Calculate_UnMixing_Fiji.m

% @matrix Image
% @matrix dSpectrum
% @matrix aSpectrum
% @matrix tSpectrum
% @double dQY
% @double aQY
% @double BackGround
% @double Threshold
% @OUTPUT double[] KDA
% @OUTPUT double[] KAD
% @OUTPUT double[] K3
% @OUTPUT double[] FDonly
% @OUTPUT double[] Eapp

%function [KDA, KAD, K3, FDonly, Eapp] = Calculate_UnMixing_Data(Image, dSpectrum, aSpectrum, tSpectrum, BackGround, Threshold)
%------------------------------------------------------------------------------
% this function performs the unmixing procedure for three components
% unmixing. 
%
% The input parameters are:
% Image      - The mixed image that contains 3-D matrix two of the dimensions
%              are x and y direcion in the image and the third dimension is the
%              wavelength dimension.
% dSpectrum  - The donors elementary spectrum
% dQY        - Donor quantum yeild
% aSpectrum  - The acceptor elementary spectrum.
% aQY        - Acceptor quantum yeild
% tSpectrum  - Third component spectrum
% BackGround - The bais value substructed from the image
% Threshold  - The value above which the signal is considered and below
%              which the signal is equal to 0;
%
% The output parameters are:
% KDA        - Emission intensity of the donor in the presence of an
%              acceptor
% KAD        - Emission intensity of the Accepto in the presence of a
%              donor
% K3         - Emission intensity of the Third component.
% FDonly     - Donor only emission calculated using the unmixed quantities
% Eapp       - Apparent FRET efficiency on the pixel level
%------------------------------------------------------------------------------
% Written by Gabriel biener, Ph.D on Thursday 10/24/2013
% Prof. Valerica Raicu Group
% Physics Department
% University of Wisconsin, Milwaukee, WI 53211
%------------------------------------------------------------------------------

NoBG_Image   = Image-BackGround; % image without the bias background
xx           = size(Image,1);    % x scale of the image
yy           = size(Image,2);    % y scale of the image
Wd           = sum(dSpectrum);   % Donor spectral integral
Wa           = sum(aSpectrum);   % Acceptor spectral integral


FMxFD        = zeros(xx,yy);     % Projection of the measurements on the donor spectrum
FMxFA        = zeros(xx,yy);     % % Projection of the measurements on the acceptor spectrum
FMxF3        = zeros(xx,yy);     % Projection of the measurements on the third component spectrum
for jj = 1:size(Image,3)
    FMxFD = FMxFD+squeeze(NoBG_Image(:,:,jj)).*dSpectrum(jj);
    FMxFA = FMxFA+squeeze(NoBG_Image(:,:,jj)).*aSpectrum(jj);
    FMxF3 = FMxF3+squeeze(NoBG_Image(:,:,jj)).*tSpectrum(jj);
end;
        
FAxFD        = sum(dSpectrum.*aSpectrum); % Integrating the cross product of the donor spectrum and the acceptor spectrum
FAxF3        = sum(aSpectrum.*tSpectrum); % Integrating the cross product of the acceptor spectrum and the third component spectrum
FDxF3        = sum(dSpectrum.*tSpectrum); % Integrating the cross product of the donor spectrum and the third component spectrum
FD2          = sum(dSpectrum.^2);         % Integrating the donor spectrum powered by 2
FA2          = sum(aSpectrum.^2);         % Integrating the acceptor spectrum powered by 2
F32          = sum(tSpectrum.^2);         % Integrating the third components spectrum powered by 2

% Calculating the unmixed quantities
KAD          = (((FMxFD*((FAxF3*FDxF3)-(F32*FAxFD)))+(FMxFA*((FD2*F32)-(FDxF3*FDxF3)))+(FMxF3*((FAxFD*FDxF3)-(FD2*FAxF3))))/((FD2*FA2*F32)+(2*FAxFD*FDxF3*FAxF3)-(FD2*FAxF3*FAxF3)-(FA2*FDxF3*FDxF3)-(F32*FAxFD*FAxFD)));
KDA          = (((FMxFD*((FA2*F32)-(FAxF3*FAxF3)))+(FMxFA*((FDxF3*FAxF3)-(F32*FAxFD)))+(FMxF3*((FAxFD*FAxF3)-(FA2*FDxF3))))/((FD2*FA2*F32)+(2*FAxFD*FDxF3*FAxF3)-(FD2*FAxF3*FAxF3)-(FA2*FDxF3*FDxF3)-(F32*FAxFD*FAxFD)));
K3           = (((FMxFD*((FAxFD*FAxF3)-(FA2*FDxF3)))+(FMxFA*((FAxFD*FDxF3)-(FD2*FAxF3)))+(FMxF3*((FD2*FA2)-(FAxFD*FAxFD))))/((FD2*FA2*F32)+(2*FAxFD*FDxF3*FAxF3)-(FD2*FAxF3*FAxF3)-(FA2*FDxF3*FDxF3)-(F32*FAxFD*FAxFD)));
% Gettrin rid of negative values
KAD(KAD < 0) = 0;
KDA(KDA < 0) = 0;
K3(K3 < 0)   = 0;

% Thresholding the unmixed images
K3(K3   <= Threshold) = 0;
KAD(KDA <= Threshold) = 0;
KAD(KAD <= Threshold) = 0;
KDA(KDA <= Threshold) = 0;
KDA(KAD <= Threshold) = 0;

FDonly                = KDA*Wd+KAD*Wa*dQY/aQY;
Eapp                  = KAD./(KAD + KDA.*aQY/dQY*Wd/Wa)*100;