%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%  [doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
%\documentclass{bmcart}

%%% Load packages
%\usepackage{amsthm,amsmath}
%\RequirePackage{natbib}
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
%\RequirePackage{hyperref}
\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails

\usepackage{bbm}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{multirow}
\usepackage[normalem]{ulem}
\usepackage{mathtools, cuted} % for long formulaes
\usepackage{url}

\newcommand{\R}{{\mathbbm{R}}}
\newcommand{\shape}{{\mathcal{S}}}
\newcommand{\shapespace}{{\textbf{S}}}
\newcommand{\object}{{\mathcal{O}}}
\newcommand{\W}{{\mathcal{W}}}
\renewcommand{\d}{{\mathrm{d}}}
\newcommand{\cov}{{\textbf{Cov}}}
\newcommand{\maha}{{d}}
\newcommand{\oldArticle}[1]{}
\newcommand{\revSP}[1]{{\color{green}{\textbf{SP: }#1}}}
\newcommand{\revBH}[1]{{\color{green}{\textbf{BH: }#1}}}
\newcommand{\revMR}[1]{{\color{green}{\textbf{MR: }#1}}}
\newcommand{\revBW}[1]{{\color{green}{\textbf{BW: }#1}}}
\newcommand{\done}{{\color{red}\textbf{DONE!}}}

\newcommand{\comment}[1]{{\color{red}{[\textbf{comment: }#1]}}}
\newcommand{\missing}[1]{{\color{red}\textbf{MISSING:} {#1}}}

\newcommand{\myCaption}[1]{\caption{}}

\DeclareRobustCommand{\citeext}[1]{\citeauthor{#1}~\cite{#1}}

\useunder{\uline}{\ul}{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\def\includegraphic{}
%\def\includegraphics{}



%%% Put your definitions there:
\startlocaldefs
\endlocaldefs


%%% Begin ...
\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Statistical shape analysis of tap roots: \\ a methodological case study on laser scanned sugar beets}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author[
addressref={aff1},                   % id's of addresses, e.g. {aff1,aff2}
corref={aff1},                       % id of corresponding address, if any
noteref={n1},                        % id's of article notes, if any
email={behrend.heeren@ins.uni-bonn.de}   % email address
]{\inits{BH}\fnm{Behrend} \snm{Heeren}}

\author[
   addressref={aff2},
	% corref={aff2},
	 noteref={n1}, 
   email={paulus@ifz-goettingen.de}
]{\inits{SP}\fnm{Stefan} \snm{Paulus}}


\author[
   addressref={aff3},
   email={h.goldbach@uni-bonn.de}
]{\inits{HG}\fnm{Heiner} \snm{Goldbach}}

\author[
addressref={aff4},
email={heiner.kuhlmann@uni-bonn.de}
]{\inits{HK}\fnm{Heiner} \snm{Kuhlmann}}

\author[
addressref={aff2},
email={mahlein@ifz-goettingen.de}
]{\inits{AKM}\fnm{Anne-Katrin} \snm{Mahlein}}

\author[
   addressref={aff1},
   email={Martin.rumpf@ins.uni-bonn.de}
]{\inits{MR}\fnm{Martin} \snm{Rumpf}}

\author[
addressref={aff5},
email={Benedikt.Wirth@uni-muenster.de}
]{\inits{BW}\fnm{Benedikt} \snm{Wirth}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\address[id=aff1]{%                           % unique id
  \orgname{Institute for Numerical Simulation, University of Bonn}, % university, etc
  \street{Endenicher Allee 60},                     %
  \postcode{53115}                                % post or zip code
  \city{Bonn},                              % city
  \cny{Germany}                                    % country
}

\address[id=aff2]{%                           % unique id
  \orgname{Institute of Sugar Beet Research, Germany}, % university, etc
  \street{Holtenser Landstr. 77},                     %
  \postcode{37079}                                % post or zip code
  \city{G\"{o}ttingen},                              % city
  \cny{Germany}                                    % country
}


\address[id=aff3]{%                           % unique id
  \orgname{INRES Plant Nutrition, University of Bonn}, % university, etc
  \street{Karlrobert-Kreiten-Strasse 13},                     %
  \postcode{53115}                                % post or zip code
  \city{Bonn},                              % city
  \cny{Germany}                                    % country
}

\address[id=aff4]{%                           % unique id
  \orgname{Institute for Geodesy and Geoinformation,  University of Bonn}, % university, etc
  \street{Nussallee 17},                     %
  \postcode{53115}                                % post or zip code
  \city{Bonn},                              % city
  \cny{Germany}                                    % country
}

\address[id=aff5]{%                           % unique id
  \orgname{Institute for Analysis and Numerics, University of Münster}, % university, etc
  \street{Einsteinstr. 62},                     %
  \postcode{48149}                                % post or zip code
  \city{Münster},                              % city
  \cny{Germany}                                    % country
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{artnotes}
%\note{Sample of title note}     % note to the article
\note[id=n1]{Equal contributor} % note, connected to author
\end{artnotes}

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract} % abstract

\parttitle{Background} %if any
The efficient and robust statistical analysis of the shape of plant organs of different cultivars is an important investigation issue in plant breeding and enables a robust cultivar description within the breeding progress. 
Laserscanning is a highly accurate and high resolution technique to acquire the 3D shape of plant surfaces. 
The computation of a shape based principal component analysis (PCA)
built on concepts from continuum mechanics has proven to be an effective tool for a qualitative and quantitative shape examination.

\parttitle{Results} 
The shape based PCA was used for a statistical analysis of 140 sugar beet roots of different cultivars. 
The calculation of the mean sugar beet root shape and the description of the main variations was possible. 
Furthermore, unknown and individual tap roots could be attributed to their cultivar by means of a robust classification tool based on the PCA results. % (classification), based on characteristic tap root shapes.

\parttitle{Conclusion} 
The method demonstrates that it is possible to identify principal modes of root shape variations automatically and to quantify associated variances out of laserscanned 3D sugar beet tap root models. 
The introduced approach is not limited to the 3D shape description by laser scanning. A transfer to 3D MRI or radar data is also conceivable.

\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The keywords begin here                  %%
%%                                          %%
%% Put each keyword in separate \kwd{}.     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{keyword}
\kwd{tap roots}
\kwd{shape variation}
\kwd{energy minimization}
\kwd{3D point clouds}
\end{keyword}

% MSC classifications codes, if any
%\begin{keyword}[class=AMS]
%\kwd[Primary ]{}
%\kwd{}
%\kwd[; secondary ]{}
%\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%% start of article main body
% <put your article body there>

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Background %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Material and Methods %%
\section{Background}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PLANTS BACKGROUND AND RELATED WORK
In breeding and precision agriculture there is a need for a precise description of the 3D architecture of a crop, a plant organ or a harvested product, in a fast and reproducible way \cite{Bradshaw2010}. 

Potential applications involve automated selection procedures in plant breeding of phenotypes with the desirable features and traits (like grain or root shape), 
assessment of crop development during the growth period, enabling an optimised crop management and sorting out special forms or non-crop contaminants (stones, weed seeds) of harvested products. 
An additional application is the adaption of harvester settings, based on archetypal 3D-geometries of harvested cultivars to avoid mechanical damage and by this quality reduction.

%%% PART ANNE
Yield formation is the combined result of genetic characteristics and environmental effects. In sugar beet the annual increase in sugar yield accounts 1.5\% \cite{Maerlaender2018}. 
This is mainly due to an increase in the root yield, while the sugar content of varieties remained stable \cite{Hoffmann2018}. 
One of the main targets of sugar beet breeding is the root development, including relevant quality parameters such as sugar content, non-sugar compounds or mark-content. 
It can be observed that the dynamic of storage root development in sugar beet has no specific growth stages. 
There is no phase of maturation, and therefore, yield development and potential can be estimated depending on the duration of the growing period \cite{Hoffmann2018}.  
Unfortunately, in sugar beet development of leaf biomass is not correlated to storage root biomass or yield formation. 
However, the shape of the sugar beet tap root plays a significant role for the entire processing chain. 
%%% PART ANNE

Besides effects of varieties, it is well known that the shape of tap roots mainly depends on the type and status of soil, management and environmental conditions \cite{Vamerali2009,Broom1997}. A strong correlation exists between basic root shape parameters like area, length, or radial variation and sugar yield and quality. %\cite{Tsialtas2010}.
Additionally, size and shape is characteristic for cultivars. 
Factors such as formation of branches or soil-adhesion are relevant breeding traits. 
On-field devices are proposed to measure and monitor basic parameters of tap roots \cite{Shimazu2010} with varying degrees of success. 
Shape models are used for the classification of crops and weeds \cite{Persson2008}. 
Image processing methodology already plays a significant role for the observation of plant growth \cite{Midtiby2012}.
Computer vision is also used for the analysis of plant root shapes. However, currently these methods are mostly based on 2D imaging modalities \cite{Moreda2012}, whereas the most precise description would result from a robust statistical analysis of the true 3D description of tap roots. 
First results in imaging the 3D shape of sugar beet have been used to extract scalar parameters of the tap root \cite{paulus2014b} such as height, width, volume and surface area.
By imaging over time the development and growth can be observed. The resulting crop growth model can be improved by a true 4D description of the crop plants \cite{BeFlHe13}.
Analysis of 3D point clouds with recent mathematical models or machine learning approaches further improve the efficiency and biological interpretability of plant sensor data \cite{Wahabzada2015}. 
This helps to assess the genotype-phenotype-environment interactions and to dissect important traits \cite{Tardieu2010} \cite{Mahlein2018}.\\[1ex]

%%%% PART TAKEN FROM DISCUSSION
Agricultural crops have to cope with adverse environments but still have to maintain productivity at the highest possible level. Here, breeding has to identify and make use of the best combinations of traits \cite{Chahal2002,Acquaah2007}. 
These traits are mostly quantitative and inherited in a complex manner \cite{Tardieu2010}. Although morphological changes  are downstream effects of altered gene expression, metabolic adaptations and environmental changes, 
their precise, unequivocal, and unbiased identification is of utmost importance in a breeding process to assess Genotype $\times$ Environment $\times$ Phenotype (G$\times$E$\times$P) interactions \cite{Iwata2002}. 
So far, breeding involves classification of desirable traits (usually around 20) by using empirical scales, which can be assessed mostly by visual assessment in traditional breeding schemes \cite{Chahal2002}. 
Considering that there are easily 10,000 individuals or more that have to be classified within a very short time span, the common method is prone to human bias (tiredness, adverse and rapidly changing light conditions). 
The same holds true for the harvested produce where shapes e.\,g.\ of tuber crops have to be assessed. Unequivocal identification of shapes or outlines is of high importance as well for purposes of precision farming where there is a need for an on-the-flight identification of shapes and their assessment, 
e.\,g.\ for the differentiation between crop and weed plants, diseased and healthy plants or plant organs or species identification (see e.\,g.\ \cite{Neto2006}). 
Further refinement of phenotyping and the increase of traits that can be simultaneously observed calls for objective and automated ways of trait identification. 
There have been various attempts so far to describe leaf or tap root shapes with elliptic Fourier analysis \cite{Neto2006,Iwata2002b,Iwata2002}. 
The challenge is now to analyze true 3D shape variations into main components, while most methods so far rely on 2D image analysis \cite{Iwata2002,Tsialtas2010}. \\[1ex]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SHAPE STATISTICS BACKGROUND AND RELATED WORK
Knowledge of the statistical variation of root shapes has various applications. 
In plant breeding such information can be used to select tendencies with favorable traits, such as an even shape, large crop size, 
or an optimal form for the use of crop harvesters and to minimize breaking root tips, soil tare and mechanical damages. 
In crop management the statistical data may help to detect deviations from the expected plant development, which can be used in a feedback mechanism to adjust growth control parameters, such as fertilizer or pesticides.

In this study we apply a computer vision tool for statistical analysis that is purely shape based. 
%Statistical models of the shape or appearance of a class of objects are widely used in computer vision and graphics to model the variability over the object class. 
%They can be used to constrain synthesis or analysis problems and their parameter space provides a compact representation which can be used for classification or intuitive editing. 
%Human faces and bodies have proven particularly amenable to statistical modelling. In this case, the source of shape variability is usually changes in identity or dynamic deformation.
Statistical models of shape have been used widely in computer vision and graphics \cite{brunton2014review}. 
In a 2D setting, PCA-based models such as Active Shape \cite{cootes1995active} 
or Appearance Models \cite{cootes2001active} provide a parametric representation of shape that can be used for segmentation, tracking and recognition. 
In a 3D setting, they are typically used for fitting to noisy or ambiguous data or for 3D reconstruction via analysis-by-synthesis. 
Essentially, the statistical model provides a constraint that significantly reduces the parameter space for many shape processing problems. 
%Such models begin by transforming raw data into a shape space. 
%This requires establishing a dense point-to-point correspondence between the training samples. 
%This problem has previously been addressed using non-rigid ICP \cite{amberg2007optimal}, template fitting \cite{salazar2014fully}, optical flow \cite{blanz2003face} and minimum description length \cite{davies20023d}.

It is a wide spread assumption that the input data is collected from a shape space that is considered as a Riemannian manifold. 
The classical treatment of shape space is due to Kendall \cite{kendall1984shape}, 
in which sets of landmarks are considered points on a shape manifold in which the effects of scale, 
rotation and translation are factored out. 
The tangent plane to Kendall's shape space enables linear principle component analysis (PCA) in which Euclidean distance approximates Procrustes distance. 
Srivastava et al.~\cite{srivastava2011shape} propose a representation for analysing shapes of curves under an elastic metric. 
%This allows them to define geodesics and geodesic distance between curves in a way that is invariant to transformations including reparameterisation.
Killian et al.~\cite{kilian2007geometric} model the space of triangulated shapes.
%They use a Riemannian metric which measures the stretching of triangle edges in order to compute geodesics. 
Statistical analysis can be performed on shape spaces in a manner that respects the Riemannian geometry of the manifold. 
This requires Riemannian notions of concepts such as distance, mean value and covariance \cite{pennec2006intrinsic}.
%Pennec \cite{pennec2006intrinsic} showed how to compute these measures for a number of geometric primitives that do not form vector spaces. 
Based on these quantities Fletcher et al.~\cite{FlLuPi04} transferred the concept of PCA to manifolds by considering a principal geodesic analysis (PGA). 
The idea of PGA was used by Tournier et al.~\cite{tournier2009motion} to build a statistical skeleton model and by
Heeren et al.~\cite{HeRuSmZh18} to perform a statistical analysis in the space of triangle meshes. 
%Tycowicz et al. \cite{TyAmMu18} presented a non-Euclidean statistical analysis of triangle meshes. % and consider medical applications such as shape-based classification of morphological disorders. Besides modelling shape variation, PGA has also been applied to distributions of probability measures using an optimal transport metric \cite{seguy2015principal}.}
However, all these approaches assume the underlying shape space to be a Riemannian manifold. 
In particular, the dissimilarity of two shapes is quantified by the length of an optimal, connecting curve. 
%where the underlying metric is frequently measuring the lack of isometry, 
%e.g. via a (linearised) elastic energy, of a shape deformation path. 
In contrast, we here consider a purely elastic model that measures shape dissimilarity by the amount of elastic deformation energy \cite{Ci97}---for details on the physical and geometrical differences we refer to \cite{RuWi10}. 
Using the physical model of 3D elasticity Rumpf and Wirth compute shape averages \cite{RuWi08} and
describe a covariance analysis \cite{RuWi09,RuWi09b} of shapes represented as boundary contours of elastic objects. 
Here the elastic average is defined as the shape that minimizes the sum of elastic deformation energies to all the given input shapes. 
With an average at hand a classical PCA is then applied to the displacements with respect to the input shapes. 
This elastic shape analysis was performed on the space of triangle meshes in \cite{ZhHeRuSm15} and 
will be applied in the present work to the shape space formed by sugar beet tap roots.\\[1ex]

%Frequently, the underlying metric is based on measuring the lack of isometry, 
%e.g. via a (linearised) elastic energy acting on the Cauchy-Green strain tensor of an associated
%infinitesimal deformation.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% DATA BACKGROUND AND RELATED WORK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In our experiment the acquisition of input shapes is realized based on 3D laser scanning---a method to recover 3D point clouds from objects. 
This method is well established in the agricultural context and has been used on plant level for 3D modelling of the canopy of tomato plants \cite{hosoi}, 
for in-field scanning of pear-trees \cite{palacin} and imaging physiological responses of leaves \cite{Omasa}. 
More detailed scans for organ specific parameterization were also possible when using close-up laserscanning. 
This enables e.g. identification of single organs \cite{Paulus2013, paulus2014b, Wahabzada2015} or tracking of growth on organ level  \cite {paulus2014c}, and it has shown to be very accurate \cite{paulus2014a}. 
Close-up scanning enables high resolution and high accuracy imaging with point to point distances below a millimeter \cite{Wagner2011}. 
This enables highly accurate 3D surface models of various plant types.\\[1ex]

 

%%%% CONTRIUBTION SUMMARIZED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Our contribution.} 
In this work, we demonstrate the extraction and mathematical description of characteristic shape features in the tap roots of different cultivars of sugar beets. 
To this end, a huge data set consisting of detailed 3D descriptions of beet root samples 
is reduced to a small number of important parameters without loosing relevant biological information. 
Our main contributions are twofold. 
First, we perform a shape based analysis of laser scanned sugar beet tap roots. 
In particular, one can compute a robust and reliable shape mean 
as well as principal modes of shape variation on large ensembles of sugar beets with a substantial variability in the shape geometry.  
In particular, the statistical tool is not based on predefined quantitative properties (such as length or volume), and it is invariant under rigid body motions.
Second, we propose an automatic classification tool based on results from the statistical analysis. 
The number of tap roots that were classified correctly without having been in the training data set is significantly higher than random classification. 
This indicates that a cultivar based assignment of a tap root to the corresponding cultivar is possible by only using geometrical shape parameters. 
%Thus, the 3D shape can be seen as a robust tool for plant description. 
The introduced approach is not limited to the 3D shape description by laser scanning 
and can be generalized for any other 3D measuring device with high precision as structure-from-motion approaches \cite{Rose2015} or volume-carving methods \cite{Jahnke2016}.


%% LEFTOVERS
% The description of the arithmetic mean as the position in the vector space which minimizes the sum of squared distances to the input positions 
% inspires the definition of a shape mean as the shape which minimizes the sum of integral deformation energies. 
% Furthermore, the covariance operator on a vector space can be represented by a matrix whose entries are scalar products of displacement vectors from the input positions to the average position.
% In our setup these are integral scalar products of displacement fields associated with the deformations of the input shapes into the average shape.
% The PCA is then based on a spectral decomposition of this representing matrix.
% The resulting eigenvectors depict modes of main variations within the data set and are used to define an extended version of the classical Mahalanobis distance.
% This distance measure allows us to define a classification tool that is able to distinguish between different cultivars.
% Conceptually, our  approach is motivated by the theory of minimizers of physical energy associated with the deformation of one shape into another shape \cite{Ci97} where its mathematical concept is given in detail in \cite{RuWi08,RuWi09b}. 
% Here, we develop this approach embedded in a processing pipeline for the shape analysis of sugar beet tap roots.\\[1ex]

%Here we demonstrate the extraction and mathematically description of shape features on a large set of sugar beet tap roots. 
%Thereby, the huge amount of measurement data consisting of detailed 3D descriptions of all individual root samples has been reduced to a small number of important parameters without losing relevant biological information. 
%To this end, we revisit a principal component analysis (PCA) in (linear) vector spaces \cite{HaBlBaAn} and generalize this concept to nonlinear shape analysis.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results}
We introduce a method for the statistical analysis of large ensembles of 3D tap roots as well as a classification tool based on information gained from the statistical analysis.
To prove the benefits of our approach, four cultivars with 35 sugar beet tap roots each were measured with a laser scanner. 
After several pre-processing steps each beet is finally represented as a characteristic function\footnote{A characteristic function assigns $1\!=\!true$ to every data point inside the sugar beet volume and $0\!=\!false$ to all other points.}. 
In order to validate our classification tool, we randomly removed 5 beets of each cultivar from the training data set.
Afterwards, we computed a \emph{shape average} as the minimizer of an elastic matching functional for each of the four training data sets separately.
In the same physical setup, a principal component analysis has been applied to the displacements between average and input shapes. 
This way, we obtained principal modes of main variation within all four training data sets.
Again, each of these modes is represented by a displacement of the shape average. 
%and it is shown in the subsequent figures as a variation of the shape average. 
%\oldArticle{ approach was tested using sugar beets within different phases of growth.
%Furthermore we calculated average and main variations for different cultivars of sugar beets within the same phase of growth.}

The results are shown in Fig.~\ref{fig:Berenika} to Fig.~\ref{fig:Pauletta} for the four different cultivars. 
In each figure we show the input shapes on the left and the five main shape variations on the right, 
where the middle shape (dark gray) represents the average beet shape. 
To better illustrate the directions of main variations, the displacements are displayed with different positive and negative magnitudes. 
Note that all input beets were scaled to have uniform volume, hence the previously most dominant mode of uniform scaling is no longer relevant. 
However, the relative order within the variances has not been effected substantially by the rescaling. 
In detail, if $(\tilde\lambda_k)_{k\geq 0}$ describe the variances (in decreasing order) before the uniform rescaling, 
where $\tilde\lambda_0$ is the dominant eigenvalue that represents the impact of uniform scaling,
and $(\lambda_k)_{k\geq 1}$  the variances after rescaling, we observe $\lambda_{k-1}/\lambda_{k} \approx \tilde \lambda_{k-1} / \tilde\lambda_{k}$ for $k > 1$.

%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
\centering
\vspace*{-2mm}
 % trim=left bottom right top, clip
\includegraphics[width=0.99\columnwidth, trim = 32 5 35 5, clip ]{figure1.pdf}
\vspace*{-7mm}
\myCaption{Input beets of sort Berenika (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).}
\label{fig:Berenika}
\end{figure}
%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
 \setlength{\unitlength}{\linewidth}
 \centering
\vspace*{-2mm}
 % trim=left bottom right top, clip
\includegraphics[width=0.99\columnwidth, trim = 32 5 35 5, clip ]{figure2.pdf}
\vspace*{-7mm}
\myCaption{Input beets of sort Cesira (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).}
\label{fig:Cesira}
\end{figure}
%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
 \setlength{\unitlength}{\linewidth}
 \centering
\vspace*{-2mm}
 % trim=left bottom right top, clip
\includegraphics[width=0.99\columnwidth, trim = 32 5 35 5, clip ]{figure3.pdf}
\vspace*{-7mm}
\myCaption{Input beets of sort Mauricia (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).}
\label{fig:Mauricia}
\end{figure}
%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
 \setlength{\unitlength}{\linewidth}
 \centering
\vspace*{-2mm}
 % trim=left bottom right top, clip
\includegraphics[width=0.99\columnwidth, trim = 32 5 35 5, clip ]{figure4.pdf}
\vspace*{-7mm}
\myCaption{Input beets of sort Pauletta (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).}
\label{fig:Pauletta}
\end{figure}
%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff


%Result describing
%Describing the shown modes of variation is rather complicated as these modes are not intuitive. 
For the variation in  cultivar \emph{Berenika} (see Fig.~\ref{fig:Berenika}) the first mode can be described as a tendency to a multiple or at least double apex (2.01), the second variation as a tendency to a long apex or a more dull one (1.00). 
For the cultivar \emph{Cesira} (see Fig.~\ref{fig:Cesira})  the first variation (1.99) is according to the second mode of variation of \emph{Berenika}, the tendency from a long apex to a dull one. The second variation is similar to the first one of \emph{Berenika} (1.63). 
The main variation for the cultivar \emph{Mauricia} (see Fig.~\ref{fig:Mauricia}) is similar to the long apex or a dull one (0.94), the second variation is the affinity to a multiple or clear and pointy apex (0.79). 
The first two main variations of \emph{Pauletta} (see Fig.~\ref{fig:Pauletta}) can be described as the affinity to a pointy or dull apex (0.93) 
and the tendency from a clear and pointy apex to a second apex at the side (0.40).

%All in all, the expression of the variation within the single cultivars differs enormous. 
%As the cultivar Berenika has a variance measurement from 2.0 to 0.51, 
%the same reaches in the cultivar Cesira from 1.99 to 0.47, in the cultivar Mauricia from 0.94 to 0.41 and in the cultivar Pauletta from 0.93 to 0.24.

%\missing{unit of variation!}

In a second step, the principal modes were used to derive a classification tool. 
Based on the classical \emph{Mahalanobis distance}, we first define a distance measure for each cultivar.
Then a given beet is classified to the cultivar that induces the lowest distance. 
For example, when computing the distance of that given beet to cultivar \emph{Pauletta},
we effectively compute the distance to the average shape of  the training data set of \emph{Pauletta}. 
In particular, the distance measure penalizes deviations in direction of dominant principal modes of this data set less, 
as these directions are supposed to be characteristic for cultivar \emph{Pauletta}.

%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
  \setlength{\unitlength}{\linewidth}
  \centering
\vspace*{-2mm}
 % trim=left bottom right top, clip
\includegraphics[width=0.99\columnwidth, trim = 32 0 35 0, clip ]{figure5.pdf}
\vspace*{-7mm}
 \myCaption{Five beets of sort Berenika (blue), Cesira (red), Mauricia (black) and Pauletta (green) to be classified. 
Vertical color bars upon each shape are proportional to the distance to the four different cultivars. 
 %Extended Mahalanobis distances (with $M_j=13$ and $\beta = 10^{-5}$) for every beet to the four training sets. 
 Overall $55\%$ beets were classified correctly (whereas $25\%$ is random).}
 \label{fig:classification}
 \end{figure}
%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

Classification results are depicted in Fig.~\ref{fig:classification} and Fig.~\ref{fig:confusion}, respectively.
In detail, we show for every beet to be classified the (extended) Mahalanobis distance %(with $\beta = 10^{-5}$) 
to each of the four cultivars 
as vertical color bars whose height is proportional to the distance (cf.~Fig.~\ref{fig:classification}). 
For example, the height of the blue bar always represents the distance to the \emph{Berenika} training set. 
%The confusion matrix shows a classification accuracy of $20\%$ for the cultivar Berenika, $60\%$ for the cultivar Cesira, $60\%$ for the cultivar Mauricia and $80\%$ for the cultivar Pauletta.
The validation of the classification  reveals a significantly better success (i.e. $55\%$) in comparison to random classification (with a $p$-value of $0.0045$).
Note that we considered only the first $13$ principal modes for each cultivar to design the distance measures (details will be explained in the "Methods" section). 
%since this was the lowest number that realized the best classification result\footnote{For $M_j < 13$ the classification success was monotonically increasing in $M_j$ while it was roughly constant for $M_j > 13$.}.
If we ignore the principal modes and simply compute the distance to the average %(i.e. $M_j = 0$) 
we still obtain a correct overall classification rate of $40\%$. 
However, we observe huge differences in classification success when distinguishing between cultivars (cf.~Fig.~\ref{fig:confusion}). % (cf.~Fig.~\ref{fig:classification}, B). 
In detail, a beet from cultivar \emph{Pauletta} was classified correctly in $80\%$ of the trials, whereas \emph{Berenika} classification success was more or less random (i.e. $20\%$). 


%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
 % trim=left bottom right top, clip
\includegraphics[width=0.9\columnwidth, trim=1.7 1.1 4.0 8.6, clip ]{figure6.png}
\myCaption{Visualization of results shown in Fig.~\ref{fig:classification} by means of a confusion matrix to demonstrate differences in classification success for the four different cultivars.}
\label{fig:confusion}
\end{figure}
%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffff


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Simplified PCA approach.} 
For comparison reasons we also computed a Euclidean principal component analysis on five characteristic parameters extracted directly from the 3D point cloud. 
In detail, these five different parameters are root length,  width, surface and volume (cf.~\cite{paulus2014b}) as well as root complexity  
%Root volume and root surface were described by using the volume and the area of the convex hull in 3D .  
which is defined as the quotient between root surface and volume. 
%Using a principal component approach on the 3D data the extension on the first component was defined to be the length. 
%According to this, width is the extension using the second principal component.
%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
\centering
\includegraphics [width=0.9\linewidth] {figure7.pdf}
\myCaption{A principal component analysis of measured root traits based on five measured traits for comparison. 
By using these parameters a differentiation between the cultivars is hard to achieve.}
\label{cultivar_distribution}
\end{figure}
%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
Fig.~\ref{cultivar_distribution} shows the distribution of the different cultivars using the first two principal components of the PCA. 
Obviously, distinguishing different cultivars or even a reliable classification is hard, if not impossible. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion}

The aim of our study was to perform a statistical shape analysis and to develop a method which allows to distinguish between different cultivars in a non-parametric way. 
%Below, we will discuss the identification of non-predefined tuber shape classifiers as a case study and proof of concept.
%The relative advantage of such a procedure is that it would become possible to dissect phenotypic traits in an unbiased and objective way. 
%It would allow as well to differentiate between environment and genotype effects when applying the method to the same lines or cultivars from different environments.
Fig.~\ref{fig:Berenika} to Fig.~\ref{fig:Pauletta} show how our method identifies the main components of shape variation. 
These are mostly in line with traits for sugar beet shapes like size, maximum length, maximum width, average radius, radial variation, circularity, the ratio width to length, and further shape factors like surface roughness or furrow formation. 
%It was found that image processing can be used to identify parameters which were related to sugar yield \cite{Tsialtas2010},
%where the cultivar effect was only seen on the width although it has been identified e.\,g.\ as having a significant impact on all quantitative and qualitative traits of sugar beet yields. 
Moreover, properties that seem to be characteristic for a particular training data set can actually be extracted in form of dominant variations. 
For example, the training data set of \emph{Pauletta} (cf.~Fig.~\ref{fig:Pauletta}, left) suggests that this cultivar tends to have a rather slim and longish body that varies predominantly in length 
whereas branching of the main tap root is usually not preferred. 
On the other hand, the clearly separated first principal mode ($\lambda_1 > 2 \lambda_2$) captures exactly this variation with respect to length 
and none of the first five modes represents a clear branching (compared to the other cultivars in Fig.~\ref{fig:Berenika} to Fig.~\ref{fig:Mauricia}). 
Besides the qualitative modes of variation, our method allows a robust quantification of the associated variance of the main components of 3D shape variation. 
A next step would be to compare the results also at different environments or sets of stress conditions which enables quantitative statements about the influence of these parameters on the plant's development.
The present outcomes suggest that it will become possible to dissect environmental from genotypic effects on the phenotype and thus allow a true G$\times$E$\times$P interaction analysis \cite{Tardieu2010}. 
Furthermore, it is envisaged that the proposed method is adopted to other plant organs and species and likewise will be useful for approaches in precision farming.


%% classification
The results of Fig.~\ref{fig:classification} have shown a classification accuracy of about $55\%$ which is significantly better than a random assignment. 
However, the robustness of the classification can still be improved. 
In particular, the dependence on the number $M_j$ of considered principal modes in the definition of the distance measure has to be further investigated. 
Furthermore, the striking differences of classification success between different cultivars (cf.~Fig.~\ref{fig:confusion}) has to be explored and explained in detail. 
We believe that increasing the number of samples in the training data sets might improve the robustness as well as the classification success of individual cultivars. 
%Furthermore, it might be interesting to determine shape variances for different environments and climate conditions, 
%to enable quantitative statements about the influence of these parameters on the plant's development.

%\missing{Why is Pauletta working so well and Berenika so bad?}

%% future improvements
Currently, the resolution of the shapes is limited by the spatial resolution of the characteristic functions which are defined on a $129^3$ voxel grid. 
The choice of spatial resolution is a trade-off between a desired level of details and limitation of computation time. 
However, the sensor provides a resolution of about $50$ microns, hence it is desirable to be able to use the full resolution in the numerical simulations as well. 
To this end, future work will focus on the reduction of the algorithm's complexity on the one hand, 
e.g. by using linearized elasticity and to use more sophisticated algorithmic features on the other hand, such as adaptive grids and advanced parallel computing techniques. 
An altenative approach is applying the method proposed in \cite{ZhHeRuSm15} to the original triangle meshes constructed directly from the laser scans. 
This \emph{shell PCA} performs the same statistical shape analysis but with a different physical model (in detail, shapes are treated as hollow objects and one studies elastic deformations of the surface only). 
Moreover, it might be interesting to investigate the potential application of machine learning approaches to tackle the classification problem.  
Finally, we aim at the replacement of a manually moved scanner by a device that is able to perform a high number of scans automatically. 
This replacement is in particular necessary since we aim at increasing the number of samples in the training data sets.


%% Übertragbarkeit auf andere Pflanzen? Kulturen? Anwendungen?
%Potatoes, Corn, Cereals,


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
We applied an established computer vision algorithm to perform a statistical shape analysis on 3D laser scanned tap roots. 
The method is applied to sugar beets, where the mean root shape and the main variations within a group of sugar beets of the same growth period have been computed. 
Our investigations can be considered as a case study for the statistical analysis of storage roots without predefined classification criteria. 
The method indeed demonstrates that it is possible to identify principal modes of shape variation automatically and to quantify the associated variances. 
%Based on our findings unknown samples could be linked to their variety using four different varieties for training.
In particular, the resulting dominant modes of variations can be used to cluster scanned tap roots into categories which forms the basis for a linking to growth conditions and environmental conditions.
Furthermore, this statistical shape analysis can be used in combination with other invasive or non-invasive sensors that access the 3D shape and can be applied to other plant organs and species as well.
Since non-invasive sensors such as MRI or radar imaging are usually affected by significant noise
the statistical analysis of large ensembles  of laser scanned root shapes %of various cultivars and grown under different environmental conditions 
could help to increase the accuracy of the classification of noisy input data. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Methods %%
\section{Methods}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Plant material.}
Sugar beets were grown during summer 2012 in a central experiment of CROP.SENSe.net. The objectives of the central experiment were the provision of ground truth data for the interpretation of sensor measurements. 
The soil was a silty loam soil (WRB: Haplic Luvisol) at the Klein Altendorf experimental research station (6$^\circ$59'N, 50$^\circ$37'E) near Bonn, Germany. 
The mean annual temperature is 9.6\,$^\circ$C with an average annual precipitation of 625\,mm.
The cultivars \textit{Pauletta}, \textit{Berenika}, \textit{Mauricia} and \textit{Cesira}---obtained commercially from KWS Saat\footnote{www.kws.de}---were chosen because of their differing phenology. % canopy.
At the end of the growing season 35 randomly chosen beets of each variety were sampled, i.e. 140 beets in total.
Rows for sampling were $2$\,m $\times$ $3$\,m within a plot of $50$\,m$^2$ each. 
Seed density was $90,000$ per hectare at a row distance of $0.2$\,m (plant to plant distance).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Data acquisition.}
To acquire the 3D shape of the sugar beet tap roots a close-up laser scanner coupled to an articulated measuring arm device \cite{Paulus2013, paulus2014b} was used (see Fig.~\ref{fig:realitaet2digital}). 
This is a well evaluated combination for 3D plant imaging \cite{paulus2014c, paulus2014a}.
Hardware details such as resolution and accuracy are given in Table \ref{hardware_details}. 
With its seven degrees of freedom it is possible to image the plant from various viewpoints to get an occlusion-free 3D model. 
The output of the laser scanner is a point cloud with XYZ coordinates with more than 3 million points that is later parsed to an automatic triangulation algorithm. 
%This was reduced in point number, smoothed and a meshing algorithm calculates a triangulation (see Fig.~\ref{fig:realitaet2digital}).

%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{table}[h]
	\myCaption{Hardware details for the measuring setup.}
	\begin{tabular}{ll}
		\hline
		\multicolumn{2}{c}{3D scanning hardware} \\ \hline
		Manufacturer supporter & Hexagon Metrology Ltd. UK \\
		Model supporter & Romer Infinite 2.0 \\
		Manufacturer scanner & Perceptron Inc. Plymouth, MI, USA \\
		Model scanner & Perceptron Scan Works V5 \\
		Resolution & $17 \mu m$ \\
		Accuracy & $45 \mu m$ \\
		Measureable Volume & spherical (1.4m radius) \\
		used Wavelength & 660nm \\
		main application field & quality management \\ \hline
	\end{tabular}
\label{hardware_details}
\end{table}
%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}%[ht]
\centering
%figure1
\includegraphics[width=0.8\linewidth]{figure8.png}
\myCaption{The laserscanner-measuring arm combination with its seven degrees of freedom enables the 3D imaging of the complete sugar beet tap root. 
It provides a point accuracy of $45\mu m$ and a point resolution of $17 \mu m$ (A). 
An RGB image of a sugar beet (B) is shown together with its laserscanned point cloud with a shaded visualization with 3 million points (C) and the down-sampled and smoothed wireframe representation (D)}
\label{fig:realitaet2digital}
\end{figure}
%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Preprocessing (cf.~Fig.~\ref{fig:dataflow}, left column).}
Before being fed to our statistical method,
outliers and background measurements were removed from the recorded point cloud data sets, and a subsequent automatic point cloud triangulation was performed (see Fig.~\ref{fig:realitaet2digital}). 
We then applied a uniform rescaling such that all shapes have the same (inscribed) volume\footnote{The rescaling proved to improve the later classification substantially while it does not qualitatively change the shape analysis. See discussion above.}. 
 
%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure*}%[ht]
 \setlength{\unitlength}{\linewidth}
 \centering
 \vspace*{-3mm}
\includegraphics[width=0.98\textwidth]{figure9.pdf}
\vspace*{-8mm}
\myCaption{Workflow of our method for data pre-processing (left), statistical analysis (middle) and classification (right). }
\label{fig:dataflow}
\end{figure*}
%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff 
 
The proposed statistical method requires the representation of each root shape $\shape$ %after a suitable and uniform rescaling
via a characteristic function $\chi_{\shape}$ on the computational domain $\Omega=[0,1]^3$. %This rescaling transforms all shapes in a way that they have a constant volume.
%\missing{Rescaling of whole data set to have constant volume!}
%This characteristic function assigns $1=true$ to every point inside the sugar beet volume and $0=false$ to all other points.
Then, the sugar beet surface $\shape$ is implicitly represented as the interface between these two regions.
The function $\chi_{\shape}$ is discretized on a regular grid with $129^3$ nodes (129 in each spatial direction).
It is obtained from the triangulated sugar beet surfaces by first computing the signed distance function $d:\Omega\to\R$ \cite{RuSm00a}, a function whose value at each point is the shortest distance to the surface with positive sign inside and negative sign outside the sugar beet. 
To this end, we first compute the signed distance on the grid nodes closest to the triangles and then we compute its values far from the triangular surface via information propagation by a fast marching algorithm \cite{Se96a}.
Since the triangulated surface may be non-closed due to local scan deficiencies or self-occlusion, the signed distance function $d$ may have an incorrect sign at a number of nodes.
As a remedy, we regularize the resulting signed distance function replacing $d$ by the minimizer $d_{\mathrm{new}}$ of the following
energy
\begin{equation*}
\int_\Omega(|d|_\epsilon-|d_{\mathrm{new}}|_\epsilon)^2\,\d x+\sigma\int_\Omega|\nabla d_{\mathrm{new}}|^2\,\d x\,
\end{equation*}
where $|d|_\epsilon=\sqrt{|d|^2+\epsilon^2}$ with $\epsilon=10^{-4}$, $\nabla$ denotes the gradient mapping, and the filter width parameter is $\sigma=0.01$ in our application.
Computations are performed using a standard finite element approach with multilinear basis function on the regular grid.
Once $d_{\mathrm{new}}$ is computed, we set $\chi_\shape(x)=1$ if $d_{\mathrm{new}}(x)<0$ and $\chi_\shape(x)=0$ else,
which in the algorithm is done for all grid nodes.
After data compression, the required memory to store the function $\chi_\shape$ is comparable to the one needed to store an explicit triangular representation of the shape $\shape$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Statistical analysis model (cf.~Fig.~\ref{fig:dataflow}, middle column).}
In the following we will generalize the standard statistical analysis of point sets in a (linear) vector space, e.\,g.\ $\R^n$,
to the space of shapes $\shapespace$, where $\shape \in \shapespace$ represents the surface of a volumetric object $ \object = \object(\shape) \subset \R^n $.
To underpin this generalization let us briefly recall the basic concepts of the mean and the principal component analysis (PCA)
on vector spaces \cite{HaBlBaAn}.  

The arithmetic mean $\bar x$ of $m$ points $x_1,\,x_2,\, \ldots, \, x_m$ in $\R^n$ is given as
$\bar x = \frac1m (x_1 + x_2 +  \ldots + x_m)$. It can also be characterized as the minimizer of the
energy $E[x] = \sum_{i=1}^m \W[x_i-x]$, where $\W[y]=\mu |y|^2$ is the elastic energy stored in a spring stretched along a vector $y$
and $\mu$ is the associated stiffness coefficient of the spring. Indeed, the arithmetic mean is the equilibrium position of the central hub of a network with $m$ springs, where one end of each spring is attached to one of the input points and the other end is connected to all other springs at a hub so that all springs pull at this hub.
If $X = [x_1 - \bar x|\ldots|x_m - \bar x]$ is the $n\times m$ (centered) \textit{data matrix}, the covariance matrix is given by $\tfrac1mXX^T$.
A PCA now involves a spectral decomposition of this $n\times n$ matrix, where $n$ might be very large.
Hence we use the fact that the $m\times m$ matrix $C= \tfrac1m X^T X$ has the same (non-trivial) eigenvalues as the covariance matrix, where $m \ll n$.
Note that $C_{ij} = \frac1m g(x_i-\bar x,x_j-\bar x)$ where $g$ denotes a suitable scalar product on $\R^n$ 
(typically the standard Euclidean inner product $(x_i-\bar x)\cdot(x_j-\bar x)$).
If $C=Q\Lambda Q^T$ is the singular value decomposition of $C$ 
with an orthogonal matrix $Q$ and a diagonal matrix $\Lambda$ with diagonal entries $\lambda_1 \geq \lambda_2\geq \ldots \geq \lambda_m$,
then the eigenvectors of the covariance matrix are obtained via $w_j = \sum_{i=1}^m \sqrt{\lambda_j^{-1}} q_j^i(x_i-\bar x)$.
These are exactly the principal directions of variations, and the eigenvalues $\lambda_j$ describe the variance in that direction.
Here, $q_j^i$ denotes the $i$\textsuperscript{th} entry in the $j$\textsuperscript{th} column $Q$.
As $g$ is supposed to be positive (semi-)definite, we have $\lambda_i \geq 0$ for $i = 1, \ldots, m$. 
Furthermore, if we assume $x_1, \ldots, x_m$ to be linearly independent we have $\lambda_1 \geq \ldots \geq \lambda_{m-1} > \lambda_m = 0$ as we consider centered data, 
i.e. $x_1-\bar x, \ldots , x_m - \bar x$ spans a $(m-1)$-dimensional subspace.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now  this general concept is transferred to shape statistics.
In the context of tap root shapes $\shape_i$ for $i=1,\ldots, m$, 
the corresponding volumetric objects $\object_i = \object(\shape_i)$ play the role of the points $x_i$, 
and the vector $x_i-x$ is replaced by the optimal deformation $\phi_i:\object_i\to\R^n$ of the object
$\object_i$ into an object $\object$. Here, optimal means that the deformation $\phi_i$ costs the least elastic energy
$\W[\psi_i, \object_i, \object] = \int_{\object_i} W(D \psi_i) \d x$ of all deformations with $\psi_i(\object_i) = \object$, where
$W(\cdot)$ denotes a hyperelastic energy density acting on the local deformation gradient $D \psi_i\in\R^{n\times n}$.   
This energy and thus the thereby defined shape mean is invariant under rigid body motions,
i.\,e.\ it does not change if the position or the orientation of the root shapes is varied.
The arithmetic mean of $m$ input objects $\object_1,\,\object_2,\, \ldots, \, \object_m$ is defined as the object 
$\bar \object$ which minimizes the energy
\begin{align*}
 E[\object] = \sum_{i=1}^m \W[\phi_i,\object_i,\object]\,.
\end{align*}
For details we refer to \cite{RuWi08}.
Different from the vectors $x_i-\bar x$ above, the deformations $\phi_i$ are strongly nonlinear and cannot be used in a (linear) PCA.
Instead, one can replace the deformation $\phi_i$ by its linear representative, the associated elastic stress (cf.~Fig.~\ref{fig:stressBalance})
$$
\sigma_i = W'((D(\phi_i^{-1}))^{-1}) (\det D(\phi_i^{-1})) (D(\phi_i^{-1}))^{-T} \nu\,,
$$
evaluated on the average shape $\bar \shape$ (the surface of the average object $\bar\object$)
with surface normal vector $\nu$ \cite{Ci97}.
Even more intuitive, one can also replace the nonlinear deformation $\phi_i$
by the displacement $u_i:\bar\shape\to\R^n$ of each point on the surface $\bar\shape$ of $\bar\object$
which is observed if the elastic stress $\sigma_i$ is applied at the surface of the average object.
The Hessian of the corresponding elastic energy
implies a natural scalar product $g(\cdot, \cdot)$ on these displacements $u_i$ \cite{Ci97,RuWi09b}.
Thus, we define the matrix 
$$	C = \frac1m  \left( g(u_i,u_j)\right)_{i,j=1,\ldots, m}	$$ as a representation of the covariance operator, % $\cov :\shapespace \rightarrow \shapespace$, 
which is defined by 
$$ \cov \, u = \frac1m \sum_{i=1}^m g(u_i, u) u_i $$
and can be regarded as the analogon of the covariance matrix in the vector space case.
Again, the whole construction is rigid body motion invariant.
Via the same spectral decomposition  $C=Q \Lambda Q^T$ of the $m\times m$ matrix $C$ as above, 
we finally obtain the principal modes of shape variation,
$$
w_j(x)= \sum_{i=1}^m \sqrt{\lambda_j^{-1}} q_j^i u_i(x)\,,
$$
and associated variances $\lambda_j$.
A comprehensive introduction to this concept can be found in \cite{RuWi09, RuWi10}.
Let us emphasize that this statistical approach does not make any assumptions on how the different shapes and their variations are configured. 
It is also not assumed that the original shapes are actually elastic---this is just a mathematical tool to define the dissimilarity between shapes.
%

%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
\begin{figure}[h]
\centering
\vspace*{-2mm}
 % trim=left bottom right top, clip
\includegraphics[width=0.99\columnwidth, trim = 65 0 50 0, clip ]{figure10.pdf}
\vspace*{-5mm}
\myCaption{Sketch of the stresses $\sigma_i$ induced by the deformations $\phi_i$ on the averaged shape $\bar \shape = \partial\bar\object$, here for $i=1,2,3$.}
\label{fig:stressBalance}
\end{figure}
%fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Classification (cf.~Fig.~\ref{fig:dataflow}, right column).}
Sample beets were available for $N = 4$ different cultivars.
The training data set of cultivar $j$ is represented as a set of displacements $(u^j_k)_{k\leq K_j}$ from the group average $\bar\shape_j$, where $K_j$ denotes the sample size.
The aim of classification is to assign an arbitrary shape $\shape$ or a corresponding shape displacement $u$, respectively, to one of these $N$ cultivars.
In detail, we need a distance measure $d_j$ for each cultivar that quantifies the distance of $u$ to the set $(u^j_k)_{k\leq K_j}$.
We say that $u$ is classified to cultivar $l$ if $d_l(u) \leq d_j(u)$ for $j = 1, \ldots, N$.
	
Our distance measure $d_j$ is an extension of the classical Mahalanobis distance. 
First, we assume that we have already performed a PCA (as described above) for each cultivar independently. 
Let  $(\lambda^j_k)_{k\leq K_j}$ and $(w^j_k)_{k\leq K_j}$ denote the eigenvalues and eigenvectors computed in the PCA of cultivar $j$, respectively.
We  define $U_j$ to be the linear span of $(u^j_k)_{k\leq K_j}$ and $U_j^\perp$ the orthogonal complement with respect to the metric $g_j(.,.)$ associated by the $j$th cultivar.  
%For $K_j = n$, 
For some truncation value $M_j \leq K_j$ the classical (squared) Mahalanobis distance is given by
\begin{equation*}
 m_j^2(u) := \sum_{k=1}^{M_j} \frac{g_j(u, w_k^j)^2}{\lambda_k^j }\, .
\end{equation*}
Note that $m_j^2(w^j_k) = (\lambda^j_k)^{-1}$, that means the Mahalanobis distance penalizes deviations in direction of dominant eigenvectors less.
On the other hand, we have $m_j(u^\perp) = 0$ for each $u^\perp \in U_j^\perp$, which requires some kind of regularization.
If $P_j u$ is a projection of $u$ onto $U_j$ we decompose $u = P_j u + (u\!-\!P_j u) \in U_j \oplus U_j^\perp$ and finally define an \emph{extended Mahalanobis distance} by
\begin{equation*}
\maha_j^2(u) = m_j^2(P_j u) + \frac1\beta g_j( u\!-\!P_j u, u\!-\!P_j u ) \, .   
\end{equation*}
Here $\beta>0$ is a regularization parameter. 
In the experiments shown in Fig.~\ref{fig:classification} we have chosen $\beta = 10^{-5}$ and $M_j = 13$ for all four cultivars.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Numerical implementation.}
The matching deformations $\phi_i$ between an input object $\object_i$ and the object $\object$, both represented by characteristic functions
$\chi_i$ and $\chi$, respectively,  are computed using a penalty approach, where one minimizes
\begin{equation*}
\int_{\object_i} W(D \phi) \d x+ \gamma \int_{[0,1]^3}(\chi_i-\chi \circ\phi)^2\,\d x
\end{equation*}
(for some penalty parameter $\gamma=50$) with respect to the deformation $\phi$.
The first term ensures that the deformation with minimum energy is found,
and the second term ensures that $\phi$ indeed deforms $\object_i$ into $\object$.
%For the sake of completeness, let us give here the concrete energy density $W(\cdot)$ used in the algorithm:
In detail, we made use of the polyconvex energy density $W(D\phi) = \bar W( \|D\phi\|, \det D\phi)$ with
\begin{equation*}
\bar W(a,d) =  \frac\mu2 a^2 +\frac\lambda4 d^2-(\mu+\frac\lambda2)\log d-\frac{3\mu}2-\frac\lambda4\, ,
\end{equation*}
where the elastic parameters are $\mu=\lambda=1$.
Consequently, the average shape $\bar\shape$, represented by $\chi$, and the deformations $\phi_1,\ldots,\phi_m$
are obtained by minimizing
\begin{equation*}
\sum_{i=1}^m \int_{\object_i} W(D \phi_i) \d x+\gamma \int_{[0,1]^3}(\chi_i-\chi \circ\phi_i)^2\,\d x
\end{equation*}
over all $\phi_i$ and $\chi$. For details we again refer to \cite{RuWi08}.
The functions $\chi_i$ and deformations $\phi_i$ are discretized by $129^3$ nodes
(i.\,e.\ a function value is assigned to each node of a grid with $129$ vertices in each space direction)
so that the number of degrees of freedom scales with $129^3\cdot3\cdot m$,
yielding a high-dimensional problem.
The principal component analysis basically involves a singular value decomposition of a symmetric $m\times m$ correlation matrix
and the solution of $m$ linear systems of equations of size $129^3$ 
in order to translate the boundary stresses into displacements or shape variations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% MANDATORY SECTION REQIRED BY THE JOURNAL (INCLUDING ALL SUBSECTIONS!)
\subsection*{Abbreviations} 
%PGA: principal geodesic analysis; 
G$\times$E$\times$P: Genotype $\times$ Environment $\times$ Phenotype;
PCA: principle component analysis

%\section*{Declarations}

\subsection*{Ethics approval and consent to participate} Not applicable.

\subsection*{Consent for publication} Not applicable.

\subsection*{Availability of data and materials} The datasets used and analysed during the current study are available from the corresponding author on reasonable request.

\subsection*{Competing interests} Not applicable.
% The authors declare that they have no competing interests.

\subsection* {Funding}
This study was funded by the German Federal Ministry of Education and Research (BMBF) within the scope of the interdisciplinary network program \emph{CROP.SENSe.net} (Funding code: 0315529). 
The BMBF played no role in the design of the study, the collection, analysis, and interpretation of data and in writing the manuscript.

%This work could be carried out due to the financial support of the German Federal Ministry of Education and Research (BMBF) within the scope of the competitive grants program "Networks of excellence in agricultural and nutrition research - CROP.SENSe.net (Funding code: 0315529).

\subsection*{Author's contributions}
SP, HK, HG, MR and BW designed the study. 
HG maintained the plant material, 
SP carried out the measuring, pre-processing was done by SP and BH. 
BH and BW implemented the main algorithm and performed all numerical simulations as well as the post-processing. 
SP, BH, AKM, BW and MR interpreted the data and drafted the manuscript.  
%SP, HG, MR, HK and BW directed the research and gave initial input. 
All authors read and approved the final manuscript.

\subsection*{Acknowledgements} Not applicable.
%KWS for the plant material.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Backmatter begins here                   %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{backmatter}

  %This work could be carried out due to the financial support of the German Federal Ministry of Education and Research (BMBF) within the scope of the competitive grants program "Networks of excellence in agricultural and nutrition research - CROP.SENSe.net (Funding code: 0315529), subprojects D 2 and GS 4. The authors also want to express their gratitude to \ldots and to \ldots for proofreading the manuscript. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if your bibliography is in bibtex format, use those commands:
\bibliographystyle{bmc-mathphys} % Style BST file (bmc-mathphys, vancouver, spbasic).
\bibliography{bmc_article}      % Bibliography file (usually '*.bib' )
% for author-year bibliography (bmc-mathphys or spbasic)
% a) write to bib file (bmc-mathphys only)
% @settings{label, options="nameyear"}
% b) uncomment next line
%\nocite{label}

% or include bibliography directly:
% \begin{thebibliography}
% \bibitem{b1}
% \end{thebibliography}
\end{backmatter}

\subsection*{Figure captions} 

\begin{enumerate}
\item Input beets of sort Berenika (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).
\item Input beets of sort Cesira (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).
\item Input beets of sort Mauricia (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).
\item Input beets of sort Pauletta (left) and first five modes of main variation (right) with average beet in dark grey and corresponding variances $\lambda_1, \ldots, \lambda_5$ (multiplied by $100$).
\item Five beets of sort Berenika (blue), Cesira (red), Mauricia (black) and Pauletta (green) to be classified. 
Vertical color bars upon each shape are proportional to the distance to the four different cultivars. 
 Overall $55\%$ beets were classified correctly (whereas $25\%$ is random).
 \item Visualization of results shown in Fig.~\ref{fig:classification} by means of a confusion matrix to demonstrate differences in classification success for the four different cultivars. 
  \item A principal component analysis of measured root traits based on five measured traits for comparison. 
By using these parameters a differentiation between the cultivars is hard to achieve.
  \item The laserscanner-measuring arm combination with its seven degrees of freedom enables the 3D imaging of the complete sugar beet tap root. 
It provides a point accuracy of $45\mu m$ and a point resolution of $17 \mu m$ (A). 
An RGB image of a sugar beet (B) is shown together with its laserscanned point cloud with a shaded visualization with 3 million points (C) and the down-sampled and smoothed wireframe representation (D).
\item Workflow of our method for data pre-processing (left), statistical analysis (middle) and classification (right).
\item Sketch of the stresses $\sigma_i$ induced by the deformations $\phi_i$ on the averaged shape $\bar \shape = \partial\bar\object$, here for $i=1,2,3$.
\end{enumerate}


\end{document}
