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

Linux Cross Reference
Tina4/src/tools/coreg/coreg_auto.c

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

  1 #include <stdio.h>
  2 #include <tina/sys.h>
  3 #include <tina/sysfuncs.h>
  4 #include <tina/math.h>
  5 #include <tina/mathfuncs.h>
  6 #include <tina/vision.h>
  7 #include <tina/visionfuncs.h>
  8 
  9 /* global variables set up by init voxchisq */
 10 static Imrect *slicexx=NULL, *sliceyx=NULL, *slicezx=NULL;
 11 static Imrect *slicexy=NULL, *sliceyy=NULL, *slicezy=NULL;
 12 static Imrect *maskx=NULL, *masky=NULL, *maskz=NULL;
 13 static Imregion *roix=NULL, *roiy=NULL, *roiz=NULL;
 14 static double xcentre, ycentre, zcentre;
 15 static double c_test1 = 0.00001, c_test2 = 0.1;
 16 static double perfect;
 17 static double smooth;
 18 static int minmask;
 19 
 20 /* access functions to coregistration images */
 21 extern Imrect *get_zcoreg_im();
 22 extern Imrect *get_ycoreg_im();
 23 extern Imrect *get_xcoreg_im();
 24 extern Imrect *coreg_slicez(float zslice, Imregion *within);
 25 extern Imrect *coreg_slicey(float yslice, Imregion *within);
 26 extern Imrect *coreg_slicex(float xslice, Imregion *within);
 27 extern void coreg_centre(double *x, double *y, double *z, Vec3*c);
 28 extern int coreg_get_vec(double *coreg_vec, int mask);
 29 extern int coreg_set_vec(double *coreg_vec, int mask);
 30 extern void store_rot_reset(), store_rot_init();
 31 
 32 static Imrect *im_prepare(Imrect *im)
 33 {
 34      Imrect *im1, *im2;
 35      im1 = imf_lsf_smooth(im,smooth);
 36      im2 = imf_lsf_smooth(im,smooth);
 37      im_free(im1);
 38      return(im2);
 39 }
 40 
 41 static double  im_compare(Imrect *im1, Imrect *im2, Imrect *im3,
 42                           Imrect *im4, Imrect *mask, Imregion *roi)
 43 {
 44     int lx, ux, ly, uy;
 45     int i,j;
 46     float *row1, *row2, *row3 , *row4, *row5 ;
 47     double dot_prod=0.0;
 48     double norm=0.0;
 49 
 50     lx = roi->lx;
 51     ux = roi->ux;
 52     ly = roi->ly;
 53     uy = roi->uy;
 54 
 55     for (i = ly; i < uy; ++i)
 56     {
 57         row1 = (float*) IM_ROW(im1,i);
 58         row2 = (float*) IM_ROW(im2,i);
 59         row3 = (float*) IM_ROW(im3,i);
 60         row4 = (float*) IM_ROW(im4,i);
 61         row5 = (float*) IM_ROW(mask,i);
 62         row1 = &row1[lx];
 63         row2 = &row2[lx];
 64         row3 = &row3[lx];
 65         row4 = &row4[lx];
 66         row5 = &row5[lx];
 67         for (j = lx; j < ux; ++j)
 68         {
 69             if (*(row5++)>0.5)
 70             {
 71                norm += (*row3 * *row3) + (*row4 * *row4);
 72                dot_prod +=  *(row1++)* *(row3++) +
 73                             *(row2++)* *(row4++);
 74 /*
 75                dot_prod += sqrt((*(row1)**(row1++)+
 76                                  *(row2)**(row2++))
 77                                *(*(row3)**(row3++)+
 78                                  *(row4)**(row4++)) );
 79 */
 80             }
 81             else
 82             {
 83                row1++;
 84                row2++;
 85                row3++;
 86                row4++;
 87             }
 88         }
 89     }
 90 
 91     return( dot_prod/sqrt(norm) );
 92 }
 93 
 94 Imrect *crop_region(Imrect *im, double border)
 95 {
 96    Imrect *im1;
 97    Imregion roi;
 98 
 99    if(im!=NULL )  roi = *(im->region);
100    else  return(NULL);
101    roi.lx += (int)border;
102    roi.ux -= (int)border;
103    roi.ly += (int)border;
104    roi.uy -= (int)border;
105    im1 = im_subim(im, &roi);
106    
107    return(im1);
108 }
109 
110 void init_voxchisq(double threshold, double border)
111 {
112     Imrect *im0,*im1,*smooth_im;
113     Vec3 dummy;
114 
115     if (maskx!=NULL) im_free(maskx);
116     if (masky!=NULL) im_free(masky);
117     if (maskz!=NULL) im_free(maskz);
118     if (slicexx!=NULL) im_free(slicexx);
119     if (sliceyx!=NULL) im_free(sliceyx);
120     if (slicezx!=NULL) im_free(slicezx);
121     if (slicexy!=NULL) im_free(slicexy);
122     if (sliceyy!=NULL) im_free(sliceyy);
123     if (slicezy!=NULL) im_free(slicezy);
124     maskx = masky = maskz = NULL;
125     slicexx = sliceyx = slicezx = NULL;
126     slicexy = sliceyy = slicezy = NULL;
127 
128     coreg_centre(&xcentre, &ycentre, &zcentre, &dummy);
129     im0 = crop_region(get_xcoreg_im(),border); 
130     im1 = imf_mod(im0);
131     maskx = imf_bthresh(threshold,im1);
132     im_free(im1);
133     smooth_im = im_prepare(im0);
134     im_free(im0);
135     slicexx = imf_diffx(smooth_im);
136     slicexy = imf_diffy(smooth_im);
137     im_free(smooth_im);
138 
139     im0 = crop_region(get_ycoreg_im(),border);
140     im1 = imf_mod(im0);
141     masky = imf_bthresh(threshold,im1);
142     im_free(im1);
143     smooth_im = im_prepare(im0);
144     im_free(im0);
145     sliceyx = imf_diffx(smooth_im);
146     sliceyy = imf_diffy(smooth_im);
147     im_free(smooth_im);
148 
149     im0 = crop_region(get_zcoreg_im(),border);
150     im1 = imf_mod(im0);
151     maskz = imf_bthresh(threshold,im1);
152     im_free(im1);
153     smooth_im = im_prepare(im0);
154     im_free(im0);
155     slicezx = imf_diffx(smooth_im);
156     slicezy = imf_diffy(smooth_im);
157     im_free(smooth_im);
158 
159     if (slicexx!=NULL) roix = slicexx->region;
160     else roix = NULL;
161     if (sliceyx!=NULL) roiy = sliceyx->region;
162     else roiy = NULL;
163     if (slicezx!=NULL) roiz = slicezx->region;
164     else roiz =NULL;
165     perfect = im_compare(slicexx,slicexy,slicexx,slicexy,maskx,roix);
166     perfect += im_compare(sliceyx,sliceyy,sliceyx,sliceyy,masky,roiy);
167     perfect += im_compare(slicezx,slicezy,slicezx,slicezy,maskz,roiz);
168 }
169 
170 static double voxchisq(int n_par, double *a)
171 {
172      double chisq = MAXDOUBLE;
173      Imrect *reslicex, *reslicey, *reslicez;
174      Imrect *smooth_im;
175      Imrect *gradx, *grady, *gradz;
176 
177      if (roix==NULL|| roiy==NULL || roiz == NULL) return(chisq);
178      if(coreg_set_vec(a,minmask)) 
179      {
180          chisq = perfect;
181          reslicex = coreg_slicex(xcentre, roix);
182          smooth_im = im_prepare(reslicex);
183          im_free(reslicex);
184          gradx = imf_diffx(smooth_im);
185          grady = imf_diffy(smooth_im);
186          im_free(smooth_im);
187          chisq -= im_compare(slicexx, slicexy, gradx, grady,
188                              maskx, roix);
189          im_free(gradx);
190          im_free(grady);
191 
192          reslicey = coreg_slicey(ycentre, roiy);
193          smooth_im = im_prepare(reslicey);
194          im_free(reslicey);
195          gradx = imf_diffx(smooth_im);
196          grady = imf_diffy(smooth_im);
197          im_free(smooth_im);
198          chisq -= im_compare(sliceyx, sliceyy, gradx, grady,
199                              masky, roiy);
200          im_free(gradx);
201          im_free(grady);
202 
203          reslicez = coreg_slicez(zcentre, roiz);
204          smooth_im = im_prepare(reslicez);
205          im_free(reslicez);
206          gradx = imf_diffx(smooth_im);
207          grady = imf_diffy(smooth_im);
208          im_free(smooth_im);
209          chisq -= im_compare(slicezx, slicezy, gradx, grady,
210                              maskz, roiz);
211          im_free(gradx);
212          im_free(grady);
213      }
214      return(chisq);
215 }
216 
217 coreg_auto(double threshold, double bscale, double border, int mask)
218 {
219     double *params;
220     double *scales;
221     double chisq;
222     int npar;
223 
224     smooth = bscale; 
225     minmask = mask;
226     params = dvector_alloc(0,9);
227     store_rot_init();
228     npar = coreg_get_vec(params,mask);
229     scales = dvector_alloc(0,9);
230     scales[0] = 1.0;
231     scales[1] = 1.0;
232     scales[2] = 1.0;
233     scales[3] = 0.02;
234     scales[4] = 0.02;
235     scales[5] = 0.02;
236     scales[6] = 0.02;
237     scales[7] = 0.02;
238     scales[8] = 0.02;
239     init_voxchisq(threshold,border);
240     chisq = simplexmin2(npar,params,scales,voxchisq,c_test1,(void(*)())format); 
241     chisq = voxchisq(npar,params);
242     dvector_free(params,0);
243     dvector_free(scales,0);
244     store_rot_reset();
245 }
246 

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