%
%
%  Spots Convex Hull Function for Imaris 9.0.0
%
%  Chi-Li Chiu 20180710
%  Modified from XTFilamentsConvexHull
%
%
%  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="Spots Functions">
%        <Item name="Spots Convex Hull" icon="Matlab" tooltip="Create a Surface which contains the convex hull of the Spots points.">
%          <Command>MatlabXT::XTSpotsConvexHull(%i)</Command>
%        </Item>
%       </Submenu>
%      </Menu>
%      <SurpassTab>
%        <SurpassComponent name="bpSpots">
%          <Item name="Convex Hull">
%            <Command>MatlabXT::XTSpotsConvexHull(%i)</Command>
%          </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
% 
%
%  Description:
%
%   Create a Surface which contains the convex hull of the Spots points.
%   Works with time lapse images. 
%   Does not work with 2D images.
%

function XTSpotsConvexHull(aImarisApplicationID)
  % connect to Imaris interface
  if ~isa(aImarisApplicationID, 'Imaris.IApplicationPrxHelper')
    javaaddpath ImarisLib.jar
    vImarisLib = ImarisLib;
    if ischar(aImarisApplicationID)
      aImarisApplicationID = round(str2double(aImarisApplicationID));
    end
    vImarisApplication = vImarisLib.GetApplication(aImarisApplicationID);
  else
    vImarisApplication = aImarisApplicationID;
  end

  % we need a dataset for generating the convex hull surface
  if isempty(vImarisApplication.GetDataSet())
    msgbox('An image is needed for this XTension');
    return;
  end

  % get the spot
  vFactory = vImarisApplication.GetFactory;
  vSpots = vFactory.ToSpots(vImarisApplication.GetSurpassSelection);

  % search the spot if not previously selected
  vSurpassScene = vImarisApplication.GetSurpassScene;
  if ~vFactory.IsSpots(vSpots)        
    for vChildIndex = 1:vSurpassScene.GetNumberOfChildren
      vDataItem = vSurpassScene.GetChild(vChildIndex - 1);
      if vFactory.IsSpots(vDataItem)
        vSpots = vFactory.ToSpots(vDataItem);
        break;
      end
    end
    % did we find the spot?
    if isequal(vSpots, [])
      msgbox('Please create some spot!');
      return;
    end
  end

  vDataSet = createWorkingDataSet(vImarisApplication);
  
  if isequal(vDataSet.GetSizeZ(),1)
      msgbox('Cannot work on 2D.');
      return;
  end
  
  vCountDegeneratedData = 0;
  % create result object
  vSurfaceHull = vFactory.CreateSurfaces;
  
  vSpotsXYZTotal = double(vSpots.Get.mPositionsXYZ);
  
  vSpotsTimeIndices = vSpots.GetIndicesT;
  maxT = max(vSpotsTimeIndices);
  
  %loop through time point
  for vIndexT = 0:maxT

    
    vSpotsT = find(vSpotsTimeIndices == vIndexT);
    if size(vSpotsT,1) < 4
        continue
    end
    % get the points coordinates
    vSpotsXYZ = vSpotsXYZTotal(vSpotsT,:);    

    try
      % yes, matlab does most of the work!
      vConvexHull = convhulln(vSpotsXYZ);
    catch er
      er.message; % suppress warning
      vCountDegeneratedData = vCountDegeneratedData + 1;
    end
    % select the necessary points
    vVertices = vSpotsXYZ(unique(vConvexHull(:)), :);
    vTriangulation = delaunayTriangulation(vVertices);
    [X, Y, Z] = getVoxelCenterCoordinates(vDataSet);
    vSimplexIndex = pointLocation(vTriangulation, X(:), Y(:), Z(:));
    vMask = ~isnan(vSimplexIndex);
    vVoxels = reshape(vMask, size(X));
    % create a dataset from convex hull mask and detect surface on it
    vOrigDataSet = vImarisApplication.GetDataSet();
    vSpotTimeIndex = vIndexT;
    vSpotTimepoint = vOrigDataSet.GetTimePoint(vSpotTimeIndex);
    vDataSet.SetTimePoint(0, vSpotTimepoint);

    for vIndexZ = 1:size(vVoxels,3)
      vDataSet.SetDataSliceFloats(vVoxels(:, :, vIndexZ), vIndexZ-1, 0, 0);
    end

    vSurface = DetectSurfaceWithoutSmoothing(vImarisApplication, vDataSet);
    vSurface.CopySurfacesToSurfaces(0, vSurfaceHull);

  end

  vSurfaceHull.SetName(['Convex Hull of ', char(vSpots.GetName)]);
  vSurfaceHull.SetColorRGBA(vSpots.GetColorRGBA);
  vSpots.GetParent.AddChild(vSurfaceHull, -1);

  if vCountDegeneratedData > 0
    msgbox(['Could not create convex hull for some of the Spots. ', ...
      'spot points must not lie on a plane to build a valid convex hull.'])
  end
  
end

function vDataSet = createWorkingDataSet(aImarisApplication)
  vOrigDataSet = aImarisApplication.GetDataSet();
  vDataSet = aImarisApplication.GetFactory().CreateDataSet;
  vSizeX = vOrigDataSet.GetSizeX();
  vSizeY = vOrigDataSet.GetSizeY();
  vSizeZ = vOrigDataSet.GetSizeZ();
  vSizeC = 1;
  vSizeT = 1;
  vType = Imaris.tType.eTypeUInt8;
  vExtentMin = [vOrigDataSet.GetExtendMinX(), vOrigDataSet.GetExtendMinY(), vOrigDataSet.GetExtendMinZ()];
  vExtentMax = [vOrigDataSet.GetExtendMaxX(), vOrigDataSet.GetExtendMaxY(), vOrigDataSet.GetExtendMaxZ()];
  vDataSet.Create(vType, vSizeX, vSizeY, vSizeZ, vSizeC, vSizeT);
  vDataSet.SetExtendMinX(vExtentMin(1));
  vDataSet.SetExtendMinY(vExtentMin(2));
  vDataSet.SetExtendMinZ(vExtentMin(3));
  vDataSet.SetExtendMaxX(vExtentMax(1));
  vDataSet.SetExtendMaxY(vExtentMax(2));
  vDataSet.SetExtendMaxZ(vExtentMax(3));
end

function [X, Y, Z] = getVoxelCenterCoordinates(aDataSet)
  vSize = [aDataSet.GetSizeX(), aDataSet.GetSizeY(), aDataSet.GetSizeZ()];
  vExtentMin = [aDataSet.GetExtendMinX(), aDataSet.GetExtendMinY(), aDataSet.GetExtendMinZ()];
  vExtentMax = [aDataSet.GetExtendMaxX(), aDataSet.GetExtendMaxY(), aDataSet.GetExtendMaxZ()];
  vVoxelSize = (vExtentMax - vExtentMin) ./ vSize;
  vX = vExtentMin(1) + (1:vSize(1)) * vVoxelSize(1) - 0.5 * vVoxelSize(1);
  vY = vExtentMin(2) + (1:vSize(2)) * vVoxelSize(2) - 0.5 * vVoxelSize(2);
  vZ = vExtentMin(3) + (1:vSize(3)) * vVoxelSize(3) - 0.5 * vVoxelSize(3);
  [X,Y,Z] = ndgrid(vX, vY, vZ);
end

function vSurface = DetectSurfaceWithoutSmoothing(aImarisApplication, aDataSet)
  vSize = [aDataSet.GetSizeX(), aDataSet.GetSizeY(), aDataSet.GetSizeZ()];
  vExtentMin = [aDataSet.GetExtendMinX(), aDataSet.GetExtendMinY(), aDataSet.GetExtendMinZ()];
  vExtentMax = [aDataSet.GetExtendMaxX(), aDataSet.GetExtendMaxY(), aDataSet.GetExtendMaxZ()];
  vVoxelSize = (vExtentMax - vExtentMin) ./ vSize;

  vROIs = [];
  vChannelIndex = 0;
  %in FilamentConvexHull vSmoothFilterWidth is half the size
  vSmoothFilterWidth = sqrt(dot(vVoxelSize, vVoxelSize)) * 0.25;
  vLocalContrastFilterWidth = 0;
  vIntensityThresholdAutomatic = false;
  vIntensityThresholdManual = 0.1;
  vSurfaceFiltersString = '';

  vImageProcessing = aImarisApplication.GetImageProcessing();
  vSurface = vImageProcessing.DetectSurfaces(aDataSet, ...
    vROIs, vChannelIndex, vSmoothFilterWidth, ...
    vLocalContrastFilterWidth, vIntensityThresholdAutomatic, ...
    vIntensityThresholdManual, vSurfaceFiltersString);
end
    