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

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

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