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

Linux Cross Reference
Tina5/tina-libs/tina/geometry/geomImg_canny.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/geometry/geomImg_canny.c,v $
 23  * Date    :  $Date: 2008/12/16 11:34:54 $
 24  * Version :  $Revision: 1.8 $
 25  * CVS Id  :  $Id: geomImg_canny.c,v 1.8 2008/12/16 11:34:54 neil Exp $
 26  *
 27  * Author  : Legacy TINA
 28  */
 29  /**
 30  * @file
 31  * @brief TINA implementation of the Canny edge detector.
 32  *
 33  * These functions implement the Canny edge detector.  The algorithm requires two thresholds:
 34  * the edge strength must be greater than the upper threshold in order to initiate detection,
 35  * but then the edge is followed until its strength falls below the lower threshold.  This 
 36  * hysteresis takes account of the fact that, if a neighbouring pixel has been confirmed as 
 37  * part of an edge, this offers additional evidence that the current pixel is part of an edge.
 38  * Two versions of the algorithm are offered here: the function canny is a standard implementation,
 39  * and the function iso_canny prioritises edges whose grey-level is consistent with a user provided 
 40  * grey level, allowing the user to pick the edge type detected (e.g. edge between grey matter and 
 41  * white matter in brain MRI).
 42  */
 43 
 44 #include "geomImg_canny.h"
 45 
 46 #if HAVE_CONFIG_H
 47   #include <config.h>
 48 #endif
 49 
 50 #include <math.h>
 51 #include <float.h>
 52 #include <tina/sys/sysDef.h>
 53 #include <tina/sys/sysPro.h>
 54 #include <tina/math/mathDef.h>
 55 #include <tina/math/mathPro.h>
 56 #include <tina/image/imgDef.h>
 57 #include <tina/image/imgPro.h>
 58 #include <tina/geometry/geom_EdgeDef.h>
 59 #include <tina/geometry/geomImg_nonmax.h>
 60 #include <tina/geometry/geomImg_thresh.h>
 61 #include <tina/geometry/geom_EdgePro.h>
 62 
 63 /**
 64  * @brief TINA implementation of the standard Canny edge detector.
 65  * @param im Image in which you wish to detect edges.
 66  * @param sigma Standard deviation of a Gaussian smoothing kernel applied to the image (0 = no amoothing).
 67  * @param precision The precision at which the Gaussian smoothing kernel is truncated.
 68  * @param lowthres Lower Canny threshold (in grey levels squared).
 69  * @param upthres Upper Canny threshold (in grey levels squared).
 70  * @param lengththres Minimum allowed length of edges.
 71  * @return edrect Image of the detected edges.
 72  *
 73  * Implementation of the standard Canny edge detector.  The image provided is smoothed using a Gaussian
 74  * kernel defined by the standard deviation "sigma" (0 = no smoothing) and truncated at the probability 
 75  * "precision".  An image of the sum-squared gradients in the x and y directions is then calculated, and 
 76  * edges stronger than the lower threshold are detected using non-maximal suppression. Finally, edges that
 77  * do not feature a voxel stronger than the upper threshold, or are shorter than the length threshold, are
 78  * removed. The edge image is then returned to the calling function.
 79  */
 80 Imrect *canny(Imrect * im, double sigma, double precision, double lowthres, double upthres, int lengththres)
 81 {
 82     Imrect *gim;
 83     Imrect *edrect;
 84     Imrect *gradx;
 85     Imrect *grady;
 86     Imrect *gradsq;
 87     unsigned int label;         /* for blocked allocation */
 88 
 89     if (im == NULL)
 90         return (NULL);
 91 
 92     label = ralloc_end_blocked();
 93 
 94     if (sigma == 0.0)
 95         gim = im_copy(im);
 96     else
 97         gim = imf_gauss(im, sigma, precision);
 98     gradx = imf_grad_h(gim);
 99     grady = imf_grad_v(gim);
100     gradsq = imf_sumsq(gradx, grady);
101     im_free(gim);
102 
103     if (label)                  /* allocation was blocked */
104         (void) ralloc_start_blocked(label);     /* re-start blocking */
105 
106     edrect = nonmaxsup(gradx, grady, gradsq, lowthres);
107     im_free(gradx);
108     im_free(grady);
109     im_free(gradsq);
110     er_find_edge_strings(edrect);
111     er_rm_edges(edrect, EDGE_GET_CONN_MASK, EDGE_NOLINK);
112     er_edge_strings_thres(edrect, lengththres, upthres);
113     er_set_row_index(edrect);
114     return (edrect);
115 }
116 
117 
118 static Imrect *get_ssqim(Imrect *im, double midpoint)
119 {
120   Imrect *im1 = NULL, *im2 = NULL, *im3 = NULL, *im4 = NULL, *im5 = NULL;
121   double sigma_x = 0.0, sigma_y = 0.0, noise = 0.0;
122   Imregion *roi = NULL;
123   double factor = 0.0;
124   
125   im1 = im_add(midpoint, im);
126   im2 = im_sqr(im1);
127   im3 = im_sqrt(im2);
128 
129   roi = im3->region;
130   sigma_x = imf_diffx_noise(im3, NULL);
131   sigma_y = imf_diffy_noise(im3, NULL);
132     
133   noise = (sigma_x + sigma_y)/2;
134   factor = -1/(noise*10);
135   
136   im4 = im_times(factor, im3);
137   im5 = im_exp(im4);
138   
139   im_free(im1);
140   im_free(im2);
141   im_free(im3);
142   im_free(im4);
143 
144   return(im5);
145 
146 }
147 
148 Imrect *iso_canny(Imrect * im, double sigma, double precision, double lowthres, double upthres, 
149                   int lengththres, double midpoint)
150 {
151     Imrect *gim;
152     Imrect *edrect;
153     Imrect *gradx;
154     Imrect *grady;
155     Imrect *ssqim;
156     unsigned int label;         /* for blocked allocation */
157 
158     if (im == NULL)
159         return (NULL);
160 
161     label = ralloc_end_blocked();
162 
163 
164 /* without gaussian smoothing to limit the gradient we can loose high gradient iso contours NAT 12/10/08 
165     tsim = im_tsmooth(im);
166 */
167 
168     if (sigma == 0.0)
169         gim = im_tsmooth(im);
170     else
171         gim = imf_gauss(im, sigma, precision);
172     gradx = imf_grad_h(gim);
173     grady = imf_grad_v(gim);
174     ssqim = get_ssqim(gim, midpoint);
175 
176     if (label)                  /* allocation was blocked */
177         (void) ralloc_start_blocked(label);     /* re-start blocking */
178 
179     edrect = nonmaxsup(gradx, grady, ssqim, lowthres);
180     im_free(gim);
181     im_free(gradx);
182     im_free(grady);
183     im_free(ssqim);
184     er_find_edge_strings(edrect);
185     er_rm_edges(edrect, EDGE_GET_CONN_MASK, EDGE_NOLINK);
186     er_edge_strings_thres(edrect, lengththres, upthres);
187     er_set_row_index(edrect);
188     return (edrect);
189 }
190 
191 void gradient_im(Imrect *im, Imrect **grad_im_x, Imrect **grad_im_y)
192 {
193    float             *row1, *row2, *row3, *row4, *row5, *swap; 
194    int               lx,ly,ux,uy;
195    int               i,j;
196    double            pixvala, pixvalb, pixvalc, mag;
197    float             theta;
198    
199    if(im==NULL) return;
200 
201    lx=im->region->lx;
202    ly=im->region->ly;
203    ux=im->region->ux;
204    uy=im->region->uy;
205    
206    row1 = fvector_alloc(lx, ux);
207    row2 = fvector_alloc(lx, ux);
208    row3 = fvector_alloc(lx, ux);
209    row4 = fvector_alloc(lx, ux);
210    row5 = fvector_alloc(lx, ux);
211 
212    *grad_im_x=im_alloc(im->height,im->width,im->region,float_v);
213    *grad_im_y=im_alloc(im->height,im->width,im->region,float_v);
214    
215     im_get_rowf(row1, im, ly, lx, ux);
216     im_get_rowf(row2, im, ly+1, lx, ux);
217     im_get_rowf(row3, im, ly+2, lx, ux);
218     im_get_rowf(row4, im, ly+3, lx, ux);
219     for(i=ly+2;i<uy-2;++i)
220     {
221         im_get_rowf(row5, im, i+2, lx, ux);
222         for(j=lx+2;j<ux-2;++j)
223         {
224             pixvala = -0.25*(2.0*row3[j]-row3[j-2]-row3[j+2]);
225             pixvalc = -0.25*(2.0*row3[j]-row1[j]-row5[j]);
226 /* cross derivative is sqrt(2) bigger than expected.. factor of 2 -> sqrt(2) here .. NAT */
227             pixvalb = 0.3535*(row2[j-1]-row2[j+1]-row4[j-1]+row4[j+1]);
228 
229             theta = atan2(pixvalb,pixvala-pixvalc);
230             if (theta < 0.0) theta += 2.0*PI;
231             theta /= 2.0;
232 /* etimate signal to noise on orientation and scale gradient to be consistent with step edge
233    estimates for later use NAT */
234             mag = 3.0*sqrt(pixvalb*pixvalb+(pixvala-pixvalc)*(pixvala-pixvalc));
235             IM_PIX_SET(*grad_im_x,i,j,mag*sin(theta));
236             IM_PIX_SET(*grad_im_y,i,j,-mag*cos(theta));
237         }
238         swap = row1;
239         row1 = row2;
240         row2 = row3;
241         row3 = row4;
242         row4 = row5;
243         row5 = swap;
244     }
245     fvector_free(row1, lx);
246     fvector_free(row2, lx);
247     fvector_free(row3, lx);
248     fvector_free(row4, lx);
249     fvector_free(row5, lx);
250 }
251 
252 Imrect *var_canny(Imrect * im, double sigma, double precision, double lowthres, double upthres, int lengththres)
253 {
254     Imrect *im1, *im2, *im3, *im4, *im5, *im6, *im7;
255     Imrect *edrect;
256     Imrect *gradx, *grady;
257     unsigned int label;         /* for blocked allocation */
258 
259     if (im == NULL)
260         return (NULL);
261 
262     label = ralloc_end_blocked();
263 
264     im1 = imf_sqr(im);
265     im2 = imf_gauss(im1, sigma, precision);
266     im3 = imf_gauss(im, sigma, precision);
267     im4 = imf_sqr(im3);
268     im5 = im_diff(im2,im4);
269     im6 = imf_thresh(0.0, im5);
270     im_free(im1);
271     im_free(im2);
272     im_free(im4);
273     im_free(im5);
274     im7 = imf_sqrt(im6);
275     gradient_im(im7, &gradx, &grady);
276     im_free(im3);
277 
278     if (label)                  /* allocation was blocked */
279         (void) ralloc_start_blocked(label);     /* re-start blocking */
280 
281     edrect = nonmaxsup(gradx, grady, im6, lowthres);
282 
283     im_free(im7);
284     im_free(gradx);
285     im_free(grady);
286     im_free(im6);
287     er_find_edge_strings(edrect);
288     er_rm_edges(edrect, EDGE_GET_CONN_MASK, EDGE_NOLINK);
289     er_edge_strings_thres(edrect, lengththres, upthres);
290     er_set_row_index(edrect);
291 
292     return (edrect);
293 }
294 
295 Imrect *rank_canny(Imrect * im, double sigma, double precision, double lowthres, double upthres, int lengththres)
296 {
297     Imrect *im1,*im2;
298     Imrect *edrect;
299     Imrect *gradx;
300     Imrect *grady;
301     unsigned int label;         /* for blocked allocation */
302 
303     if (im == NULL)
304         return (NULL);
305 
306     label = ralloc_end_blocked();
307 
308     im1 = imf_gauss(im, sigma, precision);
309     im2 = im_rank(im1, 5, lowthres*10);
310     
311     gradient_im(im2, &gradx, &grady);
312 
313     if (label)                  /* allocation was blocked */
314         (void) ralloc_start_blocked(label);     /* re-start blocking */
315 
316     edrect = nonmaxsup(gradx, grady, im2, 2.0*lowthres);
317     im_free(gradx);
318     im_free(grady);
319     im_free(im1);
320     im_free(im2);
321     er_find_edge_strings(edrect);
322     er_rm_edges(edrect, EDGE_GET_CONN_MASK, EDGE_NOLINK);
323     er_edge_strings_thres(edrect, lengththres, 2.0*upthres);
324     er_set_row_index(edrect);
325     return (edrect);
326 }
327 
328 Imrect *im_canny(Imrect * im, double sigma, double precision, double lowthres, double upthres, int lengththres)
329 {
330     Imrect *im1;
331     Imrect *edrect;
332     Imrect *gradx;
333     Imrect *grady;
334     unsigned int label;         /* for blocked allocation */
335 
336     if (im == NULL)
337         return (NULL);
338 
339     label = ralloc_end_blocked();
340 
341     im1 = imf_gauss(im, sigma, precision);
342 
343     gradient_im(im1, &gradx, &grady);
344 
345     if (label)                  /* allocation was blocked */
346         (void) ralloc_start_blocked(label);     /* re-start blocking */
347 
348     edrect = nonmaxsup(gradx, grady, im1, 2.0*lowthres);
349     im_free(gradx);
350     im_free(grady);
351     im_free(im1);
352     er_find_edge_strings(edrect);
353     er_rm_edges(edrect, EDGE_GET_CONN_MASK, EDGE_NOLINK);
354     er_edge_strings_thres(edrect, lengththres, 2.0*upthres);
355     er_set_row_index(edrect);
356     return (edrect);
357 }
358 
359 

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