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);
}

Digital Image Processing Assignment II

Main Contents

  • Image binarization;
  • Binary image erosion;
  • Binary image dilation;
  • Binary image opening;
  • Binary image closing.

Step One: Define Some Functions

Here we define two functions that generate a blank grayscale image from a given image,

// given bmHeader and bmInfo, initialize others
void init_bmp(BITMAP *bmImg) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    bmImg->bmBytesPerRow = ((bmiHeader->biWidth * bmiHeader->biBitCount + 31) >> 5) << 2;
    bmImg->bmBytesPerPel = bmiHeader->biBitCount >> 3;
    bmImg->bmData = (uint8_t *) malloc(bmImg->bmBytesPerRow * bmiHeader->biHeight);
    memset(bmImg->bmData, 255, bmImg->bmBytesPerRow * bmiHeader->biHeight);
}

// generate a blank grayscale image
void gen_gray(BITMAP *bmImg, BITMAP *bmGray) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    bmGray->bmHeader = bmImg->bmHeader;
    bmGray->bmInfoSize = sizeof(BITMAPINFOHEADER) + (sizeof(RGBQUAD) << 8); // add palette
    bmGray->bmInfo = (BITMAPINFO *) malloc(bmGray->bmInfoSize);
    bmGray->bmInfo->bmiHeader = *bmiHeader;
    bmGray->bmInfo->bmiHeader.biBitCount = 8; // one byte per pixel
    // initilize palette
    for (int16_t i = 0; i < 256; ++i) {
        RGBQUAD *rgb = &bmGray->bmInfo->bmiColors[i];
        rgb->rgbBlue = i;
        rgb->rgbGreen = i;
        rgb->rgbRed = i;
        rgb->rgbReserved = 0;
    }
    init_bmp(bmGray);
    // IMPORTANT: MODIFY byOffBits AND bfSize
    bmGray->bmHeader.bfOffBits = sizeof(BITMAPFILEHEADER) + bmGray->bmInfoSize;
    bmGray->bmHeader.bfSize = bmGray->bmHeader.bfOffBits + bmGray->bmBytesPerRow * bmiHeader->biHeight;
}

and a function that change a color image into grayscale.

void color2gray(BITMAP *bmImg, BITMAP *bmGray) {
    BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
    gen_gray(bmImg, bmGray);
    uint8_t minY = 255, maxY = 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];
            uint8_t Y = adjust(0.299 * *R + 0.587 * *G + 0.114 * *B);
            if (Y < minY) {
                minY = Y;
            }
            if (Y > maxY) {
                maxY = Y;
            }
        }
    }
    // rearrange gray indensity
    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];
            uint8_t Y = adjust(0.299 * *R + 0.587 * *G + 0.114 * *B);
            uint32_t _pos = h * bmGray->bmBytesPerRow + w;
            bmGray->bmData[_pos] = adjust(255. * (Y - minY) / (maxY - minY));
        }
    }
}

Step Two: Change Gray to Binary

In order to binarize the image, we need determine a threshold, and change pixels whose grayscale is below the threshold to black, and the others to white. But how to find the optimal threshold? There’s an excellent algorithm called Otsu’s method, which, in brief, maximize inter-class variance.

void otsu_gray2binary(BITMAP *bmGray, BITMAP *bmBinary) {
    BITMAPINFOHEADER *bmiHeader = &bmGray->bmInfo->bmiHeader;
    gen_gray(bmGray, bmBinary);
    uint8_t minG = 255, maxG = 0;
    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 * bmGray->bmBytesPerRow + w;
            ++p[bmGray->bmData[pos]];
            if (bmGray->bmData[pos] < minG) {
                minG = bmGray->bmData[pos];
            }
            if (bmGray->bmData[pos] > maxG) {
                maxG = bmGray->bmData[pos];
            }
        }
    }
    double muT = 0;
    for (int16_t k = minG; k <= maxG; ++k) {
        p[k] /= bmiHeader->biHeight * bmiHeader->biWidth;
        muT += k * p[k];
    }
    uint8_t threshold = 0;
    double omegaK = 0, muK = 0, maxB = 0;
    for (int16_t k = minG; k < maxG; ++k) {
        omegaK += p[k];
        muK += k * p[k];
        double sigmaK = (muT * omegaK - muK) * (muT * omegaK - muK) / omegaK / (1 - omegaK);
        if (maxB < sigmaK) {
            maxB = sigmaK;
            threshold = k;
        }
    }
    free(p);
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmGray->bmBytesPerRow + w;
            bmBinary->bmData[pos] = bmGray->bmData[pos] < threshold ? 0 : 255;
        }
    }
}

Otsu’s method sets the global threshold, but for some special image with different luminance, adaptive thresholding may produce better result. There are at least two ways to apply adaptive thresholding.

void pixel_adaptive_gray2binary(BITMAP *bmGray, BITMAP *bmBinary, uint32_t border, uint32_t offset) {
    BITMAPINFOHEADER *bmiHeader = &bmGray->bmInfo->bmiHeader;
    gen_gray(bmGray, bmBinary);
    border >>= 1;
    uint32_t *integral = (uint32_t *) malloc(bmGray->bmBytesPerRow * bmiHeader->biHeight << 2);
    memset(integral, 0, bmGray->bmBytesPerRow * bmiHeader->biHeight * 4);
    // get integral image
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos1 = h * bmGray->bmBytesPerRow + w;
            uint32_t pos2 = (h - 1) * bmGray->bmBytesPerRow + w;
            uint32_t pos3 = h * bmGray->bmBytesPerRow + w - 1;
            uint32_t pos4 = (h - 1) * bmGray->bmBytesPerRow + w - 1;
            integral[pos1] = bmGray->bmData[pos1]
                           + (h > 0 ? integral[pos2] : 0)
                           + (w > 0 ? integral[pos3] : 0)
                           - (h > 0 && w > 0 ? integral[pos4] : 0);
        }
    }
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmGray->bmBytesPerRow + w;
            // range checking
            uint32_t x1 = h > border ? h - border : 0;
            uint32_t x2 = h + border + 1 < bmiHeader->biHeight ? h + border + 1 : bmiHeader->biHeight;
            uint32_t y1 = w > border ? w - border : 0;
            uint32_t y2 = w + border + 1 < bmiHeader->biWidth ? w + border + 1 : bmiHeader->biWidth;
            uint32_t pos1 = (x2 - 1) * bmGray->bmBytesPerRow + y2 - 1;
            uint32_t pos2 = (x1 - 1) * bmGray->bmBytesPerRow + y2 - 1;
            uint32_t pos3 = (x2 - 1) * bmGray->bmBytesPerRow + y1 - 1;
            uint32_t pos4 = (x1 - 1) * bmGray->bmBytesPerRow + y1 - 1;
            uint64_t cnt = (x2 - x1) * (y2 - y1);
            uint64_t sum = integral[pos1]
                         - (x1 > 0 ? integral[pos2] : 0)
                         - (y1 > 0 ? integral[pos3] : 0)
                         + (x1 > 0 && y1 > 0 ? integral[pos4] : 0);
            bmBinary->bmData[pos] = bmGray->bmData[pos] * cnt * 100 < sum * (100 - offset) ? 0 : 255;
        }
    }
}

void block_adaptive_gray2binary(BITMAP *bmGray, BITMAP *bmBinary, uint32_t border, double step, uint32_t offset) {
    BITMAPINFOHEADER *bmiHeader = &bmGray->bmInfo->bmiHeader;
    gen_gray(bmGray, bmBinary);
    uint32_t *integral = (uint32_t *) malloc(bmGray->bmBytesPerRow * bmiHeader->biHeight << 2);
    memset(integral, 0, bmGray->bmBytesPerRow * bmiHeader->biHeight * 4);
    // get integral image
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos1 = h * bmGray->bmBytesPerRow + w;
            uint32_t pos2 = (h - 1) * bmGray->bmBytesPerRow + w;
            uint32_t pos3 = h * bmGray->bmBytesPerRow + w - 1;
            uint32_t pos4 = (h - 1) * bmGray->bmBytesPerRow + w - 1;
            integral[pos1] = bmGray->bmData[pos1]
                           + (h > 0 ? integral[pos2] : 0)
                           + (w > 0 ? integral[pos3] : 0)
                           - (h > 0 && w > 0 ? integral[pos4] : 0);
        }
    }
    for (uint32_t h = 0; h < bmiHeader->biHeight; h += border * step) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; w += border * step) {
            uint32_t pos = h * bmGray->bmBytesPerRow + w;
            // range checking
            uint32_t x1 = h;
            uint32_t x2 = h + border + 1 < bmiHeader->biHeight ? h + border + 1 : bmiHeader->biHeight;
            uint32_t y1 = w;
            uint32_t y2 = w + border + 1 < bmiHeader->biWidth ? w + border + 1 : bmiHeader->biWidth;
            uint32_t pos1 = (x2 - 1) * bmGray->bmBytesPerRow + y2 - 1;
            uint32_t pos2 = (x1 - 1) * bmGray->bmBytesPerRow + y2 - 1;
            uint32_t pos3 = (x2 - 1) * bmGray->bmBytesPerRow + y1 - 1;
            uint32_t pos4 = (x1 - 1) * bmGray->bmBytesPerRow + y1 - 1;
            uint64_t cnt = (x2 - x1) * (y2 - y1);
            uint64_t sum = integral[pos1]
                         - (x1 > 0 ? integral[pos2] : 0)
                         - (y1 > 0 ? integral[pos3] : 0)
                         + (x1 > 0 && y1 > 0 ? integral[pos4] : 0);
            for (uint32_t _h = x1; _h < x2; ++_h) {
                for (uint32_t _w = y1; _w < y2; ++_w) {
                    uint32_t _pos = _h * bmGray->bmBytesPerRow + _w;
                    // when overlapping, black has the priority
                    bmBinary->bmData[_pos] &= bmGray->bmData[_pos] * cnt * 100 < sum * (100 - offset) ? 0 : 255;
                }
            }
        }
    }
}

Step Three: Erosion and Dilation

Let’s say that white is the foreground and black is the background. Assume you have a solid circle, and you can only place it on white pixels, which means the circle shouldn’t cover the black pixels. We make the available area where the center can be put white and the others black, and get a new image—the result of erosion.

void erode(BITMAP *bmBinary, BITMAP *bmErosion, uint32_t border) {
    BITMAPINFOHEADER *bmiHeader = &bmBinary->bmInfo->bmiHeader;
    gen_gray(bmBinary, bmErosion);
    border >>= 1;
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmBinary->bmBytesPerRow + w;
            // range checking
            uint32_t x1 = h > border ? h - border : 0;
            uint32_t x2 = h + border + 1 < bmiHeader->biHeight ? h + border + 1 : bmiHeader->biHeight;
            uint32_t y1 = w > border ? w - border : 0;
            uint32_t y2 = w + border + 1 < bmiHeader->biWidth ? w + border + 1 : bmiHeader->biWidth;
            uint8_t flag = 255;
            for (uint32_t _h = x1; _h < x2; ++_h) {
                for (uint32_t _w = y1; _w < y2; ++_w) {
                    uint32_t _pos = _h * bmBinary->bmBytesPerRow + _w;
                    flag &= bmBinary->bmData[_pos];
                }
            }
            bmErosion->bmData[pos] = flag;
        }
    }
}

This time, you want to put the center of the circle on white pixels regardless of whether other points lie on black pixels. We make the area that can be covered by the circle white and the others black. This is what dilation do. Surprisingly, the code is very very similar with erosion.

void dilate(BITMAP *bmBinary, BITMAP *bmDilation, uint32_t border) {
    BITMAPINFOHEADER *bmiHeader = &bmBinary->bmInfo->bmiHeader;
    gen_gray(bmBinary, bmDilation);
    border >>= 1;
    for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
        for (uint32_t w = 0; w < bmiHeader->biWidth; ++w) {
            uint32_t pos = h * bmBinary->bmBytesPerRow + w;
            // range checking
            uint32_t x1 = h > border ? h - border : 0;
            uint32_t x2 = h + border + 1 < bmiHeader->biHeight ? h + border + 1 : bmiHeader->biHeight;
            uint32_t y1 = w > border ? w - border : 0;
            uint32_t y2 = w + border + 1 < bmiHeader->biWidth ? w + border + 1 : bmiHeader->biWidth;
            uint8_t flag = 0;
            for (uint32_t _h = x1; _h < x2; ++_h) {
                for (uint32_t _w = y1; _w < y2; ++_w) {
                    uint32_t _pos = _h * bmBinary->bmBytesPerRow + _w;
                    flag |= bmBinary->bmData[_pos];
                }
            }
            bmDilation->bmData[pos] = flag;
        }
    }
}

Step Four: Opening and Closing

In fact, opening operation and closing operation are just the combination of the two operations we mentioned above, the only difference is the order. Opening operation is an erosion followed by a dilation, while closing operation is the reverse.

Visually speaking, opening operation makes the gaps between black pixels disappeared, and closing operation makes the gaps between white pixels disappeared.

Note that these two operations are both idempotent, which means performing one of them multiple times is equivalent to performing it one time.

Give Up Sometimes

The expression “Never, never give up” means to keep trying and never stop working for your goals. Do you agree or disagree with this statement? Use specific reasons and examples to support your answer.

Topic

Many people believe in the saying that we should never give up casually, which tells them to work incessantly for their goals. Admittedly, this spirit-lifting saying has given countless people a pat on the back, inspired them to strive for their long-held dreams and acquire prosperity in the end. But generally speaking, I suppose we don’t have to shy away from giving up. Sometimes we need to put aside our persistence and accept the abnegation.

First of all, the initial orientation matters. We can’t guarantee that the first step we take is towards the right direction, which means we can be heading the wrong without perceiving it. If so, it would be better to just turn around than to press on. There’s no need to sacrifice your success for simply senseless insistence. One example of this is an ancient idiom from China, narrating a man who intended to go south steered his carriage to the north regardless of the prompting from a passer-by. He had not realized the sheer meaninglessness of his action. Of course, he has not arrived his destination yet. Had he listened to others’ advice and turned back, he would have reached it.

Besides, we should pay heed to the progress. Not only are we likely to move forward to where we do not expect at the beginning, but also during the process, and thus adjusting ourselves dynamically should be attached significance to. In addition, even though we know clearly what we are doing and why we are doing it, we may lack the capacity to do it well, and the progress can fall short of our expectation. When we’re patently unsuitable for tackling certain type of things, insisting may not be the optimal option. Some classmates in my high school will be rather representative instances. As soon as they entered the school, each of them picked an Olympiad in a certain subject. Yet before long, they found themselves not very good at competition, thereby opting out, focusing on regular curriculum. It turned out ultimately that they obtained fairly high scores in the college entrance examination. We can derive from it the fact that giving up does not necessarily mean losing, but somehow means wining.

What’s more, the outcome may be perplexing. When it come out, it’s time to judge whether it is right or not. We can scarcely ensure that the result will always fulfil our anticipation, so when we receive an unexpected result which is against the basic facts, probably it should not be calculated in statistics. At this time, alleging that the result has no problem seems a bit stubborn. On the other hand, if the result does not correspond with what we used to deem it should be, the chances are what deceives you is not the result, but the acknowledgement. If we look back upon history of science and technology, there are innumerable examples, such as Galileo who doubted the previous law of free fall, Bohr who suspected the previous nuclear model, etc. All of these indicate that even the most seeming truth has the possibility to be broken someday. The interesting thing is that, the result and the truth were dropped alternately, which set the pace for the development of human civilization. We wouldn’t have today’s advancement if there were no giving up in the world.

The old saying “Never, never give up” will continue to exert its positive attitude and lift people up. However, I take the viewpoint that people sometimes need to give up what they think, have or believe. This is not persuading you to give up everything in your life, but encouraging you to give up in a more sensible way, in which you will better encounter what you love, accomplish what you dream, and make what you desire to be.

Digital Image Processing Assignment I

Main Contents

  • Read a color bmp;
  • RGB->YUV;
  • Color to gray: gray=Y in YUV color space;
  • Rearrange gray intensity to lie between [0,255];
  • Write a grayscale bmp;
  • Change the luminance value Y;
  • YUV->RGB;
  • Write a color bmp.

Step One: BMP File Structure

A common BMP file is comprised of four part: image file header, image information header, palette and image data.

The image file header is a struct whose length is 14 bytes. Here gives its definition,

typedef struct tagBITMAPFILEHEADER {
    WORD bfType;  
    DWORD bfSize;  
    WORD bfReserved1;  
    WORD bfReserved2; 
    DWORD bfOffBits; 
} BITMAPFILEHEADER;

and the explanation for every variable.

  • bfType: must always be set to ‘BM’ to declare that this is a .bmp-file;
  • bfSize: specifies the size of the file in bytes;
  • bfReserved1: must always be set to zero;
  • bfReserved2: must always be set to zero;
  • bfOffBits: specifies the offset from the beginning of the file to the bitmap data.

The image information header is also a struct, while its length is 40 bytes. The definition

typedef struct tagBITMAPINFOHEADER {
    DWORD biSize;     
    LONG biWidth;    
    LONG biHeight;    
    WORD biPlanes;    
    WORD biBitCount   
    DWORD biCompression;  
    DWORD biSizeImage;   
    LONG biXPelsPerMeter;  
    LONG biYPelsPerMeter;  
    DWORD biClrUsed;  
    DWORD biClrImportant; 
} BITMAPINFOHEADER;

The explanation

  • biSize: number of bytes to define BITMAPINFORHEADER structure;
  • biWidth: image width (number of pixels);
  • biHeight: image height (number of pixels), note that if it’s a positive number, the image is inverted, otherwise upright;
  • biPlanes: number of planes, should always be 1;
  • biBitCount: bits per pixel, which may be 1, 4, 8, 16, 24, 32;
  • biCompression: compression type, only non-compression(BI_RGB) is discussed here;
  • biSizeImage: image size with bytes, when biCompression is BI_RGB, biSizeImage is 0;
  • biXPelsPerMeter: horizontal resolution, pixels per meter;
  • biYPelsPerMeter: vertical resolution, pixels per meter;
  • biClrUsed: number of color indices used in the bitmap, when it’s 0, all the palette items are used;
  • biClrImportant: number of important color indices for image display, when it’s 0, all items are important.

The palette has a series of RGBQUADs, which is defined like this.

typedef struct tagRGBQUAD {
    uint8_t rgbBlue;
    uint8_t rgbGreen;
    uint8_t rgbRed;
    uint8_t rgbReserved;
} RGBQUAD;

Note that the order of the color is blue, green, and red, not the reverse. The number of RGBQUADs is decided by biBitCount and biClrUsed.

Next we need to define BITMAPINFO.

typedef struct tagBITMAPINFO {
    BITMAPINFOHEADER bmiHeader;
    RGBQUAD bmiColors[1];
} BITMAPINFO;

As we know nothing about the number of RGBQUADs an image uses, the BITMAPINFO should be defined as a pointer so that right amount of memory can be allocated to it.

The image data contains color of all pixels, and every biBitCount bit(s) represents a pixel.

At last, we define a struct to storage a full BMP image.

typedef struct tagBITMAP {
    BITMAPFILEHEADER bmHeader;
    BITMAPINFO *bmInfo;
    uint32_t bmInfoSize;
    uint32_t bmBytesPerRow;
    uint8_t bmBytesPerPel;
    uint8_t *bmData;
} BITMAP;

When defining the two structs, we need to add a line #pragma pack(push, 1) to avoid struct padding.

Step Two: Read/Write a BMP File

A BMP file is a binary file, so we need to add "b" to the second parameter when using fopen.

Another important thing is that the number of bytes in one row must always be adjusted to fit into the border of a multiple of four, and we need to calculate how many bytes are there in one row.

For convenience, we define an initialize function, which receives bmHeader and bmInfo, and initializes other variables.

// given bmHeader and bmInfo, initialize others
void init_bmp(BITMAP *bmImg) {
    BITMAPINFOHEADER *bmiHeader = &(bmImg->bmInfo->bmiHeader);
    bmImg->bmBytesPerRow = ((bmiHeader->biWidth * bmiHeader->biBitCount + 31) >> 5) << 2;
    bmImg->bmBytesPerPel = bmiHeader->biBitCount >> 3;
    bmImg->bmData = (uint8_t *) malloc(bmImg->bmBytesPerRow * bmiHeader->biHeight);
}

The read function is showed below.

// read a BMP from file
void read_bmp(BITMAP *bmImg, char *filepath) {
    FILE *fiInImg = fopen(filepath, "rb");
    BITMAPINFOHEADER bmiHeader;
    fread(&(bmImg->bmHeader), sizeof(BITMAPFILEHEADER), 1, fiInImg);
    fread(&bmiHeader, sizeof(BITMAPINFOHEADER), 1, fiInImg);
    // if biBitCount is less than 16, use all the palette, otherwise do not use palette
    if (bmiHeader.biBitCount < 16) {
        bmImg->bmInfoSize = sizeof(BITMAPINFOHEADER) + (sizeof(RGBQUAD) << bmiHeader.biBitCount);
    } else {
        bmImg->bmInfoSize = sizeof(BITMAPINFOHEADER);
    }
    bmImg->bmInfo = (BITMAPINFO *) malloc(bmImg->bmInfoSize);
    bmImg->bmInfo->bmiHeader = bmiHeader;
    if (bmiHeader.biBitCount < 16) {
        fread(bmImg->bmInfo->bmiColors, sizeof(RGBQUAD), 1 << bmiHeader.biBitCount, fiInImg);
    }
    init_bmp(bmImg);
    fread(bmImg->bmData, bmImg->bmBytesPerRow, bmiHeader.biHeight, fiInImg);
    fclose(fiInImg);
}

The write function looks similar with the read function.

// write a BMP to file
void write_bmp(BITMAP *bmImg, char *filepath) {
    FILE *fiOutImg = fopen(filepath, "wb");
    fwrite(&(bmImg->bmHeader), sizeof(BITMAPFILEHEADER), 1, fiOutImg);
    fwrite(bmImg->bmInfo, bmImg->bmInfoSize, 1, fiOutImg);
    fwrite(bmImg->bmData, bmImg->bmBytesPerRow, bmImg->bmInfo->bmiHeader.biHeight, fiOutImg);
    fclose(fiOutImg);
}

Step Three: Change Color to Gray

To make things easier, we define a function to duplicate a BMP file,

// duplicate a BMP
void copy_bmp(BITMAP *bmDes, BITMAP *bmSrc) {
    memcpy(bmDes, bmSrc, sizeof(BITMAP));
    bmDes->bmInfo = (BITMAPINFO *) malloc(bmSrc->bmInfoSize);
    memcpy(bmDes->bmInfo, bmSrc->bmInfo, bmSrc->bmInfoSize);
    BITMAPINFOHEADER *bmiHeader = &(bmSrc->bmInfo->bmiHeader);
    bmDes->bmData = (uint8_t *) malloc(bmSrc->bmBytesPerRow * bmiHeader->biHeight);
    memcpy(bmDes->bmData, bmSrc->bmData, bmSrc->bmBytesPerRow * bmiHeader->biHeight);
}

and a function to make sure the RGB value of a pixel lie between [0, 255].

// make RGB value legal
uint8_t adjust(double val) {
    int16_t ret = (int16_t) (val + 0.5);
    return ret < 0 ? 0 : ret > 255 ? 255 : ret;
}

We need to change RGB to YUV first, as the grayscale is determined by Y value. The formula is

$$\begin{bmatrix} 0.299 & 0.587 & 0.114 \\ -0.147 & -0.289 & 0.436 \\ 0.615 & -0.515 & -0.100 \end{bmatrix} \times \begin{bmatrix} R \\ G \\ B \end{bmatrix} = \begin{bmatrix} Y \\ U \\ V \end{bmatrix}$$

    // calculate YUV value
    double *bmYUV = (double *) malloc(sizeof(double) * bmImg.bmBytesPerRow * bmiHeader->biHeight);
    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];
            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;
        }
    }

To change color to gray, we can simply make the R, G and B value of a pixel equal to its Y value. It would be more complex if we want to change it into an image with biBitCount equal to 8, because we need to set the palette manually.

One more step, rearrange gray intensity to lie between [0,255]. It’s just a math problem. So the code

    // color to gray
    BITMAP bmGray;
    bmGray.bmHeader = bmImg.bmHeader;
    bmGray.bmInfoSize = sizeof(BITMAPINFOHEADER) + (sizeof(RGBQUAD) << 8);
    bmGray.bmInfo = (BITMAPINFO *) malloc(bmGray.bmInfoSize);
    bmGray.bmInfo->bmiHeader = *bmiHeader;
    bmGray.bmInfo->bmiHeader.biBitCount = 8;
    for (int i = 0; i < 256; ++i) {
        RGBQUAD *rgb = &(bmGray.bmInfo->bmiColors[i]);
        rgb->rgbBlue = i;
        rgb->rgbGreen = i;
        rgb->rgbRed = i;
        rgb->rgbReserved = 0;
    }
    init_bmp(&bmGray);
    uint8_t min = 255, max = 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;
            double *Y = &bmYUV[pos];
            if (*Y < min) {
                min = *Y;
            }
            if (*Y > max) {
                max = *Y;
            }
        }
    }
    // rearrange gray indensity
    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];
            uint32_t _pos = h * bmGray.bmBytesPerRow + w * bmGray.bmBytesPerPel;
            bmGray.bmData[_pos] = adjust(255 * (*Y - min) / (max - min));
        }
    }

Step Four: Change the Luminance

The luminance is depend on Y value, too. What we need to do is just changing the Y value and applying the inverse formula below.

$$\begin{bmatrix} 1.000 & 0.000 & 1.140 \\ 1.000 & -0.3946 & -0.5805 \\ 1.000 & 2.032 & -0.0005 \end{bmatrix} \times \begin{bmatrix} Y \\ U \\ V \end{bmatrix} = \begin{bmatrix} R \\ G \\ B \end{bmatrix}$$

And the code is simple.

    // change luminance
    BITMAP bmLight, bmDark;
    copy_bmp(&bmLight, &bmImg);
    copy_bmp(&bmDark, &bmImg);
    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];
            bmLight.bmData[pos] = adjust(*Y + 25 + 2.032 * *U - 0.0005 * *V);
            bmLight.bmData[pos + 1] = adjust(*Y + 25 - 0.3946 * *U - 0.5805 * *V);
            bmLight.bmData[pos + 2] = adjust(*Y + 25 + 1.140 * *V);
            bmDark.bmData[pos] = adjust(*Y - 50 + 2.032 * *U - 0.0005 * *V);
            bmDark.bmData[pos + 1] = adjust(*Y - 50 - 0.3946 * *U - 0.5805 * *V);
            bmDark.bmData[pos + 2] = adjust(*Y - 50 + 1.140 * *V);
        }
    }

Step Five: Complete Source Code

bmp.h

#ifndef _BMP_H_
#define _BMP_H_

#include <stdint.h>

#pragma pack(push, 1) // avoid struct padding

typedef struct tagBITMAPFILEHEADER {
    uint16_t bfType;
    uint32_t bfSize;
    uint16_t bfReserved1;
    uint16_t bfReserved2;
    uint32_t bfOffBits;
} BITMAPFILEHEADER;

typedef struct tagBITMAPINFOHEADER {
    uint32_t biSize;
    int32_t biWidth;
    int32_t biHeight;
    uint16_t biPlanes;
    uint16_t biBitCount;
    uint32_t biCompression;
    uint32_t biSizeImage;
    int32_t biXPelsPerMeter;
    int32_t biYPelsPerMeter;
    uint32_t biClrUsed;
    uint32_t biClrImportant;
} BITMAPINFOHEADER;

typedef struct tagRGBQUAD {
    uint8_t rgbBlue;
    uint8_t rgbGreen;
    uint8_t rgbRed;
    uint8_t rgbReserved;
} RGBQUAD;

typedef struct tagBITMAPINFO {
    BITMAPINFOHEADER bmiHeader;
    RGBQUAD bmiColors[1];
} BITMAPINFO;

typedef struct tagBITMAP {
    BITMAPFILEHEADER bmHeader;
    BITMAPINFO *bmInfo;
    uint32_t bmInfoSize;
    uint32_t bmBytesPerRow;
    uint8_t bmBytesPerPel;
    uint8_t *bmData;
} BITMAP;

#pragma pack(pop)

#endif

bmp.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "bmp.h"

// given bmHeader and bmInfo, initialize others
void init_bmp(BITMAP *bmImg) {
    BITMAPINFOHEADER *bmiHeader = &(bmImg->bmInfo->bmiHeader);
    bmImg->bmBytesPerRow = ((bmiHeader->biWidth * bmiHeader->biBitCount + 31) >> 5) << 2;
    bmImg->bmBytesPerPel = bmiHeader->biBitCount >> 3;
    bmImg->bmData = (uint8_t *) malloc(bmImg->bmBytesPerRow * bmiHeader->biHeight);
}

// read a BMP from file
void read_bmp(BITMAP *bmImg, char *filepath) {
    FILE *fiInImg = fopen(filepath, "rb");
    BITMAPINFOHEADER bmiHeader;
    fread(&(bmImg->bmHeader), sizeof(BITMAPFILEHEADER), 1, fiInImg);
    fread(&bmiHeader, sizeof(BITMAPINFOHEADER), 1, fiInImg);
    // if biBitCount is less than 16, use all the palette, otherwise do not use palette
    if (bmiHeader.biBitCount < 16) {
        bmImg->bmInfoSize = sizeof(BITMAPINFOHEADER) + (sizeof(RGBQUAD) << bmiHeader.biBitCount);
    } else {
        bmImg->bmInfoSize = sizeof(BITMAPINFOHEADER);
    }
    bmImg->bmInfo = (BITMAPINFO *) malloc(bmImg->bmInfoSize);
    bmImg->bmInfo->bmiHeader = bmiHeader;
    if (bmiHeader.biBitCount < 16) {
        fread(bmImg->bmInfo->bmiColors, sizeof(RGBQUAD), 1 << bmiHeader.biBitCount, fiInImg);
    }
    init_bmp(bmImg);
    fread(bmImg->bmData, bmImg->bmBytesPerRow, bmiHeader.biHeight, fiInImg);
    fclose(fiInImg);
}

// duplicate a BMP
void copy_bmp(BITMAP *bmDes, BITMAP *bmSrc) {
    memcpy(bmDes, bmSrc, sizeof(BITMAP));
    bmDes->bmInfo = (BITMAPINFO *) malloc(bmSrc->bmInfoSize);
    memcpy(bmDes->bmInfo, bmSrc->bmInfo, bmSrc->bmInfoSize);
    BITMAPINFOHEADER *bmiHeader = &(bmSrc->bmInfo->bmiHeader);
    bmDes->bmData = (uint8_t *) malloc(bmSrc->bmBytesPerRow * bmiHeader->biHeight);
    memcpy(bmDes->bmData, bmSrc->bmData, bmSrc->bmBytesPerRow * bmiHeader->biHeight);
}

// write a BMP to file
void write_bmp(BITMAP *bmImg, char *filepath) {
    FILE *fiOutImg = fopen(filepath, "wb");
    fwrite(&(bmImg->bmHeader), sizeof(BITMAPFILEHEADER), 1, fiOutImg);
    fwrite(bmImg->bmInfo, bmImg->bmInfoSize, 1, fiOutImg);
    fwrite(bmImg->bmData, bmImg->bmBytesPerRow, bmImg->bmInfo->bmiHeader.biHeight, fiOutImg);
    fclose(fiOutImg);
}

// make RGB value legal
uint8_t adjust(double val) {
    int16_t ret = (int16_t) (val + 0.5);
    return ret < 0 ? 0 : ret > 255 ? 255 : ret;
}

int main(int argc, char *argv[]) {
    // read bmp file
    BITMAP bmImg;
    read_bmp(&bmImg, "original.bmp");
    BITMAPINFOHEADER *bmiHeader = &(bmImg.bmInfo->bmiHeader);

    // calculate YUV value
    double *bmYUV = (double *) malloc(sizeof(double) * bmImg.bmBytesPerRow * bmiHeader->biHeight);
    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];
            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;
        }
    }

    // color to gray
    BITMAP bmGray;
    bmGray.bmHeader = bmImg.bmHeader;
    bmGray.bmInfoSize = sizeof(BITMAPINFOHEADER) + (sizeof(RGBQUAD) << 8);
    bmGray.bmInfo = (BITMAPINFO *) malloc(bmGray.bmInfoSize);
    bmGray.bmInfo->bmiHeader = *bmiHeader;
    bmGray.bmInfo->bmiHeader.biBitCount = 8;
    for (int i = 0; i < 256; ++i) {
        RGBQUAD *rgb = &(bmGray.bmInfo->bmiColors[i]);
        // to prove that the palette works well, play a small trick
        rgb->rgbBlue = (i >> 4) << 4;
        rgb->rgbGreen = (i >> 4) << 4;
        rgb->rgbRed = (i >> 4) << 4;
        rgb->rgbReserved = 0;
    }
    init_bmp(&bmGray);
    uint8_t min = 255, max = 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;
            double *Y = &bmYUV[pos];
            if (*Y < min) {
                min = *Y;
            }
            if (*Y > max) {
                max = *Y;
            }
        }
    }
    // rearrange gray indensity
    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];
            uint32_t _pos = h * bmGray.bmBytesPerRow + w * bmGray.bmBytesPerPel;
            bmGray.bmData[_pos] = adjust(255 * (*Y - min) / (max - min));
        }
    }
    write_bmp(&bmGray, "gray.bmp");

    // change luminance
    BITMAP bmLight, bmDark;
    copy_bmp(&bmLight, &bmImg);
    copy_bmp(&bmDark, &bmImg);
    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];
            bmLight.bmData[pos] = adjust(*Y + 25 + 2.032 * *U - 0.0005 * *V);
            bmLight.bmData[pos + 1] = adjust(*Y + 25 - 0.3946 * *U - 0.5805 * *V);
            bmLight.bmData[pos + 2] = adjust(*Y + 25 + 1.140 * *V);
            bmDark.bmData[pos] = adjust(*Y - 50 + 2.032 * *U - 0.0005 * *V);
            bmDark.bmData[pos + 1] = adjust(*Y - 50 - 0.3946 * *U - 0.5805 * *V);
            bmDark.bmData[pos + 2] = adjust(*Y - 50 + 1.140 * *V);
        }
    }
    write_bmp(&bmLight, "light.bmp");
    write_bmp(&bmDark, "dark.bmp");

    return 0;
}

Primes

题目描述

This is an interactive problem.
For two positive integers $x, y$, we define $\pi(x, y)$ to be the number of distinct primes that divide both $x$ and $y$. For example $\pi(2, 3) = 0, \; \pi(8, 16) = 1$ and $\pi(30, 105) = 2$.
For two positive integers $a, b$, where $a \le b$, we define $S(a, b)$ to be the sum of values $\pi(x, y)$ over all pairs of integers $(x, y)$ satisfying $a \le x \lt y \le b$.
Your task is to compute the values $S(a, b)$ for many query pairs $(a, b)$. To make your task more challenging, all the queries have to be answered online.

题意概述

给定两个整数$a, b$,定义$\pi(x, y)$为所有整除$x$和$y$的质数的个数,$S(a, b)$为所有满足$a \le x \lt y \le b$的$\pi(x, y)$的和,求$S(a, b)$。有$q$组询问,强制在线。

数据范围:$1 \le q \le 5 \times 10^4, \; 1 \le a \le b \le 10^6$。

算法分析

分别考虑每个质数对答案的贡献。若在区间$[a, b]$中有$k$个数能被质数$p$整除,则$p$对答案的贡献为${k \choose 2}$。

首先筛出所有质数。但是每次询问直接枚举质数会超时。考虑对于一个整数$n$,$\lfloor {n \over i} \rfloor$只有$O(\sqrt{n})$种不同的取值。因此可以对$\lfloor {a \over i} \rfloor$和$\lfloor {b \over i} \rfloor$进行分段枚举。

代码

#include <cstdio>
#include <cstring>
#include <algorithm>

int const N = 1000005;

int tp, prime[N], rec[N], vis[N], sum[N];

void init() {
    for (int i = 2; i < N; ++i) {
        if (!vis[i]) {
            prime[tp++] = i;
        }
        for (int j = 0; j < tp && i * prime[j] < N; ++j) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                break;
            }
        }
    }
    for (int i = 2; i < N; ++i)
        sum[i] = sum[i - 1] + !vis[i];
}

int main() {
    init();
    int q;
    scanf("%d", &q);
    for (; q--;) {
        int a, b;
        scanf("%d%d", &a, &b);
        long long ans = 0;
        for (int i = 0; i < tp && prime[i] < b;) {
            int cnt = b / prime[i] - (a - 1) / prime[i];
            int nxt = 1e9;
            if (b / prime[i]) nxt = std::min(nxt, b / (b / prime[i]));
            if ((a - 1) / prime[i]) nxt = std::min(nxt, (a - 1) / ((a - 1) / prime[i]));
            ans += 1ll * cnt * (cnt - 1) / 2 * (sum[nxt] - i);
            i = sum[nxt];
        }
        printf("%lld\n", ans);
        fflush(stdout);
    }
    return 0;
}

Spring II

It was a surprise to discover a new nest opposite the window.

Through the telescope, I scrupulously watched the doves in the nest, who were reposing restfully. Narrowly I discerned that they were a family comprised of a mother and two babies. Sighing with relief, I thought, “Finally the dove has returned, and luckily I can still see the family.” I even expected to witness the first flight of the baby doves.

I started wondering where the doves came from, as it was odd that they were inclined to construct their home here between the buildings, and after a period of cogitating, I ultimately deduced that once they might live in the decrepit houses nearby which were to be pulled down owing to the modernization and urbanization.

From this viewpoint, the doves seemed to be pitiable. “It can be tough for them to find a congenial shelter in the city.” I uttered in my heart, as I observed them fluttering their gray wings. It looked as if they were enjoying the serene life. Yet my brain did not cease contemplating.

Living in the era with rapid technological advancements which largely hasten our pace of life, we ineluctably forgo the joyous time being close to nature. And not until the demise of the erstwhile happiness we attained from nature comes do we realize that it is deplorable to turn ourselves into a hardened machine.

The blood-red sun was sinking and a livid cloud in the east received its rays. A twinge of woe surged through me. If we humanity kept pursuing interests by damaging the environment, what would the attendant consequence be? Would it only be driving the doves homeless?

I couldn’t help recalling the days when the city was shrouded in thick smog, and people could scarcely distinguish the things in the distance. And the rivers in the city, which were limpid before, had become turbid in recent decades. The foregoing instances were attributed to human activities. Not to mention the global warming, the ozone depletion, and the extinction of multitudinous species. Urbanization was just all of these in miniature.

The sky was tinted with black inch by inch, and the doves quieted down again. I shifted my eyes from them.

I leaned on the windowsill, looking up at the firmament, which used to be spangled with innumerable stars that were now unseen to us. Embedded aloft in the inky night sky, the pale yellow full moon cast its glow over the city, like an infinite piece of cloth mantling, yet faraway as the neon lights were, their harsh beams penetrated the textile that the moon carefully sewed without any difficulty, as if it was nonexistent. The dazzling rays pricked not only my eyes, but also my heart. “Perhaps the doves dread the blinding light, too.”

There was a time when all the doves could flit from tree to tree jovially in the dense and emerald woodlands, where they interacted with other creatures complying with the principles of nature. It was not until human beings arose that changes occurred. Trees were felled at first, and then sewage and fumes were disgorged by pipes. As the time elapsed our compunction was attenuated unawares by the profits we gained at the price of the environment we lived. We averted our gaze from what was prime, focusing on what was subordinate. The cupidity for gain and the scant heed we paid to conservation incurred the mire we confront today.

We really should pause, ruminate for a while, and inquire ourselves, what is our original orientation, and does what we endeavor to do diverge from that aim. Is our purpose to make every hue around us grey? Or to fetter ourselves in the forest of reinforced concrete? Or, to alter the earth to meet human’s multifarious demands regardless of the sustainability?

The onus is on us to maintain an agreeable environment. We need introspection, though it may seems arduous to atone for the devastation.

Fortunately, we gradually awake to the pernicious effects, and we are making a change. Eco-civilization construction has been put forward, and the departments concerned are sparing no effort to ameliorate the imbalanced ecology.

It was late spring, and every plant was flourishing. I heard the birds warbling a mellifluous night chorus, my heart mollified.

I entertain the hope that one day, we can thoroughly relish the felicity from nature. I envision an unspoiled environment. If we esteem nature, make optimal decision, and take prompt action.

I’m looking forward to the arrival of that very day.

Spring I

I heard something twittering on the balcony.

Walking on my tiptoe, I reached the window, through which I spotted a turtle dove standing on the stainless-steel-made railing, with a sprig in its beak. Barely had I observed it attentively when it hastily fled away, leaving behind a fuzzy silhouette.

I was leaning on the windowsill, looking down, when a ball of white captivated me, which I instantly recognized as peach blossoms. Bathed in radiant sunlight, they stretched their exquisite bodies in all directions. The balmy breeze gently swayed the flowers, carrying the scent they exuded wherever it passed. Eventually here came the spring. As I was about to return to my room, my eyes rested on a half-finished bird nest that lay on the board put on the railings.

It was unbelievable that the dove was planning to build a nest on our balcony. Yet this was the case. What’s more, I witnessed two doves erecting their lovely home together. Before long, a complete nest came out.

I could see in the mild twilight, a tiny nest lying steadily behind the huge flower pots, and an adorable dove perching on the nest. The setting sun cast its silky beams on this little cute creature through the verdant plants in the pots, resembling a warm hand patting tenderly. This tranquil scene was set off by the distant golden horizon, which formed the most enchanting landscape in the world. I could picture that under its little body, eggs were being warmed up, and new life was being brewed.

But we all know that good times do not last long, which I grasped on that drizzly evening when the dove had gone away and had not returned.

I stared at the empty nest vacantly, a sense of loss welling up in my heart. There lay the sole egg placidly, reflecting the cold light which the moon shed on the earth. Why didn’t the dove return? Maybe, I unintentionally frightened it, or, it had found a superior habitat. But, did it desert its cherished egg? On this fairly chilly rainy night, where could it go? Would it be shivering all over in the freezing wind? As the moonlight suffused the ground, I perceived that, never would I know the answer.

I leaned on the windowsill, overlooking the swaying flowers. I could hear raindrops knocking the canopy, composing a rumpled melody agitating my thoughts. I could see them trickling down the eaves, descending, sinking into puddles, feebly splattering a spray of droplets. Also I could hear the moaning wind dashing through the alley, wildly rattling loose windows and doors. I could see it wobbling the leaves furiously, severing them from branches and tugging them down. But, in the meantime, I could see the blossoms that the tree threw up in the air, being hit by the rain and being swung by the wind, while they kept emanating aroma that wafted all around. I could see the myriad scattered petals which fell on the moist tarmac, all colors draining rapidly, yet somehow, they were peculiarly appealing with the lamplight flickering in dewdrops.

It suddenly occurred to me that, unlike what I considered before, what was called spring, was not a picture in a particular style. It not only meant arrival and growth, but also meant departure and wilting to some extent.

So the spring is not all about a couple of turtle doves who are busy constructing a cozy neat twiggy nest where the new life will be brought, nor about the pale pink blossoms that bloom wherever you behold, adorning the corner of the bright blue sky. The essence of the spring lies in the persistence that the doves possess, even if they get so close to human beings and lose one of the only two eggs. It lies in the patience with which the dove squats in its nest unceasingly, awaiting the arrival of the baby, no matter the icy drizzle and the roaring gale. It lies in the optimism which encourages the blossoms to bloom vigorously in the frigid sprinkle. It lies in the devotion of the blossoms as they make the world transfigured despite the adversity they are confronted with.

I lingered on the balcony, hearkening wind and rain, dismay dissipating gradually. And here came the spring. I could sense that, distinctly.

The Beginning

The bell is about to ring. The new year is in sight.

Time flies as we grow and mature, leaving naivety behind. No sooner have we relived the past year in detail than we face the coming year hurriedly. Lay down the heavy burden, repack the travelling bag, take our initial dream, step forward readily and firmly.

It’s said those who have received compliments should be rather careful, for other’s acclaim may block the consciousness of weakness, which surely leads to failure. Yet by no means shall we belittle ourselves, since nobody will have confidence in us unless we are confident. Owning a distinct understanding of oneself, cultivating the better part, remedying what have been done that turns out a mistake, are ways to adjust our mentality.

When we gaze at the vast, dark night sky, sparsely-arranged stars hang high, whispering about what they have just witnessed on this tiny, lonely, but vigorous cerulean planet, such as a large amount of teenagers burning the midnight oil, waiting for that extraordinary moment, after which they may become an adult, or may bare his or her heart to someone else, or, may just sigh with the feeling that it seems something in life, which is hard to be aware of, has vanished silently and drastically. Whatever they see, the sparkling stars keep extending their best wishes to whom they have not known and may never know. The stars, the glinting stars, the stars that split the profound night sky, for how long will we still be able to look up at you?

When we overlook the pure white snowfield, accumulated snowflakes lie serenely, mirroring the warm, golden light that streetlamps cast, illumining the bottom of passers-by’s hearts, which originally filled with some annoying things, like working overtime without gaining anything, or getting into a conflict with someone he or she considers dearest, or, biting, frigid wind blowing right into cheek, as if thousands of sharp knives are pricking the skin. But at the time, the light disperses the dust, and pushes the haze in the heart aside, bringing peace and quiet, awakening hope for the future, for the year approaching. The snow, the crystal snow, the snow that decorates the empty city streets, will you still be here years later?

No matter how we get through the year just past, it doesn’t make any sense to immerse in the memory. We shall always look forward, but don’t forget, up in the sky, down on the ground, there must be someone else, caring you, supporting you, giving you hope whenever it seems impossible.

Anyhow, it’s a brand new beginning.

Goodbye, 2018. Hello, 2019.