45 unsigned int myID = thrParams->myID;
46 unsigned int numThreads = thrParams->numThreads;
47 double sigma_s = thrParams->sigma_s;
48 double sigma_r = thrParams->sigma_r;
51 bool fast = thrParams->fast;
56 sigma_s =
CEIL( sigma_s );
61 sigma_s =
CEIL( sigma_s );
65 double inv_2_sigma_s_2 = 1.0/( 2.0* sigma_s * sigma_s);
66 double inv_2_sigma_r_2 = 1.0/( 2.0* sigma_r * sigma_r);
67 double sigma_r_2 = sigma_r * sigma_r;
68 auto _3_sigma_s = (int)(3 * sigma_s);
69 double _3_sigma_r = 3.0 * sigma_r;
76 int curr_x, curr_y, curr_z;
77 int prev_x, prev_y, prev_z;
80 int myFirstY, myLastY;
82 int numElems = y_max-y_min+1;
83 numElems /= numThreads;
85 myFirstY = myID * numElems + y_min;
87 if( myID == numThreads -1 )
90 myLastY = myFirstY + numElems - 1;
94 std::cerr <<
"Progress (over a total of " << z_max - z_min + 1 <<
" slices): ";
99 for (
int k=z_min;
k<=z_max;
k++)
104 std::cerr <<
k - z_min <<
" ";
106 std::cerr <<
k + z_min <<
" ";
109 for (
int i=myFirstY;
i<=myLastY;
i++)
111 for (
int j=x_min;
j<=x_max;
j++)
121 int xcOld, ycOld, zcOld;
141 for(
int z_j = -(
int)sigma_s ; z_j <=(int)sigma_s ; z_j++ )
144 if( z2 >= z_min && z2 <= z_max)
146 int z_j_2 = z_j * z_j;
148 for(
int y_j = -(
int)sigma_s ; y_j <= (int)sigma_s ; y_j++ )
151 if( y2 >= y_min && y2 <= y_max)
153 int y_j_2 = y_j * y_j;
155 for(
int x_j = -(
int)sigma_s ; x_j <= (int)sigma_s ; x_j++ )
158 if( x2 >= x_min && x2 <= x_max)
161 if( x_j*x_j+y_j_2+z_j_2 <= sigma_s*sigma_s )
163 double Y2 =
A3D_ELEM( input, z2, y2, x2);
166 if( dY*dY <= sigma_r_2 )
182 double num_ = 1.0/num;
184 xc = (int) (mx*num_+0.5);
185 yc = (int) (my*num_+0.5);
186 zc = (int) (mz*num_+0.5);
191 double dY = Yc-YcOld;
193 shift = dx*dx+dy*dy+dz*dz*dY*dY;
196 while(shift>3 && iters<100);
205 for (
int k=z_min;
k<=z_max;
k++)
208 std::cerr <<
k <<
" ";
209 for (
int i=myFirstY;
i<=myLastY;
i++)
211 for (
int j=x_min;
j<=x_max;
j++)
227 double x_sum = 0, y_sum = 0, z_sum = 0;
231 for(
int z_j = curr_z - _3_sigma_s ; z_j <= curr_z + _3_sigma_s ; z_j++ )
233 if( z_j >= z_min && z_j <= z_max)
235 for(
int y_j = curr_y - _3_sigma_s ; y_j <= curr_y + _3_sigma_s ; y_j++ )
237 if( y_j >= y_min && y_j <= y_max)
239 for(
int x_j = curr_x - _3_sigma_s ; x_j <= curr_x + _3_sigma_s ; x_j++ )
241 if( x_j >= x_min && x_j <= x_max)
243 double I_j =
A3D_ELEM( input, z_j, y_j, x_j);
244 double I_dist=fabs(I_j - curr_I);
246 if( I_dist <= _3_sigma_r )
249 double eucl_dist = (curr_x - x_j)*(curr_x - x_j)+(curr_y - y_j)*(curr_y - y_j)+(curr_z - z_j)*(curr_z - z_j);
250 double aux = exp(-(eucl_dist*inv_2_sigma_s_2) - (I_dist*I_dist*inv_2_sigma_r_2));
265 prev_x = (int)
round(curr_x);
266 prev_y = (int)
round(curr_y);
267 prev_z = (int)
round(curr_z);
269 double isum_denom=1.0/sum_denom;
270 curr_x = (int)
round(x_sum*isum_denom);
271 curr_y = (int)
round(y_sum*isum_denom);
272 curr_z = (int)
round(z_sum*isum_denom);
273 curr_I = I_sum*isum_denom;
275 error = fabs(prev_x - curr_x) + fabs(prev_y - curr_y) + fabs(prev_z - curr_z);
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define A3D_ELEM(V, k, i, j)