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,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
// 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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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.

Coding > Digital Image Processing
#C/C++ #Course Project
Digital Image Processing Assignment II
https://regmsif.cf/2019/10/11/coding/digital-image-processing-assignment-ii/
作者
RegMs If
发布于
2019年10月11日
许可协议