Digital Image Processing Assignment III

Main Contents

• Gamma correction;
• Visibility enhancement;
• Histogram equalization.

Step One: Define Some Functions

We add a function that generate a blank color image from a given image,

// generate a blank color image
void gen_color(BITMAP *bmImg, BITMAP *bmColor) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    bmColor->bmHeader = bmImg->bmHeader;
    bmColor->bmInfoSize = bmImg->bmInfoSize;
    bmColor->bmInfo = (BITMAPINFO *) malloc(bmColor->bmInfoSize);
    bmColor->bmInfo->bmiHeader = *bmiHeader;
    init_bmp(bmColor);
}

Step Two: Gamma Correction

Gamma encoding of images is used to optimize the usage of bits when encoding an image, or bandwidth used to transport an image, by taking advantage of the non-linear manner in which humans perceive light and color. Gamma correction is a nonlinear operation used to encode and decode luminance values in video or still image systems. Gamma correction is, in the simplest cases, defined by the following power-law expression.

$$
L_d=L_w^{1 \over \gamma}
$$

void gamma_correct(BITMAP *bmImg, BITMAP *bmGamma, double gamma) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    gen_color(bmImg, bmGamma);
    uint8_t *p = (uint8_t *) malloc(1 << 11);
    for (int16_t i = 0; i < 256; ++i) {
        p[i] = adjust(255 * pow(i / 255., 1 / gamma));
    }
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmImg->bmBytesPerRow + w * bmImg->bmBytesPerPel;
            bmGamma->bmData[pos] = p[bmImg->bmData[pos]];
            bmGamma->bmData[pos + 1] = p[bmImg->bmData[pos + 1]];
            bmGamma->bmData[pos + 2] = p[bmImg->bmData[pos + 2]];
        }
    }
    free(p);
}

Step Three: Visibility Enhancement

We use a logarithmic operator to adjust the pixel value,

$$
L_d={\log(L_w+1) \over \log(L_{max}+1)}
$$

where $L_d$ refers to display luminance, $L_w$ refers to original luminance, and $L_{max}$ is the maximal luminance in the original image.

This mapping function make sure that, no matter the dynamic range of the scene, the maximal luminance value will be 1 (white), and the others change smoothly.

void enhance(BITMAP *bmImg, BITMAP *bmEnhance) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    gen_color(bmImg, bmEnhance);
    double *bmYUV = (double *) malloc(bmImg->bmBytesPerRow * bmiHeader->biHeight << 3);
    double minD = 1, maxD = 0;
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmImg->bmBytesPerRow + w * bmImg->bmBytesPerPel;
            uint8_t *B = &bmImg->bmData[pos];
            uint8_t *G = &bmImg->bmData[pos + 1];
            uint8_t *R = &bmImg->bmData[pos + 2];
            // transform RGB into YUV
            bmYUV[pos] = 0.299 * *R + 0.587 * *G + 0.114 * *B;
            bmYUV[pos + 1] = -0.147 * *R - 0.289 * *G + 0.436 * *B;
            bmYUV[pos + 2] = 0.615 * *R - 0.515 * *G - 0.100 * *B;
            double D = log(bmYUV[pos] + 1) / log(256);
            if (D < minD) {
                minD = D;
            }
            if (D > maxD) {
                maxD = D;
            }
        }
    }
    if (minD < maxD) {
        for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
            for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
                uint32_t pos = h * bmImg->bmBytesPerRow + w * bmImg->bmBytesPerPel;
                double *Y = &bmYUV[pos];
                double *U = &bmYUV[pos + 1];
                double *V = &bmYUV[pos + 2];
                double D = log(*Y + 1) / log(256);
                // histogram stretching
                *Y = 255 * (D - minD) / (maxD - minD);
                // transform YUV into RGB
                bmEnhance->bmData[pos] = adjust(*Y + 2.032 * *U);
                bmEnhance->bmData[pos + 1] = adjust(*Y - 0.395 * *U - 0.581 * *V);
                bmEnhance->bmData[pos + 2] = adjust(*Y + 1.140 * *V);
            }
        }
    }
    free(bmYUV);
}

Step Four: Histogram Equalization

Histogram equalization is a method in image processing of contrast adjustment using the image’s histogram. Let $p_i$ be the probability of an occurrence of a pixel of level $i$ in the image, then its new level $s_i=\sum_{k=0}^i p_k$. The image should be converted to Lab color space or HSL/HSV color space before equalization, so that the algorithm can be applied to the luminance or value channel without resulting in changes to the hue and saturation of the image.

void histeq(BITMAP *bmImg, BITMAP *bmHistEq) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    gen_color(bmImg, bmHistEq);
    double *bmHSL = (double *) malloc(bmImg->bmBytesPerRow * bmiHeader->biHeight << 3);
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmImg->bmBytesPerRow + w * bmImg->bmBytesPerPel;
            double B = bmImg->bmData[pos] / 255.;
            double G = bmImg->bmData[pos + 1] / 255.;
            double R = bmImg->bmData[pos + 2] / 255.;
            double *H = &bmHSL[pos];
            double *S = &bmHSL[pos + 1];
            double *L = &bmHSL[pos + 2];
            // transform RGB into HSL
            double min = B < G ? B < R ? B : R : G < R ? G : R;
            double max = B > G ? B > R ? B : R : G > R ? G : R;
            if (min == max) {
                *H = 0;
            } else if (max == R) {
                *H = 60 * (G - B) / (max - min);
            } else if (max == G) {
                *H = 60 * (2 + (B - R) / (max - min));
            } else if (max == B) {
                *H = 60 * (4 + (R - G) / (max - min));
            }
            if (*H < 0) {
                *H += 360;
            }
            *L = (min + max) / 2;
            if (min == 1 || max == 0) {
                *S = 0;
            } else {
                *S = (max - *L) / (*L < 0.5 ? *L : 1 - *L);
            }
        }
    }
    // histogram equalization
    double *p = (double *) malloc(1 << 11);
    memset(p, 0, 1 << 11);
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmImg->bmBytesPerRow + w * bmImg->bmBytesPerPel;
            ++p[adjust(bmHSL[pos + 2] * 255)];
        }
    }
    for (int16_t i = 0; i < 256; ++i) {
        p[i] /= bmiHeader->biHeight * bmiHeader->biWidth;
        if (i > 0) {
            p[i] += p[i - 1];
        }
    }
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmImg->bmBytesPerRow + w * bmImg->bmBytesPerPel;
            double *H = &bmHSL[pos];
            double *S = &bmHSL[pos + 1];
            double *L = &bmHSL[pos + 2];
            *L = p[adjust(*L * 255)];
            // transform HSL into RGB
            double C = (1 - fabs(*L * 2 - 1)) * *S;
            double Hapos = *H / 60;
            double X = C * (1 - fabs(fmod(Hapos, 2) - 1));
            double R, G, B;
            if (Hapos <= 1) {
                R = C, G = X, B = 0;
            } else if (Hapos <= 2) {
                R = X, G = C, B = 0;
            } else if (Hapos <= 3) {
                R = 0, G = C, B = X;
            } else if (Hapos <= 4) {
                R = 0, G = X, B = C;
            } else if (Hapos <= 5) {
                R = X, G = 0, B = C;
            } else {
                R = C, G = 0, B = X;
            }
            double m = *L - C / 2;
            bmHistEq->bmData[pos] = adjust((B + m) * 255);
            bmHistEq->bmData[pos + 1] = adjust((G + m) * 255);
            bmHistEq->bmData[pos + 2] = adjust((R + m) * 255);
        }
    }
    free(p);
    free(bmHSL);
}