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/vision/visCorr_warp.c,v $
27 * Date : $Date: 2004/08/04 15:07:47 $
28 * Version : $Revision: 1.1 $
29 * CVS Id : $Id: visCorr_warp.c,v 1.1 2004/08/04 15:07:47 paul Exp $
30 *
31 * Author : Simon Crossley (see thesis Univ. Sheffield U.K.), NAT.
32 */
33 /**
34 * @file stereo correlation algorithms
35 * @brief Routines to compute warp correlation.
36 *
37 *********
38 */
39
40 #include "visCorr_warp.h"
41
42 #include <math.h>
43 #include <float.h>
44 #include <limits.h>
45 #include <tina/sys/sysDef.h>
46 #include <tina/sys/sysPro.h>
47 #include <tina/math/mathDef.h>
48 #include <tina/math/mathPro.h>
49 #include <tina/geometry/geomDef.h>
50 #include <tina/geometry/geomPro.h>
51 #include <tina/image/imgDef.h>
52 #include <tina/image/imgPro.h>
53
54 #include "vis_CorrDef.h"
55 #include "vis_CorrPro.h"
56
57 static void imf_fast_warp_correlate(Imrect *im1, Imrect *im2, Imrect *ctf_im,
58 List ***mgrid, int brows, int bcols,
59 int ctf_scale,
60 int stretch, int compress, int width, int height,
61 int range, float mthresh, Bool min_disp_cut)
62 {
63 int lx1, ly1, ux1, uy1, lx2, ly2, ux2, uy2;
64 int i, j, m, n;
65 Imregion roi1, roi2, disp_lim;
66 int pos, disp_range;
67 float **scores, **cross_partial, *auto_partial2, *adj_partial2;
68 float norm1, **norm2, **norm2_realigned, **alpha, **beta;
69 int **kappa;
70 int default_disp, norm2_calced=0;
71 Ctf_info *info_ptr;
72
73 lx1 = (im1->region->lx + PAD_DEPTH) / ctf_scale;
74 ly1 = (im1->region->ly + PAD_DEPTH) / ctf_scale;
75 ux1 = (im1->region->ux - PAD_DEPTH) / ctf_scale;
76 uy1 = (im1->region->uy - PAD_DEPTH) / ctf_scale;
77 lx2 = (im2->region->lx + PAD_DEPTH) / ctf_scale;
78 ly2 = (im2->region->ly + PAD_DEPTH) / ctf_scale;
79 ux2 = (im2->region->ux - PAD_DEPTH) / ctf_scale;
80 uy2 = (im2->region->uy - PAD_DEPTH) / ctf_scale;
81
82 default_disp = (lx2 + ux2 - lx1 - ux1)/2;
83
84 alpha = farray_alloc(-compress, -width/2, stretch+1, width/2);
85 beta = farray_alloc(-compress, -width/2, stretch+1, width/2);
86 kappa = iarray_alloc(-compress, -width/2, stretch+1, width/2);
87 setup_corLUTs(alpha, beta, kappa, compress, stretch, width);
88
89 auto_partial2 = fvector_alloc(lx2, ux2);
90 adj_partial2 = fvector_alloc(lx2, ux2);
91 norm2 = farray_alloc(lx2+width/2+stretch, -compress,
92 ux2-width/2-stretch+1, stretch+1);
93
94 for (i = ly1+height/2, m=0; i <= uy1-height/2; i+=height, m+=ctf_scale)
95 {
96 /* left and right vertical roi settings */
97 roi1.ly = i - height/2;
98 roi1.uy = i + height/2;
99 roi2.ly = roi1.ly;
100 roi2.uy = roi1.uy;
101
102 /* check roi2 is inside im2 vertically */
103 if ((roi2.uy > uy2) || (roi2.ly < ly2)) continue;
104
105 norm2_calced = 0;
106
107 for (j=lx1+width/2, n=0; j <= ux1-width/2; j+=width, n+=ctf_scale)
108 {
109 roi1.lx = j - width/2;
110 roi1.ux = j + width/2;
111
112 info_ptr = (Ctf_info *)IM_PTR(ctf_im, m, n);
113
114 if (info_ptr==NULL || info_ptr->ctf_scale!=ctf_scale ||
115 info_ptr->edge_count < 1)
116 continue;
117
118 disp_range = range/2;
119 pos = default_disp;
120
121 if (info_ptr->max_disp!=-FLT_MAX)
122 {
123 disp_range = range/2 + ((int)info_ptr->max_disp -
124 (int)info_ptr->min_disp)/2;
125 pos = (int)(info_ptr->max_disp + info_ptr->min_disp)/2;
126 }
127
128 if (min_disp_cut)
129 {
130 int min_disp, max_disp;
131
132 min_disp = pos - disp_range;
133 max_disp = pos + disp_range;
134
135 /* only allow disparities for objects in
136 front of the camera */
137 if (default_disp > 0)
138 /* positive disparities in right to left matching */
139 min_disp = MAX(min_disp, 0);
140 else
141 /* negative disparities in left to right matching */
142 max_disp = MIN(max_disp, 0);
143
144 disp_range = (max_disp - min_disp) / 2;
145 pos = (max_disp + min_disp) / 2;
146 }
147
148 scores = farray_alloc(-disp_range, -compress, disp_range+1, stretch+1);
149 cross_partial = farray_alloc(-width/2, -disp_range-width/2-stretch,
150 width/2, disp_range+width/2+stretch+1);
151
152 if (!norm2_calced)
153 {
154 roi2.lx = lx2;
155 roi2.ux = ux2;
156 norm2_calced = 1;
157 norm2 = auto_correlate(im2, &roi2, width, compress, stretch,
158 auto_partial2, adj_partial2, norm2,
159 kappa, alpha, beta, ctf_scale);
160 }
161
162 /* keep the disparity search within im2 */
163 disp_lim.lx = MAX(-disp_range, lx2-(j+pos)+width/2+stretch);
164 disp_lim.ux = MIN(disp_range, ux2-(j+pos)-width/2-stretch);
165
166 roi2.lx = (j+pos) + disp_lim.lx - width/2 - stretch;
167 roi2.ux = (j+pos) + disp_lim.ux + width/2 + stretch;
168 norm2_realigned = &norm2[j+pos];
169
170 cross_partial = partial_sum(im1, im2, &roi1, &roi2, &disp_lim,
171 width, stretch, compress,
172 cross_partial, &norm1, ctf_scale);
173
174 scores = partial_accumulate(cross_partial, kappa, alpha, beta,
175 stretch, compress, width, &disp_lim,
176 scores, norm1, norm2_realigned);
177
178 mgrid[m][n] = add_match_branch(scores, mthresh, disp_lim.lx,
179 disp_lim.ux, compress, stretch, pos);
180
181 farray_free(scores,-disp_range, -compress, disp_range+1, stretch+1);
182 farray_free(cross_partial,-width/2, -disp_range-width/2-stretch ,
183 width/2, disp_range+width/2+stretch+1);
184 }
185 }
186
187 farray_free(alpha, -compress, -width/2, stretch+1, width/2);
188 farray_free(beta, -compress, -width/2, stretch+1, width/2);
189 iarray_free(kappa, -compress, -width/2, stretch+1, width/2);
190 fvector_free(auto_partial2, lx2);
191 fvector_free(adj_partial2, lx2);
192 farray_free(norm2, lx2+width/2+stretch, -compress, ux2-width/2-stretch+1,
193 stretch+1);
194 }
195
196 static void gen_disp_ims(Imrect *im1, Imrect *edge_im1, Imrect *edge_im2,
197 Imrect *ctf_im1, Imrect **out_disp_im_ptr,
198 List ***mgrid, int blk_width, int blk_height,
199 int brows, int bcols, double Edelta)
200 {
201 List *lptr;
202 Stereo_match *match;
203 int m, n;
204 Ctf_info *info_ptr;
205 int disparity, warp;
206 Imregion roi1, roi2;
207 double edge_delta;
208
209 #define out_disp_im (*out_disp_im_ptr)
210
211 out_disp_im = im_alloc(im1->height, im1->width, im1->region, float_v);
212 FOR_IM (out_disp_im->region, m, n)
213 IM_FLOAT(out_disp_im, m, n) = NO_DISP;
214
215 for (m=0; m<brows; m++)
216 for (n=0; n<bcols; n++)
217 {
218 if ((info_ptr = (Ctf_info *)IM_PTR(ctf_im1, m, n)) == NULL)
219 continue;
220
221 if ((lptr = mgrid[m][n]))
222 {
223 match = lptr->to;
224 disparity = (match->disparity + match->pos) * info_ptr->ctf_scale;
225 warp = match->warp * info_ptr->ctf_scale;
226 edge_delta = Edelta * (double)(info_ptr->ctf_scale);
227
228 roi1 = info_ptr->pix_roi;
229 roi2.ly = roi1.ly;
230 roi2.uy = roi1.uy;
231 roi2.lx = roi1.lx + disparity - warp;
232 roi2.ux = roi1.ux + disparity + warp;
233 stretch_output_block(out_disp_im, edge_im1, edge_im2,
234 &roi1, &roi2, edge_delta);
235 }
236 }
237
238 #undef out_disp_im
239 }
240
241 #ifdef __STORE_CTF__
242 static void push_ctf_status(Imrect *ctf_im, Imrect *grey_im,
243 Imrect *disp_im)
244 {
245 Imrect *status_im;
246 Ctf_info *info_ptr;
247 int x, y, m, n;
248
249 status_im = im_alloc(grey_im->height, grey_im->width, grey_im->region, char_v);
250
251 /* mark all areas as 'unprocessed' */
252 FOR_IM (status_im->region, y, x)
253 IM_CHAR(status_im, y, x) = -50;
254
255 /* fill in the ctf level areas */
256 FOR_IM (ctf_im->region, y, x)
257 if (info_ptr = (Ctf_info *)IM_PTR(ctf_im, y, x))
258 {
259 int lx, ly, ux, uy;
260
261 FOR_IM (&(info_ptr->pix_roi), m, n)
262 IM_CHAR(status_im, m, n) = 10 * (char)info_ptr->ctf_level;
263
264 /* put a 1 pixel thick border round each block */
265
266 lx = info_ptr->pix_roi.lx;
267 ux = info_ptr->pix_roi.ux - 1;
268 ly = info_ptr->pix_roi.ly;
269 uy = info_ptr->pix_roi.uy - 1;
270
271 for (m=ly, n=lx; m<=uy; m++)
272 IM_CHAR(status_im, m, n) = -50;
273 for (m=ly, n=ux; m<=uy; m++)
274 IM_CHAR(status_im, m, n) = -50;
275 for (m=ly, n=lx; n<=ux; n++)
276 IM_CHAR(status_im, m, n) = -50;
277 for (m=uy, n=lx; n<=ux; n++)
278 IM_CHAR(status_im, m, n) = -50;
279 }
280
281 /* add in the pickup disparity locations */
282 if (disp_im != NULL)
283 FOR_IM (disp_im->region, y, x)
284 if (IM_FLOAT(disp_im, y, x) != NO_DISP)
285 IM_CHAR(status_im, y, x) = 50;
286
287 stack_push((void *)status_im, IMRECT, im_free);
288 }
289 #endif
290
291 void imf_warp_correlate(Imrect *left_im, Imrect *right_im,
292 Imrect *left_edge_im, Imrect *right_edge_im,
293 Imrect **out_disp_im_ptr, Imrect *in_disp_im,
294
295 int stretch, int compress, int width, int height,
296 int range, double Edelta, float mthresh,
297
298 Bool global_DG, double dg_limit, double dg_pass,
299 int dg_range,
300
301 int pickup_range, int max_ctf_level,
302 Bool min_disp_cut, Bool fix_ctf_max,
303 Bool lr_only)
304 {
305 static Imrect *left_ctf_im = NULL, *right_ctf_im = NULL;
306 Imrect *left_disp_im = NULL, *right_disp_im = NULL;
307 Imrect *in_right_disp_im = NULL;
308 List ***mgrid, *lptr;
309 int ctf, left_brows, left_bcols, right_brows, right_bcols;
310 Imregion roi;
311 int x, y, ctf_scale;
312
313 #define out_disp_im (*out_disp_im_ptr)
314
315 if (in_disp_im &&
316 (in_disp_im->region->lx!=left_im->region->lx ||
317 in_disp_im->region->ux!=left_im->region->ux ||
318 in_disp_im->region->ly!=left_im->region->ly ||
319 in_disp_im->region->uy!=left_im->region->uy))
320 {
321 error("imf_warp_correlate: image roi(s) changed\nignoring previous bootstrap", warning);
322 in_disp_im = NULL;
323 }
324
325 if (lr_only == false)
326 {
327 /* run the correlation matcher using right to left matching */
328
329 in_right_disp_im = imf_disp_conv(in_disp_im, right_im);
330
331 roi = *(right_im->region);
332 roi.lx += PAD_DEPTH;
333 roi.ly += PAD_DEPTH;
334 roi.ux -= PAD_DEPTH;
335 roi.uy -= PAD_DEPTH;
336 right_bcols = (roi.ux - roi.lx) / width;
337 right_brows = (roi.uy - roi.ly) / height;
338
339 mgrid = (List ***)parray_alloc(-1, -1, right_brows+1, right_bcols+1);
340
341 if (right_ctf_im == NULL || in_disp_im==NULL)
342 {
343 ctf_im_free(right_ctf_im);
344 right_ctf_im = ctf_im_alloc(right_brows, right_bcols);
345 }
346
347 if (in_disp_im == NULL)
348 /* algorithm startup */
349 init_ctf_im(right_ctf_im, right_im, max_ctf_level, width, height);
350 else
351 update_ctf_im(right_ctf_im, right_im, in_right_disp_im, pickup_range,
352 max_ctf_level, width, height);
353
354 update_edge_info(right_ctf_im, right_edge_im);
355
356 for (ctf=0; ctf<=max_ctf_level; ctf++)
357 {
358 int disp_range = range;
359
360 ctf_scale = 1 << ctf;
361
362 if (fix_ctf_max && max_ctf_level != 0)
363 {
364 /* ensure at coarsest scale, search range is 50% im width
365 and intermediate ranges grow exponentially in size */
366
367 double dw, dr, dc, n, r;
368
369 dw = (double)(left_im->region->ux - left_im->region->lx);
370 dr = 0.5 / (double)range;
371 dc = 1.0 / (double)max_ctf_level;
372
373 /* calculate the exponential factor n */
374 n = pow(dw * dr, dc);
375 /* calculate the full scale range parameter
376 for the current scale */
377 r = (double)range * pow(n, (double)ctf);
378 /* take into account change in image scale */
379 r /= (double)ctf_scale;
380
381 disp_range = (int)r;
382 }
383
384 imf_fast_warp_correlate(right_im, left_im, right_ctf_im,
385 mgrid, right_brows, right_bcols,
386 ctf_scale,
387 stretch, compress, width, height,
388 disp_range, mthresh, min_disp_cut);
389 }
390
391 if (global_DG)
392 {
393 apply_global(right_ctf_im, mgrid, dg_limit,
394 dg_pass, dg_range);
395 }
396
397 gen_disp_ims(right_im, right_edge_im, left_edge_im, right_ctf_im,
398 &right_disp_im, mgrid, width, height,
399 right_brows, right_bcols, Edelta);
400
401 roi.lx = -1;
402 roi.ly = -1;
403 roi.ux = right_bcols+1;
404 roi.uy = right_brows+1;
405
406 FOR_IM(&roi, y, x)
407 {
408 if ((lptr = mgrid[y][x]))
409 list_rm(lptr, rfree);
410 }
411 parray_free(mgrid, -1, -1, right_brows+1, right_bcols+1);
412
413 im_free(in_right_disp_im);
414 }
415
416 /* run the correlation matcher using left to right matching */
417
418 roi = *(left_im->region);
419 roi.lx += PAD_DEPTH;
420 roi.ly += PAD_DEPTH;
421 roi.ux -= PAD_DEPTH;
422 roi.uy -= PAD_DEPTH;
423 left_bcols = (roi.ux - roi.lx) / width;
424 left_brows = (roi.uy - roi.ly) / height;
425
426 mgrid = (List ***)parray_alloc(-1, -1, left_brows+1, left_bcols+1);
427
428 if (left_ctf_im==NULL || in_disp_im==NULL)
429 {
430 ctf_im_free(left_ctf_im);
431 left_ctf_im = ctf_im_alloc(left_brows, left_bcols);
432 }
433
434 if (in_disp_im == NULL)
435 /* algorithm startup */
436 init_ctf_im(left_ctf_im, left_im, max_ctf_level, width, height);
437 else
438 update_ctf_im(left_ctf_im, left_im, in_disp_im, pickup_range,
439 max_ctf_level, width, height);
440
441 update_edge_info(left_ctf_im, left_edge_im);
442
443 for (ctf=0; ctf<=max_ctf_level; ctf++)
444 {
445 int disp_range = range;
446
447 ctf_scale = 1 << ctf;
448
449 if (fix_ctf_max && max_ctf_level != 0)
450 {
451 /* ensure at coarsest scale, search range is 50% im width
452 and intermediate ranges grow exponentially in size */
453
454 double dw, dr, dc, n, r;
455
456 dw = (double)(left_im->region->ux - left_im->region->lx);
457 dr = 0.5 / (double)range;
458 dc = 1.0 / (double)max_ctf_level;
459
460 /* calculate the exponential factor n */
461 n = pow(dw * dr, dc);
462 /* calculate the full scale range parameter
463 for the current scale */
464 r = (double)range * pow(n, (double)ctf);
465 /* take into account change in image scale */
466 r /= (double)ctf_scale;
467
468 disp_range = (int)r;
469 }
470
471 imf_fast_warp_correlate(left_im, right_im, left_ctf_im,
472 mgrid, left_brows, left_bcols,
473 ctf_scale,
474 stretch, compress, width, height,
475 disp_range, mthresh, min_disp_cut);
476 }
477
478 if (global_DG)
479 {
480 apply_global(left_ctf_im, mgrid, dg_limit,
481 dg_pass, dg_range);
482 }
483
484 gen_disp_ims(left_im, left_edge_im, right_edge_im, left_ctf_im,
485 &left_disp_im, mgrid, width, height,
486 left_brows, left_bcols, Edelta);
487
488 roi.lx = -1;
489 roi.ly = -1;
490 roi.ux = left_bcols+1;
491 roi.uy = left_brows+1;
492
493 FOR_IM(&roi, y, x)
494 {
495 if ((lptr = mgrid[y][x]))
496 list_rm(lptr, rfree);
497 }
498 parray_free(mgrid, -1, -1, left_brows+1, left_bcols+1);
499
500 if (lr_only == false)
501 {
502 /* merge the two way matches to ensure uniqueness */
503 out_disp_im = edge_im_disp_merge(left_disp_im, left_edge_im, right_edge_im);
504 }
505 else
506 {
507 out_disp_im = im_copy(left_disp_im);
508 }
509
510 #ifdef __STORE_CTF__
511 push_ctf_status(left_ctf_im, left_im, out_disp_im);
512 #endif
513
514 im_free(left_disp_im);
515 im_free(right_disp_im);
516
517 return;
518
519 #undef out_disp_im
520 }
521
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.