Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 22 additions & 0 deletions MATLAB/+smi_core/@DataToPhotons/convertToPhotons.m
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,28 @@
CameraOffset = single(CameraOffset);
CameraReadNoise = single(CameraReadNoise);

% For sCMOS cameras with scalar calibration values (modern low-variation
% sensors like Orca Fusion), expand scalars to 2D arrays matching image
% dimensions. This ensures consistent behavior throughout the pipeline while
% supporting simplified calibration for uniform sensors.
if isscalar(CameraGain) && isscalar(CameraOffset) && ...
(isempty(CameraReadNoise) || isscalar(CameraReadNoise))
% Get image dimensions from RawData
[NRows, NCols, ~] = size(RawData);

% Expand scalars to 2D arrays
CameraGain = CameraGain * ones(NRows, NCols, 'single');
CameraOffset = CameraOffset * ones(NRows, NCols, 'single');
if ~isempty(CameraReadNoise)
CameraReadNoise = CameraReadNoise * ones(NRows, NCols, 'single');
end

% Update CalibrationROI to match full image size if not properly defined
if isempty(CalibrationROI) || isequal(CalibrationROI, [1, 1, 1, 1])
CalibrationROI = [1, 1, NRows, NCols];
end
end

% Compare array sizes and ROI definitions to ensure consistency.
SizeRawData = size(RawData);
SizeCameraGain = size(CameraGain);
Expand Down
36 changes: 34 additions & 2 deletions MATLAB/+smi_core/@DataToPhotons/unitTest.m
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@
% units of photons^2 for non-scalar gain and offset.
% Success(4): Same as Success(1) for scalar gain.
% Success(5): Same as Success(2) for scalar offset.
% Success(6): Test sCMOS camera with scalar calibration (Orca
% Fusion scenario) - scalars auto-expand to 2D arrays.

% Created by:
% David J. Schodt (Lidke Lab, 2020)
Expand Down Expand Up @@ -89,7 +91,7 @@
SMF.Data.CameraGain = CameraGain;
SMF.Data.CameraOffset = CameraOffset;
SMF.Data.CameraReadNoise = CameraReadNoise;
Success = zeros(1, 5, 'logical');
Success = zeros(1, 6, 'logical');
try
DTP = smi_core.DataToPhotons(SMF, RawDataBottomRight, ...
ROIBottomRight, CalibrationROI, true);
Expand Down Expand Up @@ -161,7 +163,37 @@
DTP.CameraReadNoise = CameraReadNoise;
[CorrectedData, CorrectedNoise] = DTP.convertData();
Success(4) = all(abs(CorrectedData(:)-DataWithScalarReadNoise(:)) < 0.1);
Success(5) = (abs(CorrectedNoise-ScalarReadNoisePhotons) < 0.1);
% CorrectedNoise may be expanded to 2D array by scalar expansion logic
if isscalar(CorrectedNoise)
Success(5) = (abs(CorrectedNoise-ScalarReadNoisePhotons) < 0.1);
else
% All elements should equal the scalar value
Success(5) = all(abs(CorrectedNoise(:)-ScalarReadNoisePhotons) < 0.1);
end

% Test sCMOS camera with scalar calibration (modern uniform sensors like
% Orca Fusion). This tests the auto-expansion of scalar values to 2D arrays.
SMF_sCMOS = smi_core.SingleMoleculeFitting;
SMF_sCMOS.Data.CameraType = 'SCMOS';
SMF_sCMOS.Data.CameraGain = CameraGain; % Use the scalar value already defined
SMF_sCMOS.Data.CameraOffset = CameraOffset;
SMF_sCMOS.Data.CameraReadNoise = CameraReadNoise;
SMF_sCMOS.Data.CalibrationFilePath = ''; % Empty to trigger scalar path
try
DTP_sCMOS = smi_core.DataToPhotons(SMF_sCMOS, RawData, [], [], 0, false);
[CorrectedData_sCMOS, CorrectedNoise_sCMOS] = DTP_sCMOS.convertData();
% Verify the conversion worked correctly
% Note: The scalar expansion happens internally in convertToPhotons,
% so DTP_sCMOS.CameraGain remains scalar, but the output arrays are
% properly sized and the conversion is correct.
Success(6) = all(abs(CorrectedData_sCMOS(:)-DataWithScalarReadNoise(:)) < 0.1) ...
&& (size(CorrectedNoise_sCMOS, 1) == FrameSizeFull) ...
&& (size(CorrectedNoise_sCMOS, 2) == FrameSizeFull);
catch MException
warning('sCMOS scalar test failed: %s - %s', ...
MException.identifier, MException.message)
Success(6) = false;
end


end
5 changes: 5 additions & 0 deletions MATLAB/+smi_core/@LocalizeData/genLocalizations.m
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,11 @@
SMDCandidates.X = SMDCandidates.X + SMDCandidates.XBoxCorner;
SMDCandidates.Y = SMDCandidates.Y + SMDCandidates.YBoxCorner;

% Merge duplicate localizations that may have been generated when
% FindROI detects multiple identical pixel maxima (Issue #21).
% This is rare but can occur in noise-free simulation data.
SMDCandidates = obj.mergeDuplicateLocalizations(SMDCandidates);

% Threshold localizations.
MinMax = smi_core.Threshold.setMinMax(obj.SMF);
if obj.SMF.Thresholding.AutoThreshLogL
Expand Down
89 changes: 89 additions & 0 deletions MATLAB/+smi_core/@LocalizeData/mergeDuplicateLocalizations.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,89 @@
function [SMD] = mergeDuplicateLocalizations(obj, SMD)
%mergeDuplicateLocalizations merges localizations at same position/frame.
% This method identifies and merges duplicate localizations that occur at
% the same position within the same frame. This can happen when FindROI
% detects multiple identical pixel maxima (common in noise-free simulation
% data). When duplicates are found, the localization with the higher
% LogLikelihood is kept.
%
% INPUTS:
% SMD: Single Molecule Data structure with candidate localizations
%
% OUTPUTS:
% SMD: Modified SMD structure with duplicates removed
%
% NOTE: This addresses Issue #21 where FindROI can generate multiple boxes
% for the same emitter when pixels have identical maximum values.

% Created by:
% Claude/Keith Lidke (Lidke Lab, 2025)


% Return early if no localizations to process
if isempty(SMD.X) || numel(SMD.X) < 2
return;
end

% Distance threshold for considering duplicates (pixels)
% Two localizations are considered duplicates if they are within this
% distance in the same frame. 0.1 pixels is conservative - true duplicates
% from identical pixel maxima will be < 0.01 pixels apart.
dist_threshold = 0.1;

% Find duplicates frame-by-frame
unique_frames = unique(SMD.FrameNum);
keep_indices = true(size(SMD.X));
n_duplicates_found = 0;

for ff = 1:length(unique_frames)
frame_mask = SMD.FrameNum == unique_frames(ff);
frame_indices = find(frame_mask);

if length(frame_indices) < 2
continue; % No duplicates possible with single localization
end

% Check pairwise distances within this frame
for ii = 1:length(frame_indices)-1
if ~keep_indices(frame_indices(ii))
continue; % Already marked for removal
end

for jj = ii+1:length(frame_indices)
if ~keep_indices(frame_indices(jj))
continue; % Already marked for removal
end

% Calculate Euclidean distance
dx = SMD.X(frame_indices(ii)) - SMD.X(frame_indices(jj));
dy = SMD.Y(frame_indices(ii)) - SMD.Y(frame_indices(jj));
dist = sqrt(dx^2 + dy^2);

if dist < dist_threshold
% Duplicates found - keep the one with higher LogLikelihood
% (better fit quality)
if SMD.LogLikelihood(frame_indices(ii)) >= ...
SMD.LogLikelihood(frame_indices(jj))
keep_indices(frame_indices(jj)) = false;
n_duplicates_found = n_duplicates_found + 1;
else
keep_indices(frame_indices(ii)) = false;
n_duplicates_found = n_duplicates_found + 1;
break; % Move to next ii since this one is removed
end
end
end
end
end

% Extract only the non-duplicate localizations
if n_duplicates_found > 0
if obj.Verbose > 1
fprintf(['\tLocalizeData.mergeDuplicateLocalizations(): ', ...
'Removed %d duplicate localizations\n'], n_duplicates_found);
end
SMD = smi_core.SingleMoleculeData.isolateSubSMD(SMD, keep_indices);
end


end
Loading