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

Linux Cross Reference
Tina5/tina-libs/tina/vision/visPgh_hist.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/visPgh_hist.c,v $
 23  * Date    :  $Date: 2008/12/02 22:03:07 $
 24  * Version :  $Revision: 1.6 $
 25  * CVS Id  :  $Id: visPgh_hist.c,v 1.6 2008/12/02 22:03:07 paul Exp $
 26  *
 27  * Author  : NAT
 28  */
 29 /**
 30  *   @file Functions for the construction of variants of the pairwise geometric histogram.
 31  *   @brief Straight line segments are encoded in a 2D plot of perpendicular distance versus
 32  *          relative angle. Errors are explicitly coded as convolutions with a Gaussian
 33  *          function so that resulting plots are estimates of a geometric density distributiuon
 34  *          and can then be matched using a Bhattacharrya distance.
 35  *
 36  *
 37 */
 38 
 39 #include "visPgh_hist.h"
 40 #include <stdio.h>
 41 #include <math.h>
 42 #include <tina/vision/visPgh_defs.h>
 43 #include <tina/math/mathGeom_vec2.h>
 44 #include <tina/math/math_TypPro.h>
 45 #include <tina/image/imgGen_alloc.h>
 46 #include <tina/geometry/geomImg_poly_crop.h>
 47 
 48 
 49 static double dbin_min,dbin_max,dbin_size;
 50 static int num_abin,pairs_type=-1;
 51 static double abin_size,angle_sigma,dist_ramp;
 52 static float *pair_dist=NULL;
 53 static float *pair_angle=NULL;
 54 static double erf_array[81];
 55 
 56 double norm_prob(double sigma,double mu,double y)
 57 {
 58 /* returns probability that a random number generated with a normal
 59    distribution around mu with variance sigma would be less than y */
 60     double x;
 61     double prob;
 62     int index;
 63 
 64     x = (y-mu)/sigma;
 65     index  = tina_int(x/ERFACC)+40;
 66 
 67     if (index<0) prob = 0.0; 
 68     else if (index>80) 
 69                 prob = 1.0;
 70     else 
 71                 prob = 0.5*(1.0+erf_array[index]);
 72 
 73     return(prob);
 74 }
 75 
 76 void init_erf(void)
 77 {
 78     double y;
 79 
 80     for (y=-ERFRANGE+ERFACC/2.0;y<=ERFRANGE+ERFACC;y+=ERFACC)
 81     {
 82        erf_array[tina_int(y/ERFACC)+40] = (double) erf((float)y);
 83 /*         printf("ERF %lf y %lf %d\n",(double)erf(y),y,tina_int(y/ERFACC)+40);
 84 */
 85     }
 86 }
 87 
 88 void get_gdist(float *gdist,double angle,double bin_width,double sigma,
 89                int *low_angle_bin,int *high_angle_bin)
 90 {
 91    double angle_fac, integral = 0.0;
 92    int bin_edge;
 93 
 94    *low_angle_bin = tina_int((angle - ERFRANGE*sigma)/bin_width);
 95    *high_angle_bin = tina_int((angle + ERFRANGE*sigma)/bin_width);
 96 
 97    for (bin_edge = *low_angle_bin;bin_edge<=*high_angle_bin; bin_edge++)
 98    {
 99        angle_fac = norm_prob(angle_sigma,angle,(bin_edge+1)*abin_size)- integral;
100        integral += angle_fac;
101        gdist[bin_edge] = (float)angle_fac;
102    }
103 }
104 
105 void get_pdist(float *pdist,double pdmin,double pdmax,double bin_width,double ramp,
106                int *low_bin,int *high_bin)
107 /* routine to fill a perpedicular distance array with ramped weighted values */
108 /*                                                              NAT 27/10/92 */
109 {
110     double pbmin,pbmax;
111     double pddiff,ptmin,ptmax,mean;
112     double total=0.0;
113     int i,j;
114 
115     pddiff = pdmax - pdmin;
116     mean =  (pdmax+pdmin)/2.0;
117     pbmin = mean - 0.5*pddiff - 0.5*ramp;
118     pbmax = mean + 0.5*pddiff + 0.5*ramp;
119     for (i = tina_int(pbmin/bin_width) ; i <=tina_int(pbmax/bin_width); i++)
120        pdist[i]=(float)0.0;
121 
122 
123     /* set the blurring width of the data to 'ramp' or 'pddiff' whichever is smaller */
124 
125     if (pddiff<ramp) 
126     {
127         ptmin = pbmin + pddiff;
128         ptmax = pbmax - pddiff;
129     }
130     else
131     {
132         ptmin = pbmin + ramp;
133         ptmax = pbmax - ramp; 
134     } 
135     
136     /* deal with the linear section of the trapezium */     
137     i = tina_int(ptmin/bin_width);
138     j = tina_int(ptmax/bin_width);
139         
140     if (i==j)
141     {
142         pdist[i] = (float)((ptmax-ptmin)/bin_width);
143     }
144     else
145     {
146         pdist[i] = (float)(1.0+i-ptmin/bin_width);
147         for (i = tina_int(ptmin/bin_width) +1 ; i<tina_int(ptmax/bin_width); i++)
148             pdist[i] = (float)1.0;
149         pdist[j] = (float)(ptmax/bin_width-j);
150     }
151 
152     /*  now fill leading and trailing triangles */
153 
154     i = tina_int(pbmin/bin_width);
155     j = tina_int(ptmin/bin_width);
156 
157     if (i==j)
158     {
159         pdist[i] += (float)(0.5*(ptmin-pbmin)/bin_width);
160     }
161     else
162     {
163         pdist[i] += (float)(0.5*(1.0+i-pbmin/bin_width)*((i+1.0)*bin_width-pbmin)/(ptmin-pbmin));
164 
165         for (i = tina_int(pbmin/bin_width) +1 ; i<tina_int(ptmin/bin_width); i++)
166             pdist[i] += (float)(((i+0.5)*bin_width-pbmin)/(ptmin-pbmin));
167 
168         pdist[j] += (float)(((j*bin_width+ptmin)/2.0 - pbmin)*(ptmin/bin_width-j)/(ptmin-pbmin));
169     }
170 
171     i = tina_int(ptmax/bin_width);
172     j = tina_int(pbmax/bin_width);
173 
174     if (i==j)
175     {
176         pdist[i] += (float)(0.5*(pbmax-ptmax)/bin_width);
177     }
178     else
179     {
180         pdist[i] += (float)((pbmax - ((i+1)*bin_width+ptmax)/2.0)
181                    *((i+1)-ptmax/bin_width)/(pbmax-ptmax));
182 
183         for (i = tina_int(ptmax/bin_width) +1 ; i<tina_int(pbmax/bin_width); i++)
184             pdist[i] += (float)((pbmax-(i+0.5)*bin_width)/(pbmax-ptmax)); 
185 
186         pdist[j] += (float)(0.5*(pbmax-j*bin_width)*(pbmax/bin_width-j)/(pbmax-ptmax));
187     }
188 
189    /* now normalise */
190 
191    for (i = tina_int(pbmin/bin_width) ; i <=tina_int(pbmax/bin_width); i++)
192        total += pdist[i];
193 
194    for (i = tina_int(pbmin/bin_width) ; i <=tina_int(pbmax/bin_width); i++)
195        pdist[i] /= (float)total;
196 
197    *low_bin = tina_int(pbmin/bin_width);
198    *high_bin = tina_int(pbmax/bin_width);
199 }
200 
201 
202 Vec2 dir_vec(Vec2 isct,Vec2 p1,Vec2 p2)
203 /* routine to calculate the direction of the line p1_p2 away from isct */
204 /* replacement for Aluns modify dir_vec.                  NAT 11/11/92 */
205 {
206   Vec2 tv_1;
207 
208   if (vec2_dist(p1,isct)>vec2_dist(p2,isct))
209   {
210      tv_1 = vec2_unit(vec2_diff(p1,p2));
211   }
212   else
213   {
214      tv_1 = vec2_unit(vec2_diff(p2,p1));
215   }
216   return (tv_1);
217 }
218 
219 
220 void init_pairs_entry(int new_pairs_type,double newdbin_max,double newdbin_size,
221                       int newnum_abin,double newangle_sigma,double newdist_ramp)
222 /* set up the dimensions of the pairwise histogram NAT 11/11/92 */
223 {
224     int ndbin;
225     
226     pairs_type = new_pairs_type;
227     if (pair_dist !=NULL) 
228                 rfree(pair_dist+tina_int(dbin_min/dbin_size));
229 
230     if (pair_angle !=NULL)
231                 rfree(pair_angle-tina_int(1.0+3.0*angle_sigma/abin_size));
232 
233     dbin_max = newdbin_max;
234     dbin_min = -newdbin_max;
235 
236     if (pairs_type == MIRROR) 
237                 dbin_min = -newdist_ramp;
238 
239     dbin_size = newdbin_size;
240     num_abin = newnum_abin;
241     abin_size = 3.141592/num_abin;
242 
243     if (pairs_type == DIRECTED ) 
244                 num_abin*=2;
245 
246     angle_sigma = newangle_sigma;
247     dist_ramp = newdist_ramp;
248     ndbin = (int)((dbin_max - dbin_min)/dbin_size);
249     pair_dist = (float*) ralloc(sizeof(float)*ndbin);
250     pair_dist -= tina_int(dbin_min/dbin_size);
251     pair_angle = (float*) ralloc(sizeof(float)*(2*num_abin+
252                                     2*tina_int(1.0+3.0*angle_sigma/abin_size)));
253     
254         pair_angle += tina_int(1.0+3.0*angle_sigma/abin_size);
255 
256 
257     init_erf();
258 
259     set_pairs_rotate(dbin_min,dbin_max,dbin_size,num_abin,pair_dist,
260                      pair_angle,abin_size,angle_sigma,dist_ramp);
261     set_pairs_mirror(dbin_min,dbin_max,dbin_size,num_abin,pair_dist,
262                      pair_angle,abin_size,angle_sigma,dist_ramp);
263     set_pairs_directed(dbin_min,dbin_max,dbin_size,num_abin,pair_dist,
264                      pair_angle,abin_size,angle_sigma,dist_ramp);
265 }
266 
267 Imregion im_pairs_roi(void)
268 {
269     Imregion roi;
270 
271     roi.lx = 0;
272     roi.ux = num_abin;
273     roi.ly = tina_int(-dbin_max/dbin_size);
274     if (pairs_type == MIRROR) roi.ly = 0;
275     roi.uy = tina_int(dbin_max/dbin_size);
276 
277     /* Added to prevent compiler warning PAB 02/12/2008*/
278 
279     roi.ts_id = 0;
280 
281     return(roi);
282 }
283 
284 Imrect *im_pairs_alloc(void)
285 {
286     Imregion roi;
287     Imrect *im;
288     int width,height;
289 
290     roi = im_pairs_roi();
291     width = num_abin;
292     height = roi.uy-roi.ly;
293 
294     im = im_alloc(height, width, &roi, float_v);
295 
296     return(im);
297 }
298 
299 void compare_lines(Imrect *im,Line2 *l1,Line2 *l2,int *type)
300 /* add pairwise contributions for lines l1 and l2 to histogram in im */
301 /* l1 is reference line, l2 is comparison line.         NAT 29/7/93 */
302 {
303 
304     *type=pairs_type;
305 
306     switch(pairs_type)
307     {
308       case MIRROR:
309       compare_lines_mirror(im,l1,l2);
310     break;
311       case ROTATE:
312       compare_lines_rotate(im,l1,l2);
313     break;
314       case DIRECTED:
315       compare_lines_directed(im,l1,l2);
316     break;
317     }
318 }
319 
320 /* ---------- mirror ---------- */
321 
322 void set_pairs_mirror(double new_dbin_min,double new_dbin_max,
323                       double new_dbin_size,int new_num_abin,
324                       float *new_pair_dist,float *new_pair_angle,
325                       double new_abin_size,
326                       double new_angle_sigma,double new_dist_ramp)
327 {
328      dbin_min = new_dbin_min;
329      dbin_max = new_dbin_max;
330      dbin_size = new_dbin_size;
331      num_abin = new_num_abin;
332      pair_dist = new_pair_dist;
333      pair_angle = new_pair_angle;
334      abin_size = new_abin_size;
335      angle_sigma = new_angle_sigma;
336      dist_ramp = new_dist_ramp;
337 }
338 
339 void make_entry_mirror(Imrect *im,double min_dist,double max_dist,
340                        double angle,double weight)
341 {
342    int low_angle_bin,high_angle_bin;
343    int low_dist_bin,high_dist_bin;
344    int bin_edge,dbin,gbin;
345    int angle_bin;
346    double current_contents;
347    double angle_fac;
348    double temp;
349    Imregion *roi;
350  
351    if (im==NULL) 
352            return;
353 
354    if ((roi=im->region) == NULL) 
355            return;
356 
357    if (min_dist>max_dist)
358    {
359        temp = min_dist;
360        min_dist = max_dist;
361        max_dist = temp;
362    }
363 
364    if (min_dist>dbin_max-dist_ramp) return;
365    if (max_dist<dbin_min+dist_ramp) return;
366    if (min_dist<dbin_min+dist_ramp)
367    {
368       weight *= (dbin_min+dist_ramp-min_dist)/(max_dist-min_dist);
369       min_dist = dbin_min+dist_ramp;
370    }
371    if (max_dist>dbin_max-dist_ramp)
372    {
373       weight *= (dbin_max-dist_ramp-min_dist)/(max_dist-min_dist);
374       max_dist = dbin_max-dist_ramp;
375    }
376  
377    roi = im->region;
378 
379    get_pdist(pair_dist,min_dist,max_dist,dbin_size,dist_ramp,
380                                   &low_dist_bin,&high_dist_bin);
381    get_pdist(pair_angle,angle-angle_sigma,angle+angle_sigma,abin_size,
382                    2.0*angle_sigma,&low_angle_bin,&high_angle_bin);
383  
384    for (bin_edge = low_angle_bin;bin_edge<=high_angle_bin; bin_edge++ )
385    {
386       angle_fac = pair_angle[bin_edge];
387       for (dbin = low_dist_bin;dbin<=high_dist_bin;dbin++)
388       {
389           angle_bin = bin_edge;
390           while (angle_bin<0) angle_bin = angle_bin+num_abin;
391           while (angle_bin>=num_abin) angle_bin = angle_bin-num_abin;
392           gbin = dbin;
393  
394           if (dbin<0)
395           {
396             gbin = - dbin -1;
397             angle_bin = num_abin - angle_bin -1;
398             if (angle_bin>=num_abin) angle_bin -= num_abin;
399           }
400           IM_PIX_GET(im, gbin, angle_bin, current_contents);
401           IM_PIX_SET(im, gbin, angle_bin, current_contents +
402                                           angle_fac*weight*pair_dist[dbin]);
403           if (current_contents + angle_fac*weight*pair_dist[dbin]<0.0)
404                printf("bollocks\n");
405       }
406    }
407 }
408  
409 void compare_lines_mirror(Imrect *im, Line2 *l1, Line2 *l2)
410 /* add pairwise contributions for lines l1 and l2 to histogram in im */
411 /* l1 is reference line, l2 is comparison line.         NAT 10/11/92 */
412 {
413    Vec2 v1,v2,normv1,isct;
414    Bool l1_intsct,l2_intsct;
415    double weight,min_dist,max_dist,angle;
416 
417    isct = get_intersection(l1,l2,&l1_intsct,&l2_intsct);
418  
419    if( !l1_intsct && !l2_intsct)
420    {
421        /* both lines no intersections */
422        v1 = dir_vec(isct,l1->p1,l1->p2);
423        v2 = dir_vec(isct,l2->p1,l2->p2);
424        angle = fabs(vec2_angle(v1,v2));
425        normv1 = vec2(-((l1->v).el[1]),(l1->v).el[0]);
426        min_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p1,isct)));
427        max_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p2,isct)));
428        weight = vec2_dist(l1->p1,l1->p2)*vec2_dist(l2->p1,l2->p2);
429        make_entry_mirror(im,min_dist,max_dist,angle,weight);
430    }
431    else
432    {
433        if(l2_intsct)
434        {
435           if(l1_intsct)
436           {
437              /* l1_p1 to isct and l2_p1 to isct */
438              v1 = dir_vec(isct,l1->p1,isct);
439              v2 = dir_vec(isct,l2->p1,isct);
440              angle = fabs(vec2_angle(v1,v2));
441              normv1 = vec2(-((l1->v).el[1]),(l1->v).el[0]);
442              min_dist = 0.0;
443              max_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p1,isct)));
444              weight = vec2_dist(l1->p1,isct)*vec2_dist(l2->p1,isct);
445              make_entry_mirror(im,min_dist,max_dist,angle,weight);
446  
447              /* isct to l1_p2 and l2_p1 to isct */
448              weight = vec2_dist(l1->p2,isct)*vec2_dist(l2->p1,isct);
449              make_entry_mirror(im,min_dist,max_dist,PI-angle,weight);
450  
451              /* l1-p1 to isct and isct to l2_p2 */
452              max_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p2,isct)));
453              weight = vec2_dist(l1->p1,isct)*vec2_dist(isct,l2->p2);
454              make_entry_mirror(im,min_dist,max_dist,PI-angle,weight);
455  
456              /* isct to l1_p2 and isct to l2_p2 */
457              weight = vec2_dist(l1->p2,isct)*vec2_dist(isct,l2->p2);
458              make_entry_mirror(im,min_dist,max_dist,angle,weight);
459           }
460           else
461           {
462              /* l1_p1 to l1_p2 and l2_p1 to isct */
463              v1 = dir_vec(isct,l1->p1,l1->p2);
464              v2 = dir_vec(isct,l2->p1,isct);
465              angle = fabs(vec2_angle(v1,v2));
466              normv1 = vec2(-((l1->v).el[1]),(l1->v).el[0]);
467              min_dist = 0.0;
468              max_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p1,isct)));
469              weight = vec2_dist(l1->p1,l1->p2)*vec2_dist(l2->p1,isct);
470              make_entry_mirror(im,min_dist,max_dist,angle,weight);
471  
472              /* l1_p1 to l1_p2 and isct to l2_p2 */
473              max_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p2,isct)));
474              weight = vec2_dist(l1->p1,l1->p2)*vec2_dist(isct,l2->p2);
475              make_entry_mirror(im,min_dist,max_dist,PI-angle,weight);
476           }
477        }
478        else
479        {
480            /* l1_p1 to isct and l2_p1 to l2_p2 */
481            v1 = dir_vec(isct,l1->p1,isct);
482            v2 = dir_vec(isct,l2->p1,l2->p2);
483            angle = fabs(vec2_angle(v1,v2));
484            normv1 = vec2(-((l1->v).el[1]),(l1->v).el[0]);
485            min_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p1,isct)));
486            max_dist = fabs(vec2_dot(normv1,vec2_diff(l2->p2,isct)));
487            weight = vec2_dist(l1->p1,isct)*vec2_dist(l2->p1,l2->p2);
488            make_entry_mirror(im,min_dist,max_dist,angle,weight);
489  
490            /* isct to l1_p2 and l2_p1 to l2_p2 */
491            weight = vec2_dist(isct,l1->p2)*vec2_dist(l2->p1,l2->p2);
492            make_entry_mirror(im,min_dist,max_dist,PI-angle,weight);
493        }
494    }
495 }
496 
497 /* ---------- rotate ---------- */
498 
499 void set_pairs_rotate(double new_dbin_min, double new_dbin_max,
500                       double new_dbin_size, int new_num_abin,
501                       float *new_pair_dist, float *new_pair_angle,
502                       double new_abin_size,
503                       double new_angle_sigma, double new_dist_ramp)
504 {
505      dbin_min = new_dbin_min;
506      dbin_max = new_dbin_max;
507      dbin_size = new_dbin_size;
508      num_abin = new_num_abin;
509      pair_dist = new_pair_dist;
510      pair_angle = new_pair_angle;
511      abin_size = new_abin_size;
512      angle_sigma = new_angle_sigma;
513      dist_ramp = new_dist_ramp;
514 }
515 
516 void make_entry_rotate(Imrect *im, double min_dist, double max_dist,
517                        double angle, double weight)
518 {
519    int low_angle_bin,high_angle_bin;
520    int low_dist_bin,high_dist_bin;
521    int bin_edge,dbin,gbin;
522    int angle_bin;
523    double current_contents;
524    double angle_fac;
525    double temp;
526    Imregion *roi;
527    int err=0;
528  
529 
530    if (im==NULL) 
531            return;
532    if ((roi=im->region) == NULL) 
533            return;
534    
535    if (min_dist>max_dist)
536    {
537        temp = min_dist;
538        min_dist = max_dist;
539        max_dist = temp;
540    }
541 
542    if (min_dist>dbin_max-dist_ramp) 
543            return;
544 
545    if (max_dist<dbin_min+dist_ramp) 
546            return;
547 
548    if (min_dist<dbin_min+dist_ramp)
549    {
550       weight *= (max_dist-dbin_min-dist_ramp)/(max_dist-min_dist);
551       min_dist = dbin_min+dist_ramp;
552    }
553 
554    if (max_dist>dbin_max-dist_ramp)
555    {
556       weight *= (dbin_max-dist_ramp-min_dist)/(max_dist-min_dist);
557       max_dist = dbin_max-dist_ramp;
558    }
559 
560    roi = im->region;
561  
562    get_pdist(pair_dist,min_dist,max_dist,dbin_size,dist_ramp,
563                                                  &low_dist_bin,&high_dist_bin);
564    
565    get_pdist(pair_angle,angle-angle_sigma,angle+angle_sigma,abin_size,
566                    2.0*angle_sigma,&low_angle_bin,&high_angle_bin);
567 
568                                    /* gaussian blurring version NAT 11/11/93
569    get_gdist(pair_angle,angle,abin_size,angle_sigma,
570                       &low_angle_bin,&high_angle_bin);
571 */
572 
573    for (bin_edge = low_angle_bin;bin_edge<=high_angle_bin; bin_edge++)
574    {
575       angle_fac = pair_angle[bin_edge];
576       for (dbin = low_dist_bin;dbin<=high_dist_bin;dbin++)
577       {
578           if (bin_edge<0)
579           {
580                angle_bin = bin_edge+num_abin;
581                gbin = -dbin -1;
582           }
583           else if (bin_edge>=num_abin)
584           {
585                angle_bin = bin_edge-num_abin;
586                gbin = -dbin -1;
587           }
588           else /* generally the case */
589           {
590                angle_bin = bin_edge;
591                gbin = dbin;
592           }
593  
594 
595         
596           IM_PIX_GET(im, gbin, angle_bin, current_contents);
597                   
598                   /*printf("gbin %d angle_bin %d angle_fac %lf weight %lf pair_dist[dbin] %lf current_contents %lf\n", gbin, angle_bin,
599 angle_fac,weight, pair_dist[dbin],current_contents);*/
600 
601           IM_PIX_SET(im, gbin, angle_bin, current_contents +
602                                           angle_fac*weight*pair_dist[dbin]);
603 
604 
605                   /*check the values written*/
606 
607           if ((current_contents + angle_fac*weight*pair_dist[dbin])<0.0)
608                   {
609                           err=1;
610                           printf("gbin %d angle_bin %d angle_fac %g weight %g pair_dist[dbin] %g current_contents %g\n", 
611                                   gbin, angle_bin, angle_fac,weight, pair_dist[dbin],current_contents);
612                   }
613 
614       }
615    }
616 
617    if(err==1)
618          printf("Error! negative entries in pairwise\n");
619 
620 }
621 
622 void compare_lines_rotate(Imrect *im,Line2 *l1, Line2 *l2)
623 /* add pairwise contributions for lines l1 and l2 to histogram in im */
624 /* l1 is reference line, l2 is comparison line.         NAT 10/11/92 */
625 {
626    Vec2 v1,normv1,isct;
627    Bool l1_intsct,l2_intsct;
628    double weight,min_dist,max_dist,angle;
629  
630    isct = get_intersection(l1,l2,&l1_intsct,&l2_intsct);
631  
632    if( !l1_intsct)
633    {
634        /* both lines no intersections */
635        v1 = dir_vec(isct,l1->p1,l1->p2);
636        angle = vec2_angle(l2->v,v1);
637 
638        if (angle<0.0) 
639                    angle += PI;
640 
641        normv1 = vec2(-v1.el[1],v1.el[0]);
642        min_dist = vec2_dot(normv1,vec2_diff(l2->p1,isct));
643        max_dist = vec2_dot(normv1,vec2_diff(l2->p2,isct));
644        weight = vec2_dist(l1->p1,l1->p2)*vec2_dist(l2->p1,l2->p2);
645        make_entry_rotate(im,min_dist,max_dist,angle,weight);
646    }
647    else
648    {
649        /* l1_p1 to isct and l2_p1 to l2_p2 */
650        v1 = dir_vec(isct,l1->p1,isct);
651        angle = vec2_angle(l2->v,v1);
652        if (angle<0.0) angle += PI;
653        normv1 = vec2(-v1.el[1],v1.el[0]);
654        min_dist = vec2_dot(normv1,vec2_diff(l2->p1,isct));
655        max_dist = vec2_dot(normv1,vec2_diff(l2->p2,isct));
656        weight = vec2_dist(l1->p1,isct)*vec2_dist(l2->p1,l2->p2);
657        make_entry_rotate(im,min_dist,max_dist,angle,weight);
658  
659        /* isct to l1_p2 and l2_p1 to l2_p2 */
660        weight = vec2_dist(isct,l1->p2)*vec2_dist(l2->p1,l2->p2);
661        make_entry_rotate(im,-min_dist,-max_dist,angle,weight);
662    }
663 }
664 
665 /* ---------- directed ---------- */
666 
667 /* Routines for generating line directed pairwise histograms in the range */
668 /* 0 -> 2 PI and +- distance.                        NAT 29/7/93          */
669 
670 void set_pairs_directed(double new_dbin_min, double new_dbin_max,
671                         double new_dbin_size, int new_num_abin,
672                         float *new_pair_dist, float *new_pair_angle,
673                         double new_abin_size,
674                         double new_angle_sigma, double new_dist_ramp)
675 {
676      dbin_min = new_dbin_min;
677      dbin_max = new_dbin_max;
678      dbin_size = new_dbin_size;
679      num_abin = new_num_abin;
680      pair_dist = new_pair_dist;
681      pair_angle = new_pair_angle;
682      abin_size = new_abin_size;
683      angle_sigma = new_angle_sigma;
684      dist_ramp = new_dist_ramp;
685 }
686 
687 void make_entry_directed(Imrect *im, double min_dist, double max_dist,
688                          double angle, double weight)
689 /* treat the poles in the representation at 0 PI and 2PI */
690 {
691    double adiff;
692 
693    if ((adiff = fabs(angle))<angle_sigma)
694    {
695       make_entry_direct(im,min_dist,max_dist,angle,
696                            weight*(0.5+0.5*adiff/angle_sigma));
697       make_entry_direct(im,min_dist,max_dist,angle+PI,
698                            weight*(0.5-0.5*adiff/angle_sigma));
699    }
700    else if ((adiff = fabs(angle-PI))<angle_sigma)
701    {
702       make_entry_direct(im,min_dist,max_dist,angle,
703                            weight*(0.5+0.5*adiff/angle_sigma));
704       if ((angle-PI)>=0.0)
705           make_entry_direct(im,min_dist,max_dist,angle-PI,
706                                weight*(0.5-0.5*adiff/angle_sigma));
707       else
708           make_entry_direct(im,min_dist,max_dist,angle+PI,
709                                weight*(0.5-0.5*adiff/angle_sigma));
710    }
711    else if ((adiff=fabs(angle-2.0*PI))<angle_sigma)
712    {
713       make_entry_direct(im,min_dist,max_dist,angle,
714                            weight*(0.5+0.5*adiff/angle_sigma));
715       make_entry_direct(im,min_dist,max_dist,angle-PI,
716                            weight*(0.5-0.5*adiff/angle_sigma));
717 
718    }
719    else  make_entry_direct(im,min_dist,max_dist,angle,weight);
720 }
721 
722 void make_entry_direct(Imrect *im, double min_dist, double max_dist,
723                        double angle, double weight)
724 {
725    int low_angle_bin,high_angle_bin;
726    int low_dist_bin,high_dist_bin;
727    int bin_edge,dbin,gbin;
728    int angle_bin;
729    double current_contents;
730    double angle_fac;
731    double temp;
732    Imregion *roi;
733  
734    if (im==NULL) return;
735    if ((roi=im->region) == NULL) return;
736    if (min_dist>max_dist)
737    {
738        temp = min_dist;
739        min_dist = max_dist;
740        max_dist = temp;
741    }
742 
743    if (min_dist>dbin_max-dist_ramp) return;
744    if (max_dist<dbin_min+dist_ramp) return;
745    if (min_dist<dbin_min+dist_ramp)
746    {
747       weight *= (max_dist-dbin_min-dist_ramp)/(max_dist-min_dist);
748       min_dist = dbin_min+dist_ramp;
749    }
750    if (max_dist>dbin_max-dist_ramp)
751    {
752       weight *= (dbin_max-dist_ramp-min_dist)/(max_dist-min_dist);
753       max_dist = dbin_max-dist_ramp;
754    }
755  
756    roi = im->region;
757 
758    get_pdist(pair_dist,min_dist,max_dist,dbin_size,dist_ramp,
759                                          &low_dist_bin,&high_dist_bin);
760    get_pdist(pair_angle,angle-angle_sigma,angle+angle_sigma,abin_size,
761                    2.0*angle_sigma,&low_angle_bin,&high_angle_bin);
762  
763    for (bin_edge = low_angle_bin;bin_edge<=high_angle_bin; bin_edge++ )
764    {
765       angle_fac = pair_angle[bin_edge];
766       for (dbin = low_dist_bin;dbin<=high_dist_bin;dbin++)
767       {
768           angle_bin = bin_edge;
769           while (angle_bin<0) angle_bin = angle_bin+num_abin;
770           while (angle_bin>=num_abin) angle_bin = angle_bin-num_abin;
771           gbin = dbin;
772  
773           IM_PIX_GET(im, gbin, angle_bin, current_contents);
774           IM_PIX_SET(im, gbin, angle_bin, current_contents +
775                                           angle_fac*weight*pair_dist[dbin]);
776           if (current_contents + angle_fac*weight*pair_dist[dbin]<0.0)
777                format("bollocks again");
778       }
779    }
780  
781 }
782  
783 void compare_lines_directed(Imrect *im,Line2 *l1,Line2 *l2)
784 /* add pairwise contributions for lines l1 and l2 to histogram in im */
785 /* l1 is reference line, l2 is comparison line.         NAT 10/11/92 */
786 {
787    Vec2 v1,normv1,isct;
788    Bool l1_intsct,l2_intsct;
789    double weight,min_dist,max_dist,angle;
790  
791    isct = get_intersection(l1,l2,&l1_intsct,&l2_intsct);
792  
793    if(!l1_intsct)
794    {
795 /* both lines no intersections */
796        v1 = dir_vec(isct,l1->p1,l1->p2);
797        angle = vec2_angle(l2->v,v1);
798        if (angle<0.0) angle += PI;
799 /* if intersection direction disagrees with nominal line direction */
800 /* then flip into Pi-2Pi domain of plot                            */
801        if (vec2_dot(l1->v,v1)<0.0) angle += PI;
802        normv1 = vec2(-((l1->v).el[1]),(l1->v).el[0]);
803        min_dist = vec2_dot(normv1,vec2_diff(l2->p1,isct));
804        max_dist = vec2_dot(normv1,vec2_diff(l2->p2,isct));
805        weight = vec2_dist(l1->p1,l1->p2)*vec2_dist(l2->p1,l2->p2);
806        make_entry_directed(im,min_dist,max_dist,angle,weight);
807    }
808    else
809    {
810 /* l1_p1 to isct and l2_p1 to l2_p2 */
811        v1 = dir_vec(isct,l1->p1,isct);
812        angle = vec2_angle(l2->v,v1);
813        if (angle<0.0) angle += PI;
814 /* if intersection direction disagrees with nominal line direction */
815 /* then flip into Pi-2Pi domain of plot                            */
816        if (vec2_dot(l1->v,v1)<0.0) angle += PI;
817        normv1 = vec2(-((l1->v).el[1]),(l1->v).el[0]);
818        min_dist = vec2_dot(normv1,vec2_diff(l2->p1,isct));
819        max_dist = vec2_dot(normv1,vec2_diff(l2->p2,isct));
820        weight = vec2_dist(l1->p1,isct)*vec2_dist(l2->p1,l2->p2);
821        make_entry_directed(im,min_dist,max_dist,angle,weight);
822  
823 /* isct to l1_p2 and l2_p1 to l2_p2 */
824        weight = vec2_dist(isct,l1->p2)*vec2_dist(l2->p1,l2->p2);
825        make_entry_directed(im,min_dist,max_dist,PI+angle,weight);
826    }
827 }
828 

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