发信人: zjliu (秋天的萝卜), 信区: Graphics
标 题: Canny algorithm for edge detection in c
发信站: 哈工大紫丁香 (Sun Nov 16 20:31:19 2003), 站内信件
The code implementing the Canny edge detection algorithm is shown below. It
carries out the 6 steps mentioned above as well as some additional steps tha
t are explained in my other tutorials. These can be found on my home page.
/*
FILE: canny.c -
AUTH: Bill Green
DESC: Canny edge detector
DATE: 08/17/02
Program Steps:
1) 5x5 Gaussian mask for noise removal
2) Use Sobel masks to approximate gradient
3) Find orientation of gradient
4) Implement nonmaximum suppression to assign edges
5) Hysteresis thresholding (2 thresholds)
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <alloc.h>
/*-------STRUCTURES---------*/
typedef struct {int rows; int cols; unsigned char* data;} sImage;
/*-------PROTOTYPES---------*/
long getImageInfo(FILE*, long, int);
void copyImageInfo(FILE* inputFile, FILE* outputFile);
void copyColorTable(FILE* inputFile, FILE* outputFile, int nColors);
int main(int argc, char* argv[])
{
FILE *bmpInput, *bmpOutput, *dataOutput;
sImage originalImage, gaussImage, edgeImage;
unsigned int r, c;
int I, J;
int sumX, sumY, SUM;
int leftPixel, rightPixel;
int P1, P2, P3, P4, P5, P6, P7, P8;
int nColors;
unsigned long vectorSize;
unsigned long fileSize;
int MASK[5][5];
int GX[3][3];
int GY[3][3];
unsigned char *pChar, someChar;
unsigned int row, col;
float ORIENT;
int edgeDirection;
int highThreshold, lowThreshold;
someChar = '0'; pChar = &someChar;
/* 5x5 Gaussian mask. Ref: http://www.cee.hw.ac.uk/hipr/html/gsmooth.htm
l */
MASK[0][0] = 2; MASK[0][1] = 4; MASK[0][2] = 5; MASK[0][3] = 4; MASK[0
][4] = 2;
MASK[1][0] = 4; MASK[1][1] = 9; MASK[1][2] = 12; MASK[1][3] = 9; MASK[1
][4] = 4;
MASK[2][0] = 5; MASK[2][1] = 12; MASK[2][2] = 15; MASK[2][3] = 12; MASK[2
][4] = 5;
MASK[3][0] = 4; MASK[3][1] = 9; MASK[3][2] = 12; MASK[3][3] = 9; MASK[3
][4] = 4;
MASK[4][0] = 2; MASK[4][1] = 4; MASK[4][2] = 5; MASK[4][3] = 4; MASK[4
][4] = 2;
/* 3x3 GX Sobel mask. Ref: www.cee.hw.ac.uk/hipr/html/sobel.html */
GX[0][0] = -1; GX[0][1] = 0; GX[0][2] = 1;
GX[1][0] = -2; GX[1][1] = 0; GX[1][2] = 2;
GX[2][0] = -1; GX[2][1] = 0; GX[2][2] = 1;
/* 3x3 GY Sobel mask. Ref: www.cee.hw.ac.uk/hipr/html/sobel.html */
GY[0][0] = 1; GY[0][1] = 2; GY[0][2] = 1;
GY[1][0] = 0; GY[1][1] = 0; GY[1][2] = 0;
GY[2][0] = -1; GY[2][1] = -2; GY[2][2] = -1;
if(argc < 2) {
printf("Usage: %s bmpInput.bmp\n", argv[0]);
exit(0);
};
printf("Reading filename %s\n", argv[1]);
/* open files for reading and writing to */
bmpInput = fopen(argv[1], "rb");
bmpOutput = fopen("canny.bmp", "wb");
dataOutput = fopen("canny.txt", "w");
/* start pointer at beginning of file */
fseek(bmpInput, 0L, SEEK_END);
/* retrieve and print filesize and number of cols and rows */
fileSize = getImageInfo(bmpInput, 2, 4);
originalImage.cols = (int)getImageInfo(bmpInput, 18, 4);
originalImage.rows = (int)getImageInfo(bmpInput, 22, 4);
gaussImage.rows = edgeImage.rows = originalImage.rows;
gaussImage.cols = edgeImage.cols = originalImage.cols;
printf("Width: %d\n", originalImage.cols);
printf("Height: %d\n", originalImage.rows);
printf("File size: %lu\n", fileSize);
/* retrieve and print Number of colors */
nColors = (int)getImageInfo(bmpInput, 46, 4);
printf("nColors: %d\n", nColors);
/* allocate memory for originalImage.data */
vectorSize = fileSize - (14+40+4*nColors);
originalImage.data = farmalloc(vectorSize*sizeof(unsigned char));
if(originalImage.data == NULL) {
printf("Failed to malloc originalImage.data\n");
exit(0);
}
printf("%lu bytes malloc'ed for originalImage.datt\n", vectorSize);
/* allocate memory for gaussImage.data */
printf("vectorSize: %lu\n", vectorSize);
gaussImage.data = farmalloc(vectorSize*sizeof(unsigned char));
if(gaussImage.data == NULL) {
printf("Failed to malloc gaussImage.data\n");
exit(0);
}
printf("%lu bytes malloc'ed for gaussImage.data\n", vectorSize);
/* allocate memory for edgeImage.data */
printf("vectorSize: %lu\n", vectorSize);
edgeImage.data = farmalloc(vectorSize*sizeof(unsigned char));
if(edgeImage.data == NULL) {
printf("Failed to malloc edgeImage.data\n");
exit(0);
}
printf("%lu bytes malloc'ed for edgeImage.data\n", vectorSize);
copyImageInfo(bmpInput, bmpOutput);
copyColorTable(bmpInput, bmpOutput, nColors);
fseek(bmpInput, (14+40+4*nColors), SEEK_SET);
fseek(bmpOutput, (14+40+4*nColors), SEEK_SET);
/* Read input.bmp and store it's raster data into originalImage.data */
for(row=0; row<=originalImage.rows-1; row++) {
for(col=0; col<=originalImage.cols-1; col++) {
fread(pChar, sizeof(char), 1, bmpInput);
*(originalImage.data + row*originalImage.cols + col) = *pChar;
}
}
/**********************************************
* GAUSSIAN MASK ALGORITHM STARTS HERE
***********************************************/
for(r=0; r<=(originalImage.rows-1); r++) {
for(c=0; c<=(originalImage.cols-1); c++) {
SUM = 0;
/* image boundaries */
if(r==0 || r==1 || r==originalImage.rows-2 || r==originalImage.
rows-1)
SUM = *pChar;
else if(c==0 || c==1 || c==originalImage.cols-2 ||
c==originalImage.co
ls-1)
SUM = *pChar;
/* Convolution starts here */
else {
for(I=-2; I<=2; I++) {
for(J=-2; J<=2; J++) {
SUM = SUM + (int)( (*(originalImage.data + c + I + (r +
J)*originalI
mage.cols)) * (MASK[I+2][J+2])/115);
}
}
}
if(SUM>255) SUM=255;
if(SUM<0) SUM=0;
*(gaussImage.data + c + r*originalImage.cols) = (unsigned
char)(SUM);
}
}
/***********************************************
* SOBEL GRADIENT APPROXIMATION
***********************************************/
for(r=0; r<=(originalImage.rows-1); r++) {
for(c=0; c<=(originalImage.cols-1); c++) {
sumX = 0;
sumY = 0;
/* image boundaries */
if(r==0 || r==originalImage.rows-1)
SUM = 0;
else if(c==0 || c==originalImage.cols-1)
SUM = 0;
/* Convolution starts here */
else {
/***********************************
* X gradient approximation
***********************************/
for(I=-1; I<=1; I++) {
for(J=-1; J<=1; J++) {
sumX = sumX + (int)( (*(gaussImage.data + c + I + (r +
J)*originalIm
age.cols)) * GX[I+1][J+1]);
}
}
/**************************
* Y gradient approximation
**************************/
for(I=-1; I<=1; I++) {
for(J=-1; J<=1; J++) {
sumY = sumY + (int)( (*(gaussImage.data + c + I + (r +
J)*originalI
mage.cols)) * GY[I+1][J+1]);
}
}
/***********************************************
* GRADIENT MAGNITUDE APPROXIMATION (Myler p.218)
***********************************************/
SUM = abs(sumX) + abs(sumY);
if(SUM>255) SUM=255;
if(SUM<0) SUM=0;
fprintf(dataOutput, "MAG: \t\t\t%d\n", SUM);
/***************************
* Magnitude orientation
***************************/
/* Cannot divide by zero*/
if(sumX == 0) {
if(sumY==0) ORIENT = 0.0;
else if (sumY<0) {
sumY = -sumY;
ORIENT = 90.0;
}
else ORIENT = 90.0;
}
/* Can't take invtan of angle in 2nd Quad */
else if(sumX<0 && sumY>0) {
sumX = -sumX;
ORIENT = 180 - ((atan((float)(sumY)/(float)(sumX))) *
(180/M_PI));
}
/* Can't take invtan of angle in 4th Quad */
else if(sumX>0 && sumY<0) {
sumY = -sumY;
ORIENT = 180 - ((atan((float)(sumY)/(float)(sumX))) *
(180/M_PI));
}
/* else angle is in 1st or 3rd Quad */
else ORIENT = (atan((float)(sumY)/(float)(sumX))) * (180/M_PI);
fprintf(dataOutput, "Gradient Dir: \t\t%f\n", ORIENT);
/***************************************************
* Find edgeDirection by assigning ORIENT a value of
* either 0, 45, 90 or 135 degrees, depending on which
* value ORIENT is closest to
****************************************************/
if(ORIENT < 22.5) edgeDirection = 0;
else if(ORIENT < 67.5) edgeDirection = 45;
else if(ORIENT < 112.5) edgeDirection = 90;
else if(ORIENT < 157.5) edgeDirection = 135;
else edgeDirection = 0;
fprintf(dataOutput, "Edge Dir: \t\t%d\t\t(row: %d, col %d)\n",
edgeD
irection, r, c);
/***************************************************
* Obtain values of 2 adjacent pixels in edge
* direction.
****************************************************/
if(edgeDirection == 0) {
leftPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c - 1));
rightPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c + 1));
}
else if(edgeDirection == 45) {
leftPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c + origin
alImage.cols - 1));
rightPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c - origi
nalImage.cols + 1));
}
else if(edgeDirection == 90) {
leftPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c - origin
alImage.cols));
rightPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c + origi
nalImage.cols));
}
else {
leftPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c - origin
alImage.cols - 1));
rightPixel = (int)(*(gaussImage.data + r*originalImage.cols
+
c + origi
nalImage.cols + 1));
}
fprintf(dataOutput, "Left: %d, Right: %d\n", leftPixel,
rightPixel);
/*********************************************
* Compare current magnitude to both adjacent
* pixel values. And if it is less than either
* of the 2 adjacent values - suppress it and make
* a nonedge.
*********************************************/
if(SUM < leftPixel || SUM < rightPixel) SUM = 0;
else {
/**********************
* Hysteresis
**********************/
highThreshold = 120;
lowThreshold = 40;
if(SUM >= highThreshold)
SUM = 255; /* edge */
else if(SUM < lowThreshold)
SUM = 0; /* nonedge */
/* SUM is between T1 & T2 */
else {
/* Determine values of neighboring pixels */
P1 = (int)(*(gaussImage.data + r*originalImage.cols + c -
originalIma
ge.cols - 1));
P2 = (int)(*(gaussImage.data + r*originalImage.cols + c
- originalImage.cols));
P3 = (int)(*(gaussImage.data + r*originalImage.cols + c
- originalImage.cols + 1));
P4 = (int)(*(gaussImage.data + r*originalImage.cols + c
- 1));
P5 = (int)(*(gaussImage.data + r*originalImage.cols + c
+ 1));
P6 = (int)(*(gaussImage.data + r*originalImage.cols + c
+ originalImage.cols - 1));
P7 = (int)(*(gaussImage.data + r*originalImage.cols + c
+ originalImage.cols));
P8 = (int)(*(gaussImage.data + r*originalImage.cols + c
+ originalImage.cols + 1));
/* Check to see if neighboring pixel values are edges */
if(P1 > highThreshold || P2 > highThreshold || P3 >
highThreshold ||
P4 > highThreshold || P5 > highThreshold || P6 > highThreshold || P7 > highT
hreshold || P8 > highThreshold)
SUM = 255; /* make edge */
else SUM = 0; /* make nonedge */
}
}
} /* else loop ends here (starts after b.c.) */
fprintf(dataOutput, "New MAG: \t\t%d\n\n\n\n", SUM);
*(edgeImage.data + c + r*originalImage.cols) = 255 - (unsigned
char)(S
UM);
fwrite( (edgeImage.data + c + r*originalImage.cols), sizeof(char)
,
1,
bmpOutput);
}
}
printf("See canny.bmp for results\n");
fclose(bmpInput);
fclose(bmpOutput);
farfree(gaussImage.data); /* Finished with gaussImage.data */
farfree(edgeImage.data); /* Finished with edgeImage.data */
farfree(originalImage.data); /* Finished with originalImage.data */
return 0;
}
/*----------GET IMAGE INFO SUBPROGRAM--------------*/
long getImageInfo(FILE* inputFile, long offset, int numberOfChars)
{
unsigned char *ptrC;
long value = 0L;
unsigned char dummy;
int i;
dummy = '0';
ptrC = &dummy;
fseek(inputFile, offset, SEEK_SET);
for(i=1; i<=numberOfChars; i++)
{
fread(ptrC, sizeof(char), 1, inputFile);
/* calculate value based on adding bytes */
value = (long)(value + (*ptrC)*(pow(256, (i-1))));
}
return(value);
} /* end of getImageInfo */
/*-------------COPIES HEADER AND INFO HEADER----------------*/
void copyImageInfo(FILE* inputFile, FILE* outputFile)
{
unsigned char *ptrC;
unsigned char dummy;
int i;
dummy = '0';
ptrC = &dummy;
fseek(inputFile, 0L, SEEK_SET);
fseek(outputFile, 0L, SEEK_SET);
for(i=0; i<=50; i++)
{
fread(ptrC, sizeof(char), 1, inputFile);
fwrite(ptrC, sizeof(char), 1, outputFile);
}
}
/*----------------COPIES COLOR TABLE-----------------------------*/
void copyColorTable(FILE* inputFile, FILE* outputFile, int nColors)
{
unsigned char *ptrC;
unsigned char dummy;
int i;
dummy = '0';
ptrC = &dummy;
fseek(inputFile, 54L, SEEK_SET);
fseek(outputFile, 54L, SEEK_SET);
for(i=0; i<=(4*nColors); i++) /* there are (4*nColors) bytesin color tabl
e */
{
fread(ptrC, sizeof(char), 1, inputFile);
fwrite(ptrC, sizeof(char), 1, outputFile);
}
}
--
╔═══════════════════╗
║★★★★★友谊第一 比赛第二★★★★★║
╚═══════════════════╝
※ 来源:．哈工大紫丁香 bbs.hit.edu.cn [FROM: 202.118.229.162]