%
%
%  Spots And Surfaces Distance 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="Spots And Surfaces Distance" icon="Matlab" tooltip="Compute distance between spots and surfaces.">
%          <Command>Matlab::BPSpotsAndSurfacesDistance(%i)</Command>
%        </Item>
%      </Menu>
%     <SurpassToolbar>
%        <Item name="Spots And Surfaces Distance" icon="Matlab" tooltip="Compute distance between spots and surfaces.">
%          <Command>Matlab::BPSpotsAndSurfacesDistance(%i)</Command>
%        </Item>
%     </SurpassToolbar>
%      <SurpassTab>
%        <SurpassComponent name="bpComponentGroup">
%          <Item name="Compute Distance between Spots And Surfaces" icon="Matlab" tooltip="Compute distance between spots and surfaces.">
%            <Command>Matlab::BPSpotsAndSurfacesDistance(%i)</Command>
%          </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
% 
%
%  Description:
%   
%   Compute distance between all the pairs of visible and valid objects 
%       that are in the same group. 
%   Valid objects are surfaces and spots.
%   The center of a spots object is equal to the mean of the coordinates of
%       all the single spots that belong to the object. 
%   The center of a surface object is an approximation of the center of
%       mass of the object represented.
%
%

function BPSpotsAndSurfacesDistance(aImarisApplicationID)

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

vFactory = vImarisApplication.mFactory;

vSelection = vImarisApplication.mSurpassSelection;
if vFactory.IsDataContainer(vSelection)
    vParentsList{1} = vSelection;
elseif isequal(vSelection,[])
    vParentsList{1} = vImarisApplication.mSurpassScene;
    if isequal(vParentsList{1},[])
        msgbox 'Please create a surpass scene!';
        return
    end
else
    vParentsList{1} = vSelection.GetParent;
end    

vParentsNumber = 1;

vSpotsList = {};
vSpotsNumber = 0;

vSurfacesList = {};
vSurfacesNumber = 0;

% get all the spots and surfaces (including sub directories)
while vParentsNumber
    vParent = vParentsList{vParentsNumber};
    for vChild = 1:vParent.GetNumberOfChildren
        vItem = vParent.GetChild(vChild-1);
        if ~vItem.mVisible
            continue
        end
        if vFactory.IsSpots(vItem)
            vSpotsNumber = vSpotsNumber + 1;
            vSpotsList{vSpotsNumber} = vItem;
        elseif vFactory.IsSurface(vItem)
            vSurfacesNumber = vSurfacesNumber + 1;
            vSurfacesList{vSurfacesNumber} = vItem;
        elseif ~vFactory.IsTrack(vItem) && vFactory.IsDataContainer(vItem)
            vParentsNumber = vParentsNumber + 1;
            vParentsList = [{vItem},vParentsList];
        end
    end
    vParentsNumber = vParentsNumber - 1;
    vParentsList = vParentsList(1:vParentsNumber);
end


% compute the center of the objects
vObjectsNumber = vSpotsNumber+vSurfacesNumber;
vMaxObjects = 17;
if vObjectsNumber > vMaxObjects
    msgbox(sprintf(['WARNING: %i objects found. Unable to display the ', ...
        'distance between more than %i objects.'],vObjectsNumber,vMaxObjects));
    vSpotsNumber = round(vSpotsNumber/vObjectsNumber*vMaxObjects);
    vSurfacesNumber = vMaxObjects - vSpotsNumber;
    vObjectsNumber = vMaxObjects;
elseif vObjectsNumber < 2
    msgbox('Please create at least two valid objects!');
    return
end

vObjectNames = cell(1,vObjectsNumber);

vSpotsCenter = zeros(3,vSpotsNumber);
for vSpots = 1:vSpotsNumber
    vXYZ = vSpotsList{vSpots}.GetPositionsXYZ;
    vSpotsCenter(:,vSpots) = mean(vXYZ,1)';
    
    vObjectNames{vSpots} = vSpotsList{vSpots}.mName;
end

% define the maximal size of the mask, for cubic surfaces the mask has size
%  equal to vMaxEdge^3
vMaxEdge = 100;

vSurfacesCenter = zeros(3,vSurfacesNumber);
for vSurface = 1:vSurfacesNumber
%      vXYZ = vSurfacesList{vSurface}.GetVertices;
%      vSurfacesCenter(:,vSurface) = mean(vXYZ)'; 
     
    % get the center of the mask, not the mean of the vertices
    vSurfObject = vSurfacesList{vSurface};
    vXYZ = vSurfObject.GetVertices;
    vMin = min(vXYZ);
    vMax = max(vXYZ);
    vExtends = vMax - vMin;
    vSize = ceil(vExtends/max([vExtends,0.1])*vMaxEdge);
    vMask = vSurfObject.GetMask(vMin(1),vMin(2),vMin(3), ...
        vMax(1),vMax(2),vMax(3),vSize(1),vSize(2),vSize(3));
    
    vImage = zeros(vSize);
    for vSlice = 1:vSize(3)
        vImage(:,:,vSlice) = vMask.GetDataSlice(vSlice-1,0,0);
    end
    
    % count how many cells there are at x = 1,2,3,...
    vXY = sum(vImage,3);
    vX = sum(vXY,2);
    vY = sum(vXY,1);
    vZ = sum(sum(vImage,1),2);

    % computes the weighted mean
    vSurfacesCenter(1,vSurface) = linspace(vMin(1),vMax(1),vSize(1))*vX./sum(vX);
    vSurfacesCenter(2,vSurface) = linspace(vMin(2),vMax(2),vSize(2))*vY'./sum(vY);
    vSurfacesCenter(3,vSurface) = linspace(vMin(3),vMax(3),vSize(3))*shiftdim(vZ,2)./sum(vZ);

    vObjectNames{vSpotsNumber+vSurface} = vSurfacesList{vSurface}.mName;
end


% compute the pairwise distance
vDistanceMatrix = zeros(vObjectsNumber);

vOnes3xNSpots = ones(3,vSpotsNumber);
vOnes3xNSurfaces = ones(3,vSurfacesNumber);

for vSpots = 1:vSpotsNumber
    vXYZ = diag(vSpotsCenter(:,vSpots));
    vValid = vSpots+1:vSpotsNumber;
    vDistanceMatrix(vSpots,vValid) = sqrt(sum( ...
        (vSpotsCenter(:,vValid)-vXYZ*vOnes3xNSpots(:,vValid)).^2 ...
        ));
    
    vDistanceMatrix(vSpots,vSpotsNumber+(1:vSurfacesNumber)) = sqrt(sum( ...
        (vSurfacesCenter-vXYZ*vOnes3xNSurfaces).^2 ...
        ));
end

for vSurface = 1:vSurfacesNumber
    vXYZ = diag(vSurfacesCenter(:,vSurface));
    vValid = vSurface+1:vSurfacesNumber;
    vDistanceMatrix(vSpotsNumber+vSurface,vSpotsNumber+vValid) = sqrt(sum( ...
        (vSurfacesCenter(:,vValid)-vXYZ*vOnes3xNSurfaces(:,vValid)).^2 ...
        ));
end

% display the results
clf
vObjectsNumber = vObjectsNumber - 1;
vMaxNumber = 8;

vDeltaX = 125;
vDeltaY = 20;
vBorder = 20;

% color spots with red, surfaces with blue, min and max with green
vColors = [0.8,0.8,1;1,0.8,1;1,1,1;0.8,1,0.8];

% prevent to get dummy values (zeros) as min or max
vLower = max(max(vDistanceMatrix))*tril(ones(size(vDistanceMatrix)));

[vMax,vMaxIndexY] = max(vDistanceMatrix-vLower);
[vMax,vMaxIndexX] = max(vMax);
vMaxIndexY = vMaxIndexY(vMaxIndexX);
[vMin,vMinIndexY] = min(vDistanceMatrix+vLower);
[vMin,vMinIndexX] = min(vMin);
vMinIndexY = vMinIndexY(vMinIndexX);

vMinX = vBorder+(2*vMaxNumber-vObjectsNumber)*vDeltaX;%vDeltaX*min(vObjectsNumber-vMaxNumber,8)+vBorder;
vMinY = vDeltaY*vObjectsNumber+vBorder;


if vObjectsNumber > vMaxNumber

    for vObject = 1:vObjectsNumber
        uicontrol('Style', 'text', 'String', [vObjectNames{vObject},' '], ...
            'HorizontalAlignment', 'left', ...
            'BackgroundColor', vColors(int8(vObject<=vSpotsNumber)+1,:), ...
            'FontName', 'Times New Roman', 'FontSize', 10, ...
            'Position', [vMinX+(vObjectsNumber-vMaxNumber)*vDeltaX,vMinY-vDeltaY*vObject+1,vDeltaX,vDeltaY-2]);
        if vObject > vMaxNumber
            uicontrol('Style', 'text', 'String', [vObjectNames{vObject+1},' '], ...
                'HorizontalAlignment', 'right', ...
                'BackgroundColor', vColors(int8(vObject+1<=vSpotsNumber)+1,:), ...
                'FontName', 'Times New Roman', 'FontSize', 10, ...
                'Position', [vDeltaX*(vObject-1-vMaxNumber)+vMinX+1,vMinY,vDeltaX-2,vDeltaY]);
        end
        if vObject <= vMaxNumber
            for vObject2 = vMaxNumber+1:vObjectsNumber
                uicontrol('Style', 'edit', 'String', num2str(vDistanceMatrix(vObject,vObject2+1)), ...
                    'HorizontalAlignment', 'right', ...
                    'BackgroundColor', vColors(3+int8(vObject==vMaxIndexY)*int8(vObject2+1==vMaxIndexX)+ ...
                    int8(vObject==vMinIndexY)*int8(vObject2+1==vMinIndexX),:), ...
                    'FontName', 'Times New Roman', 'FontSize', 10, ...
                    'Position', [vDeltaX*(vObject2-1-vMaxNumber)+vMinX, ...
                    vMinY-vDeltaY*vObject,vDeltaX,vDeltaY]);
            end
        else
            for vObject2 = vObject:vObjectsNumber
                uicontrol('Style', 'edit', 'String', num2str(vDistanceMatrix(vObject,vObject2+1)), ...
                    'HorizontalAlignment', 'right', ...
                    'BackgroundColor', vColors(3+int8(vObject==vMaxIndexY)*int8(vObject2+1==vMaxIndexX)+ ...
                    int8(vObject==vMinIndexY)*int8(vObject2+1==vMinIndexX),:), ...
                    'FontName', 'Times New Roman', 'FontSize', 10, ...
                    'Position', [vDeltaX*(vObject2-1-vMaxNumber)+vMinX, ...
                    vMinY-vDeltaY*(vObject),vDeltaX,vDeltaY]);
            end
        end
    end
    
    vMinY = vMinY + vDeltaY*10;
    vObjectsNumber = vMaxNumber;
end

vMinX = vDeltaX*min(vObjectsNumber,vMaxNumber)+vBorder;

for vObject = 1:vObjectsNumber
    uicontrol('Style', 'text', 'String', [vObjectNames{vObject},' '], ...
        'HorizontalAlignment', 'left', ...
        'BackgroundColor', vColors(int8(vObject<=vSpotsNumber)+1,:), ...
        'FontName', 'Times New Roman', 'FontSize', 10, ...
        'Position', [vMinX,vMinY-vDeltaY*vObject+1,vDeltaX,vDeltaY-2]);
    uicontrol('Style', 'text', 'String', [vObjectNames{vObject+1},' '], ...
        'HorizontalAlignment', 'right', ...
        'BackgroundColor', vColors(int8(vObject+1<=vSpotsNumber)+1,:), ...
        'FontName', 'Times New Roman', 'FontSize', 10, ...
        'Position', [vDeltaX*(vObject-1)+vBorder+1,vMinY,vDeltaX-2,vDeltaY]);
    for vObject2 = vObject:vObjectsNumber
        uicontrol('Style', 'edit', 'String', num2str(vDistanceMatrix(vObject,vObject2+1)), ...
            'HorizontalAlignment', 'right', ...
            'BackgroundColor', vColors(3+int8(vObject==vMaxIndexY)*int8(vObject2+1==vMaxIndexX)+ ...
            int8(vObject==vMinIndexY)*int8(vObject2+1==vMinIndexX),:), ...
            'FontName', 'Times New Roman', 'FontSize', 10, ...
            'Position', [vDeltaX*(vObject2-1)+vBorder,vMinY-vDeltaY*vObject,vDeltaX,vDeltaY]);
    end
end

vWindowPosition = get(gcf, 'Position');
vWindowPosition(3:4) = [vMinX+vDeltaX+vBorder,vMinY+vDeltaY+vBorder];
set(gcf, 'Position', vWindowPosition);