line2hist.py

# usage: python -u line2hist.py <inputimage>
import os

# os.chdir("/shm")
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import sys
import math

PHI = 1.6180339887498948482  # ppl says this is a beautiful number :)


def freeman(x, y):
    if y == 0:
        y = 1e-9  # so that we escape the divby0 exception
    if x == 0:
        x = -1e-9  # biased to the left as the text progresses leftward
    if (abs(x / y) < pow(PHI, 2)) and (abs(y / x) < pow(PHI, 2)):  # corner angles
        if (x > 0) and (y > 0):
            return 1
        elif (x < 0) and (y > 0):
            return 3
        elif (x < 0) and (y < 0):
            return 5
        elif (x > 0) and (y < 0):
            return 7
    else:  # square angles
        if (x > 0) and (abs(x) > abs(y)):
            return int(0)
        elif (y > 0) and (abs(y) > abs(x)):
            return 2
        elif (x < 0) and (abs(x) > abs(y)):
            return 4
        elif (y < 0) and (abs(y) > abs(x)):
            return 6


RESIZE_FACTOR = 2
SLIC_SPACE = 3
SLIC_SPACE = SLIC_SPACE * RESIZE_FACTOR

THREVAL = 60
RASMVAL = 160

CHANNEL = 2


def draw(img):  # draw the bitmap
    plt.figure(dpi=600)
    plt.grid(False)
    if len(img.shape) == 3:
        plt.imshow(cv.cvtColor(img, cv.COLOR_BGR2RGB))
    elif len(img.shape) == 2:
        plt.imshow(cv.cvtColor(img, cv.COLOR_GRAY2RGB))


filename = sys.argv[1]
# filename= 'topanribut.png'
imagename, ext = os.path.splitext(filename)
image = cv.imread(filename)
resz = cv.resize(
    image,
    (RESIZE_FACTOR * image.shape[1], RESIZE_FACTOR * image.shape[0]),
    interpolation=cv.INTER_LINEAR,
)
image = resz.copy()
image = cv.bitwise_not(image)
height = image.shape[0]
width = image.shape[1]

image_gray = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
image_gray = image[:, :, CHANNEL]
_, gray = cv.threshold(image_gray, 0, THREVAL, cv.THRESH_OTSU)  # less smear
# _, gray= cv.threshold(selective_eroded, 0, THREVAL, cv.THRESH_TRIANGLE) # works better with dynamic-selective erosion
# draw(gray)
render = cv.cvtColor(gray, cv.COLOR_GRAY2BGR)

# SLIC
cue = gray.copy()
slic = cv.ximgproc.createSuperpixelSLIC(
    cue, algorithm=cv.ximgproc.SLICO, region_size=SLIC_SPACE
)
slic.iterate()
mask = slic.getLabelContourMask()
result_mask = cv.bitwise_and(cue, mask)
num_slic = slic.getNumberOfSuperpixels()
lbls = slic.getLabels()

# moments calculation for each superpixels, either voids or filled (in-stroke)
moments = [np.zeros((1, 2)) for _ in range(num_slic)]
moments_void = [np.zeros((1, 2)) for _ in range(num_slic)]
# tabulating the superpixel labels
for j in range(height):
    for i in range(width):
        if cue[j, i] != 0:
            moments[lbls[j, i]] = np.append(
                moments[lbls[j, i]], np.array([[i, j]]), axis=0
            )
            render[j, i, 0] = 140 - (10 * (lbls[j, i] % 6))
        else:
            moments_void[lbls[j, i]] = np.append(
                moments_void[lbls[j, i]], np.array([[i, j]]), axis=0
            )


# moments[0][1] = [0,0] # random irregularities, not quite sure why
# some badly needed 'sanity' check
def remove_zeros(moments):
    temp = []
    v = len(moments)
    if v == 1:
        return temp
    else:
        for p in range(v):
            if moments[p][0] != 0.0 and moments[p][1] != 0.0:
                temp.append(moments[p])
        return temp


for n in range(len(moments)):
    moments[n] = remove_zeros(moments[n])

# draw(render)

######## // image preprocessing ends here

# generating nodes
scribe = nx.Graph()  # start anew, just in case

# valid superpixel
filled = 0
for n in range(num_slic):
    if (
        len(moments[n]) > SLIC_SPACE
    ):  # remove spurious superpixel with area less than 2 px
        cx = int(np.mean([array[0] for array in moments[n]]))  # centroid
        cy = int(np.mean([array[1] for array in moments[n]]))
        if cue[cy, cx] != 0:
            render[cy, cx, 1] = 255
            scribe.add_node(
                int(filled),
                label=int(lbls[cy, cx]),
                area=(len(moments[n]) - 1) / pow(SLIC_SPACE, 2),
                hurf="",
                pos_bitmap=(cx, cy),
                pos_render=(cx, -cy),
                color="#FFA500",
                rasm=True,
            )
            # print(f'point{n} at ({cx},{cy})')
            filled = filled + 1


def pdistance(point1, point2):
    x1, y1 = point1
    x2, y2 = point2
    distance = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
    return distance


# connected componentscv.circle(disp, pos[compodef line_iterator(img, point0, point1):
from dataclasses import dataclass, field
from typing import List
from typing import Optional


@dataclass
class ConnectedComponents:
    rect: (int, int, int, int)  # from bounding rectangle
    centroid: (int, int)  # centroid moment
    area: Optional[int] = field(default=0)
    nodes: List[int] = field(default_factory=list)
    mat: Optional[np.ndarray] = field(default=None, repr=False)
    node_start: Optional[int] = field(default=-1)  # right-up
    distance_start: Optional[int] = field(default=0)  # right-up
    node_end: Optional[int] = field(default=-1)  # left-down
    distance_end: Optional[int] = field(default=0)  # left-down


pos = nx.get_node_attributes(scribe, "pos_bitmap")
components = []
for n in range(scribe.number_of_nodes()):
    # fill
    seed = pos[n]
    ccv = gray.copy()
    cv.floodFill(ccv, None, seed, RASMVAL, loDiff=(5), upDiff=(5))
    _, ccv = cv.threshold(ccv, 100, RASMVAL, cv.THRESH_BINARY)
    mu = cv.moments(ccv)
    if mu["m00"] > pow(SLIC_SPACE, 2) * PHI:
        mc = (int(mu["m10"] / (mu["m00"])), int(mu["m01"] / (mu["m00"])))
        area = mu["m00"]
        pd = pdistance(seed, mc)
        node_start = n
        box = cv.boundingRect(ccv)
        # append keypoint if the component already exists
        found = 0
        for i in range(len(components)):
            if components[i].centroid == mc:
                components[i].nodes.append(n)
                # calculate the distance
                tvane = freeman(seed[0] - mc[0], mc[1] - seed[1])
                # if seed[0]>mc[0] and pd>components[i].distance_start and (tvane==2 or tvane==4): # potential node_start for long rasm
                if (
                    seed[0] > mc[0] and pd > components[i].distance_start
                ):  # potential node_start
                    components[i].distance_start = pd
                    components[i].node_start = n
                elif (
                    seed[0] < mc[0] and pd > components[i].distance_end
                ):  # potential node_end
                    components[i].distance_end = pd
                    components[i].node_end = n
                found = 1
                # print(f'old node[{n}] with component[{i}] at {mc} from {components[i].centroid} distance: {pd})')
                break
        if found == 0:
            components.append(ConnectedComponents(box, mc))
            idx = len(components) - 1
            components[idx].nodes.append(n)
            components[idx].mat = ccv.copy()
            components[idx].area = int(mu["m00"] / THREVAL)
            if seed[0] > mc[0]:
                components[idx].node_start = n
                components[idx].distance_start = pd
            else:
                components[idx].node_end = n
                components[idx].distance_end = pd
            # print(f'new node[{n}] with component[{idx}] at {mc} from {components[idx].centroid} distance: {pd})')


components = sorted(components, key=lambda x: x.centroid[0], reverse=True)
# for n in len(components):
#     for i in components[n].nodes:
#         distance= pdistance(components[n].centroid, pos[i])
#         print(f'{i}: {distance}')

# drawing the starting node (bitmap level)
disp = cv.cvtColor(gray, cv.COLOR_GRAY2BGR)
for n in range(len(components)):
    # print(f'{n} at {components[n].centroid} size {components[n].area}')
    # draw green line for rasm at edges, color the rasm brighter
    if components[n].area > 4 * PHI * pow(SLIC_SPACE, 2):
        disp = cv.bitwise_or(disp, cv.cvtColor(components[n].mat, cv.COLOR_GRAY2BGR))
        seed = components[n].centroid
        cv.circle(disp, seed, 2, (0, 0, 120), -1)
        if components[n].node_start != -1:
            cv.circle(disp, pos[components[n].node_start], 2, (0, 120, 0), -1)
        if components[n].node_end != -1:
            cv.circle(disp, pos[components[n].node_end], 2, (120, 0, 0), -1)
        r = components[n].rect[0] + int(components[n].rect[2])
        l = components[n].rect[0]
        if l < width and r < width:  # did we ever went beyond the frame?
            for j1 in range(int(SLIC_SPACE * PHI), height - int(SLIC_SPACE * PHI)):
                disp[j1, r, 1] = 120
            for j1 in range(
                int(SLIC_SPACE * pow(PHI, 3)), height - int(SLIC_SPACE * pow(PHI, 3))
            ):
                disp[j1, l, 1] = 120
    else:
        m = components[n].centroid[1]
        i = components[n].centroid[0]
        # draw blue line for shakil 'connection'
        for j2 in range(
            int(m - (2 * SLIC_SPACE * PHI)), int(m + (2 * SLIC_SPACE * PHI))
        ):
            if j2 < height and j2 > 0:
                disp[j2, i, 1] = RASMVAL / 2

draw(disp)


### several profiles to detect inter-hurf transition
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

STRAYPIXEL = 4

# vertical projection histogram
vertical_projection = np.zeros((len(components), gray.shape[1]))
for n in range(len(components)):
    vertical_projection[n] = np.sum(components[n].mat, axis=0) / RASMVAL

# Plot each component histogram with different colors
colors = cm.nipy_spectral(np.linspace(0, 1, vertical_projection.shape[0]))
plt.figure(figsize=(14, 2))
for n, row in enumerate(vertical_projection):
    # Only plot values >= 4
    x = np.arange(len(row))  # X-axis positions for each row element
    y = np.where(row >= STRAYPIXEL, row, np.nan)  # Set values <4 to NaN for skipping
    plt.plot(x, y, color=colors[n], label=f"rasm {i+1}")

# plt.xlabel("Column Index")
plt.ylabel("number of pixels")
plt.title("vertical projection for each connected component (rasm)")
plt.show()

# countour of the stroke
edge_top = np.zeros((len(components), gray.shape[1]))
edge_bot = np.zeros((len(components), gray.shape[1]))

for n in range(len(components)):
    non_zero_cols = np.unique(np.where(components[n].mat[:, :] == RASMVAL)[1])
    for j in enumerate(non_zero_cols):
        non_zero_row = np.where(components[n].mat[:, j[1]])[0]
        edge_top[n][j[1]] = -np.min(non_zero_row)
        edge_bot[n][j[1]] = -np.max(non_zero_row)
        edge_width = edge_top[n][j[1]] - edge_bot[n][j[1]]

plt.figure(figsize=(14, 2))
for n, row in enumerate(edge_top):
    x = np.arange(len(row))  # X-axis positions for each row element
    y_top = np.where(edge_top[n] < -STRAYPIXEL, edge_top[n], np.nan)
    y_bot = np.where(edge_bot[n] < -STRAYPIXEL, edge_bot[n], np.nan)
    plt.plot(x, y_top, color=colors[n], label=f"rasm {n} top")
    plt.plot(x, y_bot, color=colors[n] * 0.8, label=f"rasm {n} bot")
plt.ylabel("contour position")
plt.title("top-bottom contour for each connected component (rasm)")
plt.show()


# TODO: SLANTED HISTOGRAM