From 546457fe2990d74625cd8be5ee66528a5ff10a84 Mon Sep 17 00:00:00 2001
From: ac_fx <arc334@drexel.edu>
Date: Sun, 12 Nov 2023 13:03:50 -0500
Subject: [PATCH] SSF clip quantiles

---
 matlab/+SSF/getClipLimits.m    |  9 ++++-----
 matlab/+SSF/getClipQuantiles.m | 12 ++++++++++++
 matlab/+SSF/go_ssf_ncd.m       |  2 +-
 3 files changed, 17 insertions(+), 6 deletions(-)
 create mode 100644 matlab/+SSF/getClipQuantiles.m

diff --git a/matlab/+SSF/getClipLimits.m b/matlab/+SSF/getClipLimits.m
index c3d93db2..5eb239d7 100644
--- a/matlab/+SSF/getClipLimits.m
+++ b/matlab/+SSF/getClipLimits.m
@@ -29,10 +29,9 @@ parfor ff=1:length(flist)
 end 
 
 kymoPixels = vertcat(kymoPixels{:});
-cl = [mean(kymoPixels)-std(kymoPixels), mean(kymoPixels)+std(kymoPixels)];
-q1 = length(find(kymoPixels<cl(1)))./length(kymoPixels);
-q2 = 1 - length(find(kymoPixels>cl(2)))./length(kymoPixels);
-pxTarget = linspace(q1,q2,254);
-clipLimits = quantile(kymoPixels,pxTarget);
+
+clip_neg = SSF.getClipQuantiles(kymoPixels(kymoPixels<0),127);
+clip_pos = SSF.getClipQuantiles(kymoPixels(kymoPixels>0),127);
+clipLimits = [clip_neg,clip_pos];
 clipLimits = repmat({clipLimits},length(flist),1);
 4;
diff --git a/matlab/+SSF/getClipQuantiles.m b/matlab/+SSF/getClipQuantiles.m
new file mode 100644
index 00000000..b6cf656d
--- /dev/null
+++ b/matlab/+SSF/getClipQuantiles.m
@@ -0,0 +1,12 @@
+function clipLimits = getClipQuantiles(voxels,nQuant)
+% 
+% mx = mean(voxels);
+% sx = std(voxels);
+% cl = [mx-sx, mx+sx];
+% q1 = length(find(voxels<cl(1)))./length(voxels);
+% q2 = 1 - length(find(voxels>cl(2)))./length(voxels);
+% pxTarget = linspace(q1,q2,nQuant);
+% clipLimits = quantile(voxels,pxTarget);
+
+pxTarget = linspace(2.5,97.5,nQuant);
+clipLimits = prctile(voxels,pxTarget);
diff --git a/matlab/+SSF/go_ssf_ncd.m b/matlab/+SSF/go_ssf_ncd.m
index db4ad479..5a06be40 100644
--- a/matlab/+SSF/go_ssf_ncd.m
+++ b/matlab/+SSF/go_ssf_ncd.m
@@ -24,7 +24,7 @@ kymoPixels = {};
  
 p = ljsStartParallel();
 % 
-if ~exist('clipLimits','var')
+if ~exist('clipLimits','var') || isempty(clipLimits)
     clipLimits = SSF.getClipLimits(flist,targetChannels,nKeep,tClip);
 end
 %
-- 
GitLab