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

Linux Cross Reference
Tina6/tina-libs/tina/vision/visPgh_histfuncs.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: visPgh_histfuncs.c $
 27  * Date    :  $Date: 2012/06/22 10:00 $
 28  * Version :  $Revision: 1.8 $
 29  *
 30  * Author  : Legacy TINA modified NAT
 31  */
 32 /**
 33  * @file Wrapper routines for building and comparing Pairwise Geometric Histogrammes.
 34  * @brief Includes normalisation and comparison routines for fixed size histogrammes.
 35 */
 36 #include "visPgh_histfuncs.h"
 37 
 38 #include <stdio.h>
 39 #include <math.h>
 40 #include <tina/sys/sysDef.h>
 41 #include <tina/sys/sysPro.h>
 42 #include <tina/math/mathDef.h>
 43 #include <tina/math/mathPro.h>
 44 #include <tina/image/img_GenPro.h>
 45 #include <tina/geometry/geom_LinePro.h>
 46 #include <tina/vision/visPgh_hist.h>
 47 #include <tina/vision/visPgh_var.h>
 48 
 49 
 50 int count=1;
 51 /* ---------- Functions ---------- */
 52 
 53 void sqr_root_pairwise(Imrect *hist)
 54 {
 55     int     i,j, min_x, min_y, max_x, max_y;
 56     float  *hist_data;
 57 
 58     min_x = hist->region->lx;
 59     min_y = hist->region->ly;
 60     max_x = hist->region->ux;
 61     max_y = hist->region->uy;
 62 
 63 /* Each histogram entry is replaced by its square root */
 64 
 65     for(i=min_y;i<max_y;i++)
 66     {
 67         hist_data  = &IM_FLOAT( hist, i, min_x );
 68         for(j=min_x;j<max_x;j+=4)
 69         {
 70 /*
 71             IM_PIX_GET(hist,i,j,hist_data);
 72             if (hist_data < 0.0) hist_data = 0.0;
 73             IM_PIX_SET(hist,i,j, (double)sqrt((double)hist_data));
 74 */
 75             if (*hist_data < 0.0) *hist_data = 0.0;
 76             else *hist_data = (float)sqrt((double)*hist_data);
 77             hist_data++;
 78             if (*hist_data < 0.0) *hist_data = 0.0;
 79             else *hist_data = (float)sqrt((double)*hist_data);
 80             hist_data++;
 81             if (*hist_data < 0.0) *hist_data = 0.0;
 82             else *hist_data = (float)sqrt((double)*hist_data);
 83             hist_data++;
 84             if (*hist_data < 0.0) *hist_data = 0.0;
 85             else *hist_data = (float)sqrt((double)*hist_data);
 86             hist_data++;
 87         }
 88     }
 89 }
 90 
 91 //TODO since pghs currently directed and many of an object's edges will be at the edge of the model, for many pghs; top/bottom half will be blank,
 92 //     therefore, tag pgh dist bins if corresponding pgh angle row is blank, then just skip these rows for any calculations asociated with pghs (assuming 0.0 bin vals)
 93 double dot_pairwise( Imrect *im1, Imrect *im2)
 94 {
 95     register float dp=0.0;
 96     register float *grey1, *grey2;
 97     register int i, j, lx, ux, ly, uy;
 98 
 99     lx = im1->region->lx;
100     ux = im1->region->ux;
101     ly = im1->region->ly;
102     uy = im1->region->uy;
103  
104     for ( i = ly; i < uy; i++ )
105     {
106         grey1 = &IM_FLOAT( im1, i, lx );
107         grey2 = &IM_FLOAT( im2, i, lx );
108         for (j=lx;j<ux;j+=4)
109         {
110             dp += *grey1 * *grey2;
111             grey1++; grey2++;
112             dp += *grey1 * *grey2;
113             grey1++; grey2++;
114             dp += *grey1 * *grey2;
115             grey1++; grey2++;
116             dp += *grey1 * *grey2;
117             grey1++; grey2++;
118         }
119     }
120 
121     return ( dp );
122 }
123 
124 /* ---------- */
125 
126 void sqr_root_and_normalize_pairwise(Imrect *hist)
127 {
128 int      i,j, min_x, min_y, max_x, max_y;
129 double   sum=0.0;
130 float *gl;
131 
132         /* (1) Find the sum of the histogram entries.
133            (2) Divide each histogram entry by this value. 
134            (3) Replace each histogram entry by its sqr root.
135                This is done for statistical reasons. */
136 
137     min_x = hist->region->lx;
138     min_y = hist->region->ly;
139     max_x = hist->region->ux;
140     max_y = hist->region->uy;
141 
142     for(i=min_y;i<max_y;i++)
143     {
144         gl  = &IM_FLOAT( hist, i, min_x );
145         for(j=min_x;j<max_x;j+=4)
146         {
147 /*
148             IM_PIX_GET(hist,i,j,gl);
149             sum += gl;
150 */
151             sum+= *gl;
152             gl ++;
153             sum+= *gl;
154             gl ++;
155             sum+= *gl;
156             gl ++;
157             sum+= *gl;
158             gl ++;
159         }
160     }
161 
162     if (sum!=0.0)
163     {
164         for(i=min_y;i<max_y;i++)
165         {
166             gl  = &IM_FLOAT( hist, i, min_x );
167             for(j=min_x;j<max_x;j+=4)
168             {
169 /*
170                 IM_PIX_GET(hist,i,j,gl);
171                 IM_PIX_SET(hist,i,j,gl/(double)sum); 
172 */
173                 *gl = *gl/sum;
174                 gl++;
175                 *gl = *gl/sum;
176                 gl++;
177                 *gl = *gl/sum;
178                 gl++;
179                 *gl = *gl/sum;
180                 gl++;
181             }
182         }
183         sqr_root_pairwise(hist);
184     }
185 }
186 
187 /* ---------- Functions ---------- */
188 
189 
190 double pairwise_length( Imrect *hist )
191 {
192     int     i,j, min_x, min_y, max_x, max_y;
193     float  *gl;
194     double  length = 0.0;
195 
196     min_x = hist->region->lx;
197     min_y = hist->region->ly;
198     max_x = hist->region->ux;
199     max_y = hist->region->uy;
200 
201     for(i=min_y;i<max_y;i++)
202     {
203         gl  = &IM_FLOAT( hist, i, min_x );
204         for(j=min_x;j<max_x;j+=4)
205         {
206 /*
207             IM_PIX_GET(hist,i,j,gl);
208             length += gl*gl;
209 */
210             length += *gl * *gl;
211             gl++;
212             length += *gl * *gl;
213             gl++;
214             length += *gl * *gl;
215             gl++;
216             length += *gl * *gl;
217             gl++;
218         }
219     }
220     length = sqrt( fabs(length) );
221     return (length);
222 }
223 
224 
225 Imrect *flip_pairwise( Imrect *hist )
226 {
227     int     i,j, min_x, min_y, max_x, max_y, ishift;
228     float  *gl, *go1, *go2; 
229     Imrect *flip_hist;
230     Imregion *roi;
231 
232     roi = hist->region;
233     min_x = roi->lx;
234     max_x = roi->ux;
235     min_y = roi->ly;
236     max_y = roi->uy;
237     ishift = (int)((max_x-min_x)/2.0);
238 
239     flip_hist = im_pairs_alloc();
240 
241     for(i=min_y;i<max_y;i++)
242     {
243         gl  = &IM_FLOAT( hist, i, min_x );
244         go1  = &IM_FLOAT( flip_hist, -i-1, min_x+ishift );
245         go2  = &IM_FLOAT( flip_hist, -i-1, min_x );
246         for(j=min_x;j<ishift;j+=4)
247         {
248             *go1 = *gl;
249              gl++; go1++;
250             *go1 = *gl;
251              gl++; go1++;
252             *go1 = *gl;
253              gl++; go1++;
254             *go1 = *gl;
255              gl++; go1++;
256         }
257         for(j=ishift;j<max_x;j+=4)
258         {
259             *go2 = *gl;
260              gl++; go2++;
261             *go2 = *gl;
262              gl++; go2++;
263             *go2 = *gl;
264              gl++; go2++;
265             *go2 = *gl;
266              gl++; go2++;
267         }
268     }
269 
270     flip_hist->props = proplist_copy(hist->props);
271 
272     return ( flip_hist );
273 }
274 
275 Imrect *normalize_pairwise( Imrect *hist )
276 {
277     int     i,j, min_x, min_y, max_x, max_y;
278     float  *gl, *go; 
279     double  length;
280     Imrect *norm_hist;
281     Imregion *roi;
282 
283     length = pairwise_length( hist );
284 
285     roi = hist->region;
286     min_x = roi->lx;
287     max_x = roi->ux;
288     min_y = roi->ly;
289     max_y = roi->uy;
290 
291 
292     norm_hist = im_pairs_alloc();
293 
294     for(i=min_y;i<max_y;i++)
295     {
296         gl  = &IM_FLOAT( hist, i, min_x );
297         go  = &IM_FLOAT( norm_hist, i, min_x );
298         for(j=min_x;j<max_x;j+=4)
299         {
300 /*
301             IM_PIX_GET(hist,i,j,gl);
302             IM_PIX_SET(norm_hist,i,j,gl/length); 
303 */
304             *go = *gl/length;
305              gl++;
306              go++;
307             *go = *gl/length;
308              gl++;
309              go++;
310             *go = *gl/length;
311              gl++;
312              go++;
313             *go = *gl/length;
314              gl++;
315              go++;
316         }
317     }
318 
319     norm_hist->props = proplist_copy(hist->props);
320 
321     return ( norm_hist );
322 }
323 
324 Imrect *normalize_pairwise_param( Imrect *hist, double param )
325 /* I dont think this has ever been used NAT 15/5/09 */
326 {
327     int     i,j, min_x, min_y, max_x, max_y;
328     float  *gl;
329     Imrect *norm_hist;
330     Imregion *roi;
331 
332     roi = hist->region;
333     min_x = roi->lx;
334     max_x = roi->ux;
335     min_y = roi->ly;
336     max_y = roi->uy;
337 
338     for(i=min_y;i<max_y;i++)
339     {
340         gl  = &IM_FLOAT( hist, i, min_x ); 
341         for(j=min_x;j<max_x;j+=4)
342         {
343 /*
344             IM_PIX_GET(hist,i,j,gl);
345             IM_PIX_SET(hist,i,j,gl/param); 
346 */
347             *gl = *gl/param;
348             gl++;
349             *gl = *gl/param;
350             gl++;
351             *gl = *gl/param;
352             gl++;
353             *gl = *gl/param;
354             gl++;
355         }
356     }
357     return ( hist );
358 }
359 
360 Imrect *normalize_pairwise_param2( Imrect *hist, double param )
361 {
362     int      i,j, min_x, min_y, max_x, max_y;
363     float  *gl,*go;
364     double recparam = 1.0/param;
365     Imrect *norm_hist;
366     Imregion *roi;
367 
368     roi = hist->region;
369     min_x = roi->lx;
370     max_x = roi->ux;
371     min_y = roi->ly;
372     max_y = roi->uy;
373 
374     norm_hist = im_alloc( hist->height, hist->width, roi, float_v);
375 
376     for(i=min_y;i<max_y;i++)
377     {
378         gl  = &IM_FLOAT( hist, i, min_x ); 
379         go  = &IM_FLOAT( norm_hist, i, min_x ); 
380         for(j=min_x;j<max_x;j+=4)
381         {
382 /*
383             IM_PIX_GET(hist,i,j,gl);
384             IM_PIX_SET(norm_hist,i,j,gl/param); 
385 */
386             *go = *gl *recparam;
387             gl++;
388             go++;
389             *go = *gl *recparam;
390             gl++;
391             go++;
392             *go = *gl *recparam;
393             gl++;
394             go++;
395             *go = *gl *recparam;
396             gl++;
397             go++;
398         }
399     }
400     return ( norm_hist );
401 }
402 
403 Imrect *diff_pairwise( Imrect *hist1, Imrect *hist2 )
404 {
405     int      i,j, min_x, min_y, max_x, max_y;
406     float  *g1,*g2,*go;
407     Imrect *out_hist;
408     Imregion *roi;
409 
410     roi = hist1->region;
411     min_x = roi->lx;
412     max_x = roi->ux;
413     min_y = roi->ly;
414     max_y = roi->uy;
415 
416     out_hist = im_alloc( hist1->height, hist1->width, roi, float_v);
417 
418     for(i=min_y;i<max_y;i++)
419     {
420         g1  = &IM_FLOAT( hist1, i, min_x );
421         g2  = &IM_FLOAT( hist2, i, min_x );
422         go  = &IM_FLOAT( out_hist, i, min_x );
423         for(j=min_x;j<max_x;j+=4)
424         {
425             *go = *g1 - *g2;
426             g1++;
427             g2++;
428             go++;
429             *go = *g1 - *g2;
430             g1++;
431             g2++;
432             go++;
433            *go = *g1 - *g2;
434             g1++;
435             g2++;
436             go++;
437             *go = *g1 - *g2;
438             g1++;
439             g2++;
440             go++;
441         }
442     }
443     return ( out_hist );
444 }
445 
446 Imrect *sum_pairwise( Imrect *hist1, Imrect *hist2 )
447 {
448     int      i,j, min_x, min_y, max_x, max_y;
449     float  *g1,*g2,*go;
450     Imrect *out_hist;
451     Imregion *roi;
452 
453     roi = hist1->region;
454     min_x = roi->lx;
455     max_x = roi->ux;
456     min_y = roi->ly;
457     max_y = roi->uy;
458 
459     out_hist = im_alloc( hist1->height, hist1->width, roi, float_v);
460 
461     for(i=min_y;i<max_y;i++)
462     {
463         g1  = &IM_FLOAT( hist1, i, min_x );
464         g2  = &IM_FLOAT( hist2, i, min_x );
465         go  = &IM_FLOAT( out_hist, i, min_x );
466         for(j=min_x;j<max_x;j+=4)
467         {
468             *go = *g1 + *g2;
469             g1++;
470             g2++;
471             go++;
472             *go = *g1 + *g2;
473             g1++;
474             g2++;
475             go++;
476             *go = *g1 + *g2;
477             g1++;
478             g2++;
479             go++;
480             *go = *g1 + *g2;
481             g1++;
482             g2++;
483             go++;
484         }
485     }
486     return ( out_hist );
487 }
488 
489 double tot_pairwise(Imrect *hist)
490 {
491     int     i,j, min_x, min_y, max_x, max_y;
492     float  *gl;
493     double  sum = 0.0;
494 
495     min_x = hist->region->lx;
496     min_y = hist->region->ly;
497     max_x = hist->region->ux;
498     max_y = hist->region->uy;
499 
500     for(i=min_y;i<max_y;i++)
501     {
502         gl  = &IM_FLOAT( hist, i, min_x );
503         for(j=min_x;j<max_x;j+=4)
504         {
505 /*
506             IM_PIX_GET(hist,i,j,gl);
507             sum += gl;
508 */
509             sum += *gl;
510             gl++;
511             sum += *gl;
512             gl++;
513             sum += *gl;
514             gl++;
515             sum += *gl;
516             gl++;
517         }
518     }
519     return (sum);
520 }
521 
522 /* ---------- */
523 //note: same as build normalised_pairwise - except for call to above funcs - could condense
524 
525 Imrect *build_pairwise(Line2 *ref_line, Model_poly_header *ref_model, List *line_list,
526                        double window_r)
527 {
528     Imrect   *pairwise;
529     List     *ptr;
530     Line2    *lptr, *clip_line;
531     Hist_ref *ref;
532     double    l, l0;
533     int type = -1;
534     int i,j, *line_id, *parent_id, prev=-9999;
535 
536         /* (1) Allocate space for the histogram.
537            (2) Compile the histogram from the ref line and
538                the line list.
539            (3) Sqrt the hist to get the correct statistics.
540            (4) Add a Hist_ref to the histograms props list.
541         */
542     pairwise = im_pairs_alloc();
543     if (pairwise==NULL)
544     {
545         format("Error: Not space for pairwise histogram.\n");
546         return NULL;
547     }
548     for(i = pairwise->region->ly; i < pairwise->region->uy; i++)
549         for(j = pairwise->region->lx; j < pairwise->region->ux; j++)
550             IM_PIX_SET(pairwise,i, j, 0.0);
551 
552     for(ptr=line_list;ptr!=NULL;ptr=ptr->next)
553     {
554         lptr=ptr->to;
555 /*
556         if ( (vec2_dist(ref_line->p, lptr->p) < ( 2.0 * window_r )) )
557               if ( vec2_perp_dist2( ref_line->p1, ref_line->p2, lptr->p ) < window_r )
558 */
559         clip_line = line2_circ_inter(lptr, ref_line->p, 2.0*window_r);
560         if (clip_line!=NULL)
561         {
562             compare_lines(pairwise, ref_line, clip_line, &type);                
563             line2_free(clip_line);
564         }
565     }
566     l0 = vec2_dist(ref_line->p1, ref_line->p2);
567     l  = tot_pairwise(pairwise)/l0;
568     sqr_root_pairwise(pairwise);
569     ref = ralloc(sizeof(Hist_ref));
570     if (ref==NULL)
571     {
572         format("Error: No space for Hist_ref.\n");
573         im_free(pairwise);
574         return NULL;
575     }
576     ref->model = ref_model;
577     ref->matches = NULL;
578     ref->line_seg = ref_line; 
579     ref->l0 = l0;
580     ref->l = l;
581     ref->type = type;
582     ref->scale_factor = 1.0;
583 
584     pairwise->props = proplist_add(pairwise->props, ref, HIST_REF_TYPE, hist_ref_free);
585 
586     return (pairwise);
587 }
588 
589 
590 /* ---------- */
591 
592 Imrect *build_normalized_pairwise(Line2 *ref_line, Model_poly_header *ref_model,
593                                   List *line_list, double window_r)
594 {
595     Imrect   *pairwise;
596     List     *ptr;
597     Line2    *lptr, *clip_line;
598     Hist_ref *ref;
599     double    l, l0;
600     int type = -1, i, j, *line_id, *parent_id, prev=-999, last=-999;//, line_c=0;
601     double length;
602     Vector *findex;
603 
604 
605         /* (1) Allocate space for the histogram.
606            (2) Compile the histogram from the ref line and
607                the line list.
608            (3) Normalize the histogram.
609            (4) Sqrt the hist to get the correct statistics.
610            (5) Add a Hist_ref to the histograms props list.
611         */
612 
613 //printf("refl2 p1 = (%f %f) (%f %f)\n", vec2_x(ref_line->p1), vec2_y(ref_line->p1), vec2_x(ref_line->p2), vec2_y(ref_line->p2));
614 
615     pairwise = im_pairs_alloc();
616     if (pairwise==NULL)
617     {
618         printf("Error: No space for pairwise histogram.\n");
619         return NULL;
620     }
621 
622         for(i = pairwise->region->ly; i < pairwise->region->uy; i++)
623                 for(j = pairwise->region->lx; j < pairwise->region->ux; j++)
624                         IM_PIX_SET(pairwise,i, j, 0.0);
625 
626     for(ptr=line_list;ptr!=NULL;ptr=ptr->next)
627     {
628         lptr=ptr->to;
629 /*
630             if ( (vec2_dist(ref_line->p, lptr->p) < ( 2.0 * window_r )) )
631             {
632                if ( vec2_perp_dist2( ref_line->p1, ref_line->p2, lptr->p ) < window_r )
633 */
634         clip_line = line2_circ_inter(lptr, ref_line->p, 2.0*window_r);
635         if (clip_line!=NULL && clip_line->length > 0.5)
636         {
637             compare_lines(pairwise, ref_line, clip_line, &type);
638             line2_free(clip_line);
639         }
640     }
641 
642     l0 = vec2_dist(ref_line->p1, ref_line->p2);
643     l = tot_pairwise(pairwise)/l0;
644     if( type == ROTATE ) pairwise = balance_extreme_angs( pairwise );
645 
646 
647     sqr_root_and_normalize_pairwise(pairwise);
648 
649     ref = ralloc(sizeof(Hist_ref));
650     if (ref==NULL)
651     {
652         printf("Error: No space for Hist_ref.\n"); //(format("Error: No space for Hist_ref.\n");)
653         im_free(pairwise);
654         return NULL;
655     }
656     ref->model = ref_model;
657     ref->matches = NULL;
658     ref->line_seg = ref_line;
659     ref->l0 = l0;
660     ref->l = l;
661     ref->type = type;
662     ref->scale_factor = 1.0;
663     pairwise->props = proplist_add(pairwise->props, ref, HIST_REF_TYPE, hist_ref_free);
664 /*
665     pairwise->label = ref_line->label;
666 */
667 
668 /*
669 findex = nonzero_index(pairwise);
670 if(findex==NULL)printf("ERROR findex = NULL\n");
671 pairwise->props = proplist_add(pairwise->props, findex, FAST_INDEX1, vector_free );
672 */
673 
674     return (pairwise);
675 }
676 
677 Imrect *balance_extreme_angs( Imrect *pgh )
678 {
679     int ramp_shift, lx, ux, ly, uy, i, j; 
680     Imregion *roi;
681     double data, data2, dist_ramp, dbin_size;
682     Pairs_hist_def *hist_def;
683 
684     hist_def = pairs_hist_def_get();
685 
686     roi = pgh->region;
687     lx = roi->lx; 
688     ux = roi->ux; 
689     ly = roi->ly;
690     uy = roi->uy;
691 
692     dist_ramp = hist_def->dist_ramp;
693     dbin_size = hist_def->dbin_size;
694     ramp_shift = (int)ceil( dist_ramp / dbin_size );
695 
696     for ( i = lx; i < ux; i++ )
697     {
698         if ( (i <= lx + ramp_shift) || (i >= ux - ( ramp_shift + 1 )) )
699         {
700             for ( j = ly; j < 0; j++ )
701             {
702                 IM_PIX_GET( pgh, j, i, data );
703                 IM_PIX_GET( pgh, (j * -1)-1, i, data2 );
704 
705                 if ( data != 0.0 || data2 != 0.0 )
706                 {
707                     IM_PIX_SET( pgh, j, i, (data+data2)/2.0 );
708                     IM_PIX_SET( pgh, (j * -1)-1, i, (data+data2)/2.0 );
709                 }
710             }
711         } 
712     }
713     return(pgh);
714 }
715 
716 /* ---------- */
717 
718 void hist_ref_free(Hist_ref *hist_ref, int type)
719 {
720 
721         /* (1) Free the matches list. (Only the links,
722                not the data which they point to).
723            (2) Free the copied line seg.
724            (3) Free the hist_ref itself. */
725 
726     if (type!=HIST_REF_TYPE)
727     {
728         format("Error: Ilegal use of 'hist_ref_free'.\n");
729         return;
730     }
731 
732     list_rm_links(hist_ref->matches);
733 /*    line2_free(hist_ref->line_seg);
734       not to be done here NAT
735 */
736     rfree(hist_ref);
737 
738 }
739 
740 /* ---------- */
741 
742 Imrect *pairs_build_norm_hist_scale(Line2 *ref_line, Model_poly_header *ref_model, List *line_list, double window_r, double scale)
743 {
744     Imrect   *pairwise;
745     List     *ptr;
746     Line2    *lptr, *scaled_ref, *scaled_line, *clip_line;
747     Hist_ref *ref;
748     double    l, l0;
749     int type = -1;
750     Transform2 trans;
751     int i,j;
752         
753 
754         /* (1) Allocate space for the histogram.
755            (2) Compile the histogram from the scaled ref line and
756                the scaled line list.
757            (3) Normalize the histogram.
758            (4) Sqrt the hist to get the correct statistics.
759            (5) Add a Hist_ref to the histograms props list.
760         */
761 
762     trans.R =  rot2_with_scale(0.0, scale);
763     trans.t =  vec2_zero();
764 
765     scaled_ref = line2_copy(ref_line);
766     line2_transform(scaled_ref, trans);
767 
768     pairwise = im_pairs_alloc();
769     if (pairwise==NULL)
770     {
771         printf("Error: Not space for pairwise histogram.\n");
772         return NULL;
773     }
774 
775 
776     for(i = pairwise->region->ly; i < pairwise->region->uy; i++)
777         for(j = pairwise->region->lx; j < pairwise->region->ux; j++)
778             IM_PIX_SET(pairwise,i, j, 0.0);
779                         
780     for(ptr=line_list;ptr!=NULL;ptr=ptr->next)
781     {
782         lptr= (Line2 *) ptr->to;
783         scaled_line = line2_copy(lptr);
784         line2_transform(scaled_line, trans);
785 
786         clip_line = line2_circ_inter(scaled_line, scaled_ref->p, 2.0*window_r);
787         if (clip_line!=NULL)
788         {
789            compare_lines(pairwise, scaled_ref, clip_line, &type);
790            line2_free(clip_line);
791         }
792         line2_free(scaled_line);
793     }
794 
795     l0 = vec2_dist(scaled_ref->p1, scaled_ref->p2);
796     l = tot_pairwise(pairwise)/l0;
797 
798 
799     sqr_root_and_normalize_pairwise(pairwise);
800 
801     ref = ralloc(sizeof(Hist_ref));
802     if (ref==NULL)
803     {
804         printf("Error: No space for Hist_ref.\n");
805         im_free(pairwise);
806         return NULL;
807     }
808 
809     line2_free(scaled_ref);
810 
811     ref->model = ref_model;
812     ref->matches = NULL;
813     ref->line_seg = ref_line;
814     ref->l0 = l0;
815     ref->l = l;
816     ref->type = type;
817         ref->scale_factor = scale;
818 
819     pairwise->props = proplist_add(pairwise->props, ref, HIST_REF_TYPE, 
820                                                                hist_ref_free);
821    /*  No longer supported APA 19/6/95
822        pairwise->label = ref_line->label; */
823 
824     return (pairwise);
825 
826 }
827 
828 /* ---------- */
829 
830 Hist_ref *hist_ref_get(Imrect *hist)
831 {
832     Hist_ref *hist_ref;
833 
834     hist_ref = prop_get(hist->props, HIST_REF_TYPE);
835 
836     if (hist_ref==NULL)
837     {
838         printf("Error: Histogram has no hist_ref.\n");
839         return NULL;
840     }
841 
842     return hist_ref;
843 }
844 
845 /* ---------- */
846 
847 Line2 *ref_line_from_hist(Imrect *hist)
848 {
849     Hist_ref  *hist_ref;
850     Line2     *ref_line;
851 
852     if ((hist_ref = hist_ref_get(hist))==NULL) return NULL;
853 
854     if ((ref_line = hist_ref->line_seg)==NULL)
855     {
856         printf("Error: Histogram has no reference line.\n");
857         return NULL;
858     }
859 
860     return ref_line;
861 }    
862 
863 /* ---------- */
864 
865 List *model_geom_from_hist(Imrect *hist)
866 {
867     Hist_ref          *hist_ref;
868     Model_poly_header *model_poly_head;
869     List              *model_geom;
870 
871     if ((hist_ref = hist_ref_get(hist))==NULL) return NULL;
872 
873     if ((model_poly_head = hist_ref->model)==NULL)
874     {
875         printf("Error: Histogram has no model_poly_header.\n");
876         return NULL;
877     }
878 
879     if ((model_geom = model_poly_head->model_poly_data)==NULL)
880     {
881         printf("Error: Histogram has no model geometry.\n");
882         return NULL;
883     }
884 
885     return model_geom;
886 }
887 
888 /* ---------- */
889 
890 List *matches_list_from_hist(Imrect *hist)
891 {
892     Hist_ref          *hist_ref;
893     List              *matches_list;
894 
895     if ((hist_ref = hist_ref_get(hist))==NULL) return NULL;
896 
897     if ((matches_list = hist_ref->matches)==NULL)
898     {
899         printf("Error: Histogram has no matches.\n");
900         return NULL;
901     }
902 
903     return matches_list;
904 }
905 
906 /* ---------- */
907 
908 Hist_ref *hist_ref_copy(Hist_ref *hist_ref)
909 {
910     Hist_ref *copy;
911 
912     copy = ralloc(sizeof(Hist_ref));
913 
914     copy->line_seg = hist_ref->line_seg;
915     copy->model = hist_ref->model;
916     copy->l0 = hist_ref->l0;
917     copy->l = hist_ref->l;
918     copy->scale_factor = hist_ref->scale_factor;
919     copy->type = hist_ref->type;
920 
921     return copy; 
922 }
923 
924 /* ---------- */
925 
926 double dot_product(Imrect *im1, Imrect *im2, int type1, int type2)
927 {
928     double dp = 0.0;
929     Vector  *index;
930     Imregion *roi;
931     float *row1, *row2;
932     int lx, ux, ly, uy;
933     int i, j;
934  
935     if (im1 == NULL || im2 == NULL)
936         return (-1.0);
937     if (type1!=type2) return(-1.0); 
938 
939 
940     if ((index = prop_get(im1->props,FAST_INDEX1))) 
941     {
942         dp = fast_dot_prod(im1,im2,type1,type2,index);
943         return(dp);
944     }
945     else if ((index = prop_get(im2->props,FAST_INDEX1)))
946     {
947         dp = fast_dot_prod(im1,im2,type1,type2,index);
948         return(dp);
949     }
950 
951     roi = roi_inter(im1->region, im2->region);
952     if (roi == NULL)
953         return (-1.0);
954     lx = roi->lx;
955     ux = roi->ux;
956     ly = roi->ly;
957     uy = roi->uy;
958 
959     for (i = ly; i < uy; ++i)
960     {
961         row1 = IM_ROW(im1,i);
962         row2 = IM_ROW(im2,i);
963         for (j = lx; j < ux; ++j)
964         {
965             dp += row1[j]*row2[j];
966         }
967     }
968     rfree(roi);
969     return (dp);
970 }
971 
972 /*
973 Vector *nonzero_index( Imrect *im )
974 {
975     Imregion *roi;
976     float *row, *row0;
977     int *index;
978     Vector *new_index, *vector_alloc();
979     int lx, ux, ly, uy;
980     int width, height;
981     int i, j, k;
982 
983     if ( im == NULL ) return( NULL );
984 
985     roi = im->region;
986     if ( roi == NULL ) return( NULL );
987     lx = roi->lx;
988     ux = roi->ux;
989     ly = roi->ly;
990     uy = roi->uy;
991 
992     index = (int *) ivector_alloc( 0, (ux-lx)*(uy-ly) );
993     row0 = &IM_FLOAT( im, ly, lx );
994 
995     for ( i = ly, k = 0; i < uy; ++i )
996     {
997         row =  &IM_FLOAT( im, i, 0 );
998 
999         for ( j = lx; j < ux; ++j )
1000         {
1001             if ( row[j] > 0.0 )
1002             {
1003                 index[k++] = &row[j] - row0;
1004             }
1005         }
1006     } 
1007     new_index = vector_alloc( k, int_v );
1008 
1009     for (i=0;i<k;i++) ((float*)new_index->data)[i] = index[i];
1010 
1011     ivector_free( index, 0 );
1012     return( new_index );
1013 }
1014 */ 
1015  
1016 double fast_dot_prod(Imrect *im1, Imrect *im2, int type1, int type2, Vector *findex )
1017 {
1018     register float dp=0.0;
1019     register int *offset;
1020     register float *grey1, *grey2;
1021     register int i, n, offint;
1022 
1023     grey1 = &IM_FLOAT( im1, im1->region->ly, im1->region->lx );
1024     grey2 = &IM_FLOAT( im2, im2->region->ly, im2->region->lx );
1025 
1026     offset = (int *)findex->data;
1027 
1028     n = findex->n-8;
1029 
1030     for (i=0;i<n;)
1031     {
1032         //offint = offset[i++];
1033         dp += *(grey1 + *offset) * *(grey2 + *offset);
1034         offset++;
1035         dp += *(grey1 + *offset) * *(grey2 + *offset);
1036         offset++;
1037         dp += *(grey1 + *offset) * *(grey2 + *offset);
1038         offset++;
1039         dp += *(grey1 + *offset) * *(grey2 + *offset);
1040         offset++;
1041         dp += *(grey1 + *offset) * *(grey2 + *offset);
1042         offset++;
1043         dp += *(grey1 + *offset) * *(grey2 + *offset);
1044         offset++;
1045         dp += *(grey1 + *offset) * *(grey2 + *offset);
1046         offset++;
1047         dp += *(grey1 + *offset) * *(grey2 + *offset);
1048         offset++;
1049     }
1050 
1051     for ( ;i<n+8;i++)
1052     {
1053         dp += *(grey1 + *offset) * *(grey2 + *offset);
1054     }
1055 
1056     return ( dp );
1057 }
1058 
1059 double dot_productCP(Imrect *im1, Imrect *im2, int type1, int type2)
1060 {
1061     double dp = 0.0;
1062     Vector  *index;
1063     Imregion *roi;
1064     float *row1, *row2;
1065     int lx, ux, ly, uy;
1066     int i, j;
1067 
1068     if (im1 == NULL || im2 == NULL)
1069         return (-1.0);
1070     if (type1!=type2) return(-1.0);
1071 
1072     if ((index = prop_get(im1->props,FAST_INDEX1)))
1073     {
1074 /*  remove next bit until finished NAT 15/5/95
1075         dp = fast_dot_prod(im1,im2,type1,type2,index);
1076 */
1077         return(dp);
1078     }
1079     else if ((index = prop_get(im2->props,FAST_INDEX1)))
1080     {
1081 /*  remove next bit until finished NAT 15/5/95
1082         dp = fast_dot_prod(im1,im2,type1,type2,index);
1083 */
1084         return(dp);
1085     }
1086 
1087     roi = roi_inter(im1->region, im2->region);
1088     if (roi == NULL)
1089         return (-1.0);
1090     lx = roi->lx;
1091     ux = roi->ux;
1092     ly = roi->ly;
1093     uy = roi->uy;
1094 
1095     for (i = ly; i < uy; ++i)
1096     {
1097         row1 = IM_ROW(im1,i);
1098         row2 = IM_ROW(im2,i);
1099         for (j = lx; j < ux; ++j)
1100             dp += row1[j]*row2[j];
1101     }
1102 
1103     rfree(roi);
1104     return (dp);
1105 }
1106 
1107 double dot_product2(Imrect *im1, Imrect *im2, int type1, int type2)
1108 {
1109     double dp = 0.0;
1110     Imregion *roi;
1111     Vector *index1,*index2;
1112     float *row1, *row2;
1113     int lx, ux, ly, uy;
1114     int ishift;
1115     int i, j;
1116 
1117     if (im1 == NULL || im2 == NULL)
1118         return (-1.0);
1119     if (type1!=type2) return(-1.0);
1120 
1121     if ((index1 = prop_get(im1->props,FAST_INDEX1))
1122      && (index2 = prop_get(im1->props,FAST_INDEX2)))
1123     {
1124 /* remove until ready NAT 16/5/95
1125         dp = fast_dot_prod2(im1,im2,type1,type2,index1,index2);
1126 */
1127         return(dp);
1128     }
1129     else if ((index1 = prop_get(im2->props,FAST_INDEX1))
1130           && (index2 = prop_get(im2->props,FAST_INDEX2)))
1131     {
1132 /*        dp = fast_dot_prod2(im2,im1,type1,type2,index1,index2);
1133 */
1134         return(dp);
1135     }
1136 
1137     roi = roi_inter(im1->region, im2->region);
1138     if (roi == NULL)
1139         return (-1.0);
1140     lx = roi->lx;
1141     ux = roi->ux;
1142     ly = roi->ly;
1143     uy = roi->uy;
1144     ishift = (int)((ux-lx)/2.0);
1145     for (i = ly; i < uy; ++i)
1146     {
1147         row1 = IM_ROW(im1,i);
1148         row2 = IM_ROW(im2,-i-1);
1149         for (j = lx; j < ishift; ++j)
1150             dp += row1[j]*row2[j+ishift];
1151         for (j = ishift; j < ux; ++j)
1152             dp += row1[j]*row2[j-ishift];
1153     }
1154 
1155     rfree(roi);
1156     return (dp);
1157 }
1158 
1159 /* ---------- */
1160 
1161 /* following maths from wolfram point-line distance 2d */
1162 
1163 double vec2_perp_dist2( Vec2 l1, Vec2 l2, Vec2 p ) 
1164 {
1165     double l1_x, l1_y, l2_x, l2_y, p_x, p_y, dist;
1166     
1167     l1_x = vec2_get_x( &l1 );
1168     l1_y = vec2_get_y( &l1 );
1169     l2_x = vec2_get_x( &l2 );
1170     l2_y = vec2_get_y( &l2 ); 
1171     p_x  = vec2_get_x( &p );
1172     p_y  = vec2_get_y( &p ); 
1173 
1174     dist =  ( l2_x - l1_x ) * ( l1_y - p_y );
1175     dist -= ( l1_x - p_x ) * ( l2_y - l1_y );
1176     dist /= sqrt( (( l2_x - l1_x )*( l2_x - l1_x )) + (( l2_y - l1_y )*( l2_y - l1_y )) );
1177 
1178     return ( fabs( dist ) );
1179 }
1180 
1181 

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