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