反光检测

Automatic segmentation and inpainting of specular highlights for endoscopic imaging 论文反光检测部分pythonn实现

代码写的很乱，可直接运行，后期有时间在整理整理，写写注释

import cv2
import numpy as np
import matplotlib.pyplot as plt
import scipy

from skimage import measure, filters, data

img = cv2.imread('../data/fazhi/95.jpg')
w = img.shape[0]
h = img.shape[1]

# model1
cR = img[:, :, 0]
cG = img[:, :, 1]
cB = img[:, :, 2]

cE = 0.2989 * cR + 0.5870 * cG + 0.1140 * cB


def calc_module1_specular_mask(cE, cG, cB, T1):
    p95_cG = cG * 0.95
    p95_cE = cE * 0.95
    p95_cB = cB * 0.95
    # p95_cG = prctile(cG(:), 95)
    # p95_cB = prctile(cB(:), 95)
    # p95_cE = prctile(cE(:), 95)
    rGE = p95_cG / p95_cE
    rBE = p95_cB / p95_cE
    img_new = np.zeros((w, h, 1))
    for i in range(0, w):
        for j in range(0, h):

            if all([(cG[i][j] > rGE[i][j] * T1), (cB[i][j] > rBE[i][j] * T1), (cE[i][j] > T1)]):
                img_new[i][j] = 255
            else:
                img_new[i][j] = 0

    return img_new


module1_specular_mask = calc_module1_specular_mask(cE, cG, cB, T1=240)
cv2.imwrite('../data/New/95-model1.jpg',module1_specular_mask)


# cv2.imshow('s',img)


def calc_centroid_color_info(specular_mask_T2_abs, img):
    kernel = np.ones((2, 2), np.uint8)
    dilated_mask_1 = cv2.erode(specular_mask_T2_abs, kernel, iterations=1)
    kernel = np.ones((4, 4), np.uint8)
    dilated_mask_2 = cv2.erode(specular_mask_T2_abs, kernel, iterations=1)
    centroid_color_area = dilated_mask_2 - dilated_mask_1
    labeled_area = measure.label(centroid_color_area)  # 连通区域标记
    num_region = max(labeled_area.reshape(labeled_area.shape[0] * labeled_area.shape[1]))
    centroid_color_info = []
    for i in range(1, num_region):
        dict = {}
        [row_index, col_index] = np.where(labeled_area == i)
        num_possible_specular_points = len(row_index)
        dict['centroid_row'] = np.mean(row_index)
        dict['centroid_col'] = np.mean(col_index)
        dict['centroid_color'] = np.mean(img[row_index, col_index, :])  ## 不明确
        centroid_color_info.append(dict)

    return centroid_color_info


def calc_distance(x, y, x1, y1):
    distance_to_centroid = np.sqrt((x - x1) ** 2 + (y - y1) ** 2)
    return distance_to_centroid


def find_the_nearest_region(centroid_color_info, pixel_row, pixel_col):
    num_region = len(centroid_color_info)
    nearset_region_index = 0
    nearset_distance = 1e6
    for j in range(1, num_region):
        distance_to_centroid = calc_distance(pixel_row, pixel_col, centroid_color_info[j]['centroid_row'],
                                             centroid_color_info[j]['centroid_col'])
        if distance_to_centroid < nearset_distance:
            nearset_distance = distance_to_centroid
            nearset_region_index = j

    nearest_region = centroid_color_info[nearset_region_index]
    return nearest_region


def filling_image_using_centroid_color(specular_mask_T2_abs, img):
    # filling possible specular highlights by the centroid color
    centroid_color_info = calc_centroid_color_info(specular_mask_T2_abs, img)
    specular_mask_T2_abss = specular_mask_T2_abs.reshape(specular_mask_T2_abs.shape[0], specular_mask_T2_abs.shape[1])
    [row_index, col_index] = np.where(specular_mask_T2_abss == 255)
    num_possible_specular_points = len(row_index)
    filled_img = img
    for i in range(1, num_possible_specular_points):
        # looking for the nearst centroid color for every specular point
        # and fill it
        nearest_region = find_the_nearest_region(centroid_color_info, row_index[i], col_index[i])
        filled_img[row_index[i], col_index[i], :] = nearest_region['centroid_color']
    return filled_img


def contrast_coeffcient(c):
    mean_c = np.mean(c);
    std_c = np.std(c);
    t = 1 / ((mean_c + std_c) / mean_c)
    return t


def calc_modul2_specular_mask(filled_img, T2_rel, cR, cG, cB):
    R = filled_img[:, :, 0]
    fR = cv2.medianBlur(R, 31)
    fG = cv2.medianBlur(filled_img[:, :, 1], 31)
    fB = cv2.medianBlur(filled_img[:, :, 2], 31)
    filtered_img = np.stack((fR, fG, fB), axis=2)

    for i in range(filled_img.shape[0]):
        for j in range(filled_img.shape[1]):
            if (fR[i][j] < 2.2204e-16):
                fR[i][j] = 1e7
            if (fG[i][j] < 2.2204e-16):
                fG[i][j] = 1e7
            if (fB[i][j] < 2.2204e-16):
                fB[i][j] = 1e7
    tR = contrast_coeffcient(cR)
    tG = contrast_coeffcient(cG)
    tB = contrast_coeffcient(cB)

    max_img = np.stack(((tR * cR / fR), (tG * cG / fG), (tB * cB / fB)), axis=2)
    e_max = np.amax(max_img, 2)
    module2_specular_mask = e_max > T2_rel
    return module2_specular_mask

    # fR(fR <  2.2204e-16) = 1e7
    # fG(fG <  2.2204e-16) = 1e7
    # fB(fB <  2.2204e-16) = 1e7


# model2
specular_mask_T2_abs = calc_module1_specular_mask(cE, cG, cB, T1=190)

filled_img = filling_image_using_centroid_color(specular_mask_T2_abs, img)
plt.imshow(filled_img)
plt.show()
module2_specular_mask = calc_modul2_specular_mask(filled_img, T2_rel=1.2, cR=cR, cG=cG, cB=cB)
# module2_specular_mask = np.array(module2_specular_mask)

final_mask = np.zeros((w, h, 1))
for i in range(0, w):
    for j in range(0, h):
        if module2_specular_mask[i][j] == True or module1_specular_mask[i][j] == 255:
            final_mask[i][j][0] = 255
        else:
            final_mask[i][j][0] = 0

N_min = 5000
T3 = 5


def postprocessing(final_mask, cE, N_min, T3):
    kernel = np.ones((3, 3), np.uint8)
    final_mask = cv2.erode(final_mask, kernel, iterations=1)
    labeled_area = measure.label(final_mask)
    num_region = np.max(labeled_area[:])
    post_specular_mask = final_mask
    for i in range(1, num_region):
        index = np.where(labeled_area == i)
        if (len(index) >= N_min):
            post_specular_mask[index] = 0

    return post_specular_mask


mask = postprocessing(final_mask, cE, N_min=3000, T3=5)

mg_gray = cv2.imread('../data/fazhi/95.jpg', cv2.IMREAD_GRAYSCALE)
# 利用图像中要提取的目标区域与其背景在灰度特性上的差异，把图像看作具有不同灰度级的两类区域(目标区域和背景区域)的组合，选取一个比较合理的阈值
thresh = filters.threshold_otsu(mg_gray)

# ret,thresh = cv2.threshold(img,cv2.THRESH_BINARY)
# 根据阈值分割
TTTT = np.zeros((w, h))
dst = (mg_gray >= thresh) * 255.0

for i in range(0, w):
    for j in range(0, h):
        if mask[i][j] > 0 and dst[i][j] > 0:
            TTTT[i][j] = 255

image2 = np.concatenate([TTTT, mask, dst], axis=1)
#cv2.imwrite(r"D:\code dp\YXTX\data\New\TTTT2.jpg", TTTT)
plt.set_cmap("binary")
# plt.imshow(TTTT)
# plt.imshow(mask)
# plt.imshow(dst)
plt.imshow(image2)

plt.show()
cv2.imwrite('../data/New/95_mask.jpg',TTTT)
cv2.imwrite('../data/New/95_mask-image2.jpg',image2)

# cv2.imshow('d', module2_specular_mask)
# cv2.waitKey()
# module2_specular_mask = module2_specular_mask*255
# cv2.imshow('6',filled_img)
# plt.imshow(module2_specular_mask)
# #plt.set_camp('binary')
# plt.show()


### inpainting

decay_win_size = 10
decay_cof = 20


def InpainttingArnold2010(mask, img, decay_win_size, decay_cof):
    filled_img = filling_image_using_centroid_color(mask, img)
    cv2.imwrite('../data/New/95filled.jpg', filled_img)
    # plt.set_cmap("binary")
    # plt.imshow(filled_img)
    # plt.show()
    sig = 8
    gaussian_filtered_img = cv2.GaussianBlur(filled_img, (3, 3), sig, sig)
    cv2.imshow('s',gaussian_filtered_img)

    #mx = imfilter(double(specular_mask), ones(decay_win_size)/decay_cof)
    mx = scipy.ndimage.filters.convolve(mask, np.ones((decay_win_size,decay_win_size))/decay_cof, mode='nearest')
    cv2.imshow('mx',mx)
    cv2.waitKey()
    mx = mx + mask
    mx = (mx > 1) *1
    mx = np.array([mx])
    mx = mx.reshape(w,h,1)
    inpainted_img = mx * (gaussian_filtered_img) + (1 - mx) * img
    return inpainted_img


inpainted_img = InpainttingArnold2010(TTTT, img, decay_win_size, decay_cof)
#plt.imshow(inpainted_img)
cv2.imwrite('../data/New/95inpaint.jpg',inpainted_img)
plt.show()

自己添加一个固定阈值后做并集

import cv2
import numpy as np
import matplotlib.pyplot as plt
import scipy

from skimage import measure, filters, data

img = cv2.imread('../data/fazhi/8.jpg')
w = img.shape[0]
h = img.shape[1]

# model1
cR = img[:, :, 0]
cG = img[:, :, 1]
cB = img[:, :, 2]

cE = 0.2989 * cR + 0.5870 * cG + 0.1140 * cB


def calc_module1_specular_mask(cE, cG, cB, T1):
    p95_cG = cG * 0.95
    p95_cE = cE * 0.95
    p95_cB = cB * 0.95
    # p95_cG = prctile(cG(:), 95)
    # p95_cB = prctile(cB(:), 95)
    # p95_cE = prctile(cE(:), 95)
    rGE = p95_cG / p95_cE
    rBE = p95_cB / p95_cE
    img_new = np.zeros((w, h, 1))
    for i in range(0, w):
        for j in range(0, h):

            if all([(cG[i][j] > rGE[i][j] * T1), (cB[i][j] > rBE[i][j] * T1), (cE[i][j] > T1)]):
                img_new[i][j] = 255
            else:
                img_new[i][j] = 0

    return img_new


module1_specular_mask = calc_module1_specular_mask(cE, cG, cB, T1=240)
cv2.imwrite('../data/New/8-model1.jpg',module1_specular_mask)


# cv2.imshow('s',img)


def calc_centroid_color_info(specular_mask_T2_abs, img):
    kernel = np.ones((4, 4), np.uint8)
    dilated_mask_1 = cv2.erode(specular_mask_T2_abs, kernel, iterations=1)
    kernel = np.ones((6, 6), np.uint8)
    dilated_mask_2 = cv2.erode(specular_mask_T2_abs, kernel, iterations=1)
    centroid_color_area = dilated_mask_2 - dilated_mask_1
    labeled_area = measure.label(centroid_color_area)  # 连通区域标记
    num_region = np.max(labeled_area.reshape(labeled_area.shape[0] * labeled_area.shape[1]))
    centroid_color_info = []
    for i in range(1, num_region):
        dict = {}
        [row_index, col_index] = np.where(labeled_area == i)
        num_possible_specular_points = len(row_index)
        dict['centroid_row'] = np.mean(row_index)
        dict['centroid_col'] = np.mean(col_index)
        ii = img[row_index, col_index, :]
        dict['centroid_color'] = np.mean(ii,axis=0)  ## 不明确
        centroid_color_info.append(dict)

    return centroid_color_info


def calc_distance(x, y, x1, y1):
    distance_to_centroid = np.sqrt((x - x1) ** 2 + (y - y1) ** 2)
    return distance_to_centroid


def find_the_nearest_region(centroid_color_info, pixel_row, pixel_col):
    num_region = len(centroid_color_info)
    nearset_region_index = 0
    nearset_distance = 1e6
    for j in range(1, num_region):
        distance_to_centroid = calc_distance(pixel_row, pixel_col, centroid_color_info[j]['centroid_row'],
                                             centroid_color_info[j]['centroid_col'])
        if distance_to_centroid < nearset_distance:
            nearset_distance = distance_to_centroid
            nearset_region_index = j

    nearest_region = centroid_color_info[nearset_region_index]
    return nearest_region


def filling_image_using_centroid_color(specular_mask_T2_abs, img):
    # filling possible specular highlights by the centroid color
    centroid_color_info = calc_centroid_color_info(specular_mask_T2_abs, img)
    specular_mask_T2_abss = specular_mask_T2_abs.reshape(specular_mask_T2_abs.shape[0], specular_mask_T2_abs.shape[1])
    [row_index, col_index] = np.where(specular_mask_T2_abss > 0)
    num_possible_specular_points = len(row_index)
    filled_img = img
    for i in range(1, num_possible_specular_points):
        # looking for the nearst centroid color for every specular point
        # and fill it
        nearest_region = find_the_nearest_region(centroid_color_info, row_index[i], col_index[i])
        filled_img[row_index[i], col_index[i], :] = nearest_region['centroid_color']
    return filled_img


def contrast_coeffcient(c):
    mean_c = np.mean(c);
    std_c = np.std(c);
    t = 1 / ((mean_c + std_c) / mean_c)
    return t


def calc_modul2_specular_mask(filled_img, T2_rel, cR, cG, cB):
    R = filled_img[:, :, 0]
    fR = cv2.medianBlur(R, 31)
    fG = cv2.medianBlur(filled_img[:, :, 1], 31)
    fB = cv2.medianBlur(filled_img[:, :, 2], 31)
    filtered_img = np.stack((fR, fG, fB), axis=2)

    for i in range(filled_img.shape[0]):
        for j in range(filled_img.shape[1]):
            if (fR[i][j] < 2.2204e-16):
                fR[i][j] = 1e7
            if (fG[i][j] < 2.2204e-16):
                fG[i][j] = 1e7
            if (fB[i][j] < 2.2204e-16):
                fB[i][j] = 1e7
    tR = contrast_coeffcient(cR)
    tG = contrast_coeffcient(cG)
    tB = contrast_coeffcient(cB)

    max_img = np.stack(((tR * cR / fR), (tG * cG / fG), (tB * cB / fB)), axis=2)
    e_max = np.amax(max_img, 2)
    module2_specular_mask = e_max > T2_rel
    return module2_specular_mask

    # fR(fR <  2.2204e-16) = 1e7
    # fG(fG <  2.2204e-16) = 1e7
    # fB(fB <  2.2204e-16) = 1e7


# model2
specular_mask_T2_abs = calc_module1_specular_mask(cE, cG, cB, T1=190)

filled_img = filling_image_using_centroid_color(specular_mask_T2_abs, img)
#plt.imshow(filled_img)
plt.show()
module2_specular_mask = calc_modul2_specular_mask(filled_img, T2_rel=1.2, cR=cR, cG=cG, cB=cB)
# module2_specular_mask = np.array(module2_specular_mask)

final_mask = np.zeros((w, h, 1))
for i in range(0, w):
    for j in range(0, h):
        if module2_specular_mask[i][j] == True or module1_specular_mask[i][j] == 255:
            final_mask[i][j][0] = 255
        else:
            final_mask[i][j][0] = 0

N_min = 5000
T3 = 5


def postprocessing(final_mask, cE, N_min, T3):
    kernel = np.ones((3, 3), np.uint8)
    final_mask = cv2.erode(final_mask, kernel, iterations=1)
    labeled_area = measure.label(final_mask)
    num_region = np.max(labeled_area[:])
    post_specular_mask = final_mask
    for i in range(1, num_region):
        index = np.where(labeled_area == i)
        if (len(index) >= N_min):
            post_specular_mask[index] = 0

    return post_specular_mask


mask = postprocessing(final_mask, cE, N_min=3000, T3=5)

mg_gray = cv2.imread('../data/fazhi/8.jpg', cv2.IMREAD_GRAYSCALE)
# 利用图像中要提取的目标区域与其背景在灰度特性上的差异，把图像看作具有不同灰度级的两类区域(目标区域和背景区域)的组合，选取一个比较合理的阈值
thresh = filters.threshold_otsu(mg_gray)

# ret,thresh = cv2.threshold(img,cv2.THRESH_BINARY)
# 根据阈值分割
TTTT = np.zeros((w, h))
dst = (mg_gray >= thresh) * 255.0

for i in range(0, w):
    for j in range(0, h):
        if mask[i][j] > 0 and dst[i][j] > 0:
            TTTT[i][j] = 255


## 固定阈值

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

blurred = cv2.GaussianBlur(gray, (11, 11), 0)

# threshold the image to reveal light regions in the
# blurred image

#y = blurred[:,:,2]/(blurred[:,:,1] + blurred[:,:,2] + blurred[:,:,3])

th = cv2.threshold(blurred, 200, 255, cv2.THRESH_BINARY)[1]


for i in range(0, w):
    for j in range(0, h):
        if TTTT[i][j]>0 or th[i][j]>0:
            TTTT[i][j]=255




image2 = np.concatenate([TTTT, mask, dst], axis=1)
#cv2.imwrite(r"D:\code dp\YXTX\data\New\TTTT2.jpg", TTTT)
plt.set_cmap("binary")
# plt.imshow(TTTT)
# plt.imshow(mask)
# plt.imshow(dst)
#plt.imshow(image2)

plt.show()
cv2.imwrite('../data/New/8_mask+uzhi.jpg',TTTT)