Skip to content

Commit 3e63002

Browse files
kalidkeclaude
andcommitted
Fix sCMOS scalar calibration and non-square image bug (Issue #37)
Allow sCMOS cameras to use scalar calibration values for modern uniform sensors like Orca Fusion AND fix pre-existing bug preventing non-square sCMOS images from being processed. Issue #37 - Scalar Calibration: When CameraType='SCMOS' with scalar gain/offset/readnoise and empty CalibrationFilePath, sCMOS CUDA kernels expect 2D variance images to index into, not scalars. Pre-existing Bug - Non-Square sCMOS Images: The gauss_sCMOS() function had incorrect dimension handling after permutation, causing failures with non-square sCMOS data (e.g., 429×244 gattaquant beads). Bug was hidden with square images (256×256) where transpose doesn't change dimensions. EMCCD worked fine (different code path). Root cause: gauss_sCMOS used out-of-place pattern but referenced original Stack dimensions after permuting working arrays, unlike gaussInPlace which updates Stack reference and uses current dimensions. Changes: - FindROI.m constructor (lines 86-104): For sCMOS with scalar calibration, expand scalars to 2D variance image matching Data dimensions. - FindROI.gauss_sCMOS (lines 326-377): Refactored to match gaussInPlace pattern - use working copy and reference current dimensions after permutation, not original Stack dimensions. This fixes non-square bug and improves code consistency. - SingleMoleculeFitting.m: Update documentation to clarify sCMOS supports both pixel-wise arrays and scalars (auto-expanded). Why different filtering algorithms? - EMCCD: Uniform noise → fast recursive Gaussian (Young & van Vliet 1995) - sCMOS: Pixel-varying noise → variance-weighted convolution (Huang 2013) This algorithmic difference is necessary. Dimension handling should be consistent - now it is. Testing verified: - Non-square sCMOS + 2D arrays: ✓ (was broken, now fixed) - Non-square sCMOS + scalars: ✓ (Issue #37) - Square sCMOS: ✓ (backward compatible) - DataToPhotons: ✓ (scalars handled correctly) Fixes #37 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude <noreply@anthropic.com>
1 parent 7bef54a commit 3e63002

2 files changed

Lines changed: 51 additions & 26 deletions

File tree

MATLAB/+smi_core/@FindROI/FindROI.m

Lines changed: 46 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -83,11 +83,29 @@
8383
switch SMF.Data.CameraType
8484
case 'SCMOS'
8585
obj.IsSCMOS = 1;
86-
obj.Varim = gpuArray(single(SMF.Data.CameraReadNoise./SMF.Data.CameraGain.^2));
86+
% For sCMOS with scalar calibration (modern uniform
87+
% sensors like Orca Fusion), expand to 2D variance image
88+
% matching data dimensions. This supports Issue #37 fix.
89+
CameraReadNoise = SMF.Data.CameraReadNoise;
90+
CameraGain = SMF.Data.CameraGain;
91+
if isscalar(CameraReadNoise) && isscalar(CameraGain)
92+
% Need 2D variance image for sCMOS CUDA kernels
93+
% Expand based on Data dimensions if available
94+
if nargin > 1 && ~isempty(Data)
95+
[NRows, NCols, ~] = size(Data);
96+
CameraReadNoise = CameraReadNoise * ones(NRows, NCols, 'single');
97+
CameraGain = CameraGain * ones(NRows, NCols, 'single');
98+
else
99+
% Data not yet available - will fail if used
100+
% This case shouldn't happen in normal workflow
101+
warning('FindROI: sCMOS with scalar calibration but no Data provided');
102+
end
103+
end
104+
obj.Varim = gpuArray(single(CameraReadNoise./CameraGain.^2));
87105
otherwise
88106
obj.IsSCMOS = 0;
89107
end
90-
108+
91109
end
92110

93111
if nargin > 1
@@ -306,24 +324,24 @@ function showOverlay(obj)
306324
end
307325

308326
function [Out] = gauss_sCMOS(Stack, Varim, Sigma)
309-
%gauss_sCMOS() Out of place, weighted Gaussian filter
310-
% [SubStack] = gauss_sCMOS(SubStack, Varim, Sigma)
327+
%gauss_sCMOS() Variance-weighted Gaussian filter for sCMOS
328+
% [Out] = gauss_sCMOS(Stack, Varim, Sigma)
311329
%
312-
% Filtering is applied using direct, weighted calcualtion, but done
313-
% separably along x and y.
330+
% Filtering is applied using variance-weighted convolution, done
331+
% separably along x and y. This implements the optimal filtering
332+
% for sCMOS cameras with pixel-dependent noise (Huang et al. 2013).
314333
%
315-
% This is an out of place filter, meaning that inputs are not
316-
% modified.
334+
% This is an out-of-place filter (input Stack is not modified).
317335
%
318336
% If the input is not gpuArray, it is copied to the GPU
319-
%
337+
%
320338
% INPUTS:
321339
% Stack: Stack of 2D images to be filtered
322-
% Varim: A gain-corrected variance image
340+
% Varim: A gain-corrected variance image (2D array)
323341
% Sigma: Gaussian Kernel (Pixels)
324342
%
325343
% OUTPUTS:
326-
% Out: 2D stack of images (gpuArray)
344+
% Out: Filtered 2D stack of images (gpuArray)
327345
%
328346
% REQUIRES:
329347
% MATLAB 2014a or later versions
@@ -332,21 +350,27 @@ function showOverlay(obj)
332350
% smi_cuda_FindROI.ptx
333351
% smi_cuda_FindROI.cu
334352
%
335-
336-
Out1=gpuArray(zeros(size(Stack),'single'));
337-
353+
354+
% Create working copy (matches gaussInPlace pattern)
355+
Out = gpuArray(zeros(size(Stack),'single'));
356+
338357
%Creating GPU CUDA kernel objects from PTX and CU code
339358
K_Gauss = parallel.gpu.CUDAKernel('smi_cuda_FindROI.ptx','smi_cuda_FindROI.cu','kernel_gaussMajor_sCMOS');
340359
K_Gauss.GridSize(1) = size(Stack,3);
341360
K_Gauss.ThreadBlockSize(1) = size(Stack,2);
342-
343-
%Calling the gpu code to apply Gaussian filter along major
344-
Out1 = feval(K_Gauss,gpuArray(Stack),Varim,Out1,size(Stack,1),Sigma);
345-
346-
%Permuting and doing the other dimension
347-
Out = gpuArray(zeros(size(Stack),'single'));
348-
Out = feval(K_Gauss,permute(Out1,[2 1 3]),permute(Varim,[2 1]),Out,size(Stack,1),Sigma);
349-
361+
362+
%Calling the gpu code to apply Gaussian filter along major dimension
363+
Out = feval(K_Gauss,gpuArray(Stack),Varim,Out,size(Stack,1),Sigma);
364+
365+
%Permute for second pass (matches gaussInPlace pattern)
366+
Out = permute(Out,[2 1 3]);
367+
Varim = permute(Varim,[2 1]);
368+
369+
%Second pass along other dimension - use current Out dimensions
370+
K_Gauss.GridSize(1) = size(Out,3);
371+
K_Gauss.ThreadBlockSize(1) = size(Out,2);
372+
Out = feval(K_Gauss,Out,Varim,Out,size(Out,1),Sigma);
373+
350374
%Permute back
351375
Out = permute(Out,[2 1 3]);
352376

MATLAB/+smi_core/@SingleMoleculeFitting/SingleMoleculeFitting.m

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -391,10 +391,11 @@
391391
'precedence over the inclusion set']);
392392
obj.SMFFieldNotes.Data.CameraType.Tip = ...
393393
sprintf(['Type of camera used to collect the raw data.\n', ...
394-
'CameraGain, CameraOffset, CameraReadNoise) should be\n', ...
395-
'scalars if the CameraType is EMCCD while if SCMOS,\n', ...
396-
'square arrays taken from the file located at the ', ...
397-
'CalibrationFilePath.']);
394+
'CameraGain, CameraOffset, CameraReadNoise should be:\n', ...
395+
' EMCCD: scalars\n', ...
396+
' SCMOS: either square arrays from CalibrationFilePath\n', ...
397+
' OR scalars for uniform sensors (e.g., Orca Fusion).\n', ...
398+
'Scalar values for SCMOS will be expanded automatically.']);
398399
obj.SMFFieldNotes.Data.CameraGain.Tip = ...
399400
'Gain of the camera used to collect the raw data';
400401
obj.SMFFieldNotes.Data.CameraOffset.Tip = ...

0 commit comments

Comments
 (0)