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

Linux Cross Reference
Tina5/tina-tools/tinatool/tlmedical/tlmedCoreg_auto.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    :  $Source: /home/tina/cvs/tina-tools/tinatool/tlmedical/tlmedCoreg_auto.c,v $
 37  * Date    :  $Date: 2004/12/09 14:32:55 $
 38  * Version :  $Revision: 1.6 $
 39  * CVS Id  :  $Id: tlmedCoreg_auto.c,v 1.6 2004/12/09 14:32:55 paul Exp $
 40  *
 41  * Notes :
 42  *
 43  *********
 44 */
 45 
 46 #include "tlmedCoreg_auto.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #   include <config.h>
 50 #endif
 51 
 52 #include <stdio.h>
 53 #include <math.h>
 54 #include <float.h>
 55 #include <tina/sys/sysDef.h>
 56 #include <tina/sys/sysPro.h>
 57 #include <tina/math/mathDef.h>
 58 #include <tina/math/mathPro.h>
 59 #include <tina/image/imgDef.h>
 60 #include <tina/image/imgPro.h>
 61 #include <tina/geometry/geomDef.h>
 62 #include <tina/geometry/geomPro.h>
 63 #include <tinatool/tlbase/tlbasePro.h>
 64 #include <tinatool/tlmedical/tlmedCoreg_view.h>
 65 #include <tinatool/tlmedical/tlmedCoreg_mui.h>
 66 #include <tinatool/tlmedical/tlmedCoreg_muipab.h>
 67 #include <tinatool/tlmedical/tlmedCoreg_covar.h>
 68 
 69 
 70 /* global variables set up by init voxchisq */
 71 static Imrect *slicexx=NULL, *sliceyx=NULL, *slicezx=NULL;
 72 static Imrect *slicexy=NULL, *sliceyy=NULL, *slicezy=NULL;
 73 static Imrect *maskx=NULL, *masky=NULL, *maskz=NULL;
 74 static Imregion *roix=NULL, *roiy=NULL, *roiz=NULL;
 75 static double xcentre, ycentre, zcentre;
 76 static double c_test1 = 0.00001;
 77 static double perfect;
 78 static double smooth;
 79 static int minmask;
 80 
 81 /* Extra globals for Mutual Information coregistration  */
 82 
 83 int no_rbins=256;
 84 int no_fbins=256;
 85 static  Imrect *Rslicex=NULL, *Rslicey=NULL, *Rslicez=NULL;
 86 static  Imrect *Tslicex=NULL, *Tslicey=NULL, *Tslicez=NULL;
 87 static shistogram *hist2 = NULL;
 88 
 89 /* Coreg interface flags: see tlmedCoreg_tool for descriptions */
 90 
 91 extern int mi_switch;
 92 extern int mod_switch;
 93 extern int a_switch;
 94 extern int bin_switch;
 95 extern int bin_size_switch;
 96 extern int step_switch;
 97 
 98 Imrect *get_R_zcoreg_im()
 99 {
100   return(Rslicez);
101 }
102 
103 
104 Imrect *get_R_ycoreg_im()
105 {
106   return(Rslicey);
107 }
108 
109 
110 Imrect *get_R_xcoreg_im()
111 {
112   return(Rslicex);
113 }
114 
115 
116 static Imrect *im_prepare(Imrect *im)
117 {
118      Imrect *im1, *im2;
119      im1 = imf_lsf_smooth(im,smooth);
120      im2 = imf_lsf_smooth(im,smooth);
121      im_free(im1);
122      return(im2);
123 }
124 
125 
126 static double  im_compare(Imrect *im1, Imrect *im2, Imrect *im3,
127                           Imrect *im4, Imrect *mask, Imregion *roi)
128 {
129     int lx, ux, ly, uy;
130     int i,j;
131     float *row1, *row2, *row3 , *row4, *row5 ;
132     double dot_prod=0.0;
133     double norm=0.0;
134 
135     lx = roi->lx;
136     ux = roi->ux;
137     ly = roi->ly;
138     uy = roi->uy;
139 
140     if(mod_switch == 0)
141     {
142     /* For the data of the same modalitry */
143        for (i = ly; i < uy; ++i)
144        {
145            row1 = (float*) IM_ROW(im1,i);
146            row2 = (float*) IM_ROW(im2,i);
147            row3 = (float*) IM_ROW(im3,i);
148            row4 = (float*) IM_ROW(im4,i);
149            row5 = (float*) IM_ROW(mask,i);
150            row1 = &row1[lx];
151            row2 = &row2[lx];
152            row3 = &row3[lx];
153            row4 = &row4[lx];
154            row5 = &row5[lx];
155            for (j = lx; j < ux; ++j)
156            {
157                if (*(row5++)>0.5)
158                {
159                   norm += (*row3 * *row3) + (*row4 * *row4);
160                /* for the same modality data */
161                  dot_prod +=  *(row1++)* *(row3++) +
162                               *(row2++)* *(row4++);
163                }
164                else
165                {
166                   row1++;
167                   row2++;
168                   row3++;
169                   row4++;
170                }
171            }
172         }
173      }
174      else if(mod_switch ==1)
175      {
176     /* For the data of different modality */
177        for (i = ly; i < uy; ++i)
178        {
179            row1 = (float*) IM_ROW(im1,i);
180            row2 = (float*) IM_ROW(im2,i);
181            row3 = (float*) IM_ROW(im3,i);
182            row4 = (float*) IM_ROW(im4,i);
183            row5 = (float*) IM_ROW(mask,i);
184            row1 = &row1[lx];
185            row2 = &row2[lx];
186            row3 = &row3[lx];
187            row4 = &row4[lx];
188            row5 = &row5[lx];
189            for (j = lx; j < ux; ++j)
190            {
191                if (*(row5++)>0.5)
192                {
193                   norm += (*row3 * *row3) + (*row4 * *row4);
194 
195                /* for different modality data */
196                   dot_prod += sqrt((*(row1)**(row1++)+
197                                  *(row2)**(row2++))
198                                *(*(row3)**(row3++)+
199                                  *(row4)**(row4++)) );
200 
201                }
202                else
203                {
204                   row1++;
205                   row2++;
206                   row3++;
207                   row4++;
208                }
209            }
210        }
211     }
212     return( dot_prod/sqrt(norm) );
213 }
214 
215 
216 Imrect *crop_region(Imrect *im, double border)
217 {
218    Imrect *im1;
219    Imregion roi;
220 
221    if(im!=NULL )  roi = *(im->region);
222    else  return(NULL);
223    roi.lx += (int)border;
224    roi.ux -= (int)border;
225    roi.ly += (int)border;
226    roi.uy -= (int)border;
227    im1 = im_subim(im, &roi);
228 
229    return(im1);
230 }
231 
232 
233 void init_voxchisq(double threshold, double border)
234 {
235     Imrect *im0,*im1,*smooth_im;
236     Vec3 dummy;
237 
238     if (maskx!=NULL) im_free(maskx);
239     if (masky!=NULL) im_free(masky);
240     if (maskz!=NULL) im_free(maskz);
241     if (slicexx!=NULL) im_free(slicexx);
242     if (sliceyx!=NULL) im_free(sliceyx);
243     if (slicezx!=NULL) im_free(slicezx);
244     if (slicexy!=NULL) im_free(slicexy);
245     if (sliceyy!=NULL) im_free(sliceyy);
246     if (slicezy!=NULL) im_free(slicezy);
247     maskx = masky = maskz = NULL;
248     slicexx = sliceyx = slicezx = NULL;
249     slicexy = sliceyy = slicezy = NULL;
250 
251     coreg_centre(&xcentre, &ycentre, &zcentre, &dummy);
252     im0 = crop_region(get_xcoreg_im(),border);
253     im1 = imf_mod(im0);
254     maskx = imf_bthresh(threshold,im1);
255     im_free(im1);
256     smooth_im = im_prepare(im0);
257     im_free(im0);
258     slicexx = imf_diffx(smooth_im);
259     slicexy = imf_diffy(smooth_im);
260     im_free(smooth_im);
261 
262     im0 = crop_region(get_ycoreg_im(),border);
263     im1 = imf_mod(im0);
264     masky = imf_bthresh(threshold,im1);
265     im_free(im1);
266     smooth_im = im_prepare(im0);
267     im_free(im0);
268     sliceyx = imf_diffx(smooth_im);
269     sliceyy = imf_diffy(smooth_im);
270     im_free(smooth_im);
271 
272     im0 = crop_region(get_zcoreg_im(),border);
273     im1 = imf_mod(im0);
274     maskz = imf_bthresh(threshold,im1);
275     im_free(im1);
276     smooth_im = im_prepare(im0);
277     im_free(im0);
278     slicezx = imf_diffx(smooth_im);
279     slicezy = imf_diffy(smooth_im);
280     im_free(smooth_im);
281 
282     if (slicexx!=NULL) roix = slicexx->region;
283     else roix = NULL;
284     if (sliceyx!=NULL) roiy = sliceyx->region;
285     else roiy = NULL;
286     if (slicezx!=NULL) roiz = slicezx->region;
287     else roiz =NULL;
288     perfect = im_compare(slicexx,slicexy,slicexx,slicexy,maskx,roix);
289     perfect += im_compare(sliceyx,sliceyy,sliceyx,sliceyy,masky,roiy);
290     perfect += im_compare(slicezx,slicezy,slicezx,slicezy,maskz,roiz);
291 }
292 
293 
294 static double voxchisq(int n_par, double *a)
295 {
296      double chisq = DBL_MAX;
297      Imrect *reslicex, *reslicey, *reslicez;
298      Imrect *smooth_im;
299      Imrect *gradx, *grady;
300 
301      if (roix==NULL|| roiy==NULL || roiz == NULL) return(chisq);
302      if(coreg_set_vec(a,minmask))
303      {
304          chisq = perfect;
305          reslicex = seq_slicex(xcentre, roix, coreg_bproj);
306          smooth_im = im_prepare(reslicex);
307          im_free(reslicex);
308          gradx = imf_diffx(smooth_im);
309          grady = imf_diffy(smooth_im);
310          im_free(smooth_im);
311          chisq -= im_compare(slicexx, slicexy, gradx, grady,
312                              maskx, roix);
313          im_free(gradx);
314          im_free(grady);
315 
316          reslicey = seq_slicey(ycentre, roiy, coreg_bproj);
317          smooth_im = im_prepare(reslicey);
318          im_free(reslicey);
319          gradx = imf_diffx(smooth_im);
320          grady = imf_diffy(smooth_im);
321          im_free(smooth_im);
322          chisq -= im_compare(sliceyx, sliceyy, gradx, grady,
323                              masky, roiy);
324          im_free(gradx);
325          im_free(grady);
326 
327          reslicez = seq_slicez(zcentre, roiz, coreg_bproj);
328          smooth_im = im_prepare(reslicez);
329          im_free(reslicez);
330          gradx = imf_diffx(smooth_im);
331          grady = imf_diffy(smooth_im);
332          im_free(smooth_im);
333          chisq -= im_compare(slicezx, slicezy, gradx, grady,
334                              maskz, roiz);
335          im_free(gradx);
336          im_free(grady);
337      }
338      return(chisq);
339 }
340 
341 void coreg_auto(double threshold, double bscale, double border, int mask)
342 {
343     double *params, *scales, chisq;
344     int npar;
345 
346     smooth = bscale;
347     minmask = mask;
348     if(a_switch==0) store_rot_init();
349     params = dvector_alloc(0,10);
350     npar = coreg_get_vec(params,mask);
351     scales = dvector_alloc(0,10);
352 
353     scales[0] = 1.0;
354     scales[1] = 1.0;
355     scales[2] = 1.0;
356     scales[3] = 0.1;
357     scales[4] = 0.1;
358     scales[5] = 0.1;
359     scales[6] = 0.1;
360     scales[7] = 0.1;
361     scales[8] = 0.1;
362     scales[9] = 0.1;
363 
364     init_voxchisq(threshold,border);
365     chisq = simplexmin2(npar,params,scales,voxchisq,NULL,c_test1,(void(*)())format);
366     chisq = voxchisq(npar,params);
367     dvector_free(params,0);
368     dvector_free(scales,0);
369 
370     if(a_switch==0) store_rot_reset();
371 }
372 
373 
374 void border_detect(double border)
375 {
376         /* Search the volume for identical zero pixels and modify the
377         border to ignore them. It is important to find the first non-zero
378         pixel, not the last zero pixel, in there are zeros within the
379         data block.  The border parameter should be set by the user to
380         deal with any aliasing problems in the data (missing partial
381         slices due to rotation. PAB 16/05/2003 */
382 
383         void **images=NULL;
384         Imrect *im=NULL;
385         Imregion *roi=NULL, proi;
386         int i, j, l;
387         int bly, buy, blx, bux;
388 
389         images = pvector_alloc(0, 3);
390 
391         images[0] = get_xcoreg_im();
392         images[1] = get_ycoreg_im();
393         images[2] = get_zcoreg_im();
394 
395         for(i=0; i<3; i++)
396         {
397                 im = (Imrect *)images[i];
398                 roi = im->region;
399 
400                 blx = roi->lx;
401                 bux = roi->ux;
402                 bly = roi->ly;
403                 buy = roi->uy;
404 
405                 l = (roi->ux-roi->lx)/2;
406                 for(j=roi->ly; j<(roi->uy-roi->ly)/2; j++)
407                 {
408                         if(im_get_pixf(im, j, l)!=0.0) 
409                         {
410                                 bly = j;                        
411                                 break;
412                         }
413                 }
414                 for(j=roi->uy; j>(roi->uy-roi->ly)/2; j--)
415                 {
416                         if(im_get_pixf(im, j, l)!=0.0) 
417                         {
418                                 buy = j;        
419                                 break;
420                         }
421                 }
422 
423                 l = (roi->uy-roi->ly)/2;
424                 for(j=roi->lx; j<(roi->ux-roi->lx)/2; j++)
425                 {
426                         if(im_get_pixf(im, l, j)!=0.0) 
427                         {
428                                 blx = j;
429                                 break;
430                         }
431                 }
432                 for(j=roi->ux; j>(roi->ux-roi->lx)/2; j--)
433                 {
434                         if(im_get_pixf(im, l, j)!=0.0) 
435                         {
436                                 bux = j;
437                                 break;  
438                         }
439                 } 
440 
441                 blx = blx+border;
442                 bux = bux-border;
443                 bly = bly+border;
444                 buy = buy-border;
445 
446                 proi = *(im->region);
447 
448                 proi.lx = blx;
449                 proi.ux = bux;
450                 proi.ly = bly;
451                 proi.uy = buy;
452 
453                 images[i] = im_subim(im, &proi);
454         }
455         Rslicex = images[0];
456         Rslicey = images[1];
457         Rslicez = images[2];
458         pvector_free(images, 0);
459 }
460 
461 
462 void imlimits_finder(float *fmin, float *fmax)
463 {
464         int i, j, k, lz, uz, ly, uy, lx, ux;
465         void ***imptrs=NULL;
466         float max, min, pix;
467 
468         max = -FLT_MAX;
469         min = FLT_MAX;
470 
471         imptrs = seq_limits(&lz, &uz, &ly, &uy, &lx, &ux);
472         for(i=lz; i<uz; i++)
473         {
474                 for(j=ly; j<uy; j++)
475                 {
476                         for(k=lx; k<ux; k++)
477                         {
478                                 pix = this_pixel(imptrs, i, j, k);
479                                 if(pix<min) min = pix;
480                                 if(pix>max) max = pix;
481                         }
482                 }
483         }
484         *fmin = min;
485         *fmax = max;
486 }
487 
488 
489 void imnoise_finder(float *rnoise, float *fnoise)
490 {
491         Imrect *im;
492         Imregion *roi=NULL;
493         double k;
494         int type;
495 
496         im = im_copy(Rslicez);
497         roi = im->region;
498         stack_push(im, IMRECT, im_free);
499         imcalc_noise(&k, roi);
500         im = (Imrect *)stack_pop(&type);
501         im_free(im);
502         *rnoise = k;
503 
504         roi = Rslicez->region;
505         im = seq_slicez(zcentre, roi, coreg_bproj);
506         stack_push(im, IMRECT, im_free);
507         imcalc_noise(&k,roi);
508         im = (Imrect *)stack_pop(&type);
509         im_free(im);
510         *fnoise = k;
511 }
512 
513 
514 void init_mui_voxchisq(double threshold, double border)
515 {
516    float rmin, r1min, r2min, r3min, rmax, r1max, r2max, r3max, rnoise;
517    float fmin, fmax, fnoise, bin_mult=1.0;
518    Vec3 dummy;
519 
520    if (hist2 != NULL)
521    {
522        hfree(hist2);
523        hist2 = NULL;
524    }
525 
526    if(Rslicez!=NULL)
527    {
528       im_free(Rslicez);
529       Rslicez=NULL;
530    }
531    if(Rslicey!=NULL)
532    {
533       im_free(Rslicey);
534       Rslicey=NULL;
535    }
536    if(Rslicex!=NULL)
537    {
538       im_free(Rslicex);
539       Rslicex=NULL;
540    }
541    border_detect(border);
542 
543    if (Rslicex == NULL || Rslicey == NULL || Rslicez == NULL)
544      return;
545 
546    if (Rslicex != NULL) roix = Rslicex->region;
547    else roix = NULL;
548    if (Rslicey != NULL) roiy = Rslicey->region;
549    else roiy = NULL;
550    if (Rslicez != NULL) roiz = Rslicez->region;
551    else roiz = NULL;
552 
553    coreg_centre(&xcentre, &ycentre, &zcentre, &dummy);
554 
555    Tslicez = seq_slicez(zcentre, roiz, coreg_bproj);
556    Tslicey = seq_slicey(ycentre, roiy, coreg_bproj);
557    Tslicex = seq_slicex(xcentre, roix, coreg_bproj);
558 
559    if (Tslicex == NULL || Tslicey == NULL || Tslicez == NULL) return;
560    if (hist2 == NULL)
561    {
562       imf_minmax(Rslicex,&r1min,&r1max);
563       imf_minmax(Rslicey,&r2min,&r2max);
564       imf_minmax(Rslicez,&r3min,&r3max);
565 
566       if (r1min < r2min) rmin = r1min;
567           else rmin = r2min;
568       if (r3min < rmin) rmin = r3min;
569 
570       if (r1max > r2max) rmax = r1max;
571           else rmax = r2max;
572       if (r3max > rmax) rmax = r3max;
573 
574       imlimits_finder(&fmin,&fmax);
575 
576         imnoise_finder(&rnoise, &fnoise);
577         fmin = fmin-(int)(2.0*fnoise);
578         fmax = fmax+(int)(2.0*fnoise);
579         rmin = rmin-(int)(2.0*rnoise);
580         rmax = rmax+(int)(2.0*rnoise);
581 
582         if(bin_size_switch==0)
583         {
584                 bin_mult=1.0;
585         }
586         else
587         {
588                 bin_mult=2.0;
589         }
590 
591         no_rbins = (int)(bin_mult*(rmax-rmin)/rnoise);
592         no_fbins = (int)(bin_mult*(fmax-fmin)/fnoise);
593 
594         if(bin_switch==1)
595         {
596                 if(no_rbins>256) no_rbins = 256;
597                 if(no_rbins<64) no_rbins=64;
598                 if(no_fbins>256) no_fbins = 256;
599                 if(no_fbins<64) no_fbins=64;
600         }
601 
602         hist2 = hbook2(0,"image contents\0",
603                         rmin,rmax,no_rbins,fmin,fmax,no_fbins);
604    }
605 
606    if(Tslicez!=NULL)
607    {
608       im_free(Tslicez);
609       Tslicez=NULL;
610    }
611    if(Tslicey!=NULL)
612    {
613       im_free(Tslicey);
614       Tslicey=NULL;
615    }
616    if(Tslicex!=NULL)
617    {
618       im_free(Tslicex);
619       Tslicex=NULL;
620    }
621 }
622 
623 
624 static double mui_voxchisq_maja(int n_par, double *a)
625 {
626    double chisq = 0.0;
627    if (roix==NULL|| roiy==NULL || roiz == NULL) return(chisq);
628    if(coreg_set_vec(a,minmask))
629    {
630       hreset(hist2);
631 
632       Tslicez = seq_slicez(zcentre, roiz, coreg_bproj);
633       mui_joint_hist2(Rslicez,Tslicez,hist2);
634       im_free(Tslicez);
635 
636       Tslicey = seq_slicey(ycentre, roiy, coreg_bproj);
637       mui_joint_hist2(Rslicey,Tslicey,hist2);
638       im_free(Tslicey);
639 
640       Tslicex = seq_slicex(xcentre, roix, coreg_bproj);
641       mui_joint_hist2(Rslicex,Tslicex,hist2);
642       hist_smooth(hist2);
643       im_free(Tslicex);
644 
645       chisq = (double) mui_IFR(hist2->array, no_fbins, no_rbins);
646    }
647    return(chisq);
648 }
649 
650 
651 double mui_voxchisq_pab(int n_par, double *a)
652 {
653         double chisq = 0.0;
654 
655         if (roix==NULL|| roiy==NULL || roiz == NULL) return(chisq);
656         if(coreg_set_vec(a,minmask))
657         {
658                 Tslicez = seq_slicez(zcentre, roiz, coreg_bproj);
659                 Tslicey = seq_slicey(ycentre, roiy, coreg_bproj);
660                 Tslicex = seq_slicex(xcentre, roix, coreg_bproj);
661 
662                 hreset(hist2);
663                 mui_joint_hist2(Rslicex, Tslicex, hist2);
664                 mui_joint_hist2(Rslicey, Tslicey, hist2);
665                 mui_joint_hist2(Rslicez, Tslicez, hist2);
666 
667                 chisq = (double) mui_IFR_pab(hist2, no_rbins, no_fbins, Rslicex, Tslicex, Rslicey,
668                                                 Tslicey, Rslicez, Tslicez, roix, roiy, roiz);
669 
670                 im_free(Tslicey);
671                 im_free(Tslicez);
672                 im_free(Tslicex);
673         }
674         return(chisq);
675 }
676 
677 
678 static double mui_voxchisq(int n_par, double *a)
679 {
680     double chisq = 0.0;
681     if(mi_switch==0)
682     {
683        chisq = mui_voxchisq_maja(n_par, a);
684     }
685     else
686     {
687        chisq = mui_voxchisq_pab(n_par, a);
688     }
689    return(chisq);
690 }
691 
692 
693 void mui_coreg_auto(double threshold, double bscale, double border, int mask)
694 {
695     double *params;
696     double *scales;
697     double chisq;
698     int npar;
699 
700     smooth = bscale;
701     minmask = mask;
702     if(a_switch==0) store_rot_init();
703     params = dvector_alloc(0,10);
704     npar = coreg_get_vec(params,mask);
705     scales = dvector_alloc(0,10);
706 
707     scales[0] = 1.0;
708     scales[1] = 1.0;
709     scales[2] = 1.0;
710     scales[3] = 0.1;
711     scales[4] = 0.1;
712     scales[5] = 0.1;
713     scales[6] = 0.1;
714     scales[7] = 0.1;
715     scales[8] = 0.1;
716     scales[9] = 0.1;
717 
718     init_mui_voxchisq(threshold,border);
719     chisq = simplexmin2(npar,params,scales,mui_voxchisq,NULL,c_test1,(void(*)())format);
720     chisq = mui_voxchisq(npar,params);
721 
722     dvector_free(params,0);
723     dvector_free(scales,0);
724     if (hist2 != NULL)
725     {
726        hfree(hist2);
727        hist2 = NULL;
728     }
729     if(a_switch==0) store_rot_reset();
730 }
731 
732 
733 void coreg_covar(double threshold,double border, int mask)
734 {
735         double *step_vec=NULL, *params=NULL, **cov_mat=NULL;
736         int npar, i, j, return_flag=0;
737 
738         params = dvector_alloc(0, 10);
739         minmask=mask;
740         init_mui_voxchisq(threshold, border);
741         npar = coreg_get_vec(params, mask);
742 
743         if(step_switch!=2)
744         {
745                 for(i=0; i<npar; i++)
746                 {
747                         if(params[i]==0.0)
748                         {
749                                 format("Bad step size: param %i = 0 \n", i);
750                                 return_flag=1;
751                         }
752                 }
753                 if(return_flag==1)
754                 {
755                         dvector_free(params, 0);
756                         return;
757                 }
758         }
759 
760         if(step_switch==0)              /* For brainweb data */
761         {
762                 step_vec = dvector_alloc(0, 10);
763                 step_vec[0] = params[0]+params[0]*50.0*0.00001;
764                 step_vec[1] = params[1]+params[1]*50.0*0.00001;
765                 step_vec[2] = params[2]+params[2]*50.0*0.00001;
766                 step_vec[3] = params[3]+params[3]*50.0*0.0001;
767                 step_vec[4] = params[4]+params[4]*50.0*0.01;
768                 step_vec[5] = params[5]+params[5]*50.0*0.01;
769                 step_vec[6] = params[6]+params[6]*50.0*0.00001;
770                 step_vec[7] = params[7]+params[7]*50.0*0.00001;
771                 step_vec[8] = params[8]+params[8]*50.0*0.0001;
772         }
773         else if(step_switch==1)     /* For thacker data */
774         {
775                 step_vec = dvector_alloc(0, 10);
776                 step_vec[0] = params[0]+params[0]*50.0*0.001;
777                 step_vec[1] = params[1]+params[1]*50.0*0.001;
778                 step_vec[2] = params[2]+params[2]*50.0*0.001;
779                 step_vec[3] = params[3]+params[3]*50.0*0.001;
780                 step_vec[4] = params[4]+params[4]*50.0*0.001;
781                 step_vec[5] = params[5]+params[5]*50.0*0.001;
782                 step_vec[6] = params[6]+params[6]*50.0*0.001;
783                 step_vec[7] = params[7]+params[7]*50.0*0.001;
784                 step_vec[8] = params[8]+params[8]*50.0*0.001;
785         }
786 
787         cov_mat = coreg_covar_calc(mask,border,hist2,Rslicex, Rslicey, Rslicez, step_vec);
788         for(i=0; i<npar; i++) format("Step[%i]=%f \n", i, step_vec[i]);
789         for (i = 0; i < npar; i++)
790         {
791                 for (j = 0; j < npar;j++)
792                 {
793                         format("Cov[%i][%i]=%f ",i,j,cov_mat[i][j]);
794                 }
795                 format("\n");
796         }
797         darray_free(cov_mat, 0, 0, npar, npar);
798 
799         if (hist2 != NULL)
800         {
801                 hfree(hist2);
802                 hist2 = NULL;
803         }
804         if(step_vec!=NULL)
805         {
806                 dvector_free(step_vec, 0);
807         }
808 }
809 
810 
811 void coreg_covar_bulk(double threshold,double border, int mask, char *output_path)
812 {
813         FILE *output_file=NULL;
814         double *step_vec=NULL, *params=NULL, **cov_mat=NULL;
815         int i, j, p, return_flag=0, npar;
816 
817         params = dvector_alloc(0, 10);
818         minmask=mask;
819         init_mui_voxchisq(threshold, border);
820         npar = coreg_get_vec(params, mask);
821 
822         if(step_switch!=2)
823         {
824                 for(i=0; i<npar; i++)
825                 {
826                         if(params[i]==0.0)
827                         {
828                                 format("Bad step size: param %i = 0 \n", i);
829                                 return_flag=1;
830                         }
831                 }
832                 if(return_flag==1)
833                 {
834                         dvector_free(params, 0);
835                         return;
836                 }
837         }
838         if((output_file = fopen(output_path, "w"))==NULL)
839         {
840                 format("Could not open output file\n");
841                 return;
842         }
843 
844         for(p=0; p<100; p++)
845         {
846                 if(step_switch==0)              /* For brainweb data */
847                 {
848                         step_vec = dvector_alloc(0, 10);
849                         step_vec[0] = params[0]-params[0]*(p-50)*0.00001;
850                         step_vec[1] = params[1]-params[1]*(p-50)*0.00001;
851                         step_vec[2] = params[2]-params[2]*(p-50)*0.00001;
852                         step_vec[3] = params[3]-params[3]*(p-50)*0.0001;
853                         step_vec[4] = params[4]-params[4]*(p-50)*0.01;
854                         step_vec[5] = params[5]-params[5]*(p-50)*0.01;
855                         step_vec[6] = params[6]-params[6]*(p-50)*0.00001;
856                         step_vec[7] = params[7]-params[7]*(p-50)*0.00001;
857                         step_vec[8] = params[8]-params[8]*(p-50)*0.0001;
858                 }
859                 else if(step_switch==1)     /* For thacker data */
860                 {
861                         step_vec = dvector_alloc(0, 10);
862                         step_vec[0] = params[0]-params[0]*(p-50)*0.001;
863                         step_vec[1] = params[1]-params[1]*(p-50)*0.001;
864                         step_vec[2] = params[2]-params[2]*(p-50)*0.001;
865                         step_vec[3] = params[3]-params[3]*(p-50)*0.001;
866                         step_vec[4] = params[4]-params[4]*(p-50)*0.001;
867                         step_vec[5] = params[5]-params[5]*(p-50)*0.001;
868                         step_vec[6] = params[6]-params[6]*(p-50)*0.001;
869                         step_vec[7] = params[7]-params[7]*(p-50)*0.001;
870                         step_vec[8] = params[8]-params[8]*(p-50)*0.001;
871                 }
872                 for(i=0; i<npar; i++) fprintf(output_file, "%f \n", step_vec[i]);
873                 cov_mat = coreg_covar_calc(mask,border,hist2,Rslicex, Rslicey, Rslicez, step_vec);
874                 for (i=0; i<npar; i++)
875                 {
876                         for (j=0; j<npar; j++) fprintf(output_file,"%f ",cov_mat[i][j]);
877                         fprintf(output_file,"\n");
878                 }
879                 fprintf(output_file," \n");
880                 darray_free(cov_mat, 0, 0, npar, npar);
881                 cov_mat=NULL;
882         }
883 
884         if (hist2 != NULL)
885         {
886                 hfree(hist2);
887                 hist2 = NULL;
888         }
889         if(step_vec!=NULL)
890         {
891                 dvector_free(step_vec, 0);
892         }
893 }
894 

~ [ 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.