%
%
%  Auto Crop 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="Crop Borders" icon="Matlab" tooltip="Automatically crop borders">
%          <Command>Matlab::BPAutoCrop(%i)</Command>
%        </Item>
%      </Menu>
%    </CustomTools>
%  
%
%  Description:
%
%   This XTension computes the bounding box of voxels with an intensity
%       higher than the specified threshold, and cuts away the "useless"
%       voxels at the borders.
%   This helps to minimize the size of the dataset without loosing 
%       significant data.
%
%

function BPAutoCrop(aImarisApplicationID)

[vSelection,vOk] = listdlg('ListString',{'X','Y','Z'}, ...
    'InitialValue', (1:3), 'SelectionMode','multiple', ...
    'ListSize',[300 100], 'Name','Crop Borders', ...
    'PromptString',{'Please select cropping direction:'});
if vOk<1, return; end

vAnswer = inputdlg(['Voxels with intensity less or equal to the threshold ', ...
    'are considered empty. Please enter the threshold value:'],...
    'Crop Borders',1,{'5'});
if isempty(vAnswer), return; end;
vThreshold = str2double(char(vAnswer));

% initialize cropping directions...    
vCropDirection = [0,0,0];
vCropDirection(vSelection) = 1;

vCropDirectionX = vCropDirection(1);
vCropDirectionY = vCropDirection(2);
vCropDirectionZ = vCropDirection(3);

%Get the application
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 dataset
vImarisDataSet = vImarisApplication.mDataSet.Clone;

% Get DataSet parameters
vNumberOfSlices = vImarisDataSet.mSizeZ;
vNumberOfTimePoints = vImarisDataSet.mSizeT;
vNumberOfChannels = vImarisDataSet.mSizeC;

vSizeX = vImarisDataSet.mSizeX;
vSizeY = vImarisDataSet.mSizeY;
vSizeZ = vImarisDataSet.mSizeZ;

vCurrentFirstSliceIndex = [vSizeX;vSizeY;vSizeZ];
vCurrentLastSliceIndex = [1;1;1];

%--------------------------------------------------------------------------
% Get the DataSet volume intensities and find the new extends for
% further cropping
vCount = 0;
vTotalCount = vNumberOfSlices*vNumberOfTimePoints*vNumberOfChannels;
vProgressDisplay = waitbar(0,'Getting Volume...');
for vChannel = 1:vNumberOfChannels        
    for vTimePoint = 1:vNumberOfTimePoints 
        vDataVolume = zeros(vSizeX,vSizeY,vSizeZ);
        for vSlice = 1:vNumberOfSlices
            vDataVolume(:,:,vSlice) = double(vImarisDataSet.GetDataSlice(vSlice-1,vChannel-1,vTimePoint-1));
            vCount = vCount + 1;
            waitbar(vCount/vTotalCount);
        end

        vFirstSliceIndex = [1;1;1];
        vLastSliceIndex = [vSizeX;vSizeY;vSizeZ];
        if vCropDirectionX == 1            
            while isempty(find(vDataVolume(vFirstSliceIndex(1),:,:)>vThreshold, 1)) && ...
                    vFirstSliceIndex(1)<vCurrentFirstSliceIndex(1)
                vFirstSliceIndex(1) = vFirstSliceIndex(1) + 1;
            end
            while isempty(find(vDataVolume(vLastSliceIndex(1),:,:)>vThreshold, 1)) && ...
                    vLastSliceIndex(1)>vCurrentLastSliceIndex(1)
                vLastSliceIndex(1) = vLastSliceIndex(1) - 1;
            end
        end

        if vCropDirectionY == 1                
            while isempty(find(vDataVolume(:,vFirstSliceIndex(2),:)>vThreshold, 1)) && ...
                    vFirstSliceIndex(2)<vCurrentFirstSliceIndex(2)
                vFirstSliceIndex(2) = vFirstSliceIndex(2) + 1;
            end
            while isempty(find(vDataVolume(:,vLastSliceIndex(2),:)>vThreshold, 1)) && ...
                    vLastSliceIndex(2)>vCurrentLastSliceIndex(2)
                vLastSliceIndex(2) = vLastSliceIndex(2) - 1;
            end
        end

        if vCropDirectionZ == 1                
            while isempty(find(vDataVolume(:,:,vFirstSliceIndex(3))>vThreshold, 1)) && ...
                    vFirstSliceIndex(3)<vCurrentFirstSliceIndex(3)
                vFirstSliceIndex(3) = vFirstSliceIndex(3) + 1;
            end
            while isempty(find(vDataVolume(:,:,vLastSliceIndex(3))>vThreshold, 1)) && ...
                    vLastSliceIndex(3)>vCurrentLastSliceIndex(3)
                vLastSliceIndex(3) = vLastSliceIndex(3) - 1;
            end
        end

        vCurrentFirstSliceIndex = vFirstSliceIndex;
        vCurrentLastSliceIndex = vLastSliceIndex;            
    end        
end
clear vDataVolume;

%--------------------------------------------------------------------------
% Send back DataVolume to Imaris

waitbar(1,vProgressDisplay,'Setting New Extends...');

% prevent the last index of being greater than the first
vCurrentLastSliceIndex = max([vCurrentLastSliceIndex,vCurrentFirstSliceIndex],2);

% resize the current Imaris DataSet in order to fit the new one
vImarisDataSet.Resize(vCurrentFirstSliceIndex(1)-1,...
                      vCurrentLastSliceIndex(1)-vCurrentFirstSliceIndex(1)+1,...
                      vCurrentFirstSliceIndex(2)-1,...
                      vCurrentLastSliceIndex(2)-vCurrentFirstSliceIndex(2)+1,...
                      vCurrentFirstSliceIndex(3)-1,...
                      vCurrentLastSliceIndex(3)-vCurrentFirstSliceIndex(3)+1,...
                      0,vNumberOfChannels,0,vNumberOfTimePoints);

close(vProgressDisplay);

vImarisApplication.mDataSet = vImarisDataSet;