%UNMIXCHANNELS   Unmixes the signals in the two input channels
%   [OUTPUT,H,VAL,POINTS] = UNMIX2CHANNELS(CHAN1,CHAN2,MASK,BG1,BG2)
%
%   CHAN1 and CHAN2 are the two input channels (any dimensionality), and
%   are expected to be 12-bit images. MASK is a binary mask image
%   specifying which pixels to use to compute the cross talk.
%
%   BG1 and BG2 are the background levels for the two channels and default
%   to 0.
%
%   OUTPUT is CHAN2 with the cross talk removed.
%
%   H is the joint histogram used to determine the cross talk. POINTS are
%   the coordinates of the relevant points found in the histogram, and VAL
%   is the slope and offset of the line fitted to the points.
%
%   VAL(1) is the cross talk level for CHAN1 into CHAN2.
%   OUTPUT = CHAN2 - (CHAN1-BG1)*VAL(1).

% (c)2005-2007, Cris Luengo, Lawrence Berkeley National Laboratory.
% This code runs on MATLAB 7.0 or later with the DIPimage toolbox version
% 1.5 or later.

function [chan2,h,val,points] = unmix2channels(chan1,chan2,mask,bg1,bg2)

val = [0;0];
points = [];

if nargin<2
   error('More input arguments expected!')
end
if nargin<5
   bg2 = 0;
end
if nargin<4
   bg1 = 0;
end

% Select the pixels
if nargin>2
   c1 = chan1(mask);
   c2 = chan2(mask);
else
   c1 = chan1;
   c2 = chan2;
end
mx = 4064;
I = ~(c1>=mx | c2>=mx);
c1 = c1(I)-bg1;
c2 = c2(I)-bg2;

% Create the joint histogram
h = double(mdhistogram(newimar(c1,c2),...
                       {{'lower_abs','lower',-500,'upper_abs','upper',4000,'bins',100,'no_correction'},...
                        {'lower_abs','lower',-500,'upper_abs','upper',4000,'bins',100,'no_correction'}}));
h2 = double(gaussf(h,1));

% Detect points in the histogram
c2 = zeros(size(h,2),1);
for ii=size(h,2):-1:ceil(700/45) % Magic number! The index for ~200 above 0.
   c = h2(:,ii);
   if sum(c)>100 & sum(c(1:ii))>sum(c(ii:end)) 
      [tmp,a] = max(c);
      if sum(c(1:a))<=sum(c(a:end))
         c2(ii) = a;
      end
   end
end
c1 = find(c2>0);
if ~isempty(c1)
   c2 = (c2(c1)-1)*45-500;
   c1 = (c1-1)*45-500;
   g = corrcoef(c2,c1);
   if length(c2)>=8 & g(2) > 0.7
      a = [c1,ones(length(c1),1)]\c2;
      b = a(2);
      a = a(1);
      if a > 0.05 & a < 0.8       % Sanity check. Also, we don't care if the bleedthrough is less than 5% or so.
         val = [a;b];
         chan2 = chan2 - (chan1-bg1)*a;
         points = [c2,c1];
      end
   end
end