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

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

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