~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedCoreg_covar.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :
 37  * Date    :
 38  * Version :
 39  * CVS Id  :
 40  *
 41  * Notes :
 42  *
 43  *********
 44 */
 45 
 46 #include "tlmedCoreg_covar.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #   include <config.h>
 50 #endif
 51 
 52 #include <math.h>
 53 #include <float.h>
 54 #include <tina/sys/sysDef.h>
 55 #include <tina/sys/sysPro.h>
 56 #include <tina/math/mathDef.h>
 57 #include <tina/math/mathPro.h>
 58 #include <tina/image/imgDef.h>
 59 #include <tina/image/imgPro.h>
 60 #include <tinatool/tlmedical/tlmedCoreg_view.h>
 61 #include <tinatool/tlmedical/tlmedCoreg_auto.h>
 62 #include <tinatool/tlmedical/tlmedCoreg_mui.h>
 63 #include <tinatool/tlmedical/tlmedCoreg_muipab.h>
 64 
 65 static Imrect *Rslicex=NULL, *Rslicey=NULL, *Rslicez=NULL;
 66 static double xcentre, ycentre, zcentre;
 67 static int minmask;
 68 static double min_chisq=0.0;
 69 static int param_of_interest=0;
 70 static double *global_params_vec=NULL;
 71 shistogram *image_hist=NULL;
 72 static int x_cent, y_cent;
 73 
 74 /* Coreg interface flags: see tlmedCoreg_tool.c for descriptions */
 75 
 76 extern int mi_switch;
 77 extern int mod_switch;
 78 extern int a_switch;
 79 extern int bin_switch;
 80 extern int peak_switch;
 81 
 82 
 83 static void local_deriv(shistogram *imhist, float ref_pix, float trans_pix, float *n_max, float *j_max, float *pderiv,
 84 float *pmax)
 85 {
 86         float pix0, pix1, pix2, deriv, max;
 87 
 88         pix0 = hfill2s(imhist, ref_pix, (trans_pix - imhist->yincr), 0.0);
 89         pix1 = hfill2s(imhist, ref_pix, trans_pix, 0.0);
 90         pix2 = hfill2s(imhist, ref_pix, (trans_pix + imhist->yincr), 0.0);
 91         max = (float)interpolated_max_finder(imhist, ref_pix, n_max, j_max);
 92 
 93         deriv = (pix2-pix0)/(2*imhist->yincr);
 94 
 95         *pderiv = deriv;
 96         *pmax = max;
 97 }
 98 
 99 
100 static int blank_region_blocker(float g1, float g2, shistogram *imhist)
101 {
102         /* Find out if the bin identified by the pixel values g1, g2, lies within a certain radius of the
103            array coordinates identified by x_cent, y_cent */
104 
105         int x_bin, y_bin;
106         float radius_r, radius_f, radius;
107 
108         x_bin = (int)( (g1-imhist->xmin)/imhist->xincr );
109         y_bin = (int)( (g2-imhist->ymin)/imhist->yincr );
110 /*
111         radius_r = (int) 300 / (2*imhist->xincr);
112         radius_f = (int) 300 / (2*imhist->yincr);
113 */
114 
115         radius_r = 8;
116         radius_f = 8;
117 
118         radius = sqrt(sqr((x_bin-x_cent)*1.1/radius_r)+sqr((y_bin-y_cent)*1.1/radius_f));
119 
120         if(radius<1)
121         {
122                 return 0;
123         }
124         else
125         {
126                 return 1;
127         }
128 }
129 
130 
131 static void covar_deriv_im(float cent, Imregion *roi, int mask,Imrect *im_ref, Imrect *im_tr, double **sum_cov,
132                                                 int sw, shistogram *imhist, double *step_vector, float *n_max,
133 float *j_max)
134 {
135    Imrect **im_delta_diff, **im_delta, **im_minusdelta;
136    double *params_old,*params_delta,*params_step, *cov_vec_pix;
137    int npar,i,j,k,k1,k2;
138    void *memcpy();
139    double pix,pix1,pix2,pix3,pix5,  **cov;
140    float g1, g2;
141    float max, deriv;
142    Imrect *tmp=NULL;
143 
144    params_old = dvector_alloc(0,10);
145    params_delta = dvector_alloc(0,10);
146    params_step = dvector_alloc(0,10);
147    npar = coreg_get_vec(params_old, mask);
148    cov_vec_pix = dvector_alloc(0,npar);
149 
150    im_delta = tvector_alloc(0,npar,Imrect *);
151    im_minusdelta = tvector_alloc(0,npar,Imrect *);
152    im_delta_diff = tvector_alloc(0,npar,Imrect *);
153 
154    cov = darray_alloc(0,0,npar,npar);
155    for (k1 = 0; k1 < npar; k1++)
156      for (k2 = 0; k2 < npar;k2++)
157         cov[k1][k2] = 0.0;
158 
159    for (k = 0; k <npar; k++)
160    {
161       memcpy(params_delta,params_old,npar*sizeof(double));
162 
163       params_delta[k] = step_vector[k];
164       params_step[k] = step_vector[k]-params_old[k];
165 
166       coreg_set_vec(params_delta,mask);
167       if (sw == 0)
168          im_delta[k] = (Imrect *) seq_slicex(cent,roi, coreg_bproj);
169       else if (sw == 1)
170          im_delta[k] = (Imrect *) seq_slicey(cent,roi, coreg_bproj);
171       else if (sw == 2)
172          im_delta[k] = (Imrect *) seq_slicez(cent,roi, coreg_bproj);
173       else
174          return;
175       im_delta_diff[k] = im_diff(im_delta[k],im_tr);
176 
177       for(i=0; i<3; i++)
178       {
179         tmp = im_tsmooth(im_delta_diff[k]);
180         im_free(im_delta_diff[k]);
181         im_delta_diff[k] = tmp;
182       }
183 
184    }
185    for (i = roi->ly; i < roi->uy; i++)
186    {
187       for (j = roi->lx; j < roi->ux; j++)
188       {
189          for (k = 0; k< npar; k++) cov_vec_pix[k] =  im_get_pixf(im_delta_diff[k],i,j);
190 
191          g1 = im_get_pixf(im_ref,i,j);
192          g2 = im_get_pixf(im_tr,i,j);
193 
194         if(peak_switch==0||blank_region_blocker(g1, g2, imhist)==1)
195         {
196                 pix1 = hfill2s(imhist,g1,g2,0.0);
197                 local_deriv(imhist, g1, g2, n_max, j_max, &deriv, &max);
198                 pix3 = (pix1+1.0)/(max+1.0);
199                 if(pix3>1.0) pix3=1.0;
200                 pix = -(sqr(deriv)-0.00/(2.0*imhist->yincr*imhist->yincr));
201                 pix5 = 2 * (pix1+1.0)*(pix1+1.0)*log(pix3);
202                 pix2 = pix/pix5;
203 
204                 if(pix5!=0)
205                 {
206                     for (k1 = 0; k1 < npar; k1++)
207                     {
208                        for (k2 = 0; k2 < npar;k2++)
209                        {
210                           cov[k1][k2] = cov_vec_pix[k1]*cov_vec_pix[k2]*pix2;
211                           if(cov[k1][k2]>0.0||cov[k1][k2]<=0.0)
212                           {
213                                 sum_cov[k1][k2]+=
214                                    cov[k1][k2]/(params_step[k1]*params_step[k2]);
215                           }
216                        }
217                     }
218                 }
219            }
220       }
221    }
222    coreg_set_vec(params_old, mask);
223 
224    for (i = 0; i < npar; i++)
225    {
226      im_free(im_delta[i]);
227    }
228    tvector_free(im_delta,0,Imrect *);
229    for (i = 0; i < npar; i++)  im_free(im_delta_diff[i]);
230 
231    tvector_free(im_delta_diff,0,Imrect *);
232    dvector_free(params_old,0);
233    dvector_free(params_delta,0);
234    dvector_free(params_step,0);
235    dvector_free(cov_vec_pix,0);
236    darray_free(cov, 0, 0, npar, npar);
237 }
238 
239 
240 double step_size_cost(int npar, double *a)
241 {
242         double cost = 0.0, chisq=0.0, *params_vec=NULL;
243         int i;
244 
245         params_vec = dvector_alloc(0, 10);
246         for(i=0; i<10; i++) params_vec[i] = global_params_vec[i];
247         params_vec[param_of_interest] = a[0];
248         chisq = mui_voxchisq_pab(npar, params_vec);
249         cost = fabs(100*(chisq-min_chisq)/min_chisq-2);
250         dvector_free(params_vec, 0);
251         return cost;
252 }
253 
254 
255 double *step_size_calc(int mask)
256 {
257         int npar, i;
258         double *step_vec=NULL, chisq, *init_vec=NULL, *scales=NULL, *final_vec=NULL;
259         double c_test1 = 0.001;
260 
261         global_params_vec = dvector_alloc(0, 10);
262         npar = coreg_get_vec(global_params_vec, mask);
263         step_vec = dvector_alloc(0, npar);
264         init_vec = dvector_alloc(0, 1);
265         final_vec = dvector_alloc(0, 10);
266         scales = dvector_alloc(0, 1);
267 
268         for(i=0; i<npar; i++)
269         {
270                 scales[0] = 0.1;
271                 if(i<7) scales[0] = 0.01;
272                 if(i<3) scales[0] = 1.0;
273 
274                 init_vec[0] = global_params_vec[i];
275                 param_of_interest = i;
276                 chisq = simplexmin2(1,init_vec,scales,step_size_cost,
277                         NULL,c_test1,(void(*)())format);
278                 npar = coreg_get_vec(final_vec, mask);
279                 step_vec[i] = final_vec[i];
280                 coreg_set_vec(global_params_vec, mask);
281         }
282 
283         dvector_free(init_vec, 0);
284         dvector_free(final_vec, 0);
285         dvector_free(scales, 0);
286         dvector_free(global_params_vec, 0);
287         return step_vec;
288 }
289 
290 
291 void main_peak_remover(shistogram *imhist)
292 {
293         /* Remove a peak at 0,0 in the histogram */
294 
295         int i, j, peak_i=0, radius_r, radius_f;
296         float **array=NULL;
297         float *j_max=NULL, *n_max=NULL, peak=FLT_MIN;
298 
299         array = imhist->array;
300         n_max = (float *)nvector_alloc(0,imhist->xbins,sizeof(float));
301         j_max = (float *)nvector_alloc(0,imhist->xbins,sizeof(float));
302         array_find_max(imhist->array, imhist->xbins, imhist->ybins, n_max, j_max);
303 
304         for(i=0; i<imhist->xbins; i++)
305         {
306                 if(n_max[i]>peak)
307                 {
308                         peak_i = i;
309                         peak = n_max[i];
310                 }
311         }
312 
313         x_cent = (int)peak_i;
314         y_cent = (int)j_max[peak_i];
315 /*
316         radius_r = (int) 300 / (2*imhist->xincr);
317         radius_f = (int) 300 / (2*imhist->yincr);
318 */
319 
320         radius_r = 8;
321         radius_f = 8;
322 
323         for(i=-radius_f; i<radius_f; i++)
324         {
325                 for(j=-radius_r; j<radius_r; j++)
326                 {
327 
328                         if(sqrt(sqr(((float)i)/radius_f)+sqr(((float)j)/radius_r))<1.0)
329                           array[(y_cent+i)][(x_cent+j)] = 0.0;
330                 }
331         }
332         nvector_free((void *)n_max, 0, sizeof(float *));
333         nvector_free((void *)j_max, 0, sizeof(float *));
334 }
335 
336 
337 double **coreg_covar_calc(int mask,double border, shistogram *imhist, Imrect *slicex_ref,
338                       Imrect *slicey_ref, Imrect *slicez_ref, double *step_vec)
339 {
340         int i, j, npar, outside_step_vec=0;
341         float *j_max=NULL, *n_max=NULL;
342         double **covx=NULL,**covy=NULL, **covz=NULL;
343         double *params, **sum_cov, **cov_inv;
344         Vec3 dummy;
345         Imrect *Tslicex=NULL, *Tslicey=NULL, *Tslicez=NULL;
346 
347         minmask=mask;
348         coreg_centre(&xcentre,&ycentre,&zcentre,&dummy);
349         params = dvector_alloc(0,10);
350         npar = coreg_get_vec(params, mask);
351         covx = darray_alloc(0,0,npar,npar);
352         covy = darray_alloc(0,0,npar,npar);
353         covz = darray_alloc(0,0,npar,npar);
354         sum_cov = darray_alloc(0,0,npar,npar);
355         cov_inv = darray_alloc(0,0,npar,npar);
356 
357         if(step_vec==NULL)
358         {
359                 step_vec = dvector_alloc(0, 10);
360                 outside_step_vec=1;
361                 step_vec = step_size_calc(mask);
362         }
363         for (i = 0; i < npar; i++)
364         {
365                 for (j = 0; j < npar;j++)
366                 {
367                         sum_cov[i][j] = 0.0;
368                         cov_inv[i][j] = 0.0;
369                         covx[i][j] = 0.0;
370                         covy[i][j] = 0.0;
371                         covz[i][j] = 0.0;
372                 }
373         }
374 
375         Rslicex = slicex_ref;
376         Rslicey = slicey_ref;
377         Rslicez = slicez_ref;
378 
379         image_hist = imhist;
380         min_chisq = mui_voxchisq_pab(npar, params);
381 
382         Tslicez = seq_slicez(zcentre,Rslicez->region, coreg_bproj);
383         Tslicey = seq_slicey(ycentre,Rslicey->region, coreg_bproj);
384         Tslicex = seq_slicex(xcentre,Rslicex->region, coreg_bproj);
385 
386         hreset(imhist);
387 
388         mui_joint_hist2(Rslicex,Tslicex,imhist);
389         mui_joint_hist2(Rslicey,Tslicey,imhist);
390         mui_joint_hist2(Rslicez,Tslicez,imhist);
391 
392         if(peak_switch==1) main_peak_remover(imhist);
393         hist_smooth(imhist);
394 
395         n_max = (float *)nvector_alloc(0,imhist->xbins,sizeof(float));
396         j_max = (float *)nvector_alloc(0,imhist->xbins,sizeof(float));
397         array_find_max(imhist->array, imhist->xbins,
398                 imhist->ybins, n_max, j_max);
399 
400         for (i = 0; i < npar; i++)
401         {
402                 for (j = 0; j < npar;j++)
403                 {
404                         sum_cov[i][j] = 0.0;
405                         cov_inv[i][j] = 0.0;
406                         covx[i][j] = 0.0;
407                         covy[i][j] = 0.0;
408                         covz[i][j] = 0.0;
409                 }
410         }
411         covar_deriv_im(xcentre,Rslicex->region,mask,Rslicex,
412                 Tslicex,covx,0,imhist,step_vec, n_max, j_max);
413         covar_deriv_im(ycentre,Rslicey->region,mask,Rslicey,
414                 Tslicey,covy,1,imhist,step_vec, n_max, j_max);
415         covar_deriv_im(zcentre,Rslicez->region,mask,Rslicez,
416                 Tslicez,covz,2,imhist,step_vec, n_max, j_max);
417 
418         for (i = 0; i < npar; i++)
419                 {
420                 for (j = 0; j < npar;j++)                       {
421                         sum_cov[i][j] += covx[i][j]+covy[i][j]+covz[i][j];
422                 }
423         }
424 
425 
426         for (i = 0; i < npar; i++)
427         {
428                 for (j = 0; j < npar;j++)
429                 {
430                         sum_cov[i][j] = sum_cov[i][j]/1000000;
431                 }
432         }
433 
434         if(darray_invert(sum_cov,npar)== NULL)
435         {
436                 error("Could not find invert of covariance matrix!", warning);
437         }
438         else
439         {
440                 cov_inv = darray_invert(sum_cov, npar);
441         }
442 
443         if (Tslicez != NULL) im_free(Tslicez);
444         if (Tslicey != NULL) im_free(Tslicey);
445         if (Tslicex != NULL) im_free(Tslicex);
446 
447         darray_free(sum_cov, 0, 0, npar, npar);
448         dvector_free(params,0);
449         darray_free(covx, 0, 0, npar, npar); 
450         darray_free(covy, 0, 0, npar, npar); 
451         darray_free(covz, 0, 0, npar, npar); 
452         nvector_free((void *)n_max, 0, sizeof(float *));
453         nvector_free((void *)j_max, 0, sizeof(float *));
454         if(outside_step_vec==1) dvector_free(step_vec, 0);
455         return cov_inv;
456 }
457 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.