Part-to-whole Registration of Histology and MRI using Shape Elements
Jonas Pichat Juan Eugenio Iglesias Sotiris Nousias Tarek Yousry Sebastien Ourselin Marc Modat







Part-to-whole Registration of Histology and MRI using Shape Elements
    Jonas Pichat^1
Juan Eugenio Iglesias^1
Sotiris Nousias^1
Tarek Yousry^2
Sébastien Ourselin^1,3    Marc Modat^1

^1Translational Imaging Group, CMIC, University College London, UK

^2Department of Brain Repair & Rehabilitation, UCL Institute of Neurology, UK

^3Wellcome / EPSRC Centre for Interventional and Surgical Sciences, UCL, UK

jonas.pichat.13@ucl.ac.uk






    
===================================================================================================================================================================================================================================================================================================================================================================================

empty




Image registration between histology and magnetic resonance imaging (MRI) is a challenging task due to differences in structural content and contrast.  Too thick and wide specimens cannot be processed all at once and must be cut into smaller pieces. This dramatically increases the complexity of the problem, since each piece should be individually and manually pre-aligned. To the best of our knowledge, no automatic method can reliably locate such piece of tissue within its respective whole in the MRI slice, and align it without any prior information. We propose here a novel automatic approach to the joint problem of multimodal registration between histology and MRI, when only a fraction of tissue is available from histology. The approach relies on the representation of images using their level lines so as to reach contrast invariance. Shape elements obtained via the extraction of bitangents are encoded in a projective-invariant manner, which permits the identification of common pieces of curves between two images. We evaluated the approach on human brain histology and compared resulting alignments against manually annotated ground truths. Considering the complexity of the brain folding patterns, preliminary results are promising and suggest the use of characteristic and meaningful shape elements for improved robustness and efficiency.




















§ INTRODUCTION








Histology is concerned with the various methods of microscopic examination of a thin tissue section. Cutting through a specimen permits the investigation of its internal topography and the observation of complex differentiated structures through staining. 

MRI constitutes an invaluable resource for routine, accurate, non-invasive study of biological structures in three dimensions. Relative to histology, MRI avoids irreversible damage and distortions induced by processing, cutting, mounting and staining during the histological preparation. However, resolution-wise, it is outperformed by histology. 


One of the many benefits of combining histology and MRI is to confirm non-invasive measures with baseline information on the actual properties of tissues <cit.> by accessing simultaneously the chemical and cellular information of the former and the rich structural information of the latter.


Such combination relies on image registration and this can be achieved using iconic (intensity-based) <cit.> or geometric (feature-based) <cit.> approaches. Unfortunately, the extraction and manipulation of meaningful information from histology and clinical images is a very complicated task because each modality has, by nature, its own features and there does not always exist a mapping between their constituents: local intensity mappings are non-linear and images exhibit different structures—which is also a reason why intensity-based methods tend to get trapped in local optima.
Hence, classical feature description methods, such as SIFT <cit.>, fail to match features <cit.>. Incidentally, manual extraction of landmarks may remain the safest way to establish correspondences across modalities <cit.>.



Besides, it is common for histopathology laboratories to receive tissue samples that are: (P1) too wide or (P2) too thick, to be processed as they are. The sample is therefore cut into separate sub-blocks, each of which is processed individually. If no scan of each sub-block is available (unlike in <cit.> for example), one must keep track on which part of the sample each sub-block corresponds to and use that knowledge to initialise the registration of histological slices with the clinical image, or manually align them. As for problem (P2), attempts at using similarity measures have been made to initialise registrations, but those are ambiguous and rely on absolute measures rather than relative ones <cit.>. On that matter, it was shown in <cit.> that direct comparison of images from different modalities is non-trivial, and fails to reliably determine slice correspondences. To the best of our knowledge, no automatic method to address (P1) (see Fig. <ref>) has been proposed in the literature.








 §.§ Related work


Regarding geometric approaches, one possible strategy to align histology and clinical imaging is to simplify the images into their contours, so as to come down to a monomodal registration problem and use the shape information provided by the external boundaries. In <cit.>, contours from both histology and slices from a rat brain atlas were extracted via thresholding and represented using B-splines. Then, they were described by means of sets of affine invariants constructed from the sequence of area patches bounded by the contour and the line connecting two consecutive inflections. In <cit.>, Curvature Scale Space <cit.> was used for the registration of whole-slide images of histological sections in order to represent shape (the tissue boundary) at various scales. In <cit.>, curvature maps at different scales were used to match boundaries of full brain MRI extracted via an active contour algorithm. The main weaknesses of active contours are the number of parameters and the sensitivity to initialisation. 

An alternative to using a single contour was proposed by morphologists, observing that level lines (the boundaries of level sets) provide a complete, contrast-invariant representation of images. Furthermore, level lines fit the boundaries of structures and sub-structures of objects very well. Then, given two images, the problem is to retrieve all the level lines that are common to both images; this is however feasible only if curves have been appropriately simplified (smoothed) <cit.> (p.95). Like in <cit.>, smooth pieces of level lines (the shape elements <cit.>) can be encoded to represent shape locally in e.g., an affine-invariant manner <cit.>. The comparison of the resulting canonical curves then permits to identify portions of level lines common to two images. 

Problem (P1) being multimodal and fractional by nature, it seems natural to formulate a solution that involves contrast- and geometric-invariance, as well as locality. 


Here, we present a novel approach to (P1) based on: (i) representing both histology and MRI images using their level lines <cit.>. This allows to reach contrast invariance and to consider implicitly several structural layers of the images—as opposed to relying solely on the outer boundaries of tissues. From there, characteristic shape elements can be extracted locally along the level lines via their bitangents (<ref>). (ii) Representing those elements in a projective invariant manner (<ref>) as introduced by Rothwell in <cit.>, so as to be robust to some non-linear deformations that tissues undergo during the histological process. Combining the two procedures permits the partial matching of shape elements regardless of the orientation of the tissue on glass slides. Registration is then obtained as a result of shape recognition (<ref>). 









 §.§ Contributions

	

    
  * We address the joint problem of multimodal registration between a fraction of histology and its whole in an MRI slice as a result of shape recognition using portions of level lines. 
    
  * We introduce an efficient refinement of bitangents via ellipses.
    
  * We extend Rothwell's framework to bitangents crossing the level lines and compare the resulting canonical curves using the Fréchet distance. 



















§ PREPROCESSING


We used two standard preprocessing steps: first, smoothing, in order to simplify the image, preserve the shape of the tissue, remove unnecessary details and obtain smooth level lines (Fig. <ref>a); then, intensity correction, in order to account for inhomogeneities of the field in MR images (Fig. <ref>c) or illumination in histology. 

Smoothing is based on Affine Morphological Scale Space (AMSS) <cit.>. It is governed by the partial differential equation: ∂ u/∂ t = |Du|(u)^1/3 where u is the image, |Du| is the gradient of the image, (u) is the curvature of the level line and t is a scale parameter. AMSS smoothes homogeneous regions but enhances tissue boundaries. The sequence of updates necessary to its computation follows that presented in  <cit.> (equations of 2.3).

Image intensities correction relies on surface fitting <cit.>: the low-frequency bias of an image can be estimated using an adequate basis of smooth and orthogonal polynomial functions. It then comes down to solving the least square problem A𝖼=𝖻, where 𝖻∈ℝ^N is the vector of all the pixels values, and 𝖼 the coefficients of one linear combination of basis functions. A ∈ℝ^N× (n+1)(m+1) is the matrix of the system: its k-th row is the vectorised outer product Φ(x_k)⊗Φ(y_k) with Φ(x_k) = [P_0(x_k), P_1(x_k), …, P_m(x_k)]^T and Φ(y_k) = [P_0(y_k), P_1(y_k), …, P_n(y_k)] for pixel k≤ N. P_i(.) denotes a certain 1D polynomial of degree i.  Degrees m and n are usually taken small so as not to overfit the image intensities. The left inverse of A (it is full rank) gives the bias image, and correction is straightforward (Fig. <ref>b).

















§ FINDING BITANGENTS



Characteristic shape elements are extracted by means of bitangents of level lines. Bitangents are identified via the tangent space (<ref>) and each one is refined using two ellipses fitted in the neighbourhood of estimated bitangent points (<ref>). Since two ellipses have at most four bitangents (<ref>), one needs to be singled out which corresponds to the refined bitangent of the level line (<ref>).

In the following, a bitangent point is one of the two points where a bitangent is in contact with the level line. The length of a bitangent is defined as the number of inflections of the portion of level line that it covers. As a result, a short bitangent refers to a bitangent that covers portions with exactly two inflections and a long one, more than two.








 §.§ Dual curve



Let ℒ be a Jordan curve (level lines are plane simple curves, though closedness is not guaranteed for all of them in practice). Duality is defined as the polarity that sends any point to a line and vice versa. The image of a point with parameter t=t_0 is the line:

    ux(t_0) + vy(t_0) + 1 = 0.

If the parameter t covers the whole range of definition, the resulting set of straight lines is the envelope of ℒ: the dual ℒ^* of ℒ is the set of its tangent lines. A parametrisation of ℒ^* in homogeneous coordinates can be obtained from (<ref>) by differentiation  w.r.t the parameter t and elimination. This yields u=-ẏ(t)/ẋ(t)y(t) - ẏ(t)x(t) and v =ẋ(t)/ẋ(t)y(t) - ẏ(t)x(t) with x,y≠ 0 (dot notation is used for differentiation).

Dual curves feature the following properties: an inflection of ℒ maps onto a cusp of the dual, and two points sharing a common tangent map onto a double point of the dual curve. More generally, a set of n points sharing a common tangent line maps onto a point of multiplicity n of the dual curve. Finding the bitangents of ℒ is therefore equivalent to finding self-intersections of the polygonal curve ℒ^* (Fig. <ref>). To that end, we used the Bentley-Ottmann algorithm <cit.>, which is a line sweep algorithm that reports all intersections among line segments in the plane.








 §.§ Refining bitangents locations


The refinement of bitangents is preferable: since the slopes of tangents vary substantially in portions of high curvature, the lengths of segments of the dual curve increase on portions where a self-intersection may happen. The evaluation of that double point thus degrades, which directly affects the estimation of bitangents.













  §.§.§ Ellipse fitting


In order to cope with bitangent errors, we propose to refine their locations by fitting ellipses <cit.> around estimated bitangent points. This allows skipping the rotation part prior to the quadratic fitting in <cit.>. Beforehand, bitangents lying on almost straight edges of the level lines are removed by looking at the residual of a line fit on the portions bounded by the two bitangent points. This is intended to avoid the degenerate case of fitting an ellipse to a nearly straight line.

Let F be a general conic. It is defined as the set of points such that:

    F(𝖺,𝐱) = 𝖺.𝐱 = ax^2 + by^2 + cxy + dx + ey + f = 0,


where 𝖺=[abcdef]^T and 𝐱=[x^2y^2xyxy1]^T.

The constrained least square problem we wish to solve here is: _𝖺 = 𝖺^TS𝖺 subject to 𝖺^TC𝖺=1, where S=D^TD is the scatter matrix, D is the design matrix, made of the N points to be fitted and C is the constraint matrix which expresses the constraint 4ac - b^2 = 1 on the conic parameters to make it an ellipse. This translates in a [6× 6] matrix where C_22=-1 and  C_31=C_13=2, the rest being zeros. 

This yields the generalised eigenvalue problem (GEP):

    S𝖺 = λ C𝖺.


The ellipse coefficients, 𝖺 are the elements of the eigenvector that corresponds to the only positive eigenvalue. Although the impact of S being nearly singular and C being singular on the stability of the eigenvalues computation is discussed in <cit.>, we did not encounter any problem in our experiments. 
















  §.§.§ Bitangents of ellipses


The main goal of this section is to compute the bitangents of two ellipses efficiently. This is achieved by transforming a system of two polynomial equations into a polynomial eigenvalue problem, and for further performance, into a generalised eigenvalue problem.

Let us consider two ellipses, E_1(𝖺_1,𝐱) and E_2(𝖺_2,𝐱) defined by bivariate quadratic polynomials, like in (<ref>). The tangent line, Ty=ux+v to say E_1, is the line that intersects E_1 at exactly one point. By substitution, one gets a degree 2 polynomial in x, which has a single root if and only if its discriminant, Δ(α_1,𝗎)=0. When considering the tangent to both ellipses, this gives a system of n=2 polynomial equations in unknowns u,v: 

    (s1)
    α_11u^2 + α_12v^2 + α_13uv + α_14u + α_15v + α_16 = 0
    α_21u^2 + α_22v^2 + α_23uv + α_24u + α_25v + α_26 = 0,


α_i1 = e_i^2 - 4c_if_iα_i2 = b_i^2 - 4a_ic_iα_i3 = 4c_id_i-2b_id_i, α_i4 = 2d_ie_i-2b_if_iα_i5 = 2b_id_i-4a_ie_i α_i6 = d_i^2-4a_if_i, i={1,2}.

To start with, u is hidden in the coefficient field; (s1) becomes a system of two equations f_1(u,v) and f_2(u,v) in one variable v and coefficients from ℝ[u] i.e. f_1,f_2 ∈ (ℝ[u])[v]. The degrees of these two equations are d_1=d_2=2.

Homogenising (s1) using a new variable w gives (s2), a system of two homogeneous polynomial equations F_1(v,w) and F_2(v,w) in two unknowns v,w:

    (s2)
    0.88
    α_11u^2 + α_12v^2 + α_13uv + α_14uw + α_15vw + α_16w^2 = 0
    0.88
    α_21u^2 + α_22v^2 + α_23uv + α_24uw + α_25vw + α_26w^2 = 0,


The total degree d=∑_i=1^n(d_i-1)+1 equals 3. This gives the set 𝒮 of n+d-1d=4 possible monomials ω^δ=v^δ_2w^δ_3 in variables v,w of total degree d i.e., such that |δ|=∑_i=2^3δ_i=3: 𝒮={v^3,v^2w, vw^2, w^3}. The set 𝒮 can be partitioned into two subsets according to a modified Macaulay-based method <cit.>:

    [ 𝒮_1={ω^δ: |δ|=3, v^d_1|ω^δ},; 𝒮_2={ω^δ: |δ|=3, w^d_2|ω^δ}. ]

In other words, 𝒮_1 (resp. 𝒮_2) is the set of monomials of total degree 3 that can be divided by v^2 (resp. w^2). This gives 𝒮_1={v^3, v^2w} and 𝒮_2={vw^2, w^3}, from which the extended set of four polynomial equations: vF_1 = 0, wF_1 = 0, vF_2 = 0 and wF_2 = 0 can be derived.

After dehomogenisation (by setting w=1), the extended system can be rewritten as a polynomial eigenvalue problem (PEP):

    𝖢(u)𝗏=0,

where 𝗏=[v^3v^2v1]^T and

𝖢(u)=
[                 α_12         α_13u + α_15 α_11u^2+α_14u + α_16                    0;                    0                 α_12         α_13u + α_15 α_11u^2+α_14u + α_16;                 α_22         α_23u + α_25 α_21u^2+α_24u + α_26                    0;                    0                 α_22         α_23u + α_25 α_21u^2+α_24u + α_26 ].

Non-trivial solutions to (<ref>) are the roots of (𝖢), which gives up to 4 real solutions for u. 

For each one of them e.g., u_1, the corresponding singular value decomposition has the form: 𝖢(u_1)=𝖴𝖲𝖵^T, where the solution vector [𝗏_1𝗏_2𝗏_3𝗏_4]^T is the column of 𝖵 that corresponds to the smallest singular value. The particular solution v_1 associated with u_1 is e.g. 𝗏_3/𝗏_4, meaning that one bitangent is parametrised by T_1y=u_1x+v_1.

For the sake of completeness, the PEP (<ref>) can be further transformed into a GEP by first rewriting it as:

    0.84
    (
    [    0    0 α_11    0;    0    0    0 α_11;    0    0 α_21    0;    0    0    0 α_21 ]_𝖢_2 u^2 + 
    [    0 α_13 α_14    0;    0    0 α_13 α_14;    0 α_23 α_24    0;    0    0 α_23 α_24 ]_𝖢_1 u +
    [ α_12 α_15 α_16    0;    0 α_12 α_15 α_16; α_22 α_25 α_26    0;    0 α_22 α_25 α_26 ]_𝖢_0)𝗏 = 0,

which is equivalent to the GEP:

    𝖠𝗒=u𝖡𝗒,

with 𝖠=
[  0_4  I_4; -𝖢_0 -𝖢_1 ] and 𝖡=
[ I_4 0_4; 0_4 𝖢_2 ], 0_4 and I_4 being the [4× 4] zero and identity matrices, and 𝗒=
[  𝗏; u𝗏 ]=[𝗒_1𝗒_2…𝗒_8]^T. A particular solution v_1 is e.g., the quotient 𝗒_3/𝗒_4 (or equivalently √(𝗒_1/𝗒_4)) from the eigenvector associated with eigenvalue u_1. 

Note that the resolution of (<ref>) is two orders of magnitude faster compared to (<ref>) using linear algebra packages.

Lastly, when the two ellipses E_1 and E_2 intersect in two points, two out of the four eigenvalues obtained for u are complex. These correspond to the two internal bitangents: in that case, ellipses have only two external bitangents associated with the other two real eigenvalues. It is also worth noting that, when they exist, internal bitangents are associated with the extremal (real) eigenvalues.









  §.§.§ Selecting one bitangent


In this section, we identify the only bitangent of E_1 and E_2 that is also a bitangent of ℒ (Fig. <ref>)—referred to as the usable bitangent.

Let us consider: (i) bitangents directed from E_1 to E_2, (ii) E_1 is oriented positively and (iii) Δ is its left-most vertical tangent. Bitangents of E_1 can be cyclically ordered by considering independently the tangents below (in blue in Fig. <ref> Left), and above (in red) it, and sorting them by decreasing y-intercept with the ellipse's left-most tangent, Δ. This holds for cases where an ellipse lies above (resp. below) all of the bitangents. Lemma 1 in <cit.> states that the resulting cyclic order of the bitangent directions is 𝒞: [LL,LR,RL,RR] (L and R stand for left and right and refer to the locations of an ellipse relative to a bitangent). 

Four possible cases arise: (c1) E_2 stands to the right of E_1, (c2) is above E_1 intersecting Δ, (c3) is to the left of E_1, and (c4) is below E_1 intersecting Δ. For each case, the first bitangent encountered starting from Δ, counter-clockwise, has type LL, RR, RL and LR respectively; the next up to three bitangents for each case have their types deduced from the positive cyclic order 𝒞.



Now in order to select the usable bitangent, one has to rely on the geometry of the level line ℒ. Let us define the unit curvature vector k, at every point along ℒ as the vector pointing toward the centre of the osculating circle: k=κn=⟨ k_x, k_y ⟩, where κ is the scalar curvature and n is the normal (it is colinear to the gradient of the image along ℒ and directed toward the inside of the clockwise-oriented closed curve here). The orientation of k allows differentiating otherwise ambiguous situations; for example, two pairs of ellipses (E_1, E_2) and (E_1, E_3), all of them fitting portions with same curvature and satisfying the configuration of case (c1), can be associated with a different type of usable bitangent, RR and RL respectively. This happens when k_1 and k_3 have opposite sense, while k_1 and k_2 have the same. In the following, positiveness is defined for (c1) and (c3) as k_y>0 and as k_x>0 for (c2) and (c4), and is denoted with the superscript (+).




From there we define four patterns: (p1) (k_1^(+),k_2^(+)), (p2) (k_1^(+),k_2^(-)), (p3) (k_1^(-),k_2^(+)) and (p4) (k_1^(-),k_2^(-)). In cases (c1) and (c2), they correspond to the usable bitangent type LL, LR, RL, RR respectively. Conversely, in cases (c3) and (c4), they correspond to the type RR, RL, LR, LL respectively. Since there is a one to one correspondence between the four bitangents and the four types, it only requires identifying one of four patterns (p) and one of four cases (c) to pick the usable bitangent parameters.









We also extend the mapping to intersecting ellipses (Fig. <ref> Middle) by observing that the cyclic order of bitangents is of the form [T_e, T_i, T_i, T_e] (subscripts e and i stand for external and internal). Since only external bitangents exist in the case where E_1 and E_2 intersect in two points (<ref>), we are left with the cyclic order [T_e, _, _, T_e].

Bitangent points are straightforward to obtain for E_1 and E_2 by substitution of the tangent equation in the ellipses equations. Finally, we select the point of ℒ that is the closest to an ellipse bitangent point. Note that once all bitangents are refined, some bitangent points may collapse to similar locations. In order to reduce ineffective redundancy, only one bitangent out of those that have their end points close to each other is kept <cit.>.









§ PROJECTIVE SHAPE REPRESENTATION


We now have a set of refined bitangents. Let us consider one bitangent and its endpoints b_1 and b_2. In order to encode the shape of a portion of (oriented) level line ℒ_r = ℒ[b_1,b_2] (assuming b_1 comes before b_2) in a projective invariant manner (as opposed to affine invariant <cit.>, used in <cit.>), two more points are required: the cast points c_.. The four points b_1, c_1, c_2, b_2, invariant under projective transformation, form the vertices of a polygon—the level line frame ℱ_l—and are mapped to the unit square vertices, ℱ_c (the canonical frame) <cit.>. The resulting projection is applied to ℒ_r and provides a canonical curve that can be used for shape comparison and matching.

A cast point c_1 (resp. c_2) is defined as the contact point of the tangent to ℒ_r that intersects the level line at b_1 (resp. b_2). There exist several such points for each bitangent point in the case of long bitangents. It thus becomes critical to ensure that a candidate frame ℱ_l forms a convex polygon so as to get an acceptable projection of ℒ_r to the canonical frame. In the case of short bitangents, the construction of ℱ_l is straightforward as only two cast points exist. As for long bitangents, a single portion of curve may be associated with several canonical curves, each of which depends on the frame configuration. As noted in <cit.>, it is preferable to pick those making a wide angle between the bitangent and the cast tangents, as well as those having the cast points as far from one another as possible: unbalanced frames may give distorted canonical curves. This holds for bitangents crossing the level line. It is also worth mentioning that this step drastically prunes the set of bitangents that can lead to satisfying frames. 













 §.§ Canonical curves



The goal is here to determine the 2D homography matrix such that 𝐱_𝐢=ρ𝖳𝐗_𝐢 <cit.>, where 𝐗_𝐢 = [X_iY_i1]^T is the i-th point in F_l (which no 3 are colinear) in homogeneous coordinates, 𝐱_𝐢 = [x_iy_i1]^T is the i-th vertex of the unit square defined by (0,0,1), (0,1,1), (1,1,1) and (1,0,1), 𝖳 is a [3× 3] matrix of the transformation parameters with 𝖳_33=1 and ρ is a non-zero scalar that gives by elimination 8 equations from four correspondences, linear in the parameters. The solution we are seeking is the unit singular vector corresponding to the smallest singular value of the matrix of the system.

A normalisation step, which consists of translating and scaling, is recommended for it forces the entries of the matrix of the system to have similar magnitude. Further details can be found in <cit.> (p.108). 





 §.§ Comparing polygonal curves


Contrary to <cit.>, who relied on rays extended from an origin (1/2,0) in ℱ_c and designed a feature vector made of all the distances from every intersection point with the canonical curve to the origin, we compare canonical curves by means of the Fréchet distance (Fig. <ref>). The rationale is that we also consider bitangents that cross level lines. This means that the canonical curves may cross the base of ℱ_c one or several times with more or less complex convolutions, making the use of rays impractical.

There are (at least) two common ways of defining the similarity between polygonal curves: the Hausdorff distance <cit.> and the Fréchet distance. The latter has the advantage that it takes into account the ordering of the points along the curves, thereby capturing curves structure better <cit.>. For the sake of speed, we used the discrete Fréchet distance (see Table 1 of <cit.>), which is an approximation of the continuous Fréchet distance: it only uses the curves vertices for measurements. From there, one can also define the reachable free space, which is the set of points for which the distance between two curves is lower than a distance parameter, δ and this allows tracking local similarity <cit.>. The Fréchet distance is the minimum δ that allows reaching the top right corner of the free space starting from (0,0).









§ MATCHING AND REGISTRATION


By cross-comparisons between histology and MRI, one obtains a measure of shape similarity (<ref>). Because each level line is associated with many canonical curves (one for each shape element), matches are found when the Fréchet distance is minimum and below a certain threshold. We can then use correspondences between ℱ_l in histology and MRI to compute an affine transformation (same principle as in <ref> with only 3 points—each providing two equations—and 𝗉=[T_11T_12…T_23001]^T). In order to minimise the global alignment error, the points must be well-arranged in images, i.e. the frames should be as wide as possible (hence the advantage of using long bitangents). When considering several level lines in both modalities, each canonical curve of each level line from histology returns at most one matching canonical curve for each level line in the MRI. False matches are filtered out using random sample consensus (RANSAC) <cit.> and a single global transformation is computed.








§ RESULTS AND DISCUSSION



We evaluated the method on 7 pieces of tissue, altogether covering 3 different subjects. For each subject, we had access to T2w, PSIR and PD MRI volumes (7 slices, 0.25×0.25×2mm^3). From these volumes we selected the slice that visually looked the most similar to one piece of histology. Histological images were a series of 11 consecutive 2μ m-thick sections, stained with 11 different dyes. At this point, it is worth noting that because the histological slab was about 25μm-thick—compared to a 2mm-thick slice from MRI—projective invariance not only allowed being robust to tissue distortions in the recognition process but was also required in order to tolerate morphological variations happening within that 2mm gap.

A ground truth arrangement similar to that of Fig. <ref> was available for direct assessment of success or failure of the alignment. It was made by a histopathologist at the time of the tissue preparation and essentially consisted of reporting the cassettes locations onto a slice of a medical image in order to keep track on which part of the sample the tissue piece was cut from. In the following, we call confusing (as opposed to meaningful <cit.> i.e., the tissue outer/inner boundaries) level lines, those not providing relevant information about the tissue shape. 


We ran two experiments (Fig. <ref>): (E1) consisted of using levels multiple of 16, 12, 8 and 4 in histology and MR images to investigate two questions: what is the impact of confusing level lines as well as their number, on the matching and the alignment? Can level lines be used as they are, without any form of prior knowledge about the tissue boundaries in images? Note that level lines were computed at quantised levels 0 to 255 by steps of 1. We expect that the sparser the set of level lines, the less informative about the actual tissue shape they can be (since information is lost when quantisation is coarse). This indeed translates in higher numbers of false than true matches when using between 1/16th and 1/8th of all available level lines (except for pieces 1 and 3 when using 1/8th, but this is hardly representative). When sufficient information comes in (1/4th), recognition becomes more successful: despite finding more false than true matches for piece 2, RANSAC was able to return the correct transformation—most of the false matches being isolated and spread across the MR image domain in that case. In contrary, RANSAC was unable to deal with false matches for pieces 5 and 6, those being related to ambiguities (shape elements were small and confusing).

The second experiment (E2) investigated the question: how robust is the matching/alignment when injecting confusing information into a subset of meaningful level lines? As such, we increased the number of neighbouring level lines from ±5 to ±20 around a meaningful one. In practice, meaningful level lines are those around structural layers (contrasted boundaries) of the tissue and we manually picked the corresponding levels. We can observe that the more localised around relevant information the level lines are, the higher the ratio true/false matches and the more trustful the set of correspondences fed to RANSAC. This is where redundancy is very valuable. However, the more levels one includes, the further one goes from meaningful information, and the more confusing it can get (see the increase in false matches). Due to the complexity of the information and the sinuosity of the shape, we believe that starting from a meaningful subset of level lines is an important consideration. 







Resulting alignments are shown in Fig. <ref> for 3 pieces, considering neighbourhoods of ±10 level lines. Overall, 5 pieces were matched correctly and two incorrectly. As for piece 6, no shape element was discriminative enough to be correctly matched with an MRI portion of level line without any ambiguity (Fig. <ref>c), as only relatively short bitangents could be extracted. As for piece 7, this is due to the fact that it is close to convex (and thus was not considered in the previous experiments). As a result, a few or no bitangents could be extracted from that histological image and no match was therefore available. 

The main requirements of the approach are twofold and relate to the length of the bitangents and the threshold on the Fréchet distance. As stressed out earlier, short bitangents convey little and ambiguous information about shape. This results in false matches especially because of the tolerance of the projective-invariant setting and the sinuosity of the MRI level lines. As a matter of fact, we constrained the approach to using long bitangents: in practice, we used those covering portions of a level line with more than 6 inflections. If a histological image happened to have informative portions with more than two inflections but less than 6—as it was the case for piece 5—then the longest bitangents were used (4 and 5 inflections in that case). An upper bound was also set (we chose 10 inflections) in order to speed up the matching process and avoid aberrant comparisons with bitangents covering the whole MR image; that range was applied to both MRI and histology. The rationale for considering such a range is also that it is not guaranteed that two level lines have the exact same number of inflections on corresponding portions across modalities, but their smoothness ensures those numbers are close. Long bitangents produce characteristic canonical curves (furthermore associated with wide frames) and allow for lower thresholds on the Fréchet distance while discarding false matches better.













§ CONCLUSION


This paper stands as a proof of concept that multimodal registration between a piece of tissue from histology and its whole in an MRI—which, to the best of our knowledge, remains to be addressed—is achievable as a result of shape recognition using portions of their level lines. Such a formulation allows for contrast, projective invariant representation of shape elements and partial matching regardless of the orientation of the piece of tissue on the glass slide (flips, rotations). We also introduced a computationally efficient refinement of bitangents using ellipses, from which a single bitangent was retained according to the local geometry of the level line. All this however, is to be related with the complexity of medical images; successful alignments require subsets of meaningful level lines along with characteristic shape elements. Those were obtained via the extension of Rothwell's framework to bitangents crossing the level lines and by preferring long bitangents.

Future works include: (i) the automatic extraction of meaningful level lines <cit.>; (ii) the use of shortcut Fréchet distance <cit.>, which bypasses large dissimilarities. This could improve robustness to tissue tears: a level line in histology may be globally close in terms of its shape to part of another in MRI but because of a tear that it follows, the distance between the associated canonical curves will be large.



§ ACKNOWLEDGMENTS

The authors would like to thank Prof. Olga Ciccarelli and her group (UCL Institute of Neurology, Queen Square MS Centre), for kindly providing the data.




This research was supported by the European Research Council (Starting Grant 677697, project BUNGEE-TOOLS), the University College London Leonard Wolfson Experimental Neurology Centre (PR/ylr/18575), the Alzheimer's Society UK (AS-PG-15-025), the EPSRC Centre for Doctoral Training in Medical Imaging (EP/L016478/1), the National Institute for Health Research University College London Hospitals Biomedical Research Centre and Wellcome/EPSRC (203145Z/16/Z, NS/A000050/1).























































































































































































































































































































































































ieee



