- Research Article
- Open Access
Constrained distance transforms for spatial atlas registration
- Bill Hill^{1}Email author and
- Richard A Baldock^{1}
https://doi.org/10.1186/s12859-015-0504-5
© Hill and Baldock; licensee BioMed Central. 2015
- Received: 13 October 2014
- Accepted: 20 February 2015
- Published: 18 March 2015
Abstract
Background
Spatial frameworks are used to capture organ or whole organism image data in biomedical research. The registration of large biomedical volumetric images is a complex and challenging task, but one that is required for spatially mapped biomedical atlas systems. In most biomedical applications the transforms required are non-rigid and may involve significant deformation relating to variation in pose, natural variation and mutation. Here we develop a new technique to establish such transformations for mapping data that cannot be achieved by existing approaches and that can be used interactively for expert editorial review.
Results
This paper presents the Constrained Distance Transform (CDT), a novel method for interactive image registration. The CDT uses radial basis function transforms with distances constrained to geodesics within the domains of the objects being registered. A geodesic distance algorithm is discussed and evaluated. Examples of registration using the CDT are presented.
Conclusion
The CDT method is shown to be capable of simultaneous registration and foreground segmentation even when very large deformations are required.
Keywords
- Constrained distance transform
- Non-rigid registration
- Atlas informatics
Background
The use of spatially mapped databases has become widespread within the biomedical research community [1]. Many of these databases, such as EMAGE [2] and the Allen Brain Atlas [3], are based on volumetric atlases or reference models with assay data mapped onto the atlas models through non–linear spatial transformations or warps. Compared to the warps used in medical imaging to register images in a longitudinal study or between patients, the warps required for these atlases are complex and challenging. With the significant challenges including variations in pose, mutant phenotypes, inter–species registration and the frequently non-corresponding image values due to gene expression or other spatial signals. Yet it is often in these most challenging of cases that the biological interest is greatest.
There may be only a weak relationship between the image values of an assay and the appropriate atlas image. This may be because of differences in imaging modalities, but more critically because the assay image is of a spatially distributed signal and the presence of that signal obscures or modifies the structural image across modalities. In such cases it is not possible to acquire a reference image closely similar to the atlas model. When combined with extreme variation in pose, these problems may result in an algorithm struggling to find points of correspondence between such atlas and assay images; an expert however can often find these points relatively quickly. In such cases the time spent by an expert may be significantly less than that spent correcting correspondences found automatically by an algorithm.
Radial basis function (RBF) transforms are frequently used for interactive image registration and perform well for small deformations. However, when the deformation gradients are large these methods may produce non-diffeomorphic, mirrored or extremely distorted mappings [4]. Methods based on elasticity with a uniform homogeneous material are also unable to register images where large deformations are required because of pose, giving rise to severe image distortion. In many cases the problematic large deformation gradients required for these methods are necessary only because the methods do not respect object boundaries. Current methods for registering images where large deformations are required include articulated and fluid models. Articulated models have been used for registering hand radiographs [5], in which skeletons composed of articulated rods are registered using landmarks at the ends of the rods and the displacements away from the rods are interpolated using weighted combinations of affine transforms. Fluid models based on solving viscoelastic systems have been developed but the fluid deformation model, like the articulated model, is often inappropriate [4] and the computation time for fluid models may be prohibitive.
For an interactive registration method to be useful it must be possible to compare the registered assay image with the model in a reasonable time. RBF transforms have a high computational cost because they typically rely on the evaluation of transcendental or other expensive functions, with the number of evaluations being proportional to the image volume and number of landmarks. Large numbers of landmarks may be required to force acceptable deformations, together with 3D image volumes, these may make the computation time unacceptable for interactive use. The computational cost of basis function evaluations can however be reduced by using mesh based methods, with evaluation of the basis functions only at the nodes off the mesh and a low cost interpolation within the mesh elements [6,7].
Method
Object representation
Images are frequently represented as simple rectangular or cuboid arrays of image values, but there are many situations in which it is useful to separate the representations of regions of space or domains from the image values that the domains cover. Such representations allow multiple objects with different domains to share a single set of values or value table; alternatively they may allow multiple objects with different value tables to share a single domain. In many cases it may be only the domains (and not the values within them) that are of interest and in these cases operations on the domains may be made more efficient through an appropriate representation such as interval coding [8]. This separation of domain and values is particularly appropriate for spatial atlases, where spatial queries can then be more efficiently handled. Using this representation; an image object O will always have a domain Ω and may have a value table V. A transformation (T) may then be expressed as a mapping from a source domain to a target domain: T:Ω _{ s }→Ω _{ t }.
Radial basis functions
where O _{0} is a zero matrix, O _{1} is a zero column vector, a is a column vector of the polynomial coefficients, λ is a column vector of the basis function coefficients, R is a symmetric matrix with the radial basis function values evaluated at the landmarks, X is a matrix containing the coordinates of the landmarks and D is a column vector of the displacements for the landmarks. The design equation may then be solved using a linear system solver such as singular value decomposition [9], although in practise it is often beneficial to rescale the parameters to reduce the condition number of the design matrix [6].
Radial basis functions and their form
Radial basis function | Form |
---|---|
Thin plate spline | ϕ(r)=r ^{2} logr |
Multiquadric | \( \phi (r) = (r^{2} + c^{2})^{\frac {1}{2}}\) |
Inverse multiquadric | \(\phi (r) = (r^{2} + c^{2})^{-\frac {1}{2}}\) |
Wendland’s functions | \(\phi (r) = \left \{ \begin {array}{r@{\quad :\quad }l} p(r) & 0 \leq r \leq 1 \\ 0 & r > 1 \end {array} \right. \) |
Geodesic distances
RBFs conventionally use Euclidean distances, however in the SCDT we wish to constrain the transformations using distances evaluated along paths that are constrained to the object’s domain. The minimum distance between two points within a convex domain is always the Euclidean distance, but when the domain is non-convex and the path between the two points is constrained to the domain then the Euclidean distance is the lower limit for the geodesic distance [16]. Reviews and evaluations of geodesic distance transform algorithms, are given in [17] and [18]. It is this geodesic distance that is used in the CDT. One of the first algorithms for computing geodesic distances was that of Piper and Granum [19], which like many later algorithms is based on region growing.
An early implementation of the CDT used a region growing algorithm based on morphological operators (similar to [19]) to compute geodesic distances, but as this was computationally prohibitive for interactive landmark placement a faster mesh based algorithm was developed. Because for the CDT geodesic distances are only required at the nodes of the mesh and at the landmark points, this led us to develop an algorithm for computing the geodesic distances at only the mesh nodes directly within the mesh using a fast marching algorithm.
Mesh based transforms
When a spatial transformation is applied to an object then it is common to apply a forward transformation to the object’s domain; image values within the domain are then interpolated in the space of the transformed domain using an inverse transformation. This two step transformation is in general necessary to avoid unassigned values in the transformed object. Many functions used for image transformation, such as radial basis functions, can not be analytically inverted. A mesh provides both a fast approximation to the RBF (or any other) displacement transform and allows the transformation to be inverted very simply and efficiently. Where the transform is costly to compute, then significant savings can be made by only computing the transform displacement at the nodes of a mesh. Wolberg [21] describes algorithms for 2D image re-sampling using regular quadrilateral meshes and scanline traversal. But there are significant advantages in using non-regular meshes: Biomedical images may contain little useful information in certain regions and these may have a courser mesh; the object being transformed may be non-rectangular and occupy a small part of it’s axis aligned bounding box or only the transformed domain may be required. All these factors allow significant performance increases to be achieved through the use of non–regular meshes, at the cost of increased complexity in the data structures and algorithms.
In this work the transformation of both object domains and values are accomplished using non-regular simplical conforming meshes. The use of conforming meshes simplifies image re-sampling in situations, such as limb articulation or lumen shrinkage, which would otherwise give rise to large mesh element deformations requiring complex re-meshing schemes. Objects are transformed by first defining a mesh covering their domain. This mesh then has displacements computed through the evaluation of the RBF nodal displacement with geodesic distances. The mesh is first used to forward transform the source object’s domain and then, if the source object has values associated with it, a new value table is created for the transformed domain and the new values interpolated using a mesh sweep line algorithm.
Mesh generation
Constrained distance transforms
We have developed a registration method which we have called the Constrained Distance Transform (CDT). In this method displacements are computed within a mesh conforming to the target (or a source) domain by evaluating RBFs at the mesh nodes and using distances evaluated along paths constrained to geodesics within the mesh. Displacements within mesh elements are computed by interpolating nodal values. Using a CDT, connectivity and distance are defined by the domain conforming mesh and the problems associated with large deformations (such as to correct for pose) are significantly reduced. Regions close in Euclidean space may be easily pulled apart without distortion provided that they are distant with respect to geodesics in the conforming mesh. The CDT also allows multiple geodesically distant regions of a object to be fused.
Implementation
The algorithms and data structures described in this paper have been implemented in C within the Woolz open source image processing system. Woolz also includes a Java Native Interface (JNI) binding making it possible to use Woolz from Java as well as native applications (https://github.com/ma-tech).
Results and discussion
Applicability to atlas registration
Because CDTs are invertible the mesh may be defined either on the source or target. This is important because for atlas systems, in which the target is an atlas model and the source an assay object, the number of meshes required is reduced to the number of atlas models and the additional cost of building these meshes may not be significant when compared to the total cost of building the atlas model. By defining the mesh using the atlas model, assay images are segmented from their background through their registration to the pre-segmented assay model, since values that fall outside of the mesh (and consequently the atlas model’s domain) are not mapped. We do not claim this segmentation to be novel, however it is an intrinsic feature of our method. In most cases this is desirable, but if ignored, it could lead to incorrect conclusions being drawn from the absence of mapped source regions.
Distance evaluation errors
Errors and execution times for 2D distance evaluations
Algorithm | mesh | mesh_nlsi | c4 | c8 | oct |
---|---|---|---|---|---|
% error | -0.0005 ± | -0.007 ± | -0.07 ± | -0.07 ± | -0.03 ± |
0.0007 | 0.009 | 0.02 | 0.02 | 0.01 | |
Time (ms) | 0.9 | 0.9 | 23 | 16 | 17 |
Errors and execution times for 3D distance evaluations
Algorithm | mesh | mesh_nlsi | c6 | c26 | oct |
---|---|---|---|---|---|
% error | -0.004 ± | -0.013 ± | -0.07 ± | -0.11 ± | -0.045 ± |
0.002 | 0.006 | 0.02 | 0.02 | 0.008 | |
Time (ms) | 491 | 492 | 3832 | 2991 | 2878 |
Displacement errors
Using a mesh in which the displacements are computed at it’s nodes and approximated within it’s elements results in some displacement error. These errors were evaluated for 3D conforming meshes with varying numbers of nodes but corresponding to the same domain.
Mesh resolution
3D Displacement errors
N | \({\frac {1}{\sqrt [3]{N}}}\) | l _{ RBF } | l _{ MESH } | ε |
---|---|---|---|---|
7963 | 0.050 | 65.7 ±78.5 | 65.4 ±78.1 | 0.0134 ±0.0268 |
13899 | 0.042 | 61.1 ±75.0 | 60.9 ±74.5 | 0.0095 ±0.0235 |
40255 | 0.029 | 62.2 ±75.4 | 62.0 ±74.5 | 0.0061 ±0.0278 |
117537 | 0.020 | 58.9 ±75.0 | 58.8 ±75.0 | 0.0046 ±0.0325 |
1294844 | 0.009 | 67.1 ±84.2 | 67.0 ±82.9 | 0.0023 ±0.0259 |
Mesh defined locality
Within a CDT displacements are determined by the choice of RBF, landmarks and the constrained distance from the landmarks. A 2D assay image was selected from the e-Mouse Atlas Gene Expression database (EMAGE) in which the head and tail are close in Euclidean space (http://www.emouseatlas.org/emap/home.html and http://www.emouseatlas.org/emage/home.php). The image was segmented from its background and a conforming mesh was created for the segmented domain. Landmark points were then evenly distributed around the domain, with all the landmarks having zero displacements except for those at the tip of the tail which had a displacement set away from the head. Standard unconstrained RBF and CDT warps were applied using a MQ RBF in both warps. Figure 1 shows this image together with the landmark points and the displacement applied. The figure also shows the results of applying the CDT and an unconstrained RBF transform to the image. In the CDT warped image the head does not show any significant deformation despite the large deformation experienced by the tail, unlike the unconstrained RBF transform which resulted in large deformations to both the head and the tail. In both cases the same landmarks and displacements were used. The resulting warp shows how locality in CDTs is defined by the mesh.
2D Atlas registration
A 2D assay image was selected from the EMAGE database along with the corresponding 2D atlas image. The 2D atlas image is a projection of the 3D atlas model in a standard pose with the tail articulated away from the body to allow independent 2D spatial mapping on both parts. Using conventional mapping techniques this mapping would need to be performed separately for the body and tail, but since connectivity is defined by the mesh for CDT warps, they can be transformed simultaneously when using a mesh defined on the target. A mesh was constructed for the pre–segmented atlas (target) image and 25 landmark pairs were defined between the assay (source) and target images. Together with an inverse multiquadric RBF and δ value of 0.05 these defined the transformation from the assay (source) image to the atlas (target) image. Using a 2D CDT the assay image was registered to the atlas image and segmented in a single operation with the assay tail correctly mapped to the articulated tail. Image values in regions of the tail appear twice in the warped image illustrating that the CDT may produce one-to-many mappings. The source, target and transformed images can be seen in Figure 4.
3D Atlas registration
3D Displacement and warp computation times
N | t _{ disp } (ms) | t _{ warp } (ms) |
---|---|---|
7963 | 208 | 1764 |
13899 | 358 | 2027 |
40255 | 1336 | 2630 |
117537 | 4636 | 3678 |
1294844 | 65518 | 10749 |
Conclusions
In this paper we have described a novel image registration algorithm, the CDT and demonstrated it’s applicability for atlas registration though it’s ability to register both 2D and 3D images to an atlas whilst performing simultaneous foreground segmentation in the presence of large deformations such as those required to correct for pose. We are unaware of any other published method which is capable of performing such registration.
A novel initialisation for computing distance within meshes using fast marching has been described which improves accuracy yet has minimal impact on the computation time.
The cost and difficulty of defining large number of landmarks required for very accurate alignment may be prohibitive and a two part approach may be preferred in which large deformations are computed using CDT and remaining small deformations through other algorithms.
Declarations
Acknowledgements
We acknowledge support from the Medical Research Council, UK and NIH support under grant #1R01MH070370-01A2.
Authors’ Affiliations
References
- de Boer BA, Ruijter JM, Voorbraak FPJM, Moorman AFM. More than a decade of developmental gene expression atlases: where are we now?Nucleic Acids Res. 2009; 37(22):7349–59.View ArticlePubMedPubMed CentralGoogle Scholar
- Baldock R, Bard J, Burger A, Burton N, Christiansen J, Feng G, et al. EMAP and EMAGE: A framework for understanding spatially organised data. Neuroinformatics. 2003; 1:309–25.View ArticlePubMedGoogle Scholar
- Lein ES, Hawrylycz MJ, Ao N, Ayres M, Bensinger A, Bernard A, et al. Genome-wide atlas of gene expression in the adult mouse brain. Nature. 2007; 445(7124):168–76.View ArticlePubMedGoogle Scholar
- Modersitzki J. Numerical methods for image registration. New York: Oxford University Press; 2004.Google Scholar
- Martín-Fernández MA, Muñoz-Moreno E, Martín-Fernández M, Alberola-López C. Articulated registration: Elastic registration based on a wire-model In: Fitzpatrick JM, Reinhardt JM, editors. Medical Imaging 2005: Image Processing. San Diego, USA: SPIE Press: 2005. p. 182–91.Google Scholar
- Baldock RA, Hill B. 9. Image warping and spatial data mapping In: Baldock R, Graham J, editors. Oxford, UK: Oxford Uiversity Press: 2000.Google Scholar
- Qu L, Peng H. Littlequickwarp: An ultrafast image warping tool. Methods. 2015; 73(0):38–42.View ArticlePubMedGoogle Scholar
- Piper J, Rutovitz D. Data structures for image processing in a c language and unix environment. Pattern Recognition Lett. 1985; 3:119–129.View ArticleGoogle Scholar
- Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipies in C: The art of scientific computing. Cambridge, England: Cambridge University Press; 1992.Google Scholar
- Bookstein FL. Thin-plate splines and the atlas problem for biomedical images. In: Lecture Notes in Computer Science 511: Information Processing in Medical Imaging. Heidelberg, Germany: Springer: 1991. p. 326–42.Google Scholar
- Fogel DN. Image registration using multiquadric functions, the finite element method, bivariate mapping polynomials and the thin plate spline. Technical Report 96, National Center for Geographic Information and Analysis, USA, 1996.Google Scholar
- Buhmann MD. Radial basis functions. New York, NY, USA: Cambridge University Press; 2003.View ArticleGoogle Scholar
- Wendland H. Piecewise polynomial, positive definite and compactly supported radial basis functions of minimal degree. Adv Comput Math. 1995; 4:389–96.View ArticleGoogle Scholar
- Jacobs MA, Windham JP, Soltanian-Zadeh H, Peck DJ, Knight RA. Registration and warping of magnetic resonance images of histological sections. Med Phys. 1999; 26(8):1568–78.View ArticlePubMedGoogle Scholar
- Franke R. A critical comparison of some methods for interpolation of scattered data. Technical Report NPS-53-79-003. Monterey, California, USA: Naval Postgraduate School; 1979.Google Scholar
- Lantuejoul C, Maisonneuve F. Geodesic methods in quantitative image analysis. Pattern Recognit. 1984; 17(2):177–87.View ArticleGoogle Scholar
- Cuisenaire O. Distance transformations: Fast algorithms and applications to medical image processing. Belgium: PhD thesis, Laboratoire de Telecommunications et Teledetection, Université catholique de Louvain; 1999.Google Scholar
- Cárdenes R, Alberola-López C, Ruiz-Alzola J. Fast and accurate geodesic distance transform by ordered propagation. Image Vis. Comput. 2010; 28(3):307–16.View ArticleGoogle Scholar
- Piper J, Granum E. Computing distance transformations in convex and non-convex domains. Pattern Recognit. 1987; 20(6):599–615.View ArticleGoogle Scholar
- Qian J, Zhang YT, Zhao HK. Fast sweeping methods for eikonal equations on triangulated meshes. SIAM J Numeric Anal. 2007; 45:83–107.View ArticleGoogle Scholar
- Wolberg G. Digital image warping. Los Alamitos, California: IEE Computer Society Press; 1992.Google Scholar
- Zhang Y, Bajaj C, Sohn B-S. Adaptive and quality 3d meshing from imaging data. In: SM ’03: Proceedings of the Eighth ACM Symposium on Solid Modeling and Applications. New York, NY, USA: ACM Press: 2003. p. 286–91.Google Scholar
- Schöberl J. Netgen - an advancing front 2D/3D-mesh generator based on abstract rules. Comput Visual Sci. 1997; 1:41–52.View ArticleGoogle Scholar
Copyright
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.