blob: 2f89f9c800f6efdf42424f978f31e62c6bd182eb [file] [log] [blame]
/****
Copyright (C) 1996 McGill University.
Copyright (C) 1996 McCAT System Group.
Copyright (C) 1996 ACAPS Benchmark Administrator
benadmin@acaps.cs.mcgill.ca
This program is free software; you can redistribute it and/or modify
it provided this copyright notice is maintained.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
****/
/* ------------------------------------------------------------------------
* L_canny.c Canny Edge Detection
*
* V
* University of Illinois at Urbana-Champaign
*
* Program written by A. Tabbone (CRIN/CNRS- Lorraine, FR)
* Modified for V by C. P. Hess (Univ. Illinois at Urbana, USA)
* Original copyright appears at EOF.
*
* $Id$
* ------------------------------------------------------------------------ */
#include <math.h>
#include <ctype.h>
#include "V.h"
int GaussianMask();
int DGaussianMask();
int dfilter();
float *dnon_max();
int L_canny(sigma, image, nc, nr, edge, err)
float sigma;
float *image;
int nc;
int nr;
float **edge;
char *err;
{
float *Ix, *Iy, *gr, *gmask, *dgmask;
int tgauss, tdgauss, i;
if(GaussianMask (sigma, &tgauss, &gmask, err)) return(1);
if(DGaussianMask(sigma, &tdgauss, &dgmask, err)) return(1);
if(dfilter(image, gmask, dgmask, tgauss, tdgauss,
nc, nr, &Ix, err)) return(1);
if(dfilter(image, dgmask, gmask, tdgauss, tgauss,
nc, nr, &Iy, err)) return(1);
gr = (float *) calloc(nr*nc, FWS);
if (!gr) {
sprintf(err, "Out of memory");
return(1);
}
for (i=0; i<nr*nc; i++) gr[i] = MOD(Iy[i], Ix[i]);
*edge = dnon_max(gr, Ix, Iy, nc, nr);
return(0);
}
int GaussianMask(sigma, tg, mask, err)
float sigma;
int *tg;
float **mask;
char *err;
{
short i;
float *coeff_gauss;
float *p;
*tg = (short) (GAUSS_MASK*sigma);
if (*tg % 2 == 0) *tg += 1;
coeff_gauss = p = (float *)calloc(*tg, FWS);
if (!coeff_gauss) {
sprintf(err,"Out of memory");
return(1);
}
for (i=-(*tg/2); i<=*tg/2; i++) {
if ((i+(*tg/2)) == 0)
*coeff_gauss++ =
(float)((erf((double)(i+0.5)/(sqrt((double)2.0)*sigma))+1.0)/2.0);
else if ((i+(*tg/2)) == *tg-1)
*coeff_gauss++ =
(float)((-erf((double)(i-0.5)/(sqrt((double)2.0)*sigma))+1.0)/2.0);
else *coeff_gauss++ =
(float)0.5*(erf((double)(i+0.5)/(sqrt((double) 2.0)*sigma))-
erf((double)(i-0.5)/(sqrt((double) 2.0)*sigma)));
}
*mask = p;
return(0);
}
int DGaussianMask(sigma, tg, mask, err)
float sigma;
int *tg;
float **mask;
char *err;
{
short i;
float *coeff_dgauss;
float *p;
*tg = (short) (DGAUSS_MASK*sigma);
if (*tg % 2 == 0) *tg += 1;
coeff_dgauss = p = (float *)calloc(*tg, FWS);
if (!coeff_dgauss) {
sprintf(err,"Out of memory");
return(1);
}
for (i=-(*tg/2); i<=*tg/2; i++) {
if ((i+(*tg/2)) == 0)
*coeff_dgauss++ =
(float)(1.0/(sigma*sqrt((double)(2.0*M_PI))))*
exp(-(double)((i+0.5)*(i+0.5))/(double)(2.0*sigma*sigma));
else if ((i+(*tg/2)) == *tg-1)
*coeff_dgauss++ =
(float)-(1.0/(sigma*sqrt(2.0*M_PI)))*
exp(-(double)((i-0.5)*(i-0.5))/(double)(2.0*sigma*sigma));
else *coeff_dgauss++ =
(float)((exp(-(double)((i+0.5)*(i+0.5))/(double)(2.0*sigma*sigma))-
exp(-(double)((i-0.5)*(i-0.5))/(double)(2.0*sigma*sigma)))/
(sigma*sqrt((double) 2.0*M_PI)));
}
coeff_dgauss = p;
*mask = p;
return(0);
}
int dfilter(image, g, f, tm_g, tm_f, nc, nr, filt, err)
float *image, *g, *f;
int nc, nr, tm_g, tm_f;
float **filt;
char *err;
{
register short i, j, l, high = nr, larg = nc;
float nv, *d, *d1, *d2, *temp;
d = d1 = (float *) calloc(nc*nr, FWS);
if (!d) {
sprintf(err,"Out of memory");
return(1);
}
for (i=0; i<high; i++)
for (j=0; j<larg; j++) {
for (l=-(tm_g/2); l<=(tm_g/2); l++) {
if ((j+l) < 0) nv = (float) image[i*larg];
else if ((j+l) >= larg) nv = (float) image[((i+1)*larg)-1];
else nv = (float) image[(i*larg)+j+l];
*d = (nv * g[(tm_g/2)-l]) + *d;
}
d++;
}
d = d1;
d2 = (float *) calloc(nc*nr, FWS);
if (!d2) {
sprintf(err,"Out of memory");
return(1);
}
for (j=0; j<larg ; j++)
for (i=0; i<high; i++) {
for (l=-(tm_f/2); l<=(tm_f/2); l++) {
if ((i+l) < 0) nv = d[j];
else if ((i+l) >= high) nv = d[((high-1)*nc)+j];
else nv = d[((i+l)*larg)+j];
temp = (d2+(i*larg)+j);
*temp = (nv * f[(tm_f/2)-l]) + *temp;
}
}
*filt = d2;
return(0);
}
float *dnon_max(gr, Ix, Iy, nc, nr)
float *gr;
float *Ix;
float *Iy;
int nc;
int nr;
{
int i, j;
float *maxima;
float R, ampl1, ampl2;
maxima = (float *) calloc(nc*nr, FWS);
for (i=1;i<nr-2;i++)
for (j=1;j<nc-2;j++) {
if (ZERO(I(Ix,i,j,nc))) {
ampl1 = I(gr,i,j-1,nc); ampl2 = I(gr,i,j+1,nc);
} else {
R = I(Iy,i,j,nc) / I(Ix,i,j,nc);
if ((R >=0.) && (R < .4)) {
ampl1 = I(gr,i+1,j+1,nc) * R + (1.-R)* I(gr,i+1,j,nc);
ampl2 = I(gr,i-1,j-1,nc) * R + (1.-R)* I(gr,i-1,j,nc);
}
/*
else if ( R >= 1.) {
ampl1 = (I(gr,i+1,j+1,nc) + (R - 1.) * I(gr,i,j+1,nc)) / R;
ampl2 = (I(gr,i-1,j-1,nc) + (R - 1.) * I(gr,i,j-1,nc)) / R;
} */
else if ((R <= 0.) && (R > -.4)) {
ampl1 = - I(gr,i+1,j-1,nc) * R + ( 1. + R)* I(gr,i+1,j,nc);
ampl2 = - I(gr,i-1,j+1,nc) * R + ( 1. + R)* I(gr,i-1,j,nc);
}
/*
else if (R <= -1.) {
ampl1 = (- I(gr,i+1,j-1,nc) + (R + 1.) * I(gr,i,j-1,nc)) / R;
ampl2 = (- I(gr,i-1,j+1,nc) + (R + 1.) * I(gr,i,j+1,nc)) / R;
} */
else
{
I(maxima,i,j,nc) = 0;
continue;
}
}
if (GEPS(I(gr,i,j,nc), ampl1) && GEPS(I(gr,i,j,nc),ampl2))
I(maxima,i,j,nc) = NGMAX;
}
return(maxima);
}
/* ORIGINAL COPYRIGHT *****************************************************
* Copyright (C) 1994 CRIN *
* *
* This software was written by A. Tabbone at : *
* C.R.I.N/CNRS-INRIA lorraine *
* Boite Postale 239 *
* 54506 Vandoeuvre-les-Nancy Cedex *
* France *
* It may be distributed or copied, in whole or in part, within the *
* following restrictions: *
* (1) It may not be sold at a profit. *
* (2) This credit and notice must remain intact. *
* This software may be distributed with other software by a commercial *
* vendor, provided that it is included at no additional charge. *
* *
* Please report bugs to tabbone@loria.crin.fr *
***************************************************************************/