Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • scicomp/soans_et_al_opticalcup_surface
1 result
Show changes
Commits on Source (3)
Hide whitespace changes
Inline Side-by-side
Showing
with 2199 additions and 3 deletions
AnalyseOrientations/AnalyseOpticalCupLamininOrientation.m 0 → 100644
View file @ fc349606
% Analyse laminin orientations in projected data of zebrafish optical cup
% in data by Karen Soans, Caren Norden group, MPI-CBG
% Note that the input to this script are derotated and max-projected ROIs
% from the original image stack, obtained within Fiji by the script
% extractDerotatedSurface.py by Noreen Walker, MPI-CBG.
% Here, orientations are inferred from small image patches, and the
% alignment between orientations (or their coherency) is calculated locally
% and globally across each region of interest.
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% last edit: 2022_09_05
%
% developed using MATLAB R2019, R2020, and R2021
%
%%% usage %%%
% It is assumed that regions of interest (ROIs) of the original image stack
% were extracted, derotated by extractDerotatedSurface.py, and each saved
% into a folder structure
% head_folder/
% <Stage>/
% <Embryo>/
% <RegionOfTheEye>/
% <DerotationConditions>/
% Roi<nn>_ZMAX_of_Resliced_of_Derotated_with2ndCrop.tif
% [the pilot data had only <head_folder>/<Stage>/<Embryo>/<EyeAndChannel>/<ROI>/ ]
% Then per-ROI observables are collected in the array output_collected
% together with their relevant settings and conditions.
%
% Orientations are expressed as vector with coordinates (u,v) where
% anti-podal vectors (u,v) and (-u,-v) represent the same orientation,
% or as magnitude magn and angle theta (from the interval [0,pi]) with
% (u,v) = +- magn * (cos theta, sin theta).
% CHANGELOG
% 2022_08_23 clean-up for publication
% 2021_05_25 adapt plotting to downsampling
% 2021_05_18 include downsampling (i) of the input grayscale image and
% (ii) on the level of orientions.
% 2021_03_26 smoothing now discouraged - replaced by downsampling (on
% grayscale or on orientations)
% 2021_03_16 curate computation of theta_for_hists
% 2021_03_12 switch to orientation estimation by the gradient method (with
% gradient estimate by a Scharr operator),
% introduced kernelSize as opposed to boxRadius
% 2021_01_20 adapted to the structure of 2021 data
%% settings for orientation analysis
% Scharr operator with kernel size 5 is default and highly recommended for
% extremely low error in orientation estimation, see Tabelle B.11 in
%
% Hanno Scharr (2000): Optimale Operatoren in der digitalen
% Bildverarbeitung, PhD Thesis at University Heidelberg, Germany,
% DOI=10.11588/heidok.00000962, http://www.ub.uni-heidelberg.de/archiv/962
%
% Alternative orientation estimation based on the tensor of inertia (tensor
% of second moments of image intensities) and on the Fourier transform are
% included for completeness
method = 'gradient'; %one of "fourier", "inertia", "gradient"
kernelSize = 5; %only relevant for method="gradient" - allowed values are 5 (recommended), 3 and 2
maxImageIntensity = 256; %only relevant for method="gradient" - necessary for scaling of the gradient
boxRadius = 10; %only relevant for method="inertia" and method="fourier" - orientation estimation is performed in square boxes of side length 2*boxRadius+1
forCenterOfMass = true; %only relevant for method="inertia" - compute second moments for center of the image patch or center of mass (image intensity)
subtractBackgroud = true; %only relevant for method="inertia"
minWaveLength = 1; %only relevant for method="fourier" - wavelength to consider for radial integration in the Fourier domain
maxWaveLength = 10; %only relevant for method="fourier"
%% settings for optional downsampling of the grayscale image and/or the extracted orientation field
downSamplFactorGrayscale = 1; %needs to be integer, 1=no downsampling
downSamplFactorOrient = 1; %needs to be integer, 1=no downsampling
downSamplTotal = downSamplFactorGrayscale * downSamplFactorOrient;
%% collect prepocessing settings in a "methodInfo"
switch method
case 'fourier'
methodInfo = [ method '_radius' num2str(boxRadius) '_WaveLengths' num2str(minWaveLength) 'to' num2str(maxWaveLength) ];
case 'inertia'
methodInfo = [ method '_radius' num2str(boxRadius) '_' repmat('NOT', 1-forCenterOfMass) 'forCenterOfMass' '_' repmat('NOT', 1-subtractBackgroud) 'subtractBackgroud' ];
case 'gradient'
methodInfo = [ method '_kernel' num2str(kernelSize) ];
boxRadius = kernelSize; % use boxRadius as means to store the value of kernelSize, p.e. for usage in image erosion
end
methodInfo = [methodInfo '_downsample' num2str(downSamplFactorGrayscale) 'x' num2str(downSamplFactorOrient) ];
%% settings for analysis of coherency, that is degree of local alignment betweem orientations in a box
boxRadius_LocalAlignment = 10; %this quantity counts in ORIGINAL pixels, defining a box of side length 1+2*boxRadius_LocalAlignment
boxRadius_LocalAlignment_downSampled = floor( boxRadius_LocalAlignment / downSamplTotal );
powerOfMagn_LocalAlignment = 1;
myStrengthExponent = 1.0; % for smoothing
mySigmaForDirectySmoothing = [1.0 1.0]; %kernel radius for smoothing
requiredFracOfOrientInfo = 0.9; %fraction of orientations that must lay inside the region
% of interest (or its foreground) for a box to be included into the statistics
%% setting data folder and preparing output collection
head_folder = '/Users/khoffman/Documents/DirectionalityDefects/processedData/2021_08_06_KarenSoans/';
genericFilename = 'ZMAX_of_Resliced_of_Derotated_with2ndCrop.tif'; % according to the output of extractDerotatedSurface.py,
% specify WITHOUT leading individual filename portions like 'ROI<number>_'
filesToAnalyse = dir([head_folder '/**/Roi*_' genericFilename ]);
% the levels below head_folder ougth to be <Stage>/<Embryo>/<RegionOfTheEye>/<DerotationConditions>/
% whereby <RegionOfTheEye> is one of the ROInames = ["hollow", "inBetween", "ridge", "rim", "back"];
% prepare matrix/array for collection of statistics
output_folder_stats = [head_folder 'collectedResults' '_' methodInfo ...
'_coherencyRadius' num2str(boxRadius_LocalAlignment) 'Power' num2str(powerOfMagn_LocalAlignment) '/'];
temp = strsplit(output_folder_stats, '/');
% within the folder, specify file to hold the statistics (its name = name of the deepest subfolder)
output_collected = [output_folder_stats temp{end-1} '.csv'];
%folder for the statistics across ROIs and plots
if ~exist(output_folder_stats, 'dir')
mkdir( output_folder_stats )
end
%% analysis file by file (ROI by ROI)
all_vals_mat = [];
folders_arr = {};
for numFile = 1:numel(filesToAnalyse)
input_folder = [filesToAnalyse(numFile).folder '/' ] ;
output_folder= input_folder;
%% load the grayscale image
input_fullPath = [input_folder filesToAnalyse(numFile).name] ;
[~, input_soleFilename, ~] = fileparts(input_fullPath);
leading_ROIinfo_inFilename = input_soleFilename( 1 : strfind(input_soleFilename, genericFilename(1:end-4))-1 ); %find name part that is individual to this file
[~, GrayImageStack] = loadGrayStack(input_fullPath);
assert(numel(size(GrayImageStack)) ==2) %2D data
NumSlices = 1;
% scale image intensities to a fixed range for comparable results by the gradient
if max(GrayImageStack(:))>255
warning('image intensities are being scaled down to <=255')
GrayImageStack = GrayImageStack / max(GrayImageStack(:)) ;
end
% optional: show loaded data
% imshow(GrayImageStack_orig / 255); title('original data')
GrayImageStack_orig = GrayImageStack;
%% optional downsampling of the grayscale image (downSamplFactorGrayscale = 1 means no downsampling)
% NOTE: Identification of the foreground still happens in the original image,
% but orientation analysis is performed in the downsampled grayscale data
size_y_orig = size(GrayImageStack_orig,1);
size_x_orig = size(GrayImageStack_orig,2);
assert( int8(downSamplFactorGrayscale) == downSamplFactorGrayscale, "must be integer")
GrayImageStack_blocks = 1/(downSamplFactorGrayscale^2) ...
* conv2(GrayImageStack, ones(downSamplFactorGrayscale), 'valid') ;
GrayImageStack_downSampl = GrayImageStack_blocks(1:downSamplFactorGrayscale:end, 1:downSamplFactorGrayscale:end);
GrayImageStack = GrayImageStack_downSampl ;
size_y_grayBlocks = size(GrayImageStack_downSampl,1);
size_x_grayBlocks = size(GrayImageStack_downSampl,2);
%% load info on preprocessing
% in particular rotations performed to check for anisotropy along the z-axis in the original imaging system
% that might result from the form of the point-spread function
input_preprocessingInfoPath = [input_folder leading_ROIinfo_inFilename 'logFileParameters.txt'];
fileID = fopen(input_preprocessingInfoPath);
preprocessingInfo_raw = fread(fileID,'*char')';
fclose(fileID);
% find text location where normal vector is represented
pattern = 'Smallest eigen vector (V0) = normal vector of plane: (nx,ny,nz)= ';
textlocation = strfind(preprocessingInfo_raw, pattern);
if isempty(textlocation)
warning('no info on normal vector found')
NormalVecOrig = [NaN, NaN, NaN];
else
NormalVecOrig_raw = strsplit( preprocessingInfo_raw(textlocation+length(pattern) : end)); %split at whitespace and newline
NormalVecOrig_raw = NormalVecOrig_raw(1:3);
NormalVecOrig = str2double(NormalVecOrig_raw);
assert( abs( 1- sum(NormalVecOrig.^2) ) < 1e-5, 'Unit length violated by normal vector of the plane that was read')
end
%NormalVecOrig is the normal vector of the plane in the original data,
% where the z-axis is the axis of anisotropy
Rotation_theta_x = atan2(NormalVecOrig(2) , NormalVecOrig(3)); % equals atan(NormalVecOrig(2) / NormalVecOrig(3))
Rotation_theta_y = - atan2(NormalVecOrig(1) , sqrt(1-NormalVecOrig(1)^2) ); % named 'theta_yprime' in the preprocessingInfo_raw
AnisotropyAxis = [ sin(Rotation_theta_y), ...
-cos(Rotation_theta_y) * sin(Rotation_theta_x), ...
cos(Rotation_theta_y) * cos(Rotation_theta_x) ] ;
%AnisotropyAxis is the axis of anisotropy in the rotated data,
% where the z-axis is the normal axis of the plane (and this dimension is omitted by projection)
%AnisotropyAxis is in the left-handed coordinate systems of Fiji (and MATLAB)
% therefore distinctino between "MATLABcoords" and righthanded "NORMALcoords"
AnisotropyAxis_inThePlane_MATLABcoords = AnisotropyAxis(1:2) / sqrt(sum( AnisotropyAxis(1:2).^2)) ;
%for the angle, take the same approach as for orientations theta: use the standard right-handed coordinate layout
AngleAnisotropy_inThePlane_NORMALcoords = atan2( -AnisotropyAxis(2), AnisotropyAxis(1) );
AngleAnisotropyProjection = acos(NormalVecOrig(3)) ; % or = acos(AnisotropyAxis(3));
% in 3D: angle between axis of anisotropy and direction of projection
%% orientation analysis
%% Prepare output arrays
theta = zeros( size(GrayImageStack) );
magn = zeros( size(GrayImageStack) );
%% Orientation analysis
time_needed_total=0;
for useSlice = 1:NumSlices
tic
switch method
case 'inertia'
[directy, ~] = ComputeDirecty(GrayImageStack(:,:,useSlice), boxRadius, forCenterOfMass, subtractBackgroud) ;
case 'fourier'
[directy, ~] = ComputeDirecty_FourierAndQ(GrayImageStack(:,:,useSlice), boxRadius, minWaveLength, maxWaveLength);
case 'gradient'
directy = ComputeDirectyGradient(GrayImageStack(:,:,useSlice), kernelSize );
otherwise
warning(['method "' method '" not known'])
end
theta(:,:,useSlice) = atan2( directy(:,:,2), directy(:,:,1) );
if strcmp(method,'gradient')
magn (:,:,useSlice) = directy(:,:,3) / maxImageIntensity;
else %'inertia', 'fourier'
magn (:,:,useSlice) = directy(:,:,3);
end
time_needed=toc;
time_needed_total = time_needed_total + time_needed;
disp(['... done in ', num2str(time_needed), ' seconds.'])
end
disp(['All slices done in ', num2str(time_needed_total), ' seconds.'])
%% save result arrays (together with their settings) to file
% as checkpoint, because orientation estimation with method='inertia' and 'fourier' can take up to hours
filename_intermediate_output = [input_fullPath(1:end-4) '_' methodInfo '.mat'] ;
save( filename_intermediate_output, ...
'GrayImageStack_orig', ... %the original data
'input_fullPath', ... % its source
'NormalVecOrig', ... % the normal vector of the plane in original data, from which the applied rotation follows
'method', 'boxRadius', 'forCenterOfMass', 'subtractBackgroud', ...
'minWaveLength', 'maxWaveLength', ... %the settings for orientation analysis
'theta', 'magn', ... %the results
'downSamplFactorGrayscale', 'downSamplFactorOrient' ) %information on downsampling
%load(filename_intermediate_output)
%% find foreground = image parts sufficiently off the boundary
% (determined without downsampling in the original image data of highest resolution)
% note that the actual data is surrounded by a dim region of polygonal
% shape due to the repeated cutting and projecting in Fiji
% fg_closed = approximate the foreground region (after morphological closing for smoother shape and convex hull operation
% as the shapes du to cutting/projecting are known to be convex polygons)
% fg_inner = all those pixels, for which the complete box used for orientation estimation was on the foreground fg_closed
% fg_innermost = all those pixels, for which all those pixels taken into
% account for local alignment had their respective box used for orientation estimation on the foreground fg_closed
foreground_orig = imgaussfilt(GrayImageStack_orig/255, 5, 'FilterDomain','spatial', 'Padding', 0.0 );
foreground_orig = imbinarize(foreground_orig/max(foreground_orig(:)), 'global'); %imbinarize expects input normalised in [0,1]
% imshow(foreground)
structureElement = strel('diamond',5);
fg_closed_orig = imclose(foreground_orig, structureElement);
%fg_closed_orig = imfill(fg_closed_orig, 'holes');
fg_closed_orig = bwconvhull(fg_closed_orig);
% compare
% figure(1); imshow(fg_closed_orig)
% figure(2); imshow(GrayImageStack_orig/255+0.5)
fg_inner_orig = imerode(fg_closed_orig, strel('square', 2*boxRadius+1) );
% figure(3); imshow(fg_inner_orig)
fg_innermost_orig = imerode(fg_inner_orig, strel('square', 2*boxRadius_LocalAlignment+1) );
% figure(4); imshow(fg_innermost_orig)
%% downsample the foreground information according to the downsampling of grayscale and orientation blocks
% two variants: count a block in the downsampled image as forground, when
% (A) any of the contained pixels is foreground, or when
% (B) all of the contained pixels are foreground (" .."strict")
%% foreground downsampled to the resoution of grayscale downsampling
fg_closed_blocks = 1/(downSamplFactorGrayscale^2) ...
* conv2(fg_closed_orig, ones(downSamplFactorGrayscale), 'valid') ;
fg_closed = (0 < fg_closed_blocks(1:downSamplFactorGrayscale:end, 1:downSamplFactorGrayscale:end) );
fg_closed_strict = (1 == fg_closed_blocks(1:downSamplFactorGrayscale:end, 1:downSamplFactorGrayscale:end) );
fg_inner_blocks = 1/(downSamplFactorGrayscale^2) ...
* conv2(fg_inner_orig, ones(downSamplFactorGrayscale), 'valid') ;
fg_inner = (0 < fg_inner_blocks(1:downSamplFactorGrayscale:end, 1:downSamplFactorGrayscale:end) );
fg_inner_strict = (1 == fg_inner_blocks(1:downSamplFactorGrayscale:end, 1:downSamplFactorGrayscale:end) );
fg_innermost_blocks = 1/(downSamplFactorGrayscale^2) ...
* conv2(fg_innermost_orig, ones(downSamplFactorGrayscale), 'valid') ;
fg_innermost = (0 < fg_inner_blocks(1:downSamplFactorGrayscale:end, 1:downSamplFactorGrayscale:end) );
fg_innermost_strict = (1 == fg_innermost_blocks(1:downSamplFactorGrayscale:end, 1:downSamplFactorGrayscale:end) );
%% foreground downsampled to the (coarser) resolution which the orientation downsampling will attain
% downsample the full-resolution foreground information like "fg_closed_orig"
% by the factor downSamplTotal = downSamplFactorGrayscale * downSamplFactorOrient;
fg_closed_blocks = 1/(downSamplTotal^2) ...
* conv2(fg_closed_orig, ones(downSamplTotal), 'valid') ;
fg_closed_downSampl = (0 < fg_closed_blocks(1:downSamplTotal:end, 1:downSamplTotal:end) );
fg_closed_downSampl_strict = (1 == fg_closed_blocks(1:downSamplTotal:end, 1:downSamplTotal:end) );
fg_inner_blocks = 1/(downSamplTotal^2) ...
* conv2(fg_inner_orig, ones(downSamplTotal), 'valid') ;
fg_inner_downSampl = (0 < fg_inner_blocks(1:downSamplTotal:end, 1:downSamplTotal:end) );
fg_inner_downSampl_strict = (1 == fg_inner_blocks(1:downSamplTotal:end, 1:downSamplTotal:end) );
fg_innermost_blocks = 1/(downSamplTotal^2) ...
* conv2(fg_inner_orig, ones(downSamplTotal), 'valid') ;
fg_innermost_downSampl = (0 < fg_innermost_blocks(1:downSamplTotal:end, 1:downSamplTotal:end) );
fg_innermost_downSampl_strict = (1 == fg_innermost_blocks(1:downSamplTotal:end, 1:downSamplTotal:end) );
%% smooth orientations with Gaussian kernel, and find versions restricted to the foreground
[theta_smoothed, magn_smoothed] = smoothDirecty_ThetaStrength(theta, magn, mySigmaForDirectySmoothing, myStrengthExponent);
% alternative versions restrict to points in the foreground (_fg)
% and those for which the image patch used for orientation estimation was fully within the foreground (_inner)
% and those for which all pixels considered for local alignment had their respective box for orientation estimation fully in the foreground (_innermost)
% (the last is NOT the same as "all pixels considered in the smoothing step")
magn_smoothed_fg = NaN (size(magn_smoothed));
magn_smoothed_fg(fg_closed) = magn_smoothed(fg_closed);
magn_smoothed_fg_inner = NaN (size(magn_smoothed));
magn_smoothed_fg_inner(fg_inner) = magn_smoothed(fg_inner);
magn_smoothed_fg_innermost = NaN (size(magn_smoothed));
magn_smoothed_fg_innermost(fg_innermost) = magn_smoothed(fg_innermost);
%% optional: use the smoothed orientation field for further analysis
% basically set theta = theta_smoothed but assert values in the interval [0, pi]
% NOTE: smoothing might interfere with the subsequent analysis of alignment
%{
theta = mod(theta_smoothed, pi) ;
magn = magn_smoothed;
%}
%% optional: downsampling of the orientation field (downSamplFactorOrient = 1 means no downsampling)
assert( int8(downSamplFactorOrient) == downSamplFactorOrient, "must be integer")
[theta_downSampl, magn_downSampl] = downSamplDirecty_ThetaStrength(theta, magn, downSamplFactorOrient); % skipRows=0 and skipCols=0 fits to the grayscale downsampling in this file, strength_exponent)
size_y_grayAndOrientBlocks = size(theta_downSampl,1);
size_x_grayAndOrientBlocks = size(theta_downSampl,2);
% also downsample the grayscale intensities (again) as they will be used as weights in the check for anisotropy
GrayImageStack_blocksCoarse = 1/(downSamplFactorOrient^2) ...
* conv2(GrayImageStack_downSampl, ones(downSamplFactorOrient), 'valid') ;
GrayImageStack_downSamplTwice = GrayImageStack_blocksCoarse(1:downSamplFactorOrient:end, 1:downSamplFactorOrient:end);
assert( all([size_y_grayAndOrientBlocks, size_x_grayAndOrientBlocks] == size(GrayImageStack_downSamplTwice) ) )
% use downsampled orientation field for further analysis
theta_orig = theta;
magn_orig = magn ;
% basically set theta = theta_downSampl but assert values in the interval [0, pi]
theta = mod(theta_downSampl, pi) ;
magn = magn_downSampl;
%% alignment as a function of space
%% compute local alignment in fixed box - potentially time consuming
% use boxRadius_LocalAlignment_downSampled that is adjusted to the downsampling performed beforehand
[localOrient, localCoherency] = ComputeLocalAlignment(theta,magn,boxRadius_LocalAlignment_downSampled, powerOfMagn_LocalAlignment);
%% save result arrays of local alignment (together with their settings) to file
% as checkpoint
filename_intermediate_output1 = [filename_intermediate_output(1:end-4) ...
'_coherencyRadius' num2str(boxRadius_LocalAlignment) 'Power' num2str(powerOfMagn_LocalAlignment) '.mat'];
save( filename_intermediate_output1, ...
'GrayImageStack_orig', ... %the original data
'input_fullPath', ... % its source
'NormalVecOrig', ... % the normal vector of the plane in original data, from which the applied rotation follows
'method', 'boxRadius', 'forCenterOfMass', 'subtractBackgroud', ...
'minWaveLength', 'maxWaveLength', ... %the settings for orientation analysis
'theta', 'magn', ... %the results
'downSamplFactorGrayscale', 'downSamplFactorOrient', ... %information on downsampling
'fg_closed', 'fg_inner', 'fg_innermost', ... % different definitions of foreground
'boxRadius_LocalAlignment', 'boxRadius_LocalAlignment_downSampled', 'powerOfMagn_LocalAlignment', ... %the settings for local alignment analysis
'localOrient', 'localCoherency' ) %the results
%load(filename_intermediate_output1)
%% local alignment restricted to foreground versions
% recall: this restricts to points in the foreground (_fg)
% and those for which the image patch used for orientation estimation was fully within the foreground (_inner)
% and those for which all pixels considered for local alignment had their respective box for orientation estimation fully in the foreground (_innermost)
localCoherency_fg = NaN (size(localCoherency));
localCoherency_fg(fg_closed_downSampl) = localCoherency(fg_closed_downSampl);
localCoherency_fg_inner = NaN (size(localCoherency));
localCoherency_fg_inner(fg_inner_downSampl) = localCoherency(fg_inner_downSampl);
localCoherency_fg_innermost = NaN (size(localCoherency));
localCoherency_fg_innermost(fg_innermost_downSampl) = localCoherency(fg_innermost_downSampl);
%% compute global alignment of ALL orientations present
% in the whole image, in the foreground (_fg), and of those for which the
% image patch used for orientation estimation was fully within the foreground (_inner)
[overallOrient, overallCoherency] = ComputeQTensor(theta, magn, powerOfMagn_LocalAlignment);
[overallOrient_fg, overallCoherency_fg] = ComputeQTensor(theta, magn.*fg_closed_downSampl, powerOfMagn_LocalAlignment);
[overallOrient_fg_inner, overallCoherency_fg_inner] = ComputeQTensor(theta, magn.*fg_inner_downSampl , powerOfMagn_LocalAlignment);
%% Alternatively, limit further analysis of local orientation to points
% for which the neighbourhood considered contains not too much pixels without orientation information
OrientInfoPresent = (~isnan(theta))&(~isnan(magn));
OrientInfoPresent(1:boxRadius, :) = 0;
OrientInfoPresent(end-boxRadius+1:end, :) = 0;
OrientInfoPresent(:, 1:boxRadius ) = 0;
OrientInfoPresent(:, end-boxRadius+1:end) = 0;
LocalOrientIsReliable = requiredFracOfOrientInfo * conv2( ones(size(theta)) , ones( 2*boxRadius_LocalAlignment_downSampled+1), 'same' ) ...
<= conv2( OrientInfoPresent , ones( 2*boxRadius_LocalAlignment_downSampled+1), 'same' ) ;
LocalOrientIsReliably_fg = requiredFracOfOrientInfo * conv2( ones(size(theta)) , ones( 2*boxRadius_LocalAlignment_downSampled+1), 'same' ) ...
<= conv2( OrientInfoPresent.*fg_closed_downSampl , ones( 2*boxRadius_LocalAlignment_downSampled+1), 'same' ) ;
[overallOrient_reliable, overallCoherency_reliable] = ComputeQTensor(theta, magn.*LocalOrientIsReliable, powerOfMagn_LocalAlignment);
[overallOrient_reliably_fg, overallCoherency_reliably_fg] = ComputeQTensor(theta, magn.*LocalOrientIsReliably_fg, powerOfMagn_LocalAlignment);
%% compute alignment of orientations with axis of anisotropy in the data
% alignment1 = linear = absolute value of scalar product between two vectors
% alignment2 = quadratic like the nematic order parameter, but with the (projected) axis of anisotropy as reference
% s = average over (3/2 cos^2 beta - 1/2)
% where beta = angle between orientation and (projected) axis of anisotropy, hence
% cos beta = scalar product of orientation with (projected) axis of anisotropy (when both are expressed as normalised vectors)
% values of s range from -0.5 (all detected orientations perpendicular to axis of anisotropy)
% to +1.0 (all detected orientations along the axis of anisotropy), but
% should be around zero (indicating little bias of the orientations by the anisotropy of imaging in the original z-axis)
alignment1 = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) )) / nansum(nansum( isfinite(theta) )) ;
alignment1_fg = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* fg_closed_downSampl )) ;
alignment1_fg_inner = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* fg_inner_downSampl )) ;
alignment2 = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 )) / nansum(nansum( isfinite(theta) )) ;
alignment2_fg = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* fg_closed_downSampl )) ;
alignment2_fg_inner = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* fg_inner_downSampl )) ;
% versions that are weighted by magnitude of the orientation information and/or image intensity at that pixel
alignment1_magn = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* magn )) / nansum(nansum( isfinite(theta) .* magn )) ;
alignment1_magn_fg = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* magn .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* fg_closed_downSampl )) ;
alignment1_magn_fg_inner = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* magn .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* fg_inner_downSampl )) ;
alignment2_magn = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* magn )) / nansum(nansum( isfinite(theta) .* magn )) ;
alignment2_magn_fg = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* magn .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* fg_closed_downSampl )) ;
alignment2_magn_fg_inner = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* magn .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* fg_inner_downSampl )) ;
alignment1_intens = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* GrayImageStack_downSamplTwice )) / nansum(nansum( isfinite(theta) .* GrayImageStack_downSamplTwice )) ;
alignment1_intens_fg = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) ;
alignment1_intens_fg_inner = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) ;
alignment2_intens = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* GrayImageStack_downSamplTwice )) / nansum(nansum( isfinite(theta) .* GrayImageStack_downSamplTwice )) ;
alignment2_intens_fg = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) ;
alignment2_intens_fg_inner = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) ;
alignment1_magn_intens = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* magn .* GrayImageStack_downSamplTwice )) / nansum(nansum( isfinite(theta) .* magn .* GrayImageStack_downSamplTwice )) ;
alignment1_magn_intens_fg = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* magn .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) ;
alignment1_magn_intens_fg_inner = -2/pi+ nansum(nansum( abs( cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ) .* magn .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) ;
alignment2_magn_intens = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* magn .* GrayImageStack_downSamplTwice )) / nansum(nansum( isfinite(theta) .* magn .* GrayImageStack_downSamplTwice )) ;
alignment2_magn_intens_fg = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* magn .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* GrayImageStack_downSamplTwice .* fg_closed_downSampl )) ;
alignment2_magn_intens_fg_inner = -0.5 + nansum(nansum( (cos(theta).* cos(AngleAnisotropy_inThePlane_NORMALcoords) + sin(theta) .* sin(AngleAnisotropy_inThePlane_NORMALcoords) ).^2 .* magn .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) / nansum(nansum( isfinite(theta) .* magn .* GrayImageStack_downSamplTwice .* fg_inner_downSampl )) ;
%% collect results in structs with descriptions and values
%start with info on size of foreground
vals = struct();
descr = struct();
vals. size_x = size(theta,2);
descr.size_x = 'image size x';
vals. size_y = size(theta,1);
descr.size_y = 'image size y';
vals .pixels = sum(~isnan(theta(:)));
descr.pixels = 'pixels with orientation';
assert(all(unique(fg_closed) == [ 0; 1]))
vals .pixels_fg = sum(fg_closed(:));
descr.pixels_fg = 'foreground pixels';
assert(all(unique(fg_inner) == [ 0; 1]))
vals. pixels_fg_inner = sum(fg_inner(:));
descr.pixels_fg_inner = 'pixels for wich box was within foreground';
%% add results of coherencies in structs with descriptions and values
vals. overallCoherency = overallCoherency;
descr.overallCoherency = 'coherency - all orientations';
vals. overallTheta = atan2(overallOrient(2), overallOrient(1));
descr.overallTheta = 'theta - all orientations';
vals. overallCoherency_fg = overallCoherency_fg;
descr.overallCoherency_fg = 'coherency - all orientations on foreground';
vals. overallTheta_fg = atan2(overallOrient_fg(2), overallOrient_fg(1));
descr.overallTheta_fg = 'theta - all orientations on foreground';
vals. overallCoherency_fg_inner = overallCoherency_fg_inner;
descr.overallCoherency_fg_inner = 'coherency - all orientations for which box was within foreground';
vals. overallTheta_fg_inner = atan2(overallOrient_fg_inner(2), overallOrient_fg_inner(1));
descr.overallTheta_fg_inner = 'theta - all orientations for which box was within foreground';
vals. overallCoherency_reliab = overallCoherency_reliable;
descr.overallCoherency_reliab = 'coherency - all orientations for which box had mostly orientation information';
vals. overallTheta_reliab = atan2(overallOrient_reliable(2), overallOrient_reliable(1));
descr.overallTheta_reliab= 'theta - all orientations for which box box had mostly orientation information';
vals. overallCoherency_reliab_fg = overallCoherency_reliably_fg;
descr.overallCoherency_reliab_fg = 'coherency - all orientations for which box had mostly foreground and orientation information';
vals. overallTheta_reliab_fg = atan2(overallOrient_reliably_fg(2), overallOrient_reliably_fg(1));
descr.overallTheta_reliab_fg = 'theta - all orientations for which box had mostly foreground and orientation information';
vals. localCoherency = nanmean(localCoherency(:));
descr.localCoherency = 'mean local coherency (in boxes)';
vals. localCoherency_fg = nanmean(localCoherency_fg(:));
descr.localCoherency_fg = 'mean local coherency (in boxes) over foreground';
vals. localCoherency_fg_inner = nanmean(localCoherency_fg_inner(:));
descr.localCoherency_fg_inner = 'mean local coherency (in boxes) where box was within foreground';
vals. localCoherencyGaussian = nanmean(magn_smoothed(:));
descr.localCoherencyGaussian = 'mean local coherency (Gaussian)';
vals. localCoherencyGaussian_fg = nanmean(magn_smoothed_fg(:));
descr.localCoherencyGaussian_fg = 'mean local coherency (Gaussian) over foreground';
vals. localCoherencyGaussian_fg_inner = nanmean(magn_smoothed_fg_inner(:));
descr.localCoherencyGaussian_fg_inner = 'mean local coherency (Gaussian) where box was within foreground';
%% add info on anisotropy of data
vals. normalVec_x = NormalVecOrig(1);
descr.normalVec_x = 'normal vector in original data n_x';
vals. normalVec_y = NormalVecOrig(2);
descr.normalVec_y = 'normal vector in original data n_y';
vals. normalVec_z = NormalVecOrig(3);
descr.normalVec_z = 'normal vector in original data n_z';
vals. AngleAnisotropyProjection = AngleAnisotropyProjection;
descr.AngleAnisotropyProjection = 'angle between projection axis and anisotropy (radians)';
vals. orientationBias_linear = alignment1_fg_inner;
descr.orientationBias_linear = 'bias of orientations by anisotropy - linear, where box was within foreground';
vals. orientationBias_quadratic = alignment2_fg_inner;
descr.orientationBias_quadratic = 'bias of orientations by anisotropy - quadratic (as for nematic order parameter), where box was within foreground';
%% prepare writing of data
vals_mat = struct2array(vals);
all_vals_mat = [all_vals_mat; vals_mat] ;
folder_to_subfolders = strsplit(input_folder, filesep);
folders_arr = [ folders_arr ; ...
{input_folder} folder_to_subfolders(end-4), folder_to_subfolders(end-3), folder_to_subfolders(end-2), folder_to_subfolders(end-1), ...
leading_ROIinfo_inFilename];
%recall: the code assumes that folder_to_subfolders
% contains <Stage>/<Embryo>/<RegionOfTheEye>/<DerotationConditions>/
% and the different ROIs within one folder are distinguished by leading_ROIinfo_inFilename
%% write results to text file
filename_intermediate_output2 = [filename_intermediate_output(1:end-4) ...
'_coherencyRadius' num2str(boxRadius_LocalAlignment) 'Power' num2str(powerOfMagn_LocalAlignment) '.txt'];
dlmwrite( filename_intermediate_output2 , ['coherencies obtained from'], 'delimiter','') %header
dlmwrite( filename_intermediate_output2 , [filename_intermediate_output2] , '-append', 'delimiter','') % origin of data
dlmwrite( filename_intermediate_output2 , ['Parameters used'] , '-append', 'delimiter','')
dlmwrite( filename_intermediate_output2 , [num2str(boxRadius) , ' BoxRadius'], '-append', 'delimiter','') % parameters used
dlmwrite( filename_intermediate_output2 , [num2str(boxRadius_LocalAlignment) , ' BoxRadius for Local Alignment'], '-append', 'delimiter','') % parameters used
dlmwrite( filename_intermediate_output2 , [num2str(boxRadius_LocalAlignment_downSampled) , ' BoxRadius for Local Alignment in downsampled pixels'], '-append', 'delimiter','') % parameters used
dlmwrite( filename_intermediate_output2 , [num2str(powerOfMagn_LocalAlignment), ' Power of Magnitude for Local Alignment'], '-append', 'delimiter','') % parameters used
dlmwrite( filename_intermediate_output2 , [num2str(mySigmaForDirectySmoothing) , ' sigma for smoothing'], '-append', 'delimiter','') % parameters used
dlmwrite( filename_intermediate_output2 , [num2str(myStrengthExponent), ' Power of Magnitude for smoothing'], '-append', 'delimiter','') % parameters used
dlmwrite( filename_intermediate_output2 , [' '], '-append', 'delimiter','') % empty line
fieldnameslist = fieldnames(vals);
for kk = 1:numel(fieldnameslist)
dlmwrite( filename_intermediate_output2 , [num2str(vals.(fieldnameslist{kk})), ' ', descr.(fieldnameslist{kk})], '-append', 'delimiter','')
end
%% plotting
%% Color and other conventions for plotting
col_director = [0.4 0.4 1.0]; %NOT 'blue', as this is too dark
col_localOrient = [1.0 0.4 0.4];
show_localOrient = false; % show locally averaged orientation that corresponds to the computed alignment
col_anisotropy = [0.2 1 0.2] ;
show_anisotropy = true;
centerQuiverOnPixels = true; %show orientations centered over their location
savePlots = false;
useMagnForPlotting = false; % scale quiver plots by orientation magnitude
%spacing12 between plotted orientation [integer, =n plots only every n-th orientation]
spacing12 = ceil(4.0/downSamplTotal); % a reasonable start for ROIs of few hundred pixels
%spacing12 = 1
% {
%% plotting of orig image with orientation field
myData = GrayImageStack_orig(:,:,useSlice);
myData = myData/ max(GrayImageStack_orig(:)); %rescale for plotting
myData = 0.2+0.8* myData; %make the data more bright in plotting
assert(max(myData(:)) <= 1.0, 'Rescale data!')
theta_slice = theta(:,:,useSlice);
magn_slice = GrayImageStack_downSamplTwice .* magn (:,:,useSlice);
[x12, y12] = meshgrid( 0.5+downSamplTotal*(-0.5+(1:size(theta_slice,2)) ), ...
0.5+downSamplTotal*(-0.5+(1:size(theta_slice,1)) ) );
% rescale positions that are given as row/column in the arrays of data
% after downsampling twice, back to their positions in the original data
scaleToOrig_x = @(xInDownSamplTwice) 0.5+(xInDownSamplTwice-0.5)*downSamplTotal;
scaleToOrig_y = @(yInDownSamplTwice) 0.5+(yInDownSamplTwice-0.5)*downSamplTotal;
u12 = magn_slice.^useMagnForPlotting .* cos( theta_slice ) ;
v12 = magn_slice.^useMagnForPlotting .* -sin( theta_slice ) ; %negative sign to transform from theta in normal coord layout to coord layout used by MATLAB (y downwards in plotting that is started by imshow() )
u_localOrient = localOrient(:,:,1) .* localCoherency;
v_localOrient = localOrient(:,:,2) .* localCoherency;
%% plotting itself
Fig12 = figure('PaperOrientation', 'landscape', 'Position', [0, 0, 1200, 700]);
[axesHandles12, positions12] = tight_subplot(1,1,[.08 .05],[.05 .01],[.05 .02]) ; %gap, margin_height, margin_width)
fig12A = axesHandles12(1);
axes(fig12A) %make the existing axes current
OrigData12 = imshow(myData); % (myData * 2 - 3);
hold on
% add overlay to indicate the foreground - see https://de.mathworks.com/company/newsletters/articles/image-overlay-using-transparency.html
% show true-color all-red image on top
red = cat(3, ones(size(myData)), zeros(size(myData)), zeros(size(myData)));
foreground_plot = imshow(red);
% but set its transparency according to fg_closed_orig
set(foreground_plot, 'AlphaData', 0.1* fg_closed_orig)
% orientations plotted as head-less vectors
quiv_onOrigData1 = quiver(x12(1:spacing12:end,1:spacing12:end), y12(1:spacing12:end,1:spacing12:end), u12(1:spacing12:end,1:spacing12:end), v12(1:spacing12:end,1:spacing12:end));
if centerQuiverOnPixels
%duplicate quiver (with opposite sign) to have coherent appearance
quiv_onOrigData2 = quiver(x12(1:spacing12:end,1:spacing12:end), y12(1:spacing12:end,1:spacing12:end), -u12(1:spacing12:end,1:spacing12:end), -v12(1:spacing12:end,1:spacing12:end));
end
% modify quiver plot(s)
quiv_onOrigData1.ShowArrowHead = 'off';
quiv_onOrigData1.LineWidth = 1.5; %(default: 0.5)
quiv_onOrigData1.AutoScaleFactor = 1.0 ;
quiv_onOrigData1.Color = col_director;
if centerQuiverOnPixels
%give the duplicate quiver the same properties, especially the same color
quiv_onOrigData2.ShowArrowHead = quiv_onOrigData1.ShowArrowHead ;
quiv_onOrigData2.LineWidth = quiv_onOrigData1.LineWidth ;
quiv_onOrigData2.AutoScaleFactor = quiv_onOrigData1.AutoScaleFactor ;
quiv_onOrigData2.Color = quiv_onOrigData1.Color ;
end
if show_localOrient
quiv_localOrient1 = quiver(x12(1:spacing12:end,1:spacing12:end), y12(1:spacing12:end,1:spacing12:end), u_localOrient(1:spacing12:end,1:spacing12:end), -v_localOrient(1:spacing12:end,1:spacing12:end));
%Note negative sign on y-component to transform from theta in normal coord layout to coord layout used by MATLAB (y downwards in plotting that is started by imshow() )
if centerQuiverOnPixels
%duplicate quiver (with opposite sign) to have coherent appearance
quiv_localOrient2 = quiver(x12(1:spacing12:end,1:spacing12:end), y12(1:spacing12:end,1:spacing12:end), -u_localOrient(1:spacing12:end,1:spacing12:end), v_localOrient(1:spacing12:end,1:spacing12:end));
end
% modify quiver plot(s)
quiv_localOrient1.ShowArrowHead = 'off';
quiv_localOrient1.LineWidth = 0.5; %(default: 0.5)
quiv_localOrient1.AutoScaleFactor = 1.0 ;
quiv_localOrient1.Color = col_localOrient;
if centerQuiverOnPixels
%give the duplicate quiver the same properties, especially the same color
quiv_localOrient2.ShowArrowHead = quiv_localOrient1.ShowArrowHead ;
quiv_localOrient2.LineWidth = quiv_localOrient1.LineWidth ;
quiv_localOrient2.AutoScaleFactor = quiv_localOrient1.AutoScaleFactor ;
quiv_localOrient2.Color = quiv_localOrient1.Color ;
end
end
if show_anisotropy
% plot projection of the anisotropy outside of the plot area
quiv_anisotropy = quiver( (size_x_orig+15)*[1 1], 0.5*(size_y_orig)*[1 1], ...
[AnisotropyAxis_inThePlane_MATLABcoords(1), -AnisotropyAxis_inThePlane_MATLABcoords(1)], ...
[AnisotropyAxis_inThePlane_MATLABcoords(2), -AnisotropyAxis_inThePlane_MATLABcoords(2)] );
quiv_anisotropy.ShowArrowHead = 'off';
quiv_anisotropy.LineWidth = 5;
quiv_anisotropy.AutoScaleFactor = 10;
quiv_anisotropy.Color = col_anisotropy;
text( (size_x_orig+15), 0.5*(size_y_orig)-20, ...
'anisotropy axis (projected)')
text( (size_x_orig+15), 0.5*(size_y_orig)+20, [{'bias into anisotropy direction' , 'among orientations where box was within foreground', ['linear = ', num2str(alignment1_magn_intens_fg_inner)], ['quadratic = ', num2str(alignment2_magn_intens_fg_inner)] }])
text( (size_x_orig+15), 0.5*(size_y_orig)+60, ...
{'angle between', 'anisotropy and projection', ['= ' num2str(AngleAnisotropyProjection*180/pi) ' degree']} )
end
axis equal
title(['orig image with foreground highlighted, orientations'] )
drawnow
pause(0.1)
if show_anisotropy
%make more space for title and annotations on the right side
currAx = gca();
currAx.Position = currAx.Position - [0 0 0.05 0.05];
end
print(Fig12, [filename_intermediate_output(1:end-4) '_Orientations'],'-dpdf', '-r0', '-fillpage')
print(Fig12, [filename_intermediate_output(1:end-4) '_Orientations'],'-dpng', '-r0')
close all
%% plotting the distribution of orientations to check for impact of anisotropy in imaging
%% binning
% bin the data explicitly to allow for weights being considered
% bin into interval [0,pi], then duplicate data from theta to pi+theta to represent nematics
NbinsPolar = 128;
theta_for_hists = mod(theta, pi);
assert ( min(theta_for_hists(:))>= 0 & max(theta_for_hists(:))<=pi, 'theta falls out of range [0, pi]' )
EdgesPolar = linspace(0, pi, NbinsPolar+1);
% this ensures to have all data fall into the bins (NaN bins are only due to theta value NaN)
assert ( EdgesPolar(1) == 0 & EdgesPolar(end) == pi , 'Incorrect bins')
EdgesNematic = [EdgesPolar(1:end-1) pi+EdgesPolar];
whichBin = discretize(theta_for_hists, EdgesPolar);
Freq = zeros([1 numel(EdgesPolar)-1]);
Freq_fg = zeros([1 numel(EdgesPolar)-1]);
Freq_fg_inner = zeros([1 numel(EdgesPolar)-1]);
%Freq_fg_innermost = zeros([1 numel(EdgesPolar)-1]);
% _innermost does not add information, because it makes sure that coherency estimation is within the foreground
Freq_magn = zeros([1 numel(EdgesPolar)-1]);
Freq_magn_fg = zeros([1 numel(EdgesPolar)-1]);
Freq_magn_fg_inner = zeros([1 numel(EdgesPolar)-1]);
Freq_intens = zeros([1 numel(EdgesPolar)-1]);
Freq_intens_fg = zeros([1 numel(EdgesPolar)-1]);
Freq_intens_fg_inner = zeros([1 numel(EdgesPolar)-1]);
Freq_magn_intens = zeros([1 numel(EdgesPolar)-1]);
Freq_magn_intens_fg = zeros([1 numel(EdgesPolar)-1]);
Freq_magn_intens_fg_inner = zeros([1 numel(EdgesPolar)-1]);
for kk=1:numel(EdgesPolar)-1
Freq(kk) = sum( sum( (whichBin==kk) ));
Freq_fg(kk) = sum( sum( (whichBin==kk) .* fg_closed_downSampl));
Freq_fg_inner(kk) = sum( sum( (whichBin==kk) .* fg_inner_downSampl));
%Freq_fg_innermost(kk)= sum( sum( (whichBin==kk) .* fg_innermost_downSampl));
Freq_magn(kk) = nansum(nansum( (whichBin==kk) .* magn ));
Freq_magn_fg(kk) = nansum(nansum( (whichBin==kk) .* magn .* fg_closed_downSampl));
Freq_magn_fg_inner(kk) = nansum(nansum( (whichBin==kk) .* magn .* fg_inner_downSampl ));
Freq_intens(kk) = sum( sum( (whichBin==kk) .* GrayImageStack_downSamplTwice));
Freq_intens_fg(kk) = sum( sum( (whichBin==kk) .* GrayImageStack_downSamplTwice .* fg_closed_downSampl));
Freq_intens_fg_inner(kk) = sum( sum( (whichBin==kk) .* GrayImageStack_downSamplTwice .* fg_inner_downSampl ));
Freq_magn_intens(kk) = nansum(nansum( (whichBin==kk) .* magn .* GrayImageStack_downSamplTwice));
Freq_magn_intens_fg(kk) = nansum(nansum( (whichBin==kk) .* magn .* GrayImageStack_downSamplTwice .* fg_closed_downSampl));
Freq_magn_intens_fg_inner(kk) = nansum(nansum( (whichBin==kk) .* magn .* GrayImageStack_downSamplTwice .* fg_inner_downSampl ));
end
%% plotting of orientation distribution itself
colUse_All = 'black' ;
colUse_fg = [0 0.8 0]; %dark green
colUse_fg_inner = 'blue';
colUse_fg_innermost = 'red';
Fig13 = figure(13); %'PaperSize', [11 8.5]) %'Position', [100 100 300 500]);
Fig13.Position = [0, 0, 1200, 900] ;
Fig13A = polaraxes('Units', 'normalized', 'OuterPosition', [ 0 1/2 1/2 1/2]) ;
Fig13B = polaraxes('Units', 'normalized', 'OuterPosition', [ 1/2 1/2 1/2 1/2]) ;
Fig13C = polaraxes('Units', 'normalized', 'OuterPosition', [ 0 0 1/2 1/2]) ;
Fig13D = polaraxes('Units', 'normalized', 'OuterPosition', [ 1/2 0 1/2 1/2]) ;
subplot(Fig13A); %makes the existing subplot current
hold on
pHist = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq Freq ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_All);
pHist_fg = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_fg Freq_fg ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg);
pHist_fg_inner = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_fg_inner Freq_fg_inner ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg_inner);
% pHist_fg_innermost= polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_fg_innermost Freq_fg_innermost], ...
% 'DisplayStyle','stairs', 'EdgeColor', colUse_fg_innermost);
% add projection of the anisotropy
pAnisotropyAxis = polarplot( [AngleAnisotropy_inThePlane_NORMALcoords, AngleAnisotropy_inThePlane_NORMALcoords+pi], [0.8, 0.8]* max(Freq), 'Color', 'red' );
legend([pHist, pHist_fg, pHist_fg_inner, ... %, pHist_fg_innermost,...
pAnisotropyAxis], ...
{'all orientations', 'orientations in foreground', 'where box was within foreground', ... % 'where coherency box was within foreground', ...
'axis of anisotropy (projected)'}, ...
'Location', 'southwest')
text(-pi/4, max(Freq), [{'bias into anisotropy direction' , 'among orientations where box was within foreground', ['linear = ', num2str(alignment1_fg_inner)], ['quadratic = ', num2str(alignment2_fg_inner)] }])
text( pi/4, 1.5*max(Freq), ...
{'angle between', 'anisotropy and projection', ['= ' num2str(AngleAnisotropyProjection*180/pi) ' degree']} )
title('frequency of orientations')
subplot(Fig13B); %makes the existing subplot current
hold on
pHist_magn = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_magn Freq_magn ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_All);
pHist_magn_fg = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_magn_fg Freq_magn_fg ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg);
pHist_magn_fg_inner = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_magn_fg_inner Freq_magn_fg_inner ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg_inner);
pAnisotropyAxis = polarplot( [AngleAnisotropy_inThePlane_NORMALcoords, AngleAnisotropy_inThePlane_NORMALcoords+pi], [0.8, 0.8]* max(Freq_magn), 'Color', 'red' );
%legend([pHist_magn, pHist_magn_fg, pHist_magn_fg_inner, ...
% pAnisotropyAxis], ...
% {'all orientations', 'orientations in foreground', 'where box was within foreground', ...
% 'axis of anisotropy (projected)'}, ...
% 'Location', 'southwest')
text(-pi/4, max(Freq_magn), [{'bias into anisotropy direction' , 'among orientations where box was within foreground', ['linear = ', num2str(alignment1_magn_fg_inner)], ['quadratic = ', num2str(alignment2_magn_fg_inner)] }])
title('frequency of orientations - magnitude weighted')
subplot(Fig13C); %makes the existing subplot current
hold on
pHist_intens = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_intens Freq_intens ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_All);
pHist_intens_fg = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_intens_fg Freq_intens_fg ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg);
pHist_intens_fg_inner = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_intens_fg_inner Freq_intens_fg_inner ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg_inner);
pAnisotropyAxis = polarplot( [AngleAnisotropy_inThePlane_NORMALcoords, AngleAnisotropy_inThePlane_NORMALcoords+pi], [0.8, 0.8]* max(Freq_intens), 'Color', 'red' );
%legend([pHist_intens, pHist_intens_fg, pHist_intens_fg_inner, ...
% pAnisotropyAxis], ...
% {'all orientations', 'orientations in foreground', 'where box was within foreground', ...
% 'axis of anisotropy (projected)'}, ...
% 'Location', 'southwest')
text(-pi/4, max(Freq_intens), [{'bias into anisotropy direction' , 'among orientations where box was within foreground', ['linear = ', num2str(alignment1_intens_fg_inner)], ['quadratic = ', num2str(alignment2_intens_fg_inner)] }])
title('frequency of orientations - intensity weighted')
subplot(Fig13D); %makes the existing subplot current
hold on
pHist_magn_intens = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_magn_intens Freq_magn_intens ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_All);
pHist_magn_intens_fg = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_magn_intens_fg Freq_magn_intens_fg ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg);
pHist_magn_intens_fg_inner = polarhistogram('BinEdges',EdgesNematic,'BinCounts',[Freq_magn_intens_fg_inner Freq_magn_intens_fg_inner ], ...
'DisplayStyle','stairs', 'EdgeColor', colUse_fg_inner);
pAnisotropyAxis = polarplot( [AngleAnisotropy_inThePlane_NORMALcoords, AngleAnisotropy_inThePlane_NORMALcoords+pi], [0.8, 0.8]* max(Freq_magn_intens), 'Color', 'red' );
%legend([pHist_magn_intens, pHist_magn_intens_fg, pHist_magn_intens_fg_inner, ...
% pAnisotropyAxis], ...
% {'all orientations', 'orientations in foreground', 'where box was within foreground', ...
% 'axis of anisotropy (projected)'}, ...
% 'Location', 'southwest')
text(-pi/4, max(Freq_magn_intens), [{'bias into anisotropy direction' , 'among orientations where box was within foreground', ['linear = ', num2str(alignment1_magn_intens_fg_inner)], ['quadratic = ', num2str(alignment2_magn_intens_fg_inner)] }])
title('frequency of orientations - magnitude and intensity weighted')
drawnow
pause(0.1)
print(Fig13, [filename_intermediate_output(1:end-4) '_Anisotropy' ],'-dpdf', '-r0', '-fillpage')
print(Fig13, [filename_intermediate_output(1:end-4) '_Anisotropy' ],'-dpng', '-r0')
close all
%% plotting of local coherency and locally dominant orientation
% NOTE: this plot is only possible at the resolution of the orientation
% estimate, that is after the downsampling
[x14, y14] = meshgrid(1:size(localCoherency,2), 1:size(localCoherency,1) );
u_localOrient = localOrient(:,:,1) .* localCoherency;
v_localOrient = localOrient(:,:,2) .* localCoherency;
Fig14 = figure('PaperOrientation', 'landscape', 'Position', [0, 0, 1200, 700]);
[axesHandles14, positions14] = tight_subplot(1,1,[.08 .05],[.05 .01],[.05 .02]) ; %gap, margin_height, margin_width)
fig14A = axesHandles14(1);
axes(fig14A) %make the existing axes current
OrigData14 = imshow(localCoherency); % (myData * 2 - 3);
hold on
if true % if show_localOrient
quiv_localOrient1 = quiver(x14(1:spacing12:end,1:spacing12:end), y14(1:spacing12:end,1:spacing12:end), u_localOrient(1:spacing12:end,1:spacing12:end), -v_localOrient(1:spacing12:end,1:spacing12:end));
%Note negative sign on y-component to transform from theta in normal coord layout to coord layout used by MATLAB (y downwards in plotting that is started by imshow() )
if centerQuiverOnPixels
%duplicate quiver (with opposite sign) to have coherent appearance
quiv_localOrient2 = quiver(x14(1:spacing12:end,1:spacing12:end), y14(1:spacing12:end,1:spacing12:end), -u_localOrient(1:spacing12:end,1:spacing12:end), v_localOrient(1:spacing12:end,1:spacing12:end));
end
% modify quiver plot(s)
quiv_localOrient1.ShowArrowHead = 'off';
quiv_localOrient1.LineWidth = 0.5; %(default: 0.5)
quiv_localOrient1.AutoScaleFactor = 1.0 ;
quiv_localOrient1.Color = col_localOrient;
if centerQuiverOnPixels
%give the duplicate quiver the same properties, especially the same color
quiv_localOrient2.ShowArrowHead = quiv_localOrient1.ShowArrowHead ;
quiv_localOrient2.LineWidth = quiv_localOrient1.LineWidth ;
quiv_localOrient2.AutoScaleFactor = quiv_localOrient1.AutoScaleFactor ;
quiv_localOrient2.Color = quiv_localOrient1.Color ;
end
end
axis equal
%{
cbar = colorbar;
cbar.Label.String = 'Local Coherency';
%}
title(['locally dominant orientation (within radius = ' num2str(boxRadius_LocalAlignment) '(orig pixels) at power ' num2str(powerOfMagn_LocalAlignment) ')' ])
drawnow
pause(0.1)
print(Fig14, [filename_intermediate_output(1:end-4) '_coherencyRadius' num2str(boxRadius_LocalAlignment) 'Power' num2str(powerOfMagn_LocalAlignment) ],'-dpdf', '-r0', '-fillpage')
print(Fig14, [filename_intermediate_output(1:end-4) '_coherencyRadius' num2str(boxRadius_LocalAlignment) 'Power' num2str(powerOfMagn_LocalAlignment) ],'-dpng', '-r0')
close all
end
%% save collected coherency for each ROI to file output_collected as .mat and as .csv for subsequent statistical analysis
% basically write all_vals_mat with descr as its header, but precede by
% information on the respective path stored in folders_arr
%recall: the code assumes that folder_to_subfolders
% contains <Stage>/<Embryo>/<RegionOfTheEye>/<DerotationConditions>/
% and the different ROIs within one folder are distinguished by leading_ROIinfo_inFilename
% save as .mat file
save([ output_collected(1:end-4) '.mat'], 'descr', 'folders_arr', 'all_vals_mat')
% save as .csv
% write the header
fileID = fopen(output_collected,'w');
fprintf(fileID, ['full path' '\t' 'stage' '\t' 'embryo' '\t' 'region of the eye' '\t' ...
'derotation conditions' '\t' 'ROI'] );
fieldnameslist = fieldnames(descr);
for kk = 1:numel(fieldnameslist)
fprintf(fileID, [descr.(fieldnameslist{kk}) '\t'] );
end
fprintf(fileID, '\n');
fclose(fileID);
% write the content
for kk=1:size(folders_arr,1)
%write filenames
fileID = fopen(output_collected,'a');
fprintf(fileID, [ folders_arr{kk,1} '\t' ...
folders_arr{kk,2} '\t' ...
folders_arr{kk,3} '\t' ...
folders_arr{kk,4} '\t' ...
folders_arr{kk,5} '\t' ...
folders_arr{kk,6} '\t' ]);
fclose(fileID);
%write statistics
dlmwrite( output_collected , all_vals_mat(kk,:), '-append', 'delimiter','\t') %collected data
end
AnalyseOrientations/ComputeDirecty.m 0 → 100644
View file @ fc349606
function [Directy, subpatchesChangeDirecty] = ComputeDirecty(anArray, boxRadius, forCenterOfMass, subtractBackground, blackBackground)
% Compute an orientation estimate (also known as directionality) for input
% matrix of grayscale image intensities,
% using the eigenvectors of the tensor of inertia (ToI) of the intensities
% on square image patches of user-specified side length 2*boxRadius+1.
% For input of size (m,n), the return value Directy is of size (m,n,3).
% For each pixel (i,j), Directy(i,j,3) contains the strength of directionality
% and the corresponding normalised eigenvector is returned in Directy(i,j,1:2).
% where Directy(i,j,1) is the x-component and Directy(i,j,2) is the y-component
% in *standard* coordinate layout (while the indexing (i,j) is in MATLAB coordinate layout).
%
% The strength of the orientation is defined
% as (\lambda_max-\lambda_min) / (\lambda_max+\lambda_min)
% where \lambda_max >= \lambda_min are the two eigenvalues of the ToI
%
% Arguments:
% - boxRadius [not necessarily integer, >0.5] sets the scale: The ToI is computed from a
% roughly circular patch with radius boxRadius, and the result is assigned to the center pixel.
% [old: The ToI is computed from a box of size 2*boxRadius+1, where boxRadius
% needed to be integer, and the result was assigned to the center pixel. ]
% The return array Directy(i,j,1:3) contains NaN within distance
% boxRadius from the boundaries.
% The second return value subpatchesChangeDirecty is only experimental: it
% is the same size as the input anArray and contain logicals that
% indicate whether a consideration of subpatches changes the orientation
% estimate to a roughly perpendicular direction (which can be the case
% where the original patch touches two bright foreground objects only at
% its boundary). Even if its value is true, the returned Directy is not
% changed!
%
% Input arguments are:
% - forCenterOfMass [boolean, optional, default=true]: whether to base the
% ToI on the geometric center [false] or on the center of mass (say, center
% of intensity) [true].
% - subtractBackground [boolean, optional, default=true]: Whether to
% subtract the minimal intensity value occuring in the current box prior to
% the computation of ToI
% Recommended when forCenterOfMass=true, as the background intensity does
% not cancel out for a point of reference unequal the geometric center
% - blackBackground [boolean, optional, default=true]: Set true for
% high-intensity objects on low intensity background, false for the inverse
%
% Iteration over the small boxes is with parallelisation (parfor)
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% last edit: 2022-08-25
%
% CHANGELOG
% 2022-08-25 clean-up of experimental subpatchesChangeDirecty
% 2021-03-05 leave NaN values where no orientation can be computed,
% especially at the image boundary (instead of setting all
% components of Directy to zero)
%% preparations
% sanity checks
assert( size( size(squeeze(anArray)), 2) <=2 , 'Input array of ComputeDirecty must be essentially 2-dimensional at most' )
assert(boxRadius > 0.5, 'Input boxRadius to ComputeDirecty must be larger than 0.5 to include more than one pixel')
if min(anArray(:)) < 0
warning('Negative entries in input array of ComputeDirecty may cause unexpected results')
end
% set optional arguments if not provided
if (~exist('forCenterOfMass', 'var'))
forCenterOfMass = true;
end
if (~exist('subtractBackground', 'var'))
subtractBackground = true;
end
if (~exist('blackBackground', 'var'))
blackBackground = true;
end
% Prepare a mask of a roughly circular region with radius boxRadius:
% Determine where a circle of radius boxRadius would intersect the boundaries of the pixels.
% whereby a pixel (i,j) is covering the space [i-0.5,i+0.5]x[j-0.5,j+0.5] in \mathbb{R}^2
% If there are intersection points, connect them linearly and compute the area of the resulting polygon
% The result is equal in all four quadrants.
% prelim_mask is the mask for the upper right quadrant, with prelim_mask(1,1) representing the pixel [-0.5,0.5]x[-0.5,0.5] around the origin
% Here, assume the first coordinate (i) to count horizontally and the second (j) to count vertically [this is contrary to MATLAB default]
mask = CircularMask(boxRadius);
outerBoxRadius = (size(mask,1)-1) /2; %outerBoxRadius will preselect a square of pixels covering the [roughly] circular mask
% alternative to restore old behaviour with square image patch:
%{
outerBoxRadius = ceil(boxRadius);
mask = ones( [2*outerBoxRadius+1, 2*outerBoxRadius+1] );
%}
i_max = size(anArray,1);
j_max = size(anArray,2);
Directy = NaN( i_max, j_max, 3 );
subpatchesChangeDirecty = NaN ( i_max, j_max );
%% computation
% prepare kernel templates to multiply the actual intensities with
% distances as row and column vector
% Recall that array elements are accessed as array(row_number, column_number) in MATLAB, which is
% anImage(y,x) for pixel data. Take this into account for the computation of distances
% row distances = y-distances
% col distances = x-distances
dist_row_forOrigin = -outerBoxRadius:outerBoxRadius ;
dist_col_forOrigin = transpose(-outerBoxRadius:outerBoxRadius);
template_Ixx_forOrigin = repmat(dist_row_forOrigin .* dist_row_forOrigin, 2*outerBoxRadius+1, 1) ;
template_Iyy_forOrigin = repmat(dist_col_forOrigin .* dist_col_forOrigin, 1, 2*outerBoxRadius+1) ;
template_Ixy_forOrigin = - dist_col_forOrigin * dist_row_forOrigin;
% prepare outside parfor loop
Directy_1 = NaN(i_max,j_max);
Directy_2 = NaN(i_max,j_max);
Directy_3 = NaN(i_max,j_max);
% pointwise tensor of inertia eigendecomposition
parfor i=outerBoxRadius+1:i_max-outerBoxRadius
RowOfOrientations_1 = NaN(1,j_max);
RowOfOrientations_2 = NaN(1,j_max);
RowOfStrenths = NaN(1,j_max);
RowOfsubpatchesChangeDirecty = NaN(1,j_max);
for j=outerBoxRadius+1:j_max-outerBoxRadius
localArray = anArray(i-outerBoxRadius:i+outerBoxRadius, j-outerBoxRadius:j+outerBoxRadius);
if subtractBackground
localArray = localArray - min(min(localArray)) ;
end
% masking
maskedArray = mask .* localArray;
if forCenterOfMass
%compute the tensor of inertia for the center of mass (which minimizes the eigenvalues, see parallel axis theorem [Steiner'scher Satz])
% this option should be combined with "subtractBackground", as the background intensity does not cancel out for point of reference unequal the geometric center
totalMass = sum(sum(maskedArray)) ;
if totalMass == 0 %if there is no "mass" (image intensity), set the "center of mass" to the geometric center
centerOfMass_x = 0; % corresponds to index outerBoxRadius + 1;
centerOfMass_y = 0; % corresponds to index outerBoxRadius + 1;
else
centerOfMass_x = sum(sum( maskedArray .* repmat(dist_row_forOrigin , 2*outerBoxRadius+1, 1 ) )) / totalMass ;
centerOfMass_y = sum(sum( maskedArray .* repmat(dist_col_forOrigin , 1, 2*outerBoxRadius+1 ) )) / totalMass ;
end
dist_row_forCoM = dist_row_forOrigin - centerOfMass_x ;
dist_col_forCoM = dist_col_forOrigin - centerOfMass_y ;
template_Ixx = repmat(dist_row_forCoM .* dist_row_forCoM, 2*outerBoxRadius+1, 1) ;
template_Iyy = repmat(dist_col_forCoM .* dist_col_forCoM, 1, 2*outerBoxRadius+1) ;
template_Ixy = - dist_col_forCoM * dist_row_forCoM;
else % compute the tensor of inertia for the geometric center
template_Ixx = template_Ixx_forOrigin;
template_Iyy = template_Iyy_forOrigin;
template_Ixy = template_Ixy_forOrigin;
% set the "center of mass" to the geometric center
centerOfMass_x = 0; % corresponds to index outerBoxRadius + 1;
centerOfMass_y = 0; % corresponds to index outerBoxRadius + 1;
dist_row_forCoM = dist_row_forOrigin - centerOfMass_x ;
dist_col_forCoM = dist_col_forOrigin - centerOfMass_y ;
end
tensorOfInertia = [ sum(sum(maskedArray .* template_Ixx )) sum(sum(maskedArray .* template_Ixy )) ; ...
sum(sum(maskedArray .* template_Ixy )) sum(sum(maskedArray .* template_Iyy )) ] ;
[RightEigVecs, EigValVec, LeftEigVecs] = eig(tensorOfInertia, 'vector'); %the optional argument 'vector' makes the eigenvalues be returned as vector
[MaxEigVal, PosMaxEigVal] = max(EigValVec);
EigVecOfMaxEigVal = RightEigVecs(:, PosMaxEigVal);
MinEigVal = min(EigValVec);
EigVecOfMinEigVal = RightEigVecs(:, 3-PosMaxEigVal); % extract the Eigenvector by the position explicitly computed as "the other position out of the set {1,2}"
%{
%% experimental
% decide for each box whether you see bright object on dark background
% or dark "object" an bright background (which may occur
% p.e. when the box covers the dark space between two elongated bright objects and a bit of them)
% NOTE: This makes additional assumptions on the amount of
foreground / background pixels, but the crucial point is their spatial distribution
otsuThresh = max(maskedArray(:)) * graythresh( maskedArray/max(maskedArray(:)) ); %threshold from Otsu's method that minimizes intraclass variations - assuming there are bright (foreground) and dark (background) regions
if sum( maskedArray(:) > otsuThresh ) > sum( maskedArray(:) < otsuThresh )
%there are more pixels above the threshold => rather dark object on bright background
blackBackground = 1 - blackBackground;
end
% end of experimental
%}
% preliminary Directionality result
prelimDirecty = zeros([1,3]);
if blackBackground
prelimDirecty(1:2) = EigVecOfMaxEigVal(:); %return a vector orthogonal to the eigenvector of smallest eigenvalue
% because of a positive semi-definit, symmetric matrix, the eigenvectors are orthogonal, so we return the eigenvector to the larger eigenvalue
else
prelimDirecty(1:2) = EigVecOfMinEigVal(:); %return a vector orthogonal to the eigenvector of largest eigenvalue
% because of a positive semi-definit, symmetric matrix, the eigenvectors are orthogonal, so we return the eigenvector to the smaller eigenvalue
end
prelimDirecty(3) = (MaxEigVal - MinEigVal) / (MaxEigVal + MinEigVal) ;
%% experimental: check whether orientation estimate is testified by subpatches
% when splitting the domain orthogonal to prelimDirecty(1:2) - this helps adjust the
% Directy for example when intensity is distributed like a bar-bell
% splitting orthogonal to prelimDirecty(1:2) by scalar product <0 or >0
% NOTE that we have to use a negative sign with the y-component prelimDirecty(1)
%mask_part1 = mask .* (0 < -prelimDirecty(1) * repmat(dist_row_forCoM , 2*outerBoxRadius+1, 1 ) + prelimDirecty(2) * repmat(dist_col_forCoM , 1, 2*outerBoxRadius+1 ) );
%mask_part2 = mask .* (0 > -prelimDirecty(1) * repmat(dist_row_forCoM , 2*outerBoxRadius+1, 1 ) + prelimDirecty(2) * repmat(dist_col_forCoM , 1, 2*outerBoxRadius+1 ) );
% However, the subregions above are anisotropic, which may again introduce biases
% Therefore define two roughly circular patches of half the boxRadius
% located in each of the two halfs described above
submask_offset_y = -prelimDirecty(1)*0.5*boxRadius;
submask_offset_x = prelimDirecty(2)*0.5*boxRadius;
mask_part1 = mask .* ( (0.5*boxRadius)^2 > (repmat(dist_row_forCoM + submask_offset_y , 2*outerBoxRadius+1, 1 )).^2 + ...
(repmat(dist_col_forCoM + submask_offset_x , 1, 2*outerBoxRadius+1 )).^2 );
mask_part2 = mask .* ( (0.5*boxRadius)^2 > (repmat(dist_row_forCoM - submask_offset_y , 2*outerBoxRadius+1, 1 )).^2 + ...
(repmat(dist_col_forCoM - submask_offset_x , 1, 2*outerBoxRadius+1 )).^2 );
%% analyse tensor of inertia in each half as above for the full maskedArray
maskedArray1 = mask_part1 .* localArray;
% in the parts, always work for center of mass
%if forCenterOfMass
%compute the tensor of inertia for the center of mass (which minimizes the eigenvalues, see parallel axis theorem [Steiner'scher Satz])
% this option should be combined with "subtractBackground", as the background intensity does not cancel out for point of reference unequal the geometric center
totalMass1 = sum(sum(maskedArray1)) ;
if totalMass1 == 0 %if there is no mass, set the "center of mass" to the geometric center
centerOfMass_x1 = outerBoxRadius + 1;
centerOfMass_y1 = outerBoxRadius + 1;
else
centerOfMass_x1 = sum(sum( maskedArray1 .* repmat(dist_row_forOrigin , 2*outerBoxRadius+1, 1 ) )) / totalMass1 ;
centerOfMass_y1 = sum(sum( maskedArray1 .* repmat(dist_col_forOrigin , 1, 2*outerBoxRadius+1 ) )) / totalMass1 ;
end
dist_row_forCoM1 = dist_row_forOrigin - centerOfMass_x1 ;
dist_col_forCoM1 = dist_col_forOrigin - centerOfMass_y1 ;
template_Ixx1 = repmat(dist_row_forCoM1 .* dist_row_forCoM1, 2*outerBoxRadius+1, 1) ;
template_Iyy1 = repmat(dist_col_forCoM1 .* dist_col_forCoM1, 1, 2*outerBoxRadius+1) ;
template_Ixy1 = - dist_col_forCoM1 * dist_row_forCoM1;
%else % compute the tensor of inertia for the geometric center
% template_Ixx = template_Ixx_forOrigin;
% template_Iyy = template_Iyy_forOrigin;
% template_Ixy = template_Ixy_forOrigin;
% % set the "center of mass" to the geometric center
% centerOfMass_x = outerBoxRadius + 1;
% centerOfMass_y = outerBoxRadius + 1;
% dist_row_forCoM = dist_row_forOrigin - centerOfMass_x ;
% dist_col_forCoM = dist_col_forOrigin - centerOfMass_y ;
%end
tensorOfInertia1 = [ sum(sum(maskedArray1 .* template_Ixx1 )) sum(sum(maskedArray1 .* template_Ixy1 )) ; ...
sum(sum(maskedArray1 .* template_Ixy1 )) sum(sum(maskedArray1 .* template_Iyy1 )) ] ;
[RightEigVecs1, EigValVec1, LeftEigVecs1] = eig(tensorOfInertia1, 'vector'); %the optional argument 'vector' makes the eigenvalues be returned as vector
[MaxEigVal1, PosMaxEigVal1] = max(EigValVec1);
EigVecOfMaxEigVal1 = RightEigVecs1(:, PosMaxEigVal1);
MinEigVal1 = min(EigValVec1);
EigVecOfMinEigVal1 = RightEigVecs1(:, 3-PosMaxEigVal1); % extract the Eigenvector by the position explicitly computed as "the other position out of the set {1,2}"
maskedArray2 = mask_part2 .* localArray;
% in the parts, always work for center of mass
%if forCenterOfMass
%compute the tensor of inertia for the center of mass (which minimizes the eigenvalues, see parallel axis theorem [Steiner'scher Satz])
% this option should be combined with "subtractBackground", as the background intensity does not cancel out for point of reference unequal the geometric center
totalMass2 = sum(sum(maskedArray2)) ;
if totalMass2 == 0 %if there is no mass, set the "center of mass" to the geometric center
centerOfMass_x2 = outerBoxRadius + 1;
centerOfMass_y2 = outerBoxRadius + 1;
else
centerOfMass_x2 = sum(sum( maskedArray2 .* repmat(dist_row_forOrigin , 2*outerBoxRadius+1, 1 ) )) / totalMass2 ;
centerOfMass_y2 = sum(sum( maskedArray2 .* repmat(dist_col_forOrigin , 1, 2*outerBoxRadius+1 ) )) / totalMass2 ;
end
dist_row_forCoM2 = dist_row_forOrigin - centerOfMass_x2 ;
dist_col_forCoM2 = dist_col_forOrigin - centerOfMass_y2 ;
template_Ixx2 = repmat(dist_row_forCoM2 .* dist_row_forCoM2, 2*outerBoxRadius+1, 1) ;
template_Iyy2 = repmat(dist_col_forCoM2 .* dist_col_forCoM2, 1, 2*outerBoxRadius+1) ;
template_Ixy2 = - dist_col_forCoM2 * dist_row_forCoM2;
%else % compute the tensor of inertia for the geometric center
% template_Ixx = template_Ixx_forOrigin;
% template_Iyy = template_Iyy_forOrigin;
% template_Ixy = template_Ixy_forOrigin;
% % set the "center of mass" to the geometric center
% centerOfMass_x = outerBoxRadius + 1;
% centerOfMass_y = outerBoxRadius + 1;
% dist_row_forCoM = dist_row_forOrigin - centerOfMass_x ;
% dist_col_forCoM = dist_col_forOrigin - centerOfMass_y ;
%end
tensorOfInertia2 = [ sum(sum(maskedArray2 .* template_Ixx2 )) sum(sum(maskedArray2 .* template_Ixy2 )) ; ...
sum(sum(maskedArray2 .* template_Ixy2 )) sum(sum(maskedArray2 .* template_Iyy2 )) ] ;
[RightEigVecs2, EigValVec2, LeftEigVecs2] = eig(tensorOfInertia2, 'vector'); %the optional argument 'vector' makes the eigenvalues be returned as vector
[MaxEigVal2, PosMaxEigVal2] = max(EigValVec2);
EigVecOfMaxEigVal2 = RightEigVecs2(:, PosMaxEigVal2);
MinEigVal2 = min(EigValVec2);
EigVecOfMinEigVal2 = RightEigVecs2(:, 3-PosMaxEigVal2); % extract the Eigenvector by the position explicitly computed as "the other position out of the set {1,2}"
% check whether the Directy's from the parts coincide with the Directy from the full patch.
% For Comparison use absolute value of the scalar product
if ( abs( EigVecOfMaxEigVal(1)*EigVecOfMaxEigVal1(1) + EigVecOfMaxEigVal(2)*EigVecOfMaxEigVal1(2) ) + ...
abs( EigVecOfMaxEigVal(1)*EigVecOfMaxEigVal2(1) + EigVecOfMaxEigVal(2)*EigVecOfMaxEigVal2(2) ) ...
> ...
abs(-EigVecOfMaxEigVal(1)*EigVecOfMaxEigVal1(2) + EigVecOfMaxEigVal(2)*EigVecOfMaxEigVal1(1) ) + ... % here take a vector orthogonal to EigVecOfMaxEigVal1 in scalar prduct with prelimDirecty
abs(-EigVecOfMaxEigVal(1)*EigVecOfMaxEigVal2(2) + EigVecOfMaxEigVal(2)*EigVecOfMaxEigVal2(1) ) )
% in this case, the patches confirm the result of prelimDirecty
RowOfsubpatchesChangeDirecty(j) = false;
% DO NOT ALTER the assigned orientations
%RowOfOrientations_1 (j) = prelimDirecty(1);
%RowOfOrientations_2 (j) = prelimDirecty(2);
else % the vector orthogonal to prelimDirecty might be a better fit,
RowOfsubpatchesChangeDirecty(j) = true;
% but DO NOT ASSIGN other values to to RowOfOrientations_1 (and thereby to Directy)
%RowOfOrientations_1 (j) = prelimDirecty(2);
%RowOfOrientations_2 (j) = -prelimDirecty(1);
end
%% assign orientations and strength
%%% ALWAYS from prelimDirecty irrespective of subpatchesChangeDirecty
%%% no elementwise assignments %Directy(i, j, :) = prelimDirecty; inside the parfor loop
RowOfOrientations_1 (j) = prelimDirecty(1);
RowOfOrientations_2 (j) = prelimDirecty(2);
RowOfStrenths (j) = prelimDirecty(3);
end
Directy_1(i, :) = RowOfOrientations_1;
Directy_2(i, :) = RowOfOrientations_2;
Directy_3(i, :) = RowOfStrenths;
subpatchesChangeDirecty(i, :) = RowOfsubpatchesChangeDirecty;
end
Directy(:,:,1) = Directy_1;
Directy(:,:,2) = Directy_2;
Directy(:,:,3) = Directy_3;
AnalyseOrientations/ComputeDirectyGradient.m 0 → 100644
View file @ fc349606
function [Directy] = ComputeDirectyGradient(anArray, kernel_size)
% Compute an orientation estimate (also known as directionality) for input
% matrix of grayscale image intensities,
% using Scharr gradient estimation kernels of sizes 5 [recommended], 3 or 2
% that are optimized for low orientational error, see Tabelle B.11 in
%
% Hanno Scharr (2010): Optimale Operatoren in der digitalen
% Bildverarbeitung, PhD Thesis at University Heidelberg, Germany,
% DOI=10.11588/heidok.00000962, http://www.ub.uni-heidelberg.de/archiv/962
%
% The gradient gives the direction of steepest increase (here: of image
% intensity), therefore the direction perpendicular to it is returned as
% orientation of objeccts in the image.
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% last edit: 2022_08_25
%
%
% CHANGELOG
% 2022_08_25 clean-up
% 2021_03_15 computation of u and v simplified
if min(anArray(:)) < 0
warning('Negative entries in input array of ComputeDirecty detected, but gradient estimation not impeded.')
end
% do not enforce a normalization here, leave it to higher code levels
image_size = size(anArray);
if (kernel_size == 5)
%{
% Compare Directionality plugin in Fiji
filtering_kernel_x = [-2, -1, 0, 1, 2; -3, -2, 0, 2, 3; -4, -3, 0, 3, 4]/46;
filtering_kernel_x = [filtering_kernel_x ; flipud(filtering_kernel_x(1:2,:))];
filtering_kernel_y = rot90(filtering_kernel_x);
%}
% Use "5x5-opt" from [Scharr], see Tabelle B.11, and Abbildung 6.2.
% Note error of whopping 6.7e-05 (!)
derivative_kernel = [21.27, 85.46, 0, -85.46, -21.27]/256;
smoothing_kernel = [ 5.91, 61.77, 120.64, 61.77, 5.91]/256;
filtering_kernel_x= derivative_kernel .* smoothing_kernel' ;
filtering_kernel_y = rot90(filtering_kernel_x);
elseif (kernel_size == 2)
% The simplest kernel possible.
% See [Scharr], Tabelle B.11 and Abbildung 6.2 for error of 1.8e-02
% We here directly combine derivative ("Ableitung") and smoothing ("Glättung")
% NOTE that in the original work, the normalisation factor /2 is missing in the smoothing ("Glättung")
filtering_kernel_x = [ 1, -1; 1, -1] /2;
filtering_kernel_y = rot90(filtering_kernel_x);
elseif (kernel_size == 3)
% See [Scharr], Tabelle B.11 and Abbildung 6.2 for error of 2.2e-03
% both for "3x3-int" (which is optimised for low angle error in fixed point arithmetics)
% and for "3x3-opt" (which is optimised for low angle error in floating point arithmetics)
% We here directly combine derivative ("Ableitung") and smoothing ("Glättung")
% 3x3-int:
%filtering_kernel_x = [47, 0, -47; 162, 0, -162; 47, 0, -47]/512;
% 3x3-opt
filtering_kernel_x = [46.84, 0, -46.84; 162.32, 0, -162.32; 46.84, 0, -46.84]/512;
filtering_kernel_y = rot90(filtering_kernel_x);
else
error(['Please choose a kernel size 2, 3 or 5. Chosen kernel size was ' mat2str(kernel_size) '.']);
end
%{
% the commented version yields derivative approximations even at the very
% boundary pixels where the filter is not fully within the input image
convolved_image_x = conv2(anArray, filtering_kernel_x);
convolved_image_y = conv2(anArray, filtering_kernel_y);
radius_kernel = floor((size(filtering_kernel_x,1) - 1)/2);
%prune
convolved_image_x = convolved_image_x(radius_kernel + 1 : image_size(1) + radius_kernel, ...
radius_kernel + 1 : image_size(2) + radius_kernel);
convolved_image_y = convolved_image_y(radius_kernel + 1 : image_size(1) + radius_kernel, ...
radius_kernel + 1 : image_size(2) + radius_kernel);
%}
convolved_image_x_oversized = conv2(anArray, filtering_kernel_x);
convolved_image_y_oversized = conv2(anArray, filtering_kernel_y);
convolved_image_x = NaN(size(anArray));
convolved_image_y = NaN(size(anArray));
%leave values to NaN where the kernel was not fully within the input image
kernel_oversize_tl = ceil ((size(filtering_kernel_x,1) - 1)/2); %how many pixels the kernel has to the top or left of its center ("center" = pixel that the filtering result is assigned to)
kernel_oversize_br = floor((size(filtering_kernel_x,1) - 1)/2);
convolved_image_x(kernel_oversize_tl+1:end-kernel_oversize_br, kernel_oversize_tl+1:end-kernel_oversize_br) = ...
convolved_image_x_oversized(kernel_oversize_tl+kernel_oversize_br+1 : end-kernel_oversize_tl-kernel_oversize_br, kernel_oversize_tl+kernel_oversize_br+1 : end-kernel_oversize_tl-kernel_oversize_br);
convolved_image_y(kernel_oversize_tl+1:end-kernel_oversize_br, kernel_oversize_tl+1:end-kernel_oversize_br) = ...
convolved_image_y_oversized(kernel_oversize_tl+kernel_oversize_br+1 : end-kernel_oversize_tl-kernel_oversize_br, kernel_oversize_tl+kernel_oversize_br+1 : end-kernel_oversize_tl-kernel_oversize_br);
% L2 norm
magnitudesL2 = sqrt(convolved_image_x.^2 + convolved_image_y.^2);
% normalized components of the identified direction
%{
angles = atan2(convolved_image_y , convolved_image_x);
orthogonal_angles = angles + pi/2;
u = cos(orthogonal_angles);
v = sin(orthogonal_angles);
%}
% more direct computation of normalized components u and v
u = -convolved_image_y ./magnitudesL2 ; %normalization and rotation by pi/2 in one step
v = convolved_image_x ./magnitudesL2 ; % therefore "-y" assigned to u, and "+x" assigned to v
Directy(:,:,1) = u;
Directy(:,:,2) = v;
Directy(:,:,3) = magnitudesL2;
end
AnalyseOrientations/ComputeDirecty_FourierAndQ.m 0 → 100644
View file @ fc349606
function [Directy, subpatchesChangeDirecty] = ComputeDirecty_FourierAndQ(anArray, boxRadius, minWaveLength, maxWaveLength)
% Compute an orientation estimate (also known as directionality) for input
% matrix of grayscale image intensities,
% using the Q-tensor decomposition of the Fourier transform on square image
% patches of user-specified side length 2*boxRadius+1.
% For input of size (m,n), the return value Directy is of size (m,n,3).
% For each pixel (i,j), Directy(i,j,3) contains the strength of directionality
% and the corresponding normalised eigenvector is returned in Directy(i,j,1:2)
% where Directy(i,j,1) is the x-component and Directy(i,j,2) is the y-component
% in *standard* coordinate layout (while the indexing (i,j) is in MATLAB coordinate layout).
%
% The strength of the orientation is defined as the prefactor S in the
% Q-tensor decomposition Q = S * ( n * transpose(n) - 1/2 Id ), which is
% also known as the scalar nematic order parameter
%
% Arguments:
% - boxRadius sets the scale:
% The Fourier transform is applied to a box of size 2*boxRadius+1, where boxRadius
% needs to be integer, and the result is assigned to the center pixel.
% - minWaveLength and maxWaveLength restrict to an isotropic fraction of
% the Fourier space. (If there is no power at all in that part, the
% returned orientation and strength are NaN)
%
%
% improvement ideas:
% 1. The Fourier transform could be applied to a roughly circular image
% patch of [not necessarily integer] radius = boxRadius, and the result
% is assigned to the center pixel to reduce bias by the coordinate axis
% (already attenuated by integration only between specified wavelengths)
% 2. Polar Fourier Transform (p.e. PolarLab [Averbuch, Coifman, Donoho,
% Elad, Israeli (2004) Fast and accurate polar Fourier transform] or
% NFFT [Keiner, Kunis, Potts (2009) Using NFFT 3 - a software library
% for various nonequispaced fast Fourier transforms])
%
% Iteration over the small boxes is with parallelisation (parfor)
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
%% preparations
% sanity checks
assert( size( size(squeeze(anArray)), 2) <=2 , 'Input array of ComputeDirecty must be essentially 2-dimensional at most' )
assert(boxRadius > 0 && floor(boxRadius)==boxRadius, 'Input boxRadius to ComputeDirecty must be positive integer')
assert(minWaveLength < maxWaveLength, 'Empty selection of wavelengths to consider')
if maxWaveLength > 2*boxRadius+1
warning('maximal wavelength cannot be reached due to Fourier transform on box of smaller size')
end
if not (minWaveLength>0)
warning('minimal wavelength should be larger than 0')
% wavelength = 0 should be excluded, as otherwise the magnitude of the
% constant term would be assigned to the angle q_theta=0
minWaveLength = eps;
end
i_max = size(anArray,1);
j_max = size(anArray,2);
Directy = NaN( i_max, j_max, 3 );
subpatchesChangeDirecty = NaN( i_max, j_max );
outerBoxRadius = boxRadius;
%% run trial FFT2 on a patch in the middle to find best method, see https://de.mathworks.com/help/matlab/ref/fftw.html
fftw('dwisdom',[]); %Clear any existing double-precision transform parameters.
fftw('planner','measure'); %Set the method for optimizing Fourier transforms within the current MATLAB session to 'measure' which tries the common methods (instead of 'exhaustive' which tries ALL methods, but gives only ~5% gain [Peter Steinbach, MPI-CBG Scientific Computing Facility, based on "gearshifft", see https://github.com/mpicbg-scicomp/gearshifft, see also Steinbach, Werner (2017) gearshifft – The FFT Benchmark Suite for Heterogeneous Platforms].
fft2( anArray( ceil(0.5*i_max)-outerBoxRadius:ceil(0.5*i_max)+outerBoxRadius, ...
ceil(0.5*j_max)-outerBoxRadius:ceil(0.5*j_max)+outerBoxRadius) );
fftinfo = fftw('dwisdom'); %Assign the current double-precision transform algorithm parameters to the variable fftinfo for later use.
fftw('dwisdom',fftinfo); %Apply the parameter information stored in the variable fftinfo to future transform function calls.
%% prepare frequencies corresponding to the entries in the Fourier transform - note 0-based indexing of frequencies
% the calculation of q_x, q_y, q_theta, q_radius was moved out of the loop,
% because the size of the localArray is the same in each pass
[q_x, q_y] = meshgrid( 0:2*outerBoxRadius, 0:2*outerBoxRadius );
% shift q_x and q_y to accomodate the layout of fft2 in MATLAB
q_x = q_x - (q_x > 0.5*size(q_x,2)-0.25 ) * size(q_x,2); % negative frequencies are equivalent to frequencies > half domain size
q_y = q_y - (q_y > 0.5*size(q_y,1)-0.25 ) * size(q_y,1); % adjust here to get the correct q_theta from atan2( )
% the -0.25 is just to have the condition "q_x >= 0.5*size(q_x,2)" robust to roundoff errors
% express frequency in polar coordinates (q_radius, q_theta)
q_radius= sqrt(q_y.^2 + q_x.^2) ;
q_theta = mod( atan2(q_y, q_x), 2*pi); %force to interval [0, 2*pi)
%{
%check layout of q_theta:
imagesc(fftshift(q_theta))
%}
%% computation
Directy_1 = NaN(i_max,j_max);
Directy_2 = NaN(i_max,j_max);
Directy_3 = NaN(i_max,j_max);
parfor ii=outerBoxRadius+1:i_max-outerBoxRadius
RowOfOrientations_1 = NaN(1,j_max);
RowOfOrientations_2 = NaN(1,j_max);
RowOfStrenths = NaN(1,j_max);
for jj=outerBoxRadius+1:j_max-outerBoxRadius
localArray = anArray(ii-outerBoxRadius:ii+outerBoxRadius, jj-outerBoxRadius:jj+outerBoxRadius);
%Fourier transform
localArray_hat = fft2( localArray );
%% Q-tensor calculation
% might be further improved by a polar Fourier transform (see refs. at header)
% take only points into account within the given wavenumber limits
% the frequencies' components q_x, q_y and radial form (q_theta, q_radius) have been computed outside the loop
% normalise the power (or the absolute value of the Fourier transform) to a distribution
power_to_consider = abs( localArray_hat ) .* ( (q_radius > minWaveLength) & (q_radius < maxWaveLength)) ;
%power_to_consider = abs( localArray_hat ) .^2 .* ( (q_radius > minWaveLength) & (q_radius < maxWaveLength)) ;
% turn it into a distribution from which we compute the Q-tensor (= the second moment tensor of the distribution, except for -0.5*unitmatrix)
power_to_consider = power_to_consider / sum(power_to_consider(:));
Q_xx = sum(sum( ( cos(q_theta).^2 - 0.5 ) .* power_to_consider )) ;
Q_xy = sum(sum( cos(q_theta) .* sin(q_theta) .* power_to_consider )) ;
%% eigen decomposition of the Q-tensor [Q_xx Q_xy ; Q_xy Q_yy] where Q_yy = - Q_xx
% eigenvalues are \pm \sqrt{ Q_{xx}^2 + Q_{xy}^2 }
%eigVal_large = sqrt(Q_xx ^2 + Q_xy ^2) ;
%eigVal_small = - sqrt(Q_xx ^2 + Q_xy ^2) ;
eigVec_small = [ -Q_xy ; Q_xx+sqrt(Q_xx^2+Q_xy^2) ];
if (Q_xy==0 && Q_xx<0) %arithmetically, this case is equal to eigVec_small == [ 0; Q_xx+sqrt(Q_xx^2)] == [0 ; 0] for Q_xx<0
% but there might be roundoff errors
eigVec_small = [ 1 ; 0 ];
else
eigVec_small = eigVec_small / norm(eigVec_small) ;
end
%eigVec_large = [ eigVec_small(2) ; -eigVec_small(1) ] ;
strength = 2 * sqrt( Q_xx^2 + Q_xy^2 ); % the scalar nematic order parameter
orientation = [eigVec_small(1) -eigVec_small(2)]; % strongest fluctuations are into direction of eigVec_large,
% therefore Directy is orthogonal to that ( = eigVec_small)
% however, this is in the MATLAB coordinate system (x from left to right, y from top to bottom)
% to fit the normal coordinate system (x from left to right, y from bottom to top ) we need to inverse the y-component
%%% no elementwise assignments inside the parfor loop
%Directy(ii,jj,1:2) = orientation;
%Directy(ii,jj,3) = strength;
%Directy(ii,jj,:) = [orientation strength];
RowOfOrientations_1 (jj) = orientation(1);
RowOfOrientations_2 (jj) = orientation(2);
RowOfStrenths (jj) = strength;
end
Directy_1(ii, :) = RowOfOrientations_1;
Directy_2(ii, :) = RowOfOrientations_2;
Directy_3(ii, :) = RowOfStrenths;
end
Directy(:,:,1) = Directy_1;
Directy(:,:,2) = Directy_2;
Directy(:,:,3) = Directy_3;
AnalyseOrientations/ComputeLocalAlignment.m 0 → 100644
View file @ fc349606
function [dominantOrientationField, coherencyField] = ComputeLocalAlignment(theta, magn, boxRadius_LocalAlignment, powerOfMagn_LocalAlignment)
% Compute local degree of alignment (also known as coherency) and locally
% dominant orientations for orientation field given by magnitude magn and
% angle theta. Calculation is within square subpatches of side length
% 1+2*boxRadius, with result assigned to the centeer
%
% Iteration over the subpatches is with parallelisation (parfor)
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% last edit: 2020-10-05
%% preparations
% sanity checks
assert( all(size(theta) == size(magn)) , 'Input array of ComputeLocalAlignment must be of equal size' )
assert(boxRadius_LocalAlignment > 0 && floor(boxRadius_LocalAlignment)==boxRadius_LocalAlignment, 'Input boxRadius to ComputeLocalAlignment must be positive integer')
% set optional arguments if not provided
if (~exist('PowerOfMagn', 'var'))
powerOfMagn_LocalAlignment = 1;
end
% prepare output arrays
i_max = size(theta,1);
j_max = size(theta,2);
dominantOrientationField = zeros( i_max, j_max, 2 );
dominantOrientationField_1 = zeros( i_max, j_max );
dominantOrientationField_2 = zeros( i_max, j_max );
coherencyField = zeros( i_max, j_max );
%% calculation
parfor ii=boxRadius_LocalAlignment+1:i_max-boxRadius_LocalAlignment
dominantOrientationField_1_Row = zeros(1,j_max);
dominantOrientationField_2_Row = zeros(1,j_max);
coherencyField_Row = zeros(1,j_max);
for jj=boxRadius_LocalAlignment+1:j_max-boxRadius_LocalAlignment
%[ dominantOrientationField(ii,jj,:) , coherencyField(ii,jj) ] = ...
[ dominantOrientation , coherency ] = ...
ComputeQTensor( theta(ii-boxRadius_LocalAlignment:ii+boxRadius_LocalAlignment, jj-boxRadius_LocalAlignment:jj+boxRadius_LocalAlignment), ...
magn (ii-boxRadius_LocalAlignment:ii+boxRadius_LocalAlignment, jj-boxRadius_LocalAlignment:jj+boxRadius_LocalAlignment), ...
powerOfMagn_LocalAlignment ) ;
dominantOrientationField_1_Row(jj) = dominantOrientation(1);
dominantOrientationField_2_Row(jj) = dominantOrientation(2);
coherencyField_Row(jj) = coherency;
end
dominantOrientationField_1(ii,:) = dominantOrientationField_1_Row;
dominantOrientationField_2(ii,:) = dominantOrientationField_2_Row;
coherencyField(ii,:) = coherencyField_Row;
end
dominantOrientationField(:,:,1) = dominantOrientationField_1;
dominantOrientationField(:,:,2) = dominantOrientationField_2;
AnalyseOrientations/ComputeQTensor.m 0 → 100644
View file @ fc349606
function [dominantOrientation, coherency] = ComputeQTensor(theta, magn, PowerOfMagn)
% Compute the Q-tensor of a set of orientations, given as magnitude magn
% and angle theta. The function is similar to ComputeDirecty_FourierAndQ.m
% but considers the complete input arrays rather than sliding subpatches of
% defined size.
% The shape of the input does *not* matter, and NaN values are ignored (the
% output is NaN only when all input is NaN)
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% last edit: 2019-04-26
%% use magn to form a distribution (called power_to_consider like in ComputeDirecty_FourierAndQ.m)
power_to_consider = magn .^ PowerOfMagn ;
% turn it into a distribution from which we compute the Q-tensor (= the second moment tensor of the distribution, except for -0.5*unitmatrix)
power_to_consider = power_to_consider / nansum(power_to_consider(:)) ; %normalise
Q_xx = nansum(nansum( ( cos(theta).^2 - 0.5 ) .* power_to_consider )) ;
Q_xy = nansum(nansum( cos(theta) .* sin(theta) .* power_to_consider )) ;
%% eigen decomposition of the Q-tensor [Q_xx Q_xy ; Q_xy Q_yy] where Q_yy = - Q_xx
% eigenvalues are \pm \sqrt{ Q_{xx}^2 + Q_{xy}^2 }
%eigVal_large = sqrt(Q_xx ^2 + Q_xy ^2) ;
%eigVal_small = - sqrt(Q_xx ^2 + Q_xy ^2) ;
eigVec_small = [ -Q_xy ; Q_xx+sqrt(Q_xx^2+Q_xy^2) ];
if (Q_xy==0 && Q_xx<0) %arithmetically, this case is equal to eigVec_small == [ 0; Q_xx+sqrt(Q_xx^2)] == [0 ; 0] for Q_xx<0
% but there might be roundoff errors
eigVec_small = [ 1 ; 0 ];
else
eigVec_small = eigVec_small / norm(eigVec_small) ;
end
%eigVec_large = [ eigVec_small(2) ; -eigVec_small(1) ] ;
strength = 2 * sqrt( Q_xx^2 + Q_xy^2 ); %
dominantOrientation = [ eigVec_small(2) ; -eigVec_small(1) ] ;
%dominantOrientation = eigVec_large;
coherency = strength;
\ No newline at end of file
AnalyseOrientations/LICENSE_distinguishable_colors 0 → 100644
View file @ fc349606
Copyright (c) 2010-2011, Tim Holy
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
AnalyseOrientations/LICENSE_tight_subplot 0 → 100644
View file @ fc349606
Copyright (c) 2016, Pekka Kumpulainen
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
AnalyseOrientations/conv2in2or3D.m 0 → 100644
View file @ fc349606
function convolved = conv2in2or3D(array, kernel, optargs)
% apply 2D convolution to each of the slices of a 3D array
%
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% file generated: 2021-07-08
% last edit: 2021-07-08
% for preallocation, run the convolution once for the first slice
convolved_slice1 = conv2( array(:,:,1), kernel, optargs) ;
NumSlices = size(array,3);
% if there are more than one slice, preallocate result array and perform
% convolution in each of the remaining slices
if NumSlices > 1
convolved = zeros( [size(convolved_slice1) NumSlices] );
convolved(:,:,1) = convolved_slice1;
for kk = 2:NumSlices
convolved(:,:,kk) = conv2( array(:,:,kk), kernel, optargs) ;
end
else
convolved = convolved_slice1;
end
AnalyseOrientations/distinguishable_colors.m 0 → 100644
View file @ fc349606
function colors = distinguishable_colors(n_colors,bg,func)
% DISTINGUISHABLE_COLORS: pick colors that are maximally perceptually distinct
%
% When plotting a set of lines, you may want to distinguish them by color.
% By default, Matlab chooses a small set of colors and cycles among them,
% and so if you have more than a few lines there will be confusion about
% which line is which. To fix this problem, one would want to be able to
% pick a much larger set of distinct colors, where the number of colors
% equals or exceeds the number of lines you want to plot. Because our
% ability to distinguish among colors has limits, one should choose these
% colors to be "maximally perceptually distinguishable."
%
% This function generates a set of colors which are distinguishable
% by reference to the "Lab" color space, which more closely matches
% human color perception than RGB. Given an initial large list of possible
% colors, it iteratively chooses the entry in the list that is farthest (in
% Lab space) from all previously-chosen entries. While this "greedy"
% algorithm does not yield a global maximum, it is simple and efficient.
% Moreover, the sequence of colors is consistent no matter how many you
% request, which facilitates the users' ability to learn the color order
% and avoids major changes in the appearance of plots when adding or
% removing lines.
%
% Syntax:
% colors = distinguishable_colors(n_colors)
% Specify the number of colors you want as a scalar, n_colors. This will
% generate an n_colors-by-3 matrix, each row representing an RGB
% color triple. If you don't precisely know how many you will need in
% advance, there is no harm (other than execution time) in specifying
% slightly more than you think you will need.
%
% colors = distinguishable_colors(n_colors,bg)
% This syntax allows you to specify the background color, to make sure that
% your colors are also distinguishable from the background. Default value
% is white. bg may be specified as an RGB triple or as one of the standard
% "ColorSpec" strings. You can even specify multiple colors:
% bg = {'w','k'}
% or
% bg = [1 1 1; 0 0 0]
% will only produce colors that are distinguishable from both white and
% black.
%
% colors = distinguishable_colors(n_colors,bg,rgb2labfunc)
% By default, distinguishable_colors uses the image processing toolbox's
% color conversion functions makecform and applycform. Alternatively, you
% can supply your own color conversion function.
%
% Example:
% c = distinguishable_colors(25);
% figure
% image(reshape(c,[1 size(c)]))
%
% Example using the file exchange's 'colorspace':
% func = @(x) colorspace('RGB->Lab',x);
% c = distinguishable_colors(25,'w',func);
% Copyright 2010-2011 by Timothy E. Holy
% Parse the inputs
if (nargin < 2)
bg = [1 1 1]; % default white background
else
if iscell(bg)
% User specified a list of colors as a cell aray
bgc = bg;
for i = 1:length(bgc)
bgc{i} = parsecolor(bgc{i});
end
bg = cat(1,bgc{:});
else
% User specified a numeric array of colors (n-by-3)
bg = parsecolor(bg);
end
end
% Generate a sizable number of RGB triples. This represents our space of
% possible choices. By starting in RGB space, we ensure that all of the
% colors can be generated by the monitor.
n_grid = 30; % number of grid divisions along each axis in RGB space
x = linspace(0,1,n_grid);
[R,G,B] = ndgrid(x,x,x);
rgb = [R(:) G(:) B(:)];
if (n_colors > size(rgb,1)/3)
error('You can''t readily distinguish that many colors');
end
% Convert to Lab color space, which more closely represents human
% perception
if (nargin > 2)
lab = func(rgb);
bglab = func(bg);
else
C = makecform('srgb2lab');
lab = applycform(rgb,C);
bglab = applycform(bg,C);
end
% If the user specified multiple background colors, compute distances
% from the candidate colors to the background colors
mindist2 = inf(size(rgb,1),1);
for i = 1:size(bglab,1)-1
dX = bsxfun(@minus,lab,bglab(i,:)); % displacement all colors from bg
dist2 = sum(dX.^2,2); % square distance
mindist2 = min(dist2,mindist2); % dist2 to closest previously-chosen color
end
% Iteratively pick the color that maximizes the distance to the nearest
% already-picked color
colors = zeros(n_colors,3);
lastlab = bglab(end,:); % initialize by making the "previous" color equal to background
for i = 1:n_colors
dX = bsxfun(@minus,lab,lastlab); % displacement of last from all colors on list
dist2 = sum(dX.^2,2); % square distance
mindist2 = min(dist2,mindist2); % dist2 to closest previously-chosen color
[~,index] = max(mindist2); % find the entry farthest from all previously-chosen colors
colors(i,:) = rgb(index,:); % save for output
lastlab = lab(index,:); % prepare for next iteration
end
end
function c = parsecolor(s)
if ischar(s)
c = colorstr2rgb(s);
elseif isnumeric(s) && size(s,2) == 3
c = s;
else
error('MATLAB:InvalidColorSpec','Color specification cannot be parsed.');
end
end
function c = colorstr2rgb(c)
% Convert a color string to an RGB value.
% This is cribbed from Matlab's whitebg function.
% Why don't they make this a stand-alone function?
rgbspec = [1 0 0;0 1 0;0 0 1;1 1 1;0 1 1;1 0 1;1 1 0;0 0 0];
cspec = 'rgbwcmyk';
k = find(cspec==c(1));
if isempty(k)
error('MATLAB:InvalidColorString','Unknown color string.');
end
if k~=3 || length(c)==1,
c = rgbspec(k,:);
elseif length(c)>2,
if strcmpi(c(1:3),'bla')
c = [0 0 0];
elseif strcmpi(c(1:3),'blu')
c = [0 0 1];
else
error('MATLAB:UnknownColorString', 'Unknown color string.');
end
end
end
AnalyseOrientations/downSamplDirecty_ThetaStrength.m 0 → 100644
View file @ fc349606
function [theta_downSampl, strength_downSampl] = downSamplDirecty_ThetaStrength(theta, strength, downSamplFactor, skipRows, skipCols, strength_exponent)
% Downsampling of 2D or 3D orientation field (also known as directionality
% field) in the FIRST TWO dimensions,
% based on averaging in the domain of doubled angle, which is appropriate
% for orientations with their periodicity of pi (180 degree)
% Input and output are of the form (angle, strength).
% [For input in the form (x_component, y_component) a simple convolution
% conv2( <x or y>, 1/(downSamplFactor^2)*ones(downSamplFactor) ) suffices.]
%
% inputs:
% 1.equal-size scalar fields theta and strength give the orientation angle,
% and a strength of the orientation. Strength must be nonnegative.
% strength 0 can be used to mask out some data.
% 2.downSamplFactor - blocks of square size downSamplFactor are used
% 3.optional scalar values skipRows and skipCols (default: 0). The first
% block starts after this number of rows (skipRows) and columns
% (skipCols).
% 4.optional scalar value strength_exponent controls the weight of strength
% for the averaging within each block. Typically either 1 or 2 [from what
% is found in the literature]. Default is 1.
%
% Parts of the input array, that do not fill a complete block at the last
% rows and/or columns, are dropped. Hence the size of the output array is
% floor( (size(input arrays)-[skipRows skipCols] )/downSamplFactor )
%
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% generated: 2021-03-26
% last edit: 2022-08-25
% CHANGELOG
% 2021-07-08 include conv2in2or3D to allow 3d arrays to be processed
%% preparations
% sanity checks
assert(all(size(theta)==size(strength)), 'Input fields "theta" and "strength" must have equal shape')
assert(int8(downSamplFactor) == downSamplFactor, 'Input "downSamplFactor" must be integer')
% set optional arguments if not provided
if (~exist('skipRows', 'var'))
skipRows = 0;
else
assert(int8(skipRows) == skipRows && skipRows > 0, 'Input "skipRows" must be non-negative integer')
end
if (~exist('skipCols', 'var'))
skipCols = 0;
else
assert(int8(skipCols) == skipCols && skipCols > 0, 'Input "skipCols" must be non-negative integer')
end
if (~exist('strength_exponent', 'var'))
strength_exponent = 1;
end
%% calculations
% convert to u and v (Cartesian coordinates)
prefactor = strength.^strength_exponent;
u_doubleAngle = prefactor .* cos( 2*theta );
v_doubleAngle = prefactor .* sin( 2*theta );
u_doubleAngle_blocks = 1/(downSamplFactor^2) ...
* conv2in2or3D(u_doubleAngle, ones(downSamplFactor), 'valid') ;
v_doubleAngle_blocks = 1/(downSamplFactor^2) ...
* conv2in2or3D(v_doubleAngle, ones(downSamplFactor), 'valid') ;
u_doubleAngle_downSampl = u_doubleAngle_blocks(1+skipRows:downSamplFactor:end, 1+skipCols:downSamplFactor:end, :);
v_doubleAngle_downSampl = v_doubleAngle_blocks(1+skipRows:downSamplFactor:end, 1+skipCols:downSamplFactor:end, :);
% transform back into strength and angle theta
strength_downSampl = (u_doubleAngle_downSampl .^2 + v_doubleAngle_downSampl .^2) .^0.5 ;
theta_downSampl = 0.5*atan2( v_doubleAngle_downSampl, u_doubleAngle_downSampl );
AnalyseOrientations/loadGrayStack.m 0 → 100644
View file @ fc349606
function [InfoOnImage, ImageStack] = loadGrayStack(tifPath)
% function for loading 2D or 3D grayscale images from tif files
% argument tifPath = the path of the image stack .tif file, in string form.
% Example use:
% [InfoOnImage, ImageStack] = loadGrayStack('/Users/khoffman/Documents/ImageStack.tif')
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% last edit: 2018-03-29
InfoOnImage=imfinfo(tifPath);
mImage=InfoOnImage(1).Height;
nImage=InfoOnImage(1).Width;
NumberImages=length(InfoOnImage);
ImageStack=zeros(mImage,nImage,NumberImages,'double');
for i=1:NumberImages
ImageStack(:,:,i)=imread(tifPath,'Index',i,'Info',InfoOnImage);
end
end
% Codes modified from https://blogs.mathworks.com/steve/2009/04/02/matlab-r2009a-imread-and-multipage-tiffs/
% Consider also http://www.matlabtips.com/how-to-load-tiff-stacks-fast-really-fast/
AnalyseOrientations/smoothDirecty_ThetaStrength.m 0 → 100644
View file @ fc349606
function [theta_smoothed, strength_smoothed] = smoothDirecty_ThetaStrength(theta, strength, sigma, strength_exponent, padding)
% Smoothing of orientation fields (also known as directionality fields)
% in the domain of doubled angle, which is appropriate for orientations
% with their periodicity of pi (180 degree)
% Generally, any smoothing would work. Here, a Gaussian is implemented.
% Input and output are of the form (angle, strength).
% For input in the form (x_component, y_component) use smoothDirecty_uANDv
%
% inputs:
% 1.equal-size scalar fields theta and strength give the orientation angle,
% and a strength of the orientation. Strength must be nonnegative.
% strength 0 can be used to mask out some data.
% 2.non-negative sigma = standard deviation of a Gaussian to be used
% for smoothing. Scalar for isotropic blur or vector of 2 (or 3) nonnegative
% entries for 2- (or 3-) dimensional input fields
% 3.optional scalar value strength_exponent. Typically either 1 or 2 [from
% what is found in the literature]
% Default is 1
% 4.optional string padding. This is passed on to the Gaussian blur.
% Default is 'replicate'
%
% For each pixel, the input is converted to vector representing the doubled angle as
% strength.^strength_exponent .* cos( 2*theta )
% strength.^strength_exponent .* sin( 2*theta )
% which is smoothed with a Gaussian of standard deviation sigma, and
% converted back to orientation angle theta_smoothed and strength_smoothed
%
% Karl Hoffmann, Max Planck Institute of Molecular Cell Biology and
% Genetics, Dresden, Germany
% last edit: 2018-11-19
%% preparations
% sanity checks
assert(all(size(theta)==size(strength)), 'Input fields "theta" and "strength" must have equal shape')
assert(isscalar(sigma) || size(size(theta),2) == size(size(sigma),2) ,'Input "sigma" must be scalar or fit the dimensionality of the fields')
assert(min(sigma) >= 0.0, 'Input "sigma" must be non-negative' )
% set optional arguments if not provided
if (~exist('padding', 'var'))
padding = 'replicate';
end
if (~exist('strength_exponent', 'var'))
strength_exponent = 1;
end
%% calculations
if min(sigma) > 0
prefactor = strength.^strength_exponent;
u_doubleAngle = prefactor .* cos( 2*theta );
v_doubleAngle = prefactor .* sin( 2*theta );
%smoothing itself (other than Gaussian would work as well)
if size(size(theta),2) == 2 % two-dimensional field
u_doubleAngle_smoothed = imgaussfilt(u_doubleAngle, sigma, 'FilterDomain','spatial', 'Padding', padding );
v_doubleAngle_smoothed = imgaussfilt(v_doubleAngle, sigma, 'FilterDomain','spatial', 'Padding', padding );
else
u_doubleAngle_smoothed = imgaussfilt3(u_doubleAngle, sigma, 'FilterDomain','spatial', 'Padding', padding );
v_doubleAngle_smoothed = imgaussfilt3(v_doubleAngle, sigma, 'FilterDomain','spatial', 'Padding', padding );
end
% transform back into strength and angle theta
strength_smoothed = (u_doubleAngle_smoothed .^2 + v_doubleAngle_smoothed .^2) .^0.5 ;
theta_smoothed = 0.5*atan2( v_doubleAngle_smoothed, u_doubleAngle_smoothed );
else % no smoothing
strength_smoothed = strength;
theta_smoothed = theta;
end
end
AnalyseOrientations/tight_subplot.m 0 → 100644
View file @ fc349606
function [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
% tight_subplot creates "subplot" axes with adjustable gaps and margins
%
% [ha, pos] = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
%
% in: Nh number of axes in hight (vertical direction)
% Nw number of axes in width (horizontaldirection)
% gap gaps between the axes in normalized units (0...1)
% or [gap_h gap_w] for different gaps in height and width
% marg_h margins in height in normalized units (0...1)
% or [lower upper] for different lower and upper margins
% marg_w margins in width in normalized units (0...1)
% or [left right] for different left and right margins
%
% out: ha array of handles of the axes objects
% starting from upper left corner, going row-wise as in
% subplot
% pos positions of the axes objects
%
% Example: ha = tight_subplot(3,2,[.01 .03],[.1 .01],[.01 .01])
% for ii = 1:6; axes(ha(ii)); plot(randn(10,ii)); end
% set(ha(1:4),'XTickLabel',''); set(ha,'YTickLabel','')
% Pekka Kumpulainen 21.5.2012 @tut.fi
% Tampere University of Technology / Automation Science and Engineering
if nargin<3; gap = .02; end
if nargin<4 || isempty(marg_h); marg_h = .05; end
if nargin<5; marg_w = .05; end
if numel(gap)==1;
gap = [gap gap];
end
if numel(marg_w)==1;
marg_w = [marg_w marg_w];
end
if numel(marg_h)==1;
marg_h = [marg_h marg_h];
end
axh = (1-sum(marg_h)-(Nh-1)*gap(1))/Nh;
axw = (1-sum(marg_w)-(Nw-1)*gap(2))/Nw;
py = 1-marg_h(2)-axh;
% ha = zeros(Nh*Nw,1);
ii = 0;
for ih = 1:Nh
px = marg_w(1);
for ix = 1:Nw
ii = ii+1;
ha(ii) = axes('Units','normalized', ...
'Position',[px py axw axh], ...
'XTickLabel','', ...
'YTickLabel','');
px = px+axw+gap(2);
end
py = py-axh-gap(1);
end
if nargout > 1
pos = get(ha,'Position');
end
ha = ha(:);
LICENSE
View file @ fc349606
BSD 3-Clause License
Copyright (c) 2021, Scientific Computing
Copyright (c) 2021, Scientific Computing, and 2020-2022, Karl Hoffmann
All rights reserved.
Redistribution and use in source and binary forms, with or without
......
README.md
View file @ fc349606
......@@ -4,11 +4,19 @@ This repository contains the scripts and documentation associated to the publica
**Matrix topology guides collective cell migration in vivo**
Karen Grace Soans, Ana Patricia Ramos, Jaydeep Sidhaye, Abhijeet Krishna, Anastasia Solomatina, Karl Hoffmann, Raimund Schlüßler, Jochen Guck, Ivo. F Sbalzarini, Carl Modes, Caren Norden
## Script
## Scripts
**`extractDerotatedSurface_v1.0_threshHuang.py`**
Goal: Extract a projection of a tilted surface. This Fiji script extracts a projection of a tilted surface of derotated selected regions of the image (optical cup) to align the surface as much as possible with the xy plane. The surface is then extracted (height map + reslicing) and max-projected.
**`AnalyseOrientations/AnalyseOpticalCupLamininOrientation.m`**
Goal: Estimate locally dominant orientation of fibrous structures in two-dimensional images, and calculate their coherency (i.e. degree of alignment) both locally in a sliding window and globally across the foreground in each region of interest.
The script loops over two-dimensional images organised in a folder structure `head_folder/<Stage>/<Embryo>/<RegionOfTheEye>/<DerotationConditions>/Roi<nn>_ZMAX_of_Resliced_of_Derotated_with2ndCrop.tif` (with each single image obtained by `extractDerotatedSurface_v1.0_threshHuang.py` for the publication), and collects statistics across them.
Note: Code related to the spring network also used for the publication is available from https://git.mpi-cbg.de/krishna/laminin_network.
## Usage & Disclaimer
The scripts are provided as is, without the intention that an interested reader will be able to run them directly. The authors share their scripts here, not a ready-to-use tool. Still, they should provide enough technical detail to allow replicating the analyses.
......@@ -16,7 +24,7 @@ The scripts are provided as is, without the intention that an interested reader
The scripts are provided as is under [BSD3 License](https://opensource.org/licenses/BSD-3-Clause).
## Support
For further support please contact [Norden Lab](https://gulbenkian.pt/ciencia/research-groups/cnorden/) (mail: cnorden@igc.gulbenkian.pt) or Scientific Computing Facility (mail: scicomp@mpi-cbg.de).
For further support please contact [Norden Lab](https://gulbenkian.pt/ciencia/research-groups/cnorden/) (mail: cnorden@igc.gulbenkian.pt), Karl Hoffmann (mail: karlhoff@mpi-cbg.de) or Scientific Computing Facility (mail: scicomp@mpi-cbg.de).
## DOI
This repository is resolvable via: https://dx.doi.org/21.11101/0000-0007-F43E-1 with the [GWDG](https://www.gwdg.de/en/application-services/persistent-identifier-pid/epic-pid-search) assigned PID: 21.11101/0000-0007-F43E-1