To get this branch, use:
bzr branch
http://darksoft.org/webbzr/perf/fdk
« back to all changes in this revision
Viewing changes to orig/fdk_loop_suren.cpp
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Suren A. Chilingaryan
- Date: 2017-02-08 19:08:48 UTC
- Revision ID: csa@suren.me-20170208190848-lycd9dm701hyu608
Initial
added
removed
orig/fdk_loop_suren.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 "/opt/intel/ipp/include/ippi.h"
14
#include "/opt/intel/ipp/include/ipps.h"
15
#include "/opt/intel/ipp/include/ippcore.h"
16
17
18
/* structure to pass arguments to thread */
19
struct thread_info /* Used as argument to thread_start() */
20
{
21
pthread_t thread_id; /* ID returned by pthread_create() */
22
int thread_num; /* Application-defined thread # */
23
int n_elements;
24
int n_proj;
25
Ipp32f cor;
26
Ipp32f d;
27
Ipp32f detector_size;
28
Ipp32f pixel_size;
29
Ipp32f cor_offset;
30
Ipp32f detector_offset_u;
31
Ipp32f detector_offset_v;
32
Ipp32f detector_offset_z;
33
Ipp32f *source;
34
Ipp32f *u_detector;
35
Ipp32f *v_detector;
36
Ipp32f *n_detector;
37
Ipp32f *projections;
38
Ipp32f *slice_x;
39
Ipp32f *slice_y;
40
Ipp32f *slice_coord_z;
41
Ipp32f *out_volume;
42
};
43
44
/* global variables */
45
pthread_mutex_t mutex;
46
int counter = 0;
47
48
49
void statusinfo(IppStatus st)
50
{
51
if ((int)st != 0)
52
{
53
printf("%d : %s\n", st, ippGetStatusString(st));
54
}
55
}
56
57
58
bool is_aligned(void *p, int N)
59
{
60
return (uintptr_t)p % N == 0;
61
}
62
63
64
/* canonical multiplication of square matrices */
65
int mult(Ipp32f *a, Ipp32f *b, Ipp32f *c)
66
{
67
int i, j, k;
68
69
for (k = 0; k < 4; k++)
70
{
71
for (i = 0; i < 4; i++)
72
{
73
for (j = 0; j < 4; j++)
74
{
75
c[i * 4 + j] += a[i * 4 + k] * b[k * 4 + j];
76
}
77
}
78
}
79
80
return 0;
81
} /* mult */
82
83
84
void *process (void *args)
85
{
86
struct thread_info *t_info = (thread_info*) args;
87
88
Ipp32f *px_map = NULL, *py_map = NULL;
89
Ipp32f *tmp_2 = NULL;
90
Ipp32f *current_slice = NULL;
91
92
Ipp32f angle, z, w_x, w_y, w_z;
93
94
IppiSize im_size;
95
IppiRect im_roi_size;
96
97
int imStepBytes, ippStepBytes;
98
99
int i_angle, slice_number;
100
101
long idx;
102
103
/* image size */
104
im_size = {t_info->n_elements, t_info->n_elements};
105
106
/* region of interest = radisograph and slice size */
107
im_roi_size = {0, 0, t_info->n_elements, t_info->n_elements};
108
109
/* step in bytes */
110
imStepBytes = t_info->n_elements * sizeof(float);
111
112
/* allocate temporal array */
113
tmp_2 = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
114
if (tmp_2 == NULL)
115
{
116
printf("Cannot allocate tmp_2");
117
exit(-1);
118
}
119
120
current_slice = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
121
if (current_slice == NULL)
122
{
123
printf("Cannot allocate current_slice");
124
exit(-1);
125
}
126
127
/* allocate interpolation maps */
128
px_map = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
129
if (px_map == NULL)
130
{
131
printf("Cannot allocate px_map");
132
exit(-1);
133
}
134
135
py_map = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
136
if (py_map == NULL)
137
{
138
printf("Cannot allocate py_map");
139
exit(-1);
140
}
141
142
143
while (1)
144
{
145
/* counter increment */
146
retry:
147
int err = pthread_mutex_lock(&mutex);
148
if (err)
149
{
150
printf("Retrying\n");
151
goto retry;
152
}
153
154
slice_number = counter++;
155
//if (counter > t_info->n_elements)
156
if (counter > 2)
157
{
158
pthread_mutex_unlock(&mutex);
159
break;
160
}
161
162
pthread_mutex_unlock(&mutex);
163
164
/* z coordinate of current slice */
165
z = t_info->slice_coord_z[slice_number];
166
167
/* set current slice to zero */
168
statusinfo(ippiSet_32f_C1R((Ipp32f)0, current_slice, ippStepBytes, im_size));
169
170
for (i_angle = 0; i_angle < t_info->n_proj; i_angle++)
171
{
172
/* set temporal variable to zero */
173
statusinfo(ippiSet_32f_C1R((Ipp32f)0, tmp_2, ippStepBytes, im_size));
174
175
/* current rotation angle */
176
angle = 2*M_PI/t_info->n_proj*i_angle;
177
178
/* set some matrices to calculate forward projection matrix operator in homogeneous coordinates*/
179
Ipp32f P1[4 * 4] = {(t_info->d + t_info->detector_offset_z) / t_info->pixel_size, 0, (t_info->detector_size / 2 - t_info->detector_offset_u) / t_info->pixel_size, 0,
180
0, -(t_info->d + t_info->detector_offset_z) / t_info->pixel_size, (t_info->detector_size / 2 - t_info->detector_offset_v) / t_info->pixel_size, 0,
181
0, 0, 1, 0,
182
0, 0, 0, 1};
183
184
Ipp32f P2[4 * 4] = {t_info->u_detector[0], t_info->u_detector[1], t_info->u_detector[2], 0,
185
t_info->v_detector[0], t_info->v_detector[1], t_info->v_detector[2], 0,
186
t_info->n_detector[0], t_info->n_detector[1], t_info->n_detector[2], 0,
187
0, 0, 0, 1};
188
189
Ipp32f P3[4 * 4] = {1, 0, 0, -t_info->source[0],
190
0, 1, 0, -t_info->source[1],
191
0, 0, 1, -t_info->source[2],
192
0, 0, 0, 1};
193
194
Ipp32f P4[4 * 4] = {1, 0, 0, 0,
195
0, 1, 0, 0,
196
0, 0, 1, t_info->cor + t_info->cor_offset,
197
0, 0, 0, 1};
198
199
Ipp32f P5[4 * 4] = {(Ipp32f)(cosf(angle)), 0, (Ipp32f)(-sinf(angle)), 0,
200
0, 1, 0, 0,
201
(Ipp32f)(sinf(angle)), 0, (Ipp32f)(cosf(angle)), 0,
202
0, 0, 0, 1};
203
204
/* set to zero temporal arrays */
205
Ipp32f P_tmp_1[4 * 4] = {0};
206
Ipp32f P_tmp_2[4 * 4] = {0};
207
Ipp32f P_tmp_3[4 * 4] = {0};
208
Ipp32f P[4 * 4] = {0};
209
210
/* forward projection operator */
211
mult(P1, P2, P_tmp_1);
212
mult(P_tmp_1, P3, P_tmp_2);
213
mult(P_tmp_2, P4, P_tmp_3);
214
mult(P_tmp_3, P5, P);
215
216
for (int i = 0; i < t_info->n_elements; i++)
217
{
218
for (int j = 0; j < t_info->n_elements; j++)
219
{
220
idx = i * t_info->n_elements + j;
221
222
w_x = P[0] * t_info->slice_x[idx] + P[1] * t_info->slice_y[idx] + P[2] * z + P[3];
223
w_y = P[4] * t_info->slice_x[idx] + P[5] * t_info->slice_y[idx] + P[6] * z + P[7];
224
w_z = P[8] * t_info->slice_x[idx] + P[9] * t_info->slice_y[idx] + P[10] * z + P[11];
225
226
px_map[idx] = w_x / w_z;
227
py_map[idx] = w_y / w_z;
228
}
229
}
230
231
/* interpolate */
232
statusinfo(ippiRemap_32f_C1R(t_info->projections + (ippStepBytes / sizeof(float) * (long)t_info->n_elements * i_angle), im_size, ippStepBytes, im_roi_size, px_map, ippStepBytes, py_map, ippStepBytes, tmp_2, ippStepBytes, im_size, IPPI_INTER_LINEAR));
233
234
/* accumulate */
235
/* in-place addition of two matrices */
236
statusinfo(ippiAdd_32f_C1IR(tmp_2, ippStepBytes, current_slice, ippStepBytes, im_size));
237
}
238
239
/* write current slice to final array */
240
//memcpy(t_info->out_volume + (t_info->n_elements * t_info->n_elements * slice_number), current_slice, t_info->n_elements * t_info->n_elements * sizeof(float));
241
statusinfo(ippiCopy_32f_C1R(current_slice, ippStepBytes, t_info->out_volume + ((long)t_info->n_elements * t_info->n_elements * slice_number), imStepBytes, im_size));
242
}
243
244
/* free memory */
245
ippiFree(px_map);
246
ippiFree(py_map);
247
ippiFree(tmp_2);
248
ippiFree(current_slice);
249
} /* process */
250
251
252
253
void mexFunction( int nlhs, mxArray *plhs[],
254
int nrhs, const mxArray *prhs[])
255
{
256
/* prhs[0] int n_voxels */
257
/* prhs[1] int n_proj */
258
/* prhs[2] float cor */
259
/* prhs[3] float d */
260
/* prhs[4] float detector_size */
261
/* prhs[5] float pixel_size */
262
/* prhs[6] float cor_offset */
263
/* prhs[7] float detector_offset_u */
264
/* prhs[8] float detector_offset_v */
265
/* prhs[9] float detector_offset_z */
266
/* prhs[10] float s */
267
/* prhs[11] float u_detector */
268
/* prhs[12] float v_detector */
269
/* prhs[13] float n_detector_r */
270
/* prhs[14] float projections */
271
/* prhs[15] float slice_x */
272
/* prhs[16] float slice_y */
273
/* prhs[17] float slice_coord_z */
274
/* prhs[18] int n_threads */
275
276
int n_proj, n_threads;
277
int n_elements;
278
Ipp32f cor, d, detector_size, pixel_size, cor_offset, detector_offset_u, detector_offset_v, detector_offset_z;
279
Ipp32f *source = NULL, *u_detector = NULL, *v_detector = NULL, *n_detector = NULL;
280
Ipp32f *projections = NULL;
281
Ipp32f *slice_x = NULL, *slice_y = NULL, *slice_coord_z = NULL;
282
Ipp32f *out_volume = NULL;
283
284
int tnum, m;
285
struct thread_info *t_info;
286
void *res;
287
288
int alignment = 32;
289
290
//////////// PARSE INPUT //////////////////
291
292
/* radiograph and slice size */
293
/* check n_voxels data type */
294
if (mxGetNumberOfElements(prhs[0])!=1 || !mxIsInt32(prhs[0]))
295
{
296
mexErrMsgTxt("mex: Input n_voxels is not scalar or int32.");
297
}
298
299
n_elements = (int)mxGetScalar(prhs[0]);
300
301
if (n_elements % 16 != 0)
302
{
303
printf("n_elements has to be multiple of 16 since IPPI uses 64 bytes alignment");
304
exit(-1);
305
}
306
307
308
/* number of projections */
309
/* check n_proj data type */
310
if (mxGetNumberOfElements(prhs[1])!=1 || !mxIsInt32(prhs[1]))
311
{
312
mexErrMsgTxt("mex: Input n_proj is not scalar or int32.");
313
}
314
315
n_proj = (int)mxGetScalar(prhs[1]);
316
317
318
/* source to rotation distance */
319
/* check input cor data type */
320
if (!mxIsSingle(prhs[2]) || mxIsComplex(prhs[2]) || mxGetNumberOfElements(prhs[2])!=1)
321
{
322
mexErrMsgTxt("mex: Input cor must be scalar and have type single.");
323
}
324
325
cor = (Ipp32f)mxGetScalar(prhs[2]);
326
327
/* source to detector distance */
328
/* check input d data type */
329
if (!mxIsSingle(prhs[3]) || mxIsComplex(prhs[3]) || mxGetNumberOfElements(prhs[3])!=1)
330
{
331
mexErrMsgTxt("mex: Input d must be scalar and have type single.");
332
}
333
334
d = (Ipp32f)mxGetScalar(prhs[3]);
335
336
337
/* physical size of detector (assumed to be square) */
338
/* check input detector_size data type */
339
if (!mxIsSingle(prhs[4]) || mxIsComplex(prhs[4]) || mxGetNumberOfElements(prhs[4])!=1)
340
{
341
mexErrMsgTxt("mex: Input detector_size must be scalar and have type single.");
342
}
343
344
detector_size = (Ipp32f)mxGetScalar(prhs[4]);
345
346
347
/* physical pixel size */
348
/* check input pixel_size data type */
349
if (!mxIsSingle(prhs[5]) || mxIsComplex(prhs[5]) || mxGetNumberOfElements(prhs[5])!=1)
350
{
351
mexErrMsgTxt("mex: Input pixel_size must be scalar and have type single.");
352
}
353
354
pixel_size = (Ipp32f)mxGetScalar(prhs[5]);
355
356
357
/* offset in cor */
358
/* check input cor_offset data type */
359
if (!mxIsSingle(prhs[6]) || mxIsComplex(prhs[6]) || mxGetNumberOfElements(prhs[6])!=1)
360
{
361
mexErrMsgTxt("mex: Input cor_offset must be scalar and have type single.");
362
}
363
364
cor_offset = (Ipp32f)mxGetScalar(prhs[6]);
365
366
367
/* detector offset in U direction */
368
/* check input detector_offset_u data type */
369
if (!mxIsSingle(prhs[7]) || mxIsComplex(prhs[7]) || mxGetNumberOfElements(prhs[7])!=1)
370
{
371
mexErrMsgTxt("mex: Input detector_offset_u must be scalar and have type single.");
372
}
373
374
detector_offset_u = (Ipp32f)mxGetScalar(prhs[7]);
375
376
377
/* detector offset in V direction */
378
/* check input detector_offset_v data type */
379
if (!mxIsSingle(prhs[8]) || mxIsComplex(prhs[8]) || mxGetNumberOfElements(prhs[8])!=1)
380
{
381
mexErrMsgTxt("mex: Input detector_offset_v must be scalar and have type single.");
382
}
383
384
detector_offset_v = (Ipp32f)mxGetScalar(prhs[8]);
385
386
387
/* detector offset along magnification line */
388
/* check input detector_offset_z data type */
389
if (!mxIsSingle(prhs[9]) || mxIsComplex(prhs[9]) || mxGetNumberOfElements(prhs[9])!=1)
390
{
391
mexErrMsgTxt("mex: Input detector_offset_z must be scalar and have type single.");
392
}
393
394
detector_offset_z = (Ipp32f)mxGetScalar(prhs[9]);
395
396
397
/* source position */
398
/* check input source data type */
399
if (!mxIsSingle(prhs[10]) || mxIsComplex(prhs[10]) || mxGetNumberOfElements(prhs[10])!=3)
400
{
401
mexErrMsgTxt("mex: Input source must be type single and contain 3 elements.");
402
}
403
404
source = (Ipp32f *)mxGetData(prhs[10]);
405
406
if (source == NULL)
407
{
408
printf("Cannot get source");
409
exit(-1);
410
}
411
412
413
/* detector coordinate frame */
414
/* check input u_detector data type */
415
if( !mxIsSingle(prhs[11]) || mxIsComplex(prhs[11]) || mxGetNumberOfElements(prhs[11])!=3)
416
{
417
mexErrMsgTxt("mex: Input u_detector must be type single and contain 3 elements.");
418
}
419
420
/* check input v_detector data type */
421
if( !mxIsSingle(prhs[12]) || mxIsComplex(prhs[12]) || mxGetNumberOfElements(prhs[12])!=3)
422
{
423
mexErrMsgTxt("mex: Input v_detector must be type single and contain 3 elements.");
424
}
425
426
/* check input n_detector data type */
427
if( !mxIsSingle(prhs[13]) || mxIsComplex(prhs[13]) || mxGetNumberOfElements(prhs[13])!=3)
428
{
429
mexErrMsgTxt("mex: Input n_detector must be type single and contain 3 elements.");
430
}
431
432
u_detector = (Ipp32f *)mxGetData(prhs[11]);
433
434
if (u_detector == NULL)
435
{
436
printf("Cannot get u_detector");
437
exit(-1);
438
}
439
440
v_detector = (Ipp32f *)mxGetData(prhs[12]);
441
442
if (u_detector == NULL)
443
{
444
printf("Cannot get v_detector");
445
exit(-1);
446
}
447
448
n_detector = (Ipp32f *)mxGetData(prhs[13]);
449
450
if (n_detector == NULL)
451
{
452
printf("Cannot get n_detector");
453
exit(-1);
454
}
455
456
457
/* volume with projections, x-y-angle */
458
/* check input projections data type */
459
if (!mxIsSingle(prhs[14]) || mxIsComplex(prhs[14]))
460
{
461
mexErrMsgTxt("mex: Input projections array must be type single.");
462
}
463
464
/* check that number of elements in projections is equal to n_voxels x n_voxels x n_proj */
465
if ((mxGetNumberOfElements(prhs[14]) != (long)n_elements * n_elements * n_proj))
466
{
467
mexErrMsgTxt("mex: Input projections must be a n_elements x n_elements x n_proj array.");
468
}
469
470
projections = (Ipp32f *)mxGetData(prhs[14]);
471
472
if (projections == NULL) {
473
printf("Can not get projections");
474
exit(-1);
475
}
476
477
/* check if dataset is aligned by 32 bytes */
478
if (is_aligned(projections, (int)alignment) == 0)
479
{
480
printf("Projections have to be aligned by 32 bytes\n");
481
exit(-1);
482
}
483
484
485
/* grids */
486
/* check input slice_x data type */
487
if ( !mxIsSingle(prhs[15]) || mxIsComplex(prhs[15]))
488
{
489
mexErrMsgTxt("mex: Input slice_x array must be type single.");
490
}
491
492
/* check that number of rows and columns in slice_x is equal to n_voxels */
493
if(mxGetM(prhs[15]) != n_elements || mxGetN(prhs[15]) != n_elements)
494
{
495
mexErrMsgTxt("mex: Input slice_x must be a n_voxels x n_elements array.");
496
}
497
498
/* check input slice_y data type */
499
if( !mxIsSingle(prhs[16]) || mxIsComplex(prhs[16]))
500
{
501
mexErrMsgTxt("mex: Input slice_y array must be type single.");
502
}
503
504
/* check that number of rows and columns in slice_y is equal to n_voxels */
505
if(mxGetM(prhs[16]) != n_elements || mxGetN(prhs[16]) != n_elements)
506
{
507
mexErrMsgTxt("mex: Input slice_y must be a n_voxels x n_elements array.");
508
}
509
510
/* check input slice_coord_z data type */
511
if( !mxIsSingle(prhs[17]) || mxIsComplex(prhs[17]))
512
{
513
mexErrMsgTxt("mex: Input slice_coord_z array must be type single.");
514
}
515
516
/* check that number of elements in slice_coord_z is equal to n_voxels */
517
if(mxGetNumberOfElements(prhs[17]) != n_elements)
518
{
519
mexErrMsgTxt("mex: Input slice_coord_z must be a n_voxels array.");
520
}
521
522
slice_x = (Ipp32f *)mxGetData(prhs[15]);
523
524
if (slice_x == NULL)
525
{
526
printf("Cannot get slice_x");
527
exit(-1);
528
}
529
530
/* check if slice_x is aligned by 32 bytes */
531
if (is_aligned(slice_x, (int)alignment) == 0)
532
{
533
printf("slice_x have to be aligned by 32 bytes\n");
534
exit(-1);
535
}
536
537
slice_y = (Ipp32f *)mxGetData(prhs[16]);
538
539
if (slice_y == NULL)
540
{
541
printf("Cannot get slice_y");
542
exit(-1);
543
}
544
545
/* check if slice_y is aligned by 32 bytes */
546
if (is_aligned(slice_y, (int)alignment) == 0)
547
{
548
printf("slice_y have to be aligned by 32 bytes\n");
549
exit(-1);
550
}
551
552
slice_coord_z = (Ipp32f *)mxGetData(prhs[17]);
553
554
if (slice_coord_z == NULL) {
555
printf("Cannot get slice_coord_z");
556
exit(-1);
557
}
558
559
560
/* n_threads */
561
/* check input n_threads data type */
562
if(mxGetNumberOfElements(prhs[18])!=1 || !mxIsInt32(prhs[18]))
563
{
564
mexErrMsgTxt("mex: Input n_threads is not scalar or int32.");
565
}
566
567
n_threads = (int)mxGetScalar(prhs[18]);
568
569
570
/* create the output matrix im_size x im_size*/
571
mwSize dims[3];
572
dims[0] = n_elements;
573
dims[1] = n_elements;
574
dims[2] = n_elements;
575
plhs[0] = mxCreateNumericArray((mwSize)3,dims,mxSINGLE_CLASS,mxREAL);
576
577
/* get a pointer to the data in the output matrix */
578
out_volume = (Ipp32f *)mxGetData(plhs[0]);
579
580
if (out_volume == NULL) {
581
printf("Cannot get out_volume");
582
exit(-1);
583
}
584
585
///////// INITIALIZE THREADS //////////////
586
587
/* set global variable counter to zero */
588
counter = 0;
589
590
/* initialize mutex */
591
m = pthread_mutex_init(&mutex, NULL);
592
593
if (m != 0)
594
{
595
printf("Mutex initialization failed with error %i\n", m);
596
exit(-1);
597
}
598
599
/* allocate memory for threads */
600
t_info = (thread_info*) calloc(n_threads, sizeof(struct thread_info));
601
602
if (t_info == NULL)
603
{
604
printf("Structure t_info memory allocation error\n");
605
exit(-1);
606
}
607
608
/* fill tread_info structure and create threads */
609
for (tnum = 0; tnum < n_threads; tnum++)
610
{
611
t_info[tnum].thread_num = tnum; /* Application-defined thread # */
612
t_info[tnum].n_elements = n_elements;
613
t_info[tnum].n_proj = n_proj;
614
t_info[tnum].cor = cor;
615
t_info[tnum].d = d;
616
t_info[tnum].detector_size = detector_size;
617
t_info[tnum].pixel_size = pixel_size;
618
t_info[tnum].cor_offset = cor_offset;
619
t_info[tnum].detector_offset_u = detector_offset_u;
620
t_info[tnum].detector_offset_v = detector_offset_v;
621
t_info[tnum].detector_offset_z = detector_offset_z;
622
t_info[tnum].source = source;
623
t_info[tnum].u_detector = u_detector;
624
t_info[tnum].v_detector = v_detector;
625
t_info[tnum].n_detector = n_detector;
626
t_info[tnum].projections = projections;
627
t_info[tnum].slice_x = slice_x;
628
t_info[tnum].slice_y = slice_y;
629
t_info[tnum].slice_coord_z = slice_coord_z;
630
t_info[tnum].out_volume = out_volume;
631
632
m = pthread_create(&t_info[tnum].thread_id, NULL, &process, &t_info[tnum]);
633
634
if (m != 0)
635
{
636
printf("Cannot create thread with error %i\n", m);
637
exit(-1);
638
}
639
}
640
641
/* join threads */
642
for (tnum = 0; tnum < n_threads; tnum++)
643
{
644
m = pthread_join(t_info[tnum].thread_id, &res);
645
646
if (m != 0)
647
{
648
printf("Cannot join threads with error %i\n", m);
649
exit(-1);
650
}
651
}
652
653
/* free memory */
654
pthread_mutex_destroy(&mutex);
655
656
free(t_info);
657
658
} /* mexFunction */
b'\\ No newline at end of file'
Loggerhead is a web-based interface for Breezy
Version: 1.20.0