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

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

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

  1 #include <stdio.h>
  2 #include <sys/param.h>
  3 #include <string.h>
  4 #include <tina/sys.h>
  5 #include <tina/sysfuncs.h>
  6 #include <tina/math.h>
  7 #include <tina/mathfuncs.h>
  8 #include <tina/vision.h>
  9 #include <tina/visionfuncs.h>
 10 #include <tina/file.h>
 11 #include <tina/filefuncs.h>
 12 #include <tina/file_gen.h>
 13 #include <tina/draw.h>
 14 #include <tina/drawfuncs.h>
 15 #include <tina/tv.h>
 16 #include <tina/tvfuncs.h>
 17 #include <tina/toolsfuncs.h>
 18 #include <tina/tw_Xfuncs.h>
 19 #include <pgh/pgh_defs.h>
 20 #include <pgh/pgh_types.h>
 21 #include <pgh/pgh_funcs.h>
 22 
 23 /* ---------- Functions ---------- */
 24 
 25 static void sqr_root_pairwise(Imrect *hist)
 26 {
 27     int     i,j, min_x, min_y, max_x, max_y;
 28     double  hist_data;
 29 
 30     min_x = hist->region->lx;
 31     min_y = hist->region->ly;
 32     max_x = hist->region->ux;
 33     max_y = hist->region->uy;
 34 
 35 /* Each histogram entry is replaced by its square root */
 36 
 37     for(i=min_y;i<max_y;i++)
 38         for(j=min_x;j<max_x;j++)
 39         {
 40             IM_PIX_GET(hist,i,j,hist_data);
 41             if (hist_data < 0.0) hist_data = 0.0;
 42             IM_PIX_SET(hist,i,j,(double)sqrt((double)hist_data));
 43         }
 44 }
 45 
 46 /* ---------- */
 47 
 48 void sqr_root_and_normalize_pairwise(Imrect *hist)
 49 {
 50 int      i,j, min_x, min_y, max_x, max_y;
 51 double   sum=0.0;
 52 double   gl;
 53 
 54         /* (1) Find the sum of the histogram enties.
 55            (2) Divide each histogram entry by this value. 
 56            (3) Replace each histogram entry by its sqr root.
 57                This is done for statistical reasons. */
 58 
 59     min_x = hist->region->lx;
 60     min_y = hist->region->ly;
 61     max_x = hist->region->ux;
 62     max_y = hist->region->uy;
 63 
 64     for(i=min_y;i<max_y;i++)
 65         for(j=min_x;j<max_x;j++)
 66         {
 67             IM_PIX_GET(hist,i,j,gl);
 68             sum += gl;
 69         }
 70 
 71     if (sum!=0.0)
 72     {
 73         for(i=min_y;i<max_y;i++)
 74             for(j=min_x;j<max_x;j++)
 75             {
 76                 IM_PIX_GET(hist,i,j,gl)
 77                 IM_PIX_SET(hist,i,j,gl/(double)sum);
 78             }
 79 
 80         sqr_root_pairwise(hist);
 81     }
 82 }
 83 
 84 /* ---------- Functions ---------- */
 85 
 86 double sum_pairwise(Imrect *hist)
 87 {
 88     int     i,j, min_x, min_y, max_x, max_y;
 89     double  gl;
 90     double  sum = 0.0;
 91 
 92     min_x = hist->region->lx;
 93     min_y = hist->region->ly;
 94     max_x = hist->region->ux;
 95     max_y = hist->region->uy;
 96 
 97     for(i=min_y;i<max_y;i++)
 98         for(j=min_x;j<max_x;j++)
 99         {
100             IM_PIX_GET(hist,i,j,gl);
101             sum += gl;
102         }
103 
104   return (sum);
105 }
106 
107 /* ---------- */
108 
109 Imrect *build_pairwise(Line2 *ref_line, Model_poly_header *ref_model, List *line_list,
110                        double window_r)
111 {
112     Imrect   *pairwise, *im_pairs_alloc();
113     List     *ptr;
114     Line2    *lptr;
115     Hist_ref *ref;
116     double    l, l0;
117     double    sum_pairwise();
118     int type = -1;
119 
120         /* (1) Allocate space for the histogram.
121            (2) Compile the histogram from the ref line and
122                the line list.
123            (3) Sqrt the hist to get the correct statistics.
124            (4) Add a Hist_ref to the histograms props list.
125         */
126 
127     pairwise = im_pairs_alloc();
128     if (pairwise==NULL)
129     {
130         format("Error: Not space for pairwise histogram.\n");
131         return NULL;
132     }
133 
134     for(ptr=line_list;ptr!=NULL;ptr=ptr->next)
135     {
136         lptr=ptr->to;
137         if (vec2_dist(ref_line->p, lptr->p) < window_r)
138             compare_lines(pairwise, ref_line, lptr, &type);
139     }
140 
141     l0 = vec2_dist(ref_line->p1, ref_line->p2);
142     l  = sum_pairwise(pairwise)/l0;
143 
144     sqr_root_pairwise(pairwise);
145 
146     ref = ralloc(sizeof(Hist_ref));
147     if (ref==NULL)
148     {
149         format("Error: No space for Hist_ref.\n");
150         im_free(pairwise);
151         return NULL;
152     }
153 
154     ref->model = ref_model;
155     ref->matches = NULL;
156     ref->line_seg = ref_line; 
157     ref->l0 = l0;
158     ref->l = l;
159     ref->type = type;
160         ref->scale_factor = 1.0;
161 
162     pairwise->props = proplist_add(pairwise->props, ref, HIST_REF_TYPE, 
163                                                                hist_ref_free);
164 /*  remove fast indexing scheme for now NAT 15/5/95
165     findex = nonzero_index(pairwise);
166     pairwise->props = proplist_add(pairwise->props, findex, FAST_INDEX1, 
167                                                                vector_free);
168     if (type==DIRECTED)
169     {
170         findex = nonzero_index2(pairwise);
171         pairwise->props = proplist_add(pairwise->props, findex, FAST_INDEX2, 
172                                                                vector_free);
173     }
174     pairwise->label = ref_line->label;
175 */
176     return (pairwise);
177 
178 }
179 
180 /* ---------- */
181 
182 Imrect *build_normalized_pairwise(Line2 *ref_line, Model_poly_header *ref_model,
183                                   List *line_list, double window_r)
184 {
185     Imrect   *pairwise, *im_pairs_alloc();
186     List     *ptr;
187     Line2    *lptr;
188     Hist_ref *ref;
189     double    l, l0;
190     double    sum_pairwise();
191     int type = -1;
192 
193         /* (1) Allocate space for the histogram.
194            (2) Compile the histogram from the ref line and
195                the line list.
196            (3) Normalize the histogram.
197            (4) Sqrt the hist to get the correct statistics.
198            (5) Add a Hist_ref to the histograms props list.
199         */
200 
201     pairwise = im_pairs_alloc();
202     if (pairwise==NULL)
203     {
204         format("Error: Not space for pairwise histogram.\n");
205         return NULL;
206     }
207 
208     for(ptr=line_list;ptr!=NULL;ptr=ptr->next)
209     {
210         lptr=ptr->to;
211         if (vec2_dist(ref_line->p, lptr->p) < window_r)
212             compare_lines(pairwise, ref_line, lptr, &type);
213     }
214 
215     l0 = vec2_dist(ref_line->p1, ref_line->p2);
216     l = sum_pairwise(pairwise)/l0;
217 
218     sqr_root_and_normalize_pairwise(pairwise);
219 
220     ref = ralloc(sizeof(Hist_ref));
221     if (ref==NULL)
222     {
223         format("Error: No space for Hist_ref.\n");
224         im_free(pairwise);
225         return NULL;
226     }
227 
228     ref->model = ref_model;
229     ref->matches = NULL;
230     ref->line_seg = ref_line;
231     ref->l0 = l0;
232     ref->l = l;
233     ref->type = type;
234         ref->scale_factor = 1.0;
235 
236     pairwise->props = proplist_add(pairwise->props, ref, HIST_REF_TYPE, 
237                                                                hist_ref_free);
238 /*
239     pairwise->label = ref_line->label;
240 */
241     return (pairwise);
242 
243 }
244 
245 /* ---------- */
246 
247 void hist_ref_free(Hist_ref *hist_ref, int type)
248 {
249 
250         /* (1) Free the matches list. (Only the links,
251                not the data which they point to).
252            (2) Free the copied line seg.
253            (3) Free the hist_ref itself. */
254 
255     if (type!=HIST_REF_TYPE)
256     {
257         format("Error: Ilegal use of 'hist_ref_free'.\n");
258         return;
259     }
260 
261     list_rm_links(hist_ref->matches);
262 /*    line2_free(hist_ref->line_seg);
263       not to be done here NAT
264 */
265     rfree(hist_ref);
266 
267 }
268 
269 
270 /* ---------- */
271 
272 Imrect *pairs_build_norm_hist_scale(Line2 *ref_line, Model_poly_header *ref_model, List *line_list, double window_r, double scale)
273 {
274     Imrect   *pairwise, *im_pairs_alloc();
275     List     *ptr;
276     Line2    *lptr, *scaled_ref, *scaled_line, *clip_line;
277     Hist_ref *ref;
278     double    l, l0;
279     double    sum_pairwise();
280     int type = -1;
281         Transform2 trans;
282         Line2 *line2_circ_inter();
283         
284 
285 
286         /* (1) Allocate space for the histogram.
287            (2) Compile the histogram from the scaled ref line and
288                the scaled line list.
289            (3) Normalize the histogram.
290            (4) Sqrt the hist to get the correct statistics.
291            (5) Add a Hist_ref to the histograms props list.
292         */
293 
294         trans.R = rot2_with_scale(0.0, scale);
295         trans.t = vec2_zero();
296         
297         scaled_ref = line2_copy(ref_line);
298         line2_transform(scaled_ref, trans);
299 
300     pairwise = im_pairs_alloc();
301     if (pairwise==NULL)
302     {
303         format("Error: Not space for pairwise histogram.\n");
304         return NULL;
305     }
306 
307     for(ptr=line_list;ptr!=NULL;ptr=ptr->next)
308     {
309         lptr=ptr->to;
310                 scaled_line = line2_copy(lptr);
311                 line2_transform(scaled_line, trans);
312 
313                 clip_line = line2_circ_inter(scaled_line, scaled_ref->p, window_r);
314 
315 
316                 if (clip_line!=NULL)
317                 {
318             compare_lines(pairwise, scaled_ref, clip_line, &type);
319                 }
320 
321                 line2_free(scaled_line);
322                 line2_free(clip_line);
323     }
324 
325 
326     l0 = vec2_dist(scaled_ref->p1, scaled_ref->p2);
327     l = sum_pairwise(pairwise)/l0;
328 
329     sqr_root_and_normalize_pairwise(pairwise);
330 
331     ref = ralloc(sizeof(Hist_ref));
332     if (ref==NULL)
333     {
334         format("Error: No space for Hist_ref.\n");
335         im_free(pairwise);
336         return NULL;
337     }
338 
339         line2_free(scaled_ref);
340 
341     ref->model = ref_model;
342     ref->matches = NULL;
343     ref->line_seg = ref_line;
344     ref->l0 = l0;
345     ref->l = l;
346     ref->type = type;
347         ref->scale_factor = scale;
348 
349     pairwise->props = proplist_add(pairwise->props, ref, HIST_REF_TYPE, 
350                                                                hist_ref_free);
351    /*  No longer supported APA 19/6/95
352        pairwise->label = ref_line->label; */
353 
354     return (pairwise);
355 
356 }
357 
358 /* ---------- */
359 
360 Hist_ref *hist_ref_get(Imrect *hist)
361 {
362     Hist_ref *hist_ref;
363 
364     hist_ref = prop_get(hist->props, HIST_REF_TYPE);
365 
366     if (hist_ref==NULL)
367     {
368         printf("Error: Histogram has no hist_ref.\n");
369         return NULL;
370     }
371 
372     return hist_ref;
373 }
374 
375 /* ---------- */
376 
377 Line2 *ref_line_from_hist(Imrect *hist)
378 {
379     Hist_ref  *hist_ref;
380     Line2     *ref_line;
381 
382     if ((hist_ref = hist_ref_get(hist))==NULL) return NULL;
383 
384     if ((ref_line = hist_ref->line_seg)==NULL)
385     {
386         printf("Error: Histogram has no reference line.\n");
387         return NULL;
388     }
389 
390     return ref_line;
391 }    
392 
393 /* ---------- */
394 
395 List *model_geom_from_hist(Imrect *hist)
396 {
397     Hist_ref          *hist_ref;
398     Model_poly_header *model_poly_head;
399     List              *model_geom;
400 
401     if ((hist_ref = hist_ref_get(hist))==NULL) return NULL;
402 
403     if ((model_poly_head = hist_ref->model)==NULL)
404     {
405         printf("Error: Histogram has no model_poly_header.\n");
406         return NULL;
407     }
408 
409     if ((model_geom = model_poly_head->model_poly_data)==NULL)
410     {
411         printf("Error: Histogram has no model geometry.\n");
412         return NULL;
413     }
414 
415     return model_geom;
416 }
417 
418 /* ---------- */
419 
420 List *matches_list_from_hist(Imrect *hist)
421 {
422     Hist_ref          *hist_ref;
423     List              *matches_list;
424 
425     if ((hist_ref = hist_ref_get(hist))==NULL) return NULL;
426 
427     if ((matches_list = hist_ref->matches)==NULL)
428     {
429         printf("Error: Histogram has no matches.\n");
430         return NULL;
431     }
432 
433     return matches_list;
434 }
435 
436 /* ---------- */
437 
438 Hist_ref *hist_ref_copy(Hist_ref *hist_ref)
439 {
440     Hist_ref *copy;
441 
442     copy = ralloc(sizeof(Hist_ref));
443 
444     copy->line_seg = hist_ref->line_seg;
445     copy->model = hist_ref->model;
446     copy->l0 = hist_ref->l0;
447     copy->l = hist_ref->l;
448     copy->scale_factor = hist_ref->scale_factor;
449     copy->type = hist_ref->type;
450 
451     return copy;
452 }
453 
454 /* ---------- */
455 
456 /* Modifications to speed up histogram matching */
457 /* Installed NAT 18/3/94                        */
458 
459 /*
460 Vector *nonzero_index(Imrect *im)
461 {
462     Imregion *roi;
463     float *row,*row0;
464     int *index;
465     Vector *new_index;
466     Vector *vector_alloc();
467     int lx, ux, ly, uy;
468     int width, height;
469     int i, j, k;
470 
471     if (im == NULL)
472         return (NULL);
473 
474     roi = im->region;
475     if (roi == NULL)
476         return (NULL);
477     lx = roi->lx;
478     ux = roi->ux;
479     ly = roi->ly;
480     uy = roi->uy;
481 
482     index = (int *)ivector_alloc(0, (ux-lx)*(uy-ly));
483 
484     row0 = &(im->array.float_v[ly][lx]);
485 
486     for (i = ly, k=0; i < uy; ++i)
487     {
488         row = im->array.float_v[i];
489         for (j = lx; j < ux; ++j)
490         {
491             if (row[j]>0.0)
492             {
493                  index[k++] = &row[j] - row0;
494             }
495         }
496     }
497     new_index = vector_alloc(k, int_v);
498 
499     for (i=0;i<k;i++)
500         new_index->el.int_v[i] = index[i];
501 
502     ivector_free(index, 0);
503 
504     return (new_index);
505 }
506 */
507 
508 /*
509 Vector *nonzero_index2(Imrect *im)
510 {
511     Imregion *roi;
512     float *row,*row0;
513     int *index;
514     Vector *new_index;
515     Vector *vector_alloc();
516     int lx, ux, ly, uy;
517     int width, height;
518     int i, j, k, ishift;
519 
520     if (im == NULL)
521         return (NULL);
522 
523     roi = im->region;
524     if (roi == NULL)
525         return (NULL);
526     lx = roi->lx;
527     ux = roi->ux;
528     ly = roi->ly;
529     uy = roi->uy;
530 
531     index = (int *) ivector_alloc(0, (ux-lx)*(uy-ly));
532 
533     row0 = &(im->array.float_v[ly][lx]);
534     ishift = (ux-lx)/2.0;
535     for (i = ly, k=0; i < uy; ++i)
536     {
537         row = im->array.float_v[i];
538         for (j = lx; j < ishift; ++j)
539             if (row[j]>0.0) index[k++] = &(im->array.float_v[-i-1][j+ishift])
540                                               - row0;
541         for (j = ishift; j < ux; ++j)
542             if (row[j]>0.0) index[k++] = &(im->array.float_v[-i-1][j-ishift])
543                                               - row0;
544     }
545 
546     new_index = vector_alloc(k, int_v);
547 
548     for (i=0;i<k;i++)
549         new_index->el.int_v[i] = index[i];
550 
551     ivector_free(index, 0);
552 
553     return (new_index);
554 }
555 */
556 
557 /*
558 double fast_dot_prod(Imrect *im1,Imrect *im2,int type1,int type2,Vector *findex)
559 {
560    register double dp=0.0,dp1,dp2,dp3,dp4,dp5,dp6,dp7,dp8;
561    register int *offset;
562    register float *grey1,*grey2;
563    register int i, n, offint;
564 
565    grey1 = &(im1->array.float_v[im1->region->ly][im1->region->lx]);
566    grey2 = &(im2->array.float_v[im2->region->ly][im2->region->lx]);
567    offset = (int *)findex->el.int_v;
568    n = findex->n-8;
569 
570    for (i=0;i<n;)
571    {
572        offint = offset[i++];
573        dp1 = *(grey1 + offint)
574            * *(grey2 + offint);
575        offint = offset[i++];
576        dp2 = *(grey1 + offint)
577            * *(grey2 + offint);
578        offint = offset[i++];
579        dp3 = *(grey1 + offint)
580            * *(grey2 + offint);
581        offint = offset[i++];
582        dp4 = *(grey1 + offint)
583            * *(grey2 + offint);
584        offint = offset[i++];
585        dp5 = *(grey1 + offint)
586            * *(grey2 + offint);
587        offint = offset[i++];
588        dp6 = *(grey1 + offint)
589            * *(grey2 + offint);
590        offint = offset[i++];
591        dp7 = *(grey1 + offint)
592            * *(grey2 + offint);
593        offint = offset[i++];
594        dp8 = *(grey1 + offint)
595            * *(grey2 + offint);
596 
597        dp += dp1 +dp2 +dp3 + dp4 +dp5 + dp6 + dp7 + dp8;
598    }
599    for (;i<n+8;i++)
600    {
601        offint = offset[i++];
602        dp += *(grey1 + offint)
603            * *(grey2 + offint);
604    }
605 
606    return(dp);
607 }
608 */
609 
610 /*
611 double fast_dot_prod2(Imrect *im1,Imrect *im2,int type1,int type2,Vector *findex1,Vector *findex2)
612 {
613    register double dp=0.0,dp1,dp2,dp3,dp4,dp5,dp6,dp7,dp8;
614    register int *offset1;
615    register int *offset2;
616    register float *grey1,*grey2;
617    register int i,n;
618 
619    grey1 = &(im1->array.float_v[im1->region->ly][im1->region->lx]);
620    grey2 = &(im2->array.float_v[im2->region->ly][im2->region->lx]);
621    offset1 = (int *)findex1->el.int_v;
622    offset2 = (int *)findex2->el.int_v;
623    n = findex1->n-8;
624 
625    for (i=0;i<n;)
626    {
627        dp1 = *(grey1 + offset1[i])
628            * *(grey2 + offset2[i++]);
629        dp2 = *(grey1 + offset1[i])
630            * *(grey2 + offset2[i++]);
631        dp3 = *(grey1 + offset1[i])
632            * *(grey2 + offset2[i++]);
633        dp4 = *(grey1 + offset1[i])
634            * *(grey2 + offset2[i++]);
635        dp5 = *(grey1 + offset1[i])
636            * *(grey2 + offset2[i++]);
637        dp6 = *(grey1 + offset1[i])
638            * *(grey2 + offset2[i++]);
639        dp7 = *(grey1 + offset1[i])
640            * *(grey2 + offset2[i++]);
641        dp8 = *(grey1 + offset1[i])
642            * *(grey2 + offset2[i++]);
643        dp += dp1 + dp2 + dp3 + dp4 + dp5 + dp6 + dp7 + dp8;
644    }
645 
646    for (;i<n+8;i++)
647    {
648        dp += *(grey1 + offset1[i])
649            * *(grey2 + offset2[i]);
650    }
651 
652    return(dp);
653 }
654 */
655 
656 double dot_product(Imrect *im1, Imrect *im2, int type1, int type2)
657 {
658     double dp = 0.0;
659     Vector  *index;
660     Imregion *roi;
661     float *row1, *row2;
662     int lx, ux, ly, uy;
663     int i, j;
664 
665     if (im1 == NULL || im2 == NULL)
666         return (-1.0);
667     if (type1!=type2) return(-1.0);
668 
669     if ((index = prop_get(im1->props,FAST_INDEX1)))
670     {
671 /*  remove next bit until finished NAT 15/5/95
672         dp = fast_dot_prod(im1,im2,type1,type2,index);
673 */
674         return(dp);
675     }
676     else if ((index = prop_get(im2->props,FAST_INDEX1)))
677     {
678 /*  remove next bit until finished NAT 15/5/95
679         dp = fast_dot_prod(im1,im2,type1,type2,index);
680 */
681         return(dp);
682     }
683 
684     roi = roi_inter(im1->region, im2->region);
685     if (roi == NULL)
686         return (-1.0);
687     lx = roi->lx;
688     ux = roi->ux;
689     ly = roi->ly;
690     uy = roi->uy;
691 
692     for (i = ly; i < uy; ++i)
693     {
694         row1 = IM_ROW(im1,i);
695         row2 = IM_ROW(im2,i);
696         for (j = lx; j < ux; ++j)
697             dp += row1[j]*row2[j];
698     }
699 
700     rfree(roi);
701     return (dp);
702 }
703 
704 double dot_product2(Imrect *im1, Imrect *im2, int type1, int type2)
705 {
706     double dp = 0.0;
707     Imregion *roi;
708     Vector *index1,*index2;
709     float *row1, *row2;
710     int lx, ux, ly, uy;
711     int ishift;
712     int i, j;
713 
714     if (im1 == NULL || im2 == NULL)
715         return (-1.0);
716     if (type1!=type2) return(-1.0);
717 
718     if ((index1 = prop_get(im1->props,FAST_INDEX1))
719      && (index2 = prop_get(im1->props,FAST_INDEX2)))
720     {
721 /* remove until ready NAT 16/5/95
722         dp = fast_dot_prod2(im1,im2,type1,type2,index1,index2);
723 */
724         return(dp);
725     }
726     else if ((index1 = prop_get(im2->props,FAST_INDEX1))
727           && (index2 = prop_get(im2->props,FAST_INDEX2)))
728     {
729 /*        dp = fast_dot_prod2(im2,im1,type1,type2,index1,index2);
730 */
731         return(dp);
732     }
733 
734     roi = roi_inter(im1->region, im2->region);
735     if (roi == NULL)
736         return (-1.0);
737     lx = roi->lx;
738     ux = roi->ux;
739     ly = roi->ly;
740     uy = roi->uy;
741     ishift = (int)((ux-lx)/2.0);
742     for (i = ly; i < uy; ++i)
743     {
744         row1 = IM_ROW(im1,i);
745         row2 = IM_ROW(im2,-i-1);
746         for (j = lx; j < ishift; ++j)
747             dp += row1[j]*row2[j+ishift];
748         for (j = ishift; j < ux; ++j)
749             dp += row1[j]*row2[j-ishift];
750     }
751 
752     rfree(roi);
753     return (dp);
754 }
755 
756 /* ---------- */
757 

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