From 7959824d420d70fbf0beba7b98a743fba9bbad97 Mon Sep 17 00:00:00 2001
From: Karl Hoffmann <karlhoff@mpi-cbg.de, karl@hoffmann-freiberg.de>
Date: Mon, 5 Sep 2022 13:22:56 +0200
Subject: [PATCH] Add scripts for orientation analysis
---
.../AnalyseOpticalCupLamininOrientation.m | 1015 +++++++++++++++++
AnalyseOrientations/ComputeDirecty.m | 336 ++++++
AnalyseOrientations/ComputeDirectyGradient.m | 118 ++
.../ComputeDirecty_FourierAndQ.m | 163 +++
AnalyseOrientations/ComputeLocalAlignment.m | 52 +
AnalyseOrientations/ComputeQTensor.m | 41 +
.../LICENSE_distinguishable_colors | 24 +
AnalyseOrientations/LICENSE_tight_subplot | 24 +
AnalyseOrientations/conv2in2or3D.m | 25 +
AnalyseOrientations/distinguishable_colors.m | 152 +++
.../downSamplDirecty_ThetaStrength.m | 74 ++
AnalyseOrientations/loadGrayStack.m | 23 +
.../smoothDirecty_ThetaStrength.m | 75 ++
AnalyseOrientations/tight_subplot.m | 66 ++
LICENSE | 2 +-
README.md | 13 +-
16 files changed, 2200 insertions(+), 3 deletions(-)
create mode 100644 AnalyseOrientations/AnalyseOpticalCupLamininOrientation.m
create mode 100644 AnalyseOrientations/ComputeDirecty.m
create mode 100644 AnalyseOrientations/ComputeDirectyGradient.m
create mode 100644 AnalyseOrientations/ComputeDirecty_FourierAndQ.m
create mode 100644 AnalyseOrientations/ComputeLocalAlignment.m
create mode 100644 AnalyseOrientations/ComputeQTensor.m
create mode 100644 AnalyseOrientations/LICENSE_distinguishable_colors
create mode 100644 AnalyseOrientations/LICENSE_tight_subplot
create mode 100644 AnalyseOrientations/conv2in2or3D.m
create mode 100644 AnalyseOrientations/distinguishable_colors.m
create mode 100644 AnalyseOrientations/downSamplDirecty_ThetaStrength.m
create mode 100644 AnalyseOrientations/loadGrayStack.m
create mode 100644 AnalyseOrientations/smoothDirecty_ThetaStrength.m
create mode 100644 AnalyseOrientations/tight_subplot.m
diff --git a/AnalyseOrientations/AnalyseOpticalCupLamininOrientation.m b/AnalyseOrientations/AnalyseOpticalCupLamininOrientation.m
new file mode 100644
index 0000000..116a00e
--- /dev/null
+++ b/AnalyseOrientations/AnalyseOpticalCupLamininOrientation.m
@@ -0,0 +1,1015 @@
+% 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
diff --git a/AnalyseOrientations/ComputeDirecty.m b/AnalyseOrientations/ComputeDirecty.m
new file mode 100644
index 0000000..ecc7959
--- /dev/null
+++ b/AnalyseOrientations/ComputeDirecty.m
@@ -0,0 +1,336 @@
+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;
+
+
diff --git a/AnalyseOrientations/ComputeDirectyGradient.m b/AnalyseOrientations/ComputeDirectyGradient.m
new file mode 100644
index 0000000..0251f14
--- /dev/null
+++ b/AnalyseOrientations/ComputeDirectyGradient.m
@@ -0,0 +1,118 @@
+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
+
diff --git a/AnalyseOrientations/ComputeDirecty_FourierAndQ.m b/AnalyseOrientations/ComputeDirecty_FourierAndQ.m
new file mode 100644
index 0000000..4e18c2d
--- /dev/null
+++ b/AnalyseOrientations/ComputeDirecty_FourierAndQ.m
@@ -0,0 +1,163 @@
+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;
+
+
+
diff --git a/AnalyseOrientations/ComputeLocalAlignment.m b/AnalyseOrientations/ComputeLocalAlignment.m
new file mode 100644
index 0000000..4a5db4f
--- /dev/null
+++ b/AnalyseOrientations/ComputeLocalAlignment.m
@@ -0,0 +1,52 @@
+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;
diff --git a/AnalyseOrientations/ComputeQTensor.m b/AnalyseOrientations/ComputeQTensor.m
new file mode 100644
index 0000000..ac41122
--- /dev/null
+++ b/AnalyseOrientations/ComputeQTensor.m
@@ -0,0 +1,41 @@
+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
diff --git a/AnalyseOrientations/LICENSE_distinguishable_colors b/AnalyseOrientations/LICENSE_distinguishable_colors
new file mode 100644
index 0000000..63168b1
--- /dev/null
+++ b/AnalyseOrientations/LICENSE_distinguishable_colors
@@ -0,0 +1,24 @@
+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.
diff --git a/AnalyseOrientations/LICENSE_tight_subplot b/AnalyseOrientations/LICENSE_tight_subplot
new file mode 100644
index 0000000..342178c
--- /dev/null
+++ b/AnalyseOrientations/LICENSE_tight_subplot
@@ -0,0 +1,24 @@
+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.
diff --git a/AnalyseOrientations/conv2in2or3D.m b/AnalyseOrientations/conv2in2or3D.m
new file mode 100644
index 0000000..3bab248
--- /dev/null
+++ b/AnalyseOrientations/conv2in2or3D.m
@@ -0,0 +1,25 @@
+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
diff --git a/AnalyseOrientations/distinguishable_colors.m b/AnalyseOrientations/distinguishable_colors.m
new file mode 100644
index 0000000..f9e9e25
--- /dev/null
+++ b/AnalyseOrientations/distinguishable_colors.m
@@ -0,0 +1,152 @@
+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
diff --git a/AnalyseOrientations/downSamplDirecty_ThetaStrength.m b/AnalyseOrientations/downSamplDirecty_ThetaStrength.m
new file mode 100644
index 0000000..08ef5fa
--- /dev/null
+++ b/AnalyseOrientations/downSamplDirecty_ThetaStrength.m
@@ -0,0 +1,74 @@
+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 );
diff --git a/AnalyseOrientations/loadGrayStack.m b/AnalyseOrientations/loadGrayStack.m
new file mode 100644
index 0000000..42d33d6
--- /dev/null
+++ b/AnalyseOrientations/loadGrayStack.m
@@ -0,0 +1,23 @@
+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/
diff --git a/AnalyseOrientations/smoothDirecty_ThetaStrength.m b/AnalyseOrientations/smoothDirecty_ThetaStrength.m
new file mode 100644
index 0000000..856b98e
--- /dev/null
+++ b/AnalyseOrientations/smoothDirecty_ThetaStrength.m
@@ -0,0 +1,75 @@
+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
+
+
+
diff --git a/AnalyseOrientations/tight_subplot.m b/AnalyseOrientations/tight_subplot.m
new file mode 100644
index 0000000..73e0d1f
--- /dev/null
+++ b/AnalyseOrientations/tight_subplot.m
@@ -0,0 +1,66 @@
+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(:);
diff --git a/LICENSE b/LICENSE
index 4e28ecd..67203fa 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
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
diff --git a/README.md b/README.md
index a4bee1f..d03d6fd 100644
--- a/README.md
+++ b/README.md
@@ -4,11 +4,20 @@ 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.
+
+`Code related to the spring network`
+
+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 +25,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
--
GitLab