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