import cv2, numpy as np
from skimage.morphology import skeletonize

src = "/home/gxs/Desktop/data/2300图像文件测试数据/R4407845900002021120182/R4407845900002021120182_01_01_01_01.bmp"
img = cv2.imread(src, cv2.IMREAD_GRAYSCALE).astype(np.float32)
h, w = img.shape

strip = img[200:260, 100:300].mean(axis=0)
strip = strip - strip.mean()
ac = np.correlate(strip, strip, 'full')[len(strip)-1:]
ac = ac / ac[0]
peak = 4 + int(np.argmax(ac[4:30]))
print(f"estimated ridge period: {peak}px")

resp = np.zeros_like(img)
for freq in (0.10, 0.125, 0.15):
    for theta in range(0, 180, 15):
        k = cv2.getGaborKernel((21, 21), 0.55, theta, 1.0/freq, 0.5, 0, ktype=cv2.CV_32F)
        g = cv2.filter2D(img, -1, k)
        resp = np.maximum(resp, g)

resp_n = (resp - resp.min()) / (resp.max() - resp.min())
den = cv2.GaussianBlur(img.astype(np.uint8), (5,5), 0)
_, gotsu = cv2.threshold(den, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
print_region = (gotsu == 0).astype(np.uint8)
if print_region.mean() > 0.5:
    print_region = 1 - print_region

mask = (resp_n > 0.45).astype(np.uint8)
mask = mask * print_region
k2 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (2,2))
mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, k2)
k3 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3))
mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, k3)
print(f"gabor ridge fraction: {mask.mean():.3f}")
skel = skeletonize(mask > 0)
print(f"skeleton px: {skel.sum()}")
canvas = np.hstack([img.astype(np.uint8), (mask*255).astype(np.uint8), (255*(skel>0)).astype(np.uint8)])
cv2.imwrite('/tmp/fp_gabor.png', canvas)
print('saved')