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

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

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