%
%
%  Find Spots On Ends Of Lines 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>
%       <Submenu name="Spots Functions">
%        <Item name="Find Spots On Ends Of Lines" icon="Matlab" tooltip="Find spots on ends of lines.">
%          <Command>Matlab::BPSpotsOnEndsOfLines(%i)</Command>
%        </Item>
%       </Submenu>
%      </Menu>
%      <SurpassTab>
%        <SurpassComponent name="bpSpots">
%        <Item name="Find Spots On Ends Of Lines" icon="Matlab" tooltip="Find spots on ends of lines.">
%          <Command>Matlab::BPSpotsOnEndsOfLines(%i)</Command>
%        </Item>
%        </SurpassComponent>
%      </SurpassTab>
%    </CustomTools>
% 
% 
%  Description:
%
%   Segment the dataset and subdivide the spots in groups of spots that 
%       lie inside (or around) the same line. 
%   For each spots group discard all but the two most distant spots,
%       leaving more or less those on the line ends.
%   The image is read from channel 1 and the segmentation is returned in 
%       channel 2; All other channels are removed.
%
%   Example: 
%       - Open Tsia10tp.ims
%       - Detect spots (minimum diameter = 4)
%       - Run BPSpotsOnEndsOfLines
%
%

function BPSpotsOnEndsOfLines(aImarisApplicationID)

% connect to Imaris Com interface
if ~isa(aImarisApplicationID, 'COM.Imaris_Application')
  vImarisServer = actxserver('ImarisServer.Server');
  vImarisApplication = vImarisServer.GetObject(aImarisApplicationID);
else
  vImarisApplication = aImarisApplicationID;
end
vImarisApplication.DataSetPushUndo('Find Spots On Ends Of Lines');

vDataSet = vImarisApplication.mDataSet.Clone;

vSourceChannel = 0;
vDataSet.mSizeC = 2;

vExtendsMin = [vDataSet.mExtendMinX,vDataSet.mExtendMinY,vDataSet.mExtendMinZ];
vExtendsMax = [vDataSet.mExtendMaxX,vDataSet.mExtendMaxY,vDataSet.mExtendMaxZ];
vDataSize = [vDataSet.mSizeX,vDataSet.mSizeY,vDataSet.mSizeZ];
vRange = max(0.001,vExtendsMax-vExtendsMin);
vVoxelSize = (vExtendsMax-vExtendsMin)./max(vDataSize-1,1);

% the user has to create a scene with some spots
if isequal(vImarisApplication.mSurpassSelection, [])
    msgbox('Please select a Spots component in the Surpass scene!');
    return;
end

% search the spots
vSpots = vImarisApplication.mFactory.ToSpots(vImarisApplication.mSurpassSelection);
if ~vImarisApplication.mFactory.IsSpots(vSpots)
    msgbox('Please select a Spots component in the Surpass scene!');
    return;
end

% get the spot coordinates
vSpotRadii = vSpots.GetRadii;
vSpotTimeIndices = vSpots.GetIndicesT;
vSpotPositionsXYZ = vSpots.GetPositionsXYZ;

vSpotIndicesXYZ = vSpotPositionsXYZ;
vSpotIndicesXYZ(:,1) = round((vSpotPositionsXYZ(:,1)-vExtendsMin(1))./ ...
    vRange(1).*(vDataSize(1)-1))+1;
vSpotIndicesXYZ(:,2) = round((vSpotPositionsXYZ(:,2)-vExtendsMin(2))./ ...
    vRange(2).*(vDataSize(2)-1))+1;
vSpotIndicesXYZ(:,3) = round((vSpotPositionsXYZ(:,3)-vExtendsMin(3))./ ...
    vRange(3).*(vDataSize(3)-1))+1;

vEndPointPositionsXYZ = [];
vEndPointTimeIndices = [];
vEndPointRadii = [];

vAnswer = inputdlg({'Line width (in pixels)','N sigma',...
    'Contrast large','Contrast small','Background factor'},...
    'Please enter the segmentation parameters:',...
    1,{'2.5','5','3.5','2.5','3'});
if isempty(vAnswer), return; end;

vSegmentParam = zeros(5,1);
for vIter = 1:5
    vSegmentParam(vIter) = sscanf(char(vAnswer(vIter)),'%f');
end

vProgressDisplay = waitbar(0,'Selecting the spots on the ends of the lines');

for vTimeIndex = 1:vDataSet.mSizeT
    vImageThisTimePoint = vDataSet.GetDataVolumeAs1DArray(vSourceChannel, vTimeIndex-1);
    vImageThisTimePoint = reshape(vImageThisTimePoint, vDataSize);

    vSpotIndicesThisTimePointXYZ = vSpotIndicesXYZ( (vSpotTimeIndices+1) == vTimeIndex, :);
    
    vSpotIndicesThisTimePointXYZ(vSpotIndicesThisTimePointXYZ(:,1)>size(vImageThisTimePoint,1),1) = size(vImageThisTimePoint,1);
    vSpotIndicesThisTimePointXYZ(vSpotIndicesThisTimePointXYZ(:,2)>size(vImageThisTimePoint,2),2) = size(vImageThisTimePoint,2);
    vSpotIndicesThisTimePointXYZ(vSpotIndicesThisTimePointXYZ(:,3)>size(vImageThisTimePoint,3),3) = size(vImageThisTimePoint,3);  
    
    vSpotIndicesThisTimePoint = sub2ind(size(vImageThisTimePoint), vSpotIndicesThisTimePointXYZ(:,1),vSpotIndicesThisTimePointXYZ(:,2),vSpotIndicesThisTimePointXYZ(:,3));
    [vSegmentation, vEndPoints] = TsiavaliarisOneTimePoint(vImageThisTimePoint, vSpotIndicesThisTimePoint, vSegmentParam);
    vEndPoints(:,3) = 1;
    vNEndPoints = size(vEndPoints,1);
    vEndPointPositionsXYZ = [vEndPointPositionsXYZ; vEndPoints];
    vEndPointTimeIndices = [vEndPointTimeIndices; vTimeIndex+zeros(vNEndPoints,1)];
    vEndPointRadii = [vEndPointRadii; linspace(vSpotRadii(1),vSpotRadii(1),vNEndPoints)'];
  
    if strcmp(vDataSet.mType,'eTypeUInt8')
        vDataSet.SetDataVolume(uint8((vSegmentation > 0)*100),1,vTimeIndex-1);
    elseif strcmp(vDataSet.mType,'eTypeUInt16')
        vDataSet.SetDataVolume(uint16((vSegmentation > 0)*100),1,vTimeIndex-1);
    elseif strcmp(vDataSet.mType,'eTypeFloat')
        vDataSet.SetDataVolume(single((vSegmentation > 0)*100),1,vTimeIndex-1);
    end
    
    waitbar(vTimeIndex/vDataSet.mSizeT)
end
vImarisApplication.mDataSet = vDataSet;

% hide original spots
vSpots.mVisible = false;

% add new spots (close)
vSpotsOnEnds = vImarisApplication.mFactory.CreateSpots;

vEndPointXYZ(:,1) = (vEndPointPositionsXYZ(:,1)-1)*vVoxelSize(1) + vExtendsMin(1);
vEndPointXYZ(:,2) = (vEndPointPositionsXYZ(:,2)-1)*vVoxelSize(2) + vExtendsMin(2);
vEndPointXYZ(:,3) = (vEndPointPositionsXYZ(:,3)-1)*vVoxelSize(3) + vExtendsMin(3);

vSpotsOnEnds.Set(vEndPointXYZ, vEndPointTimeIndices-1, vEndPointRadii);
vSpotsOnEnds.mName = sprintf('%s on ends', vSpots.mName);
vSpotsOnEnds.SetColor(1.0, 0.0, 0.0, 0.0);
vImarisApplication.mSurpassScene.AddChild(vSpotsOnEnds);

close(vProgressDisplay);


%----------------------------------------------------------------%

function [aSegmentation, aEndPositions] = TsiavaliarisOneTimePoint( ...
    aImage, aSpotPositions, aSegmentParam)

%aLineWidthInPixels = 2.5;
aLineWidthInPixels = aSegmentParam(1);

%% segmentation of channel using multiscale localcontrast
vSigmaLarge = 4*aLineWidthInPixels;
vSigmaSmall = aLineWidthInPixels;
%vNSigma = 5;
vNSigma = aSegmentParam(2);
%vContrastLarge = 3.5;
vContrastLarge = aSegmentParam(3);
%vContrastSmall = 2.5;
vContrastSmall = aSegmentParam(4);
%vBackgroundFactor = 3;
vBackgroundFactor = aSegmentParam(5);
aSegmentation = TraceHighContrastRegions(aImage, vSigmaLarge, .8*vSigmaSmall, vNSigma, vContrastLarge, vContrastSmall, vBackgroundFactor);

% connected component labeling 
aSegmentation = bwlabeln(aSegmentation>0);

% spot detection
% vSigmaSmall = .6*aLineWidthInPixels;
% vSigmaLarge = 1.7*vSigmaSmall;
% vSmooth = GaussFilter(single(aImage), vSigmaSmall);
% vSmooth2 = GaussFilter(single(aImage), vSigmaLarge);
% vLaplace = vSmooth-vSmooth2;
% vSpotPositions = LocalMax(vLaplace);
vSpotPositions = aSpotPositions';

% loop over regions/labels
vLabels = setdiff(unique(aSegmentation),0)';
vSpotPositionLabels = aSegmentation(vSpotPositions);
aEndPositions = zeros(numel(2*vLabels),2);
aEndIndices= zeros(numel(2*vLabels),1);
for vLabel = vLabels
  vSpotPositionsThisLabel = vSpotPositions(vSpotPositionLabels == vLabel);
  vSpotPositionsThisLabelXY = [];
  if numel(vSpotPositionsThisLabel) > 0
  [vSpotPositionsThisLabelXY(:,1), vSpotPositionsThisLabelXY(:,2)] = ind2sub(size(aImage), vSpotPositionsThisLabel);
  if numel(vSpotPositionsThisLabel) > 2
    % find the spots within this region that are farthest apart
    vOverallMaxDist = 0;
    for vIndex = 1:numel(vSpotPositionsThisLabel)
      vPosition = vSpotPositionsThisLabelXY(vIndex,:);
      vDist = sqrt((vSpotPositionsThisLabelXY(:,1) - vPosition(1)).^2 + (vSpotPositionsThisLabelXY(:,2) - vPosition(2)).^2); 
      [vMaxDist, vMaxIndex] = max(vDist);
      if (vMaxDist > vOverallMaxDist)
        aEndPositions(2*vLabel-1,1) = vSpotPositionsThisLabelXY(vIndex,1);
        aEndPositions(2*vLabel-1,2) = vSpotPositionsThisLabelXY(vIndex,2);
        aEndPositions(2*vLabel,1) = vSpotPositionsThisLabelXY(vMaxIndex,1);
        aEndPositions(2*vLabel,2) = vSpotPositionsThisLabelXY(vMaxIndex,2);
        aEndIndices(2*vLabel-1) = vSpotPositionsThisLabel(vIndex);
        aEndIndices(2*vLabel) = vSpotPositionsThisLabel(vMaxIndex);
        vOverallMaxDist = vMaxDist;
      end
    end
  else
    if numel(vSpotPositionsThisLabel) == 1
      aEndPositions(2*vLabel-1,1) = vSpotPositionsThisLabelXY(1,1);
      aEndPositions(2*vLabel-1,2) = vSpotPositionsThisLabelXY(1,2);
    elseif numel(vSpotPositionsThisLabel) == 2
      aEndPositions(2*vLabel-1,1) = vSpotPositionsThisLabelXY(1,1);
      aEndPositions(2*vLabel-1,2) = vSpotPositionsThisLabelXY(1,2);
      aEndPositions(2*vLabel,1) = vSpotPositionsThisLabelXY(2,1);
      aEndPositions(2*vLabel,2) = vSpotPositionsThisLabelXY(2,2);
    end
  end
  end
end


%%
function aHighContrastRegions = TraceHighContrastRegions( ...
    aImage, aSigmaLarge, aSigmaSmall, aNSigma, aContrastLarge, aContrastSmall, aBGFactor)

if nargin < 7
  aBGFactor = 2;
end

vSigmaVector = zeros(aNSigma, ndims(aImage));

for vDim = 1:ndims(aImage)
  vDimL = min(vDim, numel(aSigmaLarge));
  vDimS = min(vDim, numel(aSigmaSmall));
  vSigmaVector(:,vDim) = aSigmaLarge(vDimL) - (1:aNSigma)/aNSigma * (aSigmaLarge(vDimL)-aSigmaSmall(vDimS));
end
vContrastThreshold = aContrastLarge -  (1:aNSigma)/aNSigma * (aContrastLarge-aContrastSmall);
aHighContrastRegions = (LocalHistogramContrast(single(aImage), aSigmaLarge, aBGFactor*aSigmaLarge)>aContrastLarge);
vLabels = bwlabeln(aHighContrastRegions);
for vIndex = 1:aNSigma
  vSigma = vSigmaVector(vIndex,:);
  vLocalContrast = LocalHistogramContrast(single(aImage), vSigma, aBGFactor*vSigma);
  aHighContrastRegions = aHighContrastRegions & vLocalContrast>vContrastThreshold(vIndex);  
%  show(aHighContrastRegions,1);
end
vLabels(~aHighContrastRegions) = 0;
aHighContrastRegions = vLabels;


%%
function aContrastImage = LocalHistogramContrast(aImage, aSigmaFG, aSigmaBG)

vSmoothFG = GaussFilter(single(aImage), aSigmaFG);
vSmoothBG = GaussFilter(single(aImage), aSigmaBG);
vSmoothBG2 = GaussFilter(single(aImage).*single(aImage), aSigmaBG);
aContrastImage = (vSmoothFG - vSmoothBG)./((4*pi)^(-3/4)*det(diag(aSigmaFG))^(-1/4)*sqrt(vSmoothBG2 - vSmoothBG.^2));

%%
function aOutImage = GaussFilter(aInImage, aSigma)

for vDim = (numel(aSigma)+1):3
  aSigma(vDim) = aSigma(vDim-1);
end
global vImarisApplitionIP;
if ~isa(vImarisApplitionIP, 'COM.Imaris_Application')
  vImarisApplitionIP = actxserver('Imaris.Application');
end
vImarisApplitionIP.mViewer = 0;
vImarisApplitionIP.mDataSet.Create('eTypeFloat', size(aInImage,1), size(aInImage,2), size(aInImage,3), 1, 1);
BPSetImarisImage(vImarisApplitionIP, aInImage, 1,1);
vImarisApplitionIP.mDataSet.mExtendMaxX = vImarisApplitionIP.mDataSet.mExtendMinX + vImarisApplitionIP.mDataSet.mSizeX / aSigma(1);
vImarisApplitionIP.mDataSet.mExtendMaxY = vImarisApplitionIP.mDataSet.mExtendMinY + vImarisApplitionIP.mDataSet.mSizeY / aSigma(2);
vImarisApplitionIP.mDataSet.mExtendMaxZ = vImarisApplitionIP.mDataSet.mExtendMinZ + vImarisApplitionIP.mDataSet.mSizeZ / aSigma(3);
vImarisApplitionIP.mImageProcessing.GaussFilterChannel(vImarisApplitionIP.mDataSet, 0, 1);
aOutImage = BPGetImarisImage(vImarisApplitionIP, 1, 1);

%%
function [aImage, aVoxelSize] = BPGetImarisImage(aImarisApplication, aTimeIndex, aChannelIndex)

vImarisDataSet = aImarisApplication.mDataSet;

vExtendMin = [vImarisDataSet.mExtendMinX, vImarisDataSet.mExtendMinY, ...
  vImarisDataSet.mExtendMinZ];
vExtendMax = [vImarisDataSet.mExtendMaxX, vImarisDataSet.mExtendMaxY, ...
  vImarisDataSet.mExtendMaxZ];
vDataSize = [vImarisDataSet.mSizeX, vImarisDataSet.mSizeY, ...
  vImarisDataSet.mSizeZ, vImarisDataSet.mSizeC, vImarisDataSet.mSizeT];

aVoxelSize = (vExtendMax - vExtendMin) ./ vDataSize(1:3);

vMinT = 0;
vMinC = 0;

if nargin == 3 && aChannelIndex<=vDataSize(4) && aTimeIndex<=vDataSize(5)

  vMinT = aTimeIndex-1;
  vDataSize(5) = 1;
  vMinC = aChannelIndex-1;
  vDataSize(4) = 1;
  
elseif nargin == 2 && aTimeIndex<=vDataSize(5)

  vMinT = aTimeIndex-1;
  vDataSize(5) = 1;

end

vSize = [1, prod(vDataSize(1:3)), vDataSize(4:5)];

if strcmp(vImarisDataSet.mType,'eTypeUInt8')
  aImage = zeros(vSize, 'uint8');
elseif strcmp(vImarisDataSet.mType,'eTypeUInt16')
  aImage = zeros(vSize, 'uint16');
else
  aImage = zeros(vSize, 'single');
end

for vTimeIndex = vMinT + (1:vDataSize(5))
  for vChannelIndex = vMinC + (1:vDataSize(4))
    aImage(1, :, vTimeIndex, vChannelIndex) = vImarisDataSet.GetDataVolumeAs1DArray( ...
      vChannelIndex-1, vTimeIndex-1);
  end
end
aImage = reshape(aImage, vDataSize);


%%
function BPSetImarisImage(aImarisApplication, aImage, aTimeIndex, aChannelIndex)

vImarisDataSet = aImarisApplication.mDataSet.Clone;
if vImarisDataSet.mSizeX ~= size(aImage,1)  || vImarisDataSet.mSizeY ~= size(aImage,2) || vImarisDataSet.mSizeZ ~= size(aImage,3)
   vImarisDataSet.Resize(0,size(aImage,1), 0,size(aImage,2), 0, size(aImage,3), 0, vImarisDataSet.mSizeC, 0, vImarisDataSet.mSizeT);
end
if size(aImage,1) == vImarisDataSet.mSizeX && ...
        size(aImage,2) == vImarisDataSet.mSizeY && ...
        size(aImage,3) == vImarisDataSet.mSizeZ
   if strcmp(vImarisDataSet.mType,'eTypeUInt8')
        vMin = min(reshape(aImage,[1,numel(aImage)]));
        if vMin < 0
            aImage = aImage + vMin;
        end
        vMax = max(reshape(aImage,[1,numel(aImage)]));
        if vMax > (2^8-1)
            aImage = aImage*(2^8-1)/vMax;
        end
        vImarisDataSet.SetDataVolume(uint8(aImage), aChannelIndex-1, aTimeIndex-1);
   elseif strcmp(vImarisDataSet.mType,'eTypeUInt16')
        vMin = min(reshape(aImage,[1,numel(aImage)]));
        if vMin < 0
            aImage = aImage + vMin;
        end
        vMax = max(reshape(aImage,[1,numel(aImage)]));
        if vMax > (2^16-1)
            aImage = aImage*(2^16-1)/vMax;
        end
          vImarisDataSet.SetDataVolume(uint16(aImage), aChannelIndex-1, aTimeIndex-1);
   elseif strcmp(vImarisDataSet.mType,'eTypeFloat')
        vImarisDataSet.SetDataVolume(single(aImage), aChannelIndex-1, aTimeIndex-1);
   end
end
aImarisApplication.mDataSet = vImarisDataSet;