stevesketch
3/25/2016 - 2:15 AM

#PImage

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);
      }
    }
  }
}