Here is the code:
Code for Edge-Based Resolution Sampling In the Middle, Taylor 4/29/25+
#Extending the edge to the end of the image.
import numpy as np
import matplotlib.pyplot as plt
from skimage import io, feature
from skimage.transform import probabilistic_hough_line
from skimage.color import rgb2gray
from scipy.ndimage import gaussian_filter1d
from shapely.geometry import LineString, box
---------------------
CONFIGURATION
---------------------
FRACTION_RANGE = (0.4, 0.6) # e.g. (0.4, 0.6) = center portion of the edge
NUM_SAMPLES = 3 # number of profiles to extract in that range
---------------------
LOGIC
---------------------
def load_image(path):
img = io.imread(path)
if img.ndim == 3:
img = rgb2gray(img)
return img
def detect_edge_line(image):
edges = feature.canny(image, sigma=2)
lines = probabilistic_hough_line(edges, threshold=10, line_length=50, line_gap=10)
if not lines:
raise ValueError(“No edge line detected.”)
# Use the longest line detected
lines.sort(key=lambda l: np.linalg.norm(np.subtract(l[0], l[1])), reverse=True)
(x1, y1), (x2, y2) = lines[0]
# Extend line to image bounds using Shapely
dx, dy = x2 - x1, y2 - y1
line = LineString([(x1 - 10000 * dx, y1 - 10000 * dy), (x2 + 10000 * dx, y2 + 10000 * dy)])
image_box = box(0, 0, image.shape[1], image.shape[0])
cropped = line.intersection(image_box)
if cropped.is_empty or not isinstance(cropped, LineString):
raise ValueError("Extended line does not intersect image bounds.")
x1c, y1c, x2c, y2c = *cropped.coords[0], *cropped.coords[-1]
return ((int(x1c), int(y1c)), (int(x2c), int(y2c)))
def get_perpendicular_lines(line, length=50):
(x1, y1), (x2, y2) = line
dx, dy = x2 - x1, y2 - y1
normal = np.array([-dy, dx], dtype=np.float64)
normal /= np.linalg.norm(normal)
sample_fracs = np.linspace(FRACTION_RANGE[0], FRACTION_RANGE[1], NUM_SAMPLES)
perp_lines =
for f in sample_fracs:
mx = x1 + dx * f
my = y1 + dy * f
start = (int(mx - normal[0] * length / 2), int(my - normal[1] * length / 2))
end = (int(mx + normal[0] * length / 2), int(my + normal[1] * length / 2))
perp_lines.append((start, end))
return perp_lines
def extract_profile(image, line):
(x0, y0), (x1, y1) = line
length = int(np.hypot(x1 - x0, y1 - y0))
x_vals = np.linspace(x0, x1, length)
y_vals = np.linspace(y0, y1, length)
x_vals = np.clip(x_vals.astype(int), 0, image.shape[1] - 1)
y_vals = np.clip(y_vals.astype(int), 0, image.shape[0] - 1)
profile = image[y_vals, x_vals]
return profile
def compute_lsf(esf):
if np.mean(esf[:5]) > np.mean(esf[-5:]):
esf = esf[::-1]
lsf = np.gradient(esf)
return gaussian_filter1d(lsf, sigma=1.0)
def compute_fwhm(lsf):
peak = np.max(lsf)
half_max = peak / 2
indices = np.where(lsf >= half_max)[0]
return indices[-1] - indices[0] if len(indices) > 1 else None
---------------------
VISUALIZATION
---------------------
def visualize_results(image, edge_line, perp_lines, profiles, lsfs, fwhms):
fig, axs = plt.subplots(2, 2, figsize=(12, 10))
axs[0, 0].imshow(image, cmap='gray')
axs[0, 0].set_title("Detected Edge and Sampled Lines")
axs[0, 0].plot([edge_line[0][0], edge_line[1][0]],
[edge_line[0][1], edge_line[1][1]], 'r-', linewidth=2, label='Edge')
for idx, line in enumerate(perp_lines):
axs[0, 0].plot([line[0][0], line[1][0]],
[line[0][1], line[1][1]], '--', linewidth=1.5, label=f'Perp {idx+1}')
axs[0, 0].legend()
axs[0, 1].set_title("Edge Spread Functions (ESF)")
for i, p in enumerate(profiles):
axs[0, 1].plot(p, label=f'ESF {i+1}')
axs[0, 1].legend()
axs[1, 0].set_title("Line Spread Functions (LSF)")
for i, lsf in enumerate(lsfs):
axs[1, 0].plot(lsf, label=f'LSF {i+1}')
axs[1, 0].legend()
axs[1, 1].axis('off')
text = "\n".join([f"FWHM {i+1}: {fwhm:.2f} px" if fwhm else f"FWHM {i+1}: N/A"
for i, fwhm in enumerate(fwhms)])
axs[1, 1].text(0.1, 0.6, text, fontsize=14)
plt.tight_layout()
plt.show()
---------------------
MAIN FUNCTION
---------------------
def main(image_path):
image = load_image(image_path)
edge_line = detect_edge_line(image)
perp_lines = get_perpendicular_lines(edge_line)
profiles = [extract_profile(image, line) for line in perp_lines]
lsfs = [compute_lsf(profile) for profile in profiles]
fwhms = [compute_fwhm(lsf) for lsf in lsfs]
visualize_results(image, edge_line, perp_lines, profiles, lsfs, fwhms)
---------------------
RUN
---------------------
if name == “main”:
main(“/content/green2_G_060425_1.png”) # Replace with actual image file path
