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

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

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