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

Linux Cross Reference
Tina5/tina-libs/tina/vision/visCorr_correlate.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_correlate.c,v $
 23  * Date    :  $Date: 2008/10/30 09:43:59 $
 24  * Version :  $Revision: 1.2 $
 25  * CVS Id  :  $Id: visCorr_correlate.c,v 1.2 2008/10/30 09:43:59 neil Exp $
 26  *
 27  * Author : S. Crossley, R. Lane (see thesis Univ. Sheffield U.K.), NAT.
 28  */
 29 /**
 30  * @file stereo correlation algorithms
 31  * @brief Stretch correlation algorithms (see Richard Lanes's Thesis Univ. Sheffiled).
 32  *
 33  *********
 34 */
 35 
 36 #include "visCorr_correlate.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 #include <tina/vision/visDef.h>
 50 #include <tina/vision/visPro.h>
 51 
 52 #include <tina/vision/vis_CorrDef.h>
 53 #include <tina/vision/vis_CorrPro.h>
 54 
 55 /*
 56 #define  VERT  atan2(1.0, 0.0)
 57 */
 58 #define  ANGERR 0.85
 59 static Imregion im_roi;
 60 static int count, pi, pj;
 61 
 62 void setup_corLUTs(float **alpha, float **beta, int **kappa, int compress,
 63                    int stretch, int width)
 64 {
 65   int       i, j;
 66   float     sub_pix, step;
 67 
 68   for (i=-compress; i<stretch+1; i++)
 69     {
 70       sub_pix = -((float)(width/2)) - (float)i;
 71       step = 1.0 + 2.0 * (float)i / ((float)width);
 72       for (j=-width/2; j<width/2; j++)
 73         {
 74           kappa[i][j] = tina_int(sub_pix);
 75           alpha[i][j] = 1.0 - (sub_pix - (float)kappa[i][j]);
 76           beta [i][j] = 1.0 - alpha[i][j];
 77           sub_pix += step;
 78         }
 79     }
 80 }
 81 
 82 static float IM_CTF_FLOAT(Imrect *im, int y, int x, int ctf_scale)
 83 {
 84   int i, j;
 85   float sum = 0;
 86 
 87   y *= ctf_scale;
 88   x *= ctf_scale;
 89 
 90   for (i = y; i < y + ctf_scale; i++)
 91     for (j = x; j < x + ctf_scale; j++)
 92       sum += IM_FLOAT(im, i, j);
 93 
 94   return(sum / (float)(ctf_scale * ctf_scale));
 95 }
 96 
 97 float **auto_correlate(Imrect *im, Imregion *roi, int width, int compress,
 98                        int stretch, float *auto_partial, float *adj_partial,
 99                        float **norm, int **kappa, float **alpha, float **beta,
100                        int ctf_scale)
101 {
102   int          lx, ux, ly, uy;
103   int          i, j, k, l, width_over2, j_plus1;
104   int          stretch_plus1, upper_lim;
105   float        sum, adj_sum, pix, a, b;
106   float       *col, *adj_col, *swap_temp;
107   int         *kappa_row;
108   float       *alpha_row, *beta_row;
109 
110   lx = roi->lx;
111   ly = roi->ly;
112   ux = roi->ux;
113   uy = roi->uy;
114   width_over2 = width/2;
115   upper_lim = ux-stretch-width_over2 + 1;
116   stretch_plus1 = stretch+1;
117 
118   col = fvector_alloc(ly, uy);
119   adj_col = fvector_alloc(ly, uy);
120 
121   for (k=ly; k<uy; k++)
122     {
123       adj_col[k] = IM_CTF_FLOAT(im, k, lx, ctf_scale);
124     }
125   for (j=lx; j<ux; j++)
126     {
127       swap_temp = col;
128       col = adj_col;
129       adj_col = swap_temp;
130       j_plus1 = j+1;
131 
132       sum = adj_sum = 0;
133       for (k=ly; k<uy; k++)
134         {
135           pix = col[k];
136           sum += pix * pix;
137           adj_sum += pix * (adj_col[k] = IM_CTF_FLOAT(im, k, j_plus1, ctf_scale));
138         }
139       auto_partial[j] = sum;
140       adj_partial[j] = adj_sum;
141     }
142 
143   for (l=lx+stretch+width_over2; l<upper_lim; l++)
144     {
145       for (i=-compress; i<stretch_plus1; i++)
146         {
147           kappa_row = kappa[i];
148           alpha_row = alpha[i];
149           beta_row = beta[i];
150           sum = 0;
151           for (j=-width_over2; j<width_over2; j++)
152             {
153               k = kappa_row[j];
154               a = alpha_row[j];
155               b = beta_row[j];
156 
157               sum += a*a*auto_partial[l+k] + 2.0*a*b*adj_partial[l+k] +
158                      b*b*auto_partial[l+k+1];
159             }
160           norm[l][i] = sum;
161         }
162     }
163   return(norm);
164 }
165 
166 float **partial_sum(Imrect *im1, Imrect *im2, Imregion *roi1, Imregion *roi2,
167                     Imregion *disp, int width, int stretch, int compress,
168                     float **cross_partial, float *norm1, int ctf_scale)
169 {
170   int          lx1, ux1, ly, uy;
171   int          lx2, ux2;
172   int          i, j, l, d, row_offset;
173   int          disp_min;
174   float        sum, auto1_sum, pix;
175   float       *col;
176 
177   lx1 = roi1->lx;
178   ly  = roi1->ly;
179   ux1 = roi1->ux;
180   uy  = roi1->uy;
181   lx2 = roi2->lx;
182   ux2 = roi2->ux;
183 
184   disp_min = disp->lx;
185   row_offset = lx1 + width/2;
186   col = fvector_alloc(ly, uy);
187 
188   auto1_sum = 0;
189   for (j=lx1; j<ux1; j++)
190     {
191       for (i=ly; i<uy; i++)
192         {
193           pix = col[i] = IM_CTF_FLOAT(im1, i, j, ctf_scale);
194           auto1_sum += pix*pix;
195         }
196 
197       for (l=lx2, d=disp_min-width/2-stretch; l<ux2; l++, d++)
198         {
199           sum = 0;
200           for (i=ly; i<uy; i++)
201             {
202               sum += col[i] * IM_CTF_FLOAT(im2, i, l, ctf_scale);
203             }
204           cross_partial[j-row_offset][d] = sum;
205         }
206     }
207   fvector_free(col, ly);
208 
209   *norm1 = auto1_sum;
210   return(cross_partial);
211 }
212 
213 float **partial_accumulate(float **cross_partial, int **kappa, float **alpha,
214                            float **beta, int stretch, int compress, int width,
215                            Imregion *disp_lim, float **scores,
216                            float norm1, float **norm2)
217 {
218   int       i, j, d, disp_max, disp_min;
219   int       width_over2, stretch_plus1;
220   int      *kappa_row, k;
221   float    *alpha_row, *beta_row, *scores_row, *norm2_row, sum;
222 
223   disp_min = disp_lim->lx;
224   disp_max = disp_lim->ux;
225   stretch_plus1 = stretch+1;
226   width_over2 = width/2;
227 
228   for (d=disp_min; d<=disp_max; d++)
229     {
230       for (i=-compress; i<stretch_plus1; i++)
231         {
232           kappa_row = kappa[i];
233           alpha_row = alpha[i];
234           beta_row = beta[i];
235           sum = 0;
236           for (j=-width_over2; j<width_over2; j++)
237             {
238               k = kappa_row[j];
239               sum += alpha_row[j] * cross_partial[j][d+k] +
240                      beta_row[j] * cross_partial[j][d+k+1];
241             }
242           scores[d][i] = sum;
243         }
244     }
245 
246   for (d=disp_min; d<=disp_max; d++)
247     {
248       scores_row = scores[d];
249       norm2_row = norm2[d];
250       for (j=-compress; j<stretch_plus1; j++)
251         {
252           scores_row[j] /= sqrt(norm1 * norm2_row[j]);
253         }
254     }
255 
256   return(scores);
257 }
258 
259 int im_count_edges(Imrect *im, Imregion *roi)
260 {
261   int         lx, ly, ux, uy, i, j;
262   Imregion   *valid_region=NULL;
263   int         count=0;
264 
265   valid_region = roi_inter(roi,im->region);
266 
267   if (valid_region == NULL)
268   {
269       error("No valid region found \n", non_fatal);
270       return(0);
271   } 
272   lx = valid_region->lx;
273   ly = valid_region->ly;
274   ux = valid_region->ux;
275   uy = valid_region->uy;
276 
277   for (i = ly; i < uy; i++)
278     for (j = lx; j < ux; j++)
279       {
280         if (IM_PTR(im, i, j))
281           count++;
282       }
283 
284   rfree(valid_region);
285   return(count);
286 }
287 
288 
289 static void update_maxmin(Imrect *im, float *max, float *min)
290 {
291   float pix;
292 
293   if (pi>=im_roi.ly && pi<im_roi.uy && pj>=im_roi.lx && pj<im_roi.ux)
294     {
295       pix = IM_FLOAT(im, pi, pj);
296       if (pix != NO_DISP)
297         {
298           if (pix > *max) *max = pix;
299           if (pix < *min) *min = pix;
300           count++;
301         }
302     }
303 }
304 
305 int imf_nz_maxmin(Imrect *im, float *max, float *min, int cen_i, int cen_j,
306                   int srch_range)
307 {
308   int range;
309   float pix;
310 
311   *max = -FLT_MAX;
312   *min = FLT_MAX;
313   count = 0;
314 
315   im_roi.lx = im->region->lx;
316   im_roi.ux = im->region->ux;
317   im_roi.ly = im->region->ly;
318   im_roi.uy = im->region->uy;
319 
320   pix = IM_FLOAT(im, cen_i, cen_j);
321   if (pix != NO_DISP)
322     {
323       if (pix > *max) *max = pix;
324       if (pix < *min) *min = pix;
325       count++;
326     }
327 
328   for (range = 1; range <= srch_range; range++)
329     {
330       for (pi=cen_i-range, pj=cen_j-range; pi<=cen_i+range; pi++)
331         update_maxmin(im, max, min);
332       for (pi=cen_i-range, pj=cen_j+range; pi<=cen_i+range; pi++)
333         update_maxmin(im, max, min);
334       for (pi=cen_i-range, pj=cen_j-range+1; pj<cen_j+range; pj++)
335         update_maxmin(im, max, min);
336       for (pi=cen_i+range, pj=cen_j-range+1; pj<cen_j+range; pj++)
337         update_maxmin(im, max, min);
338     }
339 
340   return(count);
341 }
342 
343 static double edgel_realign(Edgel *algn_pix, Edgel *ref_pix)
344 {
345   double    deltay;
346   double    pos;
347 
348   if (ref_pix == NULL)
349     return(vec2_x(algn_pix->pos));
350 
351   deltay = vec2_y(ref_pix->pos) - vec2_y(algn_pix->pos);
352   pos = deltay/tan(algn_pix->orient);
353 
354   return(vec2_x(algn_pix->pos)-pos);
355 }
356 
357 float imf_hnonzero(Imrect *im, int i, int j, int range)
358 {
359   int         lx, ux, ly, uy, l, col;
360   double      pix;
361 
362   lx = im->region->lx;
363   ux = im->region->ux;
364   ly = im->region->ly;
365   uy = im->region->uy;
366   if (i < ly || i >= uy || j < lx || j >= ux)
367     return(NO_DISP);
368 
369   pix = IM_FLOAT(im, i, j);
370   if (pix != NO_DISP)
371     return(pix);
372 
373   for (l=1; l<= range; l++)
374     {
375       col = j-l;
376       if (col >= lx)
377         {
378           pix = IM_FLOAT(im, i, col);
379           if (pix != NO_DISP)
380             return(pix);
381         }
382       
383       col = j+l;
384       if (col < ux)
385         {
386           pix = IM_FLOAT(im, i, col);
387           if (pix != NO_DISP)
388             return(pix);
389         }
390     }
391   return(NO_DISP);
392 }
393 
394 static Edgel *edge_im_hnonzero(Imrect *edge_im, int y, int x, double maxdist)
395 {
396   Edgel      *edge = NULL, *best = NULL;
397   double      disp, best_pos = FLT_MAX;
398   int         i, lx, ux, ly, uy;
399 
400   if (edge_im->vtype != ptr_v)
401     return(NULL);
402 
403   lx = edge_im->region->lx;
404   ux = edge_im->region->ux;
405   ly = edge_im->region->ly;
406   uy = edge_im->region->uy;
407   if (x < lx || x >= ux || y < ly || y >= uy)
408     return(NULL);
409 
410   for (i = x - tina_rint(maxdist); i <= x + tina_rint(maxdist); i++)
411     {
412       if (i < lx || i >= ux || y < ly || y >= uy)
413         continue;
414       if ((edge = IM_PTR(edge_im, y, i)) == NULL)
415         continue;
416 
417       disp = fabs(vec2_x(edge->pos)-(double)x);
418       if ((disp < best_pos) && (disp < maxdist || !maxdist))
419         {
420           best = edge;
421           best_pos = disp;
422         }
423     }
424 
425   return(best);
426 }
427 
428 void stretch_output_block(Imrect *out_disp_im, Imrect *im1_edges,
429                           Imrect *im2_edges, Imregion *roi1,
430                           Imregion *roi2, double edge_delta)
431 {
432    int       i, j1, j2, k, l;
433    int       height, width1, width2;
434    float     roi2_pos;
435    float     estimate;
436    double    pos1, pos2;
437    double    orr_diff, xscale;
438    Edgel    *edge1, *edge2;
439    Match    *match;
440 /*
441    shistogram **hists;
442    hists = (shistogram **)hist_vec();
443 */
444 
445 
446    height = roi1->uy - roi1->ly;
447    width1 = roi1->ux - roi1->lx;
448    width2 = roi2->ux - roi2->lx;
449    i = roi1->ly + (height / 2);
450    j1 = roi1->lx + (width1 / 2);
451    j2 = roi2->lx + (width2 / 2);
452    xscale = (float)width2 / (float)width1;
453 
454    for (k = -height/2; k < height/2; k++)
455    {
456       for (l = -width1/2; l < width1/2; l++)
457       {
458         /* calculate the edge disparity */
459          if ((edge1 = edge_im_hnonzero(im1_edges, i+k, j1+l, 0.0)))
460          {
461              pos1 = vec2_x(edge1->pos) - (double)j1;
462              roi2_pos = (float)j2 + ((float)pos1 * xscale);
463 
464              if ((edge2 = edge_im_hnonzero(im2_edges, i+k, (int)roi2_pos, edge_delta)))
465              {
466                  orr_diff = edge1->orient - atan(xscale*tan(edge2->orient));
467                  if (orr_diff>PI) orr_diff -= PI;
468                  if (orr_diff<-PI) orr_diff += PI;
469                  if (orr_diff> PI/2.0) orr_diff = orr_diff - PI;
470                  if (orr_diff< -PI/2.0) orr_diff = orr_diff + PI;
471                 
472 /*
473                  hfill1(hists[0],orr_diff,1.0);
474 */
475 
476                  pos2 = edgel_realign(edge2, edge1);
477                  estimate = pos2 - ((double)j1 + pos1);
478                 
479                  /* add match property to allow uniqueness checking
480                     check for existing match from previous run first */
481                  if (prop_present(edge1->props, MATCH))
482                      edge1->props = proplist_rm(edge1->props, MATCH);
483                    
484                  if (fabs(orr_diff)< ANGERR)
485                  {
486                      IM_FLOAT(out_disp_im, i+k, j1+l) = estimate;
487                      match = match_make((void *)edge1, (void *)edge2, MATCH);
488                      edge1->props = proplist_addifnp(edge1->props, match, MATCH, 
489                                                 match_free, true);
490                  }
491              }
492           }
493        }
494     }
495 }
496 

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