In our approach, the estimation of core positions on the input TMA image is subsequently refined in successive tasks by considering different image domains (i.e. patches or regions) in the input original image. Such a strategy allows not only to avoid unnecessary processing on non-content regions but also to reduce the acquisition time, storage and processing time of high resolution data. To distinguish the inputs and outputs of each task and facilitate the comprehension of the technical details, we present in Fig. 3 a diagram which illustrates a few notations which will be used throughout the paper.
TMA core detection
The detection of approximately circular TMA cores can be performed by spot detection algorithms. Spot detection is a well-known topic in image processing (see [26] for a recent review). Over past decades, number of spot detection methods have been proposed [27–31]. To produce satisfactory results, most of these methods require fine adjustment of a critical parameter: the detection scale corresponding to the size of the objects of interest. Automatic selection of the detection scale is a challenging problem since the objects of interest may have different sizes or they may have the same size as the irrelevant objects in the background. Few methods of automatic scale selection [32–34] have been proposed recently. However, in the context of tissue microarrays, the diameter of assembled TMA cores is defined by the size of the needle used for extracting cores from paraffin tissue blocks. The determination of the scale parameter used for spot detection is straightforward from this measure which is often given by the manufacturer. We propose here a fast algorithm for tissue core detection by performing directly the wavelet decomposition at the appropriate aforementioned scale and computing a locally-adaptive threshold of the wavelet coefficients.
Pre-processing
Wavelet-based detection techniques are known to be robust to non-stationary noises like Poisson noise or mixed-Poisson-Gaussian noise as in TMA images, acquired by brightfield or fluorescence. Pre-processing operations such as image denoising or variance stabilization transform are thus not mandatory. However, our algorithm is primarily designed for detecting bright spots over a dark background, and especially adapted for fluorescence images. brightfield TMA images, in which the tissue cores are darker than the background, are first inverted before further processing.
Scale selection
Nowadays, standard TMAs are manufactured with a diameter of tissue core typically from 0.6 to 1.5 mm. Given an imaging resolution (pixel size), the optimal scale of wavelet decomposition can be determined according to the core radius. If we denote r
core
the expected average radius (in pixels) of TMA cores, the optimal scale index \({\hat {\jmath }}\) of the wavelet decomposition that best fits the size of TMA cores is defined as:
$$ \hat{\jmath} = \underset{j \in \mathbb{N}^{*}}{\text{argmin}} \left| {r}_{core} - 2^{j-1}\sigma_{1} \right| \, $$
(1)
where σ1 is selected according to the pixel size.
Fast isotropic wavelet decomposition
In contrast to multiresolution approaches, our detection method requires only the wavelet decomposition at the appropriate selected scale. To compute the decomposition at a desired scale, usual wavelet transform techniques perform a sequence of successive convolutions which are used for computing iteratively the decomposition from the smallest scale to the coarsest scale. These techniques are time consuming when dealing with large images and high number of scales. Instead, to address to this computational issue, we build a dyadic isotropic wavelet frame \(\left \lbrace \psi _{j} \right \rbrace _{j \geq 1}\) by choosing the scaling function ϕ
j
as a Gaussian function whose variance \(v_{j}^{2}\) is a function of scale \(j \in \{1, \ldots, j_{\max }\}\) and \(j_{\max }\) is the maximum index of the highest scale:
$$ \phi_{j}(\mathbf{x}) = G_{v_{j}}(\mathbf{x}) = \frac{1}{{2\pi v_{j}}}\exp \left(-\frac{\left\| \mathbf{x} \right\|_{2}^{2}}{2v_{j}^{2}} \right), $$
(2)
where ∥·∥2 denotes the Euclidean norm, \(\mathbf {x} \in \Omega \subset \mathbb {R}^{2} \) is the pixel location in the rectangular domain Ω and
$$ v_{j}^{2} = \sum_{k = 1}^{j} \sigma_{k}^{2} = \sum_{k = 1}^{j} 4^{k-1}\sigma_{1}^{2} = \sigma_{j}^{2} + v_{j-1}^{2} $$
(3)
with σ
k
=2k−1 σ1. Thanks to the semi-group property of Gaussian functions, the relationship between the scaling functions at subsequent scales can be expressed as:
$$\begin{array}{@{}rcl@{}} \phi_{j}(\mathbf{x}) &=& G_{\sqrt{\sigma_{j}^{2} + v^{2}_{j-1}}} (\mathbf{x})\\ &=& G_{\sigma_{j}} \star G_{v_{j-1}}(\mathbf{x}) = G_{\sigma_{j}} \star \phi_{j-1}(\mathbf{x}), \end{array} $$
(4)
where ⋆ denotes the convolution operator. Therefore, the wavelet decomposition Ψ
j
u of \(u: \Omega \subset \mathbb {R}^{2} \rightarrow \mathbb {R} \) at the scale \(j \in \{1, \ldots, j_{\max }\}\) is obtained by convolution of u with the wavelet atom ψ
j
as:
$$\begin{array}{@{}rcl@{}} \Psi_{j} {u}(\mathbf{x}) = \psi_{j} \star {u}(\mathbf{x}) & = &\left(\phi_{j - 1} - \phi_{j} \right) \star {u}(\mathbf{x}) \\ & = & (G_{v_{j-1}} - G_{v_{j}}) \star {u}(\mathbf{x}), \end{array} $$
with the conventions \(v_{0}^{2} = 0\) and G0(·)=δ(·) (Dirac delta function). For more technical details on the proposed wavelet frame and the wavelet decomposition and reconstruction algorithms, please refer to the Additional file 1.
Locally-adaptive thresholding
While the wavelet decomposition plays the role of a filtering which reduces the noise and enhances the objects of interest, a common way to detect objects is to threshold the filtered image – the wavelet decomposition of the input TMA image in our case. As depicted in [32], a global threshold is not appropriate to handle complex situations, especially when dealing with images acquired in fluorescence context because of their inhomogeneous background. To overcome this difficulty, we propose to define an adaptive threshold according to the local distribution of the wavelet decomposition \(\Psi _{\hat {\jmath }}{u}\) previously computed. Accordingly, we consider the following statistical test at each point x of the TMA image u:
$$\begin{array}{@{}rcl@{}} \left\lbrace \!\! \begin{array}{ll} \mathcal{H}_{0}: & \mathbf{x}\, \text{belongs to the background,} \\ \mathcal{H}_{1}: & \mathbf{x}\, \text{corresponds to tissue core (foreground).} \end{array} \right. \end{array} $$
Pixels corresponding to tissue cores have strong positive responses in the wavelet decomposition. Under the null hypothesis \(\mathcal {H}_{0}\), the wavelet coefficient \(\Psi _{\hat {\jmath }}{u}(\mathbf {x})\), which follows the local distribution of the wavelet-decomposed-image background with mean μ(x) and variance ν2(x), is lower than a certain value τ(x). Let \(\mathbb {P} \left (\Psi _{\hat {\jmath }}{u}(\mathbf {x}) < \tau (\mathbf {x}) \right)\) be the probability for a pixel x to be classified as “background” class. The threshold τ(x) is used to control the number of misclassification. Given a probability of false alarm \(p_{_{\text {FA}}} > 0\), the corresponding threshold \(\tau _{_{\text {FA}}}\) is selected such that the misclassification probability \(\mathbb {P} \left (\Psi _{\hat {\jmath }} {u}(\mathbf {x}) \geq \tau _{_{\text {FA}}}(\mathbf {x}) \right)\) is lower than \(p_{_{\text {FA}}}\). By applying the conventional probabilistic Tchebychev’s inequality, we get, ∀κ(x)>0:
$$ \mathbb{P} \left(\left| \Psi_{\hat{\jmath}} {u}(\mathbf{x}) - \mu(\mathbf{x}) \right| \geq \kappa(\mathbf{x}) \right) \leq \frac{\nu^{2}(\mathbf{x})}{\kappa^{2}(\mathbf{x})}. $$
(5)
It follows that
$$\mathbb{P} \left(\Psi_{\hat{\jmath}} {u}(\mathbf{x}) \geq \mu(\mathbf{x}) + \kappa(\mathbf{x}) \right) \leq \mathbb{P} \left(\left|\Psi_{\hat{\jmath}} {u}(\mathbf{x}) - \mu(\mathbf{x}) \right| \geq \kappa(\mathbf{x}) \right). $$
Now, let us define \(\tau _{_{\text {FA}}}(\mathbf {x}) = \mu (\mathbf {x}) + \kappa (\mathbf {x})\) and assume \(\left ({\nu ^{2}(\mathbf {x})}/{\kappa ^{2}(\mathbf {x})}\right) \leq p_{_{\text {FA}}}\) such that \(\mathbb {P} \left (\Psi _{\hat {\jmath }} {u}(\mathbf {x}) \geq \tau _{_{\text {FA}}}(\mathbf {x}) \right) \leq p_{_{\text {FA}}}\). Finally,
$$ \tau_{_{\text{FA}}}(\mathbf{x}) \geq \mu(\mathbf{x}) + \frac{\nu(\mathbf{x})}{\sqrt{p_{_{\text{FA}}}}} $$
(6)
and the adaptive threshold \(\tau _{_{\text {FA}}}(\mathbf {x})\) is controlled by the p-value inferred from the significance level α of the test, and set by the user. If \(p_{_{\text {FA}}} < \alpha \), this suggests that the null hypothesis (i.e. a pixel x is classified as a “background” pixel) may be rejected. In practice, one typically sets \(p_{_{\text {FA}}} = 0.05\) (or 0.01) which corresponds to a significance level α=5% (or 1% respectively).
To determine the threshold \(\tau _{_{\text {FA}}}(\mathbf {x})\), the local mean μ(x) and the local variance ν2(x) of the image background on the wavelet decomposition \(\Psi _{\hat {\jmath }}{u}\) are required. However, prior knowledge about the image background distribution is unfortunately not available in most cases. We consider thus empirical estimations of μ and ν2 at each point x from \(\Psi _{\hat {\jmath }}{u}\):
$$\begin{array}{@{}rcl@{}} \hat{\mu}(\mathbf{x}) &=& g \star \Psi_{\hat{\jmath}}{u}(\mathbf{x}) \, \end{array} $$
(7)
$$\begin{array}{@{}rcl@{}} \hat{\nu}^{2}(\mathbf{x}) &=& g \star \left(\Psi_{\hat{\jmath}}{u}\right)^{2}\!(\mathbf{x}) - \hat{\mu}^{2}(\mathbf{x}) \, \end{array} $$
(8)
where g(·) is a weighting positive function (i.e. ∥g(·)∥1=1, ∥·∥1 is the L1 norm and g(x)≥0,∀x∈Ω) mainly used to avoid the estimation of the background distribution statistics being biased from coefficients corresponding to the foreground. By construction, \(\hat {\mu }(\mathbf {x})\) and \(\hat {\nu }^{2}(\mathbf {x})\) are weighted estimators derived from \(\Psi _{\hat {\jmath }}{u}\) which is a filtered version of u by the band-pass filter \(\psi _{\hat {\jmath }}\) in order to enhance the objects of radius r
core
. It is thus convenient to define the weighting function g according to the wavelet atom \(\psi _{\hat {\jmath }}\). By using an affine transform which implies the positivity and the normalization conditions, we propose a candidate for g(·) as follows :
$$ \hat{g}(\mathbf{x}) = \frac{\ -\psi_{\hat{\jmath}}(\mathbf{x}) + \sup \psi_{\hat{\jmath}} \ }{\ \left\| -\psi_{\hat{\jmath}} + \sup \psi_{\hat{\jmath}} \right\|_{1} \ } \, $$
(9)
where \(\sup \psi _{\hat {\jmath }} = \| \psi _{\hat {\jmath }} \|_{\infty }\) denotes the supremum (\(L_{\infty }\) norm) of \(\psi _{\hat {\jmath }}\) and \(\left \| -\psi _{\hat {\jmath }} + \sup \psi _{\hat {\jmath }} \right \|_{1}\) is the normalization factor to ensure \(\|\hat {g}(\cdot)\|_{1} = 1\). The choice of this candidate is clarified in Fig. 4 showing the wavelet atom and its derived weighting function according to a given circular spot. The proposed weighting which is constructed from the wavelet atom has the same size of the considered spot and has a hollow shape at the center (see right column in Fig. 4). This specific shape allows to reduce the impact of high wavelet coefficients corresponding to foreground pixels on the estimation of the background statistics.
By substituting the empirical estimators to μ(x) and ν2(x), we obtain the estimated detection threshold:
$$ \hat{\tau}_{_{\text{FA}}}(\mathbf{x}) = \hat{\mu}(\mathbf{x}) + \frac{\hat{\nu}(\mathbf{x})}{\sqrt{p_{_{\text{FA}}}}} \ . $$
(10)
Thresholding the wavelet decomposition \(\Psi _{\hat {\jmath }}{u}\) with respect to \(\hat {\tau }_{_{\text {FA}}}\) results in a binary image \(I_{_{\text {FA}}}: \Omega \rightarrow \{0,1\}\):
$$\begin{array}{@{}rcl@{}}I_{_{\text{FA}}} (\mathbf{x}) = \left\{ \begin{array}{ll} 1 & \text{if}~~ \Psi_{\hat{\jmath}}{u} (\mathbf{x}) \geq \tau_{_{\text{FA}}}(\mathbf{x})\\ 0 & \text{otherwise} \end{array} \right. \end{array} $$
(11)
where each connected component in \(I_{_{\text {FA}}}\) represents a region which is potentially a tissue core of the TMA image. The gravity centers of these regions (or detection position) will be used as inputs for estimating the array coordinates of TMA cores. However, the detection reliability has a great impact on the de-arraying outcome: few false detections may lead to severely inaccurate results. Removing false detections (i.e. outliers) is then crucial. To this end, the size of detected regions seems to be a relevant criterion since the core size is given in most cases by the TMA manufacturer. Although, due to the complexity of backgrounds, it may be highly different from the true core size. Instead of exploiting the imprecise information derived from the binary detection map \(I_{_{\text {FA}}}\), we perform an active-contour-based segmentation to delineate the objects at each detected position. Also, we re-use the segmentation results to confirm and improve the preliminary detection results.
Segmentation of TMA cores
As depicted in previous section, the detection binary image \(I_{_{\text {FA}}}\) does not allow us to accurately determine the size of detected objects. Active contours [35] are typically well appropriate in our context since they can evolve to closely delineate the object borders and thus yield an estimation of the TMA core size. The family of parametric active contours presented below will help to refine the detected position and the size of TMA cores and eventually to determine the orientation of the potential core if it was deformed during the manufacturing process.
Since the seminal paper of Kass, Witkin, and Terzopoulos, active contour models (or snakes) [35] have been successfully used to detect discontinuities, detect objects of interest or segment images, especially in bioimaging [36]). General purpose closed contours are generally controlled by elastic forces based on local curvature and image based potentials [35, 37–39]. The curve evolves from its initial starting position towards the target object. The optimization of the underlying energy functional is traditionally performed using variational principles and finite differences techniques, which needs an appropriate initialization to converge to a relevant solution. At the end of the nineties and beginning of the 2000’s, geodesic active contours [40] based on the theory of surfaces evolution and geometric flows have been introduced to segment an arbitrary number of highly complex objects in the image. In our TMA context, the 2D shapes of tissue cores can be actually well estimated by ellipse-shaped active contours which belong to the family of parametric deformable templates.
Application-tailored parametrized templates introduced by Yuille et al. [41] were proposed in cases where strong a priori knowledge about the shape being analyzed is available (e.g. eyes or lips in human faces [41]). The models are hand-built using simple parametrized 2D geometric representations. Another line of research focused on models of random deformations for a given initial shape (deformable template). Grenander et al. [42, 43] obtained the first promising results in image segmentation by considering statistical deformable models which describe the statistics of local deformations applied to an original template. Markov models and Monte-Carlo techniques have been introduced in this context to derive optimal random deformations estimates from image data [42–46]. In the approach initially proposed by Cootes et al. [47] and successfully applied to object tracking [45], the shape structure and the parameters describing its deformations are learned from a training set of representative shapes. Meanwhile, Staib and Duncan [48] proposed to combine parametric snakes (B-splines) to the standard decomposition on a Fourier basis to analyze deformable biomedical structures. All these methods are generally robust to noise but computationally demanding if stochastic iterative procedures are used to conduct the minimization and no initial guess close to the optimal solution is provided. Very recently “snakescules” [49] combined to fast algorithms and Markov point process [50] have been proposed along the same philosophy but dedicated to the detection of cells or nuclei in fluorescence microscopy images.
Finally, the ellipse fitting concept has been furthemore introduced by Thévenaz et al. as an extension of the simple circle-shaped active contour [49] which can be defined just by two points [24]. As a consequence, a triplet of points is necessary to parametrize the ellipse-shaped version. However, this parametrization which has an extra degree of freedom increases the complexity of the model and makes the optimization of ellipse parameters more challenging when compared to the circle-shaped model. To overcome these difficulties, an alternative way was proposed in [51]: the ellipses are configured by their center, their axes and the angle between their major axis and the horizontal. Under this configuration, the cost function introduced in [24], and defined below (see Eq. (12)) as the contrast between the core and the ring defined by the pair of ellipses (see Fig. 5) – and the derivatives with respect to the ellipse parameters could be calculated efficiently by using the Green’s theorem [48]. Nevertheless, the Green’s theorem cannot be applied with no error in the discrete setting and digitized images. In order to handle properly the ellipse parametrization described in [51] instead of [24] in the discrete setting, we propose a pixel-based smooth approximation of the underlying cost energy functional. Our approximation allows us to calculate properly the derivatives of the cost function with respect to the ellipse parameters and is not based on the Green’s theorem also used in [48] for energy minimization.
Definition of the ellipse-based energy
More formally, let Γ be the outer ellipse with parameters \(\left \lbrace \mathbf {x}_{0}, a, b, \theta \right \rbrace \) where \(\mathbf {x}_{0} = \left (x_{0}, y_{0} \right)\) is the center, a and b are the semi major and minor axes respectively, and θ is the angle of rotation. The inner ellipse \({\Gamma }^{\prime }\) is defined as a concentric and coaxial ellipse of Γ such the latter has an area (denoted |Γ|) twice larger than the former: \(|{\Gamma }| = 2|{\Gamma }^{\prime }|\) (see Fig. 5). The factor 2 ensures that the area of the elliptical outer ring is equal to the area of the elliptical inner core. Let us consider a rectangular image patch P containing a potential TMA core associated to a connected component estimated by the detection method in the early stage. The ellipse energy (or cost function) is defined as a normalized image contrast between the two domains \({\Gamma }^{\prime }\subset {\Gamma } \subset P\) where P is a rectangular domain in the image domain Ω which contains a single TMA core [24, 52]:
$$\begin{array}{@{}rcl@{}} J(u, {\Gamma}) & = & \frac{1}{ab} \left(\int_{{\Gamma} \backslash {\Gamma}^{\prime}} u(\mathbf{x}) \ d\mathbf{x} - \int_{{\Gamma}^{\prime}} {u}(\mathbf{x}) \ d\mathbf{x}\right)~~~~~~ \\ & =& \frac{1}{ab} \left(\int_{{\Gamma}} {u}(\mathbf{x}) \ d\mathbf{x} - 2\int_{{\Gamma}^{\prime}} {u}(\mathbf{x}) \ d\mathbf{x}\right). \end{array} $$
(12)
To handle discrete images, the continuous image u defined in (12) can be replaced by its sampled version as follows:
where \(u\left [\mathbf {x}\right ]\) is the discrete sample of u(x) and
denotes the set indicator function such as
if x∈Γ and 0 otherwise. However, there are two major drawbacks while considering this energy function. Firstly, the calculation of the energy gradient is not trivial in the discrete setting since the indicator functions in (13) are piecewise constant which are not differentiable at some points. Secondly, due to sampling effect, brutal switch of the membership of some points from a domain to another may happen just with an infinitesimal change in the ellipse parameters, giving rise to severe numerical instabilities. Smooth approximations of the underlying piecewise constant functions is recommended to overcome both discontinuity and sampling problems. The calculations of partial derivatives of the energy functional is facilitated if we can define a fuzzy membership to avoid abrupt domain switches (see Fig. 8a and b). Our goal is then to build an approximation which favors the computation of the partial derivative of the energy with respect to each ellipse parameter as much as possible. First, we consider the following quadratic form:
$$\begin{array}{@{}rcl@{}} \left\|\mathbf{x} \right\|_{{\Gamma}}^{2} = \left\| \left[ \begin{array}{cc} a^{-1} & 0 \\ 0 & b^{-1} \end{array} \right] \left[ \begin{array}{cc} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{array} \right] \mathbf{x} \ \right\|_{2}^{2}. \end{array} $$
(14)
For a given point x, \(\left \| \mathbf {x} - {\mathbf {x}_0} \right \|_{{{\Gamma }}}\) is a normalized metric between x and the ellipse center x0 induced by the geometry of the ellipse Γ. A pixel x belongs to the interior of the ellipse Γ if and only if \( \left \| \mathbf {x} - {\mathbf {x}_0} \right \|_{{{\Gamma }}}^{2} \leq 1\) since \(\|\mathbf {x} - {\mathbf {x}_0} \|_{{{\Gamma }}}^{2}\) is always positive. The term
can be then expressed by a function of \(\left \| \mathbf {x} - {\mathbf {x}_0} \right \|_{{{\Gamma }}}\) as
. Moreover, we need to find a smooth function which closely approximates
as investigated in [24, 38] and has simple derivative. We realized that the graph of
looks similar to the \(C^{\infty }\) S-shaped logistic curve whose the derivative is easy to compute. Let us consider therefore the following logistic function:
where ε>0 controls the steepness of the curve (see the plot of t↦S
ε
(t) in Fig. 6 for several values of ε). The smaller ε, the closer the curve S
ε
approaches the graph of the indicator function
. Thanks to the property of logistic functions, the derivative of S
ε
can be easily computed as \(S_{\epsilon }'(t)\,=\,-\epsilon ^{-1}\!S_{\epsilon }(t) \left (1 - S_{\epsilon }(t)\right)\). Finally, the energy functional has the following form [24]:
$$ J(u, {\Gamma}) = \frac{1}{ab} {\sum_{\mathbf{x} \in P \cap \mathbb{Z}^{2}}} w_{\epsilon} \left(\|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2}\right)\, {u}\left[\mathbf{x}\right] \, $$
(16)
with \(w_{\epsilon }(t) = S_{\epsilon }(t) - 2S_{\epsilon }(2t) = \frac {1}{1 + e^{\frac {t-1}{\epsilon }}} - \frac {2}{1 + e^{\frac {2t-1}{\epsilon }}}\). For illustration, we present in Fig. 7 the plot of w
ε
whose the term \(w_{\epsilon }\left (\| \mathbf {x} - {\mathbf {x}_0} \|_{{\Gamma }}^{2}\right)\) is nothing else than a smooth approximation of the piecewise constant function
. These weights are very similar to those described in [38] and based on the \(\arctan \) function.
Calculation of partial derivatives
By applying the derivation rules of composite functions, the partial derivatives of the energy with respect to each ellipse parameter \(\left \lbrace \mathbf {x}_{0}, a, b, \theta \right \rbrace \) are given by:
$$\begin{array}{@{}rcl@{}} \left\{ \begin{array}{ll} \!\!\!\frac{\partial J(u, {{\Gamma}})}{\partial x_{0}} & \!\!\!\!\!\! = \frac{1}{ab} \sum_{\mathbf{x}} {u}[\mathbf{x}] w_{\epsilon}' \left(\|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2} \right) \frac{\partial \|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2}}{\partial x_{0}}, \\ \!\!\!\frac{\partial J(u, {{\Gamma}})}{\partial y_{0}} & \!\!\!\!\!\! = \frac{1}{ab} \sum_{\mathbf{x}} u[\mathbf{x}] w_{\epsilon}' \left(\|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2} \right) \frac{\partial \|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2}}{\partial y_{0}}, \\ \!\!\!\frac{\partial J(u, {{\Gamma}})}{\partial \theta} & \!\!\!\!\!\! = \frac{1}{ab} \sum_{\mathbf{x}} u [\mathbf{x}] w_{\epsilon}' \left(\|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2} \right) \frac{\partial \|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2}}{\partial \theta}, \\ \!\!\!\frac{\partial J(u, {{\Gamma}})}{\partial a} & \!\!\!\!\!\! = \frac{1}{ab} \sum_{\mathbf{x}} u [\mathbf{x}] w_{\epsilon}' \left(\|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2} \right) \frac{\partial \|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2}}{\partial a} - \frac{J(u, {{\Gamma}})}{a},~~~ \\ \!\!\!\frac{\partial J(u, {{\Gamma}})}{\partial b} & \!\!\!\!\!\! = \frac{1}{ab} \sum_{\mathbf{x}} u[\mathbf{x}] w_{\epsilon}' \left(\|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2} \right) \frac{\partial \|\mathbf{x} - {\mathbf{x}_0} \|_{{{\Gamma}}}^{2}}{\partial b} - \frac{J(u, {{\Gamma}})}{b},~~~~ \end{array} \right. \end{array} $$
$$ \begin{array}{l l} \text{where}~ w_{\epsilon}'(t) & = \frac{4S_{\epsilon}(2t)\left(1 - S_{\epsilon}(2t) \!\right) - S_{\epsilon}(t)\left(1 - S_{\epsilon}(t)\right)}{\epsilon} \\ \\ & = \frac{1}{\epsilon} \left[ \frac{4e^{\frac{2t-1}{\epsilon}}}{\left(1 + e^{\frac{2t-1}{\epsilon}} \right)^{2}} - \frac{e^{\frac{t-1}{\epsilon}}}{\left(1 + e^{\frac{t-1}{\epsilon}} \right)^{2}}\right] \end{array} $$
and the calculation of partial derivatives of \(\|\mathbf {x} - {\mathbf {x}_0} \|_{{{\Gamma }}}^{2}\) are detailed in the Additional file 1. As depicted in Fig. 8c, for a given parametrization {x0,a,b,θ}, the term \(w_{\epsilon }'\left (\| \mathbf {x} - {\mathbf {x}_0} \|_{{{\Gamma }}}^{2} \right)\) vanishes for most of points x. Thus, the computation of the partial derivatives J(u,Γ) takes account only few points near the ellipse boundaries where \(w_{\epsilon }'\left (\|. - {\mathbf {x}_0} \|_{{{\Gamma }}}^{2} \right)\) is non-zero. Our smooth approximation which is adapted for discrete images produces similar expressions of the partial derivatives of the ellipse energy when comparing with those described in [51] for continuous images. It can be viewed as the expression of the Green’s theorem in the discrete setting and an alternative to the optimization presented in [44].
Multi-ellipse segmentation for multi-tissue core analysis
Let {c
n
}1≤n≤N be the centroids of the connected components of the binary detection map \(I_{_{\text {FA}}}\). In the original image u, we extract a rectangular patch P
n
centered at c
n
with a radius ρ larger than the given tissue core radius r
core
(for example, ρ=2r
core
). Let us define
$$ \Pi_{\rho, \mathbf{c}_{n}} u = \left\lbrace {u}\left[\mathbf{x}\right], \ \left\| \mathbf{x} - \mathbf{c}_{n} \right\|_{\infty} \leq \rho \right\rbrace, $$
(17)
where x=(x,y)∈P
n
, \(\|\mathbf {x}\|_{\infty } = \sup (|x|,|y|)\) and \(\phantom {\dot {i}\!}\Pi _{\rho, \mathbf {c}_{n}} \cdot \) denotes the patch extraction operator with center c
n
and radius ρ. In order to perform a multi-object segmentation, we consider the following multi-ellipse optimization problem:
$$\begin{array}{*{20}l} & \underset{{{{\Gamma}}}_{1}, \hdots, {{{\Gamma}}}_{N}}{\arg\min} \ \sum_{n=1}^{N}\! \left\lbrace\! \frac{1}{a_{n}b_{n}}\! \sum_{\mathbf{x}}\! w_{\epsilon} \! \left(\! \|\mathbf{x} - {\mathbf{x}_0}^{n} \|_{_{{{\Gamma}}_{n}}}^{2} \! \right) \! \Pi_{\rho, \mathbf{c}_{n}} u\left[\mathbf{x}\right] \right\rbrace \\ & \text{subject to} \left({{{\Gamma}}}_{1}, {{{\Gamma}}}_{2}, \hdots, {{{\Gamma}}}_{N} \right) \in \Upsilon \, \end{array} $$
(18)
where \(\left \lbrace {\mathbf {x}_{0,n}}, a_{n}, b_{n}, \theta _{n} \right \rbrace \) are the parameters of the ellipse Γ
n
and Υ is a set of constraints to ensure the ellipses fall into an acceptable range of configurations. In practice, we typically set
$$\begin{array}{*{20}l} \!\!\!\!\!\!\!\!\!\!\!\! \Upsilon & = \big\lbrace \left\| {\mathbf{x}_{0,n}} - {\mathbf{x}_{0,n'}} \right\|_{2} > \rho;~~~~ \\ & \quad \quad \quad \quad \ \left\| {\mathbf{x}_{0,n}} - \mathbf{c}_{n} \right\|_{\infty} \leq \rho_{\max};~~~~ \\ & \quad \quad \quad \quad \ r_{\min} \leq a_{n}, b_{n} \leq r_{\max};~~~~ \\ & \quad \quad \quad \quad \ \theta_{\min} \leq \theta_{n} \leq \theta_{\max} \quad \quad \quad \big\rbrace_{1 \leq n, n'\leq N} \, \end{array} $$
for some predefined values \(\rho _{\max }, r_{\min }, r_{\max }, \theta _{\min }, \theta _{\max }\) set according to the extracted patch positions and the allowed sizes and orientations of tissue cores. The constraint \(\left \| {\mathbf {x}_{0,n}} - {\mathbf {x}_{0,n'}} \right \|_{2} > \rho \) which prevents the distance between two ellipse centers being too close helps to avoid the overlapping of segmented tissue cores. In what follows, we denote \(\mathcal {J}(u, \Gamma _{1}, \ldots, \Gamma _{n})\) the global cost function associated with the optimization problem (18).
By construction, the function \(\mathcal {J}(u, \Gamma _{1}, \ldots, \Gamma _{n})\) is differentiable with respect to \(\left ({{\Gamma }}_{1}, \hdots, {{\Gamma }}_{N} \right)\). The common way to minimize \(\mathcal {J}(u, \Gamma _{1}, \ldots, \Gamma _{n})\) under the constraint set Υ is to use a gradient method whose performance depends on how efficient is the computation of the gradient of \(\mathcal {J}(u, \Gamma _{1}, \ldots, \Gamma _{n})\). Since \(\mathcal {J}(u, \Gamma _{1}, \ldots, \Gamma _{n})\) is a linear combination of separable functions, therefore, the gradient can be simply obtained as:
$$ \nabla \mathcal{J}(u, \Gamma_{1}, \ldots, \Gamma_{n}) = \left(\begin{array}{c} \nabla J (\Pi_{\rho,\mathbf{c}_{1}} u, \Gamma_{1}) \\ \vdots \\ \nabla J (\Pi_{\rho,\mathbf{c}_{N}} u, \Gamma_{N}) \end{array} \right), $$
(19)
where
$$J (\Pi_{\rho,\mathbf{c}_{n}} u, \Gamma_{n}) = \frac{1}{a_{n}b_{n}} \sum_{\mathbf{x} \in \Gamma_{n}}\! w_{\epsilon} \left(\|\mathbf{x} - {\mathbf{x}_0}^{n} \|_{_{{{\Gamma}}_{n}}}^{2} \right) \Pi_{\rho,\mathbf{c}_{n}} {u}[\mathbf{x}] $$
and the expressions of its partial derivatives are given in (17).
The result of the multi-ellipse optimization problem (18) is a set of ellipses \(\left \lbrace {{\Gamma }}_{n} \right \rbrace _{1 \leq n \leq N}\phantom {\dot {i}\!}\) which fits the objects located in the regions of interest \(\lbrace \Pi _{\rho,\mathbf {c}_{n}} u \rbrace _{1 \leq n \leq N}\phantom {\dot {i}\!}\). Furthermore, given the major axes of these ellipses and the TMA core radius r
core
, we discard the tiny, giant and flattened ellipses and we keep those which are most similar to the expected tissue cores. The center of the selected ellipse allows us to determine the position of the recognized TMA core. This reference position will be used to determine the array coordinates of the corresponding tissue core. In the following, we denote \(\mathcal {X}_{0} = \lbrace {\mathbf {x}_{0,n}} \rbrace _{n \in \lbrace 1, \hdots, N \rbrace }\) as the set of centers of the N reliable and selected ellipses.
Estimation of array coordinate and TMA core positions
An ideal TMA is the one which has tissue cores perfectly aligned in both horizontal and vertical directions and equally spaced according to a regular square grid. The array coordinate \(\mathbf {p} = (k, l) \in \mathbb {Z}^{2}\) of a core can be simply obtained by drawing two orthogonal lines crossed at the considered core position. However, due to the deformation of the design TMA grid, the lines passing through tissue cores and their nearest neighbors may be slightly inclined with respect to the horizontal or vertical axes. Moreover, the direction of these lines may have a large spectrum of variations which makes more challenging the tracking of tissue cores over a given direction. To deal with this deformation, existing TMA de-arraying methods use usually distance-and-angle-based criteria for the purpose of defining the neighborhood of TMA cores. Although this approach estimates robustly the average core-to-core distance and the two principal directions of the deformed core grid, it may fail for some well-detected cores whose the position is strongly distorted with respect to their neighbors. In order to avoid this failure, we introduce an algorithm for estimating iteratively the deformation of the TMA grid in a way that the grid which is warped by the estimated deformation at an iteration gets closer to the observed TMA grid. To this end, we assume that the deformation of the TMA grid can be decomposed by linear and non-linear parts. Under this assumption, we estimate the linear part of the deformation by defining an oblique grid (affine warping) which is derived from the detected core positions as the initialization of the warped grid (see Fig. 9). The latter is used to find nearby cores that will be taken into account to compute an estimator of the grid deformation by using the thin-plate interpolation [25] if we do an analogy with material deformation.
Estimation of the linear deformation
Our goal is to approximate the distorted TMA grid Λ (which is observed partially with the set of point \(\mathcal {X}_{0}\)) by an oblique grid Λ0 which minimizes the distance between them in the way that the deformation of the grid is approximated by a 2D affine transform. For this purpose, we consider the set \({\mathcal {C}_0}\) of core pairs whose each pair \(\phantom {\dot {i}\!}({\mathbf {x}_{0,n}}, {\mathbf {x}_{0,n'}})\) is formed by an element of \(\mathcal {X}_{0}\) and one of its four nearest neighbors with respect to the Euclidean distance
$$ {\mathcal{C}_0} = \left\lbrace ({\mathbf{x}_{0,n}}, {\mathbf{x}_{0,n'}}) \in \mathcal{X}_{0} \times \mathcal{X}_{0},\, {\mathbf{x}_{0,n'}} \in \mathcal{N}({\mathbf{x}_{0,n}}) \right\rbrace, $$
where \(\mathcal {N}({\mathbf {x}_{0,n}})\) denotes the 4-neighborhood of x0,n. To estimate the average core-to-core distance \(\bar {d}_{cc}\), we compute the trimmed mean (denoted TM) of the length of the segment defined by the pair \(\phantom {\dot {i}\!}({\mathbf {x}_{0,n}}, {\mathbf {x}_{0,n'}})\) of \({\mathcal {C}_0}\) by discarding the most extreme values (typically 30%):
$$ {\bar{d}} = \text{TM}_{30\%} \left\lbrace \| {\mathbf{x}_{0,n}} - {\mathbf{x}_{0,n'}}\|_{2} \right\rbrace_{({\mathbf{x}_{0,n}}, {\mathbf{x}_{0,n'}}) \in {\mathcal{C}_0}} \ . $$
(20)
Let \(\text {ang}({\mathbf {x}_{0,n}}, {\mathbf {x}_{0,n'}})\phantom {\dot {i}\!}\) be the angle between the line passing through \(({\mathbf {x}_{0,n}}, {\mathbf {x}_{0,n'}})\phantom {\dot {i}\!}\) and the horizontal axis such that \(-0.25{\pi } \leq \text {ang}({\mathbf {x}_{0,n}}, {\mathbf {x}_{0,n'}}) \leq {0.75\pi }\phantom {\dot {i}\!}\). By analogy, we define the two principal angles of the deformed TMA grid as follows:
$$\begin{array}{*{20}l} {\bar{\alpha}} & = \text{TM}_{30\%} \left\lbrace \text{ang}({\mathbf{x}_{0,n}}, {\mathbf{x}_{0,n'}}) \leq \frac{\pi}{4} \right\rbrace, \\ {\bar{\beta}} & = \text{TM}_{30\%} \left\lbrace \text{ang}({\mathbf{x}_{0,n}}, {\mathbf{x}_{0,n'}}) \geq \frac{\pi}{4} \right\rbrace. \end{array} $$
Finally, we denote \(\hat {\mathbf {t}}\) as the global translation of the distorted TMA grid with respect to the origin that minimizes the distance between the set \(\mathcal {X}_{0}\) and the linearly-estimated grid Λ0
$$ \hat{\mathbf{t}} = \arg\min_{\mathbf{t}} \sum_{n = 1}^{N} \min_{\mathbf{p} \in \mathbb{Z}^{2}} \| \mathbf{t} + \mathcal{F}(\mathbf{p}) - {\mathbf{x}_{0,n}}\|_{2}^{2} $$
(21)
where \(\mathcal {F}\) maps each array coordinates \(\mathbf {p} \in \mathbb {Z}^{2}\) to a position of Λ0 corrected by \(\hat {\mathbf {t}}\) and
$$ \mathcal{F}(\mathbf{p}) = \bar{d} \ \underset{\mathbf{M}_{\bar{\alpha}, \bar{\beta}}}{\underbrace{\left[ \!\! \begin{array}{c c} \cos\bar{\alpha} & \cos\bar{\beta} \\ \sin\bar{\alpha} & \sin\bar{\beta} \end{array} \!\! \right] }} \ \mathbf{p} \ . $$
(22)
Note that \(\mathbf {M}_{\bar {\alpha }, \bar {\beta }}\) is a change-of-basis matrix of unit column vectors and \(\bar {d}\) is a scaling factor which transforms array coordinates (elements of \(\mathbb {Z}^{2}\)) to real spatial positions (in \(\Omega \subset \mathbb {R}^{2}\)). The resulting oblique grid is parametrized with four parameters \(\{\bar {d}, \bar {\alpha }, \bar {\beta }, \hat {t}\}\) and represents the affine part of the grid deformation. We thus arrive at the affine mapping function: \(\mathcal {A}(\mathbf {p}) = \hat {\mathbf {t}} + \mathcal {F}(\mathbf {p}) \in \Lambda _{0}\). The oblique grid Λ0 will serve as initialization to estimate of the non-linear deformation of the grid. Figure 9 illustrates an example showing the oblique grid obtained from a given set of points as well as its estimated parameters.
Thin-plate-based estimation of the deformation
Let Λ⋆ be the ideal design TMA grid with (0,0) as origin and d the ideal distance between two neighboring cores along the horizontal and vertical axes. The mapping is then defined as:
$$ {\mathbf{y}_{\mathbf{p}}^{\star}} = d \mathbf{p} \in {{\Lambda^{\star}}}, \forall \mathbf{p} \in \mathbb{Z}^{2}. $$
The deformation \({\mathcal {D}}\) maps each point \({\mathbf {y}_{\mathbf {p}}^{\star }} \in {{\Lambda ^{\star }}}\) onto a point \({\mathbf {y}_{\mathbf {p}}} = \mathcal {D}({\mathbf {y}_{\mathbf {p}}^{\star }})\) in the distorted grid Λ. In order to estimate the deformation \({\mathcal {D}}\) at all points of the grid Λ⋆, we aim at approximating this set from the observed set \(\mathcal {X}_{0} = \{{\mathbf {x}_{0,n}}\} \) by using the thin-plate splines as an interpolant. Indeed, given a set of points \({\mathcal {D}}^{-1} \mathcal {X}_{0} = \left \lbrace {\mathcal {D}}^{-1} \left ({\mathbf {x}_{0,n}} \right) \right \rbrace _{n \in \lbrace 1, \ldots, N \rbrace }\), the coefficients of the interpolating thin-plate splines are the minimizers of a quadratic function which is the first approximation of the bending energy of the mapping from \({\mathcal {D}}^{-1} \mathcal {X}_{0}\) to the set of target points \(\mathcal {X}_{0}\) (see [25]). Nevertheless, unlike the usual framework [25], the correspondence between the two sets of points is not established, that is \({\mathcal {D}}^{-1} \mathcal {X}_{0}\) is unknown. Instead of investigating a matching method to determine \({\mathcal {D}}^{-1} \mathcal {X}_{0}\), we propose to build a sequence of grids {Λ(m)}m≥0 which evolves iteratively to fit \(\mathcal {X}_{0}\). We initialize this sequence with the oblique grid Λ(0)=Λ0 previously computed. The linear approximation of \({\mathcal {D}}\) is then as follows:
$$\begin{array}{@{}rcl@{}} {\mathbf{y}_{\mathbf{p}}}^{(0)}&=& {\mathcal{D}}^{(0)}({\mathbf{y}_{\mathbf{p}}^{\star}}) \\ &=& \hat{\mathbf{t}} + \frac{\ \bar{d}\ }{\ {d} \ }\!\left[ \!\! \begin{array}{c c} \cos\bar{\alpha} & \cos\bar{\beta} \\ \sin\bar{\alpha} & \sin\bar{\beta} \end{array} \!\! \right]\! \mathbf{p}. \end{array} $$
(23)
At iteration m, a core position \({\mathbf {x}_{0,n}} \in \mathcal {X}_{0}\) is associated to a position \({\mathbf {y}_{\mathbf {p}}^{\star }}\) if the former is located within a radius δ from \({\mathbf {y}_{\mathbf {p}}}^{{{(m)}}} = {\mathcal {D}}^{{{(m)}}}({\mathbf {y}_{\mathbf {p}}^{\star }})\). Pairs of associated positions establish therefore the correspondence between the ideal grid Λ⋆ and the set of observed point \(\mathcal {X}_{0}\). We also note that all the positions of Λ⋆ do not have a corresponding position in \(\mathcal {X}_{0}\) as shown in Fig. 10 mainly because the cardinality of these sets are not the same. Let \(\mathcal {P}^{{{(m)}}}\) be the set of pairs of associated positions:
$$ \mathcal{P}^{{{(m)}}} = \left\lbrace \left({\mathbf{x}_{0,n}}, {\mathbf{y}_{\mathbf{p}}^{\star}} \right),\, \| {\mathcal{D}}^{{{(m)}}}({\mathbf{y}_{\mathbf{p}}^{\star}}) - {\mathbf{x}_{0,n}} \|_{2} \leq \delta \right\rbrace.~~~~ $$
(24)
Assume that N(m) is the number of elements of \(\mathcal {P}^{{{(m)}}}\). The objective is to estimate the deformation \(\mathcal {D}^{{{(m)}}}\) from the set of N(m) associated pairs \(({\mathbf {x}_{0,n}}, {\mathbf {y}_{\mathbf {p}}^{\star }})\). According to [25], we define the Gram’s matrix \(\left \lbrace \mathbf {K}^{{{(m)}}}_{n, n'}\right \rbrace _{1 \leq n, n' \leq N^{{{(m)}}}}\) as follows:
$$\begin{array}{@{}rcl@{}} \mathbf{K}^{{{(m)}}}_{n, n'} & = \left\| {\mathbf{x}_{0,n}} - {\mathbf{x}_{0,n'}} \right\|_{2}^{2} \log \left\| {\mathbf{x}_{0,n}} - {\mathbf{x}_{0,n'}} \right\|_{2}^{2}, \end{array} $$
(25)
and the additional matrices as:
$$\begin{array}{@{}rcl@{}} \mathbf{Y}^{{{(m)}}} & = &\left[ \begin{array}{cccc} 1 & 1 & \hdots & 1 \\ \mathbf{y}_{\mathbf{p}_{1}} & \mathbf{y}_{\mathbf{p}_{2}} & \hdots & \mathbf{y}_{\mathbf{p}_{N^{(m)}}} \end{array} \right], \end{array} $$
(26)
$$\begin{array}{@{}rcl@{}} \mathbf{X}^{{{(m)}}} & = &\left[ \begin{array}{ccccccc} \mathbf{x}_{0,1} & \mathbf{x}_{0,2} & \hdots & \mathbf{x}_{0,N^{{{(m)}}}} & 0 & 0 & 0\\ \end{array} \right], \end{array} $$
(28)
$$\begin{array}{@{}rcl@{}} \mathbf{W}^{{{(m)}}} & = & \left(\left(\mathbf{L}^{{{(m)}}}\right)^{-1} \left(\mathbf{X}^{{{(m)}}}\right)^{\top}\right)^{\top}, \end{array} $$
(29)
and \(\mathbf {W}^{{{(m)}}} \,=\, \left (\!\mathbf {w}_{1}^{{{(m)}}}, \mathbf {w}_{2}^{{{(m)}}}, \hdots, \mathbf {w}_{N^{{{(m)}}}}^{{{(m)}}}\!\right)^{\top }\). By using the entries of the matrix W(m), the estimators of the deformation \({\mathcal {D}}\) and of the grid Λ at the next iteration m+1 are therefore defined as:
$$\begin{array}{@{}rcl@{}} {\mathbf{y}_{\mathbf{p}}}^{(m + 1)} & = & \mathbf{w}^{^{{{(m)}}}}_{_{N^{{{(m)}}} + 1}} + \left[\mathbf{w}^{^{{{(m)}}}}_{_{N^{{{(m)}}} + 2}} \mathbf{w}^{^{{{(m)}}}}_{_{N^{{{(m)}}} + 3}} \right] {{\mathbf{y}_{\mathbf{p}}^{\star}}}\\ && + \sum_{n=1}^{N^{{{(m)}}}} \mathbf{w}^{{{(m)}}}_{n} \left(\| {\mathbf{x}_{0,n}} - {{\mathbf{y}_{\mathbf{p}}^{\star}}}\|_{2}^{2} \log \| {\mathbf{x}_{0,n}} - {{\mathbf{y}_{\mathbf{p}}^{\star}}}\|_{2}^{2}\right). \end{array} $$
(30)
This iterative scheme will be stopped at the iteration m∗=m if there are no change between Λ(m) and Λ(m+1). At convergence, the row and column coordinates of a detected cores of position \({\mathbf {x}_{0,n}} \in \mathcal {X}_{0}\) is simply given by:
$$ \hat{\mathbf{p}} = \underset{\mathbf{p} \in \mathbb{Z}^{2}}{\arg\min} \big\| {\mathbf{x}_{0,n}} - \mathcal{D}^{(m^{*})}({\mathbf{y}_{\mathbf{p}}^{\star}})\big\|_{2}^{2} \ . $$
(31)
Moreover, since the grid \({{\Lambda }}^{(m^{*})}\) is an estimator of the deformed TMA grid Λ which is partially observed in \(\mathcal {X}_{0}\), it can be used as approximated positions to recognize tissue cores which are missed during detection and segmentation processes. Indeed, to refine de-arraying result, we perform another multi-ellipse optimization at the position of remaining nodes of the grid \({{\Lambda }}^{(m^{*})}\). If there are ellipses that meet the size and the roundness criteria of standard cores, we add them to the list of detected core position and adjust the coefficients of the thin-plate splines according to the new list. An example of TMA de-arraying is depicted in Fig. 14 showing the gain of our method in term of tissue core detection.