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 ... " :-)
Hi, I tried to experiment the above code with the ballons.pgm data file which presents in http://orion.math.iastate.edu/burkardt/data/pgm/pgm.html link. when it tries to calculate a histogram in cumulativehistogram function, it crashes. Am not sure about the reason, could you plz comment on the same? Thanks in advance, Veeru.
ReplyDeleteThanks for note. There was a problem that too big image not fitted into stack. I fixed that issue by moving image struct to static area of program data. Try it now - it should work.
ReplyDelete