INTERPOLATION ARTIFACTS IN BIOMEDICAL IMAGE REGISTRATIONG.K. Rohde1, D.M. Healy, Jr.2, C.A. Berenstein2, A. Aldroubi31Department of Biomedical Engineering, Carnegie Mellon University, Pittsburgh, PA. 152132Department of Mathematics, University of Maryland, College Park, MD, USA. 20742.3Mathematics Department, Vanderbilt University, Nashville, TN, USA. 37235.ABSTRACTWe identify a cause for interpolation artifacts in objectivefunctions by observing that the energy of interpolated andtranslated 2sequences oscillates with respect to the transla-tion parameter. Using the B-spline interpolation framework,we show that such oscillations in the energy of the signals af-fect the sum of squared differences, cross correlation, and mu-tual information objective functions. We describe several ap-proaches that can be used to avoid interpolation artifacts (suchas higher degree interpolation, as well as stochastic sampling)and explain in detail why these are effective in eliminating theartifacts.Index Terms— Image registration, interpolation, artifacts,local optima1. INTRODUCTIONImage interpolation methods have been known to introduceartifacts (local optima) in objective functions used for auto-matic intensity-based image registration. To date most con-sider such artifacts to be a feature specific to the mutual in-formation similarity measure [1, 2, 3]. A common strategyfor avoiding such artifacts has been through the introductionof a randomizing sampling operation (see [2, 4, 3] for exam-ples), though a precise cause for the oscillatory artifacts whenregular sampling is used has not been mathematically demon-strated. In [5, 6] we showed that the covariance propertiesof an image undergoing spatial transformations can cause ar-tifactual oscillations, and in some cases local optima, in thesum of squared differences (SSD), cross correlation (CC), andmutual information (MI) objective functions. Our arguments,however, were entirely stochastic in that only the effects ofadditive noise were considered.Here we provide a more general and concise explanationfor oscillation artifacts in the objective functions used aboveand show that even interpolation on signals or images withno noise can cause oscillation artifacts. We explain three ap-proaches for avoiding them based on low pass filtering theimage, higher degree interpolation, and Monte Carlo integra-tion. We begin by showing that the energy of interpolated datasequences is present in all objective functions named above.Using the generic framework of B-spline image interpolationwe then show that the energy of interpolated, translated, andsampled data sequences oscillates with respect to the transla-tion parameter. We show that Monte Carlo integration meth-ods can be used to compute objective functions without ar-tifactual oscillations. Finally results and conclusions are of-fered.2. THEORYLet 2denote the space of square summable infinite dimen-sional real valued sequences. The inner product between twosequences a, b in 2is defined by a, b 2=k∈Za(k)b(k)while b22= b, b 2. A convolution between two 2se-quences is denoted a ∗ b and thus b can be though of as a dis-crete operator characterized by its transfer function B(z)=k∈Zb(k)z−k.IfB has no zeros on the unit circle, then theinverse operator (b)−1exists and is uniquely defined by:(b)−1←→ 1/B.A sequence s ∈ 2can be translated by an arbitrary amountθ ∈ R by first fitting an interpolation (or approximation)model to s and then sampling the model at a translated set ofcoordinates normally organized on a regular grid. Let β0referto the centered normalized rectangle. Then the B-spline func-tion of degree n at value x is given by βn(x)=βn−1∗β0(x).Let˜s(x)= k∈Zs(k)ηn(x − k). (1)where ηn(x)=k∈Z(bn)−1(k)βn(x − k) and bn(k)=βn(x)|x=k. The translated data sequence can then be rep-resented by˜s(k + θ)=ηnθ∗ s(k) (2)where ηnθ(k)=ηn(k + θ). These operations can be extendedto dimensions 2 and over via tensor products of η.2.1. Objective functionsThe goal in image registration is to compute a function fθ:Rd→ Rd, where θ are parameters that determine f, such6481424406722/07/$20.00 ©2007 IEEE ISBI 2007that two signals (images) are aligned. This procedure can beautomated as an optimization problem:f∗θ=argmaxfθQ(s, t, fθ),where Q(···) is the objective function chosen for a specificproblem. An intuitive figure of merit for quantifying how welltwo sequences s and t align is the squared 2norm of its dif-ference (SSD). Let ˜sθ= {˜s(fθ(k))}k∈Z. The SSD objectivefunction, to be minimized in this case, is:Q(θ)=˜sθ− t22= t2− 2 ˜sθ,t 2+ ˜sθ22.Another objective function often used is the cross correlationbetween the signals:Q(θ)=˜sθ,t 2t2˜sθ2.Finally, the registration between two images can also becomputed by maximizing the MI between them. Originallyintroduced as a similarity measure for registration problemsin [7, 8], the MI between two signals ˜sθand t quantifies theamount of ‘information’ due to ˜sθis present in t. Let SθandT represent random variables associated with the discrete sig-nals {˜s(fθ(k))}k∈Zand {t(k)}k∈Z, respectively, with prob-ability density functions (pdf) pSθ(μ), pT(υ), and joint pdfpSθ,T(μ, υ). The mutual information is defined asQ(θ)=pSθ,T(μ, υ)logpSθ,T(μ, υ)pSθ(μ)pT(υ)dμdυUsing the Gram-Charlier series expansion for the pdfs above,the MI between two random variables can also be written as[9]:MI(θ)=IG(θ)+R(θ)whereIG(T,Sθ)=−12log1 − ρ2(θ),ρ(θ)=ξ(θ)/σTσSθξ(θ) stands for the covariance between random variables Sθand T , and σ2Tand σ2Sθrepresent their respective variances.R(θ) is a residual term. Assuming stationarity and ergodic-ity, these quantities can be estimated from a set of samples˜s(fθ(k)) and t(k), k ∈ Ω:ξ(θ)=1N k∈Ω˜s(fθ(k))t(k),σ2T=1N k∈Ω|t(k)|2, andσ2Sθ=1N k∈Ω|˜s(fθ(k))|2,Fig. 1. Frequency response for cardinal B-spline translationfilters for different translation values (a) and degrees (b).assuming the signals are zero mean, where Ω represents a setof predefined coordinates and N = |Ω|. Clearly, σ2Sθis re-lated to ˜sθ22via:σ2Sθ=1N˜sθ22−1N k/∈Ω|˜s(fθ(k))|2.2.2. Oscillation artifactsNote that the term ˜sθ2is present in all three objectivefunctions