画像処理でコンクリートのひび割れ抽出

## はじめに

以下の論文のアルゴリズムを実装してみました。

Fujita, Y.; Nakamura, H.; Hamamoto, Y. Automatic and exact crack extraction from concrete surfaces using image processing techniques. J. Jpn. Soc. Civ. Eng. 2010, 66, 459–470.

https://www.jstage.jst.go.jp/article/jscejf/66/3/66_3_459/_article/-char/ja/#citedby-wrap

## 結果

![sample02.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/50315a86-4c64-40ec-5b67-b08f145a3b73.jpeg)

![mask.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/c2aedd43-16c6-f63d-7c7c-05c76c7119ff.jpeg)

![sample01.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/d714300c-b33f-5479-49d0-ead85a6e1447.jpeg)

![mask.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211712/a8de145e-4430-333c-1379-8a3de3344b24.jpeg)


## 実装

```python

import cv2
import numpy as np
from skimage.feature import hessian_matrix


def imread(filename):
    img = cv2.imread(filename)
    if img is None:
        raise FileNotFoundError()
    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    return img


def image_correction(img):
    """濃淡差や影を補正する."""
    img = cv2.GaussianBlur(img, (7, 7), 0)
    return img


def calc_hessian(img, sigma):
    Hxx, Hxy, Hyy = hessian_matrix(
        img, sigma=sigma, mode='mirror', order='xy',
        use_gaussian_derivatives=False)
    hessian = np.stack([Hxx, Hxy, Hxy, Hyy],
                       axis=2).reshape(*img.shape, 2, 2)
    return hessian


def hessian_emphasis(img, sigma=1.414, alpha=0.25):
    hessian = calc_hessian(img, sigma)
    eig_val, _ = np.linalg.eig(hessian)

    r = np.zeros_like(img, dtype=np.float64)
    for i in range(hessian.shape[0]):
        for j in range(hessian.shape[1]):
            # 各画素の固有値を取得
            lambda1, lambda2 = eig_val[i, j]

            # 小さい方をlambda2にする
            if lambda2 > lambda1:
                temp = lambda1
                lambda1 = lambda2
                lambda2 = temp
            if lambda1 <= 0:
                # 粒状構造
                lambda12 = abs(lambda2) + lambda1
            elif lambda2 < 0 < lambda1 < (abs(lambda2) / alpha):
                # 線状構造
                lambda12 = abs(lambda2) - lambda1 * alpha
            else:
                lambda12 = 0
            r[i, j] = sigma**2 * lambda12
    return r


def multiscale_hessian_emphasis(
        img,
        sigma=1.414,
        alpha=0.25,
        s=1.414,
        num_iteration=4
):
    R = []
    for i in range(num_iteration):
        r = hessian_emphasis(img, sigma=sigma * s**i, alpha=alpha)
        R.append(r)
    R = np.stack(R, axis=0)

    out = np.max(R, axis=0)
    out_scale = np.argmax(R, axis=0) + 1  # scaleを1から始める
    return out, out_scale


def stochastic_relaxation_method(img, alpha=1.0):
    """確率的弛緩法.

    Parameters
    ----------
    img : ndarray
        画素値は0-1に規格化.
        ひびの方が背景よりも輝度が高い.
    alpha : float
        補正用の重み.

    Returns
    -------
    out
        確率的弛緩法を適用後の画像.
    residual
        確率的弛緩法適用前後の2乗誤差.
    """
    P_c = np.log1p(img) / np.log1p(img.max())
    P_b = 1 - P_c

    kernels = [
        np.array([[0, 1, 0], [0, 1, 0], [0, 1, 0]]) / 3,
        np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]]) / 3,
        np.array([[0, 0, 0], [1, 1, 1], [0, 0, 0]]) / 3,
        np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) / 3
    ]

    Q_c = []
    for kernel in kernels:
        q = cv2.filter2D(P_c, -1, kernel, borderType=cv2.BORDER_DEFAULT)
        Q_c.append(q)
    Q_c = np.stack(Q_c, axis=0)
    Q_b = 1 - Q_c

    P_c_new = np.divide(
        alpha * P_c[None, ...] * Q_c,
        alpha * P_c[None, ...] * Q_c + P_b[None, ...] * Q_b,
        where=P_c[None, ...] * Q_c + P_b[None, ...] * Q_b != 0)

    out = np.max(P_c_new, axis=0)
    return out


def stepwise_threshold_processing(
        img, candidate,
        iter_stp=3,
        k_size=21,
        element=None):

    height, width = img.shape
    if element is None:
        element = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], np.uint8)

    completed = set()
    t = b = r = l = int((k_size - 1) * 0.5)
    for _ in range(iter_stp):
        flag_updated = False
        # 候補領域の外周の位置インデックスを得る
        dilated = cv2.dilate(candidate, element, iterations=1)
        outline = dilated - candidate
        ys, xs = (outline == 1).nonzero()

        result = candidate.copy()
        for x, y in zip(xs, ys):
            if (x, y) in completed:
                # 判定済み領域は除く
                continue
            # 注目画素近傍で大津の二値化のしきい値を求める
            xmin = max(x - l, 0)
            xmax = min(x + r + 1, width)
            ymin = max(y - t, 0)
            ymax = min(y + b + 1, height)
            patch = img[ymin:ymax, xmin:xmax]
            patch_candidate = candidate[ymin:ymax, xmin:xmax]

            if crack_judge(img, (x, y), patch, patch_candidate):
                # しきい値以上であれば新たな候補領域にする
                result[y, x] = 1.0
                flag_updated = True
            else:
                # しきい値未満で背景に確定する
                completed.add((x, y))

        candidate = result.copy()
        if not flag_updated:
            break
    return result


def trans2uint(img):
    img = ((img - img.min()) * (1 / (img.max() - img.min()) * 255)).astype('uint8')
    return img


def extract_crack_mean(patch, patch_candidate):
    mu_c = np.mean(patch[patch_candidate == 1])
    mu_b = np.mean(patch[patch_candidate == 0])
    return mu_c, mu_b


def crack_judge(img, coor, patch, patch_candidate, beta=0.9):
    x, y = coor
    mu_c, mu_b = extract_crack_mean(patch, patch_candidate)
    res = (abs(img[y, x] - mu_c) / abs(img[y, x] - mu_b + 1e-9) * beta) <= 1.0
    return res


def binarization(img):
    img = ((img - img.min()) * (1 / (img.max() - img.min()) * 255)).astype('uint8')
    _, res = cv2.threshold(
        img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    return res


def crack_segmentation(
    img,
    iter_mhe=4,
    iter_srm=2,
    iter_stp=10
):
    # ヘッシアンの固有値別に画像を補正
    img, _ = multiscale_hessian_emphasis(
        img,
        num_iteration=iter_mhe
    )
    emphasis_img = img.copy()
    emphasis_img = trans2uint(emphasis_img)

    # 確率的弛緩法でノイズを消す
    kernel = np.ones((5, 5), np.float32)
    img = cv2.morphologyEx(img, cv2.MORPH_CLOSE, kernel)
    for _ in range(iter_srm):
        img = stochastic_relaxation_method(img, 2.0)

    # 抽出候補を二値画像として得る
    candidate = binarization(img)
    candidate = candidate.astype(np.float32) / 255

    # 消しすぎた候補領域を段階的閾値処理で拡張する
    mask = stepwise_threshold_processing(
        emphasis_img, candidate, iter_stp=iter_stp, k_size=27)

    # 二値化
    mask = binarization(mask)
    return mask


def main():
    img = imread('sample02.jpg')

    # 画像の前処理
    img = image_correction(img)
    img = cv2.bitwise_not(img)
    # セグメンテーション
    mask = crack_segmentation(img)
    cv2.imwrite("mask.jpg", mask)


if __name__ == '__main__':
    main()
```

## 参考文献

- ryousuke kimoto and yuusuke fujita and kei kawamura and yoshihiko hamamoto, "A Fundamental Study On Automated Crack  Measurement for Concrete structures using Image Processing," *Proceedings of the Fuzzy System Symposium.*, vol. 26, pp. 202-202, 2010, doi: 10.14864/fss.26.0.202.0.
https://www.jstage.jst.go.jp/article/fss/26/0/26_0_202/_article/-char/ja/
- Yusuke FUJITA and Yoshihiro MITANI and Yoshihiko HAMAMOTO, "A Method for Crack Extraction on Concrete Surfaces Using Image Processing Techniques.", *Journal of The Japanese Society for Non-Destructive Inspection*, vol. 56, no. 7, pp. 371-377, 2010, doi: 10.11396/jjsndi.56.371.
https://www.jstage.jst.go.jp/article/jjsndi/56/7/56_7_371/_article/-char/ja/
- Pang-jo CHUN and Nozomu KATAOKA and Chihiro MIWA and Kazuaki HASHIMOTO and Mitao OHGA, "CRACK DETECTION FROM AN IMAGE OF CONCRETE SURFACE WITH CONSIDERING STATISTICAL AND GEOMETRICAL CHARACTERISTICS," *Journal of Japan Society of Civil Engineers, Ser. F3 (Civil Engineering Informatics).*, vol. 70, no. 2, pp. I_1-I_8, 2014, doi: 10.2208/jscejcei.70.I_1.
https://www.jstage.jst.go.jp/article/jscejcei/70/2/70_I_1/_article/-char/ja/
- Y. Sato et al., "Tissue classification based on 3D local intensity structures for volume rendering," in IEEE Transactions on Visualization and Computer Graphics, vol. 6, no. 2, pp. 160-180, April-June 2000, doi: 10.1109/2945.856997.
https://ieeexplore.ieee.org/document/856997
- Sato Y, Nakajima S, Shiraga N, Atsumi H, Yoshida S, Koller T, Gerig G, Kikinis R. Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images. Med Image Anal. 1998 Jun;2(2):143-68. doi: 10.1016/s1361-8415(98)80009-1. PMID: 10646760.
https://pubmed.ncbi.nlm.nih.gov/10646760/
- Pang-jo CHUN and Yuri SHIMAMOTO and Kazuaki OKUBO and Chihiro MIWA and Mitao OHGA, "DEEP LEARNING AND RANDOM FOREST BASED CRACK DETECTION FROM AN IMAGE OF CONCRETE SURFACE," *Journal of Japan Society of Civil Engineers, Ser. F3 (Civil Engineering Informatics).*, vol. 73, no. 2, pp. I_297-I_307, 2017, doi: 10.2208/jscejcei.73.I_297.
https://www.jstage.jst.go.jp/article/jscejcei/73/2/73_I_297/_article/-char/ja/
- Barron, J.T. (2020). A Generalization of Otsu’s Method and Minimum Error Thresholding. In: Vedaldi, A., Bischof, H., Brox, T., Frahm, JM. (eds) Computer Vision – ECCV 2020. ECCV 2020. Lecture Notes in Computer Science(), vol 12350. Springer, Cham. https://doi.org/10.1007/978-3-030-58558-7_27.

以上