AnalyseOpticalCupLamininOrientation.m

% 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);