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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.