PImage sourceImage;
int kernelSize = 3;
float[] gaussianKernel = generateGaussianKernel(kernelSize,0.1);
float []theta;
void setup()
{
sourceImage = loadImage("Reindeer02.jpg");
size(1024,1024);
theta = new float[sourceImage.width*sourceImage.height];
noLoop();
}
void draw()
{
PImage gray = graying(sourceImage);
PImage blur = gaussianBlur(gray,kernelSize,gaussianKernel);
PImage gradient = calculateGradient(blur);
PImage suppressed = nonMaximumSuppression(gradient);
//image(sourceImage,0,0);
//image(gray,0,0);
//image(blur,0,0);
//image(gradient,0,0);
//image(suppressed,0,0);
calculateThresholds(suppressed,gradient);
PImage edge = edgeTracing(suppressed,gradient);
image(edge,0,0);
}
PImage graying(PImage colorfulImage)
{
int imageWidth = colorfulImage.width;
int imageHeight = colorfulImage.height;
PImage grayImage = createImage(imageWidth,imageHeight,ALPHA);
for(int y=0; y<imageHeight; y++)
{
for(int x=0; x<imageWidth; x++)
{
int pixelPosition = y*imageWidth + x;
int r = (colorfulImage.pixels[pixelPosition] >> 16) & 0xFF;
int g = (colorfulImage.pixels[pixelPosition] >> 8) & 0xFF;
int b = colorfulImage.pixels[pixelPosition] & 0xFF;
grayImage.pixels[pixelPosition] = color(0.072169*r + 0.715160*g + 0.212671*b);
}
}
return grayImage;
}
float[] generateGaussianKernel(int size, float sigma)
{
float[] gaussianKernel = new float[size];
float sigma22 = 2 * sigma * sigma;
float pi2 = 2 * (float)Math.PI;
float sqrtPi2Sigma = (float)Math.sqrt(pi2) * sigma;
int halfLength = size/2;
int index = 0;
float sum = 0;
for(int i=-halfLength;i<=halfLength;i++)
{
float distance = i*i;
gaussianKernel[index] = (float)Math.exp((-distance)/sigma22)/sqrtPi2Sigma;
sum += gaussianKernel[index];
index++;
}
for(int i=0;i<size;i++)
{
gaussianKernel[i] /= sum;
}
return gaussianKernel;
}
PImage gaussianBlur(PImage grayImage, int kernelSize, float[] gaussianKernel)
{
float totalOfKernel = 0;
int halfLength = kernelSize/2;
int imageWidth = grayImage.width;
int imageHeight = grayImage.height;
PImage blurImage = createImage(imageWidth,imageHeight,ALPHA);
for(int y=0;y<imageHeight;y++)
{
for(int x=0;x<imageWidth;x++)
{
if(y<=halfLength | y>= imageHeight-halfLength | x<=halfLength | x>= imageWidth-halfLength)
{
int pixelPosition = y*imageWidth + x;
blurImage.pixels[pixelPosition] = grayImage.pixels[pixelPosition];
}
}
}
for(int y=halfLength; y<imageHeight-halfLength; y++)
{
for(int x=halfLength; x<imageWidth-halfLength; x++)
{
float rsum = 0;
for(int i=-halfLength;i<=halfLength;i++)
{
int pixelPosition = y*imageWidth + x + i;
rsum += gaussianKernel[i+halfLength] * red(grayImage.pixels[pixelPosition]);
}
blurImage.pixels[y*imageWidth + x] = color(rsum);
}
}
for(int y=halfLength; y<imageHeight-halfLength; y++)
{
for(int x=halfLength; x<imageWidth-halfLength; x++)
{
float rsum = 0;
for(int i=-halfLength;i<=halfLength;i++)
{
int pixelPosition = (y+i)*imageWidth + x;
rsum += gaussianKernel[i+halfLength] * red(blurImage.pixels[pixelPosition]);;
}
blurImage.pixels[y*imageWidth + x] = color(rsum);
}
}
return blurImage;
}
PImage calculateGradient(PImage blurImage)
{
int imageWidth = blurImage.width;
int imageHeight = blurImage.height;
float []gradientX = new float[imageWidth*imageHeight];
float []gradientY = new float[imageWidth*imageHeight];
PImage gradientImage = createImage(imageWidth,imageHeight,ALPHA);
for(int y=0;y<imageHeight;y++)
{
for(int x=0;x<imageWidth;x++)
{
if(y<1 | y>imageHeight-2 | x<1 | x>imageWidth-2)
{
int pixelPosition = y*imageWidth + x;
gradientImage.pixels[pixelPosition] = color(0);
}
}
}
for(int y=1; y<imageHeight-1; y++)
{
for(int x=1; x<imageWidth-1; x++)
{
int pixelPosition = y*imageWidth + x;
gradientX[pixelPosition] = -2*red(blurImage.pixels[pixelPosition-1]) + 2*red(blurImage.pixels[pixelPosition+1])
-1*red(blurImage.pixels[pixelPosition-1-imageWidth]) + 1*red(blurImage.pixels[pixelPosition+1-imageWidth])
-1*red(blurImage.pixels[pixelPosition-1+imageWidth]) + 1*red(blurImage.pixels[pixelPosition+1+imageWidth]);
gradientY[pixelPosition] = -2*red(blurImage.pixels[pixelPosition+imageWidth]) + 2*red(blurImage.pixels[pixelPosition-imageWidth])
-1*red(blurImage.pixels[pixelPosition+1+imageWidth]) + 1*red(blurImage.pixels[pixelPosition+1-imageWidth])
-1*red(blurImage.pixels[pixelPosition-1+imageWidth]) + 1*red(blurImage.pixels[pixelPosition-1-imageWidth]);
gradientImage.pixels[pixelPosition] = color(sqrt(red(color(gradientX[pixelPosition]))*red(color(gradientX[pixelPosition]))
+ red(color(gradientY[pixelPosition]))*red(color(gradientY[pixelPosition]))));
theta[pixelPosition] = atan2(gradientY[pixelPosition],gradientX[pixelPosition]);
if(theta[pixelPosition]<0)theta[pixelPosition] = theta[pixelPosition]+PI;
}
}
return gradientImage;
}
PImage calculateLaplacian(PImage blurImage)
{
int imageWidth = blurImage.width;
int imageHeight = blurImage.height;
PImage laplacianImage = createImage(imageWidth,imageHeight,ALPHA);
for(int y=0;y<imageHeight;y++)
{
for(int x=0;x<imageWidth;x++)
{
if(y<1 | y>imageHeight-2 | x<1 | x>imageWidth-2)
{
int pixelPosition = y*imageWidth + x;
laplacianImage.pixels[pixelPosition] = color(0);
}
}
}
for(int y=1; y<imageHeight-1; y++)
{
for(int x=1; x<imageWidth-1; x++)
{
int pixelPosition = y*imageWidth + x;
laplacianImage.pixels[pixelPosition] = color(4*red(blurImage.pixels[pixelPosition]) -1*red(blurImage.pixels[pixelPosition+1])
-1*red(blurImage.pixels[pixelPosition-1]) -1*red(blurImage.pixels[pixelPosition-imageWidth])
-1*red(blurImage.pixels[pixelPosition+imageWidth]));
}
}
return laplacianImage;
}
PImage nonMaximumSuppression(PImage gradientImage)
{
int imageWidth = gradientImage.width;
int imageHeight = gradientImage.height;
PImage suppressedImage = createImage(imageWidth,imageHeight,ALPHA);
float g1,g2,g3,g4,weight;
float mean12=0;
float mean34=0;
for(int y=0;y<imageHeight;y++)
{
for(int x=0;x<imageWidth;x++)
{
if(y<1 | y>imageHeight-2 | x<1 | x>imageWidth-2)
{
int pixelPosition = y*imageWidth + x;
suppressedImage.pixels[pixelPosition] = color(0);
}
}
}
for(int y=1; y<imageHeight-1; y++)
{
for(int x=1; x<imageWidth-1; x++)
{
int pixelPosition = y*imageWidth + x;
if(red(gradientImage.pixels[pixelPosition]) == 0)
{
suppressedImage.pixels[pixelPosition] = color(0);
}
else
{
if(theta[pixelPosition]>=0 & theta[pixelPosition]<=PI/4)
{
g1 = red(gradientImage.pixels[pixelPosition+1-imageWidth]);
g2 = red(gradientImage.pixels[pixelPosition+1]);
g3 = red(gradientImage.pixels[pixelPosition-1]);
g4 = red(gradientImage.pixels[pixelPosition-1+imageWidth]);
weight = tan(theta[pixelPosition]);
mean12 = g1*weight + g2*(1-weight);
mean34 = g4*weight + g3*(1-weight);
} else if(theta[pixelPosition]>PI/4 & theta[pixelPosition]<=PI/2)
{
g1 = red(gradientImage.pixels[pixelPosition-imageWidth]);
g2 = red(gradientImage.pixels[pixelPosition+1-imageWidth]);
g3 = red(gradientImage.pixels[pixelPosition-1+imageWidth]);
g4 = red(gradientImage.pixels[pixelPosition+imageWidth]);
weight = 1/tan(theta[pixelPosition]);
mean12 = g1*(1-weight) + g2*(weight);
mean34 = g4*(1-weight) + g3*(weight);
} else if(theta[pixelPosition]>PI/2 & theta[pixelPosition]<=3*PI/4)
{
g1 = red(gradientImage.pixels[pixelPosition-imageWidth-1]);
g2 = red(gradientImage.pixels[pixelPosition-imageWidth]);
g3 = red(gradientImage.pixels[pixelPosition+imageWidth]);
g4 = red(gradientImage.pixels[pixelPosition+imageWidth+1]);
weight = abs(1/tan(theta[pixelPosition]));
mean12 = g1*(weight) + g2*(1-weight);
mean34 = g4*(weight) + g3*(1-weight);
} else if(theta[pixelPosition]>3*PI/4 & theta[pixelPosition]<=PI)
{
g1 = red(gradientImage.pixels[pixelPosition-imageWidth-1]);
g2 = red(gradientImage.pixels[pixelPosition-1]);
g3 = red(gradientImage.pixels[pixelPosition+1]);
g4 = red(gradientImage.pixels[pixelPosition+imageWidth+1]);
weight = abs(tan(theta[pixelPosition]));
mean12 = g1*(weight) + g2*(1-weight);
mean34 = g4*(weight) + g3*(1-weight);
}
if(red(gradientImage.pixels[pixelPosition]) >= mean12 & red(gradientImage.pixels[pixelPosition]) >= mean34)
{
suppressedImage.pixels[pixelPosition] = color(128);
}
else
{
suppressedImage.pixels[pixelPosition] = color(0);
}
}
}
}
return suppressedImage;
}
int highThreshold;
int lowThreshold;
void calculateThresholds(PImage suppressedImage, PImage gradientImage)
{
int[] hist = new int[256];
int edgePoints = 0;
int maxGradientMag = 0;
int totalPixels = suppressedImage.width * suppressedImage.height;
for(int i=0;i<256;i++)
{
hist[i]=0;
}
for(int i=0;i<totalPixels;i++)
{
if(128 == red(suppressedImage.pixels[i]))
{
hist[(int)red(gradientImage.pixels[i])]++;
edgePoints++;
maxGradientMag = (red(gradientImage.pixels[i]) > maxGradientMag)?(int)red(gradientImage.pixels[i]):maxGradientMag;
}
}
int highCount = (int)(edgePoints * 0.79 + 0.5);
edgePoints = 0;
highThreshold = 0;
while(edgePoints < highCount)
{
edgePoints += hist[highThreshold];
highThreshold++;
}
highThreshold--;
lowThreshold = (int)(highThreshold * 0.5 + 0.5);
}
PImage edgeTracing(PImage suppressedImage, PImage gradientImage)
{
int imageWidth = suppressedImage.width;
int imageHeight = suppressedImage.height;
int totalPixels = imageWidth*imageHeight;
PImage edgeImage = createImage(imageWidth,imageHeight,ALPHA);
for(int y=1;y<imageHeight-1;y++)
{
for(int x=1;x<imageWidth-1;x++)
{
int pixelPosition = y*imageWidth + x;
if(128 == red(suppressedImage.pixels[pixelPosition]) && red(gradientImage.pixels[pixelPosition])>highThreshold)
{
edgeImage.pixels[pixelPosition] = color(255);
trace(x,y,suppressedImage,gradientImage,edgeImage);
}
}
}
return edgeImage;
}
void trace(int x,int y,PImage suppressedImage,PImage gradientImage,PImage edgeImage)
{
int imageWidth = suppressedImage.width;
int imageHeight = suppressedImage.height;
for (int i=-1;i<=1;i++)
{
for(int j=-1;j<=1;j++)
{
int newy = y+i;
int newx = x+j;
int pixelPosition = newy*imageWidth + newx;
if(128 == red(suppressedImage.pixels[pixelPosition]) && red(gradientImage.pixels[pixelPosition])>lowThreshold
&& newy!=y && newx!=x && 255!=red(edgeImage.pixels[pixelPosition]))
{
edgeImage.pixels[pixelPosition] = color(255);
trace(newx,newy,suppressedImage,gradientImage,edgeImage);
}
}
}
}