Digital Image Processing Assignment IV

Main Contents

  • Image translate;
  • Image mirror;
  • Image shear;
  • Image rotate;
  • Image scale.

Step One: Define Some Functions

We modify our gen_color function in order to generate a blank image with given scale.

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

Step Two: Translate, Mirror and Shear

Normally, we enumerate each pixel in the new image, map it back to the original image, and use the neighboring pixels to interpolate. However, for translating, mirroring and shearing, there is a one-to-one match between the pixels in two images, so interpolation is unnecessary.

void translate(BITMAP *bmImg, BITMAP *bmTranslate, uint32_t disX, uint32_t disY) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    uint32_t height = bmiHeader->biHeight + disX;
    uint32_t width = bmiHeader->biWidth + disY;
    gen_color(bmImg, bmTranslate, height, width);
    for (uint32_t h = 0; h < height; ++h) {
        for (uint32_t w = 0; w < width; ++w) {
            int32_t x = h;
            int32_t y = w - disY;
            uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
            uint32_t _pos = h * bmTranslate->bmBytesPerRow + w * bmTranslate->bmBytesPerPel;
            if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
                for (uint8_t k = 0; k < 3; ++k) {
                    bmTranslate->bmData[_pos + k] = bmImg->bmData[pos + k];
                }
            }
        }
    }
}

void mirror(BITMAP *bmImg, BITMAP *bmMirror) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    uint32_t height = bmiHeader->biHeight << 1;
    uint32_t width = bmiHeader->biWidth;
    gen_color(bmImg, bmMirror, height, width);
    for (uint32_t h = 0; h < height; ++h) {
        for (uint32_t w = 0; w < width; ++w) {
            int32_t x = bmiHeader->biHeight - 1 - h;
            int32_t y = w;
            uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
            uint32_t _pos = h * bmMirror->bmBytesPerRow + w * bmMirror->bmBytesPerPel;
            if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
                for (uint8_t k = 0; k < 3; ++k) {
                    bmMirror->bmData[_pos + k] = bmImg->bmData[pos + k];
                }
            }
        }
    }
}

void shear(BITMAP *bmImg, BITMAP *bmShear, uint32_t disX) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    uint32_t height = bmiHeader->biHeight + disX;
    uint32_t width = bmiHeader->biWidth;
    gen_color(bmImg, bmShear, height, width);
    for (uint32_t h = 0; h < height; ++h) {
        for (uint32_t w = 0; w < width; ++w) {
            int32_t x = h + disX * w / (width - 1) - disX;
            int32_t y = w;
            uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
            uint32_t _pos = h * bmShear->bmBytesPerRow + w * bmShear->bmBytesPerPel;
            if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
                for (uint8_t k = 0; k < 3; ++k) {
                    bmShear->bmData[_pos + k] = bmImg->bmData[pos + k];
                }
            }
        }
    }
}

Step Three: Interpolate

There are several methods of interpolation, such as nearest-neighbor interpolation, bilinear interpolation, and bicubic interpolation. Among these, nearest-neighbor interpolation has the highest speed, bicubic interpolation has the highest quality.

The bicubic interpolation can be summarized as solving a system of linear equations with $16$ variables, i.e., find $a_{ij}$, s.t. $p(x, y)=\sum_{i=0}^3 \sum_{j=0}^3 a_{ij} x^i y^j$. An interpolator with similar properties can be obtained by applying a convolution with the following kernel in both dimensions,

$$\begin{equation}
W(x)=\begin{cases}
(a+2)|x|^3-(a+3)|x|^2+1 & ,\ |x|\le 1 \\
a|x|^3-5a|x|^2+8a|x|-4a & ,\ 1<|x|\le 2 \\
0 & ,\ |x|>2
\end{cases}
\end{equation}$$

where $a$ is usually set to $-0.5$. At this time, the equation can be expressed in a more friendly manner,

$$\begin{equation}
p(t)=\frac{1}{2}
\left[\begin{matrix}
1 & t & t^2 & t^3
\end{matrix}\right]
\left[\begin{matrix}
0 & 2 & 0 & 0 \\
-1 & 0 & 1 & 0 \\
2 & -5 & 4 & -1 \\
-1 & 3 & -3 & 1
\end{matrix}\right]
\left[\begin{matrix}
f_{-1} \\ f_0 \\ f_1 \\ f_2
\end{matrix}\right]
\end{equation}$$

for $t$ between $0$ and $1$ for one dimension. Note that for $1$-dimensional cubic convolution interpolation $4$ sample points are required. For each inquiry two samples are located on its left and two samples on the right. These points are indexed from $−1$ to $2$ in this text. The distance from the point indexed with $0$ to the inquiry point is denoted by $t$ here. For two dimensions first applied once in $y$ and again in $x$.

double interpolate(double t1, double f0, double f1, double f2, double f3) {
    double t2 = t1 * t1, t3 = t2 * t1;
    double ret = (-t1 + 2 * t2 - t3) * f0;
    ret += (2 - 5 * t2 + 3 * t3) * f1;
    ret += (t1 + 4 * t2 - 3 * t3) * f2;
    ret += (-t2 + t3) * f3;
    return ret / 2;
}

Step Four: Rotate

First we need to calculate the size of the new image. Then we use nearest-neighbor interpolation and bicubic interpolation separately.

// nearest-neighbor interpolation
void simple_rotate(BITMAP *bmImg, BITMAP *bmRotate, double theta) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    double minX = fmin(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
    double maxX = fmax(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
    double minY = fmin(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
    double maxY = fmax(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
    minX = fmin(minX, fmin(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
    maxX = fmax(maxX, fmax(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
    minY = fmin(minY, fmin(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
    maxY = fmax(maxY, fmax(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
    uint32_t height = ceil(maxX) - floor(minX);
    uint32_t width = ceil(maxY) - floor(minY);
    gen_color(bmImg, bmRotate, height, width);
    for (uint32_t h = 0; h < height; ++h) {
        for (uint32_t w = 0; w < width; ++w) {
            int32_t x = round((h + minX) * cos(theta) + (w + minY) * sin(theta));
            int32_t y = round((h + minX) * -sin(theta) + (w + minY) * cos(theta));
            uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
            uint32_t _pos = h * bmRotate->bmBytesPerRow + w * bmRotate->bmBytesPerPel;
            if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
                for (uint8_t k = 0; k < 3; ++k) {
                    bmRotate->bmData[_pos + k] = bmImg->bmData[pos + k];
                }
            }
        }
    }
}

// bicubic interpolation
void rotate(BITMAP *bmImg, BITMAP *bmRotate, double theta) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    double minX = fmin(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
    double maxX = fmax(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
    double minY = fmin(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
    double maxY = fmax(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
    minX = fmin(minX, fmin(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
    maxX = fmax(maxX, fmax(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
    minY = fmin(minY, fmin(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
    maxY = fmax(maxY, fmax(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
    uint32_t height = ceil(maxX) - floor(minX);
    uint32_t width = ceil(maxY) - floor(minY);
    gen_color(bmImg, bmRotate, height, width);
    for (uint8_t k = 0; k < 3; ++k) {
        for (uint32_t h = 0; h < height; ++h) {
            for (uint32_t w = 0; w < width; ++w) {
                double x = (h + minX) * cos(theta) + (w + minY) * sin(theta);
                double y = (h + minX) * -sin(theta) + (w + minY) * cos(theta);
                uint32_t _pos = h * bmRotate->bmBytesPerRow + w * bmRotate->bmBytesPerPel;
                if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
                    double g[4];
                    for (int8_t u = -1; u < 3; ++u) {
                        double f[4];
                        int32_t _x = max(0, min(bmiHeader->biHeight - 1, (int32_t) x + u));
                        for (int8_t v = -1; v < 3; ++v) {
                            int32_t _y = max(0, min(bmiHeader->biWidth - 1, (int32_t) y + v));
                            f[v + 1] = bmImg->bmData[_x * bmImg->bmBytesPerRow + _y * bmImg->bmBytesPerPel + k];
                        }
                        // interpolate horizontally
                        g[u + 1] = interpolate(y - (int32_t) y, f[0], f[1], f[2], f[3]);
                    }
                    // interpolate vertically
                    bmRotate->bmData[_pos + k] = adjust(interpolate(x - (int32_t) x, g[0], g[1], g[2], g[3]));
                }
            }
        }
    }
}

Step Five: Scale

Similar to rotating, interpolation is vital for scaling. A slight difference is that, we can storage the results of all horizontal interpolations before performing the vertical, since they will be reused often.

// nearest-neighbor interpolation
void simple_scale(BITMAP *bmImg, BITMAP *bmScale, double ratioH, double ratioW) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    uint32_t height = ceil(ratioH * bmiHeader->biHeight);
    uint32_t width = ceil(ratioW * bmiHeader->biWidth);
    gen_color(bmImg, bmScale, height, width);
    for (uint32_t h = 0; h < height; ++h) {
        for (uint32_t w = 0; w < width; ++w) {
            int32_t x = round(h / ratioH);
            int32_t y = round(w / ratioW);
            x = max(0, min(bmiHeader->biHeight - 1, x));
            y = max(0, min(bmiHeader->biWidth - 1, y));
            uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
            uint32_t _pos = h * bmScale->bmBytesPerRow + w * bmScale->bmBytesPerPel;
            if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
                for (uint8_t k = 0; k < 3; ++k) {
                    bmScale->bmData[_pos + k] = bmImg->bmData[pos + k];
                }
            }
        }
    }
}

// bicubic interpolation
void scale(BITMAP *bmImg, BITMAP *bmScale, double ratioH, double ratioW) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    uint32_t height = ceil(ratioH * bmiHeader->biHeight);
    uint32_t width = ceil(ratioW * bmiHeader->biWidth);
    gen_color(bmImg, bmScale, height, width);
    double *p = (double *) malloc(bmiHeader->biHeight * width << 3);
    for (uint8_t k = 0; k < 3; ++k) {
        for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
            for (uint32_t w = 0; w < width; ++w) {
                double y = w / ratioW;
                double f[4];
                for (int8_t u = -1; u < 3; ++u) {
                    int32_t _y = max(0, min(bmiHeader->biWidth - 1, (int32_t) y + u));
                    f[u + 1] = bmImg->bmData[h * bmImg->bmBytesPerRow + _y * bmImg->bmBytesPerPel + k];
                }
                // interpolate horizontally
                p[h * width + w] = interpolate(y - (int32_t) y, f[0], f[1], f[2], f[3]);
            }
        }
        for (uint32_t h = 0; h < height; ++h) {
            for (uint32_t w = 0; w < width; ++w) {
                double x = h / ratioH;
                uint32_t _pos = h * bmScale->bmBytesPerRow + w * bmScale->bmBytesPerPel;
                double f[4];
                for (int8_t u = -1; u < 3; ++u) {
                    int32_t _x = max(0, min(bmiHeader->biHeight - 1, (int32_t) x + u));
                    f[u + 1] = p[_x * width + w];
                }
                // interpolate vertically
                bmScale->bmData[_pos + k] = adjust(interpolate(x - (int32_t) x, f[0], f[1], f[2], f[3]));
            }
        }
    }
    free(p);
}

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^{\frac{1}{\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=\frac{\log(L_w+1)}{\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);
}