Registration600.445 Computer-Integrated SurgeryRussell H. TaylorWhat needs registering?• Preoperative Data– 2D & 3D medical images– Models– Preoperative positions• Intraoperative Data– 2D & 3D medical images– Models– Intraoperative positioning information• The PatientA typical registration problemIntraoperativeRealityPreoperativeModelFramework• Definition of coordinate system relations• Segmentation of reference features• Definition of disparity function between features• Optimization of disparity functionDefinitionsOverall Goal: Given two coordinate systems, and coordinatesassociated with homologous features in the two coordinate systems, the general goal is to determine a transformation function T that transforms one set of coordinates into the other:Ref RefAB & xxAB & xTxAB= ()Definitions• Rigid Transformation: Essentially, our old friends 2D & 3D coordinate transformations:T(x) = R•x + pThe key assumption is that deformations may be neglected.• Elastic Transformation: Cases where must take deformations into account. Many different flavors, depending on what is being deformedUses of Rigid Transformations• Register (approximately) multiple image data sets• Transfer coordinates from preoperative data to reality (especially in orthopaedics & neurosurgery)• Initialize non-rigid transformationsUses of Elastic Transformations• Register different patients to common data base (e.g., for statistical analysis)• Overlay atlas information onto patient data• Study time-varying deformations• Assist segmentationTypical Features• Point fiducials• Point anatomical landmarks• Ridge curves• Contours• Surfaces• Line fiducialsDistance FunctionsGiven two (possibly distributed) features Fi andFj, need to define a distance metric distance (Fi,Fj) between them. Some choices include:– Minimum distance between points– Maximum of minimum distances– Area between line features– Volume between surface features– Area between point and line–etc.Disparity Functions Between Feature SetsOptimization• Global vs Local• Numerical vs Direct Solution• Local MinimaSampled 3D data to surface modelsA typical registration problemIntraoperativeRealityPreoperativeModelWhat the computer knowsFind homologous points & pull!Find homologous points & pull!Find homologous points & pull!Iterate this until convergeFind new point pairs every iterationKey challenge is finding point pairs efficiently.Head in Hat AlgortithmHead in Hat AlgorithmHead-in-hat algorithm: step 0Head-in-hat algorithm: step1Head-in-hat algorithm: step1Head-in-hat algorithm: step 2Head-in-hat algorithm: step 2Head-in-hat algorithm: step 3Head in Hat Algorithm• Strengths– Moderately straightforward to implement– Slow step is intersecting rays with surface model– Works reasonably well for original purpose (registration of skin of head) if have adequate initial guess• Weaknesses– Local minima– Assumptions behind use of centroid– Requires good initial guess and close matches during convergenceMinimizing Rigid Registration ErrorsTypically, given a set of points { in one coordinate systemand another set of points { in a second coordinate systemGoal is to find [ imizeswhere This is tricky, because of iiabRpeeeRapbR}}, ] that min().η=•=•+−∑iiiiiiMinimizing Rigid Registration Errors11211()(,)NNiiiiii iiiiiNNFrame=====− =−⋅−=− ⋅=∑∑∑aa bbaaa bbbRRa bppbRaFRprrrr%%%%rrrStep 1: Compute Step 2: Find that minimizesStep 3: Find Step 4: Desired transformation isSolving for R: iteration method(){}()0121,,, argmin,()ii i iikikiiiikk−+=−=∆∆−=∆∑∑ab R Ra bRRbRbRRa bRRR%%%%LL(%(%Given , want to find Step 0: Make an initial guess Step 1: Given compute Step 2: Compute that minimizesStep 3: Set Step 4: Iterate Steps 1-3 until residual error is sufficiently small (or other termination condition)Iterative method: Solving for ∆R22())..()min( ).ii ii iiiskew∆α∆+α∆•≈+α×∆• − ≈ − +α×α∑∑RRIRv v vvRa b a b a((%%%Approximate as ( I.e., for any vector Then, our least squares problem becomesmin This is linear least squares problem in Then com().∆αRputeDirect Techniques to solve for R,, ,, ,,,, ,, ,,,, ,, ,,Step 1: Compute Step 2: Compute the SVD of = Step 3: Step 4: Verify ( ) 1. If not, then algorithm ix ix ix iy ix iziy ix iy iy iy iziiz ix iz iy iz izab ab abab ab abab ab abDet===∑ttHHUSVRVUR may fail.• Method due to K. Arun, et. al., IEEE PAMI, Vol 9, no 5, pp 698-700, Sept 1987 • Failure is rare, and mostly fixable. The paper has details.Quarternion Technique to solve for R• B.K.P. Horn, “Closed form solution of absolute orientation using unit quaternions”, J. Opt. Soc. America, A vol. 4, no. 4, pp 629-642, Apr. 1987.• Method described as reported in Besl and McKay, “A method for registration of 3D shapes”, IEEE Trans. on Pattern Analysis and Machine Intelligence, vol. 14, no. 2, February 1992. • Solves a 4x4 eigenvalue problem to find a unit quaternion corresponding to the rotation• This quaternion may be converted in closed form to get a more conventional rotation matrixDigression: quaternions[][]012301 2 31,,,,qqqqssqq q q==+==+ + +λ+qqvvqijkqrrrrrInvented by Hamilton as a way to express the ratioof vectors. Can be thought of as 4 elements: scalar & vector: Properties:Linearity:[][][][]21212*12 12 1212 21 1 2*22 2220123,,,0,ssssss s ssqqqqµ=λ+µλ+µ=−= −=− ++×==+= +++qvvqvvqq vv v v vvqp q p qqvvrrrrrrrr rrro&rrooorr&Conjugate:Product: Transform vector: Norm:Digression continued: unit quaternions(, cos ,sin2(( , cos , sin 02RotRotθθθθ⇔θθθ=nnnnp nrrrrr r&oWe can associate a rotation by angle about an axis with the unit quaternion:) 2Exercise: Demonstrate this relationship. I.e., show) 2[],cos,sin2θθ−pnrro 2Rotation matrix from unit quaternion[]012322330123 1203 130222 3312 03 0 1 2 3 23 01223313 02 23 01 0 1 2 3,,, 12( ) 2( )() 2( ) 2( )2( ) 2( )qqqqqqqq qqqq qqqqqq qq q q q q qq qqqq qq qq qq q q q q==+−− − +=+ −+− −−+−−+qqRq;Unit quaternion from rotation matrix01230123000001231;1() ;1;1max{ } max{ } max{ } max{ }2444xx yx zxxx yy zz xx yy zzxy yy zyxx yy zz xx yy zzxz yz zzkkkkyz zy xy yxzx xzrrrarrrarrrrrrarrrarrrrrraaaaaaaarr rrarrqqqqqqqq=+ + + =+ −