%
%
%  Find Spots Close To Surface 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="Spots Functions">
%        <Item name="Find Spots Close To Surface" icon="Matlab" tooltip="Find spots close to surface.">
%          <Command>Matlab::BPSpotsCloseToSurface(%i)</Command>
%        </Item>
%       </Submenu>
%      </Menu>
%      <SurpassTab>
%        <SurpassComponent name="bpSpots">
%          <Item name="Find Spots Close To Surface" icon="Matlab" tooltip="Find spots close to surface.">
%            <Command>Matlab::BPSpotsCloseToSurface(%i)</Command>
%          </Item>
%        </SurpassComponent>
%        <SurpassComponent name="bpSurfaceObject">
%          <Item name="Find Spots Close To Surface" icon="Matlab" tooltip="Find spots close to surface.">
%            <Command>Matlab::BPSpotsCloseToSurface(%i)</Command>
%          </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
% 
%
%  Description:
%   
%   Find the spots that are closer to a surface than a threshold, and 
%       create two new spots objects: One collects the spots close to the
%       surface, the other the farer. 
%   If multiple spots or surfaces are selected, each pair [spots, surface] 
%       is analyzed: The number of spots objects created is 2 times 
%       'number of spots' times 'number of surfaces'. 
%
%

function BPSpotsCloseToSurface(aImarisApplicationID)

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

% the user has to create a scene with some spots and surface
vSurpassScene = vImarisApplication.mSurpassScene;
if isequal(vSurpassScene, [])
    msgbox('Please create some Spots and Surface in the Surpass scene!');
    return;
end

% get the spots and the surface object
vSpots = vImarisApplication.mFactory.ToSpots(vImarisApplication.mSurpassSelection);
vSurface = vImarisApplication.mFactory.ToSurface(vImarisApplication.mSurpassSelection);

vSpotsSelected = ~isequal(vSpots,[]);
vSurfaceSelected = ~isequal(vSurface,[]);
if vSpotsSelected
    vParent = vSpots.GetParent;
elseif vSurfaceSelected
    vParent = vSurface.GetParent;
else
    vParent = vSurpassScene;
end

% get the spots and surfaces
vSpotsSelection = 1;
vSurfaceSelection = 1;
vNumberOfSpots = 0;
vNumberOfSurfaces = 0;
vSpotsList = [];
vSurfacesList = [];
vSpotsName = {};
vSurfacesName = {};
for vIndex = 1:vParent.GetNumberOfChildren
    vItem = vParent.GetChild(vIndex-1);
    if vImarisApplication.mFactory.IsSpots(vItem)
        vNumberOfSpots = vNumberOfSpots + 1;
        vSpotsList(vNumberOfSpots) = vIndex;
        vSpotsName{vNumberOfSpots} = vItem.mName;
        
        if vSpotsSelected && strcmp(vItem.mName,vSpots.mName)
            vSpotsSelection = vNumberOfSpots;
        end
    elseif vImarisApplication.mFactory.IsSurface(vItem)
        vNumberOfSurfaces = vNumberOfSurfaces + 1;
        vSurfacesList(vNumberOfSurfaces) = vIndex;
        vSurfacesName{vNumberOfSurfaces} = vItem.mName;
        
        if vSurfaceSelected && strcmp(vItem.mName,vSurface.mName)
            vSurfaceSelection = vNumberOfSurfaces;
        end
    end
end

if min(vNumberOfSpots,vNumberOfSurfaces) == 0
    msgbox('Please create some spots AND a surface object!');
    return;
end

if vNumberOfSpots>1
    [vSpotsSelection,vOk] = listdlg('ListString',vSpotsName, ...
        'InitialValue', vSpotsSelection, 'SelectionMode','multiple', ...
        'ListSize',[300 300], 'Name','Find Spots Close To Surface', ...
        'PromptString',{'Please select the spots:'});
    if vOk<1, return; end
end
if vNumberOfSurfaces>1
    [vSurfaceSelection,vOk] = listdlg('ListString',vSurfacesName, ...
        'InitialValue', vSurfaceSelection, 'SelectionMode','multiple', ...
        'ListSize',[300 300], 'Name','Find Spots Close To Surface', ...
        'PromptString',{'Please select the surface:'});
    if vOk<1, return; end
end

vAnswer = inputdlg({'Please enter the threshold:'}, ...
    'Find Spots Close To Surface',1,{'5'});
if isempty(vAnswer), return; end;
vThreshold = abs(str2double(vAnswer{1}));

vProgressDisplay = waitbar(0,'Finding Spots Close To Surface');

% compute the distances and create new spots objects
vNumberOfSurfacesSelected = length(vSurfaceSelection);
vNumberOfSpotsSelected = length(vSpotsSelection);

for vSurfaceIndex = 1:vNumberOfSurfacesSelected
    vSurface = vParent.GetChild(vSurfacesList( ...
        vSurfaceSelection(vSurfaceIndex)) - 1);

    vSurfaceVertices = vSurface.GetVertices;
    vNumberOfVertices = size(vSurfaceVertices,1);
    
    % limit the memory usage to 3*1000*vNumberOfSpots
    if vNumberOfVertices>1000
        vIndices = round(((1:1000)-0.5)*vNumberOfVertices/1000);
        vSurfaceVertices = vSurfaceVertices(vIndices,:);
        vNumberOfVertices = 1000;
    end
    
    vSparseVertices = 1:vNumberOfVertices;
    
    vSurfX = sparse(vSparseVertices,vSparseVertices,double(vSurfaceVertices(:,1)));
    vSurfY = sparse(vSparseVertices,vSparseVertices,double(vSurfaceVertices(:,2)));
    vSurfZ = sparse(vSparseVertices,vSparseVertices,double(vSurfaceVertices(:,3)));

    for vSpotsIndex = 1:vNumberOfSpotsSelected
        vSpots = vParent.GetChild(vSpotsList( ...
            vSpotsSelection(vSpotsIndex)) - 1);

        [vSpotsPosition,vSpotsTime,vSpotsRadius] = vSpots.Get;
        vNumberOfSpots = size(vSpotsPosition,1);

        
        % build the distance matrix
        vOnes = ones(vNumberOfVertices,vNumberOfSpots);
        vDistX = vOnes*diag(vSpotsPosition(:,1)) - vSurfX*vOnes;
        vDistY = vOnes*diag(vSpotsPosition(:,2)) - vSurfY*vOnes;
        vDistZ = vOnes*diag(vSpotsPosition(:,3)) - vSurfZ*vOnes;
        
        vDistances = sqrt( min( vDistX.^2+vDistY.^2+vDistZ.^2 ) );
        
        vSpots.mVisible = false;
        
        vSpotsClose = find(vDistances<=vThreshold);
        vNewSpotsClose = vImarisApplication.mFactory.CreateSpots;
        vNewSpotsClose.Set(vSpotsPosition(vSpotsClose,:), ...
            vSpotsTime(vSpotsClose), vSpotsRadius(vSpotsClose));
        vNewSpotsClose.SetColor(1,0,1,0);
        vNewSpotsClose.mName = sprintf('%s close to %s [%.2f]', ...
            vSpots.mName,vSurface.mName,vThreshold);
        vParent.AddChild(vNewSpotsClose);
            
        vSpotsFar = find(vDistances>vThreshold);
        vNewSpotsFar = vImarisApplication.mFactory.CreateSpots;
        vNewSpotsFar.Set(vSpotsPosition(vSpotsFar,:), ...
            vSpotsTime(vSpotsFar), vSpotsRadius(vSpotsFar));
        vNewSpotsFar.SetColor(0,1,1,0);
        vNewSpotsFar.mName = sprintf('%s far from %s [%.2f]', ...
            vSpots.mName,vSurface.mName,vThreshold);
        vParent.AddChild(vNewSpotsFar);
                
        waitbar((vNumberOfSpotsSelected*(vSurfaceIndex-1)+vSpotsIndex) / ...
            (vNumberOfSpotsSelected*vNumberOfSurfacesSelected));
    end
end

close(vProgressDisplay);