%
%
%  Normalize Layers 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="Normalize Layers" icon="Matlab" tooltip="Normalizes the voxel intensities layer by layer.">
%          <Command>Matlab::BPNormalizeLayers(%i)</Command>
%        </Item>
%      </Menu>
%    </CustomTools>
%  
%
%  Description:
%   
%   Normalizes the voxel intensities layer by layer (same mean and stddev).
%
%

function BPNormalizeLayers(aImarisApplicationID, aExcludeZero)

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

vThreshold = 0;
if nargin == 1 
   vAnswer = inputdlg({'Exclude Voxels less than'}, 'Normalize Layers', 1, {'1'});
   if length(vAnswer) > 0
       vThreshold = str2double(vAnswer(1)); 
   else 
       return;
   end
elseif nargin == 2
    if aExcludeZero == true
        vThreshold = 1;
    end
end

% get the data set 
vImarisDataSet = vImarisApplication.mDataSet.Clone;
vImarisApplication.DataSetPushUndo('Normalize Layers');

vNumberOfTimePoints = vImarisDataSet.mSizeT;
vNumberOfChannels = vImarisDataSet.mSizeC;
vNumberOfChannelsSelected = 0;
for vChannelIndex = 1:vNumberOfChannels
    if vImarisApplication.GetChannelVisibility(vChannelIndex-1) == true
        vNumberOfChannelsSelected = vNumberOfChannelsSelected + 1;
    end
end
vNumberOfSlices = vImarisDataSet.mSizeZ;

vTotalCount = vNumberOfSlices*vNumberOfChannelsSelected*vNumberOfTimePoints;
vProgressDisplay = waitbar(0,'Normalize Layers');
vProgressCount = 0;

for vChannelIndex = 1:vNumberOfChannels
    if vImarisApplication.GetChannelVisibility(vChannelIndex-1) == false
        continue;
    end
    
    for vTimeIndex = 1:vNumberOfTimePoints
        vImage = zeros(vImarisDataSet.mSizeX,vImarisDataSet.mSizeY,vImarisDataSet.mSizeZ);
        for vSlice = 1:vImarisDataSet.mSizeZ
            vImage(:,:,vSlice) = vImarisDataSet.GetDataSlice( ...
                vSlice-1,vChannelIndex-1,vTimeIndex-1);
        end
        
        vDataVector = vImage(vImage>=vThreshold);
        vCount = numel(vDataVector);
        if vCount > 0
            vMean = mean(vDataVector);
            vVariance = var(vDataVector);
            for vSlice = 1:size(vImage,3)
                vImageSlice = vImage(:,:,vSlice);
                vDataVectorSlice = vImageSlice(vImageSlice>=vThreshold);
                vCountSlice = numel(vDataVectorSlice);
                if vCountSlice > 0
                    vMeanSlice = mean(vDataVectorSlice);
                    vVarianceSlice = var(vDataVectorSlice);
                    if vVarianceSlice == 0;
                        vVarianceSlice = vVariance;
                    end
                    vImageSlice(vImageSlice>=vThreshold) = vMean + ...
                      (vDataVectorSlice-vMeanSlice)*sqrt(vVariance/vVarianceSlice);
                    vImage(:,:,vSlice) = vImageSlice;
                end
            vProgressCount = vProgressCount + 1;
            waitbar(vProgressCount/vTotalCount)
            end
        else
            vProgressCount = vProgressCount + vNumberOfSlices;
            waitbar(vProgressCount/vTotalCount)
        end
        
        if strcmp(vImarisApplication.mDataSet.mType,'eTypeUInt8') 
            vImage = uint8(vImage);
        elseif strcmp(vImarisApplication.mDataSet.mType,'eTypeInt16') 
            vImage = uint16(vImage);
        else strcmp(vImarisApplication.mDataSet.mType,'eTypeFloat')    
            vImage = single(vImage);
        end
        vImarisDataSet.SetDataVolume(vImage,vChannelIndex-1,vTimeIndex-1);
    end
end
close(vProgressDisplay);

vImarisApplication.mDataSet = vImarisDataSet;