%
%
%  Distance Transformation 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>
%       <Submenu name="Surface Functions">
%        <Item name="Distance Transformation" icon="Matlab" tooltip="Voxel intensity of new channel indicates the shortest distance to a surface object.">
%          <Command>Matlab::BPDistanceTransform(%i)</Command>
%        </Item>
%       </Submenu>
%      </Menu>
%      <SurpassTab>
%        <SurpassComponent name="bpSurfaceObject">
%          <Item name="Distance Transformation" icon="Matlab" tooltip="Distance Transformation">
%            <Command>Matlab::BPDistanceTransform(%i)</Command>
%          </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
%  
%
%  Description:
%
%   Create and compute a new channel, where the voxel intensity indicates
%       the shortest distance to a surface object.
%
%

function BPDistanceTransform(aImarisApplicationID)

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

% Get Surpass splitted SurfaceObject
vImarisSurfaceObject = vImarisApplication.mSurpassSelection;            
% Check if there is a selection
if isempty(vImarisSurfaceObject)
    msgbox('A SurfaceObject must be selected');
    return;
end
% Check if the selection is a Splitted SurfaceObject
if ~ismethod(vImarisSurfaceObject, 'GetMask')
    msgbox('Your selection is not a SurfaceObject');
    return;
end

% if the DataSet is not of type float, display a warning. As the result of
% DistanceTransform will very probably be floating values, a DataSet of
% type int will not be able to display it correctly.
if strcmp(vImarisDataSet.mType,'eTypeUInt8') || strcmp(vImarisDataSet.mType,'eTypeUInt16')
    vWarningAnswer = warndlg('Due to the DataSet type the results of the Distance Transform function will probably be truncated! It is best to use it with a float DataSet!','!Warning!');
    waitfor(vWarningAnswer);
end

vImarisApplication.DataSetPushUndo('Distance Transform');

% Draw the input dialog and ask for a threshold.    
[vSelection,vOk] = listdlg('ListString',{'Inside SurfaceObject','Outside SurfaceObject'},...
                       'SelectionMode','single',...
                       'ListSize',[200 50],...
                       'Name','BPDistanceTransform - Selection Box',...
                       'PromptString',{'Please select Distance Transform mode:'});
% Go on if the threshold value exists.
if vOk < 1, return; end;
        
% vNumberOfTimePoints = vImarisDataSet.mSizeT;
vNumberOfChannels = vImarisDataSet.mSizeC;
vNumberOfSlices = vImarisDataSet.mSizeZ;

% Create the 'weigths' matrix - corresponds to the typical distance
% between a voxel and its neighbours
vDeltaX = (vImarisDataSet.mExtendMaxX-vImarisDataSet.mExtendMinX)/vImarisDataSet.mSizeX;
vDeltaY = (vImarisDataSet.mExtendMaxY-vImarisDataSet.mExtendMinY)/vImarisDataSet.mSizeY;
vDeltaZ = (vImarisDataSet.mExtendMaxZ-vImarisDataSet.mExtendMinZ)/vImarisDataSet.mSizeZ;

vXYDiag = sqrt(vDeltaX^2 + vDeltaY^2);
vXYZDiag = sqrt(vXYDiag^2 + vDeltaZ^2);
vXZDiag = sqrt(vDeltaX^2 + vDeltaZ^2);
vYZDiag = sqrt(vDeltaY^2 + vDeltaZ^2);

vWeigths(:,:,1) = [vXYZDiag vYZDiag vXYZDiag;vXZDiag vDeltaZ vXZDiag;vXYZDiag vYZDiag vXYZDiag];
vWeigths(:,:,2) = [vXYDiag vDeltaY vXYDiag;vDeltaX 0 vDeltaX;vXYDiag vDeltaY vXYDiag];
vWeigths(:,:,3) = vWeigths(:,:,1);

vProgressDisplay = waitbar(0,'Distance Transform: Getting Mask');

% Get the mask DataSet
vMaskDataSet = vImarisSurfaceObject.GetMask( ...
    vImarisDataSet.mExtendMinX,vImarisDataSet.mExtendMinY,vImarisDataSet.mExtendMinZ, ...
    vImarisDataSet.mExtendMaxX,vImarisDataSet.mExtendMaxY,vImarisDataSet.mExtendMaxZ, ...
    vImarisDataSet.mSizeX,vImarisDataSet.mSizeY,vImarisDataSet.mSizeZ);

% Get the volume and compute the distance transformation for each timepoint
waitbar(0,vProgressDisplay,'Distance Transform: Getting Volume');
% Create a new channel where the result will be sent
vImarisDataSet.mSizeC = vNumberOfChannels + 1;
vImarisDataSet.SetChannelName(vNumberOfChannels,['Distance to ',vImarisSurfaceObject.mName]);
vImarisDataSet.SetChannelColor(vNumberOfChannels,0,0,1,0);

vDataVolume = zeros(vImarisDataSet.mSizeX,vImarisDataSet.mSizeY,vImarisDataSet.mSizeZ);

vTime = vImarisSurfaceObject.GetTimeIndex;

for vSlice = 1:vNumberOfSlices
    %Get the volume
    if vSelection == 1
        vDataVolume(:,:,vSlice) = vMaskDataSet.GetDataSlice(vSlice-1,0,0)~=1;
    else
        vDataVolume(:,:,vSlice) = vMaskDataSet.GetDataSlice(vSlice-1,0,0)==1;
    end
end

waitbar(0.2,vProgressDisplay,'Distance Transform: Processing Transformation');
% Creates a connectivity matrix (define which voxel are
% considered as neighbours - typically a 3x3x3 matrix filled with ones)
vConnectivity = ones(3,3,3);
% Use the matlab function with customized weight matrix

vDistanceTransformImage = ddist(int8(vDataVolume),vConnectivity,vWeigths);
%vDistanceTransformImage = bwdist(vDataVolume);

waitbar(0.8,vProgressDisplay,'Distance Transform: Setting Volume');
% Set Distance transform image back to Imaris

if strcmp(vImarisDataSet.mType,'eTypeUInt8')
    vImarisDataSet.SetDataVolume(uint8(vDistanceTransformImage),vNumberOfChannels,vTime);
elseif strcmp(vImarisDataSet.mType,'eTypeUInt16')
    vImarisDataSet.SetDataVolume(uint16(vDistanceTransformImage),vNumberOfChannels,vTime);
elseif strcmp(vImarisDataSet.mType,'eTypeFloat')
    vImarisDataSet.SetDataVolume(single(vDistanceTransformImage),vNumberOfChannels,vTime);
end

close(vProgressDisplay);

vImarisApplication.mDataSet = vImarisDataSet;