From 3e6300203117dfe877f6f2e766e7dc9318a44de2 Mon Sep 17 00:00:00 2001 From: kalidke Date: Tue, 4 Nov 2025 16:47:35 -0700 Subject: [PATCH] Fix sCMOS scalar calibration and non-square image bug (Issue #37) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- MATLAB/+smi_core/@FindROI/FindROI.m | 68 +++++++++++++------ .../SingleMoleculeFitting.m | 9 +-- 2 files changed, 51 insertions(+), 26 deletions(-) diff --git a/MATLAB/+smi_core/@FindROI/FindROI.m b/MATLAB/+smi_core/@FindROI/FindROI.m index 398f8dfb..1f060ad0 100644 --- a/MATLAB/+smi_core/@FindROI/FindROI.m +++ b/MATLAB/+smi_core/@FindROI/FindROI.m @@ -83,11 +83,29 @@ switch SMF.Data.CameraType case 'SCMOS' obj.IsSCMOS = 1; - obj.Varim = gpuArray(single(SMF.Data.CameraReadNoise./SMF.Data.CameraGain.^2)); + % For sCMOS with scalar calibration (modern uniform + % sensors like Orca Fusion), expand to 2D variance image + % matching data dimensions. This supports Issue #37 fix. + CameraReadNoise = SMF.Data.CameraReadNoise; + CameraGain = SMF.Data.CameraGain; + if isscalar(CameraReadNoise) && isscalar(CameraGain) + % Need 2D variance image for sCMOS CUDA kernels + % Expand based on Data dimensions if available + if nargin > 1 && ~isempty(Data) + [NRows, NCols, ~] = size(Data); + CameraReadNoise = CameraReadNoise * ones(NRows, NCols, 'single'); + CameraGain = CameraGain * ones(NRows, NCols, 'single'); + else + % Data not yet available - will fail if used + % This case shouldn't happen in normal workflow + warning('FindROI: sCMOS with scalar calibration but no Data provided'); + end + end + obj.Varim = gpuArray(single(CameraReadNoise./CameraGain.^2)); otherwise obj.IsSCMOS = 0; end - + end if nargin > 1 @@ -306,24 +324,24 @@ function showOverlay(obj) end function [Out] = gauss_sCMOS(Stack, Varim, Sigma) - %gauss_sCMOS() Out of place, weighted Gaussian filter - % [SubStack] = gauss_sCMOS(SubStack, Varim, Sigma) + %gauss_sCMOS() Variance-weighted Gaussian filter for sCMOS + % [Out] = gauss_sCMOS(Stack, Varim, Sigma) % - % Filtering is applied using direct, weighted calcualtion, but done - % separably along x and y. + % Filtering is applied using variance-weighted convolution, done + % separably along x and y. This implements the optimal filtering + % for sCMOS cameras with pixel-dependent noise (Huang et al. 2013). % - % This is an out of place filter, meaning that inputs are not - % modified. + % This is an out-of-place filter (input Stack is not modified). % % If the input is not gpuArray, it is copied to the GPU - % + % % INPUTS: % Stack: Stack of 2D images to be filtered - % Varim: A gain-corrected variance image + % Varim: A gain-corrected variance image (2D array) % Sigma: Gaussian Kernel (Pixels) % % OUTPUTS: - % Out: 2D stack of images (gpuArray) + % Out: Filtered 2D stack of images (gpuArray) % % REQUIRES: % MATLAB 2014a or later versions @@ -332,21 +350,27 @@ function showOverlay(obj) % smi_cuda_FindROI.ptx % smi_cuda_FindROI.cu % - - Out1=gpuArray(zeros(size(Stack),'single')); - + + % Create working copy (matches gaussInPlace pattern) + Out = gpuArray(zeros(size(Stack),'single')); + %Creating GPU CUDA kernel objects from PTX and CU code K_Gauss = parallel.gpu.CUDAKernel('smi_cuda_FindROI.ptx','smi_cuda_FindROI.cu','kernel_gaussMajor_sCMOS'); K_Gauss.GridSize(1) = size(Stack,3); K_Gauss.ThreadBlockSize(1) = size(Stack,2); - - %Calling the gpu code to apply Gaussian filter along major - Out1 = feval(K_Gauss,gpuArray(Stack),Varim,Out1,size(Stack,1),Sigma); - - %Permuting and doing the other dimension - Out = gpuArray(zeros(size(Stack),'single')); - Out = feval(K_Gauss,permute(Out1,[2 1 3]),permute(Varim,[2 1]),Out,size(Stack,1),Sigma); - + + %Calling the gpu code to apply Gaussian filter along major dimension + Out = feval(K_Gauss,gpuArray(Stack),Varim,Out,size(Stack,1),Sigma); + + %Permute for second pass (matches gaussInPlace pattern) + Out = permute(Out,[2 1 3]); + Varim = permute(Varim,[2 1]); + + %Second pass along other dimension - use current Out dimensions + K_Gauss.GridSize(1) = size(Out,3); + K_Gauss.ThreadBlockSize(1) = size(Out,2); + Out = feval(K_Gauss,Out,Varim,Out,size(Out,1),Sigma); + %Permute back Out = permute(Out,[2 1 3]); diff --git a/MATLAB/+smi_core/@SingleMoleculeFitting/SingleMoleculeFitting.m b/MATLAB/+smi_core/@SingleMoleculeFitting/SingleMoleculeFitting.m index 68d21358..8ef5169c 100644 --- a/MATLAB/+smi_core/@SingleMoleculeFitting/SingleMoleculeFitting.m +++ b/MATLAB/+smi_core/@SingleMoleculeFitting/SingleMoleculeFitting.m @@ -391,10 +391,11 @@ 'precedence over the inclusion set']); obj.SMFFieldNotes.Data.CameraType.Tip = ... sprintf(['Type of camera used to collect the raw data.\n', ... - 'CameraGain, CameraOffset, CameraReadNoise) should be\n', ... - 'scalars if the CameraType is EMCCD while if SCMOS,\n', ... - 'square arrays taken from the file located at the ', ... - 'CalibrationFilePath.']); + 'CameraGain, CameraOffset, CameraReadNoise should be:\n', ... + ' EMCCD: scalars\n', ... + ' SCMOS: either square arrays from CalibrationFilePath\n', ... + ' OR scalars for uniform sensors (e.g., Orca Fusion).\n', ... + 'Scalar values for SCMOS will be expanded automatically.']); obj.SMFFieldNotes.Data.CameraGain.Tip = ... 'Gain of the camera used to collect the raw data'; obj.SMFFieldNotes.Data.CameraOffset.Tip = ...