%
%
%  Stereology Box for Imaris
%
%  Copyright Bitplane AG 2006
%
%
%  Installation:
%
%  - Copy this file into the XTensions subdirectory of the Imaris
%    Installation directory (create this directory if necessary)
%  - Start Imaris
%  - You will now find this module as a button in surpass.
%
%  Do not edit the xml below:
%
%    <CustomTools>
%      <Menu>
%        <Item name="Stereology Box" icon="Matlab" tooltip="Create a surface object with the shape of a cube.">
%          <Command>Matlab::BPStereologyBox(%i)</Command>
%        </Item>
%      </Menu>
%      <SurpassToolbar>
%        <Item name="Stereology Box" tooltip="Create a surface object with the shape of a cube." icon="Matlab">
%          <Command>Matlab::BPStereologyBox(%i)</Command>
%        </Item>
%      </SurpassToolbar>
%    </CustomTools>
%  
% 
%  Description:
% 
%   Create a surface object with the shape of a cube. 
%   The edge length is defined by the user.
%
%

function BPStereologyBox(aImarisApplicationID)

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

vDataSet = vImarisApplication.mDataSet.Clone;

if vDataSet.mSizeX == 0
  msgbox('Please Load a DataSet first');
  return;
end

vCreateBoxChannel = false;

[vExtendsMin, vExtendsMax] = BPGetExtends_(vDataSet);
vVoxelSize = BPGetVoxelSize_(vDataSet);
vCenter = (vExtendsMin + vExtendsMax)/2;
vCenter(3) = vExtendsMax(3)+vVoxelSize(3)/2;
vLength = min(vCenter-vExtendsMin);
vAnswer = inputdlg({'Length of Box'}, 'Stereology Box', 1, {num2str(vLength)});
if length(vAnswer) > 0
  vLength = str2double(vAnswer(1));
  vLength = [vLength, vLength, vLength/2];

  if vDataSet.mSizeT > 1
    vGroup = vImarisApplication.mFactory.CreateDataContainer();
    vImarisApplication.mSurpassScene.AddChild(vGroup);
    vGroup.mName = 'StereologyBox';
  else 
    vGroup = vImarisApplication.mSurpassScene;
  end
  
  if vCreateBoxChannel
    vChannelIndexStereologyBox = GetChannelIDByName(vImarisApplication, 'StereologyBox');
    if vChannelIndexStereologyBox < 0
      vDataSet.mSizeC = vDataSet.mSizeC + 1;
      vChannelIndexStereologyBox = vDataSet.mSizeC;
      vDataSet.SetChannelName(vChannelIndexStereologyBox-1,'StereologyBox');
    end
    % also color channel
    vImage = vDataSet.GetDataVolumeAs1DArray(vChannelIndexStereologyBox-1, 0);
    vImage = reshape(vImage, [vDataSet.mSizeX, vDataSet.mSizeY, vDataSet.mSizeZ]);

    vCenterIndex = vExtendsMin + vCenter./vVoxelSize;
    vLower = floor(vCenterIndex - 0.5 * vLength./vVoxelSize);
    vLower = max(vLower, [1,1,1]);
    vUpper = ceil(vCenterIndex + 0.5 * vLength./vVoxelSize);
    vSize = [1,1,1];
    vSize(1:ndims(vImage)) = size(vImage);
    vUpper = min(vUpper, vSize);
    vRange{1} = vLower(1):vUpper(1);
    vRange{2} = vLower(2):vUpper(2);
    vRange{3} = vLower(3):vUpper(3);
    vRangeInner{1} = (vLower(1)+1):(vUpper(1)-1);
    vRangeInner{2} = (vLower(2)+1):(vUpper(2)-1);
    vRangeInner{3} = (vLower(3)+1):(vUpper(3));
    vImage = uint8(zeros(size(vImage)));
    vImage(vRange{:}) = 1;
    vImage(vRangeInner{:}) = 0;
  end
  
  for vTimeIndex = 0:(vDataSet.mSizeT-1)
    vSurfaceObject = vImarisApplication.mFactory.CreateSurface();
    vVertices = [[0,0,0];[1,0,0];[0,1,0];[1,1,0]; [0,0,1];[1,0,1];[0,1,1];[1,1,1]; ...
                 [0,0,0];[1,0,0];[0,1,0];[1,1,0]; [0,0,1];[1,0,1];[0,1,1];[1,1,1]; ... 
                 [0,0,0];[1,0,0];[0,1,0];[1,1,0]; [0,0,1];[1,0,1];[0,1,1];[1,1,1]];
    vNormals = [ [0,0,1];[0,0,1];[0,0,1];[0,0,1];[0,0,1];[0,0,1];[0,0,1];[0,0,1]; ...
                 [0,-1,0];[0,-1,0];[0,-1,0];[0,-1,0];[0,-1,0];[0,-1,0];[0,-1,0];[0,-1,0]; ...
                 [1,0,0];[1,0,0];[1,0,0];[1,0,0];[1,0,0];[1,0,0];[1,0,0];[1,0,0]];
    % vSquares = [[1,2,3,4]; [5,6,7,8]; [1,2,5,6]; [3,4,7,8]; [1,3,5,7]; [2,4,6,8]];
    vTriangles = [[1,2,3]; [2,4,3]; ...  % [5,6,7];[6,8,7];
      [1,2,5]+8;[2,6,5]+8;[3,4,7]+8;[4,8,7]+8;[1,3,5]+16;[3,7,5]+16;[2,4,6]+16;[4,8,6]+16];

    vVertices = vVertices - 0.5;
    vVertices(:,3) = (vVertices(:,3)-0.5);
    vVertices(:,1) = vVertices(:,1) * vLength(1);
    vVertices(:,2) = vVertices(:,2) * vLength(2);
    vVertices(:,3) = vVertices(:,3) * vLength(3);
    vVertices(:,1) = vVertices(:,1) + vCenter(1);
    vVertices(:,2) = vVertices(:,2) + vCenter(2);
    vVertices(:,3) = vVertices(:,3) + vCenter(3);

    vSurfaceObject.SetSurface(vVertices, vTriangles-1, vNormals, vTimeIndex);
    vSurfaceObject.mName = strcat('StereologyBox', num2str(vTimeIndex));
    vGroup.AddChild(vSurfaceObject);
  
    if vCreateBoxChannel
      vDataSet.SetDataVolume(vImage, vChannelIndexStereologyBox-1, vTimeIndex);
    end
  end
end

vImarisApplication.mDataSet = vDataSet;

%%
function aChannelID = GetChannelIDByName(aImarisApplication, aChannelName)
aChannelID = -1;
for vChannelID = 1:aImarisApplication.mDataSet.mSizeC
  if strcmp(aImarisApplication.mDataSet.GetChannelName(vChannelID-1), aChannelName)
    aChannelID = vChannelID;
  end
end

%%
% function GetExtends
function [aExtendMin, aExtendMax] = BPGetExtends_(aImarisDataSet)

aExtendMin = [0,0,0];
aExtendMin(1)= aImarisDataSet.mExtendMinX;
aExtendMin(2) = aImarisDataSet.mExtendMinY;
aExtendMin(3) = aImarisDataSet.mExtendMinZ;

aExtendMax = [0,0,0];
aExtendMax(1)= aImarisDataSet.mExtendMaxX;
aExtendMax(2) = aImarisDataSet.mExtendMaxY;
aExtendMax(3) = aImarisDataSet.mExtendMaxZ;

%%
function aVoxelSize = BPGetVoxelSize_(aImarisDataSet)
aVoxelSize = [1,1,1];
aVoxelSize(1) = (aImarisDataSet.mExtendMaxX-aImarisDataSet.mExtendMinX)/aImarisDataSet.mSizeX;
aVoxelSize(2) = (aImarisDataSet.mExtendMaxY-aImarisDataSet.mExtendMinY)/aImarisDataSet.mSizeY;
aVoxelSize(3) = (aImarisDataSet.mExtendMaxZ-aImarisDataSet.mExtendMinZ)/aImarisDataSet.mSizeZ;