Ellipse detection in image by using Hough transform

How we can detect ellipses in images ? One way is to use Hough transform. I will use Hough transform algorithm variant created by Yonghong Xie and Qiang Ji. That algorithm pseudo-code:

1. Store all edge pixels in a one dimensional array.
2. Clear the accumulator array.
3. For each pixel (x1, y1), carry out the following steps from (4) to (14).
   4. For each other pixel (x2, y2), if the distance between
   (x1, y1) and (x2, y2) is greater than the required least
   distance for a pair of pixels to be considered then
   carry out the following steps from (5) to (14).
      5. From the pair of pixels (x1, y1) and (x2, y2) calculate the center,
      orientation and length of major axis for the assumed ellipse.
      6. For each third pixel (x, y), if the distance between
      (x, y) and (x0, y0) is greater than the required least
      distance for a pair of pixels to be considered then
      carry out the following steps from (7) to (9).
         7. Calculate the length of minor axis.
         8. Increment the accumulator for this length of minor axis by 1.
         9. Loop until all pixels are computed for this pair of pixels.
      10. Find the maximum element in accumulator array. The
      related length is the possible length of minor axis
      for assumed ellipse. If the vote is greater than the
      required least number for assumed ellipse, one
      ellipse is detected.
      11. Output ellipse parameters.
      12. Remove the pixels on the detected ellipse from edge pixel array.
      13. Clear accumulator array.
      14. Loop until all pairs of pixels are computed.

Proof-of-concept algorithm implementation in Python:

---

import sys
from PIL import Image,ImageFilter, ImageDraw
from math import *

# some global constants
EL_COVERAGE_RATIO = 0.9
EL_VERIFICATION_DISTANCE = 1.
EL_PATH_POINTS = 51
MIN_MINOR_FREQUENCY = 30
MIN_HALF_MAJOR = 8
MIN_HALF_MINOR = 6

def distance(p1,p2):
 x1,y1 = p1
 x2,y2 = p2
 return sqrt((x1-x2)**2 + (y1-y2)**2)

def nonnegative(v):
 return v if v >= 0 else 0

def parametricEllipse(center, a, b, angle):
 xc,yc = center
 elx = lambda t: xc + a * cos(t) * cos(angle) - b * sin(t) * sin(angle)
 ely = lambda t: yc + a * cos(t) * sin(angle) + b * sin(t) * cos(angle)
 return [(int(elx(2.*pi*x/float(EL_PATH_POINTS-1))),int(ely(2.*pi*x/float(EL_PATH_POINTS-1)))) for x in range(EL_PATH_POINTS)]

assert len(sys.argv)!=2, "missing input and/or output file !"

im = Image.open(sys.argv[1])
width, height = im.size
io = Image.new('RGB',(width, height),(255,255,255))
draw = ImageDraw.Draw(io)

# converting image to grayscale
im = im.convert('L')

# detecting edge pixels
im = im.filter(ImageFilter.FIND_EDGES)

# converting to binary image
im = im.convert('1')
pixels = im.load()
pxy = []

# extracting binary pixels coordinates
for x in range(width):
 for y in range(height):
  if pixels[x,y]==255:
   pxy.append((x,y))

# applying Hough transform for ellipses detection.
# algorithm is taken from this paper:
# http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.1.8792&rep=rep1&type=pdf
cIx = -1
colors = [(255,0,0),(0,200,0),(0,0,255)]
for x1,y1 in pxy:
  for x2,y2 in pxy:
    bbins = {}
    dist = distance((x1,y1),(x2,y2))
    if dist >= 2*MIN_HALF_MAJOR:
     cent = ((x1+x2)/2.,(y1+y2)/2.)
     a = dist/2. # semi-length of major axis
     alfa = atan2((y2 - y1),(x2 - x1))
     for rx,ry in pxy:
      d = distance((rx,ry),cent)
      if d >= MIN_HALF_MINOR:
       f = distance((rx,ry),(x2,y2))
       cost = (a**2. + d**2. - f**2.)/(0.00001+2.*a*d)
       b = sqrt(nonnegative((a**2. * d**2. * (1.-cost**2.))/(0.00001 + a**2. - d**2. * cost**2.)))  # semi-length of minor axis
       b = int(b)
       if bbins.has_key(b):
        bbins[b]+=1
       elif b > 0:
        bbins[b]=1
     bbins_rev = dict([(v,k) for k,v in bbins.iteritems()])
     max_freq = max(bbins_rev.keys())
     bmax = bbins_rev[max_freq]
     # Did we found probable ellipse ?
     if max_freq >= MIN_MINOR_FREQUENCY and alfa >=0.0 and bmax >= MIN_HALF_MINOR:
      elData = parametricEllipse(cent, a, bmax, alfa)
      supported = []
      supportRatio = 0.0
      # counting how much pixels lies on ellipse path
      for i in range(EL_PATH_POINTS):
       elx,ely = elData[i]
       added = False
       for x,y in pxy:
        if distance((elx,ely),(x,y)) <= EL_VERIFICATION_DISTANCE:
         supported.append((x,y))
         if not added:
          supportRatio += 1./float(EL_PATH_POINTS)
          added = True
      supported = list(set(supported))
      # if number of pixels on ellipse path is big enough
      if supportRatio >= EL_COVERAGE_RATIO:
       cIx = (cIx+1)%3
       print "coverage %.2f" % supportRatio,"frequency ", max_freq, "center ", cent, "angle %.2f" % alfa, "axes (%.2f,%.2f)" % (a, bmax)
       # removing founded ellipse pixels from further analysis
       for p in supported:
        pxy.remove(p)
       # drawing founded ellipse
       for i in range(EL_PATH_POINTS):
        elx,ely = elData[i]
        if i < EL_PATH_POINTS-1:
         draw.line(elData[i] + elData[i+1], fill=colors[cIx])
io.save(sys.argv[2])
print "***************************************************************"
print "************************** DONE *******************************"
print "***************************************************************"

---

(Prototype algorithm is slow, tested only on 50x50 images). So, by running this algo on this image:

we will get such algorithm output:

Have fun in computer vision !

Algorithm to determine image contrast

How to determine image contrast ? More specifically - How to determine that image contrast is low and it needs automatic contrast adjustment through histogram equalization method ?

Algorithm is this:
1. Calculate cumulative histogram of image.
2. Make linear regression of cumulative histogram in the form freq(x) = A*x + B.
3. Calculate RMSE of real_cumulative_frequency(x)-freq(x).
4. If that RMSE is close to zero - image is already equalized and should be in good contrast. (That means for equalized images cumulative histograms must be linear)

Because picture is worth a thousand words - this is how cumulative histograms looks for image with different levels of contrast:

So as you see - image with most aggressive contrast has cumulative histogram which is almost a perfect line.

Now fun part - C code which opens image in PGM format with ASCII encoding and detects if image needs contrast adjustment or not :

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define MAX_COLS 1000
#define MAX_ROWS 1000

typedef struct {
 int Width;
 int Height;
 int max_value;
 int data[MAX_ROWS][MAX_COLS];
} PGMdata;

void readPgmFile(char* fileName, PGMdata * pgmOut) {
 FILE * pFile;
 char line[1000];
 char* res;
 int lineNum = 0;

 if ((pFile = fopen(fileName , "r")) == NULL)
  printf("Error opening file %s \n", fileName);
 else {
  while ((res = fgets(line, 100, pFile)) != NULL) {
  lineNum++;
  // max value of pixel
  if (lineNum==4)
   sscanf(line,"%i",&pgmOut->max_value);
  // width and height of image
  if (lineNum==3) {
   sscanf(line,"%i %i",&pgmOut->Width, &pgmOut->Height);
  }
  // load real pixels
  if (lineNum > 4) {
   int ix = lineNum - 5;
   int row = ix/pgmOut->Width;
   int col = ix - row*pgmOut->Width;
   sscanf(line,"%i", &pgmOut->data[row][col]);
  }
  }
  fclose(pFile);
    }
}

void CumulativeHistogram(PGMdata* image, double* hist) {
 int i,j;
 double dp = 1.0/((double) image->Width*image->Height);

 // initializing histogram bins to zero
 for (i=0; i<256; i++) {
  hist[i] = 0.0;
 }

 // calculating histogram
 for (i=0; i<image->Width; i++) {
  for (j=0; j<image->Height; j++) {
   hist[image->data[j][i]] += dp;
  }
 }

 // making histogram cumulative
 for (i=0; i<255; i++) {
  hist[i+1] += hist[i];
 }
}

void CreateXmatrix(double mat[256][2]) {
 int i;
 for (i=0; i<256; i++) {
   mat[i][0] = 1;
   mat[i][1] = i;
 }
}

void TransposeMatrix(double mat[256][2], double tmat[2][256]) {
 int i,j;
 for (i=0; i<2; i++) {
  for (j=0; j<256; j++) {
   tmat[i][j] = mat[j][i];
  }
 }
}

double MultiplyMatrixes(double A[2][256], double B[256][2], int row, int col) {
 int i;
 double sum = 0.0;

 for (i=0; i<256; i++) {
  sum += A[row][i]*B[i][col];
 }

 return sum;
}

double MultiplyMatrixAndVector(double A[2][256], double Y[256], int row) {
 int i;
 double sum = 0.0;

 for (i=0; i<256; i++) {
  sum += A[row][i]*Y[i];
 }

 return sum;
}

double HistogramPredicted(double c0, double c1, double level) {
 return c0 + c1*level;
}

double RootMeanSquare(double* hist, double c0, double c1) {
 double rms = 0.0;
 int i;

 for (i=0; i<256; i++) {
  rms += pow(hist[i]-HistogramPredicted(c0,c1,i),2.0);
 }
 rms /= 256.0;

 return sqrt(rms);
}

void LinearLeastSquares(double* hist, double* c0, double* c1) {
 double X[256][2];
 double tX[2][256];
 double a1,a2,a3,a4;
 double b1, b2;

 // create matrix X composed of x'es
 CreateXmatrix(X);

 // transpose X matrix
 TransposeMatrix(X,tX);

 // calculate tX*X matrix which is
 // [a1 a2]
 // [a3 a4]

 a1 = MultiplyMatrixes(tX,X,0,0);
 a2 = MultiplyMatrixes(tX,X,0,1);
 a3 = MultiplyMatrixes(tX,X,1,0);
 a4 = MultiplyMatrixes(tX,X,1,1);

 // calculate tX*Y  (Y=HISTOGRAM) which is
 // [b1]
 // [b2]

 b1 = MultiplyMatrixAndVector(tX,hist,0);
 b2 = MultiplyMatrixAndVector(tX,hist,1);

 // solve matrix equation by using elimination of variables method
 // [a1 a2]   [c0]   [b1]
 //         *      =
 // [a3 a4]   [c1]   [b2]

 *c1 = (a1*b2-a3*b1)/(a1*a4-a2*a3);
 *c0 = (b1-a2*(*c1))/a1;

}

int main () {
 double hist[256];
 static PGMdata image = {0};
 double c0=0.0, c1=0.0;
 double rms;
 char pgm_file[20] = {0};

 while (1) {
   // get PGM file name
   printf("PGM filename: ");
   scanf("%s", pgm_file);

   // read grayscale image from PGM format
   memset(&image, 0, sizeof(image));
   readPgmFile(pgm_file, &image);
   if (image.Width == 0)
     continue;

   // create cumulative histogram
   CumulativeHistogram(&image, hist);

   // least mean squares method, to find c0,c1 coefficents in equation:
   // c0 + c1*gray_level = Frequency
   LinearLeastSquares(hist,&c0,&c1);

   // calculate RMS of Frequency[bin]-predicted_Frequency(bin)
   rms = RootMeanSquare(hist,c0,c1);

   // Low RMS shows that histogram frequencies are distributed in linear fashion
   // between bins. This means that histogram equalization was performed on image
   // and/or that global image contrast is good.
   if (rms <= 0.01)
    printf("\n==> %s contrast is OK (rms=%f)\n\n", pgm_file ,rms);
   else
    printf("\n==> %s contrast is BAD (rms=%f)\n\n", pgm_file ,rms);
 }

 return 0;
}

If you want to test this algorithm on ASCII PGM image - you better covert image to PGM format with GIMP 2.6.11, because it was only tested in this way.
So by running this C code on below image (after converting it to PGM)

we get algorithm output
"image contrast is BAD ... " :-)

Detecting copy-move forgery in images

This time I'll talk about digital image forensic and how to detect copy-move forgery in images. I`ve implemented some ad-hoc algorithm for detection of this kind of forgeries. Algorithm is robust, so it can detect forgeries in lossy image formats- such as JPEG. Please keep in mind - this algorithm is experimental toy, if you want more general solutions - you should read this paper or maybe this to get some ideas in the field. BTW, this algorithm is only tested with Python 2.6, so you better install Python 2.6.

Robust detection algorithm steps
1. Blur image for eliminating image details
2. Convert image to degraded palette
3. Decompose image into small NxN pixel blocks
4. Alphabetically order these blocks by their pixel values
5. Extract only these adjacent blocks which have small absolute color difference
6. Cluster these blocks into clusters by intersection area among blocks
7. Extract only these clusters which are bigger than block size
8. Extract only these clusters which have similar cluster, by using some sort of similarity function (in this case Hausdorff distance between clusters)
9. Draw discovered similar clusters on image

Copy-move detection script usage instructions
Run script from the command line - cmd in Windows or terminal in Linux.
1. At first try to execute script with default parameters-
python detect_copymove.py image_file.jpg

2. Then try to lower block color deviation threshold-
python detect_copymove.py image_file.jpg --blcoldev=0.05

3. Finally - run script in manual mode and try to spot similar regions by eyes-
python detect_copymove.py image_file.jpg --blcoldev=0.05 --imauto=0

If by trying all 3 steps - no copy-move tamperings revealed - there is a good chance that really there are no such tamperings in image. BTW, script has more parameters, full list of them -
python detect_copymove.py --help

Experiments - object cloning forgery
This type of forgery tries to deceive viewer that there were more objects than really was.
Original picture:

Doctored image:

Copy-move script output:

As you see script detected cloning of one tank, but did not detected cloning of tank in further row. This can be because of small size of that tank and because that tank is on boundary of image.

Experiments - object hiding forgery
In this type of forgery viewer is deceived that there were less objects than in reality was.
Original image:

Doctored image:

Copy-move script output:


Experiments - mixed type forgery
In this type of forgery some objects was hidden, some - cloned.
Original picture:

Doctored image:

Copy-move script output:

As you see script detected 2 forgeries - cloning of dog, and cloning of background (which was done with intention to hide other dog)

Conclusions
Copy-move forgery detection script works not bad. However it is also not perfect of course (may not detect all tamperings, so called false negative rate is not zero of course). And I`ve not tested very well, but I suppose script also have false positive rate - i.e. can detect areas which actually was not tampered. But all in all - I think this algorithm is pretty obvious and can be simply implemented and used as first line of defense against copy-move forgeries.

Have fun in copy-move forgeries hunting !!

Detecting similar entropy zones in image

As I thought about information entropy one idea came to me - to write application which looks similar entropy zones in image. So after some time, I came with this algorithm (pseudo code):

1. Split image into 5x5 pixel image blocks
2. Calculate information entropy of these blocks (actually sum of entropy in 3 color channels)
3. Find similar entropy blocks.
4. [...] filter out small groups of blocks (seems like noise, huh ? ). Blah, blah...
5. Picture these similar entropy blocks on top of original image as red color layer.

Now real part (as you know already) - Python code which does the job [you need PIL module to run this]:

---

from PIL import Image
from math import *

def entropysum(pixels):
 """
 Calculating information entropy for image
 region and returning entropy sum for
 all 3 color channels
 """
 cr = [r for (r,g,b) in pixels]
 cg = [g for (r,g,b) in pixels]
 cb = [b for (r,g,b) in pixels]
 
 er = 0.0
 eg = 0.0
 eb = 0.0
 
 for r in set(cr):
  p = float(cr.count(r))/len(cr)
  if p > 0.0: er += -p * log(p,2)

 for g in set(cg):
  p = float(cg.count(g))/len(cg)
  if p > 0.0: eg += -p * log(p,2)

 for b in set(cb):
  p = float(cb.count(b))/len(cb)
  if p > 0.0: eb += -p * log(p,2)

 return er + eg + eb

def decompose(image, block_len):
 """
 Decomposing given image into some number of
 smaller images of size block_len*block_len
 """
 parts = []
 w, h = image.size
 
 for x in range(0, w, block_len):
  for y in range(0, h, block_len):
   locim = image.crop((x,y,x+block_len,y+block_len))
   acc = entropysum(list(locim.getdata()))
   parts.append((acc,x,y,locim))
 
 parts.sort()
 
 return parts

def similarparts(imagparts):
 """
 Detecting similar image blocks by comparing
 entropy of given images. Two images considered
 being equal if entropy difference is not big.
 """
 dupl = []
 
 for i in range(len(imagparts)-1):
  acc1, x1, y1, im1 = imagparts[i]
  acc2, x2, y2, im2 = imagparts[i+1]
  
  if acc1 == acc2 == 0:
   gain = 0.0
  else:
   gain = 100.0 * (1.0 - acc1 / acc2)

  if 0.01 < gain < 0.1 :
   if imagparts[i] not in dupl:
    dupl.append(imagparts[i])
   if imagparts[i+1] not in dupl:
    dupl.append(imagparts[i+1])

 return dupl

def clusterparts(parts):
 """
 Grouping nearest images into groups.
 This is done, because we need to
 filter out very small groups. We
 want to know only big differences.
 """
 
 filtparts = []
 clust = {}
 belongs = {}
 w,h = parts[0][3].size
 
 # assign all parts to clusters
 for i in range(len(parts)):
  acc, x, y, im = parts[i]
  sides = []
  sides.append(str(x)+str(y)+str(x+w)+str(y))
  sides.append(str(x+w)+str(y)+str(x+w)+str(y+h))
  sides.append(str(x)+str(y+h)+str(x+w)+str(y+h))
  sides.append(str(x)+str(y)+str(x)+str(y+h))
  
  # detect side already in cluster
  fc = None
  for s in sides:
   if belongs.has_key(s):
    fc = belongs[s]
    break
  
  # if this is new cluster
  if fc == None:
   fc = len(clust) + 1
   clust[fc] = 1
  else:
   clust[fc] += 1
  
  # set cluster for rectangle sides
  for s in sides:
   if not belongs.has_key(s):
    belongs[s] = fc

 # filter out small clusters
 for i in range(len(parts)):
  acc, x, y, im = parts[i]
  side = str(x)+str(y)+str(x+w)+str(y)
  cl = belongs[side]
  if clust[cl] > 2:
   filtparts.append(parts[i])
 
 return filtparts

def marksimilar(image, dparts):
 """
 Mark found similar image blocks on
 original image, by applying red layer
 on similar parts of image.
 """
 if dparts:
  colormask = Image.new('RGB', dparts[0][3].size,(255,0,0))
  for (acc,x,y,im) in dparts:
   im = Image.blend(im, colormask, 0.4)
   image.paste(im,(x,y))
 
 return image

if __name__ == '__main__':
 im = Image.open("1.jpg")
 ls = decompose(im, 5)
 dparts = similarparts(ls)
 cparts = clusterparts(dparts)
 im = marksimilar(im, cparts)
 im.show()

---

So these are the results after running this script on several images:





Conclusion

So this algorithm is an interesting tool for exploration of information entropy in image. Maybe in some cases it could be a tool for analyzing very similar texture zones. BTW information entropy may be used for hashing image. Hashing image is useful, because it lets us to search similar images in database (for example) by its hash.

Have fun !