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

Linux Cross Reference
Tina6/tina-libs/tina/medical/medNorm_base.c

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

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  * 
 23  **********
 24  * 
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/medical/medNorm_base.c,v $
 27  * Date    :  $Date: 2008/10/13 15:34:03 $
 28  * Version :  $Revision: 1.11 $
 29  * CVS Id  :  $Id: medNorm_base.c,v 1.11 2008/10/13 15:34:03 paul Exp $
 30  *  
 31  * Author  :  paul.bromiley@manchester.ac.uk
 32  *
 33  * Notes :  
 34  * 
 35  * Modifications NAT 25/5/2001:
 36  * 
 37  * - Use of a single limit (STAT_LIM) on all statistical NULL hypotheses.
 38  * - A separate limit for low information INF_LIM.
 39  * - Normalisation software modified to improve matching of statistical
 40  *     tests to algorithm stability and eliminate edge threshold.
 41  * - Regularisation for indeterminate regions modified to reduce bias
 42  *     towards null gain change.
 43  * - Second differential integration routine addded to reduce 
 44  *     reconstruction errors
 45  * 
 46  * Modifications NAT 25/12/2001:
 47  *
 48  * - Slope estimation routine enf_integ() added to enforce integrability.
 49  *     routine added to maximise information content im_corscale(), based
 50  *     upon Bhattacharrya overlap with scale invariant histogram distribution.
 51  *     NAT 25/12/2001 
 52  *
 53  * Modifications PAB 13/10/2008:
 54  *
 55  * - Added complex image processing, where two images are corrected jointly
 56  *     (see TINA Memo No. 2007-003
 57  *
 58  *
 59  *********
 60 */  
 61 
 62 #include "medNorm_base.h"
 63 
 64 #if HAVE_CONFIG_H
 65 #include <config.h>
 66 #endif
 67     
 68 #include <stdio.h>
 69 #if HAVE_SYS_PARAM_H
 70 #include <sys/param.h>
 71 #endif
 72 
 73 #include <string.h>
 74 #include <math.h>
 75 #include <float.h>
 76 
 77 #include <tina/sys/sysDef.h>
 78 #include <tina/sys/sysPro.h>
 79 #include <tina/math/mathDef.h>
 80 #include <tina/math/mathPro.h>
 81 #include <tina/image/imgDef.h>
 82 #include <tina/image/imgPro.h>
 83 #include <tina/geometry/geomPro.h>
 84     
 85 
 86 #define STAT_LIM 2.5
 87 #define INF_LIM 200.0
 88 #define CORSCAN 26
 89 
 90 
 91 /* mask of intensities < 3 sigma */
 92 static Imrect * thresh_slice(Imrect * im, double sig_noise)
 93 {
 94     Imrect *im1,*im2,*im3,*imcut;
 95 
 96     im1 = im_thresh((STAT_LIM*(sig_noise)),im);
 97     im2 = im_thresh((-STAT_LIM*(sig_noise)),im);
 98     im3 = im_diff(im2,im1);
 99     im_free(im1); im_free(im2);
100     /* fudge for inversion recovery borders */
101     im1 = im_thresh(4096.0,im);
102     imcut = im_sum(im1,im3);
103     im_free(im1); im_free(im3);
104 
105     return (imcut);
106 }
107 
108 /* mask of edges > STAT_LIM  */
109 static Imrect *edge_slice(Imrect *image, double sig_noise) 
110 {
111     Imrect *im,*im2,*im3,*im4;
112     
113     im = imf_gauss(image, 1.0, 0.01);
114 
115     im2 = imf_diffx(im);
116     im3 = im_sqr(im2);
117     im_free(im2); 
118    
119     im2 = imf_diffy(im);
120     im4 = im_sqr(im2);
121     im_free(im); im_free(im2);
122 
123     im = im_sum(im3,im4);
124     im2 = im_sqrt(im);
125     im_free(im); im_free(im3); im_free(im4);
126 
127     im = im_thresh((STAT_LIM*sig_noise/sqrt(2)),im2);
128     im_free(im2);
129 
130     return (im);
131 }
132 
133 static Imrect *zero_cut(double thresh, Imrect * cut) 
134 {
135     Imrect *im,*im1,*imcut;
136 
137     im1 = im_mod(cut);
138     im = im_bthresh(thresh, im1);
139     imcut = im_cast(im, uchar_v); 
140     im_free(im); im_free(im1); 
141 
142     return (imcut);
143 }
144 
145 static Imrect *zero_mask(Imrect * image)  /* mask of voxels = 0 */
146 {
147     Imrect *im_new,*im1,*im2;
148     
149     im1 = zero_cut(FLT_MIN,image);
150     im2 = im_minus(im1);
151     im_new = im_add(1.0,im2);    
152     im_free(im1);  im_free(im2);
153     
154     return (im_new);
155     
156 }
157 
158 static Imrect * edge_remove(Imrect * imcut, Imrect * image)
159 {
160     Imrect *im,*new_image;
161 
162     im = im_prod(image,imcut);
163     new_image = im_diff(image,im);
164     im_free(im);
165     return (new_image);
166 }
167 
168 /*                                    */
169 /* mask combining two slices of edges */
170 /*                                    */
171 void edge_mask(double sig_noise, Imrect *im1, Imrect *im2, Imrect **cut,
172                 Imrect **image1, Imrect **image2)
173 {
174     Imrect *cut1,*cut2,*cut3;
175     Imrect *cut_temp1,*cut_temp2;
176     Imrect *im1_temp,*im2_temp;
177 
178     im1_temp = imf_median(im1);
179     im2_temp = imf_median(im2);
180     cut1 = thresh_slice(im1_temp, sig_noise);
181     cut2 = edge_slice(im1_temp, sig_noise);
182     cut3 = im_sum(cut1,cut2);
183     im_free(cut1);  im_free(cut2);
184     cut1 = zero_mask(im1_temp);
185     cut_temp1 = im_sum(cut1,cut3);
186     im_free(cut1);  im_free(cut3);
187     /* and again... */
188     cut1 = thresh_slice(im2_temp, sig_noise);
189     cut2 = edge_slice(im2_temp, sig_noise);
190     cut3 = im_sum(cut1,cut2);
191     im_free(cut1);  im_free(cut2);
192     cut1 = zero_mask(im2_temp);
193     cut_temp2 = im_sum(cut1,cut3);
194     im_free(cut1);  im_free(cut3);
195 
196     cut1 = im_sum(cut_temp1,cut_temp2);
197     im_free(cut_temp1);    im_free(cut_temp2);
198     cut2 = zero_cut(FLT_MIN,cut1);
199     im_free(cut1);
200 
201     *image1 = edge_remove(cut2,im1_temp);
202     *image2 = edge_remove(cut2,im2_temp);
203     *cut = cut2;
204     im_free(im1_temp); im_free(im2_temp);
205     return;
206 }
207 
208 static shistogram *imf_hist_ratio(double sig_noise, double range,
209                            Imrect *ima, Imrect* imb, Imrect *immask,
210                            Imregion *roi,int nbin)
211 {
212     shistogram *imhist;
213     float *rowa, *rowb, *rowmask;
214     int lx, ux, ly, uy;
215     int width, height;
216     int i, j;
217     float shift = range/(float)nbin;
218     double ratio, log_ratio, errweight;
219 
220     if (ima == NULL || imb == NULL)
221         return(NULL);
222     if (roi == NULL)
223         return(NULL);
224     lx = roi->lx;
225     ux = roi->ux;
226     ly = roi->ly;
227     uy = roi->uy;
228 
229     width = ima->width;
230     height = ima->height;
231 
232     imhist = hbook1(0,"image contents\0",(float)(-range),
233                 (float)(range), nbin);
234 
235     rowa = fvector_alloc(lx, ux);
236     rowb = fvector_alloc(lx, ux);
237     rowmask = fvector_alloc(lx, ux);
238 
239     for (i = ly; i < uy; ++i)
240     {
241         im_get_rowf(rowa, ima, i, lx, ux);
242         im_get_rowf(rowb, imb, i, lx, ux);
243         im_get_rowf(rowmask, immask, i, lx, ux);
244         for (j = lx; j < ux; ++j)
245         {
246             if (rowb[j]==0.0) ratio = FLT_MAX;
247             else ratio = rowa[j]/rowb[j];
248             if (ratio>0.0) log_ratio = log(ratio); 
249             else log_ratio = -FLT_MAX;
250 /* weight computed based one error on log-ratio NAT 20/1/2001 */ 
251             errweight = rowa[j]*rowa[j]*rowb[j]*rowb[j]
252                       /(rowa[j]*rowa[j]+rowb[j]*rowb[j]);
253             errweight /= (sig_noise*sig_noise);
254  
255             if(rowmask[j]!=0.0)
256             {
257                hfill1s(imhist,log_ratio-shift,errweight);
258                hfill1s(imhist,log_ratio+shift,errweight);
259             }
260         }
261     }
262  
263     fvector_free(rowa, lx);
264     fvector_free(rowb, lx);
265     fvector_free(rowmask, lx);
266     return(imhist);
267 }
268 
269 double im_fraction(Imrect *cut, Imrect *image1, Imrect *image2, int nbin, 
270                 double sig_noise) 
271 {
272     Imrect *immask;
273     double normal;
274     shistogram *hist=NULL;
275     Imregion *roi=NULL;
276     float maxx=0.0;
277 
278     immask   = zero_mask(cut);
279     roi = roi_inter(image1->region,image2->region);
280     hist = imf_hist_ratio(sig_noise,1.0,image1,image2,immask,roi,nbin);
281     /* removed during promotion to tina-libs
282        a.lacey@man.ac.uk 23.05.03
283        graph_hfit(imcalc_graph_tv_get(), hist); */
284 
285     normal = (double)hnearmax1(hist,&maxx);
286     if (maxx>0.0 || maxx <=0.0)
287        normal = (double)maxx;
288     else
289     {
290        normal = 0.0;
291        format("bad ratio calculated in fraction()");
292     }
293     im_free(immask);
294     rfree(roi);
295     hfree(hist);
296     return(normal);
297 }
298 
299 static void err_map(Imrect *im, Imrect *imx, Imrect *imy, double sig_noise,
300                        double range, Imrect **imerrx, Imrect **imerry)
301 {
302     Imrect *imrx,*imry;
303     Imregion *roi;
304     float *rowi,*rowx,*rowy,*rowfinx,*rowfiny;
305     int i,j,lx,ux,ly,uy;
306     double weight;
307     double thresh=sig_noise*range*sqrt(2);
308  
309     roi = im->region;
310  
311     lx = roi->lx;
312     ux = roi->ux;
313     ly = roi->ly;
314     uy = roi->uy;
315  
316     if (lx == ux || ly == uy)
317         return;
318 
319     imrx = im_alloc(im->height, im->width, roi, float_v);
320     imry = im_alloc(im->height, im->width, roi, float_v);
321     rowi = fvector_alloc(lx, ux);
322     rowx = fvector_alloc(lx, ux);
323     rowy = fvector_alloc(lx, ux);
324     rowfinx = fvector_alloc(lx, ux);
325     rowfiny = fvector_alloc(lx, ux);
326  
327     for (i = ly+1; i < uy-1; ++i)
328     {
329         im_get_rowf(rowi, im, i, lx, ux);
330         im_get_rowf(rowx, imx, i, lx, ux);
331         im_get_rowf(rowy, imy, i, lx, ux);
332         for (j = lx+1 ; j < ux-1; ++j)
333         {
334            if ((fabs(rowi[j]) <= 0.01*sig_noise) ||
335                (fabs(rowx[j]) <= 0.01*sig_noise) ||
336                (fabs(rowy[j]) <= 0.01*sig_noise) || 
337                (fabs(rowi[j]-rowx[j]) > thresh) ||
338                (fabs(rowi[j]-rowy[j]) > thresh) )
339            { 
340               rowfinx[j] = 0.0;
341               rowfiny[j] = 0.0;
342            }
343            else 
344            {
345               rowfinx[j] = (16*sig_noise*sig_noise
346                             *(rowi[j]*rowi[j]+ rowx[j]*rowx[j]))
347                             /(pow((rowi[j]+rowx[j]),4));
348               rowfiny[j] = (16*sig_noise*sig_noise
349                             *(rowi[j]*rowi[j] + rowy[j]*rowy[j]))
350                             /(pow((rowi[j]+rowy[j]),4));
351 /* weighting to eliminate any edges */
352               weight = (rowi[j]-rowx[j])*(rowi[j]-rowx[j])
353                       +(rowi[j]-rowy[j])*(rowi[j]-rowy[j]);
354               weight = 1.0 - (sqrt(weight)/thresh);
355               if (weight <= 0.0) weight = 0.0;
356               weight = sqrt(weight);
357               rowfinx[j] = weight/rowfinx[j];
358               rowfiny[j] = weight/rowfiny[j];
359             }
360         }
361         im_put_rowf(rowfinx, imrx, i, lx, ux);
362         im_put_rowf(rowfiny, imry, i, lx, ux);
363     }
364                 
365     fvector_free(rowi, lx);
366     fvector_free(rowx, lx);
367     fvector_free(rowy, lx);
368     fvector_free(rowfinx, lx);
369     fvector_free(rowfiny, lx);
370                 
371     *imerrx = imrx;
372     *imerry = imry;
373 
374     return;
375 }
376 
377 static Imrect* lsf_iter(Imrect *im,int iter, double sigma)
378 {
379     Imrect *imsm=NULL,*im1=NULL;
380     int     i;
381     
382     im1 = im_copy(im);
383     for(i = 1; i<= iter; ++i)
384     {
385        if (i>1)
386        {
387           im1 = imsm;
388        }
389        imsm = imf_lsf_smooth(im1,sigma);
390        im_free(im1);
391     }
392     return(imsm);
393 }
394 
395 void smooth_slopes(double sig_noise, double nsmear, Imrect *image,
396                    Imrect **imdx, Imrect **imdy, Imrect **corr, 
397                    Imrect **imxs, Imrect **imys)
398 /* routine to calculate smooth slopes imdx and imdy from image  */
399 /* using noise estimate sig_noise and regional smoothing nsmear */
400 {
401     Imrect *im3=NULL, *im4=NULL, *imdxdy=NULL, *imdydx=NULL;
402     Imrect *imfx=NULL, *imfy=NULL;
403     Imrect *im=NULL, *im0=NULL, *im1=NULL, *im2=NULL;
404     Imrect *imdx1=NULL, *imdy1=NULL, *imdx2=NULL, *imdy2=NULL;
405     Imrect *imx=NULL, *imy=NULL, *imx_temp=NULL, *imy_temp=NULL;
406     Imrect *imdxs=NULL, *imdys=NULL;
407     int niter = 2;
408 
409 /* smooth image to remove outliers NAT 20/1/2001 */
410     im = im_median(image);
411     im0 = im_median(im);
412     sig_noise *= 0.5;
413     im_free(im);
414     im1 = im_bshift(im0,0,1);
415     im2 = im_bshift(im0,1,0);
416 
417     imx_temp = im_sum(im1,im0);
418     imy_temp = im_sum(im2,im0);
419     imx = im_times(0.5,imx_temp);
420     imy = im_times(0.5,imy_temp);
421     im_free(imx_temp);im_free(imy_temp);
422     imdx1 = im_diff(im0,im1);
423     imdy1 = im_diff(im0,im2);
424     err_map(im0,im1,im2,sig_noise,STAT_LIM,&imfx,&imfy);
425     im_free(im0); im_free(im1); im_free(im2);
426 
427     /* scale slopes */
428     imdx2 = imf_div(imdx1,imx,STAT_LIM*sig_noise,STAT_LIM*sig_noise);
429     imdy2 = imf_div(imdy1,imy,STAT_LIM*sig_noise,STAT_LIM*sig_noise);
430     im_free(imdy1);im_free(imdx1); im_free(imx);im_free(imy);
431     *imdx = imf_prod(imdx2,imfx);
432     *imdy = imf_prod(imdy2,imfy);
433 /*
434     stack_push((void*) imfx,IMRECT,im_free);
435     stack_push((void*) imfy,IMRECT,im_free);
436 */
437     im_free(imdy2);im_free(imdx2);
438 
439     *imxs = lsf_iter(imfx,niter,nsmear);
440     *imys = lsf_iter(imfy,niter,nsmear);
441     im_free(imfx); im_free(imfy);
442     imdxs = lsf_iter(*imdx,niter,nsmear);
443     imdys = lsf_iter(*imdy,niter,nsmear);
444     im_free(*imdx); im_free(*imdy);
445 /*  restrict statistical weighting to bias towards zero image slope
446     when there is no accurate information NAT 25/4/2001 */
447     *imdx = imf_div(imdxs,*imxs,INF_LIM/nsmear,INF_LIM/nsmear);
448     *imdy = imf_div(imdys,*imys,INF_LIM/nsmear,INF_LIM/nsmear);
449     im_free(imdxs); im_free(imdys);
450     im3 = im_bshift(*imdy,0,1);
451     im4 = im_bshift(*imdx,1,0);
452     imdydx = im_diff(*imdy, im3);
453     imdxdy = im_diff(*imdx, im4);
454     im_free(im3); im_free(im4);
455     *corr = im_diff(imdxdy,imdydx);
456     im_free(imdxdy); im_free(imdydx);
457 }
458 
459 void enf_integ(Imrect **imdx, Imrect **imdy, Imrect *corr, Imrect *imxs, Imrect *imys)
460 {
461 /*  enforce integrability  by averaging of error on cross derivatives
462                                                           NAT 25/12/2001 */
463     int lx,ly,ux,uy,i,j;
464     Imregion roi;
465     double corrpix, pixval1, pixval2, integx, meanx, integy, meany;
466     double weight, norm;
467     roi = *(Imregion *)(corr->region);
468     lx = roi.lx;
469     ux = roi.ux;
470     ly = roi.ly;
471     uy = roi.uy;
472     for (i=lx;i<ux;i++)
473     {
474       IM_PIX_SET(corr, ly, i, 0.0);
475     }
476 
477     for (j=ly;j<uy;j++)
478     {
479       IM_PIX_SET(corr, j, lx, 0.0);
480     }
481 
482     for (i=lx;i<ux;i++)
483     {
484        integx =0.0;
485        meanx = 0.0;
486        norm = 0.0;
487        for (j=ly; j<uy; j++)
488        {
489           IM_PIX_GET(corr, j, i, pixval1);
490           IM_PIX_GET(imxs, j, i, weight);
491           norm += weight;
492           integx += pixval1;
493           meanx += integx*weight;
494        }
495        meanx /= norm;
496        integx = 0.0;
497        for (j=ly; j<uy; j++)
498        {
499           IM_PIX_GET(corr, j, i, pixval1);
500           IM_PIX_GET(*imdx, j, i, pixval2);
501           integx  += pixval1;
502           corrpix = pixval2 - 0.5*(integx - meanx);
503           IM_PIX_SET(*imdx, j, i, corrpix);
504        }
505     }
506     for (j=ly; j<uy; j++)
507     {
508        integy =0.0;
509        meany = 0.0;
510        norm = 0.0;
511        for (i=lx;i<ux;i++)
512        {
513           IM_PIX_GET(corr, j, i, pixval1);
514           IM_PIX_GET(imys, j, i, weight);
515           norm += weight;
516           integy += pixval1;
517           meany += integy*weight;
518        }
519        meany /= norm;
520        integy = 0.0;
521        for (i=lx; i<ux; i++)
522        {
523           IM_PIX_GET(corr, j, i, pixval1);
524           IM_PIX_GET(*imdy, j, i, pixval2);
525           integy  += pixval1;
526           corrpix = pixval2 + 0.5*(integy - meany);
527           IM_PIX_SET(*imdy, j, i, corrpix);
528        }
529     }
530 }
531 
532 static void fill_loop(int step, int midx, int midy, float **arrX, float **arrY,
533                  float **arrN, int correct)
534 /* reconstruct integral image loop with integrability enforced by linear
535    constraint if not already guaranteed by use of end_integ() NAT 25/12/25 */
536 {
537     int i,j;
538     int linex, liney;
539     float corner1, corner2, corner3, corner4;
540     float biasweight, bias, cumweight;
541 
542     linex = midx-step;
543     liney = midy-step;
544     corner1 = arrN[liney+1][linex+1] + 0.5 *
545             (-arrY[liney+1][linex+1]
546              -arrX[liney+1][linex+1]
547              -arrY[liney+1][linex]
548              -arrX[liney][linex+1]);
549     linex = midx+step;
550     liney = midy-step;
551     corner2 = arrN[liney+1][linex-1] + 0.5 *
552             (-arrY[liney+1][linex-1]
553              +arrX[liney+1][linex]
554              -arrY[liney+1][linex]
555              +arrX[liney][linex]);
556     linex = midx+step;
557     liney = midy+step;
558     corner3 = arrN[liney-1][linex-1] + 0.5 *
559             (+arrY[liney][linex-1]
560              +arrX[liney-1][linex]
561              +arrY[liney][linex]
562              +arrX[liney][linex]);
563     linex = midx-step;
564     liney = midy+step;
565     corner4 = arrN[liney-1][linex+1] + 0.5 *
566             (+arrY[liney][linex+1]
567              -arrX[liney-1][linex+1]
568              +arrY[liney][linex]
569              -arrX[liney][linex+1]);
570     biasweight = 2*step+1;
571 
572     liney = midy-step;
573     arrN[liney][midx-step] = corner1;
574     for (j=midx-step+1;j<=midx+step;j++)
575     {
576         arrN[liney][j] = arrN[liney][j-1] + arrX[liney][j];
577     }
578     if (correct)
579     {
580         bias = corner2 - arrN[liney][midx+step];
581         cumweight = 0.0;
582         for (j=midx-step+1;j<=midx+step;j++)
583         {
584             cumweight += 1.0;
585             arrN[liney][j] += bias*cumweight/biasweight;
586         }
587         for (j=midx-step+1;j<midx+step;j++)
588         {
589             arrN[liney][j] = (arrN[liney][j] 
590                             + (arrN[liney+1][j] - arrY[liney+1][j]))
591                             /2.0;
592         }
593     }
594 
595     linex = midx+step;
596     arrN[midy-step][linex] = corner2;
597     for (i=midy-step+1;i<=midy+step;i++)
598     {
599         arrN[i][linex] = arrN[i-1][linex] + arrY[i][linex];
600     }
601     if (correct)
602     {
603         bias = corner3 - arrN[midy+step][linex];
604         cumweight = 0.0;
605         for (i=midy-step+1;i<=midy+step;i++)
606         {
607             cumweight += 1.0;
608             arrN[i][linex] += bias*cumweight/biasweight;
609         }
610         for (i=midy-step+1;i<midy+step;i++)
611         {
612             arrN[i][linex] = (arrN[i][linex] 
613                             + (arrN[i][linex-1] + arrX[i][linex]))
614                             /2.0;
615         }
616     }
617 
618     liney = midy+step;
619     arrN[liney][midx+step] = corner3;
620     for (j=midx+step-1;j>=midx-step;j--)
621     {
622         arrN[liney][j] = arrN[liney][j+1] - arrX[liney][j+1];
623     }
624     if (correct)
625     {
626         bias = corner4 - arrN[liney][midx-step];
627         cumweight = 0.0;
628         for (j=midx+step-1;j>=midx-step;j--)
629         {
630             cumweight += 1.0;
631             arrN[liney][j] += bias*cumweight/biasweight;
632         }
633         for (j=midx+step-1;j>midx-step;j--)
634         {
635             arrN[liney][j] = (arrN[liney][j] 
636                             + (arrN[liney-1][j] + arrY[liney][j]))
637                             /2.0;
638         }
639     }
640 
641     linex = midx-step;
642     arrN[midy+step][linex] = corner4;
643     for (i=midy+step-1;i>=midy-step;i--)
644     {
645         arrN[i][linex] = arrN[i+1][linex] - arrY[i+1][linex];
646     }
647     if (correct)
648     {
649         bias = corner1 - arrN[midy-step][linex];
650         cumweight = 0.0;
651         for (i=midy+step-1;i>=midy-step;i--)
652         {
653             cumweight += 1.0;
654             arrN[i][linex] += bias*cumweight/biasweight;
655         }
656         for (i=midy+step-1;i>midy-step;i--)
657         {
658             arrN[i][linex] = (arrN[i][linex]
659                               + (arrN[i][linex+1] - arrX[i][linex+1]))
660                               /2.0;
661         }
662     }
663 }
664 
665 static Imrect *integrate2(Imregion *roi, Imrect *imdx, Imrect *imdy)
666 /* integrate from two sets of derivatives constraining diagonals */
667 {
668     Imrect *im_new;
669     float  **arrN,**arrY,**arrX;
670     int     lx,ly,ux,uy;
671     int     sqmin,sqmax;
672     int     i, j, midx, midy, step;
673 
674     lx = roi->lx;
675     ux = roi->ux;
676     ly = roi->ly;
677     uy = roi->uy;
678     midx = (ux + lx)/2;
679     midy = (uy + ly)/2;
680     sqmin = imin((ux - lx)/2,(uy - ly)/2);
681     sqmax = imax((ux - lx)/2,(uy - ly)/2);
682 
683     im_new = im_alloc(imdx->height, imdy->width, roi, float_v);
684     arrN = farray_alloc(ly,lx,uy,ux);
685     arrY = farray_alloc(ly,lx,uy,ux);
686     arrX = farray_alloc(ly,lx,uy,ux);
687     for (i = ly; i < uy; ++i)
688     {
689         im_get_rowf(arrX[i], imdx, i, lx, ux);
690         im_get_rowf(arrY[i], imdy, i, lx, ux);
691         for (j = lx; j< ux; ++j) arrN[i][j] = 0.0;
692     }
693 
694 /* centre */
695     arrN[midy][midx] = 0.0;
696 /* make loops */
697     for (step = 1; step < sqmin; ++step)
698     {
699        fill_loop(step, midx, midy, arrX, arrY, arrN, 0.0);
700     }
701 /* tidy up arround the edges */
702     fill_loop(sqmin-1, lx+sqmin-1, ly+sqmin-1, arrX, arrY, arrN, 0.0);
703     fill_loop(sqmin-1, midx, ly+sqmin-1, arrX, arrY, arrN, 0.0);
704     fill_loop(sqmin-1, lx+sqmin-1, midy, arrX, arrY, arrN, 0.0);
705 
706     for (i = ly; i < uy; ++i)
707     {
708         im_put_rowf(arrN[i],im_new,i,lx,ux);
709     }
710     farray_free(arrN,ly,lx,uy,ux);
711     farray_free(arrX,ly,lx,uy,ux);
712     farray_free(arrY,ly,lx,uy,ux);
713 
714     return (im_new);
715 }
716 
717 
718 /* Static but not used: comment out for now PAB */
719 /* integrate from two sets of derivatives constraining along horizontal
720    and vertical axes debugged but no longer used NAT 25/12/2001 */
721 /*
722 static Imrect *integrate(Imregion *roi, Imrect *imdx, Imrect *imdy)
723 {
724     Imrect *im_new;
725     float  **arrN,**arrY,**arrX;
726     float  *func,*val;
727     float  valave,funcave;
728     int    lx,ly,ux,uy;
729     int    sqmin,sqmax;
730     int    step,run;
731     int    i, j, midx, midy;
732     float  end,diff,c;
733 
734     lx = roi->lx;
735     ux = roi->ux;
736     ly = roi->ly;
737     uy = roi->uy;
738     midx = (ux + lx)/2;
739     midy = (uy + ly)/2;
740     sqmin = imin((ux - lx)/2,(uy - ly)/2);
741     sqmax = imax((ux - lx)/2,(uy - ly)/2);
742 
743     im_new = im_alloc(imdx->height, imdy->width, roi, float_v);
744     arrN = farray_alloc(ly,lx,uy,ux);
745     arrY = farray_alloc(ly,lx,uy,ux);
746     arrX = farray_alloc(ly,lx,uy,ux);
747     func  = fvector_alloc(0,2*(ux-lx)+2*(uy-ly));
748     val   = fvector_alloc(0,2*(ux-lx)+2*(uy-ly));
749 
750     for (i = ly; i < uy; ++i)
751     {
752         im_get_rowf(arrX[i], imdx, i, lx, ux);
753         im_get_rowf(arrY[i], imdy, i, lx, ux);
754         for (j = lx; j< ux; ++j) arrN[i][j] = 0.0;
755     }
756 
757 
758     arrN[midy][midx] = 0.0;
759 
760     for (step = 1; step < sqmin; ++step)
761     {
762         arrN[midy][midx+step] = 
763                      arrN[midy][midx+step-1] + arrX[midy][midx+step];
764         arrN[midy][midx-step] = 
765                      arrN[midy][midx-step+1] - arrX[midy][midx-step+1];
766         arrN[midy+step][midx] = 
767                      arrN[midy+step-1][midx] + arrY[midy+step][midx];
768         arrN[midy-step][midx] = 
769                      arrN[midy-step+1][midx] - arrY[midy-step+1][midx];
770 
771         run = 1;
772         funcave = 0.0;
773         valave  = 0.0;
774         c = 0.0;
775         func[run] = arrN[midy][midx+step] + arrY[midy+1][midx+step];
776         val[run] = arrN[midy+1][midx+step-1] + arrX[midy+1][midx+step];
777         if (step != 1)
778         {
779            for (i = midy+2; i < midy+step; ++i)
780            {
781               ++run;
782               func[run] = func[run-1] + arrY[i][midx+step];
783               funcave += func[run];
784               val[run] = arrN[i][midx+step-1] + arrX[i][midx+step];
785               valave += val[run];
786            }
787            ++run;
788            val[run] = 0.0;
789            func[run] = func[run-1] + arrY[midy+step][midx+step];
790 
791            for (j = midx+step-1; j > midx; --j)
792            {
793               ++run;
794               func[run] = func[run-1] - arrX[midy+step][j+1];
795               funcave += func[run];
796               val[run] = arrN[midy+step-1][j] + arrY[midy+step][j];
797               valave += val[run];
798            }
799        }
800        end = func[run] - arrX[midy+step][midx+1];
801        diff = (arrN[midy+step][midx] - end)/(2*step);
802        funcave = 0.0;
803        for (i = 1; i<= run; ++i)
804        {
805            func[i] = func[i] + i*diff;
806            if (i != step) funcave += func[i];
807        }
808        if (run != 1) c = (valave - funcave)/((float)run-1);
809        run = 0;
810        for (i = midy+1; i < midy+step; ++i)
811        {
812            ++run;
813            arrN[i][midx+step] = (func[run] - c + val[run])/2.0;
814        }
815        ++run;
816        arrN[midy+step][midx+step] = func[run] - c;
817        for (j = midx+step-1; j > midx; --j)
818        {
819           ++run;
820           arrN[midy+step][j] = (func[run] - c + val[run])/2.0;
821        }
822 
823        run = 1;
824        valave  = 0.0;
825        c = 0.0;
826        func[run] = arrN[midy+step][midx] - arrX[midy+step][midx];
827        val[run] = arrN[midy+step-1][midx-1] + arrY[midy+step][midx-1];
828        if (step != 1)
829        {
830            for (j = midx-2; j> midx-step; --j)
831            {
832                ++run;
833                func[run] = func[run-1] - arrX[midy+step][j+1];
834                val[run] = arrN[midy+step-1][j] + arrY[midy+step][j];
835                valave += val[run];
836            }
837            ++run;
838            val[run] = 0.0;
839            func[run] = func[run-1] - arrX[midy+step][midx-step+1];
840            for (i = midy+step-1; i > midy; --i)
841            {
842               ++run;
843               func[run] = func[run-1] - arrY[i+1][midx-step];
844               val[run] = arrN[i][midx-step+1] - arrX[i][midx-step+1];
845               valave += val[run];
846            }
847        }
848        end = func[run] - arrY[midy+1][midx-step];
849        diff = (arrN[midy][midx-step] - end)/(2*step);
850        funcave = 0.0;
851        for (i = 1; i<= run; ++i)
852        {
853            func[i] = func[i] + i*diff;
854            if (i != step) funcave += func[i];
855        }
856        if (run != 1) c = (valave - funcave)/((float)run-1);
857        run = 0;
858        for (j = midx-1; j> midx-step; --j)
859        {
860           ++run;
861           arrN[midy+step][j] = (func[run] - c + val[run])/2.0;
862        }
863        ++run;
864        arrN[midy+step][midx-step] = func[run] - c + val[run];
865        for (i = midy+step-1; i > midy; --i)
866        {
867           ++run;
868           arrN[i][midx-step] = (func[run] - c + val[run])/2.0;
869        }
870 
871        run = 1;
872        valave  = 0.0;
873        c = 0.0;
874        func[run] = arrN[midy][midx-step] - arrY[midy][midx-step];
875        val[run] = arrN[midy-1][midx-step+1] - arrX[midy-1][midx-step+1];
876        if (step != 1)
877        {
878            for (i = midy-2; i > midy-step; --i)
879            {
880                ++run;
881                func[run] = func[run-1] - arrY[i+1][midx-step];
882                val[run] = arrN[i][midx-step+1] - arrX[i][midx-step+1];
883                valave += val[run];
884            }
885            ++run;
886            val[run] = 0.0;
887            func[run] = func[run-1] - arrY[midy-step+1][midx-step];
888            for (j = midx-step+1; j < midx; ++j)
889            {
890               ++run;
891               func[run] = func[run-1] + arrX[midy-step][j];
892               val[run] = arrN[midy-step+1][j] - arrY[midy-step+1][j];
893               valave += val[run];
894            }
895        }
896        end = func[run] + arrX[midy-step][midx];
897        diff = (arrN[midy-step][midx] - 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 = midy-1; i > midy-step; --i)
907        {
908            ++run;
909            arrN[i][midx-step] = (func[run] - c + val[run])/2.0;
910        }
911        ++run;
912        arrN[midy-step][midx-step] = func[run] - c + val[run];
913        for (j = midx-step+1; j < midx; ++j)
914        {
915           ++run;
916           arrN[midy-step][j] = (func[run] - c + val[run])/2.0;
917        }
918 
919        run = 1;
920        valave  = 0.0;
921        c = 0.0;
922        func[run] = arrN[midy-step][midx] + arrX[midy-step][midx+1];
923        val[run] = arrN[midy-step+1][midx+1] - arrY[midy-step+1][midx+1];
924        if (step!=1)
925        {
926            for (j = midx+2; j < midx+step; ++j)
927            {
928                ++run;
929                val[run] = 0.0;
930                func[run] = func[run-1] + arrX[midy-step][j];
931                val[run] = arrN[midy-step+1][j] - arrY[midy-step+1][j];
932                valave += val[run];
933            }
934            ++run;
935            val[run] = 0.0;
936            func[run] = func[run-1] + arrX[midy-step][midx+step];
937 
938            for (i = midy-step+1; i < midy; ++i)
939            {
940               ++run;
941               func[run] = func[run-1] + arrY[i][midx+step];
942               val[run] = arrN[i][midx+step-1] + arrX[i][midx+step];
943               valave += val[run];
944            }
945        }
946        end = func[run] + arrY[midy][midx+step];
947        diff = (arrN[midy][midx+step] - end)/(2*step);
948        funcave = 0.0;
949        for (i = 1; i<= run; ++i)
950        {
951            func[i] = func[i] + i*diff;
952            if (i != step) funcave += func[i];
953        }
954        if (run != 1) c = (valave - funcave)/((float)run-1);
955        run = 0;
956        for (i = midx+1; i < midx+step; ++i)
957        {
958           ++run;
959           arrN[midy-step][i] = (func[run] - c + val[run])/2.0;
960        }
961        ++run;
962        arrN[midy-step][midx+step] = func[run] - c + val[run];
963        for (i = midy-step+1 ; i < midy; ++i)
964        {
965            ++run;
966            arrN[i][midx+step] = (func[run] - c + val[run])/2.0;
967        }
968     }
969     for (i = ly; i < uy; ++i)
970     {
971         im_put_rowf(arrN[i],im_new,i,lx,ux);
972     }
973     fvector_free(func,0);
974     fvector_free(val,0);
975     farray_free(arrN,ly,lx,uy,ux);
976     farray_free(arrY,ly,lx,uy,ux);
977     farray_free(arrX,ly,lx,uy,ux);
978 
979     return (im_new);
980 }
981 */
982 
983 double im_corscale(Imrect *im0, Imrect *corr, double noise)
984 /* routine to find the point of maximum difference to 1/x distribution */
985 /* using Bhattacharrya overap NAT 02/01/2002 */
986 {
987     Imregion roi;
988     int lx, ux, ly, uy, i, j, k, maxk;
989     float immax, immin;
990     float *err;         /* (cannot call 'error' since this conflicts with the error() function) */
991     float *rowimpix, *rowmod;
992     double diffx, impix, mod, fac, conts, range;
993     double scale = 2.6/CORSCAN, x;
994     double centre, h0, hm1, h1;
995     shistogram *hists[CORSCAN+1];
996 
997     if (noise <= 0.0) return (1.0);
998     err = fvector_alloc(0,CORSCAN);
999     imf_minmax(im0, &immin, &immax);
1000     range = 0.5*(immax-immin);
1001 
1002     for (k=0;k<CORSCAN;k++)
1003     {
1004         hists[k] = hbook1(k,"scale hist",immin-range, immax+range,
1005                           2.0*(immax-immin)/noise); 
1006         hreset(hists[k]);
1007     }
1008     hists[CORSCAN] = hbook1(CORSCAN, "performance", 0.0, 2.6, CORSCAN);
1009 
1010     roi = *(Imregion *)(im0->region); 
1011     lx = roi.lx; ux = roi.ux; ly = roi.ly; uy = roi.uy;
1012     rowimpix = fvector_alloc(lx, ux);
1013     rowmod = fvector_alloc(lx, ux);
1014     
1015     for (k = 0; k<CORSCAN; k++)
1016     {
1017         err[k] = 0.0;
1018     }
1019     for (j = ly+1; j< uy-1; j++)
1020     {
1021         im_get_rowf(rowimpix, im0, j, lx, ux);
1022         im_get_rowf(rowmod, corr, j, lx, ux);
1023         for (i=lx+1; i<ux-1; i++)
1024         {
1025              impix = rowimpix[i];
1026              mod = rowmod[i];
1027              for (k = 0; k<CORSCAN; k++)
1028              {
1029                  fac = (k+0.5)*scale; 
1030                  diffx = impix*exp(-fac*mod);
1031                  hfill1s(hists[k],diffx,1.0);
1032              }
1033         }
1034     } 
1035     fvector_free(rowimpix, lx);
1036     fvector_free(rowmod, lx);
1037 
1038     for (k = 0; k<CORSCAN; k++)
1039     {
1040         fac = (k+0.5)*scale; 
1041         for (j=0,x=hists[k]->xmin+0.5*hists[k]->xincr; j< hists[k]->xbins;
1042              j++,x+=hists[k]->xincr)
1043         {
1044             conts = hists[k]->array[0][j]; 
1045             err[k] -= sqrt(conts/(fabs(x)+0.5*noise));
1046 /*
1047    alternative standard information measure doesn't work due to lack
1048    of required invariances !!!
1049             err[k] += conts*log(1.0+conts);
1050 */
1051         }
1052         hfill1(hists[CORSCAN],fac,err[k]-err[0]);
1053     }
1054 // removed during promotion to tina-libs
1055 // a.lacey@man.ac.uk 23.05.03
1056 //    graph_hfit(imcalc_graph_tv(), hists[CORSCAN]);
1057 
1058 /* initialise to defaults for cases where no local maxima are found */
1059     maxk =0;
1060     for (k=1; k<CORSCAN;k++)
1061     {
1062         if (err[k]>err[0]) maxk = k;
1063     }
1064 
1065 /* now find nearest local maxima to fac= 1.0 */
1066     centre = 0.0;
1067     for (k = 0; k<(CORSCAN/2-1); k++)
1068     {
1069         if (err[CORSCAN/2-k]>err[CORSCAN/2-k-1]
1070           &&err[CORSCAN/2-k]>err[CORSCAN/2-k+1])
1071         {
1072             maxk = CORSCAN/2-k;
1073             h0  = err[maxk];
1074             h1  = err[maxk+1];
1075             hm1 = err[maxk-1];
1076             centre = (hm1-h1)/(2.0*(h1+hm1)-4.0*h0);
1077             break;
1078         }
1079         if (err[CORSCAN/2+k]>err[CORSCAN/2+k-1]
1080           &&err[CORSCAN/2+k]>err[CORSCAN/2+k+1])
1081         {
1082             maxk = CORSCAN/2+k;
1083             h0  = err[maxk];
1084             h1  = err[maxk+1];
1085             hm1 = err[maxk-1];
1086             centre = (hm1-h1)/(2.0*(h1+hm1)-4.0*h0);
1087             break;
1088         }
1089     }
1090 /*
1091     basic debug stuff
1092     hprint(stdout, hists[maxk]);
1093 */
1094     fac = (centre + maxk + 0.5)*scale;
1095     format(" scale factor %f = %f \n", fac, err[maxk]);
1096 
1097     for (k = 0; k<CORSCAN+1; k++) hfree(hists[k]);
1098     fvector_free(err,0);
1099     return (fac);
1100 }
1101 
1102 double im_corscale2(Imrect *im0, Imrect *corr, double noise)
1103 /* routine to determine best scaling factor of correction image based
1104    upon robust least squares of x and y image derivatives NAT 24/12/2001 */
1105 {
1106     Imrect *im1, *im2, *im3, *im4, *im5, *im6, *im7, *im8;
1107     Imrect *moddiffx, *moddiffy;
1108     Imregion *roi;
1109     int lx, ux, ly, uy, i, j, k, mink;
1110     double modx, mody, diffx, diffy, imfx, imfy, sumx, sumy, fac, 
1111            errfacx, errfacy;
1112     double thresh = STAT_LIM*noise*sqrt(4.0), weight, weightsum=0.0;
1113     double scale = 2.6/CORSCAN;
1114     shistogram *hists[CORSCAN+1];
1115     float *err;
1116 
1117     hists[CORSCAN] = hbook1(CORSCAN, "performance", 0.0, 2.6, CORSCAN);
1118     err = fvector_alloc(0,CORSCAN);
1119     im1 = im_bshift(im0,0,1);
1120     im2 = im_bshift(im0,1,0);
1121 
1122     im3 = im_sum(im1,im0);
1123     im4 = im_sum(im2,im0);
1124 
1125     im5 = im_diff(im0,im1);
1126     im6 = im_diff(im0,im2);
1127     im_free(im1); im_free(im2);
1128 
1129 
1130     im7 = im_bshift(corr,0,1);
1131     im8 = im_bshift(corr,1,0);
1132 
1133     moddiffx = im_diff(corr,im7);
1134     moddiffy = im_diff(corr,im8);
1135 
1136     roi = roi_inter(im0->region,corr->region); 
1137     lx = roi->lx; ux = roi->ux; ly = roi->ly; uy = roi->uy;
1138     
1139     for (k = 0; k<CORSCAN; k++)
1140     {
1141         err[k] = 0.0;
1142     }
1143     for (i=lx+1; i<ux-1; i++)
1144     {
1145         for (j = ly+1; j< uy-1; j++)
1146         {
1147              IM_PIX_GET(moddiffx, j, i, modx);
1148              IM_PIX_GET(moddiffy, j, i, mody);
1149              IM_PIX_GET(im3, j, i, sumx);         
1150              IM_PIX_GET(im4, j, i, sumy);         
1151              IM_PIX_GET(im5, j, i, imfx);
1152              IM_PIX_GET(im6, j, i, imfy);
1153              weight = imfx*imfx + imfy*imfy;
1154              weight = 1.0 - (sqrt(weight)/thresh);
1155              if (weight <= 0.0) weight = 0.0;
1156              else weight = 1.0;
1157              weightsum += weight;
1158 
1159              for (k = 0; k<CORSCAN; k++)
1160              {
1161                  fac = k*scale; 
1162                  diffx = (imfx - 0.5*sumx*fac*modx)/noise;
1163                  diffx *= diffx;
1164                  errfacx = 2.0*(1.0 + 0.25*fac*fac*modx*modx); 
1165                  if (diffx/errfacx >= 25.0) diffx = 25.0*errfacx;
1166 
1167                  diffy = (imfy - 0.5*sumy*fac*mody)/noise;
1168                  diffy *= diffy;
1169                  errfacy = 2.0*(1.0 + 0.25*fac*fac*mody*mody); 
1170                  if (diffy/errfacy >= 25.0) diffy = 25.0*errfacy;
1171 
1172                  err[k] += weight*(diffy/errfacy + diffx/errfacx); 
1173              }
1174         }
1175     } 
1176     mink = 0;
1177     for (k = 0; k<CORSCAN; k++)
1178     {
1179         if (err[k]<err[mink]) mink = k;
1180         hfill1(hists[CORSCAN],scale*(k+0.5),err[k]-err[0]);
1181     }
1182     
1183 // removed during promotion to tina-libs
1184 // a.lacey@man.ac.uk 23.05.03
1185 //    graph_hfit(imcalc_graph_tv(), hists[CORSCAN]);
1186 
1187     fac = mink*scale;
1188     format(" scale factor %f = %f \n", fac, 
1189           sqrt(err[mink]/2.0)/sqrt(weightsum));
1190     im_free(moddiffx);
1191     im_free(moddiffy);
1192     im_free(im3);
1193     im_free(im4);
1194     im_free(im5);
1195     im_free(im6);
1196     fvector_free(err,0);
1197     return (fac);
1198 }
1199 
1200 Imrect *im_integrate(Imrect *imdxf, Imrect *imdyf)
1201 {
1202     Imrect *im_new;
1203     Imregion *roi;
1204     
1205     if (imdxf == NULL || imdyf == NULL)
1206     {
1207         error("no images provided to imf_integrate",warning);
1208         return(NULL);
1209     }
1210 
1211     roi = imdxf->region;
1212 
1213     im_new = integrate2(roi,imdxf,imdyf);
1214 
1215     return(im_new);
1216 }
1217 
1218 
1219 void combine_corrs(Imrect **imd, Imrect **ims, Imrect *imd_r, Imrect *ims_r, Imrect *imd_i, Imrect *ims_i)
1220 {
1221    Imregion *roi=NULL;
1222    int i, j, lx, ly, ux, uy;
1223    float *row_d_r=NULL, *row_s_r=NULL, *row_d_i=NULL, *row_s_i=NULL;
1224    float *row_d=NULL, *row_s=NULL;
1225 
1226    roi = imd_r->region;
1227 
1228    lx = roi->lx;
1229    ux = roi->ux;
1230    ly = roi->ly;
1231    uy = roi->uy;
1232 
1233    *imd = im_alloc(imd_r->height, imd_r->width, roi, float_v);
1234    *ims = im_alloc(imd_r->height, imd_r->width, roi, float_v);
1235 
1236    row_d_r = fvector_alloc(lx, ux);
1237    row_s_r = fvector_alloc(lx, ux);
1238    row_d_i = fvector_alloc(lx, ux);
1239    row_s_i = fvector_alloc(lx, ux);
1240    row_d = fvector_alloc(lx, ux);
1241    row_s = fvector_alloc(lx, ux);
1242 
1243    for (i=ly; i<uy; i++)
1244    {
1245       im_get_rowf(row_d_r, imd_r, i, lx, ux);
1246       im_get_rowf(row_s_r, ims_r, i, lx, ux);
1247       im_get_rowf(row_d_i, imd_i, i, lx, ux);
1248       im_get_rowf(row_s_i, ims_i, i, lx, ux);
1249 
1250       for (j = lx; j < ux; ++j)
1251       {
1252          /* Have you got the errors the wrong way up? i.e. are they 1/s rather than s? */
1253 
1254 /*
1255          row_d[j] = (row_d_r[j]/(row_s_r[j]*row_s_r[j]) + row_d_i[j]/(row_s_i[j]*row_s_i[j]))/
1256                     (1/(row_s_r[j]*row_s_r[j]) + 1/(row_s_i[j]*row_s_i[j]));
1257          row_s[j] = sqrt(1/(1/(row_s_r[j]*row_s_r[j]) + 1/(row_s_i[j]*row_s_i[j])));
1258 */
1259          row_d[j] = (row_d_r[j]*(row_s_r[j]*row_s_r[j]) + row_d_i[j]*(row_s_i[j]*row_s_i[j]))/
1260                     (row_s_r[j]*row_s_r[j] + row_s_i[j]*row_s_i[j]);
1261          row_s[j] = sqrt(row_s_r[j]*row_s_r[j] + row_s_i[j]*row_s_i[j]);
1262       }
1263 
1264       im_put_rowf(row_d, *imd, i, lx, ux);
1265       im_put_rowf(row_s, *ims, i, lx, ux);
1266    }
1267 
1268    fvector_free(row_d, lx);
1269    fvector_free(row_s, lx);
1270    fvector_free(row_d_r, lx);
1271    fvector_free(row_s_r, lx);
1272    fvector_free(row_d_i, lx);
1273    fvector_free(row_s_i, lx);
1274 }
1275 
1276 
1277 Imrect *xy_normf(Imrect *im, double constant, double sigma, double thresh)
1278 {
1279    Imrect *im1,*im2,*im3,*imdx,*imdy;
1280    Imrect *imxs,*imys,*corr;
1281    double fac;
1282 
1283    im1 = im_square(im);
1284 
1285    smooth_slopes(constant,sigma,im1,&imdx,&imdy,&corr,&imxs,&imys);
1286    im_free(im1);
1287    enf_integ(&imdx,&imdy,corr,imxs,imys);
1288    im_free(corr);
1289    im2 = im_integrate(imdx,imdy);
1290 
1291    im_free(imxs); im_free(imys);
1292    im_free(imdx);
1293    im_free(imdy);
1294    if (im2 == NULL)
1295    {
1296       error("imcalc_xy_norm:  image not in calculator ",warning);
1297             return(NULL);
1298    }
1299 /* find maximum improvement in image histogram */
1300 
1301    fac = im_corscale2(im, im2, constant);
1302    im3 = imf_times(-fac,im2);
1303    im_free(im2);
1304    return(im3);
1305 }
1306 
1307 
1308 Imrect *xy_normz(Imrect *im, double constant, double sigma, double thresh)
1309 {
1310    Imrect *im_real=NULL, *im_cmplx=NULL, *im2=NULL, *im3=NULL, *im4=NULL, *corr=NULL;
1311    Imrect *im1_r=NULL, *imdx_r=NULL, *imdy_r=NULL, *imxs_r=NULL, *imys_r=NULL, *corr_r=NULL;
1312    Imrect *im1_i=NULL, *imdx_i=NULL, *imdy_i=NULL, *imxs_i=NULL, *imys_i=NULL, *corr_i=NULL;
1313    Imrect *imdx=NULL, *imdy=NULL, *imxs=NULL, *imys=NULL, *imdydx=NULL, *imdxdy=NULL;
1314 
1315    im_real = im_re(im);
1316    im_cmplx = im_im(im);
1317 
1318    im1_r = im_square(im_real);
1319    im1_i = im_square(im_cmplx);
1320    smooth_slopes(constant, sigma, im1_r, &imdx_r, &imdy_r, &corr_r, &imxs_r, &imys_r);
1321    smooth_slopes(constant, sigma, im1_i, &imdx_i, &imdy_i, &corr_i, &imxs_i, &imys_i);
1322    im_free(im1_r);
1323    im_free(im1_i);
1324    im_free(corr_r);
1325    im_free(corr_i);
1326 
1327    combine_corrs(&imdx, &imxs, imdx_r, imxs_r, imdx_i, imxs_i);
1328    combine_corrs(&imdy, &imys, imdy_r, imys_r, imdy_i, imys_i);
1329 
1330    /* re-est_corr */
1331 
1332    im3 = im_bshift(imdy,0,1);
1333    im4 = im_bshift(imdx,1,0);
1334    imdydx = im_diff(imdy, im3);
1335    imdxdy = im_diff(imdx, im4);
1336    im_free(im3); 
1337    im_free(im4);
1338    corr = im_diff(imdxdy,imdydx);
1339    im_free(imdxdy); 
1340    im_free(imdydx);
1341 
1342    enf_integ(&imdx, &imdy, corr, imxs, imys);
1343    im_free(corr);
1344 
1345    im2 = im_integrate(imdx,imdy);
1346 
1347    im_free(imxs);
1348    im_free(imys);
1349    im_free(imdx);
1350    im_free(imdy);
1351    im_free(imxs_r);
1352    im_free(imys_r);
1353    im_free(imdx_r);
1354    im_free(imdy_r);
1355    im_free(imxs_i);
1356    im_free(imys_i);
1357    im_free(imdx_i);
1358    im_free(imdy_i);
1359 
1360    if (im2 == NULL)
1361    {
1362       error("imcalc_xy_norm:  image not in calculator ",warning);
1363             return(NULL);
1364    }
1365 
1366    im3 = imf_times(-1.0,im2);
1367 
1368 
1369    im_free(im_real);
1370    im_free(im_cmplx);
1371    return(im3);
1372 }
1373 
1374 
1375 Imrect *xy_norm(Imrect *im, double constant, double sigma, double thresh)
1376 {
1377      if(im == NULL) return(NULL);
1378      switch(im->vtype)
1379      {
1380          case uchar_v:
1381          case char_v:
1382          case short_v:
1383          case ushort_v:
1384          case int_v:
1385          case float_v:
1386              return(xy_normf(im, constant, sigma, thresh));
1387          case complex_v:
1388              return(xy_normz(im, constant, sigma, thresh));
1389          default:
1390              return(NULL);
1391      }
1392 }
1393 

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