%
%
%  Normalize TimePoints 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 Time Points" icon="Matlab" tooltip="Normalizes the voxel intensities for each time point of a 4D dataset.">
%          <Command>Matlab::BPNormalizeTimePoints(%i)</Command>
%        </Item>
%      </Menu>
%    </CustomTools>
%  
%
%  Description:
%   
%   Normalizes the voxel intensities for each time point of a 4D dataset
%       (same mean and stddev).
%
%

function BPNormalizeTimePoints(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

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

if vNumberOfTimePoints == 1
    msgbox('There is only one time point', 'Info');
    return;
end

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

elseif nargin == 2
    if aExcludeZero == true
        vThreshold = 1;
    end
end

vImarisApplication.DataSetPushUndo('Normalize TimePoints');

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

vProgressTotalCount = 2*vNumberOfChannelsSelected*vNumberOfTimePoints;
vProgressDisplay = waitbar(0,'Normalize TimePoints');
vProgressCount = 0;

for vChannelIndex = 1:vNumberOfChannels
    if vImarisApplication.GetChannelVisibility(vChannelIndex-1) == false
        continue;
    end
    
    % get total mean and variance
    vCount = 0;
    vSum = 0;
    vSumOfSquares = 0;
    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);
        vCountVolume = numel(vDataVector);
        if vCountVolume > 0
            vCount= vCount + vCountVolume;
            vSum = vSum + sum(vDataVector);
            vSumOfSquares = vSumOfSquares + sum(vDataVector.^2);
        end
        vProgressCount = vProgressCount + 1;
        waitbar(vProgressCount/vProgressTotalCount)
    end
    vMean = vSum/vCount;
    vVariance = vSumOfSquares/vCount - vMean*vMean;

    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
            vMeanVolume = mean(vDataVector);
            vVarianceVolume = var(vDataVector);
            if vVarianceVolume == 0
                vVarianceVolume = vVariance;
            end
            vImage(vImage>=vThreshold) = (vDataVector-vMeanVolume)*sqrt(vVariance/vVarianceVolume) + vMean;
            
            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
        vProgressCount = vProgressCount + 1;
        waitbar(vProgressCount/vProgressTotalCount)
    end    
end
close(vProgressDisplay);

vImarisApplication.mDataSet = vImarisDataSet;