New routines for feature detection and analysis.

This commit is contained in:
Karthik 2025-09-12 16:46:46 +02:00
parent 22004110e0
commit 85ade404f2
7 changed files with 429 additions and 27 deletions

View File

@ -75,7 +75,8 @@ function [patchProps, patchCentroidsGlobal, imgCropped, xStart, yStart] = detect
%% --- Step 5: Detect lattice patches ---
opts = struct('sigmaSmall', params.dogGaussianSmallSigma, ...
'sigmaLarge', params.dogGaussianLargeSigma, ...
'minPeakProminence', params.minPeakProminence, ...
'adaptiveSensitivity', params.adaptiveSensitivity, ...
'adaptiveNeighborhoodSize', params.adaptiveNeighborhoodSize, ...
'minPatchArea', params.minimumPatchArea, ...
'minPeakFraction', params.minPeakFraction);
@ -95,20 +96,20 @@ function patchProps = detectPatchesCore(I, opts)
% thresholding, and connected component analysis.
% Returns a struct array with fields:
% Centroid, Area, Orientation, MajorAxisLength, MinorAxisLength
% Step 1: Difference-of-Gaussians
G1 = imgaussfilt(I, opts.sigmaSmall);
G2 = imgaussfilt(I, opts.sigmaLarge);
dogResponse = mat2gray(G1 - G2);
% Step 2: Threshold
th = max(graythresh(dogResponse), opts.minPeakProminence);
binaryMask = dogResponse > th;
% Step 2: Adaptive threshold
T = adaptthresh(dogResponse, opts.adaptiveSensitivity, 'NeighborhoodSize', opts.adaptiveNeighborhoodSize, 'Statistic','gaussian');
binaryMask = imbinarize(dogResponse, T);
% Step 3: Remove small specks, close small gaps
binaryMask = bwareaopen(binaryMask, 5);
binaryMask = imclose(binaryMask, strel('disk',1));
% Step 4: Connected components
CC = bwconncomp(binaryMask);
if CC.NumObjects == 0
@ -116,10 +117,10 @@ function patchProps = detectPatchesCore(I, opts)
'MajorAxisLength',[],'MinorAxisLength',[]);
return;
end
% Step 5: Extract patch properties
patchPropsAll = regionprops(CC, dogResponse, 'Centroid','Area','Orientation','MajorAxisLength','MinorAxisLength');
% Step 6: Filter by minimum area and adaptive per-patch peak fraction
maxDog = max(dogResponse(:)); % Global maximum DoG response
keepIdx = false(1,numel(patchPropsAll));

View File

@ -0,0 +1,117 @@
function objective = evaluateFeatureDetectionPipeline(x, img, gtMask, params, doPlot)
% evaluateFeatureDetectionPipeline
% Objective wrapper for Bayesian optimization.
% Runs the detection + shape extraction pipeline for a candidate
% parameter set and returns a score relative to ground truth.
%
% INPUTS:
% x - struct or table from bayesopt with tunable params
% img - 2D image (double or uint)
% gtMask - logical mask of ground truth region(s)
% params - base parameter struct (fixed params)
% doPlot - optional boolean flag to display results (default=false)
%
% OUTPUT:
% objective - scalar (to be MINIMIZED by bayesopt)
if nargin < 5
doPlot = false;
end
% --- Copy params and overwrite tunable fields ---
p = params;
% Clip / enforce constraints
p.dogGaussianSmallSigma = min(max(x.dogGaussianSmallSigma,0.01), x.dogGaussianLargeSigma-0.01);
p.dogGaussianLargeSigma = max(x.dogGaussianLargeSigma, p.dogGaussianSmallSigma+0.01);
p.adaptiveSensitivity = min(max(x.adaptiveSensitivity,0),1);
% Use nearestOdd helper for neighborhood size
p.adaptiveNeighborhoodSize = nearestOdd(x.adaptiveNeighborhoodSize, 3);
p.minPeakFraction = min(max(x.minPeakFraction,0),1);
p.minimumPatchArea = max(round(x.minimumPatchArea),1);
p.intensityThreshFraction = min(max(x.intensityThreshFraction,0),1);
p.edgeSigma = max(x.edgeSigma,0.01);
lowT = min(max(x.edgeThresholdLow,0),1);
highT = min(max(x.edgeThresholdHigh,0),1);
if highT <= lowT, highT = min(lowT+0.05,1); end
p.edgeThresholdLow = lowT;
p.edgeThresholdHigh = highT;
% --- Run detectPatches safely ---
try
[patchProps, ~, imgCropped, xStart, yStart] = Analyzer.detectPatches(img, p);
catch ME
warning('detectPatches failed: %s', ME.message);
objective = 0;
return;
end
% --- Run extractAndClassifyShapes safely ---
try
results = Analyzer.extractAndClassifyShapes(imgCropped, patchProps, xStart, yStart, p);
catch ME
warning('extractAndClassifyShapes failed: %s', ME.message);
objective = 0;
return;
end
% --- Build predicted mask ---
predictedMask = false(size(img));
for k = 1:numel(results)
BW = results(k).BW;
if isempty(BW), continue; end
[h,w] = size(BW);
yIdx = (yStart:yStart+h-1);
xIdx = (xStart:xStart+w-1);
predictedMask(yIdx, xIdx) = predictedMask(yIdx, xIdx) | BW;
end
% --- Compute Dice score ---
if any(gtMask(:))
score = dice(predictedMask, gtMask);
else
score = 0;
end
% --- Optional plotting ---
if doPlot
p.gtMask = gtMask; % overlay
Plotter.plotDetectedPatches(img, patchProps, results, p, xStart, yStart, 'Optimization Candidate');
drawnow;
end
% bayesopt minimizes return negative Dice
objective = -score;
end
function val = nearestOdd(x, minVal)
%NEARESTODD Round to the nearest odd integer, enforcing a minimum.
%
% val = NEARESTODD(x) rounds x to the nearest odd integer.
% val = NEARESTODD(x, minVal) enforces that the result is at least minVal.
if nargin < 2
minVal = 1;
end
% Round to nearest integer
n = round(x);
if mod(n,2) == 0
lowerOdd = n - 1;
upperOdd = n + 1;
% Choose whichever odd is closer to original x
if abs(lowerOdd - x) <= abs(upperOdd - x)
n = lowerOdd;
else
n = upperOdd;
end
end
% Enforce minimum
val = max(n, minVal);
end

View File

@ -36,7 +36,8 @@ function runInteractiveFeatureDetectorGUI(od_imgs, scan_parameter_values, file_l
'boundingBoxPadding',...
'dogGaussianSmallSigma',...
'dogGaussianLargeSigma',...
'minPeakProminence',...
'adaptiveSensitivity', ... % .
'adaptiveNeighborhoodSize', ... % Larger smoother masks, less sensitive to noise.
'minPeakFraction',...
'minimumPatchArea',...
'shapeMinArea', ...
@ -54,8 +55,9 @@ function runInteractiveFeatureDetectorGUI(od_imgs, scan_parameter_values, file_l
'Bounding Box Padding (px)',...
'DoG Small Kernel Sigma',...
'DoG Large Kernel Sigma',...
'Minimum Peak Prominence',...
'Minimum Peak Fraction',...
'Adaptive threshold sensitivity',...
'Window size for adaptive threshold', ...
'Minimum Patch Peak Fraction',...
'Minimum Patch Area (px)',...
'Minimum Shape Area (px)', ...
'Morphological Closing Radius (px)', ...
@ -72,8 +74,9 @@ function runInteractiveFeatureDetectorGUI(od_imgs, scan_parameter_values, file_l
'Number of pixels to pad around detected cloud bounding box', ...
'Sigma of small Gaussian in Difference-of-Gaussians (DoG)', ...
'Sigma of large Gaussian in DoG', ...
'Minimum prominence for DoG peak detection', ...
'Minimum fraction of max DoG response for adaptive thresholding', ...
'Global threshold that accounts for spatial variation of intensity', ...
'Local window size over which threshold is computed', ...
'Minimum fraction of max DoG response for detected patches', ...
'Minimum area (pixels) for detected patches', ...
'Minimum area (pixels) for internal shapes within patches', ...
'Radius (pixels) used for morphological closing to smooth shapes', ...

View File

@ -0,0 +1,110 @@
function plotDetectedPatches(img, patchProps, results, params, xStart, yStart, figTitle)
% plotDetectedPatches
% Plots the image, DoG-detected patch ellipses, and extracted shape boundaries
% with optional text overlay for Dice score and patch count.
%
% INPUTS:
% img - original 2D image
% patchProps - struct array from detectPatches
% results - struct array from extractAndClassifyShapes
% params - parameter struct
% xStart, yStart - offsets for cropped regions
% figTitle - title string (e.g., 'BayesOpt Candidate')
if nargin < 7
figTitle = '';
end
[Ny, Nx] = size(img);
dx = params.pixelSize / params.magnification;
dy = dx;
% Physical axes (μm)
xAxis = ((1:Nx)-(Nx+1)/2) * dx * 1e6;
yAxis = ((1:Ny)-(Ny+1)/2) * dy * 1e6;
%% --- Create Figure ---
% Try to find an existing figure by a unique tag
hFig = findobj('Type','figure','Tag','OptimizerViewer');
if isempty(hFig)
% If not found, create a new figure
hFig = figure('Name','Bayesian Optimization Live Viewer', ...
'NumberTitle','off', ...
'Position', [50 50 1000 800], ...
'KeyPressFcn',@keyPressCallback, ...
'Tag','OptimizerViewer'); % <-- unique tag
else
% If figure exists, bring it to front
figure(hFig);
clf;
end
imagesc(xAxis, yAxis, img);
axis(gca,'equal','tight');
colormap(Colormaps.inferno());
set(gca,'FontSize',14,'YDir','normal');
xlabel(gca,'x (\mum)','FontName', 'Bahnschrift', 'FontSize', 14, 'FontWeight', 'bold');
ylabel(gca,'y (\mum)','FontName', 'Bahnschrift', 'FontSize', 14, 'FontWeight', 'bold');
hold on;
% Draw patch ellipses
for k = 1:numel(patchProps)
a = patchProps(k).MajorAxisLength/2 * dx * 1e6;
b = patchProps(k).MinorAxisLength/2 * dy * 1e6;
phi = deg2rad(-patchProps(k).Orientation);
theta = linspace(0,2*pi,100);
R = [cos(phi) -sin(phi); sin(phi) cos(phi)];
ellipseXY = [a*cos(theta(:)) b*sin(theta(:))]*R';
cx_um = ((patchProps(k).Centroid(1)+xStart-1)-(Nx+1)/2)*dx*1e6;
cy_um = ((patchProps(k).Centroid(2)+yStart-1)-(Ny+1)/2)*dy*1e6;
ellipseXY(:,1) = ellipseXY(:,1) + cx_um;
ellipseXY(:,2) = ellipseXY(:,2) + cy_um;
plot(ellipseXY(:,1), ellipseXY(:,2),'g-','LineWidth',1);
end
% Overlay shapes
for k = 1:numel(results)
for n = 1:numel(results(k).boundaries)
bnd = results(k).boundaries{n};
bx = ((bnd(:,2)-(Nx+1)/2)) * dx * 1e6;
by = ((bnd(:,1)-(Ny+1)/2)) * dy * 1e6;
plot(bx, by, 'c-', 'LineWidth', 1.5);
end
end
% Compute Dice if gtMask is available
if isfield(params,'gtMask') && ~isempty(params.gtMask)
predictedMask = false(size(img));
for k = 1:numel(results)
BW = results(k).BW;
if isempty(BW), continue; end
[h,w] = size(BW);
yIdx = (yStart:yStart+h-1);
xIdx = (xStart:xStart+w-1);
predictedMask(yIdx, xIdx) = predictedMask(yIdx, xIdx) | BW;
end
diceScore = dice(predictedMask, params.gtMask);
% Convert logical mask to coordinates for overlay
[maskY, maskX] = find(params.gtMask);
maskX_um = (maskX - (Nx+1)/2) * dx * 1e6;
maskY_um = (maskY - (Ny+1)/2) * dy * 1e6;
% GT mask Overlay
overlayColor = [0.94 0.94 0.94]; % RGB from middle of coolwarm
scatter(maskX_um, maskY_um, 15, overlayColor, 'filled', ...
'MarkerFaceAlpha', 0.35, 'MarkerEdgeAlpha', 0);
else
diceScore = NaN;
end
% Add text overlay in top-right
txtStr = sprintf('Detected patches: %d\nDice score: %.3f', numel(patchProps), diceScore);
text(0.975, 0.975, txtStr, 'Units','normalized', 'Color','w', ...
'FontName', 'Bahnschrift', 'FontWeight','bold', 'FontSize',14, 'HorizontalAlignment','right', 'VerticalAlignment','top');
title(figTitle,'FontName', 'Bahnschrift', 'FontSize',16,'FontWeight','bold');
hold off;
end

View File

@ -93,21 +93,22 @@ end
%% Minimal pipeline: preprocess + detect patches + extract shapes
% ---------- user-tunable parameters ----------
params.backgroundDiskFraction = 1/8; % Fraction of image size used for morphological opening
params.backgroundDiskFraction = 0.1250; % Fraction of image size used for morphological opening
params.boundingBoxPadding = 12; % Padding around detected cloud
params.dogGaussianSmallSigma = 1.2; % Sigma for small Gaussian in Difference-of-Gaussians
params.dogGaussianLargeSigma = 3.8; % Sigma for large Gaussian in Difference-of-Gaussians
params.minPeakProminence = 0.02; % Threshold for DoG response
params.minPeakFraction = 0.80; % Fraction of max DoG response to reject weak patches
params.minimumPatchArea = 100; % Minimum area (pixels) for detected patches to be kept
params.dogGaussianSmallSigma = 1.1498; % Sigma for small Gaussian in Difference-of-Gaussians
params.dogGaussianLargeSigma = 3.9232; % Sigma for large Gaussian in Difference-of-Gaussians
params.adaptiveSensitivity = 0.4552; % Adaptive threshold sensitivity Higher more pixels marked foreground (lower threshold).
params.adaptiveNeighborhoodSize = 13; % Window size for adaptive threshold Defines the local window size over which threshold is computed. Larger smoother masks, less sensitive to noise.
params.minPeakFraction = 0.7190; % Fraction of max DoG response to reject weak patches
params.minimumPatchArea = 81; % Minimum area (pixels) for detected patches to be kept
params.shapeMinArea = 20; % Minimum shape area
params.shapeCloseRadius = 3; % Morphological closing radius (fills holes)
params.shapeFillHoles = false; % Ensures shapes are solid regions
params.intensityThreshFraction = 0.45; % Fraction of max intensity to keep
params.edgeSigma = 1.0; % Gaussian smoothing for Canny
params.edgeThresholdLow = 0.25; % Low Canny threshold fraction
params.edgeThresholdHigh = 0.55; % High Canny threshold fraction
params.intensityThreshFraction = 0.4499; % Fraction of max intensity to keep
params.edgeSigma = 1.1749; % Gaussian smoothing for Canny
params.edgeThresholdLow = 0.3383; % Low Canny threshold fraction
params.edgeThresholdHigh = 0.6412; % High Canny threshold fraction
params.pixelSize = 5.86e-6; % Physical pixel size (meters/pixel)
params.magnification = 23.94; % Magnification factor of imaging system
@ -116,5 +117,3 @@ params.magnification = 23.94; % Magnification factor of im
Analyzer.runInteractiveFeatureDetectorGUI(od_imgs, scan_parameter_values, file_list, options, params)

View File

@ -0,0 +1,77 @@
%% generateImagesForOptimization.m
% Converts a cell array of images into PNG files for Bayesian optimization
% Launches Image Segmenter for manually labeling ground truth masks
% ----------------- USER INPUTS -----------------
outputFolder = 'OptimizationImages'; % folder to save PNGs
maskFolder = 'OptimizationMasks'; % folder to save masks
startIdx = 1; % first image to export
endIdx = 10; % last image to export (set = 1 for a single test image)
createMasks = true; % set true to generate masks manually
% -----------------------------------------------
% Create folders if they don't exist
if ~exist(outputFolder,'dir')
mkdir(outputFolder);
end
if createMasks && ~exist(maskFolder,'dir')
mkdir(maskFolder);
end
% Loop over selected images
for idx = startIdx:endIdx
img = od_imgs{idx};
% Convert to double if not already
if ~isa(img,'double')
img = im2double(img);
end
% Normalize to [0,1] for PNG
img = img - min(img(:));
img = img / max(img(:));
% Save image
imgName = sprintf('image_%03d.png', idx);
imwrite(img, fullfile(outputFolder, imgName));
% Optional: create ground truth mask using Image Segmenter
if createMasks
fprintf('\n[INFO] Segment image %d/%d using Image Segmenter...\n', idx, endIdx);
% Launch Image Segmenter with the image
imageSegmenter(img);
% --- Instructions for user ---
% 1. Segment the regions you want as ground truth.
% 2. Click "Export" "Export to Workspace", export the binary mask
% 3. Close the Image Segmenter when done.
% ------------------------------------------------
% Pause script until user confirms they exported the mask
input('\n[INFO] After exporting mask from Image Segmenter, press Enter to continue...','s');
% Look for any BW variables in base workspace
bwVars = evalin('base', 'who(''BW*'')');
if ~isempty(bwVars)
% Pick the last one (e.g. BW, BW1, BW2, ...)
latestBW = bwVars{end};
gtMask = evalin('base', latestBW);
% Save mask
maskName = sprintf('mask_%03d.png', idx);
imwrite(gtMask, fullfile(maskFolder, maskName));
fprintf('\n[INFO] Mask saved: %s (from %s)\n', maskName, latestBW);
% Clear masks from workspace to avoid confusion
for v = 1:numel(bwVars)
evalin('base', sprintf('clear %s', bwVars{v}));
end
else
warning('No mask found. Skipping mask save for image %d.', idx);
end
end
end
fprintf('\n[INFO] Image and mask export complete.\n');

View File

@ -0,0 +1,95 @@
%% optimizePipeline.m
% Bayesian optimization of patch detection + shape extraction pipeline
% Joint optimization across all training images
% ------------------ USER INPUTS ------------------
imageFolder = 'OptimizationImages';
maskFolder = 'OptimizationMasks';
% Load all images and ground truth masks
imgFiles = dir(fullfile(imageFolder, '*.png'));
maskFiles = dir(fullfile(maskFolder, '*.png'));
nImages = min(numel(imgFiles), numel(maskFiles));
imgs = cell(1, nImages);
masks = cell(1, nImages);
for i = 1:nImages
imgs{i} = im2double(imread(fullfile(imageFolder, imgFiles(i).name)));
masks{i} = imread(fullfile(maskFolder, maskFiles(i).name)) > 0;
end
% ------------------ BASE PARAMETERS ------------------
params = struct();
params.backgroundDiskFraction = 1/8; % Fraction of image size used for morphological opening
params.boundingBoxPadding = 12; % Padding around detected cloud
% Initial (starting) values for optimizable parameters
params.dogGaussianSmallSigma = 1.2412;
params.dogGaussianLargeSigma = 3.9609;
params.adaptiveSensitivity = 0.4146;
params.adaptiveNeighborhoodSize = 17;
params.minPeakFraction = 0.7610;
params.minimumPatchArea = 80;
params.intensityThreshFraction = 0.4879;
params.edgeSigma = 1.0244;
params.edgeThresholdLow = 0.2384;
params.edgeThresholdHigh = 0.6127;
% Fixed shape extraction parameters
params.shapeMinArea = 20;
params.shapeCloseRadius = 3;
params.shapeFillHoles = false;
% Fixed imaging system metadata
params.pixelSize = 5.86e-6; % meters/pixel
params.magnification = 23.94;
% ------------------ OPTIMIZATION VARIABLES ------------------
optimVars = [
optimizableVariable('dogGaussianSmallSigma',[1.0,1.5])
optimizableVariable('dogGaussianLargeSigma',[3.5,4.0])
optimizableVariable('adaptiveSensitivity',[0.3,0.6])
optimizableVariable('adaptiveNeighborhoodSize',[11,19],'Type','integer')
optimizableVariable('minPeakFraction',[0.7,0.85])
optimizableVariable('minimumPatchArea',[80,150],'Type','integer')
optimizableVariable('intensityThreshFraction',[0.35,0.55])
optimizableVariable('edgeSigma',[0.8,1.2])
optimizableVariable('edgeThresholdLow',[0.2,0.35])
optimizableVariable('edgeThresholdHigh',[0.45,0.65])
];
% ------------------ OBJECTIVE FUNCTION ------------------
objFcn = @(x) mean(cellfun(@(im,mask) ...
Analyzer.evaluateFeatureDetectionPipeline(x, im, mask, params, true), ...
imgs, masks));
% ------------------ RUN BAYESIAN OPTIMIZATION ------------------
results = bayesopt(objFcn, ...
optimVars, ...
'MaxObjectiveEvaluations', 100, ...
'AcquisitionFunctionName','expected-improvement-plus', ...
'Verbose',1, ...
'PlotFcn',{@plotMinObjective,@plotAcquisitionFunction});
% ------------------ EXTRACT BEST PARAMETERS ------------------
bestX = results.XAtMinObjective;
params.dogGaussianSmallSigma = bestX.dogGaussianSmallSigma;
params.dogGaussianLargeSigma = bestX.dogGaussianLargeSigma;
params.adaptiveSensitivity = bestX.adaptiveSensitivity;
params.adaptiveNeighborhoodSize = bestX.adaptiveNeighborhoodSize;
params.minPeakFraction = bestX.minPeakFraction;
params.minimumPatchArea = bestX.minimumPatchArea;
params.intensityThreshFraction = bestX.intensityThreshFraction;
params.edgeSigma = bestX.edgeSigma;
params.edgeThresholdLow = bestX.edgeThresholdLow;
params.edgeThresholdHigh = bestX.edgeThresholdHigh;
disp('Best parameters found (joint optimization across all images):');
disp(params);
%% --- Evaluate with best parameters, plot results on first image ---
fprintf('Evaluating with best parameters and plotting results...\n');
Analyzer.runInteractiveFeatureDetectorGUI(od_imgs, scan_parameter_values, file_list, options, params)