Programming for Image ProcessingAnalysis and Visualization using The Visualization Toolkit

1
Programming for Image Processing/Analysis and
Visualization using The Visualization Toolkit
Week 7 Case Study 2 Intensity--Based Linear
Image Registration
http//noodle.med.yale.edu/seminar/seminar.html

• Xenios Papademetris
• xenios.papademetris_at_yale.edu
• TAC (ex CAB) N119 5-6148

2
Schedule Part 2

• Review of Part 1 and Course Overview
• C Pointers/Classes, Object Oriented Programming
• VTK an object library
• Adding new VTK Commands/Cmake
• Image-to-image filters
• Case Study I -- Iterative Closest Point surface
  matching
• Case Study II Intensity-Based Linear Image
  Registration

3
Allocating/De-allocating Arrays of Pointers

• As before allocation is performed by the new
  command e.g.
• int anew int10 // Allocate an array of 10
  integers
• a is a pointer which contains the memory location
  of the first element of the array. To access the
  values stored at the memory locations use
• a0 1 a32 etc.
• When a is no longer needed the memory can be
  released using the delete command e.g.
• delete a
• Do not use the delete operator to delete arrays,
  it causes lots of problems use delete !!!

4
Talk Outline

• A review of the VTK Transformation Classes
• Image Manipulation
• Reslicing -- vtkImageReslice
• Smoothing -- vtkImageGaussianSmooth
• Intensity Scaling -- vtkImageShiftScale
• Resampling -- vtkImageResample
• An Outline of the Algorithm and its
  Implementation
• vtkpxSimpleRegistration
• vtkpxSimpleLinearRegistration
• Two Key Helper Classes
• vtkpxJointHistogram
• vtkpxLinearTransform

5
The vtkAbstractTransform Branch

6
vtkAbstractTransform

• Abstract parentclass of all VTK transformations
• Defines interface (but not the implementation)
  for most common operations e.g.
• TransformPoint()
• Identity()
• Inverse()
• TranformDerivative()
• TransformNormalAtPoint()
• vtkHomogeneousTransform provides a generic
  interface for homogeneous transformations, i.e.
  transformations which can be represented by
  multiplying a 4x4 matrix with a homogeneous
  coordinate.

7
The vtkHomogeneousTransform Branch

8
Explicit 4x4 Matrix Definition

• We can explicitly specify an arbitrary 4x4 matrix
  transformation using vtkMatrix4x4 (in this case
  we specify a translation of (10,0,0))
• vtkMatrix4x4 mat
• mat Identity
• mat SetElement 0 3 10.0
• vtkTransform tr
• tr SetMatrix mat
• mat Delete

9
Talk Outline

• A review of the VTK Transformation Classes
• Image Manipulation
• Reslicing -- vtkImageReslice
• Smoothing -- vtkImageGaussianSmooth
• Intensity Scaling -- vtkImageShiftScale
• Resampling -- vtkImageResample
• An Outline of the Algorithm and its
  Implementation
• vtkpxSimpleRegistration
• vtkpxSimpleLinearRegistration
• Two Key Helper Classes
• vtkpxJointHistogram
• vtkpxLinearTransform

10
Reslicing an Image
F
r
r
Reference Image R (Static)
Transform Image T (Moving)

• In reslicing images we compute a transformation T
  which maps positions on the reference image to
  positions in the transform image. (e.g rFr)
• The standard image reslicing algorithms
  (vtkImageReslice) then takes the intensity at r
  T(r) and places it at position r in the new
  resampled image.
• Hence while F maps points from R to T, it maps
  intensities from T to R!!

11
vtkImageReslice

• vtkImageReslice is one of the most powerful
  classes in VTK.
• It can be used to reslice, resample and reorient
  images using different transformations and
  interpolation strategies.
• In its simplest complete usage it takes a number
  of parameters the most useful being
• The input image T to be resliced
• The information input image R that defines the
  output geometry Note that the output geometry
  can also be defined explicitly without using an
  information input
• The reslice transformation that maps R to T
  identity if not specified
• The interpolation strategy used to interpolate T
  (nearestneighbor, linear, cubic)

12
Image Reslice Example(From Lecture 7 of part I)

• This is the reslice, note that we use the image
  itself to define the reference
• geometry (not always a good idea!)
• vtkImageReslice resl
• resl SetInput tr GetOutput
• resl SetInformationInput tr GetOutput
• resl SetResliceTransform xform
• resl SetInterpolationModeToLinear
• resl Update

13
vtkImageShiftScale

• Performs Intensity Scaling
• Iout ( Iin shift ) scale
• Key Methods
• virtual void SetShift (float)
• virtual void SetScale (float)
• void SetOutputScalarType (int)
• void SetOutputScalarTypeToShort ()
• void SetOutputScalarTypeToFloat ()

14
vtkImageResample

• Resamples Image to a given resolution (can also
  be achieved using vtkImageReslice Identiy
  Transformation)
• Key Methods (axis0(x), 1(y), 2(z))
• SetAxisOutputSpacing (int axis, float spacing)
• virtual void SetInterpolate (int)
• virtual void InterpolateOn ()

15
vtkImageGaussianSmooth

• Implements simple smoothing with a Gaussian
  kernel
• Can specify kernel shape (in voxels)
• void SetStandardDeviations (float, float, float)
• And kernel size (also in voxels)
• SetRadiusFactors ( float, float,float)
• The radius factors determine how far out the
  gaussian kernel will go before being clamped to
  zero.
• Smoothing is a key step in multi-resolution
  algorithms (smooth resample)

16
Talk Outline

• A review of the VTK Transformation Classes
• Image Manipulation
• Reslicing -- vtkImageReslice
• Smoothing -- vtkImageGaussianSmooth
• Intensity Scaling -- vtkImageShiftScale
• Resampling -- vtkImageResample
• An Outline of the Algorithm and its
  Implementation
• vtkpxSimpleRegistration
• vtkpxSimpleLinearRegistration
• Two Key Helper Classes
• vtkpxJointHistogram
• vtkpxLinearTransform

17
The Registration Algorithm

• Given two images R (reference) and T (transform)
  estimate the transformation G that maximizes
• Similairity ( G(R), T )
• Key implementation choices
• Similarity Metric Correlation, Mutual
  Information
• Optimization Strategy Multi-resolution Hill
  Climbing
• Transformation Class Linear Affine

18
vtkpxSimpleRegistration

• This is the parent class for (potentially) a
  number of different approaches
• Derived from vtkProcessObject fairly minimal
  overhead and inherited functionality
• Typical use (note that a derived class is used!)
• Create i.e. regvtkpxSimpleLinearRegistrationNew
  ()
• Set Parameters
• reg-gtSetReference(img)
• Run the registration
• reg-gtRun()
• Get the Result
• Transformreg-gtGetTransformation()

19
vtkpxSimpleLinearRegistration

• Concrete derived class of vtkpxSimpleRegistration
  for estimating linear transformations
• Implements specific optimization strategy for
  linear transformation estimation
• Optimize()
• Uses vtkpxLinearTransform to store the
  transformation.
• Key Methods
• virtual void SetTransformModeToRigid()
• virtual void SetTransformModeToAffine()
• virtual vtkAbstractTransform GetTransformation()
  
• vtkSetObjectMacro(InitialTransform,vtkpxLinearTran
  sform)

20
The Multiresolution Approach

• Typical choice for both
• Computational Efficiency
• Avoidance of Local Minima
• If lowest resolution is 1.5 mm, and three levels
  are used
• Step 1 Blur both images and resample at 6 mm
• Optimize to Compute Tranformation G
• Step 2 Blur both images and resample at 3mm
• Use G as initial guess
• Optimize to Compute Transformation G
• Step 3 Resample at 1.5 mm (or use original
  resolution)
• Use G as initial guess
• Optimize to Compute Final Transformation G

21
A High-level View of the Implementation

• Set Parameters e.g. Images, Resolution etc.
• InitializeAll()
• For each level
• InitializeLevel (level)
• ComputeRegistration(level)
• Optimize()
• FinalizeLevel(level)
• FinalizeAll()

22
SetParameters

• Typical Minimum
• SetReferenceImage(vtkImageData )
• SetTargetImage(vtkImageData )
• SetSimilarityMeasure(int )
• e.g. SetSimilarityMeasureToNormalizedMutualInforma
  tion()
• Additional Parameters
• SetResolution(float)
• SetNumberOfLevels(int)
• SetNumberOfIterations(int)
• Registration Process is invoked using
• registration-gtRun()

23
InitializeAll()

• Nothing Special
• Check that Reference and Target Images exist
• Do some preliminary memory allocation

24
InitializeLevel(level)

• Compute current resolution and blurring kernel
  sigma
• Resample and Smooth Reference and Transform
  Images using the method
• vtkpxSimpleRegistrationResampleAndSmoothImage()
• To compute SmoothReference and SmoothTransform
  images respectively
• Given the maximum number of bins (e.g. 64) for
  the histograms adjust intensities of both
  SmoothReference and SmoothTransform in the range
  0..63 using the method
• vtkpxSimpleRegistrationCalculateNumberOfBinsAndA
  djustImage()
• Allocate memory for the Joint Histogram class
  (of type vtkpxJointHistogram)

25
ComputeRegistration(level)

• Calls the optimization routine at this level
• Implemented in the Optimize() method
• Key Sub-methods are
• Evaluate()
• Given current transformation evaluate the
  similarity between the images
• This in turn calls
• FillHistogram(img1,img2)
• ComputeSimilarity(histogram)

26
Finalize(level) FinalizeAll()

• Nothing Special
• Clean up code and tidying up!

27
Key Code Fragments

• ResampleAndSmoothImage()
• CalculateNumberOfBinsAndAdjustImage()
• Optimize() and Slow Hill Climbing
• Evaluate()
• Similarity(vtkpxJointHistogram)

28
ResampleAndSmoothImage() -- 1
int vtkpxSimpleRegistrationResampleAndSmoothImag
e( vtkImageData destination,
vtkImageData source, float baseresolution3,
float sigma, float res) if (sourceNULL
destinationNULL) return 0 float
sp3 source-gtGetSpacing(sp) float ori3
source-gtGetOrigin(ori) int dim3
source-gtGetDimensions(dim) float
gaussian3,resolution3 for (int
j0jlt2j) gaussianjsigmabaseresol
utionj resolutionjbaseresolutionjre
s // Compute Gaussian smoothing std in
voxels from mm float gp3 for (int
ia0ialt2ia) gpia(gaussiania/spia)
29
ResampleAndSmoothImage() -- 2
// Smooth then resample pipeline vtkImageGaussianS
mooth smoothvtkImageGaussianSmoothNew()
smooth-gtSetInput(source)
smooth-gtSetStandardDeviations(gp) vtkImageResamp
le resampvtkImageResampleNew() resamp-gtSetIn
put(smooth-gtGetOutput()) smooth-gtDelete() res
amp-gtInterpolateOn() resamp-gtSetDimensionality(3
) for (int ib0iblt2ib)
resamp-gtSetAxisOutputSpacing(ib,resolutionib) r
esamp-gtUpdate() destination-gtShallowCopy(resamp-
gtGetOutput()) resamp-gtDelete() return 1
30
CalculateNumberOfBinsAndAdjustImage() -- 1
// Semi-bad practice in this method we overwrite
the input with the output int vtkpxSimpleRegistrat
ionCalculateNumberOfBinsAndAdjustImage( vtkImag
eData image, int maxbin, short minv,short
maxv) float r2 image-gtGetPointData()-gtGetSc
alars()-gtGetRange(r,0) int range
int(r1 - r0 1), int nbins range, width
1 // Calculate number of bins to use if
(maxbin gt 0) while ((int)(ceil(range
/ (double)width)) gt maxbin) width nbins
(int)(ceil(range/(double)width))
minvint(r0) maxvint(r1) // Compute
Shift and Scale Parameters to force intensity to
0..nbins-1 float scalefloat(nbins)/range
float shift-r0/scale
31
CalculateNumberOfBinsAndAdjustImage() -- 2
vtkImageShiftScale scvtkImageShiftScaleNe
w() sc-gtSetInput(image)
sc-gtSetOutputScalarTypeToShort()
sc-gtSetShift(shift) sc-gtSetScale(scale)
sc-gtUpdate() image-gtShallowCopy(sc-gtGetOutput()
) sc-gtDelete() return nbins
32
Slow Hill Climbing(Studholme et al)

• If Optimization has 6 parameters pi(3
  translations, 3 rotations)
• Calculate similarity metric S at current estimate
  of p
• Compute Similarity S at pi dp and pi-dp for all
  parameters
• Take the best result and continue until best
  result is current position
• Reduce dp and repeat if desired
• Transformation stored as vtkpxLinearTransform.
  Parameters are stored as a column vector
  (translation_x,t_y,t_z,rotation_x,r_y,r_z,scale_x,
  s_y,s_z)
• Transform-gtGet(i) returns the value of
  parameter i
• Transform-gtPut(i,v) sets the value of parameter
  i to v

33
Optimize()
float vtkpxSimpleLinearRegistrationOptimize(int
maxdofs,float step,float old_similarity,int
level,int st) float similarity0.0,new_similari
ty0.0 for (int iter1iterltthis-gtNumberOfItera
tionsiter) int j 0,k 0
old_similarity new_similarity
this-gtEvaluate() // Evaluate Similarity using
current xformation for (int i 0 i lt
maxdofs i) float origthis-gtLinearTran
sform-gtGet(i) this-gtLinearTransform-gtPut(i
, orig step) similarity
this-gtEvaluate() if (similarity gt
new_similarity) new_similarity
similarity j i k 1
this-gtLinearTransform-gtPut(i, orig - step)
similarity this-gtEvaluate() if
(similarity gt new_similarity)
new_similarity similarity j i k -1
this-gtLinearTransform-gtPut(i, orig)
if (new_similarity gt
old_similarity) this-gtLinearTransform-gtPut(j,
this-gtLinearTransform-gtGet(j) k step)
old_similaritynew_similarity
else iterthis-gtNumberOfIterations1 // Break
Out return old_similarity
34
Evaluate()

• Evaluate() in vtkpxSimpleLinearRegistration
  simply calls the more complex Evaluate() in
  vtkpxLinearRegistration with the appropriate
  parameters e.g.
• float vtkpxSimpleLinearRegistrationEvaluate()
• this-gtLastSimilarityvtkpxSimpleRegistrationEv
  aluate(
• this-gtSampledReferenceImage,this-gtSampledTransfo
  rmImage,
• this-gtLinearTransform,1)
• return this-gtLastSimilarity

35
Evaluate() -- 2
The key to evaluate is to resample the target to
match the reference both in terms of position (as
defined by the current transformation) and number
of voxels. This is achieved using
vtkImageReslice. Then the vtkpxJointHistogram is
used to compute the appropriate similarity
metric. float vtkpxSimpleRegistrationEvaluate(v
tkImageData ref,vtkImageData targ,vtkAbstractTra
nsform trans, int interp,int reset)
vtkImageReslice reslvtkImageResliceNew()
resl-gtSetInput(targ) resl-gtSetInformationInput(
ref) resl-gtSetResliceTransform(trans)
resl-gtSetInterpolationMode(interp)
resl-gtUpdate() int dim3
resl-gtGetOutput()-gtGetDimensions(dim)
this-gtHistogram-gtFillHistogram(ref,resl-gtGetOutput
(),reset) float valthis-gtSimilarity(this-gtHist
ogram) resl-gtDelete() return val
36
Similarity(vtkpxJointHistogram)
float vtkpxSimpleRegistrationSimilarity(vtkpxJoi
ntHistogram histo) if (histoNULL)
histothis-gtHistogram // Note the use of
,- signes to ensure that optimum maximum
switch (this-gtSimilarityMeasure) case
0 return histo-gtCrossCorrelation()
break case 1 return
histo-gtMutualInformation() break
case 2 return histo-gtNormalizedMutualInfor
mation() break case 3 return
-histo-gtSumsOfSquaredDifferences() break

37
Talk Outline

• A review of the VTK Transformation Classes
• Image Manipulation
• Reslicing -- vtkImageReslice
• Smoothing -- vtkImageGaussianSmooth
• Intensity Scaling -- vtkImageShiftScale
• Resampling -- vtkImageResample
• An Outline of the Algorithm and its
  Implementation
• vtkpxSimpleRegistration
• vtkpxSimpleLinearRegistration
• Two Key Helper Classes
• vtkpxJointHistogram
• vtkpxLinearTransform

38
vtkpxJointHistogram

• Custom class to compute various intensity
  similarity metrics between pre-aligned images
  based on their joint histogram
• Histogram stored as a 2D matrix (see code for how
  this is allocated/deallocated and managed)
• Key methods
• int FillHistogram(vtkImageData ref,vtkImageData
  targ,int reset)
• float CrossCorrelation()
• float NormalizedMutualInformation()
• float SumOfSquareDifferences()
• See the code for all the details

39
vtkpxJointHistogram 2
int vtkpxJointHistogramFillHistogram(vtkImageDat
a ref,vtkImageData targ) int np1
ref-gtGetPointData()-gtGetScalars()-gtGetNumberOfTupl
es() int np2targ-gtGetPointData()-gtGetScalars()
-gtGetNumberOfTuples() // Images Must have the
same number of voxels! if (np1!np2) return
0 // Images Must of be of type short short
img1 (short) ref-gtGetPointData()-gtGetScalars(
)-gtGetVoidPointer(0) short img2 (short)
targ-gtGetPointData()-gtGetScalars()-gtGetVoidPointer
(0) this-gtReset() for (int i0iltnp1i)
// Add to Bin increments the value of
bin(i,j) by 1 this-gtAddToBin(img1,img2)
img1 img2 return 1
40
vtkpxJointHistogram 3
// See code for implementation of EntropyX()
etc float vtkpxJointHistogramMutualInformation(
) if (this-gtCheckHistogram()0) return
0.0 return this-gtEntropyX()
this-gtEntropyY() - this-gtJointEntropy() float
vtkpxJointHistogramNormalizedMutualInformation()
if (this-gtCheckHistogram()0) return
0.0 return (this-gtEntropyX()
this-gtEntropyY()) / this-gtJointEntropy()-1.0 f
loat vtkpxJointHistogramCrossCorrelation()
if (this-gtCheckHistogram()0) return 0.0
return this-gtCovariance() / (sqrt(this-gtVarianceX(
)) sqrt(this-gtVarianceY()))
41
vtkpxLinearTransform

• Given parameters such as 3 translations, 3
  rotations, 3 scales and 3 shears compute 4x4
  matrix
• Key methods
• Void Put(int i, float v)
• Set parameter(i) v
• float Get(i)
• Return the value of parameter i
• Derived from vtkTransform can function as a
  standard VTK transformation
• Again See the code for all the details

42
vtkpxLinearTransform -- 2
// From the header file vtkpxLinearTransform.h //
The transformation is parameterized as ///
Scaling along the x-axis, y-axis, z-axis as a
percentage i.e. 100 identity float
_sx,_sy,_sz /// Skew angle in the x-y ...
plane (in degrees) float _sxy,
_syx,_syz,_szy,_szx,_sxz /// Translation along
the x-axis (in mm) float _tx,_ty, _tz ///
Rotation around the x-axis (in degrees) float
_rx,_ry,_rz // Additional options for
mirroring of an image int FlipX int FlipY
int FlipZ
43
vtkpxLinearTransform -- 3
void vtkpxLinearTransformUpdateParameters()
// Update affine transformation Add shearing
vtkMatrix4x4 skewx vtkMatrix4x4New()
skewx-gtIdentity() skewx-gtSetElement(2,
1,tan(_szy(pi/180.0))) skewx-gtSetElement(1,
2,tan(_syz(pi/180.0))) vtkMatrix4x4 skewy
vtkMatrix4x4New() skewy-gtIdentity()
skewy-gtSetElement(2, 0,tan(_szx(pi/180.0)))
skewy-gtSetElement(0, 2,tan(_sxz(pi/180.0)))
vtkMatrix4x4 skewz vtkMatrix4x4New()
skewz-gtIdentity() skewz-gtSetElement(1,
0,tan(_sxy(pi/180.0))) skewz-gtSetElement(0,
1,tan(_syx(pi/180.0))) // Scales stored as
percentage multiply by 0.01 to convert to scale
factors and check for flip float sx_sx0.01
if (this-gtFlipX) sx-1.0sx float sy_sy0.01
if (this-gtFlipY) sy-1.0sy float sz_sz0.01
if (this-gtFlipZ) sz-1.0sz vtkTransform
trvtkTransformNew() tr-gtIdentity()
tr-gtPostMultiply() tr-gtConcatenate(skewx)
tr-gtConcatenate(skewy) tr-gtConcatenate(skewz)
tr-gtScale(sx,sy,sz) tr-gtRotateX(_rx)
tr-gtRotateY(_ry) tr-gtRotateZ(_rz)
tr-gtTranslate(_tx,_ty,_tz) skewx-gtDelete()
skewy-gtDelete() skewz-gtDelete()
this-gtSetMatrix(tr-gtGetMatrix()) tr-gtDelete()
this-gtvtkTransformInternalUpdate()
44
The End

• All Code and Lecture notes is online.
• The only way to learn programming is by
  programming!
• Hopefully this seminar provided some pointers as
  to where to start.