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
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.