%
%
%  Time Subtract Average 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="TimeSubtractAverage" icon="Matlab" tooltip="For each voxel, subtract the temporal average intensity.">
%          <Command>Matlab::BPTimeSubtractAverage(%i)</Command>
%        </Item>
%      </Menu>
%    </CustomTools>
%  
%
%  Description:
%   
%   For each voxel, subtract the temporal average intensity.
%   This function is used to eliminate motionless objects.
%   For 2D time series try to use 'Image processing - Swap Time and Z' 
%     before and after starting the XTension to speed up processing.
%
%

function BPTimeSubtractAverage(aImarisApplicationID)

% get the application object
if isa(aImarisApplicationID, 'COM.Imaris_Application')
  % called from workspace
  vImarisApplication = aImarisApplicationID;
else
  % connect to Imaris Com interface
  vImarisServer = actxserver('ImarisServer.Server');
  if ischar(aImarisApplicationID)
    aImarisApplicationID = round(str2double(aImarisApplicationID));
  end
  vImarisApplication = vImarisServer.GetObject(aImarisApplicationID);
end

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

try
  vImarisApplication.DataSetPushUndo('TimeSubtractAverage');
catch
  % nothing to do
end

try % dataset has been loaded?
  vSize = [vImarisDataSet.mSizeX, vImarisDataSet.mSizeY, vImarisDataSet.mSizeZ];
catch
  msgbox('Invalid dataset')
  return
end

vTotalCount = vImarisDataSet.mSizeC * vImarisDataSet.mSizeT * 3;

vIsTime = 1;
if vImarisDataSet.mSizeT==1
  vAnswer = questdlg(['This is not a time dataset. Would you like to ', ...
    'subtract the average of the slices instead?'], 'TimeSubtractAverage');
  if ~isequal(vAnswer, 'Yes')
    return
  end
  vIsTime = 0;

  vAtOnce = inputdlg('Please enter the number of slices to process at once:', ...
    'TimeSubtractAverage', 1, {'100'});
  vAtOnce = str2double(vAtOnce{1});
  if isnan(vAtOnce)
    msgbox('Invalid entry')
    return
  end
  vTotalCount = vTotalCount * numel(0:vAtOnce:vSize(3)-1);
end

[vDiscardMode, vOk] = listdlg('ListString', {'Set to zero', 'Get absolute value'},...
                       'SelectionMode','single',...
                       'ListSize',[200 50],...
                       'Name','TimeSubtractAverage',...
                       'PromptString', {'Neagative values:'});
if vOk<1
    return;
end

vProgressDisplay = waitbar(0, 'TimeSubtractAverage');
vCount = 0;

try
  for vChannel = 0:vImarisDataSet.mSizeC-1
    if vImarisApplication.GetChannelVisibility(vChannel)
      if vIsTime
        % compute mean volume
        vMeanVolume = double(zeros(1, prod(vSize)));
        for vTime = 0:vImarisDataSet.mSizeT-1
          vVolume = double(vImarisDataSet.GetDataVolumeAs1DArray(vChannel, vTime));
          vMeanVolume = (vMeanVolume*vTime + vVolume) / (vTime+1);
          vCount = vCount + 1;
          waitbar(vCount/vTotalCount)
        end

        % subtract mean volume
        vMin = zeros(vImarisDataSet.mSizeT, 1) + inf;
        vMax = zeros(vImarisDataSet.mSizeT, 1) - inf;
        for vTime = 0:vImarisDataSet.mSizeT-1
          vVolume = double(vImarisDataSet.GetDataVolumeAs1DArray(vChannel, vTime));
          vVolume = reshape(vVolume - vMeanVolume, vSize);

          if vDiscardMode
            vVolume(vVolume<0) = 0;
          else
            vVolume = abs(vVolume);
          end

          vMin(vTime+1) = min(vVolume(:));
          vMax(vTime+1) = max(vVolume(:));
          if strcmp(vImarisDataSet.mType,'eTypeUInt8')
            vImarisDataSet.SetDataVolume(uint8(vVolume), vChannel, vTime);
          elseif strcmp(vImarisDataSet.mType,'eTypeUInt16')
            vImarisDataSet.SetDataVolume(uint16(vVolume), vChannel, vTime);
          elseif strcmp(vImarisDataSet.mType,'eTypeFloat')
            vImarisDataSet.SetDataVolume(single(vVolume), vChannel, vTime);
          end
          vCount = vCount + 2;
          waitbar(vCount/vTotalCount)
        end
        
      else % vIstime == false
        
        % compute mean volume
        vMeanVolume = double(zeros(vSize(1:2)));
        for vBlock = 0:vAtOnce:vSize(3)-1
          vMinZ = vBlock;
          vSizeZ = min(vBlock + vAtOnce, vSize(3)) - vBlock;
          vVolume = double(vImarisDataSet.GetDataSubVolumeAs1DArray( ...
            0, 0, vMinZ, vChannel, 0, vSize(1), vSize(2), vSizeZ));
          
          vVolume = reshape(vVolume, vSize(1), vSize(2), []);

          vMeanVolume = (vMeanVolume*vMinZ + mean(vVolume, 3)*vSizeZ) / ...
            (vMinZ + vSizeZ);
          
          vCount = vCount + 1;
          waitbar(vCount/vTotalCount)
        end
        
        % subtract mean volume
        vMin = zeros(numel(0:vAtOnce:vSize(3)-1), 1) + inf;
        vMax = zeros(numel(0:vAtOnce:vSize(3)-1), 1) - inf;
        for vBlock = 0:vAtOnce:vSize(3)-1
          vMinZ = vBlock;
          vSizeZ = min(vBlock + vAtOnce, vSize(3)) - vBlock;
          vVolume = double(vImarisDataSet.GetDataSubVolumeAs1DArray( ...
            0, 0, vMinZ, vChannel, 0, vSize(1), vSize(2), vSizeZ));
          
          vVolume = reshape(vVolume, [vSize(1:2), vSizeZ]);
          for vSlice = 1:vSizeZ
            vVolume(:, :, vSlice) = vVolume(:, :, vSlice) - vMeanVolume;
          end

          if vDiscardMode
            vVolume(vVolume<0) = 0;
          else
            vVolume = abs(vVolume);
          end

          vMin(vBlock+1) = min(vVolume(:));
          vMax(vBlock+1) = max(vVolume(:));
          if strcmp(vImarisDataSet.mType,'eTypeUInt8')
            vVolume = uint8(vVolume);
          elseif strcmp(vImarisDataSet.mType,'eTypeUInt16')
            vVolume = uint16(vVolume);
          elseif strcmp(vImarisDataSet.mType,'eTypeFloat')
            vVolume = single(vVolume);
          end
          
          vImarisDataSet.SetDataSubVolume(vVolume, 0, 0, vMinZ, vChannel, 0);
          
          vCount = vCount + 2;
          waitbar(vCount/vTotalCount)
        end

      end
      vImarisDataSet.SetChannelRange(vChannel, min(vMin), max(vMax));
    end
  end
  close(vProgressDisplay)
catch
  close(vProgressDisplay);
  rethrow(lasterror)
end

vImarisApplication.mDataSet = vImarisDataSet;