%
%
%  Attenuation Correction Function for Imaris
%
%  Copyright Bitplane AG 2006
%
%
%  Installation:
%
%  - Copy this file into the XTensions folder in the Imaris installation directory
%  - You will find this function in the Image Processing menu
%
%    <CustomTools>
%      <Menu>
%        <Item name="Attenuation Correction" icon="Matlab" tooltip="Attenuation correction">
%          <Command>Matlab::BPAttenuationCorrection(%i,'%c')</Command>
%        </Item>
%      </Menu>
%    </CustomTools>
%  
%
%  Description:
%
%   The layers to the top of the data set (close to objective lense) show
%       often much higher intensities than layers deep into the tissue .
%       This XTension can correct this attenuation effect.
%   The correction is performed along the Z axis.
%
%

function BPAttenuationCorrection(aImarisApplicationID, aChannelSelection)

% convert channel selection string into array
vChannelSelection = str2double(aChannelSelection);
vChannelIndex = find(vChannelSelection, 1, 'first');

if isa(aImarisApplicationID, 'COM.Imaris_Application')
    vImarisApplication = aImarisApplicationID;
else
    % connect to Imaris Com interface
    vImarisServer = actxserver('ImarisServer.Server');
    vImarisApplication = vImarisServer.GetObject(aImarisApplicationID);
end

% get the data set
vImarisDataSet = vImarisApplication.mDataSet.Clone;

vAnswer = inputdlg({'Channel', 'Intensity Front', 'Intensity Back'}, ...
    'Attenuation Correction', 1, {num2str(vChannelIndex), '256', '128'});
if isempty(vAnswer), return; end

vNumberOfTimePoints = vImarisDataSet.mSizeT;
vNumberOfSlices = vImarisDataSet.mSizeZ;

vDeltaZ = vNumberOfSlices;

vA = str2double(vAnswer(2));
vB = str2double(vAnswer(3));    
if vB <=0 || vA <= 0
    msgbox('BPAttenuationCorrection: Intensities must be greater than zero.');
    return;
end
vAlpha= (log(vB)-log(vA))/vDeltaZ;
vChannelIndex = round(str2double(vAnswer(1)));
if vChannelIndex <1 || vChannelIndex > vImarisDataSet.mSizeC
    msgbox('Channel index out of bound');
    return;
end

vImarisApplication.DataSetPushUndo('Attenuation Correction');

vTotalCount = vNumberOfSlices*vNumberOfTimePoints;
vProgressDisplay = waitbar(0,'Attenuation Correction');
vCount = 0;

try
    for vTime = 1:vNumberOfTimePoints
        for vSlice = 1:vNumberOfSlices
            vImageSlice = double(vImarisDataSet.GetDataSlice(vSlice-1, vChannelIndex-1, vTime-1));
            vImageSlice = vImageSlice*exp(-vAlpha*(vSlice-1));
            if strcmp(vImarisDataSet.mType,'eTypeUInt8')
                vImarisDataSet.SetDataSlice(uint8(vImageSlice), vSlice-1, vChannelIndex-1,vTime-1);
            elseif strcmp(vImarisDataSet.mType,'eTypeUInt16')
                vImarisDataSet.SetDataSlice(uint16(vImageSlice), vSlice-1, vChannelIndex-1,vTime-1);
            elseif strcmp(vImarisDataSet.mType,'eTypeFloat')
                vImarisDataSet.SetDataSlice(single(vImageSlice), vSlice-1, vChannelIndex-1,vTime-1);
            end
            vCount = vCount + 1;
            waitbar(vCount/vTotalCount)
        end
    end
    close(vProgressDisplay)
catch
    close(vProgressDisplay);
    rethrow(lasterror)
end

vImarisApplication.mDataSet = vImarisDataSet;