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

Linux Cross Reference
Tina5/tina-libs/tina/vision/visShape_hough2.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 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/visShape_hough2.c,v $
 23  * Date    :  $Date: 2008/12/09 00:15:22 $
 24  * Version :  $Revision: 1.6 $
 25  * CVS Id  :  $Id: visShape_hough2.c,v 1.6 2008/12/09 00:15:22 paul Exp $
 26  *
 27  * Author  : Legacy TINA
 28  *
 29  * Notes :
 30  *
 31  *********
 32 */
 33 
 34 #include "visShape_hough2.h"
 35 
 36 #if HAVE_CONFIG_H
 37   #include <config.h>
 38 #endif
 39 
 40 
 41 #include <stdio.h>
 42 #include <stdlib.h>
 43 #include <math.h>
 44 #include <tina/sys/sysDef.h>
 45 #include <tina/sys/sysPro.h>
 46 #include <tina/math/mathDef.h>
 47 #include <tina/image/imgDef.h>
 48 #include <tina/image/imgPro.h>
 49 #include <tina/geometry/geomDef.h>
 50 #include <tina/geometry/geomPro.h>
 51 #include <tina/vision/vis_ShapeDef.h>
 52 
 53 
 54 /*#include "extra.h"
 55 */
 56 
 57 /* ---------- Functions ---------- */
 58 
 59 Imrect *hough2_alloc(double min_x, double min_y, double max_x, double max_y, double binsize_x, double binsize_y, Vartype type)
 60 {
 61     Imrect       *hough2;
 62     Imregion      roi;
 63     Hough2_region *hough2_roi;
 64 
 65 /* Define the physical size of the hough space array. Need to use the 'ceil' 
 66    function because the hough space may not be exactly divisable by the bin
 67    size. */
 68 
 69     roi.lx = 0;
 70     roi.ux = (int)ceil((max_x-min_x)/binsize_x);
 71     roi.ly = 0;
 72     roi.uy = (int)ceil((max_y-min_y)/binsize_y);
 73   
 74     hough2 = im_alloc(roi.uy, roi.ux, &roi, type);
 75 
 76 /* Add the logical hough coords to the props list */
 77 
 78     hough2_roi = ralloc(sizeof(Hough2_region));
 79 
 80     if (hough2_roi!=NULL)
 81     {
 82       hough2_roi->min_x = min_x;
 83       hough2_roi->min_y = min_y;
 84       hough2_roi->max_x = max_x;
 85       hough2_roi->max_y = max_y;
 86       hough2_roi->binsize_x = binsize_x;
 87       hough2_roi->binsize_y = binsize_y;
 88 
 89       hough2->props = proplist_add(hough2->props, hough2_roi, HOUGH2_ROI_TYPE,  
 90                                                                         rfree);
 91     }
 92 
 93     return (hough2);
 94 }
 95 
 96 /* ---------- */
 97 
 98 /* This function extends the passed line and then clips it to the logical 
 99    coords of the hough space. If the extended line doesn't pass through the
100    hough space the passed line is unaltered. The function returns TRUE if 
101    line has been altered otherwise it returns FALSE */
102 
103 int hough2_extend_line(Line2 *line, Imrect *hough2)
104 {
105     double         min_x, min_y, max_x, max_y; 
106     Vec2           start, end;
107     double         bottom_edge_intersect; 
108     double         top_edge_intersect;     
109     double         left_edge_intersect;   
110     double         right_edge_intersect;
111     double         m, c;
112     Hough2_region *hough2_roi=NULL;
113     Line2         *new_line;
114 
115 /* Extract the logical coords of the hough space */
116 
117     hough2_roi = prop_get(hough2->props, HOUGH2_ROI_TYPE);
118 
119     if (hough2_roi==NULL)
120     {
121         format("Error: hough2_roi not defined.\n");
122         return FALSE;
123     }
124 
125     min_x = hough2_roi->min_x;
126     min_y = hough2_roi->min_y;
127     max_x = hough2_roi->max_x;
128     max_y = hough2_roi->max_y;
129  
130     if (vec2_x(line->v)==0.0) /* Special case for vertical lines */
131     {
132         if ((vec2_x(line->p1)>=min_x) &&
133             (vec2_x(line->p1)<=max_x))
134             {           
135                 vec2_x(start) = vec2_x(line->p1);
136                 vec2_y(start) = (float)min_y;
137                 vec2_x(end)   = vec2_x(line->p1);
138                 vec2_y(end)   = (float)max_y;
139                 new_line = line2_make(start, end, 0);
140                 *line = *new_line;
141                 line2_free(new_line);
142                 return TRUE;
143             }
144             else
145             {
146                 return FALSE;
147             }
148     }/*endif vertical*/
149 
150     if (vec2_y(line->v)==0.0) /* Special case for horizontal lines */
151     {
152         if ((vec2_y(line->p1)>=min_y) &&
153             (vec2_y(line->p1)<=max_y))
154             {           
155                 vec2_x(start) = (float)min_x;
156                 vec2_y(start) = vec2_y(line->p1);
157                 vec2_x(end)   = (float)max_x;
158                 vec2_y(end)   = vec2_y(line->p1);
159                 new_line = line2_make(start, end, 0);
160                 *line = *new_line;
161                 line2_free(new_line);
162                 return TRUE;
163             }
164             else
165             {
166                 return FALSE;
167             }
168     }/*endif horizontal*/
169 
170 /* The line is neither horizontal or vertical - determine where it intersects
171    the lines which define the hough space. */
172 
173     m = vec2_y(line->v)/vec2_x(line->v);
174     c = vec2_y(line->p1)-m*vec2_x(line->p1);
175 
176     left_edge_intersect   = m*min_x+c;
177     right_edge_intersect  = m*max_x+c;
178     bottom_edge_intersect = (min_y-c)/m;
179     top_edge_intersect    = (max_y-c)/m;
180 
181 /* Return if this line doesn't cross the hough space */
182 
183     if (((left_edge_intersect<min_y) && (right_edge_intersect<min_y)) ||
184         ((left_edge_intersect>max_y) && (right_edge_intersect>max_y)) ||
185         ((bottom_edge_intersect<min_x) && (top_edge_intersect<min_x)) ||
186         ((bottom_edge_intersect>max_x) && (top_edge_intersect>max_x))) return FALSE;
187   
188       
189 /* Initially assume that the left and right intersects are the line endpoints */
190 
191     vec2_x(start) = (float)min_x;
192     vec2_y(start) = (float)left_edge_intersect;
193     vec2_x(end)   = (float)max_x;
194     vec2_y(end)   = (float)right_edge_intersect;
195 
196 /* Only change these assumptions if necessary */
197 
198     if (left_edge_intersect<min_y)
199     {
200         vec2_x(start) = (float)bottom_edge_intersect;
201         vec2_y(start) = (float)min_y;
202     }
203 
204     if (left_edge_intersect>max_y)
205     {
206         vec2_x(start) = (float)top_edge_intersect;
207         vec2_y(start) = (float)max_y;
208     }
209 
210 
211     if (right_edge_intersect<min_y)
212     {
213         vec2_x(end) = (float)bottom_edge_intersect;
214         vec2_y(end) = (float)min_y;
215     }
216 
217     if (right_edge_intersect>max_y)
218     {
219         vec2_x(end) = (float)top_edge_intersect;
220         vec2_y(end) = (float)max_y;
221     }
222 
223 /* Now change the old line to the new line */
224 
225     new_line = line2_make(start, end, 0);
226     *line = *new_line;
227     line2_free(new_line);
228     return TRUE;
229 }
230 
231 /* ---------- */
232 
233 void hough2_plot_line(Line2 *line, Imrect *hough2, double weight)
234 {
235     int            x_diff, y_diff;
236     double         temp, min_x, min_y, max_x, max_y, binsize_x, binsize_y;
237     double            x0, y0, x1, y1;
238     Hough2_region *hough2_roi;
239 
240         /* (1) Transform the line coords from logical coords
241                to hough space array indices.
242            (2) Swap start and end points such that x0 is 
243                always less than x1.
244            (3) The problem is now limited to 4 cases:
245                    (i)   y' positive and greater than x'
246                    (ii)  y' positive and less than x'
247                    (iii) y' negative and abs less than x'
248                    (iv)  y' negative and abs greater than x'
249                Test for case and execute appropriately.  */
250 
251 /* Extract the logical coords of the hough space and the binsizes */
252 
253     hough2_roi = prop_get(hough2->props, HOUGH2_ROI_TYPE);
254 
255     if (hough2_roi==NULL)
256     {
257         format("Error: hough2 roi not set.\n");
258         return;
259     }
260 
261     min_x = hough2_roi->min_x;
262     min_y = hough2_roi->min_y;
263     max_x = hough2_roi->max_x;
264     max_y = hough2_roi->max_y;
265     binsize_x = hough2_roi->binsize_x;
266     binsize_y = hough2_roi->binsize_y;
267 
268 /* The array indices can now be determined */
269 
270     x0 = (vec2_x(line->p1)-min_x)/binsize_x;
271     y0 = (vec2_y(line->p1)-min_y)/binsize_y;
272     x1 = (vec2_x(line->p2)-min_x)/binsize_x;
273     y1 = (vec2_y(line->p2)-min_y)/binsize_y;
274 
275 /* Swap start and end points if necessary */
276 
277     if (x0>x1)
278     {
279         temp = x0;
280         x0 = x1;
281         x1 = temp;
282 
283         temp = y0;
284         y0 = y1;
285         y1 = y0;
286     }
287 
288 /* Determine which case applies */
289 
290     x_diff = (int)(x1-x0);
291     y_diff = (int)(y1-y0);
292 
293     if (y_diff>=0) /* +ve gradient */
294     {
295         if (y_diff>=x_diff)
296             y_based_plot(y0, y1, x0, (double)x_diff/(double)y_diff, 
297                                                             hough2, weight);
298         else
299             x_based_plot(x0, x1, y0, (double)y_diff/(double)x_diff, 
300                                                             hough2, weight);
301     }
302     else  /* -ve gradient */
303     {
304         if ((-1*y_diff)>=x_diff)
305             y_based_plot(y1, y0, x1, (double)x_diff/(double)y_diff, 
306                                                             hough2, weight);
307         else
308             x_based_plot(x0, x1, y0, (double)y_diff/(double)x_diff, 
309                                                             hough2, weight);
310     }
311 }
312 
313 /* ---------- */
314 
315 void y_based_plot(double y0, double y1, double x0, double grad, Imrect *hough2, double weight)
316 {
317     int       y;
318     double    x;
319     Vartype   vtype;
320     int       lx, ly, ux, uy;  /* max and min image array indices */
321     double    gl;
322 
323     vtype =   hough2->vtype;
324 
325     lx = hough2->region->lx;
326     ly = hough2->region->ly;
327     ux = hough2->region->ux;
328     uy = hough2->region->uy;
329 
330     x = x0 + grad *(y0+0.5-tina_rint(y0));
331 
332     if (y0<ly) y0 = ly;
333     if (y1>uy) y1 = uy;
334 
335     for(y=tina_rint(y0);y<y1;y++)
336     { 
337         if ( (tina_rint(x-1)>=lx) && (tina_rint(x+1)<ux) )
338         {
339             IM_PIX_GET(hough2,y, tina_rint(x-1), gl);
340             IM_PIX_SET(hough2,y, tina_rint(x-1), gl+weight);
341             IM_PIX_GET(hough2,y, tina_rint(x), gl);
342             IM_PIX_SET(hough2,y, tina_rint(x), gl+weight);
343             IM_PIX_GET(hough2,y, tina_rint(x+1), gl);
344             IM_PIX_SET(hough2,y, tina_rint(x+1), gl+weight);
345         }
346 
347         x+=grad;
348     }
349 }
350 
351 /* ---------- */
352 
353 void x_based_plot(double x0, double x1, double y0, double grad, Imrect *hough2, double weight)
354 {
355     int    x;
356     double y;
357     double gl;
358     Vartype   vtype;
359     int       lx, ly, ux, uy;  /* max and min image array indices */
360 
361     vtype =   hough2->vtype;
362 
363     lx = hough2->region->lx;
364     ly = hough2->region->ly;
365     ux = hough2->region->ux;
366     uy = hough2->region->uy;
367 
368     y = y0 + grad *(x0+0.5-tina_rint(x0));
369 
370     if (x0<lx) x0 = lx;
371     if (x1>ux) x1 = ux;
372 
373     for(x=tina_rint(x0);x<x1;x++)
374     {
375       if ( (tina_rint(y-1)>=ly) && (tina_rint(y+1)<uy) )
376       {
377           IM_PIX_GET(hough2,tina_rint(y-1),x,gl);
378           IM_PIX_SET(hough2,tina_rint(y-1),x,gl+weight);
379           IM_PIX_GET(hough2,tina_rint(y),x,gl);
380           IM_PIX_SET(hough2,tina_rint(y),x,gl+weight);
381           IM_PIX_GET(hough2,tina_rint(y+1),x,gl);
382           IM_PIX_SET(hough2,tina_rint(y+1),x,gl+weight);
383       }
384       y+=grad;    
385     }
386 }
387 
388 /* ---------- */
389 
390 double im_get_quadcovar(Imrect * image, double x, double y, float *px, float *py
391 , Mat2 *C)
392 {
393     Imregion       *region = image->region;
394     double           pixval[3][3];
395     double           a, b, c, d, e, f;
396     double          inter;
397     double           temp;
398     short           i, j, n, m;
399     double           xs, ys;
400     double det;
401 
402     i = tina_int(x - 1);
403     j = tina_int(y - 1);
404 
405     if (j < region->ly + 1 || j > region->uy - 1
406         || i < region->lx + 1 || i > region->ux - 1)
407         return (0);
408 
409     for (n = 0; n < 3; n++)
410     {
411         for (m = 0; m < 3; m++)
412         {
413             pixval[n][m] = im_get_pixf(image, j + n, i + m);
414         }
415     }
416 
417     a = pixval[1][1];
418     b = (pixval[0][2] - pixval[0][0]
419         + pixval[1][2] - pixval[1][0]
420         + pixval[2][2] - pixval[2][0]) / 6.0;
421     c = (pixval[2][0] - pixval[0][0]
422         + pixval[2][1] - pixval[0][1]
423         + pixval[2][2] - pixval[0][2]) / 6.0;
424     d = (pixval[0][0] - 2.0 * pixval[0][1] + pixval[0][2]
425         + 3.0 * pixval[1][0] - 6.0 * pixval[1][1] + 3.0 * pixval[1][2]
426         + pixval[2][0] - 2.0 * pixval[2][1] + pixval[2][2]) / 10.0;
427     e = (pixval[0][0] - pixval[2][0]
428         + pixval[2][2] - pixval[0][2]) / 4.0;
429     f = (pixval[0][0] + 3.0 * pixval[0][1] + pixval[0][2]
430         - 2.0 * pixval[1][0] - 6.0 * pixval[1][1] - 2.0 * pixval[1][2]
431         + pixval[2][0] + 3.0 * pixval[2][1] + pixval[2][2]) / 10.0;
432 
433     temp = 4.0 * d * f - e * e;
434     *px = (float)1.5 + (float) i + (e * c - 2 * f * b) / temp;
435     *py = (float)1.5 + (float) j + (e * b - 2 * d * c) / temp;
436 
437     xs = (e * c - 2 * f * b) / temp;
438     ys = (e * b - 2 * d * c) / temp;
439 
440     inter = a + b * xs + c * ys + d * xs * xs + e * xs * ys + f * ys * ys;
441 
442 
443     /* set the covariace matrix */
444     det = d*f-e*e/4.0;
445     mat2_xx(*C) = (float)(-f/det);
446     mat2_yy(*C) = (float)(-d/det);
447     mat2_xy(*C) = mat2_yx(*C) = (float)((e/2.0)/det);
448 
449     return (inter);
450 }
451 
452 List *hough2_locate_peaks(Imrect *hough2, double thres)
453 {
454     double         min_x, min_y, max_x, max_y; /* logical coords */
455     int            x0, x1, y0, y1;             /* array indices  */
456     int            i, j;
457     double         binsize_x, binsize_y;
458     Vec2           p;
459     float          x, y;
460     List          *peaks_list=NULL;
461     double         top_lef, top_cen, top_rig;
462     double         mid_lef, mid_cen, mid_rig;
463     double         bot_lef, bot_cen, bot_rig;
464     double        *peak_val;
465     Point2        *peak;
466     Hough2_region *hough2_roi;
467     Vartype        vtype;
468     Mat2 *C;
469 
470 
471     if (hough2==NULL)
472     {
473         format("Error: hough2 not defined.\n");
474         return NULL;
475     }
476 
477 /* Extract the physical coords of the hough space */
478 
479     x0 = hough2->region->lx;
480     y0 = hough2->region->ly;
481     x1 = hough2->region->ux;
482     y1 = hough2->region->uy;
483 
484 /* Extract the logical coords of the hough space and the binsizes */
485 
486     if ((hough2_roi = prop_get(hough2->props, HOUGH2_ROI_TYPE)) == NULL)
487     {
488         format("Error: hough2_roi no defined.\n");
489         return NULL;
490     }
491 
492     min_x = hough2_roi->min_x;
493     min_y = hough2_roi->min_y;
494     max_x = hough2_roi->max_x;
495     max_y = hough2_roi->max_y;
496     binsize_x = hough2_roi->binsize_x;
497     binsize_y = hough2_roi->binsize_y;
498 
499 /* If the hough2 is less than 3x3 bins cannot find peaks */
500 
501     if ((x1-x0 < 3) || (y1-y0 < 3))
502     {
503         format("Error: hough2 too small to find peaks.\n");
504         return NULL;
505     }
506 
507 /* Access the image data */ 
508 
509     vtype =   hough2->vtype;
510 
511 /* Pass 3x3 window over the hough2. Add peaks to the peaks_list. */
512 
513     for(j=(y0+1);j<(y1-1);j++)
514         for(i=(x0+1);i<(x1-1);i++)
515         {
516             IM_PIX_GET(hough2, j  , i  , mid_cen);
517             IM_PIX_GET(hough2, j-1, i-1, top_lef);
518             IM_PIX_GET(hough2, j-1, i  , top_cen);
519             IM_PIX_GET(hough2, j-1, i+1, top_rig);
520             IM_PIX_GET(hough2, j  , i-1, mid_lef);
521             IM_PIX_GET(hough2, j  , i+1, mid_rig);
522             IM_PIX_GET(hough2, j+1, i-1, bot_lef);
523             IM_PIX_GET(hough2, j+1, i  , bot_cen);
524             IM_PIX_GET(hough2, j+1, i+1, bot_rig);
525 
526             if ( (mid_cen>=thres ) &&
527                  (mid_cen>=top_lef) &&
528                  (mid_cen>=top_cen) &&
529                  (mid_cen>top_rig) &&
530                  (mid_cen>=mid_lef) &&
531                  (mid_cen>mid_rig) &&
532                  (mid_cen>=bot_lef) &&
533                  (mid_cen>bot_cen) &&
534                  (mid_cen>bot_rig))
535                 {
536                     /* Interpolate surrounding pixels to determine peak pos
537                        and store the covariance matrix on the props list */
538 
539                     C = mat2_alloc();
540                     peak_val = (double *)ralloc(sizeof(double));
541                     *peak_val = im_get_quadcovar(hough2, i, j, &x, &y, C);
542                     
543                     vec2_x(p) = (float)(((double)x)*binsize_x+min_x);
544                     vec2_y(p) = (float)(((double)y)*binsize_y+min_y);
545                     peak = point2_make(p, HOUGH_PEAK_TYPE);
546                     peaks_list = ref_addtostart(peaks_list, peak, POINT2);
547                     peak->props = proplist_add(peak->props, peak_val, HOUGH2_PEAK,rfree);
548                     peak->props = proplist_add(peak->props, C, COVAR2, mat2_free);
549                 }
550         }/*endfor(i)*/
551 
552     return (peaks_list);                           
553 }
554 
555 /* ---------- */
556 
557 Vec2 *hough2_locate_max_peak(Imrect *hough2, double thres)
558 {
559     List *peak_list;
560     Point2 *p;
561     double peak_func();
562     Vec2 *vec2_copy();
563 
564     peak_list = hough2_locate_peaks(hough2, thres);
565 
566     if (peak_list==NULL) return NULL;
567 
568     peak_list = sort_list(peak_list, peak_func, NULL);
569     p = peak_list->to;
570 
571     return vec2_copy(&p->p);
572 
573     list_rm(peak_list, rfree);  /* ?? not reached */
574 }
575 
576 /* ---------- */
577 
578 double peak_func(Point2 *p)
579 {
580     double peak_val;
581     peak_val = *(double *)prop_get(p->props,HOUGH2_PEAK);
582     return (-peak_val);
583 }
584 
585 /* ---------- */
586 
587 void hough2_fill(Imrect *hough2, double val)
588 {
589     int lx, ly, ux, uy;             /* array indices  */
590         int i, j;
591 
592     if (hough2==NULL)
593     {
594         format("Error: hough2 not defined.\n");
595                 return;
596     }
597 
598 /* Extract the physical coords of the hough space */
599 
600     lx = hough2->region->lx;
601     ly = hough2->region->ly;
602     ux = hough2->region->ux;
603     uy = hough2->region->uy;
604 
605         for(i=lx;i<ux;i++)
606         {
607                 for(j=ly;j<uy;j++)
608                 {
609                         /*IM_PIX_SET(hough2, j, i, val);*/
610                         ((double**)hough2->data)[j][i] = val;
611                 }
612         }
613 }
614 
615 /* ---------- */
616 
617 void hough2_plot_ellipse(Imrect *hough2, Vec2 *c, Mat2 *C, 
618                          double weight) 
619 {
620     Hough2_region *hough2_roi;    
621     double minx, miny, maxx, maxy, binsizex, binsizey;
622     double ellipse_bound;
623     int ellipse_minx, ellipse_maxx, ellipse_miny, ellipse_maxy;
624     int x, y;
625     double px, py, dx, dy;
626     double chi2, f, g;
627     int lx, ly, ux, uy;
628     Mat2 invC;
629     double det;
630 
631     if (hough2==NULL)
632     {
633         format("Error: hough2 not defined.\n");
634         return;
635     }
636 
637 /* Extract the physical coords of the hough space */
638 
639     lx = hough2->region->lx;
640     ly = hough2->region->ly;
641     ux = hough2->region->ux;
642     uy = hough2->region->uy;
643 
644 /* Extract the logical coords of the hough space and the binsizes */
645 
646     if ((hough2_roi = prop_get(hough2->props, HOUGH2_ROI_TYPE)) == NULL)
647     {
648         format("Error: hough2_roi not defined.\n");
649         return;
650     }
651 
652     minx = hough2_roi->min_x;
653     miny = hough2_roi->min_y;
654     maxx = hough2_roi->max_x;
655     maxy = hough2_roi->max_y;
656     binsizex = hough2_roi->binsize_x;
657     binsizey = hough2_roi->binsize_y;
658 
659     if (mat2_xx(*C) > mat2_yy(*C))
660     {
661         ellipse_bound = 3.0*sqrt(mat2_xx(*C));
662     }
663     else
664     {
665         ellipse_bound = 3.0*sqrt(mat2_yy(*C));
666     }
667 
668     ellipse_minx = (int)((vec2_x(*c) - ellipse_bound - minx)/binsizex);
669     ellipse_maxx = (int)((vec2_x(*c) + ellipse_bound - minx)/binsizex);
670     ellipse_miny = (int)((vec2_y(*c) - ellipse_bound - miny)/binsizey);
671     ellipse_maxy = (int)((vec2_y(*c) + ellipse_bound - miny)/binsizey);
672 
673     /* Determine the inverse of the covariance matrix */
674 
675     det = mat2_xx(*C)*mat2_yy(*C) - mat2_xy(*C)*mat2_yx(*C);
676 
677     mat2_xx(invC) = (float)(mat2_yy(*C)/det);
678     mat2_xy(invC) = (float)(-mat2_xy(*C)/det);
679     mat2_yx(invC) = (float)(-mat2_yx(*C)/det);
680     mat2_yy(invC) = (float)(mat2_xx(*C)/det);
681 
682     for(x = ellipse_minx; x <= ellipse_maxx; x++)
683     {
684         for(y = ellipse_miny; y <= ellipse_maxy; y++)
685         {
686             if ( (x>=lx) && (x<ux) && (y>=ly) && (y<uy) )
687             {
688                 /* Get the position of this pixel centre */
689                 px = ((double)x+0.5)*binsizex + minx;
690                 py = ((double)y+0.5)*binsizey + miny;
691 
692                 dx = px-vec2_x(*c);
693                 dy = py-vec2_y(*c);
694 
695                 chi2 = dx*dx*mat2_xx(invC) + dx*dy*mat2_xy(invC) + 
696                        dx*dy*mat2_yx(invC) + dy*dy*mat2_yy(invC);
697 
698                 f = weight*(9.0-chi2);
699   
700                 if (f>0.0)
701                 {
702                     g = ((double**)hough2->data)[y][x];
703                     ((double**)hough2->data)[y][x] = (double)f+g;
704                 }
705             }   
706         }
707     }
708 }
709 
710 /* ---------- */
711 
712 void sqrt_hough2(Imrect *hough2)
713 {
714     int     i,j, min_x, min_y, max_x, max_y;
715     double  hough2_data;
716 
717     min_x = hough2->region->lx;
718     min_y = hough2->region->ly;
719     max_x = hough2->region->ux;
720     max_y = hough2->region->uy;
721 
722 /* Each histogram entry is replaced by its square root */
723 
724     for(i=min_y;i<max_y;i++)
725         for(j=min_x;j<max_x;j++)
726         {
727             IM_PIX_GET(hough2,i,j,hough2_data);
728             IM_PIX_SET(hough2,i,j,(double)sqrt((double)hough2_data));
729         }
730 }
731 
732 /* ---------- */
733 

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