\documentclass{InsightArticle} \usepackage{bm} \usepackage[dvips]{graphicx} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % hyperref should be the last package to be loaded. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \usepackage[dvips, bookmarks, bookmarksopen, backref, colorlinks,linkcolor={blue},citecolor={blue},urlcolor={blue}, ]{hyperref} \title{Helper classes for the \doxygen{BSplineDeformableTransform}} % % NOTE: This is the last number of the "handle" URL that % The Insight Journal assigns to your paper as part of the % submission process. Please replace the number "1338" with % the actual handle number that you get assigned. % \newcommand{\IJhandlerIDnumber}{????} % Increment the release number whenever significant changes are made. % The author and/or editor can define 'significant' however they like. \release{0.00} % At minimum, give your name and an email address. You can include a % snail-mail address if you like. \author{Marius Staring$^{1}$ and Stefan Klein$^{2}$} \authoraddress{$^{1}$Leiden University Medical Center, The Netherlands\\% $^{2}$Biomedical Imaging Group Rotterdam, Erasmus MC, Rotterdam, The Netherlands.} \begin{document} % % Add hyperlink to the web location and license of the paper. % The argument of this command is the handler identifier given % by the Insight Journal to this paper. % \IJhandlefooter{\IJhandlerIDnumber} \ifpdf \else % % Commands for including Graphics when using latex % \DeclareGraphicsExtensions{.eps,.jpg,.gif,.tiff,.bmp,.png} \DeclareGraphicsRule{.jpg}{eps}{.jpg.bb}{`convert #1 eps:-} \DeclareGraphicsRule{.gif}{eps}{.gif.bb}{`convert #1 eps:-} \DeclareGraphicsRule{.tiff}{eps}{.tiff.bb}{`convert #1 eps:-} \DeclareGraphicsRule{.bmp}{eps}{.bmp.bb}{`convert #1 eps:-} \DeclareGraphicsRule{.png}{eps}{.png.bb}{`convert #1 eps:-} \fi \maketitle \ifhtml \chapter*{Front Matter\label{front}} \fi % The abstract should be a paragraph or two long, and describe the % scope of the document. \begin{abstract} \noindent This document describes two new classes that facilitate the use of the \doxygen{BSplineDeformableTransform} in multiresolution registration algorithms. The first class, \doxygen{GridScheduleComputer}, defines the B-spline grid, based on a user-specified grid spacing and an input image (the fixed image in registration). A different grid is determined for each resolution level of the registration, given a multiresolution schedule similar to that of the image pyramids. The second class, \doxygen{UpsampleBSplineParametersFilter}, resamples (usually upsamples) the B-spline coefficients to a new control point grid having a different spacing/origin/size. This can be used to transfer the result from one resolution level to that of the next level. Summarising, suitable grid definitions follow from the first class, and are used as input for the second class, which performs the actual resampling. The classes are implemented using the Insight Toolkit \url{www.itk.org}. This paper is accompanied with the source code, including a test program. \end{abstract} \IJhandlenote{\IJhandlerIDnumber} \tableofcontents \section{Introduction} Nonrigid image registration often requires a multiresolution strategy not only for the images, but also for the transformation. In case of a transformation modelled by B-splines, it is for example common to start with a coarse B-spline grid, which takes care of the more global deformations. The grid is gradually refined, to enable matching of finer structures. The B-spline grid is defined by the spacing between the control points and the position of the `first' control point (the grid origin). Some care has to be taken to place a band of control points outside the fixed image domain, such that the transformation is valid for all voxels. The grid definition is therefore also related to the support region of the B-spline basis functions. Since the grid needs to be defined for every resolution level, the functionality is needed regularly. A dedicated class simplifies the procedure. The new class \doxygen{GridScheduleComputer} takes care of computing the required information for all resolution levels. Another aspect of multiresolution registration is the initialisation of the transformation in one resolution level, based on the result of the previous resolution level. This requires evaluation of the transformation at the control point locations of the new (finer-scale) grid, followed by a B-spline decomposition. The new class \doxygen{UpsampleBSplineParametersFilter} implements this procedure. Included with the source code is a program that tests all functionality of the two classes. \section{The GridScheduleComputer class} Figure \ref{fig:gridplacement} illustrates the grid placement strategy implemented by the \doxygen{GridScheduleComputer} class. We chose to place the B-spline grid such that the image is in the centre. Note the band of control points outside the image. More details on the exact implementation can be found in the \code{ComputeBSplineGrid()}-function of the class. \begin{figure} \centering \includegraphics[width=0.3\textwidth, angle=90]{grid.eps} \caption{Illustration of the control point grid placement. The square represents the image. The dots are the B-spline control points. The example is for a 3rd order B-spline basis function. Note that we chose to center the grid over the image.} \label{fig:gridplacement} \end{figure} In many registration problems, the nonrigid B-spline registration is preceded by a rigid or affine registration. The transformations could be simply summed, in which case there is no problem. In case the transformations are combined by \emph{composition}, the grid placement algorithm needs to be modified. Imagine, for example, that the result of rigid registration is a global translation of 10\,mm in the $x$-direction. Consider a point close to the rightmost boundary of the fixed image. The voxel is first translated 10\,mm to the right and thus possibly falls outside the valid B-spline grid region. Therefore, the B-spline grid should have been placed 10\,mm to the right. To facilitate this, the \doxygen{GridScheduleComputer} has an optional input: \code{SetInitialTransform()}. The \doxygen{GridScheduleComputer} requires four user inputs: the \code{BSplineOrder}, the \code{FinalGridSpacing}, the \code{GridSpacingSchedule}, and the domain of the fixed image (given by \code{ImageOrigin}, \code{ImageSpacing}, and \code{ImageRegion}). The \code{BSplineOrder} is used to compute the width of the band of control points outside the image. The \code{FinalGridSpacing} denotes the control point grid spacing at the finest resolution level. The \code{GridSpacingSchedule} defines the downsampling scheme used to determine the grid at lower resolution levels. The schedule has a similar meaning as in the \doxygen{MultiResolutionPyramidImageFilter}. The schedule can be supplied in two ways. The first way has the most freedom: the user supplies a vector of a vector of downsampling factors. For every resolution and for every dimension the downsampling factor can thus be specified. The second, simplified way supplies the number of resolution levels and a single downsample factor. For example, a factor 2.0 will halve the grid spacing every resolution, for each dimension. Both schedules end with the final grid spacing. An example of usage is given below: \small \begin{verbatim} typedef GridScheduleComputer GridScheduleComputerType; GridScheduleComputerType::Pointer gridScheduleComputer = GridScheduleComputerType::New(); gridScheduleComputer->SetImageOrigin( fixedImage()->GetOrigin() ); gridScheduleComputer->SetImageSpacing( fixedImage()->GetSpacing() ); gridScheduleComputer->SetImageRegion( fixedImage()->GetLargestPossibleRegion() ); gridScheduleComputer->SetInitialTransform( initialTransform ); // optional gridScheduleComputer->SetFinalGridSpacing( finalGridSpacing ); gridScheduleComputer->SetSchedule( gridSchedule ); /** This function does the work. */ gridScheduleComputer->ComputeBSplineGrid(); /** Obtain the grid definition for a certain resolution level. */ gridScheduleComputer->GetBSplineGrid( resolutionLevel, gridRegion, gridSpacing, gridOrigin ); \end{verbatim} \normalsize \section{The UpsampleBSplineParametersFilter class} This filter is based on the code in \code{DeformableRegistration6.cxx}, found in the ITK repository (\code{/Examples/Registration}). The filter resamples a B-spline coefficient image on a new grid. Usually, it will be used to \emph{upsample} the grid, from one resolution level to the next. The filter takes as input the old (low-resolution) B-spline grid (origin, spacing, region) and the corresponding coefficients (obtained by \code{bSplineTransform->GetParameters()}). The coefficients are resampled on a new, user-specified ``required grid''. Computation is done by resampling the old grid, using the \doxygen{BSplineResampleImageFunction}, followed by B-spline decomposition using the \doxygen{BSplineDecompositionImageFilter}. More details can be found in the \code{UpsampleParameters()}-function of the class. An example of usage is given below: \small \begin{verbatim} typedef UpsampleBSplineParametersFilter< ParametersType, CoefficientImageType > GridUpsamplerType; GridUpsamplerType::Pointer gridUpsampler = GridUpsamplerType::New(); OriginType currentGridOrigin = bSplineTransform->GetGridOrigin(); SpacingType currentGridSpacing = bSplineTransform->GetGridSpacing(); RegionType currentGridRegion = bSplineTransform->GetGridRegion(); gridUpsampler->SetCurrentGridOrigin( currentGridOrigin ); gridUpsampler->SetCurrentGridSpacing( currentGridSpacing ); gridUpsampler->SetCurrentGridRegion( currentGridRegion ); /** Use the previously described itk::GridScheduleComputer to * get the grid domain definition of the next resolution level. */ OriginType requiredGridOrigin; SpacingType requiredGridSpacing; RegionType requiredGridRegion; gridScheduleComputer->GetBSplineGrid( newResolutionLevel, requiredGridRegion, requiredGridSpacing, requiredGridOrigin ); gridUpsampler->SetRequiredGridOrigin( requiredGridOrigin ); gridUpsampler->SetRequiredGridSpacing( requiredGridSpacing ); gridUpsampler->SetRequiredGridRegion( requiredGridRegion ); /** Get the parameters resulting from the last resolution level. * It is assumed in this example that the * itk::MultiResolutionImageRegistrationMethod is used. */ ParametersType latestParameters = registration->GetLastTransformParameters(); ParametersType upsampledParameters; gridUpsampler->UpsampleParameters( latestParameters, upsampledParameters ); bSplineTransform->SetGridOrigin( requiredGridOrigin ); bSplineTransform->SetGridSpacing( requiredGridSpacing ); bSplineTransform->SetGridRegion( requiredGridRegion ); registration->SetInitialTransformParametersOfNextLevel( upsampledParameters ); \end{verbatim} \normalsize \section{Conclusion} Two classes were introduced that facilitate multiresolution image registration, using the B-spline transformation. The classes implement commonly used functionality, thereby simplifying integration in the registration framework. Future work includes integration with the \doxygen{BSplineDeformableTransformInitializer}, see \url{http://hdl.handle.net/1926/1338}. That class requires less user interaction, since it directly acts on the transform object. It also takes into account the direction cosines, which we have neglected. An advantage of our class is that it centers the control point grid on the image domain. \end{document}