exp-hough-circle.py

# -*- coding: utf-8 -*-
# %%
# # !jupytext exp-hough-circle.ipynb --to py:percent -o- | tee exp-hough-circle_.py
# # !jupyter nbconvert exp-hough-circle.ipynb --to html 

# %% [markdown]
# ## utility 함수 로딩

# %%
# %run ../ml-workspace-common/ml-workspace-common.py

# %% [markdown]
# ## 2020.12.11_2 데이터에서 100개를 샘플로 가져와서 로컬에 1/2 사이즈로 저장

# %%
# # !rm -fvr data

# %%
if not isdir('data'):

    makedirs('data', exist_ok=True)

    src_dir = '/G/Train_ImageSet/Zebrafish/2020.12.11_2/20201211/2차/'
    files = glob(src_dir+'*.BMP')

    np.random.seed(1)
    choices = np.random.choice(files, size=(100,), replace=False)

    itr = tqdm(choices)
    for f in itr:
        im = cv2.imread(f, flags=cv2.IMREAD_UNCHANGED)
        im_conv = cv2.resize(im, None, fx=0.5, fy=0.5,
                             interpolation=cv2.INTER_AREA)
        nf = 'data/' + splitext(basename(f))[0] + '.png'
        cv2.imwrite(nf, im_conv)
        itr.set_postfix_str(nf)


# %%
# # !ls -s data

# %%
# def show_gray_image(im,vmin=0,vmax=255):
#     print('show_gray_image: shape:',im.shape)
#     assert im.ndim == 2
#     fig, axs = plt.subplots(1,2,figsize=(11,4))
#     axs[0].imshow(im,vmin=vmin,vmax=vmax)
#     axs[1].hist(im.flatten())
#     fig.show()    

# %% [markdown]
# ## `find_circles()` - houghcircles 를 이용한 탐지 함수

# %% [raw]
# """
# Python:
#
# circles = cv.HoughCircles( image, method, dp, minDist[, circles[, param1[, \
#             param2[, minRadius[, maxRadius]]]]] )
#
#
#
# Parameters
#     image 
#         8-bit, single-channel, grayscale input image.
#     circles 
#         Output vector of found circles. Each vector is encoded as 3 or 
#         4 element floating-point vector (x,y,radius) or (x,y,radius,votes) .
#     method 
#         Detection method, see HoughModes. Currently, the only implemented 
#         method is HOUGH_GRADIENT
#     dp 
#         Inverse ratio of the accumulator resolution to the image resolution. 
#         For example, if dp=1 , the accumulator has the same resolution as 
#         the input image. If dp=2 , the accumulator has half as big width 
#         and height.
#     minDist 
#         Minimum distance between the centers of the detected circles. 
#         If the parameter is too small, multiple neighbor circles may be 
#         falsely detected in addition to a true one. If it is too large, 
#         some circles may be missed.
#     param1 
#         First method-specific parameter. In case of HOUGH_GRADIENT , it is 
#         the higher threshold of the two passed to the Canny edge detector 
#         (the lower one is twice smaller).
#     param2 
#         Second method-specific parameter. In case of HOUGH_GRADIENT , it is 
#         the accumulator threshold for the circle centers at the detection 
#         stage. The smaller it is, the more false circles may be detected. 
#         Circles, corresponding to the larger accumulator values, will be 
#         returned first.
#     minRadius 
#         Minimum circle radius.
#     maxRadius 
#         Maximum circle radius. If <= 0, uses the maximum image dimension. 
#         If < 0, returns centers without finding the radius.
# """;

# %%
def find_circles(im,
                 sf=8.0,  # 이미지를 1/8 으로 축소한 후 탐색
                 channel=0,  # R 채널을 gray로 사용
                 mb_ksize=3, # median blur kernel size
                 eq_clip=2.0, # clahe clip limit
                 eq_grid=8, # clahe tile grid size
                 hc_dp=3.65,  # image / accumulator ratio for hough circle
                 hc_min_dist_r=1.5,  # 한 화면에 한개 서클
                 hc_param1=100,  # param1 value
                 hc_param2=20,  # param2 value
                 hc_min_dia_r=0.2,  # 최소 0.2 크기
                 hc_max_dia_r=0.4,  # 최대 0.4 크기
                 circ_pad=0.0,  # 내부 패딩 비율
                 _debug=None,
                 ):
    resized = cv2.resize(im, None, fx=1/sf, fy=1/sf,
                         interpolation=cv2.INTER_AREA)
    image_size = min(*resized.shape[:2])
    gray = resized[:, :, channel]
    blur = cv2.medianBlur(gray, ksize=mb_ksize)
#     norm = cv2.normalize(blur, gray.astype(float), 0,
#                          255, norm_type=cv2.NORM_MINMAX)
#     norm = cv2.equalizeHist(blur)
    clahe = cv2.createCLAHE(clipLimit=eq_clip, tileGridSize=(eq_grid,eq_grid))
    norm = clahe.apply(blur)
    method = cv2.HOUGH_GRADIENT
    minDist = hc_min_dist_r * image_size
    minRadius = int(0.5 * hc_min_dia_r * image_size)
    maxRadius = int(0.5 * hc_max_dia_r * image_size)
    circles = cv2.HoughCircles(norm, method,
                               dp=hc_dp,
                               minDist=minDist,
                               param1=hc_param1,
                               param2=hc_param2,
                               minRadius=minRadius,
                               maxRadius=maxRadius)
    if circles is not None:
        circles = np.round(circles*[[[sf, sf, sf/(1.0-circ_pad)]]]).astype(int)
    if isinstance(_debug, dict):
        canny = cv2.Canny(norm, hc_param1, hc_param1//2)
        preview = im.copy()
        if circles is not None:
            image_size = min(*preview.shape[:2])
            for x, y, r in circles[0]:
                cv2.circle(preview, (x, y), r, (255, 128, 128), 15)
        _debug.update(locals())
    return circles


def _imshow(im, vmin=0, vmax=255, title=None, ax=None, **kwargs):
    ax = ax or plt.gca()
    ax.imshow(im, vmin=vmin, vmax=vmax, **kwargs)
    if title:
        ax.set_title(title)


# %%

# %%

# %%

# %% [markdown]
# ## 100개 샘플 중 20개에 대해서 시험

# %%
np.random.seed(0)
files = sorted(glob('data/*.png'))
samples = np.random.choice(files, size=(20,), replace=False)


for im_name in samples:
    im = cv2.imread(im_name, flags=cv2.IMREAD_UNCHANGED)
    im = im[:, :, ::-1]
    print(im_name, im.shape, flush=True)
    params = dict(
#         sf=8.0,  # 이미지를 1/8 으로 축소한 후 탐색
#         channel=0,  # R 채널을 gray로 사용
        mb_ksize=3,  # median blur kernel size
#         eq_clip=2.0,  # clahe clip limit
        eq_grid=8,  # clahe tile grid size
#         hc_dp=3.65,  # image / accumulator ratio for hough circle
#         hc_min_dist_r=1.5,  # 한 화면에 한개 서클
#         hc_param1=100,  # param1 value
#         hc_param2=20,  # param2 value
#         hc_min_dia_r=0.2,  # 최소 0.2 크기
#         hc_max_dia_r=0.4,  # 최대 0.4 크기
#         circ_pad=0.0,  # 내부 패딩 비율
    )
    _debug = dict()
    %time circles = find_circles(im, _debug=_debug, **params)
    if circles is not None:
        for x, y, r in circles[0,:]:
            image_size = min(*im.shape[:2])
            print('circle: ({:d},{:d},{:d}) rel: {:.1%}'.format(
                x, y, r, r*2/image_size))
    _, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(11, 2.5))
    _imshow(_debug['preview'], title=im_name, ax=ax1)
    _imshow(_debug['norm'], title='norm', ax=ax2)
    ax3.hist(_debug['gray'].flatten(), bins=32, alpha=0.2, label='gray')
    ax3.hist(_debug['norm'].flatten(), bins=32, alpha=0.2, label='norm')
    ax3.set_title('hist')
    _imshow(_debug['canny'], title='canny', ax=ax4)
    plt.tight_layout()
    plt.show()

# %%

# %%

# %%

# %% [markdown]
# ## 100개 샘플중 20개에 대해서 이근호 팀장 코드 시험

# %%
import math


def crop_circle_img(img, x,  y, r):
    rectX = x - r
    if rectX < 0:
        rectX = 0
    rectY = y - r
    if rectY < 0:
        rectY = 0
    maxX = x + r
    if maxX >= img.shape[1]:
        maxX = img.shape[1]
    maxY = y + r
    if maxY >= img.shape[0]:
        maxY = img.shape[0]
    crop_img = img[rectY:maxY, rectX:maxX]
    return crop_img

def average_color(crop_img, mode, r):
    if mode == 0:
        STANDARD_MEAN = [0, 255, 255, 0]
    else:
        STANDARD_MEAN = [0, 0, 0, 0]
    if crop_img.shape[0] != crop_img.shape[1]:
        return (99999999999)
    circle_img = np.zeros((crop_img.shape[0], crop_img.shape[1]), np.uint8)
    y = int(crop_img.shape[0] / 2)
    x = int(crop_img.shape[1] / 2)
    cv2.circle(circle_img, (x, y), r, (255, 255, 255), -1)
    datos_rgb = cv2.mean(crop_img, mask = circle_img)[::-1]
    list_datos_rgb = np.array(list(datos_rgb)) - np.array(list(STANDARD_MEAN))
    return np.sum(list_datos_rgb ** 2)

def detect_circle(img, cropsize):
    dst = img.copy()
    gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    blur_img = cv2.blur(gray_img, (9, 9))
    sobel_img_x = cv2.Sobel(blur_img, cv2.CV_64F, 1, 0, 3)
    sobel_img_x = cv2.convertScaleAbs(sobel_img_x)
    sobel_img_y = cv2.Sobel(gray_img, cv2.CV_64F, 0, 1, 3)
    sobel_img_y = cv2.convertScaleAbs(sobel_img_y)
    sobel_img = cv2.addWeighted(sobel_img_x, 1, sobel_img_y, 1, 0)
    ret, sobel_img = cv2.threshold(sobel_img, 20, 255, cv2.THRESH_BINARY)
    inner_circles = cv2.HoughCircles(sobel_img, cv2.HOUGH_GRADIENT, 1, 1000, param1 = 50, param2 = 20, minRadius = 140, maxRadius = 200)
    outter_circles = cv2.HoughCircles(sobel_img, cv2.HOUGH_GRADIENT, 1, 1000, param1 = 50, param2 = 20, minRadius = 240, maxRadius = 300)
    hh, ww = img.shape[:2]
    inner_circle_lst = []
    outter_circle_lst = []
    if inner_circles is None and outter_circles is None:
        return None
    if inner_circles is not None:
        for i in inner_circles[0]:
            a = (hh / 2) - int(i[1])
            b = (ww / 2) - int(i[0])
            crop_circle = crop_circle_img(img, int(i[0]), int(i[1]), int(i[2]))
            circle_mean = average_color(crop_circle, 0, int(i[2]))
            inner_circle_lst.append((math.sqrt((a * a) + (b * b)), circle_mean))
    if outter_circles is not None:
        for i in outter_circles[0]:
            a = (hh / 2) - int(i[1])
            b = (ww / 2) - int(i[0])
            crop_circle = crop_circle_img(img, int(i[0]), int(i[1]), int(i[2]))
            circle_mean = average_color(crop_circle, 1, int(i[2]))
            outter_circle_lst.append((math.sqrt((a * a) + (b * b)), circle_mean))
    inner_ar = np.array(inner_circle_lst)[:, 1].tolist()
    idx = inner_ar.index(min(inner_ar))
    xx = int(inner_circles[0][idx][0])
    yy = int(inner_circles[0][idx][1])
    if inner_circle_lst[idx][0] > 900:
        if len(outter_circle_lst) > 0:
            outter_ar = np.array(outter_circle_lst)[:, 1].tolist()
            idx = outter_ar.index(min(outter_ar))
            xx = int(outter_circles[0][idx][0])
            yy = int(outter_circles[0][idx][1])
    crop_half = cropsize // 2
    crop_rest = cropsize - crop_half
    minY = yy - crop_half if yy - crop_half > 0 else 0
    maxY = yy + crop_rest if yy + crop_rest < hh else hh
    minX = xx - crop_half if xx - crop_half > 0 else 0
    maxX = xx + crop_rest if xx + crop_rest < ww else ww
    img_cropped = img[int(minY):int(maxY), int(minX):int(maxX)]
#     if len(img_cropped) <= 0:
#         img_cropped = common.center_crop_img(img, 900)
    return img_cropped


# %%
np.random.seed(0)
files = sorted(glob('data/*.png'))
samples = np.random.choice(files, size=(20,), replace=False)
for im_name in samples:
    im = cv2.imread(im_name,flags=cv2.IMREAD_UNCHANGED)
    im = im[:,:,::-1]
    print(im_name,im.shape,flush=True)
    %time im_crop = detect_circle(im,450)
    _, ax1 = plt.subplots(1,1,figsize=(3.5,2.5))
    _imshow(im_crop,title=im_name,ax=ax1)
    plt.tight_layout()
    plt.show()


# %%

# %%

# %%

# %% [markdown]
# ## `data/` 아래 모든 샘플에 대해서 변환 후 결과를 `preview/` `blur/` `canny/` 로 저장

# %%
# !mkdir -p data/preview data/blur data/canny data/csv data/debug

# %%
def write_circle_csv(csv_name,circles):
    with open(csv_name,'w') as fout:
        if circles is not None:
            for x, y, r in circles[0,:]:
                image_size = min(*im.shape[:2])
                rel = r * 2 / image_size
                print('{:d},{:d},{:d},{:.3f}'.format(x,y,r,rel),file=fout)


# %%
from importlib import reload
import tqdm
reload(tqdm)

# %%
from tqdm import tqdm

files = sorted(glob('data/*.png'))
itr = tqdm(files)
for im_name in itr:
    im = cv2.imread(im_name, flags=cv2.IMREAD_UNCHANGED)[:, :, ::-1]
    itr.set_postfix_str(repr((im_name, im.shape)))
    
    _debug = dict()
    circles = find_circles(im, _debug=_debug)
    
    csv_name = 'data/csv/' + splitext(basename(im_name))[0] + '.csv'
    write_circle_csv(csv_name,circles)
    
    new_name = splitext(basename(im_name))[0] + '.png'
    
    pv_name = 'data/preview/' + new_name
    cv2.imwrite(pv_name, _debug['preview'][:, :, ::-1])
    
#     mb_name = 'data/blur/' + new_name
#     cv2.imwrite(mb_name, _debug['blur'])
#     cn_name = 'data/canny/' + new_name
#     cv2.imwrite(cn_name, _debug['canny'])
    
    debug_name = 'data/debug/' + new_name
    fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(22, 5), dpi=100)
    _imshow(_debug['preview'], title=im_name, ax=ax1)
    _imshow(_debug['norm'], title='norm', ax=ax2)
    ax3.hist(_debug['gray'].flatten(), bins=32, alpha=0.2, label='gray')
    ax3.hist(_debug['norm'].flatten(), bins=32, alpha=0.2, label='norm')
    ax3.set_title('hist')
    _imshow(_debug['canny'], title='canny', ax=ax4)
    fig.tight_layout()
    fig.savefig(debug_name)
    plt.clf()
    plt.close()


# %%

# %%

# %%

# %% [markdown]
# ## break

# %%
assert 0

# %% [markdown]
# ## 전체 Zebrafish 데이터 2020.12 월 분을 `data/` 에 resize 하여 저장

# %%
# !ls /G/Train_ImageSet/Zebrafish/

# %%
# directories:
#   
# /G/Train_ImageSet/Zebrafish/2020.12.09~10/
# /G/Train_ImageSet/Zebrafish/2020.12.11_1/
# /G/Train_ImageSet/Zebrafish/2020.12.11_2/

# %%
# # !ls /G/Train_ImageSet/Zebrafish/2020.12.09~10/*/*/*/*.BMP

# %%
# # !ls /G/Train_ImageSet/Zebrafish/2020.12.11_1/*/*/*.BMP

# %%
# # !ls /G/Train_ImageSet/Zebrafish/2020.12.11_2/*/*/*.BMP

# %%
# 4143it [1:17:52,  1.13s/it, ('/G/Train_ImageSet/Zebrafish/2020.12.11_2/A_plate1210/7차/H9_20201210225407.BMP', (2448, 2048, 3))]      

# %%
from tqdm import tqdm
from itertools import chain
from glob import iglob


files_glob_1 = '/G/Train_ImageSet/Zebrafish/2020.12.09~10/*/*/*/*.BMP'
files_glob_2 = '/G/Train_ImageSet/Zebrafish/2020.12.11_1/*/*/*.BMP'
files_glob_3 = '/G/Train_ImageSet/Zebrafish/2020.12.11_2/*/*/*.BMP'

itr = chain(iglob(files_glob_1),
                 iglob(files_glob_2),
                 iglob(files_glob_3))
num_files = 0
for _ in itr:
    num_files += 1
    
print('total files:', num_files)

# %%
num_converted_files = 0
for _ in iglob('data/*.png'):
    num_converted_files += 1

print('converted files:', num_converted_files)

# %%
if num_converted_files < num_files:
    itr = tqdm(chain(iglob(files_glob_1),
                     iglob(files_glob_2),
                     iglob(files_glob_3)),
               total=num_files)
    for im_name in itr:
        im = cv2.imread(im_name, flags=cv2.IMREAD_UNCHANGED)
        # 1/2x1/2 로 축소
        im = cv2.resize(im, None, fx=0.5, fy=0.5, interpolation=cv2.INTER_AREA)
        new_name = splitext(basename(im_name))[0] + '.png'
        cv2.imwrite('data' + '/' + new_name, im)
        itr.set_postfix_str(repr((im_name, im.shape)))        

# %%

# %%

# %%