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

Linux Cross Reference
Tina4/src/pgh/pgh_hist.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

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

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