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

# Gabor bank -> ridge signal (dark periodic lines)
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)   # high = periodic dark ridge

resp_n = (resp - resp.min()) / (resp.max() - resp.min())

# print region
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.55).astype(np.uint8) * print_region
print(f"gabor ridge fraction: {mask.mean():.3f}")

# orientation-aware closing: for each 32x32 block, close along local ridge dir
# (ridge direction = perpendicular to gradient of resp_n)
gx = cv2.Sobel(resp_n, cv2.CV_32F, 1, 0)
gy = cv2.Sobel(resp_n, cv2.CV_32F, 0, 1)
ang = np.arctan2(gy, gx)  # gradient direction; ridge runs perpendicular
closed = mask.copy()
B = 32
for by in range(0, h, B):
    for bx in range(0, w, B):
        a = ang[by:by+B, bx:bx+B].mean()
        # ridge direction = a + 90deg
        L = 5
        k = np.zeros((2*L+1, 2*L+1), np.uint8)
        for t in range(-L, L+1):
            px = L + int(round(t * np.cos(a + np.pi/2)))
            py = L + int(round(t * np.sin(a + np.pi/2)))
            if 0 <= px < 2*L+1 and 0 <= py < 2*L+1:
                k[py, px] = 1
        closed[by:by+B, bx:bx+B] = cv2.morphologyEx(
            mask[by:by+B, bx:bx+B], cv2.MORPH_CLOSE, k, iterations=1)
print(f"after orientation closing: {closed.mean():.3f}")

skel = skeletonize(closed > 0)
print(f"skeleton px: {skel.sum()}")
canvas = np.hstack([img.astype(np.uint8), (mask*255).astype(np.uint8), (closed*255).astype(np.uint8), (255*(skel>0)).astype(np.uint8)])
cv2.imwrite('/tmp/fp_gabor2.png', canvas)
print('saved')