m1_4kandClaude Sonnet 4.5
58485d8844
Implements Kovesi (1999) phase congruency algorithm using log-Gabor filters across multiple scales and orientations for illumination-invariant and contrast-invariant edge detection. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
235 lines
8.3 KiB
Python
235 lines
8.3 KiB
Python
import torch
|
||
import numpy as np
|
||
from PIL import Image, ImageOps
|
||
import math
|
||
|
||
|
||
class PhaseCongruencyEdge:
|
||
"""
|
||
ComfyUI Custom Node: Phase Congruency Edge Detection
|
||
|
||
Detects edges using phase congruency — an illumination/contrast invariant
|
||
feature detector based on the coherence of Fourier phase components
|
||
across multiple scales and orientations (Kovesi, 1999).
|
||
"""
|
||
|
||
@classmethod
|
||
def INPUT_TYPES(s):
|
||
return {
|
||
"required": {
|
||
"image": ("IMAGE",),
|
||
"nscale": ("INT", {
|
||
"default": 4,
|
||
"min": 2,
|
||
"max": 8,
|
||
"step": 1,
|
||
"tooltip": "Number of wavelet scales"
|
||
}),
|
||
"norient": ("INT", {
|
||
"default": 6,
|
||
"min": 2,
|
||
"max": 8,
|
||
"step": 1,
|
||
"tooltip": "Number of filter orientations"
|
||
}),
|
||
"min_wavelength": ("INT", {
|
||
"default": 3,
|
||
"min": 2,
|
||
"max": 20,
|
||
"step": 1,
|
||
"tooltip": "Minimum wavelength of log-Gabor filter"
|
||
}),
|
||
"mult": ("FLOAT", {
|
||
"default": 2.1,
|
||
"min": 1.5,
|
||
"max": 4.0,
|
||
"step": 0.1,
|
||
"tooltip": "Scaling factor between successive filter scales"
|
||
}),
|
||
"sigma_on_f": ("FLOAT", {
|
||
"default": 0.55,
|
||
"min": 0.1,
|
||
"max": 1.0,
|
||
"step": 0.05,
|
||
"tooltip": "Bandwidth of log-Gabor filter"
|
||
}),
|
||
"k": ("FLOAT", {
|
||
"default": 2.0,
|
||
"min": 0.5,
|
||
"max": 10.0,
|
||
"step": 0.5,
|
||
"tooltip": "Noise threshold (std devs above mean)"
|
||
}),
|
||
"cutoff": ("FLOAT", {
|
||
"default": 0.5,
|
||
"min": 0.1,
|
||
"max": 0.9,
|
||
"step": 0.05,
|
||
"tooltip": "Fractional measure of frequency spread below which phase congruency values get penalized"
|
||
}),
|
||
},
|
||
}
|
||
|
||
RETURN_TYPES = ("IMAGE",)
|
||
RETURN_NAMES = ("phase_congruency_edge",)
|
||
FUNCTION = "compute_phase_congruency"
|
||
CATEGORY = "image/filters"
|
||
|
||
def compute_phase_congruency(self, image, nscale, norient, min_wavelength,
|
||
mult, sigma_on_f, k, cutoff):
|
||
# image: [B, H, W, C] torch tensor, values in [0, 1]
|
||
batch_results = []
|
||
|
||
for b in range(image.shape[0]):
|
||
img_np = image[b].cpu().numpy() # [H, W, C]
|
||
|
||
# Convert to grayscale float in [0, 1]
|
||
if img_np.shape[2] == 1:
|
||
gray = img_np[:, :, 0].astype(np.float64)
|
||
else:
|
||
# Luminance conversion
|
||
gray = (0.299 * img_np[:, :, 0] +
|
||
0.587 * img_np[:, :, 1] +
|
||
0.114 * img_np[:, :, 2]).astype(np.float64)
|
||
|
||
pc = self._phase_congruency(
|
||
gray, nscale, norient, min_wavelength,
|
||
mult, sigma_on_f, k, cutoff
|
||
)
|
||
|
||
# Normalize to [0, 1]
|
||
pc_min, pc_max = pc.min(), pc.max()
|
||
if pc_max - pc_min > 1e-10:
|
||
pc = (pc - pc_min) / (pc_max - pc_min)
|
||
else:
|
||
pc = np.zeros_like(pc)
|
||
|
||
# Convert to [H, W, 3] RGB tensor
|
||
pc_rgb = np.stack([pc, pc, pc], axis=-1).astype(np.float32)
|
||
batch_results.append(torch.from_numpy(pc_rgb))
|
||
|
||
output = torch.stack(batch_results, dim=0) # [B, H, W, 3]
|
||
return (output,)
|
||
|
||
def _phase_congruency(self, img, nscale, norient, min_wavelength,
|
||
mult, sigma_on_f, k, cutoff):
|
||
"""
|
||
Computes phase congruency for a 2D grayscale image [0,1].
|
||
Based on Kovesi's algorithm using log-Gabor filters.
|
||
|
||
For each orientation:
|
||
1. Apply log-Gabor filters at each scale
|
||
2. Accumulate complex responses (energy vector) and amplitudes
|
||
3. Estimate noise threshold using Rayleigh statistics
|
||
4. PC_orient = max(|energy| - T, 0) / (sum_amplitude + epsilon)
|
||
Final PC = mean over all orientations.
|
||
"""
|
||
rows, cols = img.shape
|
||
epsilon = 1e-4 # avoid division by zero
|
||
|
||
# Frequency coordinate grids (centered, range -0.5..0.5)
|
||
cx = (np.arange(cols) - cols // 2) / cols
|
||
cy = (np.arange(rows) - rows // 2) / rows
|
||
x, y = np.meshgrid(cx, cy)
|
||
|
||
radius = np.sqrt(x**2 + y**2)
|
||
radius[rows // 2, cols // 2] = 1.0 # avoid log(0) at DC
|
||
|
||
theta = np.arctan2(-y, x) # orientation of each frequency component
|
||
|
||
# Butterworth low-pass filter to suppress high-frequency noise
|
||
lp = self._lowpass_filter(rows, cols, 0.45, 15)
|
||
|
||
# FFT of input image (shifted so DC is at center)
|
||
img_fft = np.fft.fftshift(np.fft.fft2(img))
|
||
|
||
# Accumulate per-orientation PC maps, then average
|
||
pc_sum = np.zeros((rows, cols))
|
||
|
||
for orient in range(norient):
|
||
angle = orient * np.pi / norient # filter orientation
|
||
|
||
# Angular spread function (raised cosine)
|
||
ds = np.sin(theta) * np.cos(angle) - np.cos(theta) * np.sin(angle)
|
||
dc = np.cos(theta) * np.cos(angle) + np.sin(theta) * np.sin(angle)
|
||
dtheta = np.abs(np.arctan2(ds, dc))
|
||
dtheta = np.minimum(dtheta * norient / 2.0, np.pi)
|
||
spread = (np.cos(dtheta) + 1.0) / 2.0
|
||
|
||
sum_an = np.zeros((rows, cols))
|
||
energy_real = np.zeros((rows, cols))
|
||
energy_imag = np.zeros((rows, cols))
|
||
an_finest = None # amplitude at the finest (smallest) scale
|
||
|
||
for scale in range(nscale):
|
||
wavelength = min_wavelength * (mult ** scale)
|
||
fo = 1.0 / wavelength
|
||
|
||
# Log-Gabor radial component
|
||
log_gabor = np.exp(
|
||
-(np.log(radius / fo)) ** 2 /
|
||
(2.0 * (np.log(sigma_on_f) ** 2))
|
||
)
|
||
log_gabor[rows // 2, cols // 2] = 0.0 # remove DC
|
||
|
||
# Combined filter = radial × angular × low-pass
|
||
filt = log_gabor * spread * lp
|
||
|
||
# Complex filter response
|
||
resp = np.fft.ifft2(np.fft.ifftshift(img_fft * filt))
|
||
re, im = resp.real, resp.imag
|
||
an = np.sqrt(re**2 + im**2)
|
||
|
||
sum_an += an
|
||
energy_real += re
|
||
energy_imag += im
|
||
|
||
if scale == 0:
|
||
an_finest = an # save finest scale amplitude for noise est.
|
||
|
||
# Rayleigh noise threshold from finest-scale amplitude statistics
|
||
# E[amplitude] for Rayleigh: mean = tau * sqrt(pi/2)
|
||
tau = np.sqrt(np.mean(an_finest**2) / 2.0 + epsilon)
|
||
mean_noise = tau * math.sqrt(math.pi / 2.0)
|
||
noise_std = tau * math.sqrt((4.0 - math.pi) / 2.0)
|
||
T = mean_noise + k * noise_std
|
||
|
||
# Phase congruency for this orientation
|
||
energy = np.sqrt(energy_real**2 + energy_imag**2)
|
||
pc_orient = np.maximum(energy - T, 0.0) / (sum_an + epsilon)
|
||
pc_sum += pc_orient
|
||
|
||
# Average over all orientations
|
||
pc = pc_sum / norient
|
||
return pc
|
||
|
||
def _lowpass_filter(self, rows, cols, cutoff, n):
|
||
"""
|
||
Butterworth low-pass filter in the frequency domain.
|
||
cutoff: cutoff frequency as fraction of Nyquist (0..0.5)
|
||
n: filter order
|
||
"""
|
||
if cols % 2 == 0:
|
||
cx = np.arange(-cols // 2, cols // 2) / cols
|
||
else:
|
||
cx = np.arange(-(cols - 1) // 2, (cols - 1) // 2 + 1) / cols
|
||
|
||
if rows % 2 == 0:
|
||
cy = np.arange(-rows // 2, rows // 2) / rows
|
||
else:
|
||
cy = np.arange(-(rows - 1) // 2, (rows - 1) // 2 + 1) / rows
|
||
|
||
x, y = np.meshgrid(cx, cy)
|
||
radius = np.sqrt(x**2 + y**2)
|
||
lp = 1.0 / (1.0 + (radius / cutoff) ** (2 * n))
|
||
return lp
|
||
|
||
|
||
NODE_CLASS_MAPPINGS = {
|
||
"PhaseCongruencyEdge": PhaseCongruencyEdge
|
||
}
|
||
|
||
NODE_DISPLAY_NAME_MAPPINGS = {
|
||
"PhaseCongruencyEdge": "Phase Congruency Edge"
|
||
}
|