%
%
%  Datasets Stitcher Function for Imaris
%
%  Copyright Bitplane AG 2008
%
%
%  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="Stitch Datasets">
%          <Command>Matlab::BPStitchDatasets(%i)</Command>
%        </Item>
%      </Menu>
%      <SurpassToolbar>
%        <Item name="Stitch DataSets" icon="S" tooltip="Load and align a number of datasets">
%          <Command>Matlab::BPStitchDatasets(%i)</Command>
%        </Item>
%      </SurpassToolbar>
%    </CustomTools>
% 
%
%  Description:
%   
%   Load and stitch datasets. 
%   Voxelsize, channel colors and other parameters (dataset's name, 
%       description, date, ...) are read from the first file loaded. 
%   Refinement is the average shift computed from all channels, 
%       all time points.
%   Adjustment (of the intensities) is computed independently for 
%       each channel, each timepoint.
% 
% 

function BPStitchDatasets(aImarisApplicationID)

% callbacks need (generally) persistent or global variables
persistent mImarisApplication mUIControls mOptions

% button callbacks
if isa(aImarisApplicationID, 'numeric') && aImarisApplicationID < 0
  if aImarisApplicationID == -1
    clear mImarisApplication mUIControls mOptions
  else
    [mUIControls, mOptions] = Process(-aImarisApplicationID, ...
      mImarisApplication, mUIControls, mOptions);
  end
  return
end % button callbacks

% get the application object
if isa(aImarisApplicationID, 'COM.Imaris_Application')
  % called from workspace
  mImarisApplication = aImarisApplicationID;
else
  % connect to Imaris Com interface
  vImarisServer = actxserver('ImarisServer.Server');
  if ischar(aImarisApplicationID)
    aImarisApplicationID = round(str2double(aImarisApplicationID));
  end
  mImarisApplication = vImarisServer.GetObject(aImarisApplicationID);
end

% builds the GUI
mUIControls = BuildGUI;

% default values for persistent variables
mOptions.mFilePath = 'C:\Program Files\Bitplane\';
mOptions.mFileName = '';
mOptions.mLayoutFromFileNames = 1;
mOptions.mFlipY = 0;
mOptions.mEditOverlapX = 0;
mOptions.mEditOverlapY = 0;
mOptions.mEditOverlapZ = 0;
mOptions.mEditShiftX = 0;
mOptions.mEditShiftY = 0;
mOptions.mEditShiftZ = 0;
mOptions.mAdjust = 0; 
mOptions.mFiles = {};
mOptions.mIndexSize = [0, 0, 0];
mOptions.mGridX = 1;
mOptions.mGridY = 1;
mOptions.mCoords = [];
mOptions.mImarisSlave = [];
% mOptions.mImarisSlave = mImarisApplication;

% create "slave" Imaris instance to load single datasets
mOptions.mImarisSlave = actxserver('Imaris.Application');
mOptions.mImarisSlave.mView = 'eViewerSlice';

% throw resize callback to refresh window
BPStitchDatasets(-3);


%-------------------------------------------------------------------------%
%% GUI Callbacks
%-------------------------------------------------------------------------%
function [aUIControls, aOptions] = Process(aEventNumber, ...
  aImarisApplication, aUIControls, aOptions)

% which action has been requested? 
switch aEventNumber

%-------------------------------------------------------------------------%
% Stitch command
%-------------------------------------------------------------------------%
  case 2
    
    % open waitbar
    vBar = waitbar(0, 'ImarisXT - Estimating blocks position');

    % one should not use the window while work is in progress
    vFields = fieldnames(aUIControls);
    vPanelText = 'mPanel';
    vPanelTextLenght = numel(vPanelText);
    for vIndex = 1:numel(vFields)
      vFieldName = vFields{vIndex};
      if numel(vFieldName) < vPanelTextLenght || ...
          ~isequal(vFieldName(1:vPanelTextLenght), vPanelText)
        set(aUIControls.(vFieldName), 'Enable', 'off');
      end
    end

    % copy option values (possible changes would not be persistent) 
    vFileList = aOptions.mFiles;
    for vIndex = 1:numel(vFileList)
      vPath = fileparts(vFileList{vIndex});
      if isempty(vPath)
        vFileList{vIndex} = fullfile(aOptions.mFilePath, vFileList{vIndex});
      end
    end
    vRefineStitch = [aOptions.mEditShiftX, aOptions.mEditShiftY, ...
      aOptions.mEditShiftZ];
    vAdjustIntensities = aOptions.mAdjust;
    vDimensionX = aOptions.mIndexSize(1);
    vDimensionY = aOptions.mIndexSize(2);
    vOverlap = [aOptions.mEditOverlapX, aOptions.mEditOverlapY, ...
      aOptions.mEditOverlapZ];
    vCoords = aOptions.mCoords;
        
    % ExtMode: 1 - Coords from Imaris, 2 - from text (already read) 
    %          3 - Append volumes,     4 - from filenames (already read, append)
    vExtendMode = 1 + 2*aOptions.mLayoutFromFileNames;
    if ~isempty(vCoords)
      vExtendMode = 2;
    end
            
    % need a backup?
    if ~isempty(aImarisApplication.mDataSet)
      % we need a try statement since in fact mDataSet is never empty, 
      %   but it might be invalid :/
      try
        aImarisApplication.DataSetPushUndo('Stitch DataSets');
      catch
        warning('IMARIS:pushUndo', 'Could not push undo dataset')
      end
    end
    
    % create "slave" Imaris instance to load single datasets
%     vImarisSlave = actxserver('Imaris.Application');
%     vImarisSlave.mView = 'eViewerSlice';
    vImarisSlave = aOptions.mImarisSlave;
    
    vDataSet = vImarisSlave.mDataSet;
    
    % reserve some space for the result
    vListIndex = zeros(1, numel(vFileList));
    vExtendList = zeros(numel(vFileList), 6);
    vMaxSizeCT = [0, 0];
    
    vNumFiles = 0;
    for vIndex = 1:numel(vFileList)
      vName = vFileList{vIndex};
      
      waitbar(vIndex/numel(vFileList)/4, vBar)

      % it exists?
      if ~exist(vName, 'file')
        warning('MATLAB:fileNotFound', 'File not found: %s', vName)
        continue
      end
      
      % Imaris can read it?
      vSuccess = ImarisFileOpen(vImarisSlave, vName);
      if ~vSuccess 
        continue
      end
      
      % add file to the list of the valid files
      vNumFiles = vNumFiles + 1;
      vListIndex(vNumFiles) = vIndex;
      
      % get extends and size
      vCurrentExtends = [vDataSet.mExtendMinX, vDataSet.mExtendMinY, ...
        vDataSet.mExtendMinZ, vDataSet.mExtendMaxX, ...
        vDataSet.mExtendMaxY, vDataSet.mExtendMaxZ];
      vSize = [vDataSet.mSizeX, vDataSet.mSizeY, vDataSet.mSizeZ];
      vSizeCT = [vDataSet.mSizeC, vDataSet.mSizeT];
      vMaxSizeCT = max(vMaxSizeCT, vSizeCT);      
      
      % compute extends
      switch vExtendMode 
        case 1 % coords from Imaris
          
          % current extends already computed
          
        case 2 % user defined (text)
          
          vCurrentExtends = [vCoords(vIndex, :), vCoords(vIndex, :) + ...
            (vCurrentExtends(4:6) - vCurrentExtends(1:3))];
          
        case 3 % append
          
          if vNumFiles == 1
            vIndexXYZ = [0, 0, 0]; % current index
            vXYZ = vCurrentExtends(1:3);
            vMax = vCurrentExtends(4:6); 
            vNextYZ = vCurrentExtends(5:6); 
          else
            vIndexXYZ(1) = vIndexXYZ(1) + 1;
            vXYZ(1) = vMax(1) - vOverlap(1) + vFirstPosition(1);
            vNextYZ = max(vMax(2:3), vNextYZ);
            if vIndexXYZ(1) >= vDimensionX % end of row
              vIndexXYZ(1) = 0;
              vIndexXYZ(2) = vIndexXYZ(2) + 1;
              vXYZ(1:2) = [vFirstPosition(1), vNextYZ(1) - vOverlap(2)];
              if vIndexXYZ(2) >= vDimensionY % end of column
                vIndexXYZ(2) = 0;
                vIndexXYZ(3) = vIndexXYZ(3) + 1;
                vXYZ = [vFirstPosition(1:2), vNextYZ(2) - vOverlap(3)];
                vNextYZ(1) = 0;
              end
            end
            vMax = vXYZ + vSize.*vVoxelSize;
            vCurrentExtends = [vXYZ, vMax];
          end
          
        otherwise % should not happen!
          
          ErrorBox('Invalid options')
          close(vBar)
          return
          
      end % of switch vExtendMode

      if vNumFiles == 1
        % voxelsize must be equal for all datasets
        vVoxelSize = (vCurrentExtends(4:6) - vCurrentExtends(1:3)) ./ vSize; 
        vFirstPosition = vCurrentExtends(1:3);
        
        vCurrentExtends = [0, 0, 0, vSize];
        
        vOverlap = round(vOverlap .* vVoxelSize);
      else
        % compute subblockposition
        vCurrentIndices = (vCurrentExtends(1:3) - vFirstPosition) ./ vVoxelSize;
        vCurrentIndices = round(vCurrentIndices);
 
        vCurrentExtends = [vCurrentIndices, vCurrentIndices + vSize];
      end
            
      % store extends
      vExtendList(vIndex, :) = vCurrentExtends;
    end % of loop on filenames
    
    % list of valid files shall not be empty
    if vNumFiles == 0
      ErrorBox('Please enter some VALID image filename!')
      close(vBar)
      return
    end
    
    if isequal(aImarisApplication.mView, 'eViewerSurpass') && ...
        isempty(aImarisApplication.mSurpassScene)
      aImarisApplication.mView = 'eViewerSlice';
    end
    
    % load the first file of the list
    vName = vFileList{vListIndex(1)};
    ImarisFileOpen(aImarisApplication, vName);

    % this dataset will store the result
    vMegaSet = aImarisApplication.mDataSet.Clone;

    % erase loaded volume (might be in the wrong place or have wrong intensity)
    vSize = vExtendList(vListIndex(1), 4:6);
    vVolume = 0;
    if strcmp(vMegaSet.mType, 'eTypeUInt8')
      vVolume = zeros(vSize, 'uint8');
    elseif strcmp(vMegaSet.mType, 'eTypeUInt16')
      vVolume = zeros(vSize, 'uint16');
    elseif strcmp(vMegaSet.mType, 'eTypeFloat')
      vVolume = zeros(vSize, 'single');
    end
    for vChannel = 0:vSizeCT(1)-1
      for vTime = 0:vSizeCT(2)-1
        vMegaSet.SetDataVolume(vVolume, vChannel, vTime);
      end
    end
    
    % discard unused indices
    vListIndex = vListIndex(1:vNumFiles);

    if any(vRefineStitch) || vAdjustIntensities 
      % files have to be sorted properly
      vListIndex = SortByPosition(vListIndex, vExtendList);
    end

    % discard unused extends
    vExtendList = vExtendList(vListIndex, :);
    
    % reserve space for the result
    vMin = min(vExtendList(:, 1:3), [], 1);
    vMax = max(vExtendList(:, 4:6), [], 1);
   
    vMegaSet.Resize(vMin(1), vMax(1), vMin(2), vMax(2), vMin(3), vMax(3), ...
      0, vMaxSizeCT(1), 0, vMaxSizeCT(2));

    for vIndex = 1:3
      vExtendList(:, vIndex) = vExtendList(:, vIndex) - vMin(vIndex);
      vExtendList(:, 3+vIndex) = vExtendList(:, 3+vIndex) - vMin(vIndex);
    end
    vMaxSize = vMax - vMin;
    
    % update volume extends
    vMin = vFirstPosition + vMin.*vVoxelSize;
    vMegaSet.mExtendMinX = vMin(1);
    vMegaSet.mExtendMinY = vMin(2);
    vMegaSet.mExtendMinZ = vMin(3);
    vMax = vFirstPosition + vMax.*vVoxelSize;
    vMegaSet.mExtendMaxX = vMax(1);
    vMegaSet.mExtendMaxY = vMax(2);
    vMegaSet.mExtendMaxZ = vMax(3);

    % initialize volumes position list - might be modified after refinement
    vPositions = zeros(vNumFiles, 6);
    
    % open waitbar
    waitbar(0.25, vBar, 'ImarisXT - Stitching');
    for vIndex = 1:vNumFiles
      vName = vFileList{vListIndex(vIndex)};
      
      % open dataset
      ImarisFileOpen(vImarisSlave, vName);    
      
      % compute subblockposition
      vCurrentIndices = vExtendList(vIndex, 1:3);
      vSize = vExtendList(vIndex, 4:6) - vCurrentIndices;
      
      if any(vRefineStitch)
        % compute overlaps
        vOverlaps = GetOverlaps([vCurrentIndices, vCurrentIndices] + ...
          [1, 1, 1, vSize], vPositions(1:vIndex-1, :));
        
        vTranslation = [0, 0, 0];
        % compute translations from intensities at channel 0 and time point 0 
        for vOverlapIndex = 1:size(vOverlaps, 1)
          for vChannel = 0:vSizeCT(1)-1
            for vTime = 0:vSizeCT(2)-1
              vOverlap = vOverlaps(vOverlapIndex, :);
              vOverSize = vOverlap(4:6) - vOverlap(1:3) + [1, 1, 1];
              vOverIndices = vOverlap(1:3) - vCurrentIndices - 1;
              
              vVolume1 = vDataSet.GetDataSubVolumeAs1DArray(vOverIndices(1), ...
                vOverIndices(2), vOverIndices(3), vChannel, vTime, ...
                vOverSize(1), vOverSize(2), vOverSize(3));
              vVolume1 = reshape(vVolume1, vOverSize);
              vVolume2 = vMegaSet.GetDataSubVolumeAs1DArray(vOverlap(1)-1, ...
                vOverlap(2)-1, vOverlap(3)-1, vChannel, vTime, ...
                vOverSize(1), vOverSize(2), vOverSize(3));
              vVolume2 = reshape(vVolume2, vOverSize);

              % compute refinement
              vSuggestTranslation = GetBestMatch(vVolume1, vVolume2);
              vTranslation = vTranslation + vSuggestTranslation;
            end
          end
        end
        % performs translations
        if size(vOverlaps, 1) > 0
          vTranslation = round(vTranslation / (size(vOverlaps, 1) * prod(vSizeCT)));
          vInd = vRefineStitch >= 0;
          vTranslation(vInd) = min(vTranslation(vInd), vRefineStitch(vInd));
        end
        % recompute extends
        vCurrentIndices = vCurrentIndices + vTranslation;
      
        % append at the left (negative indices)
        if min(vCurrentIndices) < 0
          vDelta = min(vCurrentIndices, 0);
          vNewSize = vMaxSize - vDelta;

          vMegaSet.Resize(vDelta(1), vNewSize(1), vDelta(2), vNewSize(2), ...
            vDelta(3), vNewSize(3), 0, vMaxSizeCT(1), 0, vMaxSizeCT(2));

          % update volume positions (all volumes have been moved!)
          for vBlock = 1:vIndex-1
            vPositions(vBlock, :) = vPositions(vBlock, :) - [vDelta, vDelta];
          end

          vCurrentIndices = max(vCurrentIndices, 0);
        end % of append at left

      end % of refine

      % update volumes position list
      vLastPos = [vCurrentIndices, vCurrentIndices] + [1, 1, 1, vSize];
      vPositions(vIndex, :) = vLastPos;
      
      % compute overlapping region for blending
      vOverlaps = GetOverlaps(vLastPos, vPositions(1:vIndex-1, :));
      vNumOverlaps = size(vOverlaps, 1);
      
      % compute new extends
      vMaxSize = max(vMaxSize, vLastPos(4:6));

      % add some space if needed
      vMegaSet.Resize(0, vMaxSize(1), 0, vMaxSize(2), 0, vMaxSize(3), ...
        0, vMaxSizeCT(1), 0, vMaxSizeCT(2));

      % copy the "slave" dataset to the main Imaris instance
      for vChannel = 0:vSizeCT(1)-1
        for vTime = 0:vSizeCT(2)-1
          vVolume = vDataSet.GetDataVolumeAs1DArray(vChannel, vTime);
          vVolume = reshape(vVolume, vSize);
          
          % adjust intensities
          if vAdjustIntensities && vNumOverlaps>0
            vAdjustMean = zeros(1, vNumOverlaps);
            vAdjustStd = ones(1, vNumOverlaps);
            
            % compare overlapping regions
            for vOverlapIndex = 1:vNumOverlaps 
              vOverlap = vOverlaps(vOverlapIndex, :);
              vOverSize = vOverlap(4:6) - vOverlap(1:3) + [1, 1, 1];
              vOverIndices = vOverlap(1:3) - vCurrentIndices - 1;

              vVolume2 = vMegaSet.GetDataSubVolumeAs1DArray(vOverlap(1)-1, ...
                vOverlap(2)-1, vOverlap(3)-1, vChannel, vTime, ...
                vOverSize(1), vOverSize(2), vOverSize(3));

              vVolume1 = vVolume(vOverIndices(1)+(1:vOverSize(1)), ...
                vOverIndices(2)+(1:vOverSize(2)), vOverIndices(3)+(1:vOverSize(3)));

              vMean = mean(double(vVolume1(:)));
              vAdjustMean(vOverlapIndex) = mean(double(vVolume2)) - vMean;
              vStdVolume = std(double(vVolume1(:)));
              vAdjustStd(vOverlapIndex) = std(double(vVolume2)) / ...
                (vStdVolume+(vStdVolume==0)); 
            end % of loop on overlaps
            
            vMean = mean(double(vVolume(:)));
            vAdjustMean = mean(vAdjustMean);
            vAdjustStd = mean(vAdjustStd);
            % rescale intensities
            vVolume = (double(vVolume) - vMean)*vAdjustStd + vMean + vAdjustMean;
            
            if strcmp(vMegaSet.mType, 'eTypeUInt8')
              vVolume = uint8(vVolume);
            elseif strcmp(vMegaSet.mType, 'eTypeUInt16')
              vVolume = uint16(vVolume);
            elseif strcmp(vMegaSet.mType, 'eTypeFloat')
              vVolume = single(vVolume);
            end
          end % of adjust intensities

          % blend overlapping regions
          for vOverlapIndex = 1:vNumOverlaps
            vOverlap = vOverlaps(vOverlapIndex, :);
            vOverSize = vOverlap(4:6) - vOverlap(1:3) + [1, 1, 1];
            vOverIndices = vOverlap(1:3) - vCurrentIndices - 1;

            vVolume2 = vMegaSet.GetDataSubVolumeAs1DArray(vOverlap(1)-1, ...
              vOverlap(2)-1, vOverlap(3)-1, vChannel, vTime, ...
              vOverSize(1), vOverSize(2), vOverSize(3));
            vVolume2 = reshape(vVolume2, vOverSize);
          
            vVolume1 = vVolume(vOverIndices(1)+(1:vOverSize(1)), ...
              vOverIndices(2)+(1:vOverSize(2)), vOverIndices(3)+(1:vOverSize(3)));
                      
            if ~strcmp(vMegaSet.mType, 'eTypeFloat')
              % int division rounds by excess
              % 5/2 + 5/2 -> would return 3 + 3 = 6
              vBothOdd = (mod(vVolume1, 2) == 1) & (mod(vVolume2, 2) == 1);
              vVolume1(vBothOdd) = vVolume1(vBothOdd) - 1;
            end
            
            % average volumes
            vVolume(vOverIndices(1)+(1:vOverSize(1)), ...
              vOverIndices(2)+(1:vOverSize(2)), ...
              vOverIndices(3)+(1:vOverSize(3))) = ...
              vVolume1 / 2 + vVolume2 / 2;
          end % of blend overlapping regions
          
          vMegaSet.SetDataSubVolume(vVolume, vCurrentIndices(1), ...
            vCurrentIndices(2), vCurrentIndices(3), vChannel, vTime);
        end % of loop on channels
      end % of loop on timepoints
      
      waitbar(0.25 + vIndex/vNumFiles*0.75, vBar)
    end % of loop on files
    
    if any(vRefineStitch)
      % discard useless border
      vMin = min(vPositions(:, 1:3), [], 1) - 1;
      vSize = max(vPositions(:, 4:6), [], 1) - vMin;
      vMegaSet.Resize(vMin(1), vSize(1), vMin(2), vSize(2), vMin(3), vSize(3), ...
        0, vMaxSizeCT(1), 0, vMaxSizeCT(2));
    end
    
    % time to show the result
    aImarisApplication.mDataSet = vMegaSet;
    
    % close the waitbar
    close(vBar)

    % free some memory
    if ~isequal(vImarisSlave, aImarisApplication)
      vDataSet.Resize(0, 1, 0, 1, 0, 1, 0, 1, 0, 1);
      vImarisSlave.mDataSet = vDataSet;
    end
    
    % restore enable / disable uicontrols
    for vIndex = 1:numel(vFields)
      vFieldName = vFields{vIndex};
      if numel(vFieldName) < vPanelTextLenght || ...
          ~isequal(vFieldName(1:vPanelTextLenght), vPanelText)
        set(aUIControls.(vFieldName), 'Enable', 'on');
      end
    end
    aOptions = UpdateToggleLayout(aUIControls, aOptions);
    UpdateRadioEnable(aUIControls, aOptions);
    
%-------------------------------------------------------------------------%
%% Resize Layout table to fit the window
%-------------------------------------------------------------------------%
  case 3 % Resize window
    
    [aUIControls, aOptions] = ResizeWindow(aUIControls, aOptions);
    UpdateSliders(aUIControls, aOptions);
    UpdateGrid(aUIControls, aOptions);
    
%-------------------------------------------------------------------------%
%% Open - Show dialog
%-------------------------------------------------------------------------%
  case 4 % Open...
    
    vFileName = [aOptions.mFilePath, '\'];
    [vFileName, vPathName] = uigetfile({'*.*', 'Image-files'}, ...
      'Add image file', vFileName);
    if ~isequal(vFileName, 0)
      vFileName = fullfile(vPathName, vFileName);
      aOptions.mFilePath = vPathName;
      aOptions.mFileName = vFileName;
      set(aUIControls.mEditFileName, 'String', vFileName);
      
      [aUIControls, aOptions] = Process(5, ...
        aImarisApplication, aUIControls, aOptions);
    end

%-------------------------------------------------------------------------%
% Edit filename - Test if valid
%-------------------------------------------------------------------------%
  case 5 
    
    aOptions.mFileName = get(aUIControls.mEditFileName, 'String');
    vValid = exist(aOptions.mFileName, 'file');

    vFiles = {};
    vIndexSize = [];
    vCoords = [];
    if vValid
      [vFiles, vCoords, vIndexSize] = ReadTextFile(aOptions.mFileName);
      if isempty(vIndexSize)
        vIndexSize = [1, numel(vFiles), 1];
      end
      
      aOptions.mTextFile = ~isempty(vFiles);
      UpdateRadioEnable(aUIControls, aOptions);
      
      if ~aOptions.mTextFile
        [vPath, vName, vExt, vVersion] = fileparts(aOptions.mFileName);
        aOptions.mFilePath = vPath;
        vFiles = GetFileNames(vPath, [vName, vExt, vVersion]);
        vIndexSize = [1, numel(vFiles), 1];
        if aOptions.mLayoutFromFileNames
          [vFiles, vIndexSize] = SortByIndex(vFiles);
        end
      else
        set(aUIControls.mRadioFileNames, 'Value', isempty(vCoords));
        aOptions = UpdateToggleLayout(aUIControls, aOptions);
      end
    end

    vStates = {'on', 'off'};
    set(aUIControls.mEditFileName, 'ForegroundColor', [isempty(vIndexSize), 0, 0]);
    set(aUIControls.mButtonStitch, 'Enable', vStates{1+isempty(vIndexSize)});
    
    if aOptions.mLayoutFromFileNames && ~isempty(vIndexSize) && aOptions.mFlipY
      vIndices = 1:prod(vIndexSize);
      vIndices = reshape(vIndices, vIndexSize);
      vIndices = flipdim(vIndices, 2);
      vFiles = vFiles(vIndices(:));
    end

    if isempty(vIndexSize)
      vIndexSize = [1, numel(vFiles), 1];
    end
    
    if numel(vFiles) > 0
      % create "slave" Imaris instance to load single datasets
      vImarisSlave = aOptions.mImarisSlave;
      if isempty(vImarisSlave)
        vImarisSlave = actxserver('Imaris.Application');
        vImarisSlave.mView = 'eViewerSlice';
        aOptions.mImarisSlave = vImarisSlave;
      end
      vDataSet = vImarisSlave.mDataSet;
    
      vName = fullfile(aOptions.mFilePath, vFiles{1});
      % Imaris can read it?
      vSuccess = ImarisFileOpen(vImarisSlave, vName);
      if vSuccess
        aOptions.mEditShiftX = round(vDataSet.mSizeX / 5);
        aOptions.mEditShiftY = round(vDataSet.mSizeY / 5);
        aOptions.mEditShiftZ = round(vDataSet.mSizeZ / 5);
        set(aUIControls.mEditShiftX, 'String', sprintf('%d', aOptions.mEditShiftX));
        set(aUIControls.mEditShiftY, 'String', sprintf('%d', aOptions.mEditShiftY));
        set(aUIControls.mEditShiftZ, 'String', sprintf('%d', aOptions.mEditShiftZ));

        if ~isequal(vImarisSlave, aImarisApplication)
          vDataSet.Resize(0, 1, 0, 1, 0, 1, 0, 1, 0, 1);
          vImarisSlave.mDataSet = vDataSet;
        end
      end
    end
    
    aOptions.mFiles = vFiles;
    aOptions.mIndexSize = vIndexSize;
    aOptions.mCoords = vCoords;
    
    [aUIControls, aOptions] = ResizeWindow(aUIControls, aOptions);
    UpdateSliders(aUIControls, aOptions);
    UpdateGrid(aUIControls, aOptions);

%-------------------------------------------------------------------------%
% Toggle layout from ...
%-------------------------------------------------------------------------%
  case 6 
    
    aOptions = UpdateToggleLayout(aUIControls, aOptions);
    [aUIControls, aOptions] = Process(5, ...
        aImarisApplication, aUIControls, aOptions);
      
%-------------------------------------------------------------------------%
% Flip Y
%-------------------------------------------------------------------------%
  case 7
    
    aOptions.mFlipY = get(aUIControls.mCheckFlipY, 'Value');
    [aUIControls, aOptions] = Process(5, ...
        aImarisApplication, aUIControls, aOptions);
    
%-------------------------------------------------------------------------%
% Adjust intensities
%-------------------------------------------------------------------------%
  case 8
    
    aOptions.mAdjust = get(aUIControls.mCheckAdjust, 'Value');

%-------------------------------------------------------------------------%
% Edit overlap or max shift
%-------------------------------------------------------------------------%
  case 9
    
    vFields = {'mEditOverlapX', 'mEditOverlapY', 'mEditOverlapZ', ...
      'mEditShiftX', 'mEditShiftY', 'mEditShiftZ'};
    for vIndex = 1:numel(vFields)
      vField = vFields{vIndex};
      
      % read value
      vValue = str2double(get(aUIControls.(vField), 'String'));
      vValue = round(vValue);
      if ~isnan(vValue)
        aOptions.(vField) = vValue;
      else % invalid
        % restore
        vValue = aOptions.(vField);
      end
      set(aUIControls.(vField), 'String', num2str(vValue));
    end
    
%-------------------------------------------------------------------------%
% Scrolls
%-------------------------------------------------------------------------%
  case 10
    
    UpdateGrid(aUIControls, aOptions);

%-------------------------------------------------------------------------%
% Someone did a joke!
%-------------------------------------------------------------------------%
  otherwise
    
    % should not happen

end


%-------------------------------------------------------------------------%
% Display error
%-------------------------------------------------------------------------%
function ErrorBox(aText)
msgbox(aText)


%-------------------------------------------------------------------------%
% Try to open a file
%-------------------------------------------------------------------------%
function aSuccess = ImarisFileOpen(aImarisApplication, aFileName)

[vPath, vName, vExt, vVersion] = fileparts(aFileName);

vTempFileName = aFileName;
vNeedRename = ~isequal(vExt, '.ims');
if vNeedRename
  vTempFileName = fullfile(vPath, ['imarisstitchtempfile', vExt, vVersion]);
  copyfile(aFileName, vTempFileName);
end

aSuccess = 1;
try
  aImarisApplication.FileOpen(vTempFileName);
catch
  warning('MATLAB:invalidFile', 'Could not load %s', aFileName)
  aSuccess = 0;
end

if vNeedRename
  delete(vTempFileName);
end


%-------------------------------------------------------------------------%
%% Compute overlapping regions
%-------------------------------------------------------------------------%
function aOverlaps = GetOverlaps(aRange, aRangeList)

% aRange = [minX, minY, minZ, maxX, maxY, maxZ];
% aRangeList = [minX1, minY1, minZ1, maxX1, maxY1, maxZ1; ...
%               minX2, minY2, minZ2, maxX2, maxY2, maxZ2; ...
% aOverlaps has the same structure as aRangeList (does not contain empty entries)

aOverlaps = aRangeList;
vNumOverlaps = 0;
for vIndex = 1:size(aRangeList, 1)
  vRange2 = aRangeList(vIndex, :);
  vOverlap = [max(aRange(1:3), vRange2(1:3)), min(aRange(4:6), vRange2(4:6))];
  if min(vOverlap(4:6) - vOverlap(1:3)) >= 0
    vNumOverlaps = vNumOverlaps + 1;
    aOverlaps(vNumOverlaps, :) = vOverlap;
  end
end
aOverlaps = aOverlaps(1:vNumOverlaps, :);


%-------------------------------------------------------------------------%
%% Compute shift so that VolumeShift(aVolume1) == aVolume2
%-------------------------------------------------------------------------%
function aDisplace = GetBestMatch(aVolume1, aVolume2)

aDisplace = [0, 0, 0];

vSize = size(aVolume1);
if ~isequal(size(aVolume2), vSize)
  warning('MATLAB:invalidMatrixSize', 'Invalid matrix: Skipping refinement');
  return
end

% calculate correlation in Fourier space
vCorrelation = fourierCorrND(aVolume1, aVolume2);

% calculate the shift
[vMax, vIndices] = max(vCorrelation(:));

% computes "[vD1, vD2, vD3] = ind2sub(vSize, vIndices);" 3-4 times faster
vD1 = mod(vIndices-1, vSize(1)) + 1;
vIndices = (vIndices - vD1) / vSize(1) + 1;
vD2 = mod(vIndices-1, vSize(2)) + 1;
vD3 = (vIndices - vD2) / vSize(2) + 1;

% offset for the center of the spectrum
vOffset = [fix(vSize/2) + 1, ones(1, 3 - numel(vSize))];

% result
aDisplace = [vD1, vD2, vD3] - vOffset;

% debug message
% disp(sprintf('Refined %i, %i, %i', aDisplace(1), aDisplace(2), aDisplace(3)))


%-------------------------------------------------------------------------%
%% Sort the tiles
%-------------------------------------------------------------------------%
function aListIndex = SortByPosition(aListIndex, aExtendList)

% discard useless data
aExtendList = aExtendList(aListIndex, :);

vNumEl = numel(aListIndex);

% routine will return the elements of aListIndex in order vOrder
vOrder = 1:vNumEl;

vCenters = (aExtendList(:, 1:3) + aExtendList(:, 4:6)) / 2; 
vMeanCenter = median(vCenters, 1);

% compute overlap size only once
vOverlaps = ones(vNumEl);

vDist = zeros(vNumEl, 1);
for vDim = 1:3
  % compute the distance from the center of the scene
  vDist = sqrt(vDist.^2 + (vCenters(:, vDim) - vMeanCenter(vDim)).^2);

  % overlap size in a single dimension
  vThisOvers = zeros(vNumEl);
  for vIndex = 1:vNumEl
    % 2 segments: A-B, C-D
    vA = aExtendList(:, vDim);
    vB = aExtendList(:, 3 + vDim);
    vC = aExtendList(vIndex, vDim) * ones(vNumEl, 1);
    vD = aExtendList(vIndex, 3 + vDim) * ones(vNumEl, 1);
    
    vMinOfMaxes = min(vB, vD);
    vMaxOfMins = max(vA, vC);
    
    % zero if there is no overlap
    vThisOvers(vIndex, :) = max(vMinOfMaxes - vMaxOfMins, 0); 
  end
  
  % multiply each dimension to get the volume
  vOverlaps = vOverlaps .* vThisOvers;
end

% overlap between tile A and tile A must be zero
vOverlaps = vOverlaps - diag(diag(vOverlaps));

% get first tile
[vDist, vCurrent] = min(vDist);

% swap indices
vOrder(vCurrent) = 1;
vOrder(1) = vCurrent;

% swap overlaps
vOverlaps = SwapLineColumn(vOverlaps, 1, vCurrent);

% selection sort with different scoring function at each selection
for vIndex = 2:vNumEl
  % sum all already placed overlaps
  vScore = sum(vOverlaps(1:vIndex-1, vIndex:vNumEl), 1);
  
  % get the new tile to place
  [vScore, vCurrent] = max(vScore);
  vCurrent = vIndex - 1 + vCurrent;

  % swap indices
  vTemp = vOrder(vCurrent);
  vOrder(vCurrent) = vOrder(vIndex);
  vOrder(vIndex) = vTemp;
  
  % swap overlaps
  vOverlaps = SwapLineColumn(vOverlaps, vIndex, vCurrent);
end

aListIndex = aListIndex(vOrder);


%-------------------------------------------------------------------------%
%% Swap line A with line B and column A with column B
%-------------------------------------------------------------------------%
function aMatrix = SwapLineColumn(aMatrix, aIndexA, aIndexB)
% aMatrix has to be simmetric
if ndims(aMatrix)~=2 || size(aMatrix, 1)~=size(aMatrix, 2) || ...
    any(any(aMatrix - aMatrix'))
  return
end

% swap lines
vLine = aMatrix(aIndexA, :);
aMatrix(aIndexA, :) = aMatrix(aIndexB, :);
aMatrix(aIndexB, :) = vLine;

% swap columns
vLine = aMatrix(:, aIndexA); 
aMatrix(:, aIndexA) = aMatrix(:, aIndexB);
aMatrix(:, aIndexB) = vLine;

% PS: aMatrix(aIndexA, aIndexB) remains unchanged (w.a.i.)


%-------------------------------------------------------------------------%
%% Sort filenames according to indices _X _Y _Z
%-------------------------------------------------------------------------%
function [aFileList, aSize] = SortByIndex(aFileList)

% try some naming conventions
vIndices = zeros(numel(aFileList), 3, 3);
for vIndex = 1:numel(aFileList)
  vName = aFileList{vIndex};
  vIndices(vIndex, :, 1) = ExtractIndices(vName, {'_X', '_Y', '_Z'});
  vIndices(vIndex, :, 2) = ExtractIndices(vName, {'_x', '_y', '_z'});
  vIndices(vIndex, :, 3) = ExtractIndices(vName, {'_', '_', '_'});
end

% choose the best convention
vMissing = squeeze(sum(sum(vIndices == -1, 1), 2));
[vValue, vIndex] = min(vMissing);
vIndices = max(vIndices(:, :, vIndex), 0);

aSize = max(vIndices, [], 1) - min(vIndices, [], 1) + 1;
vPerfectIndices = size(vIndices, 1) == prod(aSize);
if vPerfectIndices
  vPerfectIndices = PerfectIndices(vIndices);
end

% convert 3 indices to one, priority is z, y, x
vIndices = vIndices(:, 1) + aSize(1) * vIndices(:, 2) + ...
  aSize(1) * aSize(2) * vIndices(:, 3);

% sort the file list
[vValue, vIndices] = sort(vIndices);
aFileList = aFileList(vIndices);

if ~vPerfectIndices
  aSize = [];
end


%-------------------------------------------------------------------------%
%% Extract an index for each tag, return -1 if not found
%-------------------------------------------------------------------------%
function aIndices = ExtractIndices(aFileName, aTags)

% ExtractIndices('a_x101_5_y2_x3_y.ims', {'_z', '_y','_x', '_x', 'argh'})
%   -> -1, 2, 101, 3, -1

aIndices = zeros(1, numel(aTags));
% process in reverse order, since the string is processed from right
for vIndex = numel(aTags):-1:1
  vTag = aTags{vIndex};
  vElem = (1:numel(vTag)) - 1;
  
  % do not restart if we search again the same tag 
  if vIndex == numel(aTags) || ~isequal(vTag, aTags{vIndex+1})
    vPosition = numel(aFileName) - numel(vTag);
  end
  
  vValue = -1;
  while vPosition > 0 && vValue == -1

    if isequal(aFileName(vPosition + vElem), vTag)
      % read the value
      vHere = vPosition + numel(vTag);
      vInteger = str2double(aFileName(vHere));
      while ~isnan(vInteger)
        vValue = max(vValue * 10, 0);
        vValue = vValue + vInteger;
        vHere = vHere + 1;
        vInteger = str2double(aFileName(vHere));
      end
    end

    vPosition = vPosition - 1;
  end
  aIndices(vIndex) = vValue;
end

aIndices = [aIndices, zeros(3-numel(aIndices), 1)];


%-------------------------------------------------------------------------%
%% Return true if the indices are on a grid
%-------------------------------------------------------------------------%
function aResult = PerfectIndices(aIndices)

vX = unique(aIndices(:, 1));
vY = unique(aIndices(:, 2));
vZ = unique(aIndices(:, 3));

aResult = size(aIndices, 1) == numel(vX) * numel(vY) * numel(vZ);
if aResult
  % number is ok, now test the content
  for vIndex = 1:size(aIndices, 1)-1
    for vSecond = vIndex+1:size(aIndices, 1)
      if isequal(aIndices(vIndex, :), aIndices(vSecond, :))
        aResult = 0;
        return
      end
    end
  end
end
  

%-------------------------------------------------------------------------%
%% Return the file names of the dataset
%-------------------------------------------------------------------------%
function aFileNames = GetFileNames(aPath, aReferenceFile)

vFiles = dir(aPath);
aFileNames = cell(numel(vFiles), 1);
vNumFiles = 0;
vBaseReference = RemoveNumbers(aReferenceFile);
for vIndex = 1:numel(vFiles)
  if vFiles(vIndex).isdir
    continue      
  end
  vName = vFiles(vIndex).name;
  vBase = RemoveNumbers(vName);
  if isequal(vBase, vBaseReference)
    vNumFiles = vNumFiles + 1;
    aFileNames{vNumFiles} = vName;
  end
end
aFileNames = aFileNames(1:vNumFiles);


%-------------------------------------------------------------------------%
%% Removes characters from '0' to '9' from a string
%-------------------------------------------------------------------------%
function aString = RemoveNumbers(aString)

aString = aString(aString > '9' | aString < '0');


%-------------------------------------------------------------------------%
%% Read filenames and coords from .txt file
%-------------------------------------------------------------------------%
function [aFiles, aCoords, aIndexSize] = ReadTextFile(aFileName)

vFile = fopen(aFileName, 'rt');

aFiles = {};
aCoords = [];
aIndexSize = [];
% read xml file
vXml = fread(vFile);
fclose(vFile);

vXml = Trim(vXml(vXml ~= sprintf('\n')));
vValue = ReadXmlTags(vXml, 'imaris_stitch');
if isempty(vValue)
  return
end
vXml = vValue{1}.mValue;
  
vFileTags = ReadXmlTags(vXml, 'file');

vFiles = cell(numel(vFileTags), 1);
vPosition = zeros(numel(vFileTags), 3);

vPrefix = {'coord', 'index'};
vSuffix = {'x', 'y', 'z'};

for vIndex = 1:numel(vFileTags)
  vAttrs = vFileTags{vIndex}.mAttributes;
  if numel(vAttrs) ~= 4
    return
  end
  vValue = ReadValueFromAttrs(vAttrs, 'name');
  if isnan(vValue)
    return
  end
  vFiles{vIndex} = vValue;
  if iscell(vPrefix)
    vValue = str2double(ReadValueFromAttrs(vAttrs, [vPrefix{1}, vSuffix{1}]));    
    vPrefix = vPrefix{1 + isnan(vValue)};
  end
  for vDim = 1:numel(vSuffix)
    vValue = str2double(ReadValueFromAttrs(vAttrs, [vPrefix, vSuffix{vDim}]));    
    if isnan(vValue)
      return
    end
    vPosition(vIndex, vDim) = vValue;
  end
end

if isequal(vPrefix, 'coord')
  aFiles = vFiles;
  aCoords = vPosition;
  aIndexSize = [1, numel(vFileTags), 1];
else
  vIndices = vPosition;
  
  vSize = max(vIndices, [], 1) - min(vIndices, [], 1) + 1;
  vPerfectIndices = size(vIndices, 1) == prod(vSize);
  if vPerfectIndices
    vPerfectIndices = PerfectIndices(vIndices);
  end
  if ~vPerfectIndices
    return
  end
  
  % convert 3 indices to one, priority is z, y, x
  vIndices = vIndices(:, 1) + vSize(1) * vIndices(:, 2) + ...
    vSize(1) * vSize(2) * vIndices(:, 3);

  % sort the file list
  [vFileTags, vIndices] = sort(vIndices);

  aFiles = vFiles(vIndices);
  aIndexSize = vSize;
  % aCoords = [];
end
  

%-------------------------------------------------------------------------%
%% Return value from attributes
%-------------------------------------------------------------------------%
function aValue = ReadValueFromAttrs(aAttrs, vName)

aValue = NaN;
for vIndex = 1:numel(aAttrs)
  if isequal(aAttrs{vIndex}.mName, vName)
    aValue = aAttrs{vIndex}.mValue;
    return
  end
end


%-------------------------------------------------------------------------%
%% Return list of struct .mValue (string), .mAttributes (list of struct
%    .mName, .mValue)
%-------------------------------------------------------------------------%
function aValue = ReadXmlTags(aXml, aTag)

aValue = {};

vStart = aXml == '<';
vEnd = aXml == '>';
vNumberTags = sum(vStart);
if sum(vEnd) ~= vNumberTags
  return
end

vOpenStart = 1:vNumberTags;
vOpenIndex = 1:vNumberTags;
vNumTags = 0;
vNumTagsOpen = 0;

vStart = find(vStart);
vEnd = find(vEnd);
for vIndex = 1:vNumberTags
  vTag = Trim(aXml(vStart(vIndex)+1:vEnd(vIndex)-1));
  vType = 0;
  if vTag(numel(vTag)) == '/' 
    % short tag
    vType = 1;
    vTag = Trim(vTag(1:numel(vTag)-1));
  elseif vTag(1) == '/'
    % close tag
    vType = 2;   
    vTag = Trim(vTag(2:numel(vTag)));
  end
  vSep = find(vTag == ' ', 1);
  if isempty(vSep)
    vSep = numel(vTag) + 1;
  end
  vName = Trim(vTag(1:vSep-1));
  if isequal(vName, aTag)
    if vType < 2
      vNumTags = vNumTags + 1;
      aValue{vNumTags} = struct('mValue', '', ...
        'mAttributes', {ReadAttributes(vTag(vSep+1:numel(vTag)))});
    end
    
    if vType == 0
      vNumTagsOpen = vNumTagsOpen + 1;
      vOpenIndex(vNumTagsOpen) = vNumTags;
      vOpenStart(vNumTagsOpen) = vEnd(vIndex) + 1;
    elseif vType == 2
      if vNumTagsOpen < 1
        return
      end
      vTagIndex = vOpenIndex(vNumTagsOpen);
      vTagStart = vOpenStart(vNumTagsOpen);
      aValue{vTagIndex}.mValue = Trim(aXml(vTagStart:vStart(vIndex)-1));
      vNumTagsOpen = vNumTagsOpen - 1;
    end
  end
end


%-------------------------------------------------------------------------%
%% Read from xml a list of struct .mName, .mValue
%-------------------------------------------------------------------------%
function aAttributes = ReadAttributes(aXml)

aXml = Trim(aXml);
aAttributes = {};
vNum = 0;

while ~isempty(aXml)
  vIndex = find(aXml == '=', 1);
  if isempty(vIndex)
    return
  end
  vName = Trim(aXml(1:vIndex-1));
  aXml = Trim(aXml(vIndex+1:numel(aXml)));
  
  vIndex = find(aXml == '"', 1);
  if isempty(vIndex)
    return
  end
  aXml = Trim(aXml(vIndex+1:numel(aXml)));

  vIndex = find(aXml == '"', 1);
  if isempty(vIndex)
    return
  end
  vValue = Trim(aXml(1:vIndex-1));
  aXml = Trim(aXml(vIndex+1:numel(aXml)));
  
  vNum = vNum + 1;
  aAttributes{vNum} = struct('mName', vName, 'mValue', vValue);
end


%-------------------------------------------------------------------------%
%% Remove white spaces at begin and end of string
%-------------------------------------------------------------------------%
function aString = Trim(aString)

vBool = aString ~= ' ';
vMin = find(vBool, 1, 'first');
vMax = find(vBool, 1, 'last');
if ~isempty(vMin) && ~isempty(vMax)
  aString = char(aString(vMin:vMax));
else
  aString = '';
end
if size(aString, 1) > size(aString, 2)
  aString = aString';
end


%-------------------------------------------------------------------------%
%% Layout from text file?
%-------------------------------------------------------------------------%
function UpdateRadioEnable(aUIControls, aOptions)

vStates = {'on', 'off'};
vState = vStates{1 + aOptions.mTextFile};
set(aUIControls.mLabelLayout, 'Enable', vState);
set(aUIControls.mRadioFileNames, 'Enable', vState);
set(aUIControls.mRadioFilesContent, 'Enable', vState);

%-------------------------------------------------------------------------%
%% Layout from file names or files content?
%-------------------------------------------------------------------------%
function aOptions = UpdateToggleLayout(aUIControls, aOptions)

vFromNames = get(aUIControls.mRadioFileNames, 'Value');
if aOptions.mLayoutFromFileNames ~= vFromNames
  set(aUIControls.mRadioFilesContent, 'Value', ~vFromNames);
else
  vFromContent = get(aUIControls.mRadioFilesContent, 'Value');
  set(aUIControls.mRadioFileNames, 'Value', ~vFromContent);
end
aOptions.mLayoutFromFileNames = get(aUIControls.mRadioFileNames, 'Value');
vChildren = {'mCheckFlipY', 'mLabelOverlap', ...
  'mEditOverlapX', 'mEditOverlapY', 'mEditOverlapZ'};
vStates = {'off', 'on'};
vEnable = vStates{aOptions.mLayoutFromFileNames+1};
for vIndex = 1:numel(vChildren)
  vChild = vChildren{vIndex};
  set(aUIControls.(vChild), 'Enable', vEnable);
end


%-------------------------------------------------------------------------%
%% Filenames in the grid
%-------------------------------------------------------------------------%
function UpdateGrid(aUIControls, aOptions)

vX = 0;
vY = 0;
if isequal(get(aUIControls.mSliderLayoutH, 'Visible'), 'on') 
  vX = round(get(aUIControls.mSliderLayoutH, 'Value'));
end
if isequal(get(aUIControls.mSliderLayoutV, 'Visible'), 'on') 
  vY = get(aUIControls.mSliderLayoutV, 'Max') - round(get(aUIControls.mSliderLayoutV, 'Value'));
end

if aOptions.mLayoutFromFileNames
  for vIndex = 3:min(aOptions.mGridX, aOptions.mIndexSize(1)+2)
    set(aUIControls.mPanelLayoutLabels(vIndex), ...
      'String', sprintf('x = %d', vIndex+vX-2));  
  end
  for vIndex = 2:min(aOptions.mGridY, prod(aOptions.mIndexSize(2:3)) + 1)
    vMod = mod(vIndex+vY-2, aOptions.mIndexSize(2));
    if vIndex == 2 || vMod == 0
      vNum = (vIndex+vY-2-vMod)/aOptions.mIndexSize(2);
      set(aUIControls.mPanelLayoutLabels((vIndex-1)*aOptions.mGridX+1), ...
        'String', sprintf('z = %d', vNum+1));
    else
      set(aUIControls.mPanelLayoutLabels((vIndex-1)*aOptions.mGridX+1), ...
        'String', '');
    end
    set(aUIControls.mPanelLayoutLabels((vIndex-1)*aOptions.mGridX+2), ...
      'String', sprintf('y = %d', vMod+1));      
  end
else
  for vIndex = 1:min(aOptions.mGridY-1, prod(aOptions.mIndexSize(2:3)))
    set(aUIControls.mPanelLayoutLabels(vIndex*aOptions.mGridX+2), ...
      'String', sprintf('n = %d', vIndex+vY));      
  end  
end

for vIndexY = 2:min(aOptions.mGridY, prod(aOptions.mIndexSize(2:3)) + 1)
  for vIndexX = 3:min(aOptions.mGridX, aOptions.mIndexSize(1) + 2)
    vFile = aOptions.mFiles(vIndexX+vX-2 + (vIndexY+vY-2)*aOptions.mIndexSize(1));
    set(aUIControls.mPanelLayoutLabels(vIndexX + (vIndexY-1)*aOptions.mGridX), ...
      'String', vFile);
  end
end

vPosition = get(aUIControls.mPanelLayoutLabels(1), 'Position');
vSizeY = vPosition(4);
for vIndexY = 2:min(aOptions.mGridY, prod(aOptions.mIndexSize(2:3)) + 1)
  vIndex = (vIndexY-1)*aOptions.mGridX + 1;
  vPosition = get(aUIControls.mPanelLayoutLabels(vIndex), 'Position');
  vPosition(4) = vSizeY + (mod(vIndexY+vY-2, aOptions.mIndexSize(2)) ~= 0);
  set(aUIControls.mPanelLayoutLabels(vIndex), 'Position', vPosition);
end


%-------------------------------------------------------------------------%
%% Update sliders
%-------------------------------------------------------------------------%
function UpdateSliders(aUIControls, aOptions)

vIndexSize = aOptions.mIndexSize;

set(aUIControls.mSliderLayoutH, 'Visible', 'off');
vN = vIndexSize(1) - aOptions.mGridX + 2;
if vN > 0
  set(aUIControls.mSliderLayoutH, 'Max', vN);    
  set(aUIControls.mSliderLayoutH, 'SliderStep', 1/vN*[1, 3]);    
  set(aUIControls.mSliderLayoutH, 'Value', 0);
  set(aUIControls.mSliderLayoutH, 'Visible', 'on');    
end
set(aUIControls.mSliderLayoutV, 'Visible', 'off');
vN = prod(vIndexSize(2:3)) - aOptions.mGridY + 1;
if vN > 0
  set(aUIControls.mSliderLayoutV, 'Max', vN);    
  set(aUIControls.mSliderLayoutV, 'SliderStep', 1/vN*[1, 3]);    
  set(aUIControls.mSliderLayoutV, 'Value', vN);
  set(aUIControls.mSliderLayoutV, 'Visible', 'on');    
end

vStates = {'off', 'on'};
for vIndexY = 1:aOptions.mGridY
  for vIndexX = 1:aOptions.mGridX
    vIndex = vIndexX + (vIndexY-1)*aOptions.mGridX;
    set(aUIControls.mPanelLayoutLabels(vIndex), 'String', '');
    vIsVisible = (vIndexX < aOptions.mIndexSize(1)+3) && ...
      (vIndexY < prod(aOptions.mIndexSize(2:3))+2) && (vIndex > 2);
    set(aUIControls.mPanelLayoutLabels(vIndex), ...
      'Visible', vStates{1+vIsVisible});
  end
end


%-------------------------------------------------------------------------%
%% Builds the GUI, returns the list of the objects handles
%-------------------------------------------------------------------------%
function aUIControls = BuildGUI()

% create the window
%vWindow = ...
figure('Name', 'Stitcher XTension', 'NumberTitle', 'off', ...
  'MenuBar', 'none', 'ToolBar', 'none', ...
  'ResizeFcn', 'BPStitchDatasets(-3)', ...
  'DeleteFcn', 'BPStitchDatasets(-1)');
%aUIControls.mWindow = vWindow;

aUIControls.mPanelStitch = uipanel('Title', 'Stitching');
aUIControls.mPanelLayout = uipanel('Title', 'Layout');
aUIControls.mPanelSource = uipanel('Title', 'Source');

% stitching
aUIControls.mLabelShift = uicontrol('Style', 'text', 'String', 'Max shift', ...
  'Parent', aUIControls.mPanelStitch);
aUIControls.mEditShiftX = uicontrol('Style', 'edit', 'String', '0', ...
  'Parent', aUIControls.mPanelStitch, 'Callback', 'BPStitchDatasets(-9)');
aUIControls.mEditShiftY = uicontrol('Style', 'edit', 'String', '0', ...
  'Parent', aUIControls.mPanelStitch, 'Callback', 'BPStitchDatasets(-9)');
aUIControls.mEditShiftZ = uicontrol('Style', 'edit', 'String', '0', ...
  'Parent', aUIControls.mPanelStitch, 'Callback', 'BPStitchDatasets(-9)');
aUIControls.mCheckAdjust = uicontrol('Style', 'check', 'String', 'Adjust intensities', ...
  'Value', 0, ...
  'Parent', aUIControls.mPanelStitch, 'Callback', 'BPStitchDatasets(-8)');
aUIControls.mButtonStitch = uicontrol('Style', 'pushbutton', 'String', 'Stitch!', ...
  'FontWeight', 'bold', 'BackgroundColor', [1, 1, 1], ...
  'Parent', aUIControls.mPanelStitch, 'Callback', 'BPStitchDatasets(-2)');

% layout
aUIControls.mPanelLayoutLabels = [];
aUIControls.mLabelOverlap = uicontrol('Style', 'text', 'String', 'Overlap', ...
  'Parent', aUIControls.mPanelLayout);
aUIControls.mEditOverlapX = uicontrol('Style', 'edit', 'String', '0', ...
  'Parent', aUIControls.mPanelLayout, 'Callback', 'BPStitchDatasets(-9)');
aUIControls.mEditOverlapY = uicontrol('Style', 'edit', 'String', '0', ...
  'Parent', aUIControls.mPanelLayout, 'Callback', 'BPStitchDatasets(-9)');
aUIControls.mEditOverlapZ = uicontrol('Style', 'edit', 'String', '0', ...
  'Parent', aUIControls.mPanelLayout, 'Callback', 'BPStitchDatasets(-9)');
aUIControls.mLabelUnits = uicontrol('Style', 'text', 'String', '[Voxels]', ...
  'Parent', aUIControls.mPanelLayout);
aUIControls.mSliderLayoutH = uicontrol('Style', 'slider', ...
  'Min', 0, 'Max', 1, 'Value', 0, 'SliderStep', [1, 1], 'Visible', 'off', ...
  'Parent', aUIControls.mPanelLayout, 'Callback', 'BPStitchDatasets(-10)');
aUIControls.mSliderLayoutV = uicontrol('Style', 'slider', ...
  'Min', 0, 'Max', 1, 'Value', 1, 'SliderStep', [1, 1], 'Visible', 'off', ...
  'Parent', aUIControls.mPanelLayout, 'Callback', 'BPStitchDatasets(-10)');

% source
aUIControls.mButtonOpen = uicontrol('Style', 'pushbutton', 'String', 'Open', ...
  'Parent', aUIControls.mPanelSource, 'Callback', 'BPStitchDatasets(-4)', ...
  'BackgroundColor', [1, 1, 1]);
aUIControls.mEditFileName = uicontrol('Style', 'edit', 'String', 'FileName', ...
  'HorizontalAlignment', 'Left', ...
  'Parent', aUIControls.mPanelSource, 'Callback', 'BPStitchDatasets(-5)');
aUIControls.mLabelLayout = uicontrol('Style', 'text', 'String', 'Layout from', ...
  'Parent', aUIControls.mPanelSource);
aUIControls.mRadioFileNames = uicontrol('Style', 'radio', 'String', 'File names', ...
  'Parent', aUIControls.mPanelSource, 'Value', 1, 'Callback', 'BPStitchDatasets(-6)');
aUIControls.mRadioFilesContent = uicontrol('Style', 'radio', 'String', 'Files content', ...
  'Parent', aUIControls.mPanelSource, 'Value', 0, 'Callback', 'BPStitchDatasets(-6)');
aUIControls.mCheckFlipY = uicontrol('Style', 'check', 'String', 'Flip Y', ...
  'Parent', aUIControls.mPanelSource, 'Callback', 'BPStitchDatasets(-7)');


%-------------------------------------------------------------------------%
%% Callback, rebuild table
%-------------------------------------------------------------------------%
function [aUIControls, aOptions] = ResizeWindow(aUIControls, aOptions)

vDelta = 5;
vSizeX = 50;
vSizeY = 20;

vWindow = get(aUIControls.mPanelStitch, 'Parent');
vWindowCoords = get(vWindow, 'Position');
vMove = max([vSizeX*6.5, vSizeY*12.5] - vWindowCoords(3:4), 0);
vWindowCoords = vWindowCoords + [0, -vMove(2), vMove];
set(vWindow, 'Position', vWindowCoords);

vDeltaX = vDelta / vWindowCoords(3);
vDeltaY = vDelta / vWindowCoords(4);
% vSizeXWin = vSizeX / vWindowCoords(3);
vSizeYWin = vSizeY / vWindowCoords(4);

% panels
vPosition = [vDeltaX, vDeltaY, 1-2*vDeltaX, 2*vDeltaY+3*vSizeYWin];
set(aUIControls.mPanelStitch, 'Position', vPosition);
vPosition = [vDeltaX, 4*vDeltaY+3*vSizeYWin, 1-2*vDeltaX, 1-8*vDeltaY-6*vSizeYWin];
set(aUIControls.mPanelLayout, 'Position', vPosition);
vPosition = [vDeltaX, 1-3*vDeltaY-3*vSizeYWin, 1-2*vDeltaX, 2*vDeltaY+3*vSizeYWin];
set(aUIControls.mPanelSource, 'Position', vPosition);

% stitching
vPosition = [vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mLabelShift, 'Position', vPosition);
vPosition = [vDelta + vSizeX + vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mEditShiftX, 'Position', vPosition);
vPosition = [vDelta + vSizeX*2 + vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mEditShiftY, 'Position', vPosition);
vPosition = [vDelta + vSizeX*3 + vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mEditShiftZ, 'Position', vPosition);
vPosition = [vDelta + vSizeX*4 + vDelta*2, vDelta, 0, vSizeY * 2 + vDelta];
vPosition(3) = vWindowCoords(3) - vPosition(1) - vDelta * 4;
set(aUIControls.mButtonStitch, 'Position', vPosition);
vPosition = [vDelta, vDelta + vSizeY + vDelta, vSizeX*4, vSizeY];
set(aUIControls.mCheckAdjust, 'Position', vPosition);

% layout
vPosition = [vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mLabelOverlap, 'Position', vPosition);
vPosition = [vDelta + vSizeX + vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mEditOverlapX, 'Position', vPosition);
vPosition = [vDelta + vSizeX*2 + vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mEditOverlapY, 'Position', vPosition);
vPosition = [vDelta + vSizeX*3 + vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mEditOverlapZ, 'Position', vPosition);
vPosition = [vDelta + vSizeX*4 + vDelta, vDelta, vSizeX, vSizeY];
set(aUIControls.mLabelUnits, 'Position', vPosition);

vPosition = get(aUIControls.mPanelLayout, 'Position');
vPosition = [vPosition(3)*vWindowCoords(3) - vDelta - vSizeY, vDelta + vSizeY, ...
  vSizeY, vPosition(4)*vWindowCoords(4) - vSizeY*2 - vDelta];
set(aUIControls.mSliderLayoutV, 'Position', vPosition);

vPosition = get(aUIControls.mPanelLayout, 'Position');
vPosition = [vDelta + vSizeX*4 + vDelta*2, vDelta, ...
  vPosition(3)*vWindowCoords(3) - vSizeX*4 - vDelta*4 - vSizeY, vSizeY];
set(aUIControls.mSliderLayoutH, 'Position', vPosition);

% source
vPosition = [vDelta, vDelta, vSizeX*1.5, vSizeY];
set(aUIControls.mLabelLayout, 'Position', vPosition);
vPosition = [vDelta + vSizeX*1.5 + vDelta, vDelta, vSizeX*1.5, vSizeY];
set(aUIControls.mRadioFileNames, 'Position', vPosition);
vPosition = [vDelta + vSizeX*3 + vDelta,  vDelta, vSizeX*2, vSizeY];
set(aUIControls.mRadioFilesContent, 'Position', vPosition);
vPosition = [vWindowCoords(3) - 4*vDelta - vSizeX, vDelta, vSizeX, vSizeY];
set(aUIControls.mCheckFlipY, 'Position', vPosition);
vPosition = [vDelta, vDelta + vSizeY + vDelta, vSizeX, vSizeY];
set(aUIControls.mButtonOpen, 'Position', vPosition);
vPosition = [vDelta + vSizeX + vDelta, vDelta + vSizeY + vDelta, 0, vSizeY];
vPosition(3) = vWindowCoords(3) - vPosition(1) - vDelta * 4;
set(aUIControls.mEditFileName, 'Position', vPosition);

% table
vDelta = 5;
vSizeX = 100;
for vIndex = 1:numel(aOptions.mFiles)
  vSizeX = max(vSizeX, 6*numel(aOptions.mFiles{vIndex}));
end
vSizeY = 20;
vSizeHeadX = 40;
vPosition = get(aUIControls.mSliderLayoutV, 'Position');
vMinX = vDelta;
vMaxY = vPosition(2) + vPosition(4) - vSizeY;

vPosition = get(aUIControls.mLabelOverlap, 'Position');
vPosition2 = get(aUIControls.mSliderLayoutH, 'Position');
aOptions.mGridX = fix((vPosition2(1)+vPosition2(3)-vPosition(1)-2*vSizeHeadX)/vSizeX)+2;
aOptions.mGridX = max(aOptions.mGridX, 3);
vPosition = get(aUIControls.mSliderLayoutV, 'Position');
aOptions.mGridY = fix(vPosition(4)/vSizeY);

for vIndex = numel(aUIControls.mPanelLayoutLabels)+1:aOptions.mGridX*aOptions.mGridY
  aUIControls.mPanelLayoutLabels(vIndex) = uicontrol('Style', 'text', ...
    'BackgroundColor', [1, 1, 1], ... 
    'Parent', aUIControls.mPanelLayout);
end

vX = vMinX;
vY = vMaxY;
for vIndex = 1:aOptions.mGridX*aOptions.mGridY
  vPosition = [vX, vY, vSizeX-1, vSizeY-1];
  set(aUIControls.mPanelLayoutLabels(vIndex), 'HorizontalAlignment', 'Center');
  if mod(vIndex, aOptions.mGridX) == 0
    vX = vMinX;
    vY = vY - vSizeY;
  elseif mod(vIndex, aOptions.mGridX) < 3
    vX = vX + vSizeHeadX;
    vPosition(3) = vSizeHeadX-1;
  else
    vX = vX + vSizeX;
  end
  if mod(vIndex-1, aOptions.mGridX) > 1 && vIndex > aOptions.mGridX
    set(aUIControls.mPanelLayoutLabels(vIndex), 'HorizontalAlignment', 'Left');
  end
  set(aUIControls.mPanelLayoutLabels(vIndex), 'String', sprintf('%d', vIndex));
  set(aUIControls.mPanelLayoutLabels(vIndex), 'Position', vPosition);
  set(aUIControls.mPanelLayoutLabels(vIndex), 'Visible', 'on');
end

for vIndex = aOptions.mGridX*aOptions.mGridY+1:numel(aUIControls.mPanelLayoutLabels)
  set(aUIControls.mPanelLayoutLabels(vIndex), 'Visible', 'off');
end


function fcorr = fourierCorrND( in1, in2 )
% fourierCorrND calculates ND cross-correlation in Fourier space
%
% SYNOPSIS   fcorr = fourierCorrND( in1, in2 )
%
% INPUT      in1  : input image 1
%            in2  : input image 2
%
% OUTPUT     fcor : ND correlation function, DC in the middle
%
% REMARK     The correlation is NOT NORMALIZED.
%
% Aaron Ponti, 2007/03/20

%--------------------------------------------------------------------------
%
%   The contents of this file are subject to the Mozilla Public License
%   Version 1.1 (the "License"); you may not use this file except in
%   compliance with the License. You may obtain a copy of the License at 
%   http://www.mozilla.org/MPL/
% 
%   Software distributed under the License is distributed on an "AS IS"
%   basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
%   License for the specific language governing rights and limitations 
%   under the License.
% 
%   The Original Code is "Qu for MATLAB".
% 
%   The Initial Developer of the Original Code is Aaron Ponti.
%   All Rights Reserved.
%
%--------------------------------------------------------------------------

f1    = fftn( single( in1 ) );
f2    = fftn( single( in2 ) );
fcorr = fftshift( ifftn( f2 .* conj( f1 ) ) );