704 KiB
704 KiB
In [47]:
import cv2
import numpy as np
import matplotlib.pyplot as plt
import seaborn as snsIn [30]:
threshold = 5In [ ]:
bkg = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg1_50_50.jpg")
curr_img = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/1002_0_1_50_50_curr_image.tiff")In [139]:
def flatfield_correction(raw_image, flat_image, dark_image=None):
"""
Apply flat-field correction to remove uneven illumination
Args:
raw_image: Your actual image with crystal
flat_image: Image of uniform illumination (no sample)
dark_image: Dark frame (camera with no light), optional
"""
# Convert to float for calculations
raw = raw_image.astype(np.float32)
flat = flat_image.astype(np.float32)
if dark_image is not None:
dark = dark_image.astype(np.float32)
# Subtract dark frame from both
raw_corrected = raw - dark
flat_corrected = flat - dark
else:
raw_corrected = raw
flat_corrected = flat
# Avoid division by zero
flat_corrected[flat_corrected == 0] = 1
# Apply correction
mean_flat = np.mean(flat_corrected)
corrected = (raw_corrected / flat_corrected) * mean_flat
# Convert back to original dtype
return np.clip(corrected, 0, 255).astype(raw_image.dtype)In [154]:
def create_ring_mask(image, center=None, inner_radius=50, outer_radius=100):
"""
Create a mask to exclude a ring region
Args:
image_shape: (height, width) of your image
center: (cx, cy) center of ring, None for image center
inner_radius: Inner radius of ring to mask out
outer_radius: Outer radius of ring to mask out
Returns:
mask: Binary mask (0 = exclude ring, 255 = use for correction)
"""
if len(image.shape) == 3:
h, w, _ = image.shape
else:
h, w = image.shape
if center is None:
center = (w // 2, h // 2)
# Create coordinate arrays
cx, cy = center
w = int(w)
h = int(h)
x = np.arange(w)
y = np.arange(h)
X, Y = np.meshgrid(x, y)
# Calculate distances from center
distances = np.sqrt((X - cx)**2 + (Y - cy)**2)
# Create mask: exclude the ring region
mask = np.ones((h, w), dtype=np.uint8) * 255
ring_region = (distances >= inner_radius) & (distances <= outer_radius)
mask[ring_region] = 0
return mask
def flatfield_correction_with_mask(raw_image, flat_image, mask, dark_image=None):
"""
Apply flat-field correction with mask to exclude certain regions
Args:
raw_image: Your camera image
flat_image: Flat-field reference
mask: Binary mask (0 = ignore, 255 = use)
dark_image: Dark frame (optional)
"""
# Convert to float
raw = raw_image.astype(np.float32)
flat = flat_image.astype(np.float32)
if dark_image is not None:
dark = dark_image.astype(np.float32)
raw_corrected = raw - dark
flat_corrected = flat - dark
else:
raw_corrected = raw
flat_corrected = flat
# Apply mask to flat field - only use unmasked regions for correction
mask_float = mask.astype(np.float32) / 255.0
# Calculate mean only from unmasked regions
if len(flat_corrected.shape) == 3:
# Color image
mean_flat = []
for channel in range(flat_corrected.shape[2]):
masked_flat = flat_corrected[:,:,channel] * mask_float
valid_pixels = masked_flat[mask > 0]
mean_flat.append(np.mean(valid_pixels) if len(valid_pixels) > 0 else 1.0)
mean_flat = np.array(mean_flat).reshape(1, 1, -1)
else:
# Grayscale
masked_flat = flat_corrected * mask_float
valid_pixels = masked_flat[mask > 0]
mean_flat = np.mean(valid_pixels) if len(valid_pixels) > 0 else 1.0
# Avoid division by zero
flat_corrected[flat_corrected <= 0] = 1
# Apply correction
corrected = (raw_corrected / flat_corrected) * mean_flat
# In masked regions, keep original values
mask_3d = mask_float
if len(corrected.shape) == 3:
mask_3d = np.stack([mask_float] * corrected.shape[2], axis=2)
# Blend: use corrected where mask=1, original where mask=0
final_result = corrected * mask_3d + raw_corrected * (1 - mask_3d)
return np.clip(final_result, 0, 255).astype(np.uint8)
def create_interactive_ring_mask(image):
"""
Interactively create a ring mask by clicking on the image
Click center, then inner radius point, then outer radius point
"""
points = []
def mouse_callback(event, x, y, flags, param):
if event == cv2.EVENT_LBUTTONDOWN:
points.append((x, y))
cv2.circle(image, (x, y), 3, (0, 255, 0), -1)
cv2.imshow('Select Ring', image)
if len(points) == 1:
print("Click on inner radius of ring")
elif len(points) == 2:
print("Click on outer radius of ring")
cv2.imshow('Select Ring', image)
print("Click on center of ring")
cv2.setMouseCallback('Select Ring', mouse_callback)
cv2.waitKey(0)
cv2.destroyAllWindows()
if len(points) == 3:
center = points[0]
inner_radius = np.sqrt((points[1][0] - center[0])**2 + (points[1][1] - center[1])**2)
outer_radius = np.sqrt((points[2][0] - center[0])**2 + (points[2][1] - center[1])**2)
mask = create_ring_mask(image.shape, center, inner_radius, outer_radius)
return mask
else:
print("Need 3 points!")
return None
In [354]:
min_contour_area = 4
mask = create_ring_mask(curr_img, inner_radius=1000,outer_radius=2000)
diff_image = cv2.absdiff(curr_img, bkg)
diff_image = cv2.cvtColor(diff_image, cv2.COLOR_RGB2GRAY)
thresh_value, binary_image = cv2.threshold(diff_image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
thresh_value, thresh_normal = cv2.threshold(diff_image, 50, 255, cv2.THRESH_BINARY)
ff = flatfield_correction(curr_img, bkg, dark_image=None)
ff = cv2.cvtColor(ff, cv2.COLOR_RGB2GRAY)
ff_masked = flatfield_correction_with_mask(curr_img, bkg, mask, dark_image=None)
ff_masked = cv2.cvtColor(ff_masked, cv2.COLOR_RGB2GRAY)
thresh_value, thresh_mask = cv2.threshold(ff_masked, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
kernel = np.ones((3,3), np.uint8)
thresh = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel)
thresh = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel)
contours, _ = cv2.findContours(thresh_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
valid_contours = [cnt for cnt in contours if cv2.contourArea(cnt) > min_contour_area]
largest_contour = max(valid_contours, key=cv2.contourArea)
# Step 6: Find leftmost point
leftmost_point = tuple(largest_contour[largest_contour[:,:,0].argmin()][0])
print(leftmost_point)
leftmost_point = None
for contour in contours:
# Loop through each point in the current contour
for point in contour:
x, y = point[0] # Point is a nested array [ [x, y] ]
# Check if this is the leftmost point
if leftmost_point is None or x < leftmost_point[0]:
leftmost_point = (x, y)
plot_contour = contour
print(leftmost_point)
c_img = cv2.drawContours(thresh, contours, -1, (0, 255, 0), 1)
plt.subplot(2, 2, 1)
plt.imshow(diff_image)
plt.subplot(2, 2, 2)
plt.imshow(thresh_normal)
#plt.subplot(2, 2, 4)
#plt.imshow(ff_masked)
result_image = cv2.cvtColor(thresh, cv2.COLOR_GRAY2BGR) if len(thresh.shape) == 2 else thresh.copy()
# Draw contour
result_image = cv2.drawContours(result_image, [largest_contour], -1, (0, 255, 0), 2)
# Mark leftmost point
result_image = cv2.circle(result_image, leftmost_point, 100, (0, 0, 255), -1)
result_image = cv2.putText(result_image, f'Leftmost: {leftmost_point}',
(leftmost_point[0] + 10, leftmost_point[1] - 10),
cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 0, 255), 2)
print(f'Leftmost point: {leftmost_point}')
plt.subplot(2, 2, 3)
plt.imshow(binary_image)
Out [354]:
(np.int32(24), np.int32(979)) (np.int32(24), np.int32(979)) Leftmost point: (np.int32(24), np.int32(979))
<matplotlib.image.AxesImage at 0x7f97e92802c0>
In [350]:
# Load image and convert
img = cv2.imread('/home/leonarski_f/aaredaq/daq/src/aaredaq/991_0_1_50_50_curr_image.tiff', cv2.IMREAD_COLOR)
#img = cv2.imread('/home/leonarski_f/aaredaq/daq/src/aaredaq/0_280_100_50_curr_image_colour.jpg', cv2.IMREAD_COLOR)
bkg = cv2.imread('/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg1_50_50.jpg')
# Convert to grayscale & enhance
gray_bkg = cv2.cvtColor(bkg, cv2.COLOR_BGR2GRAY)
gray = cv2.subtract(gray, gray_bkg)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
gray = cv2.GaussianBlur(gray, (5, 5), 0)
gray = cv2.equalizeHist(gray)
# Mask out camera-case ring
#gray = cv2.subtract(gray, gray_bkg)
mask = np.zeros_like(gray, dtype=np.uint8)
h, w = gray.shape
radius = min(h, w)//2 - 20
center = (w//2, h//2)
norm=gray
print(radius)
#norm = flatfield_correction_with_mask(gray, gray_bkg, mask)
#radius = 500#radius#-factor
print(gray.shape)
# # gray = cv2.bitwise_and(gray, gray, mask=mask)
# # x, y = center[0] - radius, center[1] - radius
# gray_cropped = gray[y:y+2*radius, x:x+2*radius]
# #gray_bkg_cropped = gray[y:y+2*radius, x:x+2*radius]
# #cv2.circle(mask, center, radius, 93, -1)
# print(gray_cropped.shape)
#
# #masked = cv2.bitwise_and(gray, gray, mask=mask)
# #norm = gray_cropped
# masked=gray_cropped
# background= cv2.medianBlur(masked, 51)
#
# norm = cv2.subtract(masked, background)
#norm = cv2.absdiff(masked,gray_bkg)
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (25,25))
tophat = cv2.morphologyEx(norm, cv2.MORPH_TOPHAT, kernel)
enhanced = cv2.addWeighted(norm, 0.7, tophat, 0.3, 0)
# Adaptive thresholding
thresh = cv2.adaptiveThreshold(enhanced, 255, cv2.ADAPTIVE_THRESH_MEAN_C,
cv2.THRESH_BINARY_INV, 41, 2)
# Find contours
contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
best_point = None
min_x = w + 1
print(len(contours))
for cnt in contours:
x,y,wc,hc = cv2.boundingRect(cnt)
# Skip very large contours (likely mask boundary)
if cv2.contourArea(cnt) > 0.2 * h * w:
print("Skipping contour with area", cv2.contourArea(cnt))
continue
# Compute distance of contour centroid from image center
M = cv2.moments(cnt)
if M["m00"] != 0:
cx = int(M["m10"]/M["m00"])
cy = int(M["m01"]/M["m00"])
dist = np.sqrt((cx-center[0])**2 + (cy-center[1])**2)
# If centroid is near mask radius, skip (mask edge)
#if abs(dist - radius) < 60:
# print("Skipping contour near mask radius")
# continue
# Skip vertical pin
if hc > 3*wc:
print("Skipping vertical pin")
continue
leftmost = tuple(cnt[cnt[:, :, 0].argmin()][0])
if leftmost[0] < min_x:
min_x = leftmost[0]
best_point = leftmost
# Draw result
if best_point is not None:
print(f"best_point = {best_point}")
#cv2.circle(thresh, best_point, 100, (0, 0, 255), -1)
print(radius)
corrected_point = (best_point[0]+factor, best_point[1] +factor)
print(f"corrected_point {corrected_point}")
cv2.circle(img, corrected_point, 100, (0, 0, 255), -1)
print("Loop tip (leftmost point):", best_point)
else:
print("No suitable contour found")
plt.figure(figsize=(10,8))
plt.figure(figsize=(14,10))
plt.subplot(2,4,1); plt.title("Original"); plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))#; plt.axis("off")
plt.subplot(2,4,2); plt.title("Cropped"); plt.imshow(gray_cropped, cmap="gray"); plt.axis("off")
plt.subplot(2,4,3); plt.title("Background"); plt.imshow(norm, cmap="gray"); plt.axis("off")
plt.subplot(2,4,4); plt.title("Top-hat"); plt.imshow(tophat, cmap="gray"); plt.axis("off")
plt.subplot(2,4,5); plt.title("Enhanced"); plt.imshow(enhanced, cmap="gray"); plt.axis("off")
plt.subplot(2,4,6); plt.title("Threshold"); plt.imshow(thresh, cmap="gray")#; plt.axis("off")
plt.subplot(2,4,7); plt.title("image with contours"); #plt.imshow(img)#; plt.axis("off")
contour_vis = cv2.cvtColor(masked, cv2.COLOR_GRAY2BGR)
cv2.drawContours(contour_vis, contours, -1, (0,255,0), 1)
cv2.circle(contour_vis, best_point, 100, (0, 0, 255), -1)
plt.imshow(contour_vis)
plt.show()
cv2.imwrite('masked.tiff', enhanced)
cv2.imwrite('thresh.tiff', thresh)Out [350]:
1003 (2046, 2046) 10291 Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping contour with area 2054203.5 Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin Skipping vertical pin best_point = (np.int32(0), np.int32(1421)) 1003 corrected_point (np.int32(250), np.int32(1671)) Loop tip (leftmost point): (np.int32(0), np.int32(1421))
<Figure size 1000x800 with 0 Axes>
True
In [275]:
x#bkg_old = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_020925.jpg")
#bkg_new = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_03092025.jpg")
bkg_old = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg280_50_50_old.jpg")
bkg_new = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg280_50_50.jpg")
#$bkg_newest = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg280.jpg")
grey_old = cv2.cvtColor(bkg_old, cv2.COLOR_RGB2GRAY)
grey_new = cv2.cvtColor(bkg_new, cv2.COLOR_RGB2GRAY)
#grey_newest = cv2.cvtColor(bkg_newest, cv2.COLOR_RGB2GRAY)
diff_image = cv2.absdiff(grey_old, grey_new)
_, thresh_diff = cv2.threshold(diff_image, threshold, 255, cv2.THRESH_BINARY)
# Simple subtraction (may result in negative values)
diff_simple = cv2.subtract(bkg_new, bkg_old)
# Weighted difference for motion detection
diff_weighted = cv2.addWeighted(bkg_new, 0.5, bkg_old, -0.5, 128)
# Using numpy for difference
diff_numpy = np.abs(bkg_new.astype(np.int16) - bkg_old.astype(np.int16)).astype(np.uint8)
plt.subplot(1,2,1)
plt.imshow(bkg_old)
plt.subplot(1,2,2)
plt.imshow(bkg_new)Out [275]:
<matplotlib.image.AxesImage at 0x7f97bc9dc7a0>
In [196]:
bkg_old_np = bkg_old.astype(np.float32)
bkg_old_np_normalized = bkg_old.astype(np.float32) / 255.0
bkg_new_np = bkg_old.astype(np.float32)
bkg_new_np_normalized = bkg_old.astype(np.float32) / 255.0
diff_np=bkg_old_np-bkg_new_np
diff_np_normalized=bkg_old_np_normalized-bkg_new_np_normalizedIn [200]:
print("mean: ", np.mean(grey_old.astype(np.float32)-grey_new.astype(np.float32)))
print("std: ", np.std(grey_old.astype(np.float32)-grey_new.astype(np.float32)))
diff_abs = np.abs(grey_old.astype(np.float32)-grey_new.astype(np.float32))
print(np.percentile(diff_abs, 99))mean: -6.184164 std: 5.93317 21.0
In [ ]:
threshold_old = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_020925.jpg")
threshold_new = cv2.imread("/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_03092025.jpg")In [189]:
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
# Example showing the difference
image = np.random.rand(50, 50) * 100
# Add a bright spot
image[20:30, 35:45] = 200
# Without thresholding - all pixels contribute
com_no_threshold = ndimage.center_of_mass(image)
print(f"No threshold COM: ({com_no_threshold[1]:.2f}, {com_no_threshold[0]:.2f})")
# With thresholding - only bright pixels contribute
threshold = 150
thresholded = image.copy()
thresholded[thresholded < threshold] = 0
com_with_threshold = ndimage.center_of_mass(thresholded)
print(f"With threshold COM: ({com_with_threshold[1]:.2f}, {com_with_threshold[0]:.2f})")
No threshold COM: (26.11, 24.33) With threshold COM: (39.50, 24.50)
In [190]:
def center_of_mass_with_threshold(image, threshold):
"""
Apply threshold before center of mass calculation
"""
# Create binary mask
binary = (image >= threshold).astype(float)
# Calculate center of mass on thresholded image
com = ndimage.center_of_mass(binary)
# Return as (x, y) coordinates
return (com[1], com[0]) if com[0] is not np.nan else None
In [191]:
def weighted_center_of_mass_threshold(image, threshold):
"""
Keep original intensities above threshold, zero below
"""
weighted_image = image.copy()
weighted_image[weighted_image < threshold] = 0
com = ndimage.center_of_mass(weighted_image)
return (com[1], com[0]) if com[0] is not np.nan else None
In [ ]: