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

Linux Cross Reference
Tina4/src/tools/imcalc/norm.c

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

  1 #include <stdio.h>
  2 #include <sys/param.h>
  3 #include <string.h>
  4 #include <tina/sys.h>
  5 #include <tina/sysfuncs.h>
  6 #include <tina/math.h>
  7 #include <tina/mathfuncs.h>
  8 #include <tina/vision.h>
  9 #include <tina/visionfuncs.h>
 10 #include <tina/file.h>
 11 #include <tina/filefuncs.h>
 12 #include <tina/file_gen.h>
 13 #include <tina/tv.h>
 14 #include <tina/toolsfuncs.h>
 15 #include <tina/hist_funcs.h>
 16 
 17 #define STAT_LIM 2.5 
 18 #define INF_LIM 200.0 
 19 #define CORSCAN 26
 20 /* Modifications */
 21 /* Use of a single limit (STAT_LIM) on all statistical NULL hypotheses.
 22 /* A separate limit for low information INF_LIM.
 23 /* Normalisation software modified to improve matching of statistical
 24    tests to algorithm stability and eliminate edge threshold.
 25    Regularisation for indeterminate regions modified to reduce bias
 26    towards null gain change.
 27    Second differential integration routine addded to reduce 
 28    reconstruction errors NAT 25/5/2001  */
 29 /* slope estimation routine enf_integ() added to enforce integrability.
 30    routine added to maximise information content im_corscale(), based
 31    upon Bhattacharrya overlap with scale invariant histogram distribution.
 32    NAT 25/12/2001 */
 33 
 34 
 35 /* mask of intensities < 3 sigma */
 36 static Imrect * thresh_slice(Imrect * im, double sig_noise)
 37 {
 38     Imrect *im1,*im2,*im3,*imcut;
 39 
 40     im1 = im_thresh((STAT_LIM*(sig_noise)),im);
 41     im2 = im_thresh((-STAT_LIM*(sig_noise)),im);
 42     im3 = im_diff(im2,im1);
 43     im_free(im1); im_free(im2);
 44     /* fudge for inversion recovery borders */
 45     im1 = im_thresh(4096.0,im);
 46     imcut = im_sum(im1,im3);
 47     im_free(im1); im_free(im3);
 48 
 49     return (imcut);
 50 }
 51 
 52 /* mask of edges > STAT_LIM  */
 53 static Imrect *edge_slice(Imrect *image, double sig_noise) 
 54 {
 55     Imrect *im,*im2,*im3,*im4;
 56     Imrect *imcut;
 57     
 58     im = imf_gauss(image, 1.0, 0.01);
 59 
 60     im2 = imf_diffx(im);
 61     im3 = im_sqr(im2);
 62     im_free(im2); 
 63    
 64     im2 = imf_diffy(im);
 65     im4 = im_sqr(im2);
 66     im_free(im); im_free(im2);
 67 
 68     im = im_sum(im3,im4);
 69     im2 = im_sqrt(im);
 70     im_free(im); im_free(im3); im_free(im4);
 71 
 72     im = im_thresh((STAT_LIM*sig_noise/sqrt(2)),im2);
 73     im_free(im2);
 74 
 75     return (im);
 76 }
 77 
 78 static Imrect *zero_cut(double thresh, Imrect * cut) 
 79 {
 80     Imrect *im,*im1,*imcut;
 81 
 82     im1 = im_mod(cut);
 83     im = im_bthresh(thresh, im1);
 84     imcut = im_cast(im, uchar_v); 
 85     im_free(im); im_free(im1); 
 86 
 87     return (imcut);
 88 }
 89 
 90 static Imrect *zero_mask(Imrect * image)  /* mask of voxels = 0 */
 91 {
 92     Imrect *im_new,*im1,*im2;
 93     
 94     im1 = zero_cut(MINFLOAT,image);
 95     im2 = im_minus(im1);
 96     im_new = im_add(1.0,im2);    
 97     im_free(im1);  im_free(im2);
 98     
 99     return (im_new);
100     
101 }
102 
103 static Imrect * edge_remove(Imrect * imcut, Imrect * image)
104 {
105     Imrect *im,*new_image;
106 
107     im = im_prod(image,imcut);
108     new_image = im_diff(image,im);
109     im_free(im);
110     return (new_image);
111 }
112 
113 /*                                    */
114 /* mask combining two slices of edges */
115 /*                                    */
116 void edge_mask(double sig_noise, Imrect *im1, Imrect *im2, Imrect **cut,
117                 Imrect **image1, Imrect **image2)
118 {
119     Imrect *cut1,*cut2,*cut3;
120     Imrect *cut_temp1,*cut_temp2;
121     Imrect *im1_temp,*im2_temp;
122 
123     im1_temp = imf_median(im1);
124     im2_temp = imf_median(im2);
125     cut1 = thresh_slice(im1_temp, sig_noise);
126     cut2 = edge_slice(im1_temp, sig_noise);
127     cut3 = im_sum(cut1,cut2);
128     im_free(cut1);  im_free(cut2);
129     cut1 = zero_mask(im1_temp);
130     cut_temp1 = im_sum(cut1,cut3);
131     im_free(cut1);  im_free(cut3);
132     /* and again... */
133     cut1 = thresh_slice(im2_temp, sig_noise);
134     cut2 = edge_slice(im2_temp, sig_noise);
135     cut3 = im_sum(cut1,cut2);
136     im_free(cut1);  im_free(cut2);
137     cut1 = zero_mask(im2_temp);
138     cut_temp2 = im_sum(cut1,cut3);
139     im_free(cut1);  im_free(cut3);
140 
141     cut1 = im_sum(cut_temp1,cut_temp2);
142     im_free(cut_temp1);    im_free(cut_temp2);
143     cut2 = zero_cut(MINFLOAT,cut1);
144     im_free(cut1);
145 
146     *image1 = edge_remove(cut2,im1_temp);
147     *image2 = edge_remove(cut2,im2_temp);
148     *cut = cut2;
149     im_free(im1_temp); im_free(im2_temp);
150     return;
151 }
152 
153 static shistogram *imf_hist_ratio(double sig_noise, double range,
154                            Imrect *ima, Imrect* imb, Imrect *immask,
155                            Imregion *roi,int nbin)
156 {
157     shistogram *imhist;
158     Imrect *im2;
159     float *rowa, *rowb, *rowmask;
160     int lx, ux, ly, uy;
161     int width, height;
162     int i, j, ndata;
163     float imin,imax;
164     float shift = range/(float)nbin;
165     double ratio, log_ratio, errweight;
166 
167     if (ima == NULL || imb == NULL)
168         return(NULL);
169     if (roi == NULL)
170         return(NULL);
171     lx = roi->lx;
172     ux = roi->ux;
173     ly = roi->ly;
174     uy = roi->uy;
175 
176     width = ima->width;
177     height = ima->height;
178 
179     imhist = hbook1(0,"image contents\0",(float)(-range),
180                 (float)(range), nbin);
181 
182     rowa = fvector_alloc(lx, ux);
183     rowb = fvector_alloc(lx, ux);
184     rowmask = fvector_alloc(lx, ux);
185 
186     for (i = ly; i < uy; ++i)
187     {
188         im_get_rowf(rowa, ima, i, lx, ux);
189         im_get_rowf(rowb, imb, i, lx, ux);
190         im_get_rowf(rowmask, immask, i, lx, ux);
191         for (j = lx; j < ux; ++j)
192         {
193             if (rowb[j]==0.0) ratio = MAXFLOAT;
194             else ratio = rowa[j]/rowb[j];
195             if (ratio>0.0) log_ratio = log(ratio); 
196             else log_ratio = -MAXFLOAT;
197 /* weight computed based one error on log-ratio NAT 20/1/2001 */ 
198             errweight = rowa[j]*rowa[j]*rowb[j]*rowb[j]
199                       /(rowa[j]*rowa[j]+rowb[j]*rowb[j]);
200             errweight /= (sig_noise*sig_noise);
201  
202             if(rowmask[j]!=0.0)
203             {
204                hfill1s(imhist,log_ratio-shift,errweight);
205                hfill1s(imhist,log_ratio+shift,errweight);
206             }
207         }
208     }
209  
210     fvector_free(rowa, lx);
211     fvector_free(rowb, lx);
212     fvector_free(rowmask, lx);
213     return(imhist);
214 }
215 
216 double im_fraction(Imrect *cut, Imrect *image1, Imrect *image2, int nbin, 
217                 double sig_noise) 
218 {
219     Imrect *immask;
220     double normal;
221     shistogram *hist=NULL;
222     Imregion *roi=NULL;
223     float maxx=0.0,err,nerr;
224     void graph_hfit(Tv * tv, shistogram * hist);
225 
226     immask   = zero_mask(cut);
227     roi = roi_inter(image1->region,image2->region);
228     hist = imf_hist_ratio(sig_noise,1.0,image1,image2,immask,roi,nbin);
229     graph_hfit(imcalc_graph_tv_get(), hist);
230 
231     normal = (double)hnearmax1(hist,&maxx);
232     if (maxx>0.0 || maxx <=0.0)
233        normal = (double)maxx;
234     else
235     {
236        normal = 0.0;
237        format("bad ratio calculated in fraction()");
238     }
239     im_free(immask);
240     rfree(roi);
241     hfree(hist);
242     return(normal);
243 }
244 
245 static Imrect *err_map(Imrect *im, Imrect *imx, Imrect *imy, double sig_noise,
246                        double range, Imrect **imerrx, Imrect **imerry)
247 {
248     Imrect *imrx,*imry;
249     Imregion *roi;
250     float *rowi,*rowx,*rowy,*rowfinx,*rowfiny;
251     int i,j,lx,ux,ly,uy;
252     double weight;
253     double thresh=sig_noise*range*sqrt(2);
254  
255     roi = im->region;
256  
257     lx = roi->lx;
258     ux = roi->ux;
259     ly = roi->ly;
260     uy = roi->uy;
261  
262     if (lx == ux || ly == uy)
263         return (NULL);
264 
265     imrx = im_alloc(im->height, im->width, roi, float_v);
266     imry = im_alloc(im->height, im->width, roi, float_v);
267     rowi = fvector_alloc(lx, ux);
268     rowx = fvector_alloc(lx, ux);
269     rowy = fvector_alloc(lx, ux);
270     rowfinx = fvector_alloc(lx, ux);
271     rowfiny = fvector_alloc(lx, ux);
272  
273     for (i = ly+1; i < uy-1; ++i)
274     {
275         im_get_rowf(rowi, im, i, lx, ux);
276         im_get_rowf(rowx, imx, i, lx, ux);
277         im_get_rowf(rowy, imy, i, lx, ux);
278         for (j = lx+1 ; j < ux-1; ++j)
279         {
280            if ((fabs(rowi[j]) <= 0.01*sig_noise) ||
281                (fabs(rowx[j]) <= 0.01*sig_noise) ||
282                (fabs(rowy[j]) <= 0.01*sig_noise) || 
283                (fabs(rowi[j]-rowx[j]) > thresh) ||
284                (fabs(rowi[j]-rowy[j]) > thresh) )
285            { 
286               rowfinx[j] = 0.0;
287               rowfiny[j] = 0.0;
288            }
289            else 
290            {
291               rowfinx[j] = (16*sig_noise*sig_noise
292                             *(rowi[j]*rowi[j]+ rowx[j]*rowx[j]))
293                             /(pow((rowi[j]+rowx[j]),4));
294               rowfiny[j] = (16*sig_noise*sig_noise
295                             *(rowi[j]*rowi[j] + rowy[j]*rowy[j]))
296                             /(pow((rowi[j]+rowy[j]),4));
297 /* weighting to eliminate any edges */
298               weight = (rowi[j]-rowx[j])*(rowi[j]-rowx[j])
299                       +(rowi[j]-rowy[j])*(rowi[j]-rowy[j]);
300               weight = 1.0 - (sqrt(weight)/thresh);
301               if (weight <= 0.0) weight = 0.0;
302               weight = sqrt(weight);
303               rowfinx[j] = weight/rowfinx[j];
304               rowfiny[j] = weight/rowfiny[j];
305             }
306         }
307         im_put_rowf(rowfinx, imrx, i, lx, ux);
308         im_put_rowf(rowfiny, imry, i, lx, ux);
309     }
310                 
311     fvector_free((void *) rowi, lx);
312     fvector_free((void *) rowx, lx);
313     fvector_free((void *) rowy, lx);
314     fvector_free((void *) rowfinx, lx);
315     fvector_free((void *) rowfiny, lx);
316                 
317     *imerrx = imrx;
318     *imerry = imry;
319 
320     return;
321 }
322 
323 static Imrect* lsf_iter(Imrect *im,int iter, double sigma)
324 {
325     Imrect *imsm,*im1;
326     int     i;
327     
328     im1 = im_copy(im);
329     for(i = 1; i<= iter; ++i)
330     {
331        if (i>1)
332        {
333           im1 = imsm;
334        }
335        imsm = imf_lsf_smooth(im1,sigma);
336        im_free(im1);
337     }
338     return(imsm);
339 }
340 
341 void smooth_slopes(double sig_noise, double nsmear, Imrect *image,
342                    Imrect **imdx, Imrect **imdy, Imrect **corr, 
343                    Imrect **imxs, Imrect **imys)
344 /* routine to calculate smooth slopes imdx and imdy from image  */
345 /* using noise estimate sig_noise and regional smoothing nsmear */
346 {
347     Imrect *im3, *im4, *imdxdy, *imdydx;
348     Imrect *imfx, *imfy;
349     Imrect *im,*im0,*im1,*im2;
350     Imrect *imdx1,*imdy1,*imdx2,*imdy2;
351     Imrect *imx,*imy,*imx_temp,*imy_temp;
352     Imrect *imfx1,*imfy1;
353     Imrect *imdxs, *imdys;
354     int niter = 2;
355 
356 /* smooth image to remove outliers NAT 20/1/2001 */
357     im = im_median(image);
358     im0 = im_median(im);
359     sig_noise *= 0.5;
360     im_free(im);
361     im1 = im_bshift(im0,0,1);
362     im2 = im_bshift(im0,1,0);
363 
364     imx_temp = im_sum(im1,im0);
365     imy_temp = im_sum(im2,im0);
366     imx = im_times(0.5,imx_temp);
367     imy = im_times(0.5,imy_temp);
368     im_free(imx_temp);im_free(imy_temp);
369     imdx1 = im_diff(im0,im1);
370     imdy1 = im_diff(im0,im2);
371     err_map(im0,im1,im2,sig_noise,STAT_LIM,&imfx,&imfy);
372     im_free(im0); im_free(im1); im_free(im2);
373 
374     /* scale slopes */
375     imdx2 = imf_div(imdx1,imx,STAT_LIM*sig_noise,STAT_LIM*sig_noise);
376     imdy2 = imf_div(imdy1,imy,STAT_LIM*sig_noise,STAT_LIM*sig_noise);
377     im_free(imdy1);im_free(imdx1); im_free(imx);im_free(imy);
378     *imdx = imf_prod(imdx2,imfx);
379     *imdy = imf_prod(imdy2,imfy);
380 /*
381     stack_push((void*) imfx,IMRECT,im_free);
382     stack_push((void*) imfy,IMRECT,im_free);
383 */
384     im_free(imdy2);im_free(imdx2);
385 
386     *imxs = lsf_iter(imfx,niter,nsmear);
387     *imys = lsf_iter(imfy,niter,nsmear);
388     im_free(imfx); im_free(imfy);
389     imdxs = lsf_iter(*imdx,niter,nsmear);
390     imdys = lsf_iter(*imdy,niter,nsmear);
391     im_free(*imdx); im_free(*imdy);
392 /*  restrict statistical weighting to bias towards zero image slope
393     when there is no accurate information NAT 25/4/2001 */
394     *imdx = imf_div(imdxs,*imxs,INF_LIM/nsmear,INF_LIM/nsmear);
395     *imdy = imf_div(imdys,*imys,INF_LIM/nsmear,INF_LIM/nsmear);
396     im_free(imdxs); im_free(imdys);
397     im3 = im_bshift(*imdy,0,1);
398     im4 = im_bshift(*imdx,1,0);
399     imdydx = im_diff(*imdy, im3);
400     imdxdy = im_diff(*imdx, im4);
401     im_free(im3); im_free(im4);
402     *corr = im_diff(imdxdy,imdydx); 
403     im_free(imdxdy); im_free(imdydx);
404 }
405 
406 void enf_integ(Imrect **imdx, Imrect **imdy, Imrect *corr, Imrect *imxs, Imrect *imys)
407 {
408 /*  enforce integrability  by averaging of error on cross derivatives
409                                                           NAT 25/12/2001 */
410     int lx,ly,ux,uy,i,j;
411     Imregion roi;
412     double corrpix, pixval1, pixval2, integx, meanx, integy, meany;
413     double weight, norm;
414     roi = *(Imregion *)(corr->region);
415     lx = roi.lx;
416     ux = roi.ux;
417     ly = roi.ly;
418     uy = roi.uy;
419     for (i=lx;i<ux;i++)
420     {
421       IM_PIX_SET(corr, ly, i, 0.0);
422     }
423 
424     for (j=ly;j<uy;j++)
425     {
426       IM_PIX_SET(corr, j, lx, 0.0);
427     }
428 
429     for (i=lx;i<ux;i++)
430     {
431        integx =0.0;
432        meanx = 0.0;
433        norm = 0.0;
434        for (j=ly; j<uy; j++)
435        {
436           IM_PIX_GET(corr, j, i, pixval1);
437           IM_PIX_GET(imxs, j, i, weight);
438           norm += weight;
439           integx += pixval1;
440           meanx += integx*weight;
441        }
442        meanx /= norm;
443        integx = 0.0;
444        for (j=ly; j<uy; j++)
445        {
446           IM_PIX_GET(corr, j, i, pixval1);
447           IM_PIX_GET(*imdx, j, i, pixval2);
448           integx  += pixval1;
449           corrpix = pixval2 - 0.5*(integx - meanx);
450           IM_PIX_SET(*imdx, j, i, corrpix);
451        }
452     }
453     for (j=ly; j<uy; j++)
454     {
455        integy =0.0;
456        meany = 0.0;
457        norm = 0.0;
458        for (i=lx;i<ux;i++)
459        {
460           IM_PIX_GET(corr, j, i, pixval1);
461           IM_PIX_GET(imys, j, i, weight);
462           norm += weight;
463           integy += pixval1;
464           meany += integy*weight;
465        }
466        meany /= norm;
467        integy = 0.0;
468        for (i=lx; i<ux; i++)
469        {
470           IM_PIX_GET(corr, j, i, pixval1);
471           IM_PIX_GET(*imdy, j, i, pixval2);
472           integy  += pixval1;
473           corrpix = pixval2 + 0.5*(integy - meany);
474           IM_PIX_SET(*imdy, j, i, corrpix);
475        }
476     }
477 }
478 
479 static fill_loop(int step, int midx, int midy, float **arrX, float **arrY,
480                  float **arrN, int correct)
481 /* reconstruct integral image loop with integrability enforced by linear
482    constraint if not already guaranteed by use of end_integ() NAT 25/12/25 */
483 {
484     int i,j;
485     int linex, liney;
486     float corner1, corner2, corner3, corner4;
487     float biasweight, bias, cumweight;
488 
489     linex = midx-step;
490     liney = midy-step;
491     corner1 = arrN[liney+1][linex+1] + 0.5 *
492             (-arrY[liney+1][linex+1]
493              -arrX[liney+1][linex+1]
494              -arrY[liney+1][linex]
495              -arrX[liney][linex+1]);
496     linex = midx+step;
497     liney = midy-step;
498     corner2 = arrN[liney+1][linex-1] + 0.5 *
499             (-arrY[liney+1][linex-1]
500              +arrX[liney+1][linex]
501              -arrY[liney+1][linex]
502              +arrX[liney][linex]);
503     linex = midx+step;
504     liney = midy+step;
505     corner3 = arrN[liney-1][linex-1] + 0.5 *
506             (+arrY[liney][linex-1]
507              +arrX[liney-1][linex]
508              +arrY[liney][linex]
509              +arrX[liney][linex]);
510     linex = midx-step;
511     liney = midy+step;
512     corner4 = arrN[liney-1][linex+1] + 0.5 *
513             (+arrY[liney][linex+1]
514              -arrX[liney-1][linex+1]
515              +arrY[liney][linex]
516              -arrX[liney][linex+1]);
517     biasweight = 2*step+1;
518 
519     liney = midy-step;
520     arrN[liney][midx-step] = corner1;
521     for (j=midx-step+1;j<=midx+step;j++)
522     {
523         arrN[liney][j] = arrN[liney][j-1] + arrX[liney][j];
524     }
525     if (correct)
526     {
527         bias = corner2 - arrN[liney][midx+step];
528         cumweight = 0.0;
529         for (j=midx-step+1;j<=midx+step;j++)
530         {
531             cumweight += 1.0;
532             arrN[liney][j] += bias*cumweight/biasweight;
533         }
534         for (j=midx-step+1;j<midx+step;j++)
535         {
536             arrN[liney][j] = (arrN[liney][j] 
537                             + (arrN[liney+1][j] - arrY[liney+1][j]))
538                             /2.0;
539         }
540     }
541 
542     linex = midx+step;
543     arrN[midy-step][linex] = corner2;
544     for (i=midy-step+1;i<=midy+step;i++)
545     {
546         arrN[i][linex] = arrN[i-1][linex] + arrY[i][linex];
547     }
548     if (correct)
549     {
550         bias = corner3 - arrN[midy+step][linex];
551         cumweight = 0.0;
552         for (i=midy-step+1;i<=midy+step;i++)
553         {
554             cumweight += 1.0;
555             arrN[i][linex] += bias*cumweight/biasweight;
556         }
557         for (i=midy-step+1;i<midy+step;i++)
558         {
559             arrN[i][linex] = (arrN[i][linex] 
560                             + (arrN[i][linex-1] + arrX[i][linex]))
561                             /2.0;
562         }
563     }
564 
565     liney = midy+step;
566     arrN[liney][midx+step] = corner3;
567     for (j=midx+step-1;j>=midx-step;j--)
568     {
569         arrN[liney][j] = arrN[liney][j+1] - arrX[liney][j+1];
570     }
571     if (correct)
572     {
573         bias = corner4 - arrN[liney][midx-step];
574         cumweight = 0.0;
575         for (j=midx+step-1;j>=midx-step;j--)
576         {
577             cumweight += 1.0;
578             arrN[liney][j] += bias*cumweight/biasweight;
579         }
580         for (j=midx+step-1;j>midx-step;j--)
581         {
582             arrN[liney][j] = (arrN[liney][j] 
583                             + (arrN[liney-1][j] + arrY[liney][j]))
584                             /2.0;
585         }
586     }
587 
588     linex = midx-step;
589     arrN[midy+step][linex] = corner4;
590     for (i=midy+step-1;i>=midy-step;i--)
591     {
592         arrN[i][linex] = arrN[i+1][linex] - arrY[i+1][linex];
593     }
594     if (correct)
595     {
596         bias = corner1 - arrN[midy-step][linex];
597         cumweight = 0.0;
598         for (i=midy+step-1;i>=midy-step;i--)
599         {
600             cumweight += 1.0;
601             arrN[i][linex] += bias*cumweight/biasweight;
602         }
603         for (i=midy+step-1;i>midy-step;i--)
604         {
605             arrN[i][linex] = (arrN[i][linex]
606                               + (arrN[i][linex+1] - arrX[i][linex+1]))
607                               /2.0;
608         }
609     }
610 }
611 
612 static Imrect *integrate2(Imregion *roi, Imrect *imdx, Imrect *imdy)
613 /* integrate from two sets of derivatives constraining diagonals */
614 {
615     Imrect *im_new;
616     float  **arrN,**arrY,**arrX;
617     int     lx,ly,ux,uy;
618     int     sqmin,sqmax;
619     int     i, j, midx, midy, step;
620 
621     lx = roi->lx;
622     ux = roi->ux;
623     ly = roi->ly;
624     uy = roi->uy;
625     midx = (ux + lx)/2;
626     midy = (uy + ly)/2;
627     sqmin = imin((ux - lx)/2,(uy - ly)/2);
628     sqmax = imax((ux - lx)/2,(uy - ly)/2);
629 
630     im_new = im_alloc(imdx->height, imdy->width, roi, float_v);
631     arrN = farray_alloc(ly,lx,uy,ux);
632     arrY = farray_alloc(ly,lx,uy,ux);
633     arrX = farray_alloc(ly,lx,uy,ux);
634     for (i = ly; i < uy; ++i)
635     {
636         im_get_rowf(arrX[i], imdx, i, lx, ux);
637         im_get_rowf(arrY[i], imdy, i, lx, ux);
638         for (j = lx; j< ux; ++j) arrN[i][j] = 0.0;
639     }
640 
641 /* centre */
642     arrN[midy][midx] = 0.0;
643 /* make loops */
644     for (step = 1; step < sqmin; ++step)
645     {
646        fill_loop(step, midx, midy, arrX, arrY, arrN, 0.0);
647     }
648 /* tidy up arround the edges */
649     fill_loop(sqmin-1, lx+sqmin-1, ly+sqmin-1, arrX, arrY, arrN, 0.0);
650     fill_loop(sqmin-1, midx, ly+sqmin-1, arrX, arrY, arrN, 0.0);
651     fill_loop(sqmin-1, lx+sqmin-1, midy, arrX, arrY, arrN, 0.0);
652 
653     for (i = ly; i < uy; ++i)
654     {
655         im_put_rowf(arrN[i],im_new,i,lx,ux);
656     }
657     farray_free((void *) arrN,ly,lx,uy,ux);
658     farray_free((void *) arrX,ly,lx,uy,ux);
659     farray_free((void *) arrY,ly,lx,uy,ux);
660 
661     return (im_new);
662 }
663 
664 static Imrect *integrate(Imregion *roi, Imrect *imdx, Imrect *imdy)
665 /* integrate from two sets of derivatives constraining along horizontal
666    and vertical axes debugged but no longer used NAT 25/12/2001 */
667 {
668     Imrect *im_new;
669     float  **arrN,**arrY,**arrX;
670     float  *func,*val;
671     float  valave,funcave;
672     int    lx,ly,ux,uy;
673     int    sqmin,sqmax;
674     int    step,run,chi;
675     int    i, j, midx, midy;
676     float  end,diff,c;
677 
678     lx = roi->lx;
679     ux = roi->ux;
680     ly = roi->ly;
681     uy = roi->uy;
682     midx = (ux + lx)/2;
683     midy = (uy + ly)/2;
684     sqmin = imin((ux - lx)/2,(uy - ly)/2);
685     sqmax = imax((ux - lx)/2,(uy - ly)/2);
686 
687     im_new = im_alloc(imdx->height, imdy->width, roi, float_v);
688     arrN = farray_alloc(ly,lx,uy,ux);
689     arrY = farray_alloc(ly,lx,uy,ux);
690     arrX = farray_alloc(ly,lx,uy,ux);
691     func  = fvector_alloc(0,2*(ux-lx)+2*(uy-ly));
692     val   = fvector_alloc(0,2*(ux-lx)+2*(uy-ly));
693 
694     for (i = ly; i < uy; ++i)
695     {
696         im_get_rowf(arrX[i], imdx, i, lx, ux);
697         im_get_rowf(arrY[i], imdy, i, lx, ux);
698         for (j = lx; j< ux; ++j) arrN[i][j] = 0.0;
699     }
700 
701     /* centre */ 
702     arrN[midy][midx] = 0.0;
703     /* make loops */
704     for (step = 1; step < sqmin; ++step)
705     {
706         /* step out along the axes*/
707         arrN[midy][midx+step] = 
708                      arrN[midy][midx+step-1] + arrX[midy][midx+step];
709         arrN[midy][midx-step] = 
710                      arrN[midy][midx-step+1] - arrX[midy][midx-step+1];
711         arrN[midy+step][midx] = 
712                      arrN[midy+step-1][midx] + arrY[midy+step][midx];
713         arrN[midy-step][midx] = 
714                      arrN[midy-step+1][midx] - arrY[midy-step+1][midx];
715 
716     /* right bottom */
717         run = 1;
718         funcave = 0.0;
719         valave  = 0.0;
720         c = 0.0;
721         func[run] = arrN[midy][midx+step] + arrY[midy+1][midx+step];
722         val[run] = arrN[midy+1][midx+step-1] + arrX[midy+1][midx+step];
723         if (step != 1)
724         {
725            for (i = midy+2; i < midy+step; ++i)
726            {
727               ++run;
728               func[run] = func[run-1] + arrY[i][midx+step];
729               funcave += func[run];
730               val[run] = arrN[i][midx+step-1] + arrX[i][midx+step];
731               valave += val[run];
732            }
733            ++run;
734            val[run] = 0.0;
735            func[run] = func[run-1] + arrY[midy+step][midx+step];
736          /* bottom right */
737            for (j = midx+step-1; j > midx; --j)
738            {
739               ++run;
740               func[run] = func[run-1] - arrX[midy+step][j+1];
741               funcave += func[run];
742               val[run] = arrN[midy+step-1][j] + arrY[midy+step][j];
743               valave += val[run];
744            }
745        }
746        end = func[run] - arrX[midy+step][midx+1];
747        diff = (arrN[midy+step][midx] - end)/(2*step);
748        funcave = 0.0;
749        for (i = 1; i<= run; ++i)
750        {
751            func[i] = func[i] + i*diff;
752            if (i != step) funcave += func[i];
753        }
754        if (run != 1) c = (valave - funcave)/((float)run-1);
755        run = 0;
756        for (i = midy+1; i < midy+step; ++i)
757        {
758            ++run;
759            arrN[i][midx+step] = (func[run] - c + val[run])/2.0;
760        }
761        ++run;
762        arrN[midy+step][midx+step] = func[run] - c;
763        for (j = midx+step-1; j > midx; --j)
764        {
765           ++run;
766           arrN[midy+step][j] = (func[run] - c + val[run])/2.0;
767        }
768 
769        /* bottom left */
770        run = 1;
771        valave  = 0.0;
772        c = 0.0;
773        func[run] = arrN[midy+step][midx] - arrX[midy+step][midx];
774        val[run] = arrN[midy+step-1][midx-1] + arrY[midy+step][midx-1];
775        if (step != 1)
776        {
777            for (j = midx-2; j> midx-step; --j)
778            {
779                ++run;
780                func[run] = func[run-1] - arrX[midy+step][j+1];
781                val[run] = arrN[midy+step-1][j] + arrY[midy+step][j];
782                valave += val[run];
783            }
784            ++run;
785            val[run] = 0.0;
786            func[run] = func[run-1] - arrX[midy+step][midx-step+1];
787            /* left bottom */
788            for (i = midy+step-1; i > midy; --i)
789            {
790               ++run;
791               func[run] = func[run-1] - arrY[i+1][midx-step];
792               val[run] = arrN[i][midx-step+1] - arrX[i][midx-step+1];
793               valave += val[run];
794            }
795        }
796        end = func[run] - arrY[midy+1][midx-step];
797        diff = (arrN[midy][midx-step] - end)/(2*step);
798        funcave = 0.0;
799        for (i = 1; i<= run; ++i)
800        {
801            func[i] = func[i] + i*diff;
802            if (i != step) funcave += func[i];
803        }
804        if (run != 1) c = (valave - funcave)/((float)run-1);
805        run = 0;
806        for (j = midx-1; j> midx-step; --j)
807        {
808           ++run;
809           arrN[midy+step][j] = (func[run] - c + val[run])/2.0;
810        }
811        ++run;
812        arrN[midy+step][midx-step] = func[run] - c + val[run];
813        for (i = midy+step-1; i > midy; --i)
814        {
815           ++run;
816           arrN[i][midx-step] = (func[run] - c + val[run])/2.0;
817        }
818 
819        /* left top */
820        run = 1;
821        valave  = 0.0;
822        c = 0.0;
823        func[run] = arrN[midy][midx-step] - arrY[midy][midx-step];
824        val[run] = arrN[midy-1][midx-step+1] - arrX[midy-1][midx-step+1];
825        if (step != 1)
826        {
827            for (i = midy-2; i > midy-step; --i)
828            {
829                ++run;
830                func[run] = func[run-1] - arrY[i+1][midx-step];
831                val[run] = arrN[i][midx-step+1] - arrX[i][midx-step+1];
832                valave += val[run];
833            }
834            ++run;
835            val[run] = 0.0;
836            func[run] = func[run-1] - arrY[midy-step+1][midx-step];
837            /* top left */
838            for (j = midx-step+1; j < midx; ++j)
839            {
840               ++run;
841               func[run] = func[run-1] + arrX[midy-step][j];
842               val[run] = arrN[midy-step+1][j] - arrY[midy-step+1][j];
843               valave += val[run];
844            }
845        }
846        end = func[run] + arrX[midy-step][midx];
847        diff = (arrN[midy-step][midx] - end)/(2*step);
848        funcave = 0.0;
849        for (i = 1; i<= run; ++i)
850        {
851            func[i] = func[i] + i*diff;
852            if (i != step) funcave += func[i];
853        }
854        if (run != 1) c = (valave - funcave)/((float)run-1);
855        run = 0;
856        for (i = midy-1; i > midy-step; --i)
857        {
858            ++run;
859            arrN[i][midx-step] = (func[run] - c + val[run])/2.0;
860        }
861        ++run;
862        arrN[midy-step][midx-step] = func[run] - c + val[run];
863        for (j = midx-step+1; j < midx; ++j)
864        {
865           ++run;
866           arrN[midy-step][j] = (func[run] - c + val[run])/2.0;
867        }
868        /* top right */
869        run = 1;
870        valave  = 0.0;
871        c = 0.0;
872        func[run] = arrN[midy-step][midx] + arrX[midy-step][midx+1];
873        val[run] = arrN[midy-step+1][midx+1] - arrY[midy-step+1][midx+1];
874        if (step!=1)
875        {
876            for (j = midx+2; j < midx+step; ++j)
877            {
878                ++run;
879                val[run] = 0.0;
880                func[run] = func[run-1] + arrX[midy-step][j];
881                val[run] = arrN[midy-step+1][j] - arrY[midy-step+1][j];
882                valave += val[run];
883            }
884            ++run;
885            val[run] = 0.0;
886            func[run] = func[run-1] + arrX[midy-step][midx+step];
887           /* right top */
888            for (i = midy-step+1; i < midy; ++i)
889            {
890               ++run;
891               func[run] = func[run-1] + arrY[i][midx+step];
892               val[run] = arrN[i][midx+step-1] + arrX[i][midx+step];
893               valave += val[run];
894            }
895        }
896        end = func[run] + arrY[midy][midx+step];
897        diff = (arrN[midy][midx+step] - end)/(2*step);
898        funcave = 0.0;
899        for (i = 1; i<= run; ++i)
900        {
901            func[i] = func[i] + i*diff;
902            if (i != step) funcave += func[i];
903        }
904        if (run != 1) c = (valave - funcave)/((float)run-1);
905        run = 0;
906        for (i = midx+1; i < midx+step; ++i)
907        {
908           ++run;
909           arrN[midy-step][i] = (func[run] - c + val[run])/2.0;
910        }
911        ++run;
912        arrN[midy-step][midx+step] = func[run] - c + val[run];
913        for (i = midy-step+1 ; i < midy; ++i)
914        {
915            ++run;
916            arrN[i][midx+step] = (func[run] - c + val[run])/2.0;
917        }
918     }
919 /* hack for non-square images removed see im_square() NAT 20/1/2001 */
920     for (i = ly; i < uy; ++i)
921     {
922         im_put_rowf(arrN[i],im_new,i,lx,ux);
923     }
924     fvector_free((void *) func,0);
925     fvector_free((void *) val,0);
926     farray_free((void *) arrN,ly,lx,uy,ux);
927     farray_free((void *) arrY,ly,lx,uy,ux);
928     farray_free((void *) arrX,ly,lx,uy,ux);
929 
930     return (im_new);
931 }
932 
933 double im_corscale(Imrect *im0, Imrect *corr, double noise)
934 /* routine to find the point of maximum difference to 1/x distribution */
935 /* using Bhattacharrya overap NAT 02/01/2002 */
936 {
937     Imregion roi;
938     int lx, ux, ly, uy, i, j, k, maxk, mink;
939     float immax, immin;
940     float *error;
941     float *rowimpix, *rowmod;
942     double diffx, impix, mod, fac, conts, range;
943     double scale = 2.6/CORSCAN, x;
944     double centre, h0, hm1, h1;
945     shistogram *hists[CORSCAN+1];
946 
947     if (noise <= 0.0) return (1.0);
948     error = fvector_alloc(0,CORSCAN);
949     imf_minmax(im0, &immin, &immax);
950     range = 0.5*(immax-immin);
951 
952     for (k=0;k<CORSCAN;k++)
953     {
954         hists[k] = hbook1(k,"scale hist",immin-range, immax+range,
955                           2.0*(immax-immin)/noise); 
956         hreset(hists[k]);
957     }
958     hists[CORSCAN] = hbook1(CORSCAN, "performance", 0.0, 2.6, CORSCAN);
959 
960     roi = *(Imregion *)(im0->region); 
961     lx = roi.lx; ux = roi.ux; ly = roi.ly; uy = roi.uy;
962     rowimpix = fvector_alloc(lx, ux);
963     rowmod = fvector_alloc(lx, ux);
964     
965     for (k = 0; k<CORSCAN; k++)
966     {
967         error[k] = 0.0;
968     }
969     for (j = ly+1; j< uy-1; j++)
970     {
971         im_get_rowf(rowimpix, im0, j, lx, ux);
972         im_get_rowf(rowmod, corr, j, lx, ux);
973         for (i=lx+1; i<ux-1; i++)
974         {
975              impix = rowimpix[i];
976              mod = rowmod[i];
977              for (k = 0; k<CORSCAN; k++)
978              {
979                  fac = (k+0.5)*scale; 
980                  diffx = impix*exp(-fac*mod);
981                  hfill1s(hists[k],diffx,1.0);
982              }
983         }
984     } 
985     fvector_free(rowimpix, lx);
986     fvector_free(rowmod, lx);
987 
988     for (k = 0; k<CORSCAN; k++)
989     {
990         fac = (k+0.5)*scale; 
991         for (j=0,x=hists[k]->xmin+0.5*hists[k]->xincr; j< hists[k]->xbins;
992              j++,x+=hists[k]->xincr)
993         {
994             conts = hists[k]->array[0][j]; 
995             error[k] -= sqrt(conts/(fabs(x)+0.5*noise));
996 /*
997    alternative standard information measure doesn't work due to lack
998    of required invariances !!!
999             error[k] += conts*log(1.0+conts);
1000 */
1001         }
1002         hfill1(hists[CORSCAN],fac,error[k]-error[0]);
1003     }
1004     graph_hfit(imcalc_graph_tv(), hists[CORSCAN]);
1005 
1006 /* initialise to defaults for cases where no local maxima are found */
1007     maxk =0;
1008     for (k=1; k<CORSCAN;k++)
1009     {
1010         if (error[k]>error[0]) maxk = k;
1011     }
1012 
1013 /* now find nearest local maxima to fac= 1.0 */
1014     centre = 0.0;
1015     for (k = 0; k<(CORSCAN/2-1); k++)
1016     {
1017         if (error[CORSCAN/2-k]>error[CORSCAN/2-k-1]
1018           &&error[CORSCAN/2-k]>error[CORSCAN/2-k+1])
1019         {
1020             maxk = CORSCAN/2-k;
1021             h0  = error[maxk];
1022             h1  = error[maxk+1];
1023             hm1 = error[maxk-1];
1024             centre = (hm1-h1)/(2.0*(h1+hm1)-4.0*h0);
1025             break;
1026         }
1027         if (error[CORSCAN/2+k]>error[CORSCAN/2+k-1]
1028           &&error[CORSCAN/2+k]>error[CORSCAN/2+k+1])
1029         {
1030             maxk = CORSCAN/2+k;
1031             h0  = error[maxk];
1032             h1  = error[maxk+1];
1033             hm1 = error[maxk-1];
1034             centre = (hm1-h1)/(2.0*(h1+hm1)-4.0*h0);
1035             break;
1036         }
1037     }
1038 /*
1039     basic debug stuff
1040     hprint(stdout, hists[maxk]);
1041 */
1042     fac = (centre + maxk + 0.5)*scale;
1043     format(" scale factor %f = %f \n", fac, error[maxk]);
1044 
1045     for (k = 0; k<CORSCAN+1; k++) hfree(hists[k]);
1046     fvector_free((void *)error,0);
1047     return (fac);
1048 }
1049 
1050 double im_corscale2(Imrect *im0, Imrect *corr, double noise)
1051 /* routine to determine best scaling factor of correction image based
1052    upon robust least squares of x and y image derivatives NAT 24/12/2001 */
1053 {
1054     Imrect *im1, *im2, *im3, *im4, *im5, *im6, *im7, *im8;
1055     Imrect *moddiffx, *moddiffy;
1056     Imregion roi;
1057     int lx, ux, ly, uy, i, j, k, mink;
1058     double modx, mody, diffx, diffy, imfx, imfy, sumx, sumy, fac, 
1059            errfacx, errfacy;
1060     double thresh = STAT_LIM*noise*sqrt(4.0), weight, weightsum=0.0;
1061     double scale = 2.0/CORSCAN;
1062     shistogram *hists[CORSCAN+1];
1063     float *error;
1064 
1065     hists[CORSCAN] = hbook1(CORSCAN, "performance", 0.0, CORSCAN, CORSCAN);
1066     error = fvector_alloc(0,CORSCAN);
1067     im1 = im_bshift(im0,0,1);
1068     im2 = im_bshift(im0,1,0);
1069 
1070     im3 = im_sum(im1,im0);
1071     im4 = im_sum(im2,im0);
1072 
1073     im5 = im_diff(im0,im1);
1074     im6 = im_diff(im0,im2);
1075     im_free(im1); im_free(im2);
1076 
1077 
1078     im7 = im_bshift(corr,0,1);
1079     im8 = im_bshift(corr,1,0);
1080 
1081     moddiffx = im_diff(corr,im7);
1082     moddiffy = im_diff(corr,im8);
1083 
1084     roi = *(Imregion *)(moddiffy->region); 
1085     lx = roi.lx; ux = roi.ux; ly = roi.ly; uy = roi.uy;
1086     
1087     for (k = 0; k<CORSCAN; k++)
1088     {
1089         error[k] = 0.0;
1090     }
1091     for (i=lx+1; i<ux-1; i++)
1092     {
1093         for (j = ly+1; j< uy-1; j++)
1094         {
1095              IM_PIX_GET(moddiffx, j, i, modx);
1096              IM_PIX_GET(moddiffy, j, i, mody);
1097              IM_PIX_GET(im3, j, i, sumx);         
1098              IM_PIX_GET(im4, j, i, sumy);         
1099              IM_PIX_GET(im5, j, i, imfx);
1100              IM_PIX_GET(im6, j, i, imfy);
1101              weight = imfx*imfx + imfy*imfy;
1102              weight = 1.0 - (sqrt(weight)/thresh);
1103              if (weight <= 0.0) weight = 0.0;
1104              else weight = 1.0;
1105              weightsum += weight;
1106 
1107              for (k = 0; k<CORSCAN; k++)
1108              {
1109                  fac = k*scale; 
1110                  diffx = (imfx - 0.5*sumx*fac*modx)/noise;
1111                  diffx *= diffx;
1112                  errfacx = 2.0*(1.0 + 0.25*fac*fac*modx*modx); 
1113                  if (diffx/errfacx >= 25.0) diffx = 25.0*errfacx;
1114 
1115                  diffy = (imfy - 0.5*sumy*fac*mody)/noise;
1116                  diffy *= diffy;
1117                  errfacy = 2.0*(1.0 + 0.25*fac*fac*mody*mody); 
1118                  if (diffy/errfacy >= 25.0) diffy = 25.0*errfacy;
1119 
1120                  error[k] += weight*(diffy/errfacy + diffx/errfacx); 
1121              }
1122         }
1123     } 
1124     mink = 0;
1125     for (k = 0; k<CORSCAN; k++)
1126     {
1127         if (error[k]<error[mink]) mink = k;
1128         hfill1(hists[CORSCAN],k+0.5,error[k]-error[0]);
1129     }
1130     graph_hfit(imcalc_graph_tv(), hists[CORSCAN]);
1131 
1132     fac = mink*scale;
1133     format(" scale factor %f = %f \n", fac, 
1134           sqrt(error[mink]/2.0)/sqrt(weightsum));
1135     im_free(moddiffx);
1136     im_free(moddiffy);
1137     im_free(im3);
1138     im_free(im4);
1139     im_free(im5);
1140     im_free(im6);
1141     fvector_free((void *)error,0);
1142     return (fac);
1143 }
1144 
1145 Imrect *im_integrate(Imrect *imdxf, Imrect *imdyf)
1146 {
1147     Imrect *im_new, *im_new1, *im_new2, *im_new3;
1148     Imregion *roi;
1149     
1150     if (imdxf == NULL || imdyf == NULL)
1151     {
1152         error("no images provided to imf_integrate",warning);
1153         return(NULL);
1154     }
1155 
1156     roi = imdxf->region;
1157 
1158     im_new = integrate2(roi,imdxf,imdyf);
1159 
1160     return(im_new);
1161 }
1162 
1163 

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