%
%
%  Stereology SpaceBall 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 Space Ball" icon="Matlab" tooltip="Create a surface object with the shape of a half shsphere.">
%          <Command>Matlab::BPStereologySpaceBall(%i)</Command>
%        </Item>
%      </Menu>
%      <SurpassToolbar>
%        <Item name="Stereology Space Ball" tooltip="Create a surface object with the shape of a half sphere." icon="Matlab">
%          <Command>Matlab::BPStereologySpaceBall(%i)</Command>
%        </Item>
%      </SurpassToolbar>
%    </CustomTools>
%  
%
%  Description:
% 
%   Create a surface object with the shape of a half sphere. 
%   The radius size is defined by the user.
%
%

function BPStereologySpaceBall(aImarisApplicationID)

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

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

vSpots = vImarisApplication.mFactory.CreateSpots();
[vExtendsMin, vExtendsMax] = BPGetExtends_(vImarisApplication.mDataSet);
vPosition = (vExtendsMin + vExtendsMax)/2;
vPosition(3) = vExtendsMax(3);
vRadius = min(vPosition-vExtendsMin);
vAnswer = inputdlg({'Radius of SpaceBall'}, 'Stereology Box', 1, {num2str(vRadius)});
if length(vAnswer) > 0
  vRadius= str2double(vAnswer(1));
  
  %%
  vCreateSpot = false;
  if vCreateSpot
    vSizeT= vImarisApplication.mDataSet.mSizeT;
    vPositions = repmat(vPosition, [vSizeT, 1]);
    vNSpots = vSizeT;
    vRadii = repmat(vRadius, [vNSpots, 1]);
    vTimeIndices = floor((0:(vNSpots-1))/(vNSpots/vSizeT));
    vSpots.Set(vPositions, vTimeIndices, vRadii);
    vSpots.mName = 'StereologySpaceBall';
    vImarisApplication.mSurpassScene.AddChild(vSpots);
  end

  vNTheta = 20;
  vVectors = CreateRadialVectors_(vNTheta);
  vVectors = vVectors(vVectors(:,3) <= 0,:);
  vNPhi = size(vVectors,1)/vNTheta;
  vVertices = vVectors * vRadius;
  vVertices(:,1) = vVertices(:,1) + vPosition(1);
  vVertices(:,2) = vVertices(:,2) + vPosition(2);
  vVertices(:,3) = vVertices(:,3) + vPosition(3);
  vNormals = vVectors;
  vTriangles = [[1,2,vNTheta+1]; [2,vNTheta+2,vNTheta+1]]; 
  vTriangles = repmat(vTriangles, [vNTheta-1,1]) + repmat(floor((0:(2*(vNTheta-1)-1))/2)',[1,3]);
  vTriangles = [vTriangles ; [vNTheta, 1, 2*vNTheta]; [1, vNTheta+1, 2*vNTheta]];
  vTriangles = repmat(vTriangles, [vNPhi-1,1]) + repmat(vNTheta*floor((0:(2*(vNTheta)*(vNPhi-1)-1))/(2*(vNTheta)))',[1,3]);
  if vImarisApplication.mDataSet.mSizeT > 1
    vGroup = vImarisApplication.mFactory.CreateDataContainer();
    vImarisApplication.mSurpassScene.AddChild(vGroup);
    vGroup.mName = 'StereologySpaceBall';
  else
    vGroup = vImarisApplication.mSurpassScene;
  end
  for vTimeIndex = 0:(vImarisApplication.mDataSet.mSizeT-1)
    vSurfaceObject = vImarisApplication.mFactory.CreateSurface();
    vSurfaceObject.SetSurface(vVertices, vTriangles-1, vNormals, vTimeIndex);
    vSurfaceObject.mName = strcat('SpaceBall ', num2str(vTimeIndex));
    vGroup.AddChild(vSurfaceObject);
  end
end


%%
% function aVectors = CreateRadialVectors(vNSamplesEquator)
%
function aVectors = CreateRadialVectors_(vNSamplesEquator)

if nargin == 0
  vNSamplesEquator = 20;
end

vPhi = (0:vNSamplesEquator)/vNSamplesEquator*pi;

vCoord = []; 
vNTheta = vNSamplesEquator;
for vCount = 1:numel(vPhi); 
  vNewTheta = (1:vNTheta)/vNTheta*2*pi; 
  vNew = [repmat(vPhi(vCount),[vNTheta,1]),reshape(vNewTheta,[vNTheta,1])];
  vCoord = [vCoord; vNew]; 
end;

[aVectors(:,1), aVectors(:,2), aVectors(:,3)] = sph2cart(vCoord(:,2), vCoord(:,1)-pi/2, ones(size(vCoord,1),1));

%%
% 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;