Added more robust rolling-ball baseline correction. Added ability tweak

baseline to avoid negative values.

Author: dpfoose

fix_spectra/fix_click_menu.m

@@ -33,14 +33,14 @@ function fix_click_menu(hObject, eventdata, handles) %#ok<INUSD>
 %
 % Dan Homer ????
 %
-% Paul Anderson ????
+% % Paul Anderson ????
 %
 % Eric Moyer 2011-2012 (eric_moyer@yahoo.com)
 
 str = {'Get collection(s)','Load collections', ...
 ...%'','Set noise regions', ...
 	'','Load regions','Create signal map','Create new region', ...
-       'Edit region','Delete region','Clear regions','Fix baseline', 'Fix baseline (rolling ball)', ...
+       'Edit region','Delete region','Clear regions','Fix baseline', 'Rolling ball baseline', 'Fix negative points to 0', ...
     '','Crop','Set reference','Normalize to reference','Zero regions', 'Sum normalize','Normalize to weight',...     
     '','Save regions','Finalize','Save collections', ...
        'Merge processing logs','Merge locally' ...
@@ -116,8 +116,10 @@ function fix_click_menu(hObject, eventdata, handles) %#ok<INUSD>
     set_noise_regions;
 elseif strcmp(str{s},'Fix baseline')
     fix_baseline
-elseif strcmp(str{s},'Fix baseline (rolling ball)')
+elseif strcmp(str{s},'Rolling ball baseline')
     fix_rolling_ball_baseline
+elseif strcmp(str{s},'Fix negative points to 0')
+    fix_negative_points
 elseif strcmp(str{s},'Edit region')
     set_edit
 elseif strcmp(str{s},'Create new region')

fix_spectra/fix_negative_points.m

@@ -0,0 +1,14 @@
+function fix_negative_points
+%fix_negative_points adjust baseline so that negative points are set to 0.
+collections = getappdata(gcf,'collections');
+for c = 1:length(collections)
+    for s = 1:collections{c}.num_samples
+        inds = collections{c}.Y_fixed(:,s) < 0;
+        collections{c}.Y_fixed(inds,s) = 0;
+        collections{c}.Y_baseline(inds,s) = collections{c}.Y(inds,s) - collections{c}.Y_fixed(inds,s);
+    end
+end
+setappdata(gcf,'collections',collections);
+plot_all
+end
+

fix_spectra/fix_rolling_ball_baseline.m

@@ -0,0 +1,30 @@
+function fix_rolling_ball_baseline
+%rolling_ball_baseline Summary of this function goes here
+%   Detailed explanation goes here
+prompt={'Enter min/max window:', 'Enter smoothing window:'};
+name='Rolling ball';
+answer = inputdlg(prompt, name, [1, 35], {'55', '35'});
+
+if(isempty(answer))
+    return
+end
+
+wm = str2double(answer{1});
+ws = str2double(answer{2});
+
+collections = getappdata(gcf,'collections');
+for c = 1:length(collections)
+    for s = 1:collections{c}.num_samples
+        y = collections{c}.Y(:,s);
+        baseline = rolling_ball_baseline(y, wm, ws);
+        collections{c}.Y_baseline(:,s) = baseline;
+        collections{c}.Y_fixed(:,s) = y - baseline;
+    end
+end
+%setappdata(gcf,'add_processing_log','Fixed baseline.');
+setappdata(gcf,'add_processing_log',sprintf('Fix baseline (rolling ball, wm: %f, ws: %f).',wm,ws));
+setappdata(gcf,'temp_suffix','_fixed_baseline');
+setappdata(gcf,'collections',collections);
+
+plot_all
+end

fix_spectra/rolling_ball_baseline.m

@@ -0,0 +1,101 @@
+function baseline = rolling_ball_baseline(y, wm, ws)
+%rolling_ball_baseline Summary of this function goes here
+% y is a column vector representing a spectrum
+% wm is the width of min/max window
+% ws is the width of the smoothing window
+% Based on the rolling ball baseline implemented in the R package baseline
+    [m, ~] = size(y);
+    T1 = zeros(size(y)); % minimizers
+    T2 = zeros(size(y)); % maximizers
+    baseline = zeros(size(y));
+
+    % minimize
+    u1 = ceil((wm+1) / 2) + 1;
+    T1(1) = min(y(1:u1));
+    % start of spectrum
+    for i=2:wm
+        u2 = u1 + 1 + mod(i, 2);
+        T1(i) = min(min(y(u1+1:u2)), T1(i-1));
+        u1 = u2;
+    end
+    % middle of spectrum
+    for i = wm+1:m-wm
+        if ((y(u1+1) <= T1(i-1)) && (y(u1-wm) ~= T1(i-1)))
+            T1(i) = y(u1+1); % next is smaller
+        else
+            T1(i) = min(y(i-wm:i+wm));
+        end
+        u1 = u1 + 1;
+    end
+    % end of spectrum
+    u1 = m - 2*wm - 1;
+    for i = m-wm+1:m
+        u2 = u1 + 1 + mod(i+1, 2);
+        if min(y(u1:u2-1)) > T1(i-1)
+            T1(i) = T1(i-1); % removed is larger
+        else
+            T1(i) = min(y(u2:m));
+        end
+        u1 = u2;
+    end
+
+    % maximize
+    u1 = ceil((wm+1)/2) + 1;
+    % start of spectrum
+    T2(1) = max(T1(1:u1));
+    for i = 2:wm
+        u2 = u1 + 1 + mod(i, 2);
+        T2(i) = max(max(T1(u1+1:u2)), T2(i-1));
+        u1 = u2;
+    end
+    % middle of spectrum
+    for i = wm+1:m-wm
+        if (T1(u1+1) >= T2(i-1)) && (T1(u1-wm) ~= T2(i-1))
+            T2(i) = T1(u1+1);
+        else
+            T2(i) = max(T1(i-wm:i+wm));
+        end
+        u1 = u1 + 1;
+    end
+    % end of spectrum
+    u1 = m - 2*wm - 1;
+    for i = m-wm+1:m
+        u2 = u1 + 1 + mod(i+1, 2);
+        if max(T1(u1:u2-1)) < T2(i-1)
+            T2(i) = T2(i-1);
+        else
+            T2(i) = max(T1(u2:m));
+        end
+        u1 = u2;
+    end
+
+    % Smoothing
+    u1 = ceil(ws/2);
+    % start of spectrum
+    v = sum(T2(1:u1));
+    for i = 1:ws
+        u2 = u1 + 1 + mod(i, 2);
+        v = v + sum(T2(u1+1:u2));
+        baseline(i) = v/u2;
+        u1 = u2;
+    end
+    % middle of spectrum
+    v = sum(T2(1:2*ws+1));
+    baseline(ws+1) = v/(2*ws+1);
+    for i = ws+2:m-ws
+        v = v - T2(i-ws-1) + T2(i+ws);
+        baseline(i) = v/(2*ws+1);
+    end
+    u1 = m - 2*ws+1;
+    v = v - T2(u1); % sum so far
+    baseline(m-ws+1) = v/(2*ws); % mean so far
+    % end of spectrum
+    for i = m-ws+2:m
+        u2 = u1+ 1 + mod(i+1,2);
+        v = v - sum(T2(u1:u2-1));
+        baseline(i) = v/(m-u2+1);
+        u1 = u2;
+    end
+
+end
+