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

Linux Cross Reference
Tina5/tina-libs/tina/vision/visCorr_warp.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/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 

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