% 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