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

Linux Cross Reference
Tina6/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 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 

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