% % % % MatlabXT::UConnHanDewarp(%i) % % % function UConnHanDewarp(aImarisApplicationID) % Get parameters from user prompt = {'Optical slices to keep per section','Slices per section (enter as array, or leave blank if only processing one section)'}; defaultans = {'7', ''}; answer = inputdlg(prompt, 'Parameters', 1, defaultans); outputSizeZ = str2double(answer{1}); % Number of optical slices to keep per section slicesPerSection = str2num(answer{2}); % Number of slices per section in the stack, or blank if only processing one section % Connect to Imaris 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 vDataSet = vImarisApplication.GetDataSet; % Preparation sizeX = vDataSet.GetSizeX; sizeY = vDataSet.GetSizeY; fullSizeZ = vDataSet.GetSizeZ; if isempty(slicesPerSection) slicesPerSection = fullSizeZ; else if sum(slicesPerSection) ~= fullSizeZ error('Array of slices per section does not add up to total number of slices'); end end sizeS = numel(slicesPerSection); sizeC = vDataSet.GetSizeC; sizeT = vDataSet.GetSizeT; newSizeX = sizeX; newSizeY = sizeY; % Giant loop over sections currentReadingZ = 1; currentWritingZ = 1; for s = 1 : sizeS sizeZ = slicesPerSection(s); multiChStack = uint16(zeros(sizeX, sizeY, sizeZ, sizeC)); for c = 1 : sizeC for z = 1 : sizeZ fprintf('Reading section %d/%d, channel %d/%d, z = %d\n', s, sizeS, c, sizeC, currentReadingZ + z - 1) multiChStack(:, :, z, c) = vDataSet.GetDataSliceShorts(currentReadingZ + z - 2, c - 1, 0); end end % Loop over channels slicesToKeep4D = false(sizeX, sizeY, sizeZ, sizeC); for c = 1 : sizeC % Reduce resolution for calculations lowResStack = zeros(newSizeX, newSizeY, sizeZ); for z = 1 : sizeZ lowResStack(:, :, z) = multiChStack(:, :, z, c); end % Calculate total brightnesses of all possible contiguous sets of z-slices sizePossBrsZ = sizeZ - outputSizeZ + 1; possBrs = zeros(newSizeX, newSizeY, sizePossBrsZ); for z = 1 : sizePossBrsZ possBrs(:, :, z) = sum(lowResStack(:, :, z:(z + outputSizeZ - 1)), 3); end % Determine which sets maximize brightness [maxBrs, slicesToKeep] = max(possBrs, [], 3); % Interpolate to full resolution slicesToKeep(slicesToKeep > sizePossBrsZ) = sizePossBrsZ; slicesToKeep(slicesToKeep < 1) = 1; % Fill slices-to-keep 4D matrix with logical values for z = 1 : sizeZ slicesToKeep4D(:, :, z, c) = (slicesToKeep > (z - outputSizeZ)) & (slicesToKeep < (z + 1)); end end % Dewarp slicesToKeep4D = permute(slicesToKeep4D, [3, 1, 2, 4]); % Note dim 3 is first multiChStack = permute(multiChStack, [3, 1, 2, 4]); stack2 = uint16(zeros(outputSizeZ, sizeX, sizeY, sizeC)); stack2(:) = multiChStack(slicesToKeep4D); stack2 = permute(stack2, [2, 3, 1, 4]); % Send back to Imaris for c = 1 : sizeC for z = 1 : outputSizeZ fprintf('Writing section %d/%d, channel %d/%d, z = %d\n', s, sizeS, c, sizeC, currentWritingZ + z - 1) vDataSet.SetDataSliceShorts(stack2(:, :, z, c), currentWritingZ + z - 2, c - 1, 0); end end currentReadingZ = currentReadingZ + slicesPerSection(s); currentWritingZ = currentWritingZ + outputSizeZ; end % Crop fprintf('Cropping\n') vDataSet.Crop(0, sizeX, 0, sizeY, 0, sizeS * outputSizeZ, 0, sizeC, 0, sizeT);