/perf/fdk

To get this branch, use:
bzr branch http://darksoft.org/webbzr/perf/fdk

« back to all changes in this revision

Viewing changes to fdk.cpp

  • Committer: Suren A. Chilingaryan
  • Date: 2017-02-08 19:08:48 UTC
  • Revision ID: csa@suren.me-20170208190848-lycd9dm701hyu608
Initial

Show diffs side-by-side

added added

removed removed

Lines of Context:
fdk.cpp
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
 
 
4
#include <math.h>
 
5
#include <pthread.h>
 
6
#include <stdint.h>
 
7
#include <string.h>
 
8
 
 
9
#include <limits.h>
 
10
 
 
11
#include <mex.h>
 
12
 
 
13
#include <ippi.h>
 
14
#include <ipps.h>
 
15
#include <ippcore.h>
 
16
 
 
17
#include "process.h"
 
18
 
 
19
extern int counter;
 
20
 
 
21
int is_aligned(void *p, int N)
 
22
{
 
23
    return (uintptr_t)p % N == 0;
 
24
}
 
25
 
 
26
void mexFunction( int nlhs, mxArray *plhs[],
 
27
        int nrhs, const mxArray *prhs[])
 
28
{
 
29
    /* prhs[0]      int   n_voxels */
 
30
    /* prhs[1]      int   n_proj */
 
31
    /* prhs[2]      float   cor */
 
32
    /* prhs[3]      float   d */
 
33
    /* prhs[4]      float   detector_size */
 
34
    /* prhs[5]      float   pixel_size */
 
35
    /* prhs[6]      float   cor_offset */
 
36
    /* prhs[7]      float   detector_offset_u */
 
37
    /* prhs[8]      float   detector_offset_v */
 
38
    /* prhs[9]      float   detector_offset_z */
 
39
    /* prhs[10]      float   s */
 
40
    /* prhs[11]      float   u_detector */
 
41
    /* prhs[12]      float   v_detector */
 
42
    /* prhs[13]      float   n_detector_r */
 
43
    /* prhs[14]      float   projections */
 
44
    /* prhs[15]      float   slice_x */
 
45
    /* prhs[16]      float   slice_y */
 
46
    /* prhs[17]      float   slice_coord_z */
 
47
    /* prhs[18]      int     n_threads */
 
48
    
 
49
    int n_proj, n_threads;
 
50
    int n_elements;
 
51
    Ipp32f cor, d, detector_size, pixel_size, cor_offset, detector_offset_u, detector_offset_v, detector_offset_z;
 
52
    Ipp32f *source = NULL, *u_detector = NULL, *v_detector = NULL, *n_detector = NULL;
 
53
    Ipp32f *projections = NULL;
 
54
    Ipp32f *slice_x = NULL, *slice_y = NULL, *slice_coord_z = NULL;
 
55
    Ipp32f *out_volume = NULL;
 
56
    
 
57
    int tnum, m;
 
58
    struct thread_info *t_info;
 
59
    void *res;
 
60
    
 
61
    int alignment = 32;
 
62
    
 
63
    /*/////////// PARSE INPUT /////////////////*/
 
64
    
 
65
    /* radiograph and slice size */
 
66
    /* check n_voxels data type */
 
67
    if (mxGetNumberOfElements(prhs[0])!=1 || !mxIsInt32(prhs[0]))
 
68
    {
 
69
        mexErrMsgTxt("mex: Input n_voxels is not scalar or int32.");
 
70
    }
 
71
     
 
72
    n_elements = (int)mxGetScalar(prhs[0]);
 
73
    
 
74
    if (n_elements % 16 != 0)
 
75
    {
 
76
        printf("n_elements has to be multiple of 16 since IPPI uses 64 bytes alignment");
 
77
        exit(-1);
 
78
    }
 
79
 
 
80
    
 
81
    /* number of projections */
 
82
    /* check n_proj data type */
 
83
    if (mxGetNumberOfElements(prhs[1])!=1 || !mxIsInt32(prhs[1]))
 
84
    {
 
85
        mexErrMsgTxt("mex: Input n_proj is not scalar or int32.");
 
86
    }
 
87
    
 
88
    n_proj = (int)mxGetScalar(prhs[1]);
 
89
    
 
90
    
 
91
    /* source to rotation distance */
 
92
    /* check input cor data type */
 
93
    if (!mxIsSingle(prhs[2]) || mxIsComplex(prhs[2]) || mxGetNumberOfElements(prhs[2])!=1)
 
94
    {
 
95
        mexErrMsgTxt("mex: Input cor must be scalar and have type single.");
 
96
    }
 
97
    
 
98
    cor = (Ipp32f)mxGetScalar(prhs[2]);
 
99
    
 
100
    /* source to detector distance */
 
101
    /* check input d data type */
 
102
    if (!mxIsSingle(prhs[3]) || mxIsComplex(prhs[3]) || mxGetNumberOfElements(prhs[3])!=1)
 
103
    {
 
104
        mexErrMsgTxt("mex: Input d must be scalar and have type single.");
 
105
    }
 
106
    
 
107
    d = (Ipp32f)mxGetScalar(prhs[3]);
 
108
    
 
109
    
 
110
    /* physical size of detector (assumed to be square) */
 
111
    /* check input detector_size data type */
 
112
    if (!mxIsSingle(prhs[4]) || mxIsComplex(prhs[4]) || mxGetNumberOfElements(prhs[4])!=1)
 
113
    {
 
114
        mexErrMsgTxt("mex: Input detector_size must be scalar and have type single.");
 
115
    }
 
116
    
 
117
    detector_size = (Ipp32f)mxGetScalar(prhs[4]);
 
118
    
 
119
    
 
120
    /* physical pixel size */
 
121
    /* check input pixel_size data type */
 
122
    if (!mxIsSingle(prhs[5]) || mxIsComplex(prhs[5]) || mxGetNumberOfElements(prhs[5])!=1)
 
123
    {
 
124
        mexErrMsgTxt("mex: Input pixel_size must be scalar and have type single.");
 
125
    }
 
126
    
 
127
    pixel_size = (Ipp32f)mxGetScalar(prhs[5]);
 
128
   
 
129
    
 
130
    /* offset in cor */
 
131
    /* check input cor_offset data type */
 
132
    if (!mxIsSingle(prhs[6]) || mxIsComplex(prhs[6]) || mxGetNumberOfElements(prhs[6])!=1)
 
133
    {
 
134
        mexErrMsgTxt("mex: Input cor_offset must be scalar and have type single.");
 
135
    }
 
136
    
 
137
    cor_offset = (Ipp32f)mxGetScalar(prhs[6]);
 
138
    
 
139
    
 
140
    /* detector offset in U direction */
 
141
    /* check input detector_offset_u data type */
 
142
    if (!mxIsSingle(prhs[7]) || mxIsComplex(prhs[7]) || mxGetNumberOfElements(prhs[7])!=1)
 
143
    {
 
144
        mexErrMsgTxt("mex: Input detector_offset_u must be scalar and have type single.");
 
145
    }
 
146
    
 
147
    detector_offset_u = (Ipp32f)mxGetScalar(prhs[7]);
 
148
    
 
149
    
 
150
    /* detector offset in V direction */
 
151
    /* check input detector_offset_v data type */
 
152
    if (!mxIsSingle(prhs[8]) || mxIsComplex(prhs[8]) || mxGetNumberOfElements(prhs[8])!=1)
 
153
    {
 
154
        mexErrMsgTxt("mex: Input detector_offset_v must be scalar and have type single.");
 
155
    }
 
156
    
 
157
    detector_offset_v = (Ipp32f)mxGetScalar(prhs[8]);
 
158
    
 
159
    
 
160
    /* detector offset along magnification line */
 
161
    /* check input detector_offset_z data type */
 
162
    if (!mxIsSingle(prhs[9]) || mxIsComplex(prhs[9]) || mxGetNumberOfElements(prhs[9])!=1)
 
163
    {
 
164
        mexErrMsgTxt("mex: Input detector_offset_z must be scalar and have type single.");
 
165
    }
 
166
    
 
167
    detector_offset_z = (Ipp32f)mxGetScalar(prhs[9]);
 
168
    
 
169
    
 
170
    /* source position */
 
171
    /* check input source data type */
 
172
    if (!mxIsSingle(prhs[10]) || mxIsComplex(prhs[10]) || mxGetNumberOfElements(prhs[10])!=3)
 
173
    {
 
174
        mexErrMsgTxt("mex: Input source must be type single and contain 3 elements.");
 
175
    }
 
176
    
 
177
    source = (Ipp32f *)mxGetData(prhs[10]);
 
178
    
 
179
    if (source == NULL) 
 
180
    {
 
181
        printf("Cannot get source");
 
182
        exit(-1);
 
183
    }
 
184
    
 
185
    
 
186
    /* detector coordinate frame */
 
187
    /* check input u_detector data type */
 
188
    if( !mxIsSingle(prhs[11]) || mxIsComplex(prhs[11]) || mxGetNumberOfElements(prhs[11])!=3)
 
189
    {
 
190
        mexErrMsgTxt("mex: Input u_detector must be type single and contain 3 elements.");
 
191
    }
 
192
    
 
193
    /* check input v_detector data type */
 
194
    if( !mxIsSingle(prhs[12]) || mxIsComplex(prhs[12]) || mxGetNumberOfElements(prhs[12])!=3)
 
195
    {
 
196
        mexErrMsgTxt("mex: Input v_detector must be type single and contain 3 elements.");
 
197
    }
 
198
    
 
199
    /* check input n_detector data type */
 
200
    if( !mxIsSingle(prhs[13]) || mxIsComplex(prhs[13]) || mxGetNumberOfElements(prhs[13])!=3)
 
201
    {
 
202
        mexErrMsgTxt("mex: Input n_detector must be type single and contain 3 elements.");
 
203
    }
 
204
    
 
205
    u_detector = (Ipp32f *)mxGetData(prhs[11]);
 
206
    
 
207
    if (u_detector == NULL) 
 
208
    {
 
209
        printf("Cannot get u_detector");
 
210
        exit(-1);
 
211
    }
 
212
    
 
213
    v_detector = (Ipp32f *)mxGetData(prhs[12]);
 
214
    
 
215
    if (u_detector == NULL) 
 
216
    {
 
217
        printf("Cannot get v_detector");
 
218
        exit(-1);
 
219
    }
 
220
    
 
221
    n_detector = (Ipp32f *)mxGetData(prhs[13]);
 
222
    
 
223
    if (n_detector == NULL) 
 
224
    {
 
225
        printf("Cannot get n_detector");
 
226
        exit(-1);
 
227
    }
 
228
    
 
229
    
 
230
    /* volume with projections, x-y-angle */
 
231
    /* check input projections data type */
 
232
    if (!mxIsSingle(prhs[14]) || mxIsComplex(prhs[14]))
 
233
    {
 
234
        mexErrMsgTxt("mex: Input projections array must be type single.");
 
235
    }
 
236
    
 
237
    /* check that number of elements in projections is equal to n_voxels x n_voxels x n_proj */
 
238
    if ((mxGetNumberOfElements(prhs[14]) != (long)n_elements * n_elements * n_proj))
 
239
    {
 
240
        mexErrMsgTxt("mex: Input projections must be a n_elements x n_elements x n_proj array.");
 
241
    }
 
242
    
 
243
    projections = (Ipp32f *)mxGetData(prhs[14]);
 
244
    
 
245
    if (projections == NULL) {
 
246
        printf("Can not get projections");
 
247
        exit(-1);
 
248
    }
 
249
    
 
250
    /* check if dataset is aligned by 32 bytes */
 
251
    if (is_aligned(projections, (int)alignment) == 0)
 
252
    {
 
253
        printf("Projections have to be aligned by 32 bytes\n");
 
254
        exit(-1);
 
255
    }
 
256
    
 
257
    
 
258
    /* grids */
 
259
    /* check input slice_x data type */
 
260
    if ( !mxIsSingle(prhs[15]) || mxIsComplex(prhs[15]))
 
261
    {
 
262
        mexErrMsgTxt("mex: Input slice_x array must be type single.");
 
263
    }
 
264
    
 
265
    /* check that number of rows  and columns in slice_x is equal to n_voxels */
 
266
    if(mxGetM(prhs[15]) != n_elements || mxGetN(prhs[15]) != n_elements)
 
267
    {
 
268
        mexErrMsgTxt("mex: Input slice_x must be a n_voxels x n_elements array.");
 
269
    }
 
270
    
 
271
    /* check input slice_y data type */
 
272
    if( !mxIsSingle(prhs[16]) || mxIsComplex(prhs[16]))
 
273
    {
 
274
        mexErrMsgTxt("mex: Input slice_y array must be type single.");
 
275
    }
 
276
    
 
277
    /* check that number of rows  and columns in slice_y is equal to n_voxels */
 
278
    if(mxGetM(prhs[16]) != n_elements || mxGetN(prhs[16]) != n_elements)
 
279
    {
 
280
        mexErrMsgTxt("mex: Input slice_y must be a n_voxels x n_elements array.");
 
281
    }
 
282
    
 
283
    /* check input slice_coord_z data type */
 
284
    if( !mxIsSingle(prhs[17]) || mxIsComplex(prhs[17]))
 
285
    {
 
286
        mexErrMsgTxt("mex: Input slice_coord_z array must be type single.");
 
287
    }
 
288
    
 
289
    /* check that number of elements in slice_coord_z is equal to n_voxels */
 
290
    if(mxGetNumberOfElements(prhs[17]) != n_elements)
 
291
    {
 
292
        mexErrMsgTxt("mex: Input slice_coord_z must be a n_voxels array.");
 
293
    }
 
294
    
 
295
    slice_x = (Ipp32f *)mxGetData(prhs[15]);
 
296
    
 
297
    if (slice_x == NULL) 
 
298
    {
 
299
        printf("Cannot get slice_x");
 
300
        exit(-1);
 
301
    }
 
302
    
 
303
    /* check if slice_x is aligned by 32 bytes */
 
304
    if (is_aligned(slice_x, (int)alignment) == 0)
 
305
    {
 
306
        printf("slice_x have to be aligned by 32 bytes\n");
 
307
        exit(-1);
 
308
    }
 
309
    
 
310
    slice_y = (Ipp32f *)mxGetData(prhs[16]);
 
311
    
 
312
    if (slice_y == NULL) 
 
313
    {
 
314
        printf("Cannot get slice_y");
 
315
        exit(-1);
 
316
    }
 
317
    
 
318
    /* check if slice_y is aligned by 32 bytes */
 
319
    if (is_aligned(slice_y, (int)alignment) == 0)
 
320
    {
 
321
        printf("slice_y have to be aligned by 32 bytes\n");
 
322
        exit(-1);
 
323
    }
 
324
    
 
325
    slice_coord_z = (Ipp32f *)mxGetData(prhs[17]);
 
326
    
 
327
    if (slice_coord_z == NULL) {
 
328
        printf("Cannot get slice_coord_z");
 
329
        exit(-1);
 
330
    }
 
331
    
 
332
    
 
333
    /* n_threads */
 
334
    /* check input n_threads data type */
 
335
    if(mxGetNumberOfElements(prhs[18])!=1 || !mxIsInt32(prhs[18]))
 
336
    {
 
337
        mexErrMsgTxt("mex: Input n_threads is not scalar or int32.");
 
338
    }
 
339
    
 
340
    n_threads = (int)mxGetScalar(prhs[18]);
 
341
    
 
342
    
 
343
    /* create the output matrix im_size x im_size*/
 
344
    mwSize dims[3];
 
345
    dims[0] = n_elements;
 
346
    dims[1] = n_elements;
 
347
    dims[2] = n_elements;
 
348
    plhs[0] = mxCreateNumericArray((mwSize)3,dims,mxSINGLE_CLASS,mxREAL);
 
349
    
 
350
    /* get a pointer to the data in the output matrix */
 
351
    out_volume = (Ipp32f *)mxGetData(plhs[0]);
 
352
    
 
353
    if (out_volume == NULL) {
 
354
        printf("Cannot get out_volume");
 
355
        exit(-1);
 
356
    }
 
357
    
 
358
    /*//////// INITIALIZE THREADS /////////////*/
 
359
    
 
360
    /* set global variable counter to zero */
 
361
    counter = 0;
 
362
    
 
363
    /* allocate memory for threads */
 
364
    t_info = (struct thread_info*) calloc(n_threads, sizeof(struct thread_info));
 
365
    
 
366
    if (t_info == NULL)
 
367
    {
 
368
        printf("Structure t_info memory allocation error\n");
 
369
        exit(-1);
 
370
    }
 
371
    
 
372
    /* fill tread_info structure and create threads */
 
373
    for (tnum = 0; tnum < n_threads; tnum++)
 
374
    {
 
375
        t_info[tnum].thread_num = tnum;                      /* Application-defined thread # */
 
376
        t_info[tnum].n_elements = n_elements;
 
377
        t_info[tnum].n_proj = n_proj;
 
378
        t_info[tnum].cor = cor;
 
379
        t_info[tnum].d = d;
 
380
        t_info[tnum].detector_size = detector_size; 
 
381
        t_info[tnum].pixel_size = pixel_size;
 
382
        t_info[tnum].cor_offset = cor_offset;
 
383
        t_info[tnum].detector_offset_u = detector_offset_u;
 
384
        t_info[tnum].detector_offset_v = detector_offset_v;
 
385
        t_info[tnum].detector_offset_z = detector_offset_z;
 
386
        t_info[tnum].source = source;
 
387
        t_info[tnum].u_detector = u_detector; 
 
388
        t_info[tnum].v_detector = v_detector; 
 
389
        t_info[tnum].n_detector = n_detector;
 
390
        t_info[tnum].projections = projections;
 
391
        t_info[tnum].slice_x = slice_x;
 
392
        t_info[tnum].slice_y = slice_y; 
 
393
        t_info[tnum].slice_coord_z = slice_coord_z;
 
394
        t_info[tnum].out_volume = out_volume;
 
395
        
 
396
        m = pthread_create(&t_info[tnum].thread_id, NULL, &process, &t_info[tnum]);
 
397
        
 
398
        if (m != 0)
 
399
        {
 
400
            printf("Cannot create thread with error %i\n", m);
 
401
            exit(-1);
 
402
        }
 
403
    }
 
404
    
 
405
    /* join threads */
 
406
    for (tnum = 0; tnum < n_threads; tnum++)
 
407
    {
 
408
        m = pthread_join(t_info[tnum].thread_id, &res);
 
409
        
 
410
        if (m != 0)
 
411
        {
 
412
            printf("Cannot join threads with error %i\n", m);
 
413
            exit(-1);
 
414
        }
 
415
    }
 
416
 
 
417
    free(t_info);
 
418
}
 
 
b'\\ No newline at end of file'