/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 ippprocess.c

  • Committer: Suren A. Chilingaryan
  • Date: 2017-02-09 00:44:25 UTC
  • Revision ID: csa@suren.me-20170209004425-4dt67qhxz9ibdehy
Intel compiler

Show diffs side-by-side

added added

removed removed

Lines of Context:
ippprocess.c
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
 
 
4
#ifndef __USE_BSD
 
5
# define __USE_BSD
 
6
#endif
 
7
#include <math.h>
 
8
 
 
9
#include <pthread.h>
 
10
#include <stdint.h>
 
11
#include <string.h>
 
12
 
 
13
#include <limits.h>
 
14
 
 
15
#include <ippi.h>
 
16
#include <ipps.h>
 
17
#include <ippcore.h>
 
18
 
 
19
#include "process.h"
 
20
 
 
21
/* global variables */
 
22
int counter = 0;
 
23
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
 
24
 
 
25
 
 
26
void statusinfo(IppStatus st) 
 
27
{
 
28
    if ((int)st != 0)
 
29
    {
 
30
        printf("%d : %s\n", st, ippGetStatusString(st));
 
31
    }
 
32
}  
 
33
 
 
34
 
 
35
 
 
36
 
 
37
/* canonical multiplication of square matrices */
 
38
int mult(Ipp32f *a, Ipp32f *b, Ipp32f *c) 
 
39
{
 
40
    int i, j, k;
 
41
        
 
42
    for (k = 0; k < 4; k++) 
 
43
    {
 
44
        for (i = 0; i < 4; i++) 
 
45
        {
 
46
            for (j = 0; j < 4; j++) 
 
47
            {
 
48
                c[i * 4 + j] += a[i * 4 + k] * b[k * 4 + j];
 
49
            }
 
50
        }
 
51
    }
 
52
 
 
53
    return 0;
 
54
} /* mult */
 
55
 
 
56
void *process (void *args) 
 
57
{
 
58
    struct thread_info *t_info = (struct thread_info*) args;
 
59
    
 
60
    Ipp32f *px_map = NULL, *py_map = NULL;
 
61
    Ipp32f *tmp_1 = NULL, *tmp_2 = NULL, *w_x = NULL, *w_y = NULL, *w_z = NULL;
 
62
    Ipp32f *current_slice = NULL;
 
63
    
 
64
    Ipp32f angle, z;
 
65
    
 
66
    IppiSize im_size = {t_info->n_elements, t_info->n_elements};
 
67
    IppiRect im_roi_size = {0, 0, t_info->n_elements, t_info->n_elements};
 
68
    
 
69
    int imStepBytes, ippStepBytes;
 
70
    
 
71
    int i_angle, slice_number;
 
72
    
 
73
    /* step in bytes */
 
74
    imStepBytes = t_info->n_elements * sizeof(float);
 
75
 
 
76
    /* allocate temporal arrays */
 
77
    tmp_1 = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
78
    if (tmp_1 == NULL) 
 
79
    {
 
80
        printf("Cannot allocate tmp_1");
 
81
        exit(-1);
 
82
    }
 
83
    
 
84
    tmp_2 = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
85
    if (tmp_2 == NULL) 
 
86
    {
 
87
        printf("Cannot allocate tmp_2");
 
88
        exit(-1);
 
89
    }
 
90
    
 
91
    w_x = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
92
    if (w_x == NULL) 
 
93
    {
 
94
        printf("Cannot allocate w_x");
 
95
        exit(-1);
 
96
    }
 
97
    
 
98
    w_y = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
99
    if (w_y == NULL) 
 
100
    {
 
101
        printf("Cannot allocate w_y");
 
102
        exit(-1);
 
103
    }
 
104
    
 
105
    w_z = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
106
    if (w_z == NULL) 
 
107
    {
 
108
        printf("Cannot allocate w_z");
 
109
        exit(-1);
 
110
    }
 
111
    
 
112
    current_slice = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
113
    if (current_slice == NULL) 
 
114
    {
 
115
        printf("Cannot allocate current_slice");
 
116
        exit(-1);
 
117
    }
 
118
    
 
119
    /* allocate interpolation maps */
 
120
    px_map = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
121
    if (px_map == NULL) 
 
122
    {
 
123
        printf("Cannot allocate px_map");
 
124
        exit(-1);
 
125
    }
 
126
    
 
127
    py_map = ippiMalloc_32f_C1(t_info->n_elements, t_info->n_elements, &ippStepBytes);
 
128
    if (py_map == NULL) 
 
129
    {
 
130
        printf("Cannot allocate py_map");
 
131
        exit(-1);
 
132
    }
 
133
        
 
134
    
 
135
    while (1)
 
136
    {
 
137
        /* counter increment */
 
138
retry: 
 
139
        if (pthread_mutex_lock(&mutex))
 
140
        {
 
141
            printf("Retrying\n");
 
142
            goto retry;
 
143
        }
 
144
 
 
145
        slice_number = counter;
 
146
 
 
147
        if (slice_number >= t_info->n_elements) 
 
148
        {
 
149
            pthread_mutex_unlock(&mutex);
 
150
            break;
 
151
        }
 
152
 
 
153
        counter += t_info->slices_per_iter;
 
154
        
 
155
        pthread_mutex_unlock(&mutex);
 
156
        
 
157
        /* z coordinate of current slice */
 
158
        z = t_info->slice_coord_z[slice_number];
 
159
    
 
160
        /* set current slice to zero */ 
 
161
        statusinfo(ippiSet_32f_C1R((Ipp32f)0, current_slice, ippStepBytes, im_size));
 
162
        
 
163
        for (i_angle = 0; i_angle < t_info->n_proj; i_angle++)
 
164
        {    
 
165
            /* set temporal variable to zero */
 
166
            statusinfo(ippiSet_32f_C1R((Ipp32f)0, tmp_2, ippStepBytes, im_size));
 
167
            
 
168
            /* current rotation angle */
 
169
            angle = 2*M_PI/t_info->n_proj*i_angle;
 
170
            
 
171
            /* set some matrices to calculate forward projection matrix operator in homogeneous coordinates*/
 
172
            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,
 
173
                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,
 
174
                0, 0, 1, 0,
 
175
                0, 0, 0, 1};
 
176
        
 
177
            Ipp32f P2[4 * 4] = {t_info->u_detector[0], t_info->u_detector[1], t_info->u_detector[2], 0,
 
178
                t_info->v_detector[0], t_info->v_detector[1], t_info->v_detector[2], 0,
 
179
                t_info->n_detector[0], t_info->n_detector[1], t_info->n_detector[2], 0,
 
180
                0, 0, 0, 1};
 
181
        
 
182
            Ipp32f P3[4 * 4] = {1, 0, 0, -t_info->source[0],
 
183
                0, 1, 0, -t_info->source[1],
 
184
                0, 0, 1, -t_info->source[2],
 
185
                0, 0, 0, 1};
 
186
        
 
187
            Ipp32f P4[4 * 4] = {1, 0, 0, 0,
 
188
                0, 1, 0, 0,
 
189
                0, 0, 1, t_info->cor + t_info->cor_offset,
 
190
                0, 0, 0, 1};
 
191
        
 
192
            Ipp32f P5[4 * 4] = {(Ipp32f)(cosf(angle)), 0, (Ipp32f)(-sinf(angle)), 0,
 
193
                0, 1, 0, 0,
 
194
                (Ipp32f)(sinf(angle)), 0, (Ipp32f)(cosf(angle)), 0,
 
195
                0, 0, 0, 1};
 
196
            
 
197
            /* set to zero temporal arrays */
 
198
            Ipp32f P_tmp_1[4 * 4] = {0};
 
199
            Ipp32f P_tmp_2[4 * 4] = {0};
 
200
            Ipp32f P_tmp_3[4 * 4] = {0};
 
201
            Ipp32f P[4 * 4] = {0};
 
202
            
 
203
            /* forward projection operator */
 
204
            mult(P1, P2, P_tmp_1);
 
205
            mult(P_tmp_1, P3, P_tmp_2);
 
206
            mult(P_tmp_2, P4, P_tmp_3);
 
207
            mult(P_tmp_3, P5, P);
 
208
            
 
209
            /* forward projetion in homogeneous coordinates *
 
210
             * in MATLAB notation this is equivalent to following
 
211
             * w_x = P(1,1)*slice_x+P(1,2)*slice_y+P(1,3)*z(i)+P(1,4);
 
212
             * w_y = P(2,1)*slice_x+P(2,2)*slice_y+P(2,3)*z(i)+P(2,4);
 
213
             * w_z = P(3,1)*slice_x+P(3,2)*slice_y+P(3,3)*z(i)+P(3,4);
 
214
             * px_map = w_x./w_z;
 
215
             * py_map = w_y./w_z; */
 
216
            
 
217
            /* calculate w_x */
 
218
            /* not-in-place mutiplication of matrix by constant */
 
219
            statusinfo(ippiMulC_32f_C1R(t_info->slice_x, ippStepBytes, P[0], tmp_1, ippStepBytes, im_size));
 
220
            statusinfo(ippiMulC_32f_C1R(t_info->slice_y, ippStepBytes, P[1], w_x, ippStepBytes, im_size));
 
221
            
 
222
            /* in-place addition of two matrices*/
 
223
            statusinfo(ippiAdd_32f_C1IR(tmp_1, ippStepBytes, w_x, ippStepBytes, im_size));
 
224
            
 
225
            /* in-place addition of matrix and constant */
 
226
            statusinfo(ippiAddC_32f_C1IR(P[2]*z + P[3], w_x, ippStepBytes, im_size));
 
227
            
 
228
            /* calculate w_y */
 
229
            statusinfo(ippiMulC_32f_C1R(t_info->slice_x, ippStepBytes, P[4], tmp_1, ippStepBytes, im_size));
 
230
            statusinfo(ippiMulC_32f_C1R(t_info->slice_y, ippStepBytes, P[5], w_y, ippStepBytes, im_size));
 
231
            
 
232
            statusinfo(ippiAdd_32f_C1IR(tmp_1, ippStepBytes, w_y, ippStepBytes, im_size));
 
233
 
 
234
            statusinfo(ippiAddC_32f_C1IR(P[6]*z + P[7], w_y, ippStepBytes, im_size));
 
235
            
 
236
            /* calculate w_z */
 
237
            statusinfo(ippiMulC_32f_C1R(t_info->slice_x, ippStepBytes, P[8], tmp_1, ippStepBytes, im_size));
 
238
            statusinfo(ippiMulC_32f_C1R(t_info->slice_y, ippStepBytes, P[9], w_z, ippStepBytes, im_size));
 
239
            
 
240
            statusinfo(ippiAdd_32f_C1IR(tmp_1, ippStepBytes, w_z, ippStepBytes, im_size));
 
241
 
 
242
            statusinfo(ippiAddC_32f_C1IR(P[10]*z + P[11], w_z, ippStepBytes, im_size));
 
243
            
 
244
            /* calculate px_map and py_map = dehomogenize */
 
245
            /* not-in-place devision of two matrices */
 
246
            statusinfo(ippiDiv_32f_C1R(w_z, ippStepBytes, w_x, ippStepBytes, px_map, ippStepBytes, im_size));
 
247
            statusinfo(ippiDiv_32f_C1R(w_z, ippStepBytes, w_y, ippStepBytes, py_map, ippStepBytes, im_size));
 
248
            
 
249
            /* interpolate */
 
250
            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));
 
251
 
 
252
            /* accumulate */
 
253
            /* in-place addition of two matrices */
 
254
            statusinfo(ippiAdd_32f_C1IR(tmp_2, ippStepBytes, current_slice, ippStepBytes, im_size));
 
255
        }
 
256
        
 
257
        /* write current slice to final array */
 
258
        //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));
 
259
        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));
 
260
    }
 
261
    
 
262
    /* free memory */
 
263
    ippiFree(px_map);
 
264
    ippiFree(py_map);
 
265
    ippiFree(w_x);
 
266
    ippiFree(w_y);
 
267
    ippiFree(w_z);
 
268
    ippiFree(tmp_1);
 
269
    ippiFree(tmp_2);
 
270
    ippiFree(current_slice); 
 
271
} /* process */