Katholieke Universiteit Leuven
Contents
1. 6 75 seconds 0 0000000000 0 0000000000 1 0000000000 0 0000000000 OUTPUT TO FILE MIRIToutput Figure 6 3 Typical MIRIT screen output in histogram mode MIRIT Version 97 08 histogram Type double Size 256 Header 128 bytes lt binary double histogram image data gt Thu Apr 17 18 14 06 1997 256 Figure 6 4 Header format of the MIRIT output file in histogram mode 26 s s s s s s pling factors as specified in MIRITimages The values in the floating image are then overwritten by the interpolated values at the transformed position in the reference image The transformation is defined by the transformation parameters specified in MIRITparams By default trilinear interpolation is used but nearest neighbour interpolation can be used as well which is to be specified in MIRITparams The reformatted reference image is written to the specified output file using 2 bytes per voxel short A 128 byte header describes the image size figure 6 6 REFORMATTING HAS STARTED Position 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 Interpolation 2 trilinear Done 9 13 seconds OUTPUT TO FILE MIRIToutput Figure 6 5 Typical MIRIT screen output in reformatting mode MIRIT Version 97 08 reformat Thu Apr 17 18 15 57 1997 Type short Size 256 256 100 Header 128 bytes lt binary2. into the set n Pr Ti2 pr in the reference image R Only those samples Pp Sra C Sp are retained for which pg falls inside the volume of R Let f pp and r pg denote the image intensity values in image F at position pr and in image R at the transformed position pg respectively The joint image intensity histogram Ha f r of the overlapping part of F and R is computed by binning the image intensity pairs f Pr r pn for all Pr Spa The number of bins np and ng to be used for the floating and reference image intensities respectively are parameters to be specified by the user Typically we use np ng 256 In order to do the binning efficiently the floating and the reference image intensities are first linearly rescaled between 1 np 1 and 1 ng 1 respectively The histogram bins H 0 r r 0 1 nr 1 and H f 0 f 0 1 npr 1 are not being considered when computing the mutual information criterion This allows to use the value 0 in either image for special purposes such as to mark regions in the image that one wants to exclude when computing the registration measure A rectangular region of interest and sampling factors along each axis separately can be specified in either image The original images are first resampled using trilinear interpolation within these regions Only the voxels inside the overlapping volume of these regions are considered when computing the joint image intensity histogram and the re
3. RefFile images registration MR img RefHeader 512 RefROI O 255 O 255 O 179 RefSampling 1 000000 1 000000 1 000000 RefRange 1 4094 4 FLOATING IMAGE FltOrientation 1 2 3 FltDimensions 256 256 100 FltSizes 0 937500 0 937500 1 550000 FltGantry 0 000000 FltFile images registration CT img FltHeader 512 F1tROI 0 255 O 255 0 99 FltSampling 1 000000 1 000000 1 000000 FltRange 1 4093 CRITERION Criterion 4 I X Y Interpolation 3 partial volume Criterion 4 RefNrOfBins 256 FloatNrOfBins 256 Interpolation 3 ResolutionLevels 2 2 2 NrOfResolutionLevels 2 ResolutionSteps 3 3 2 4 OPTIMISATION 4 Optimisation method 1 Powell OptimisationMethod 1 4 Line minimisation parameters LinminMaxBracketStep 1 00e 01 LinminAbsTolerance 1 00e 03 LinminFracTolerance 1 00e 02 LinminMaxNrOfIterations 20 Convergence parameters OptimFracTolerance 1 00e 05 OptimMaxNrOflterations 20 Powell parameters OptimOrder 563124000000 OptimBracketSteps 1 00 1 00 1 00 1 00 1 00 1 00 0 10 0 10 0 10 0 10 0 10 0 10 OptimAcceptNewDirection 0 Figure 7 1 Typical MIRIT output file in registration mode 29 Optimisation result InitialPosition 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 00000
4. f MakefileIBMPW2 x1C compiler for IBM Power 2 architecture or gt make f MakefileSUN CC compiler or gt make f MakefileHP CC compiler or gt make f MakefileSGI CC compiler in the Mirit src directory This creates the library Mirit lib libMiritX a and the executable Mirit bin MiritX with X the architecture name IBM IBMPW2 SUN HP or SGI The compiler flags have been set to generate optimized code 03 for each platform For SGI some flags ip have to be adapted to match the processor type See the manual pages for CC on your platform 4 4 Starting the program The MIRIT program has no fancy user interface but is designed to run off line reading its input from a number of input files and writing its output to an output file Assuming that the program executable Mirit resides in the directory home registration Mirit bin the program can be run from the working directory home registration Mirit experiment by the commands gt cd home registration Mirit experiment gt home registration Mirit bin Mirit MIRITinput The file MIRITinput which should be in the working directory specifies all input to the program All input filenames mentioned in this document are arbitrary Examples of the MIRIT input files can be found in the Mirit examples directory The format of the MIRITinput file is discussed in chapter 5 17 Chapter 5 Passing input to the program 5 1 Specifying the operation mode MIRITi
5. short reformatted image data gt Figure 6 6 Header format of the MIRIT output file in reformatting mode 27 Chapter 7 Interpreting the output file The output file formats for histogram and reformatting mode have been described in sections 6 3 and 6 4 above In this section we discuss the output file in case of registration and tracing mode 7 1 Registration mode The MIRIT output file corresponding to the above example files MIRITimages and MIRITparams for registration mode is shown in figures 7 1 7 2 and 7 3 This file is in Matlab 4 format which allows easy analysis of the results The contents of the file are Reference image The parameters of the reference image as specified by the user in the MIRITimages file The intensity range of the image RefRange as determined when rescaling the image is also displayed Floating image The parameters of the floating image as specified by the user in the MIRITimages file The intensity range of the image FltRange as determined when rescaling the image is also displayed Criterion The criterion parameters as specified by the user in the MIRITparams file Optimisation The optimization parameters as specified by the user in the MIRITparams file and the optimization result This includes FinalPosition the computed registered position rotation translation scaling and skew CriterionValues the values of the information theoretic measures and of the registration criterion at the
6. 0 0 0000000000 0 0000000000 1 4 World to Ref coordinate transformation RefToFltWorldToImageTransform 1 0666666667 0 0000000000 0 0000000000 0 0000000000 1 0666666667 0 0000000000 0 0000000000 0 0000000000 0 6451612903 0 0000000000 0 0000000000 0 0000000000 1 Image coordinate transformation RefToFltImageTransform 1 0393901484 0 0474539674 0 0510532017 0 0379654176 1 0246266937 0 1881715562 0 0347252798 0 1098039724 0 6342917952 0 0000000000 0 0000000000 0 0000000000 1 Figure 7 3 Typical MIRIT output file in registration mode cont 31 127 127 89 0000000000 124 124 89 0000000000 127 127 49 0000000000 5000000000 5000000000 5000000000 6436522953 1034553716 59 0000000000 0929722627 5116550000 5116550000 5000000000 9651062228 8885192662 4813921137 0000000000 5000000000 5000000000 5000000000 7653814991 10 37 0000000000 8752806350 6200741141 7 2 Tracing mode The MIRIT output file corresponding to the above example files MIRITimages and MIRITparams for tracing mode is shown in figure 7 4 This file is in Matlab 4 format which allows easy analysis of the results The file summarizes all user specified image related input parameters together with the image intensity range as determined when rescaling the images the trace parameters and the computed criterion values along the trace These values are the number of floating image sam
7. 0 workstation under AIX 4 1 3 using IBM s x1C compiler and has been ported to SUN HP and SGI using the CC compiler Porting to other platforms such as PC should not pose too many problems due to the fact that the code is completely self contained does not require external libraries except for the standard C library and does not use the template class mechanism 4 2 Installation The program is distributed as a compressed tar file Mirit tar gz First move this file to the directory in which the program has to be installed keeping in mind that extracting the code from this file will itself create a subdirectory Mirit The code can be extracted from this file by gt gunzip Mirit tar gz gt tar xvf Mirit tar This will create the following directory structure Mirit include MIChrono h MIImageGrid h MIManager h MIMeasure h MIOptimisation h MIParams h MIPrint h MIReformat h MIShrtImage h MITrace h MITransform h MITypes h Mirit src MIChrono C MIImageGrid C MIManager C MIMeasure C MIMirit C MIOptimisation C MIParams C MIReformat C MIShrtImage C 16 MITrace C MITransform C MakefileHP MakefileIBM MakefileIBMPW2 MakefileSGI MakefileSUN Mirit bin Mirit lib Mirit man MiritManual ps gz this document Mirit examples MIRITinput MIRITimages MIRITparams 4 3 Compilation The program can be compiled by running depending on your platform gt make f MakefileIBM x1C compiler or gt make
8. 000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 1 NumberOfSteps 11 Trace values 4 Npoints H float H ref H float ref TraceValues 4587506 4 902859046003 6 025869005183 10 375254270044 4589021 4 902324662568 6 028251727764 10 374780710069 sls Criterion 7 446526218858 7 444204319737 Figure 7 4 Typical MIRIT output file in tracing mode 33 Bibliography 1 F Maes A Collignon D Vandermeulen G Marchal P Suetens Multimodality image registration by maximization of mutual information IEEE Trans Medical Imaging vol 16 no 2 pp 187 198 April 1997 2 J West M Fitzpatrick et al Comparison and evaluation of retrospective intermodality brain image registration techniques Jourbal of Computed Assisted Tomography vol 21 no 4 pp 554 566 1997 3 W H Press B P Flannery S A Teukolsky and W T Vetterling Numerical Recipes in C Second Edition Cambridge England Cambridge University Press 1992 chapter 10 pp 412 419 4 Matlab The Language of Scientific Computing The MathWorks Inc 24 Prime Park Way Natick MA 01760 1500 USA 34
9. 00000 des FinalPosition 10 0522495309 3 1595130675 2 0918915841 6 9018589962 1 0915816053 18 2176739384 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 l CriterionValues Initial Final 4589046 4372616 Npoints 4 902308020691 4 982842932259 H float 4 339177262774 4 039391170488 H float ref 6 041564058499 5 981530181951 H ref 5 478433300581 5 038078420180 H ref float 10 380741321272 10 020921352439 H float ref 0 563130757918 0 943451761771 I float ref 0 102912525629 0 172094063554 ECC float ref 0 093209432601 0 157727493312 UC float ref 9 817610563355 9 077469590668 Rho float ref 7 436869242082 7 056548238229 Criterion 1 NrOfIterations 7 29 NrOfEvaluations 595 133 728 TimeInSeconds 275 21 889 69 1164 90 ResolutionFactors 18 0 1 0 TRANSFORMATION FltToRef 4 Transformation parameters FltToRefParameters 10 0522495309 3 1595130675 2 0918915841 6 9018589962 1 0915816053 18 2176739384 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 1 Flt to world coordinate transformation FltToRefImageToWorldTransform 0 9375000000 0 0000000000 0 0000000000 119 5312500000 0 0000000000 0 9375000000 0 0000000000 119 5312500000 0 0000000000 0 0000000000 1 5500000000 76 7250000000 0 0000000000 0 0000000000 0 0000000000 1 0000000000 1 4 World coordinate transformation FltToRefWorldTransform 0 99
10. 0683506925 3 1541507013 2 0849665884 6 8953720978 1 1058616719 18 2236738732 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 OUTPUT TO FILE MIRIToutput m 7 421546973235 7 403797011246 7 403797011246 7 403797011246 Figure 6 1 Typical MIRIT screen output in registration mode 25 TRACING HAS STARTED Initial position 5 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 Step sizes 1 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 Number of steps 11 0 N 1039651 H 9 05637407 0 59871847 ECC 0 12402128 80 1 N 1039351 H 9 05759336 0 59697297 ECC 0 12366645 74 2 N 1040029 H 9 05825980 0 59340902 ECC 0 12296506 T1 8 N 1040029 H 9 04950813 0 58294867 ECC 0 12103842 89 9 N 1039351 H 9 05225641 0 58010088 ECC 0 12044837 86 10 N 1039651 H 9 05171481 0 57742469 ECC 0 11993277 88 OUTPUT TO FILE MIRIToutput m HISTOGRAM COMPUTATION HAS STARTED Reference image Floating image Interpolation Figure 6 2 Typical MIRIT screen output in tracing mode 256 bins 256 bins 3 partial volume Position 0 0000000000 0 0000000000 1 0000000000 0 0000000000 Done 0 0000000000 0 0000000000 1 0000000000 0 0000000000
11. 56 180 lt reference image dimensions gt 0 9765625 0 9765625 1 0 lt reference image voxel sizes gt 0 0 lt reference image gantry gt images registration MR img lt reference image file name gt 512 lt reference image file header gt O 255 0 255 O 179 lt reference image ROI gt 1 0 1 0 1 0 lt reference image sampling factors gt 123 lt floating image orientation gt 256 256 100 lt floating image dimensions gt 0 9375 0 9375 1 55 lt floating image voxel sizes gt 0 0 lt floating image gantry gt images registration CT img lt floating image file name gt 512 lt floating image file header gt 0 255 0 255 O 99 lt floating image ROI gt 1 0 1 0 1 0 lt floating image sampling factors gt The number of images specified should always be 2 The first image specified is taken to be the reference image the second image specified is the floating image Samples are taken from the floating image and transformed into the reference image to compute the joint image intensity histogram of the overlapping volume of both images The speed of the method is increased by selecting the smaller image of the two to be the floating image Each of the images is specified by the following parameters in this order Dimensions along the x column y row and z slice axis respectively Orientation of the x y and z axis respectively as defined in section 1 1 Voxel sizes in mm along the x column y row and z slice
12. 78145425 0 0364468009 0 0551159642 6 9018589962 0 0455558087 0 9836416260 0 1742808650 1 0915816053 0 0478623766 0 1764108339 0 9831522825 18 2176739384 0 0000000000 0 0000000000 0 0000000000 1 0000000000 l Figure 7 2 Typical MIRIT output file in registration mode cont 30 4 World to Ref coordinate transformation FltToRefWorldToImageTransform 1 0240005243 0 0000000000 0 0000000000 0 0000000000 1 0240005243 0 0000000000 0 0000000000 0 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 1 4 Image coordinate transformation FltToRefImageTransform 0 9579019608 0 0349889289 0 0874800583 0 0437335763 0 9442959609 0 2766185890 0 0448709781 0 1653851568 1 5238860379 0 0000000000 0 0000000000 0 0000000000 1 TRANSFORMATION RefToFlt 4 Transformation parameters RefToFltParameters 10 1725523080 5 9651062228 1 0000000000 0 0000000000 1 4 Flt to world coordinate transformation RefToFltImageToWorldTransform 2 7433602007 1 8885192662 1 0000000000 0 0000000000 2 6140571097 18 4813921137 1 0000000000 0 0000000000 0 9765620000 0 0000000000 0 0000000000 0 0000000000 0 9765620000 0 0000000000 0 0000000000 0 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 l 4 World coordinate transformation RefToFltWorldTransform 0 9978145425 0 0455558087 0 0478623766 0 0364468009 0 9836416260 0 1764108339 0 0551159642 0 1742808650 0 9831522825 0 000000000
13. Katholieke Universiteit Leuven Departement Elektrotechniek Afdeling ESAT MI2 Laboratorium voor Medische Beeldvorming ESAT Radiologie Kardinaal Mercierlaan 94 B 3001 Heverlee Belgium TECHNISCH RAPPORT TECHNICAL REPORT Multimodality Image Registration using Information Theory MIRIT Version 97 08 Manual Pages Frederik Maes Dirk Vandermeulen and Paul Suetens September 1997 Nr KUL ESAT MI2 9708 Cr ij CA E Ej 57 ER ES e ED CEA Dj EJ AO E E ES Boo O O E D Multimodality Image Registration using Information Theory M LR I T Version 97 08 Manual Pages Frederik Maes Dirk Vandermeulen Paul Suetens Katholieke Universiteit Leuven Medical Image Computing ESAT Radiology UZ Gasthuisberg Herestraat 49 B 3000 Leuven Belgium E mail Frederik MaesQuz kuleuven ac be Frederik Maes is Aspirant of the Fund for Scientific Research Flanders Belgium Abstract This document describes the MIRIT software developed at the Laboratory for Medical Imaging Research of the K U Leuven for automated affine registration of multimodal medical images MIRIT implements the registration algorithm described in 1 The method is based on maximization of the mutual information between corresponding voxel intensities in the two images to be registered Becaus
14. an be determined such that the elements Roo cos cos dz and R22 cos x cos y of R are positive All rotation parameters Qi can then be assumed to be in the range 7 2 7 2 These conditions can be expressed as 1 Ro Ego 21 EE E gt 0 Sx Sy Sz 1 gz dy Rog Y Lp 8 p 899995 Toy gt 0 Sy Sy Sz Substituting the expressions for g these conditions are equivalent to 1 g Jy Gz sign ss 1 7 Loo A 2 m LAN D gt 0 el syl sz Ju Ja 9 1 9 9 Gz sign s D Lao TAE Loa A no 20 Sz E s2 from which the appropriate signs can be derived to be sign s sign S sign s sign 51 sign 2 sign D sign s sign S3 with 1 Jy Jy Jz 1 Loo Loi Loz Sz syl sz Ju 91 9 1 9 9 Gz LE ds CE ECE SEG pos Sz syl s2 The sign S1 sign S2 sign D factor in the expression for sy assures that s sy sz has the same sign as D Finally the rotation parameters can be derived from R Py asin Ro2 oo R be asin a Z a y Roa rco Chapter 2 The registration measure 2 1 Histogram computation The mutual information of the images to be registered at the current registration position is evaluated from the joint histogram of the image intensities of the overlapping volume of the images at this position Let Sp pr be the set of samples pr taken from the floating image F and transformed by the image to image transformation T2
15. ansformed position po of the origin pi o of image 1 in image 2 and the differential changes Ap in position in image 2 for differential changes Ap in position in image 1 Pio 0 0 0 1 P2 0 T12 P1 0 Ti2 4 P2 0 Ap Tis pi o Api Ap Tis Api Api 1 00 1 T Aspa Ti2 1 AyP1 0 10 1 gt AyP2 Tia 2 Api 0 0 1 1 gt A 2 Tia 3 Given the transformation P2 Tis pio Ti2 4 of the origin of the floating image into the reference image the transformed position p of all other points pi of the floating image is computed incrementally while scanning the floating image and incrementing pa with Tia 1 Tia 2 or Tia 3 for a step in the floating image along the x y or z image axis respectively 1 6 Computing the parameters of the affine transformation Given the 12 affine transformation parameters t gi Si the affine transformation matrix A can be straightforwardly computed by the matrix multiplications described above Here we are interested in the inverse problem decomposing a given affine transformation matrix A into A TxRxGxS with T R G and S of the form defined above i e computing from A This is useful if the transformation A of image 1 into image 2 is known and we wish to use the parameters of the inverse transformation A to reformat image 1 along the grid of image 2 The translation parameters tz ty t of the transformation are simply given by ti Asa
16. asured in radians by 1 0 0 0 pa 0 cos y 0 0 10 sin y 1 0 0 0 0 1 assuming that y is positive for rotation of the top of the gantry caudally from patient s head towards patient s feet as is the case on Siemens scanners The image to world transformation matrix T w is thus given by w Tiw P Ti 0Ox Tx Vx C Va 0 0 Uz Co E 0 cos y vy 0 0Uy Cy DOS 0 sin y vy v Uz Cz 1 1 0 0 0 1 1 3 The affine world to world transformation The affine transformation that transforms world coordinates 4 in image 1 into world coordinates Wa in image 2 is defined by the 4 x 4 matrix A Wy Aw with A of the form A Tx Rx GxS 1 2 with S G Rand T being 4x 4 matrices respresenting scaling skew rotation and translation respectively The affine transformation is defined by 12 parameters 3 translation distances t 3 rotation angles i 3 skew factors g and 3 scale factors s from which the matrices T R G and S are computed as described below Alternative definitions of A are possible such as for instance A T x Gx Sx R but definition 1 2 has the advantage that S can be interpreted directly as correction factors to be applied to the voxelsizes of image 1 such that S V represents the correctly calibrated voxelsizes of this image This interpretation is important when the voxelsizes of image 1 are not precisely known Moreover if image 1 is a CT image with orientation 1 2 3 for which a gantry tilt was applied whe
17. axis respectively Gantry tilt in degrees This is the angle between the imaging plane and the plane orthogonal to the slice axis for CT images The gantry tilt is assumed to be positive if the top of the gantry is tilted towards the patient s feet caudally and is negative if the top of the gantry is tilted towards the patient s head cranially This number should always be 0 0 for MR and PET images File name including directory path and Header length in bytes The image data are read from the specified file name skipping the specified number of bytes for the file header MIRIT requires that the image data are stored contiguously row by row and slice by slice using 16 bits short per voxel The program exits if the specified file can not be opened or if the number of voxels read from file does not match the specified dimensions 19 Region of interest The program allows for a rectangular region of interest ROI to be defined in either image Only voxels inside this ROI will be taken into account when computing the registration criterion The ROI is specified by the first and last image coordinates Emin Umar Ymin Ymar 2min Zmaz Of the region along each of the axes These coordinates have to be specified even if the full image is to be considered Sampling The program allows to first resample either image within the specified ROI before registration This allows to artificially increase the resolution of the reference image for higher precisi
18. d s representing scaling along the x y and z axis respectively Ss 0 0 0 0 sy 00 S 0 0 s 0 1 6 0 0 0 1 The overall expression for A A j as function of the affine transformation parameters ti di gi si is given by Aoo Sm COSy cos gy cos sinz Ao sy cos sin gz SiNy gr Jz COSy COS Ago Sz Jx COSy COS Siny Apa tz Aio Sa Sin Siny cos cos sin g sin siny sin cos cos Aj sy sin Siny Sin cos cos gz sin cosy gz gz sinz SINy cos cos sin Ai s sin cosy gz sing siny cos cos SIN A13 ty A20 s 4 cos sin cos sing sin gy cos sin SIN sing cos A21 Sy c0S Siny sin sing COS gz COS COSy 9x gz COSz SINy cos sin sin A22 S C0Sz COS 9 COS SIN cos Sin sin Az tz A31 0 A32 0 1 4 The image to image transformation The affine geometrical transformation of image coordinates p measured in voxels in image 1 into image coordinates pa in image 2 is defined by W A W Tiw 2 P2 A Ti w 1 P1 Pa Qe A Tiw 1 Pi pa Ti2 Pi 1 7 with Ti w and T w 2 being the image to world coordinate transformation matrices for images 1 and 2 respectively A the affine world to world transformation matrix mapping world coordinates in image 1 into world coordinates in image 2 and Ti2 the affine i
19. e initial position as specified in the MIRITparams file The trace is computed symmetrically around this initial position such that there are 2 times the number of steps plus 1 values in the trace and that the center value in the trace corresponds to the specified initial position During program execution the initial position and the current step along the trace are displayed together with values for the number of floating image sample points in the overlapping volume N the joint entropy H the mutual information I and the entropy correlation coefficient ECC as well as the CPU time required to compute these values figure 6 2 6 3 Histogram mode In this mode the joint histogram of the overlapping parts of both images is computed at the specified registration position and written to file The input images are treated as in registration mode The position for which the histogram is to be computed is defined by the registration parameters specified in the MIRITparams file The histogram image data are stored contiguously using 8 bytes per entry double The number of columns in the histogram image equals the number of bins specified for the reference image while the number of rows in the histogram image equals the number of bins specified for the floating image The number of bins is specified in the MIRITparams file A header of 128 bytes is added to the histogram output file describing the histogram image size figure 6 4 6 4 Reforma
20. e minimization algorithm is based on the routines linmin mnbrak and brent in 3 chapter 10 The parameters that determine convergence and that have to be specified by the user are e the initial step taken when bracketing the minimum i e 1 0 in the routine linmin in 3 pp 419 This is a parameter to be specified by the user as for instance a step of 1 0 has a different impact for scaling and skew than for translation and rotation Default values are 1 0 mm for translation 1 0 degree for rotation and 0 1 for scaling and skew These initial step values are decreased by a factor of 2 after each Powell iteration in order to try to minimize the number of evaluations required for the line minimization as the optimization nears the optimum e the maximal step taken when bracketing the minimum i e GLIMIT in the routine mnbrak in 3 pp 400 401 Default value is 10 for all directions e the fractional precision to which the minimum is isolated i e tol in the routine brent in 3 pp 404 405 Default value is 10 for all directions 13 e the absolute precision to which the minimum is isolated This parameter is not included in the routine brent in 3 pp 404 405 It is included here to declare convergence when the bracket is smaller than say 0 01 mm or 0 1 degrees Default value is 107 for all directions e the maximum allowed number of iterations i e ITMAX in the routine brent in 3 pp 404 405 Default value is 20 3 2 T
21. e no limiting assumptions are made regarding the image content and the nature of the relation between corresponding voxel intensities this criterion is very general and powerful allowing for completely automated and robust registration in a variety of applications without prior segmentation or other pre processing The accuracy of the method has been validated to be subvoxel 1 2 for rigid body registration of CT MR and PET brain images compared to the stereotactic reference registration solution MIRIT also allows to compute one dimensional traces of the registration criterion through registration space and to reformat the images at the registered position The MIRIT software and this manual is distributed without warranty and the software is to be used for experimentation only All questions or feedback relating to this software or requests to obtain the MIRIT software should be addressed to Frederik Maes at the address mentioned above This document has two parts Part one describes some implementation specifics of the program the second part describes how to install and use it Part I Implementation issues Chapter 1 The geometrical transformation 1 1 Coordinate frames and axis orientation With each of the images to be registered is associated an image coordinate frame in voxel units with its origin in the first voxel of the image the x coordinate corresponding to the column number the y coordinate to the row number and the z coord
22. for convergence gt 20 lt Brent maximum number of iterations gt 1e 5 lt Powell relative tolerance for convergence gt 20 lt Powell maximum number of iterations gt 0 lt Powell update direction matrix if new direction is found gt The values given here are the default values for these parameters For most cases it is not needed to change them However if more than one resolution level is used the maximum number of Powell iterations is better set to 5 instead of 20 to increase speed performance 22 Chapter 6 Functionality When the program is started using gt home registration Mirit bin Mirit MIRITinput the files MIRITimages and MIRITparams as specified in MIRITinput are read and the images are loaded from file The program mirrors its input on screen MIRIT input file MIRITinput MIRIT images file MIRITimages MIRIT parameter file MIRITparams MIRIT mode 1 registration MIRIT output file MIRIToutput m LOADING IMAGES Reference image Orientation 1 2 3 Dimensions 256 256 180 Voxelsizes 0 976562 0 976562 1 000000 Gantry 0 000000 File images registration MR img Header 512 ROI O 255 O 255 O 179 Sampling 1 000000 1 000000 1 000000 Floating image Orientation 1 2 3 Dimensions 256 256 100 Voxelsizes 0 937500 0 937500 1 550000 Gantry 0 000000 File images registration CT img Header 512 bytes ROI 0 255 0 255 O 99 Sampling 1 000000 1 000000 1 000000 The progra
23. from which the translation matrix T can be constructed Once T has been found the scale and skew parameters can be most easily computed from L T71x A R G x S using the fact that the rotation matrix R is orthonormal RT R The product P of L and its transpose LT is therefore independent of R L TUxA Rx Gx S Plis DUELE R Gx S Rx GxS ST x Gl xR R G S S GUxG S I The 6 skew and scale parameters can be computed from P by solving for g and s the set of 6 non linear equations S G xGxS P with s2 1 9 Sq Sy 9y 9x 92 Sx Sz Jx 0 S Gl Gag Sesu Iy 9092 sy 1 92 1 9 Sy S2 92 1 9 0 Sx Sz Ja sigo 92 sx 92 0 0 0 0 1 using the additional relation that Sq Sy Sz det L Solving these equations yields the following expressions for g and s with D det L P2 2 12 P Sy 11 Pro Po1 P22 Po2 Pi2 s2 Poo Coto Pei 25 92 P2 2 02 S P2 a Ss vs ls 57 Is Vs 2 s Jr n EH Pyi P 5 Po2 Pi2 Iy Pa D E etus CER 2 185 de P2 D Sx sign s s Sy sign sy Syl sz sign sz sz Y sign sy Ja Jy sign sz 9y gz sign Sz gz in which the sign of the scale factors has yet to be determined From the computed t s and gi the matrices T S and G can be reconstructed and the rotation matrix R can be derived as RT x Ax G 8 L 87 G7 using expressions 1 6 and 1 5 for S and G The sign of the scale factors c
24. gether by assigning them the same order number For instance the parameter order specification 5 6 3 1 2 4 7 7 8 0 0 O will force the scaling parameters in the x and y direction to vary simultaneously such that if they had the same value initially they will have the same value after optimization Experiments have indicated that robustness is increased by first optimizaing the in plane parameters Le tz ty and for an axial image before optimizing the out of plane parameters tz s and dy and that the non rigid parameters s and g are best optimized after the rigid parameters t and Qi A reasonable specification for the parameter order in most cases is thereforeb 6 312400 0 0 0 O for rigid body registration and 5 6 3 1 2 4 7 8 9 10 11 12 for full affine registration 20 Step size The first parameter on this line the number of steps parameter is used in tracing mode only while the next 12 parameters the step size parameters one for each affine parameter in the same order as used to specify the initial position are used both in tracing and registration mode In registration mode the step sizes determine the size of the first step in the bracketing phase of the line minimization routine see section 3 1 Default bracketing step values are 1 0 for rotation and translation and 0 1 for scaling and skewing and these are decreased during optimization If the specified step size parameter is zero the defualt value is used If the specif
25. gistration criterion Subsampling the floating image aims at increas ing speed performance see also section 3 3 while supersampling the floating or the reference image aims at increasing registration precision 2 2 Interpolation While all samples pr are taken at grid points of the floating image F their transformed positions pg will in general not coincide with a grid point of the reference image R and interpolation of the reference image is needed to obtain the image intensity value r pg at this position Three interpolation methods have been implemented see figure 2 1 Nearest neighbour NN interpolation The nearest neighbour of pg on the grid of R is found If this neighbour is inside the volume of R its intensity value is assigned to r pg and the histogram is updated by adding 1 to the entry corresponding to f Pr r pn Nearest neighbour interpolation is generally insufficient to guarantee subvoxel accuracy as it is insensitive to translations up to one voxel Trilinear TRI interpolation The 8 nearest neighbours of pg on the grid of R are found If all these neighbours are inside the volume of R their intensities are linearly interpolated and this 10 value is assigned to r pr The histogram is updated by adding 1 to the entry corresponding to Fr r pr Linear interpolation may introduce new intensity values which are originally not present in the image R leading to unpredictable changes in the marginal histogram of t
26. he Powell algorithm The implementation of the Powell algorithm is based on the routine powel1 in 3 chapter 10 pp 417 418 The parameters that determine convergence and that have to be specified by the user are e the order in which the parameters are optimized The direction matrix is initialized with unit vectors in each direction The order in which these directions are considered may greatly influence optimization robustness due to the fact that the image resolution is generally not the same along all axes e the maximum allowed number of iterations i e ITMAX in the routine powell in 3 pp 417 418 Default value is 20 e the fractional tolerance in the function value i e ftol in the routine powell in 3 pp 417 418 Default value is 1075 e a flag indicating whether or not the direction matrix should be updated when a new direction has been determined New directions are considered according to the criteria discussed in 3 pp 416 417 Default behavior is not to update the direction matrix 3 3 Multiresolution optimization Because the time required to evaluate the registration criterion is proportional to the number of samples taken from the floating image and transformed into the reference image to construct the joint image intensity histogram large speedups are possible by starting the optimization at a lower resolution by subsampling the floating image and by switching to a higher resolution by including more samples a
27. he reference image for small changes in the registration parameters Partial volume PV interpolation The 8 nearest neighours of pg on the grid of R are found For all neighbours n i 1 2 8 that fall inside the volume of R the histogram is updated by adding wi to entry f pi r n wi being the linear interpolation weight associated with neighbour i 0 w 1 To avoid having to test for each neighbour separately whether it falls inside the volume of R R is expanded into Ry by appending a one voxel wide border to the reference volume at all sides which is assigned the value zero The histogram is updated for each sample for which Pr falls inside the volume of Ry Because histogram entries with reference image zero are not taken into account when computing the registration criterion only those neighbours that fall inside the volume of R are effectively being considered Because the histogram is build up as the sum of fractions w that vary smoothly with the registration parameters the histogram and thus the mutual information measure itself is a continuous function of nA n3 arg minn d pg ni na NN PR r Pr r na Ha f gr r Pg 1 ni n2 na na Y wii 1 TRI PR r BR D wi r ni ue TS Ha f Pr r PR 1 n n2 na na PV Pr 25 wir 1 W3 Wa Vi Ha f Pr r n wi ni n2 Figure 2 1 Graphical illustration of NN TRI and PV interpolation in 2 D NN and TRI interpolation find the reference image inte
28. ied step size parameter is non zero the initial bracketing step is set to the step size parameter if this is larger than 0 001 and smaller than 10 0 otherwise the default value is used In tracing mode the step size parameters define the direction of a one dimensional line in parameter space Values of the registration criterion are computed along this trace direction at intervals of a step size in each parameter The trace is computed symmetrically around the initial position and has a total length of 2 times the number of steps plus 1 Criterion parameters These parameters specify the registration measure 4 lt registration criterion gt 256 lt number of reference image histogram bins gt 256 lt number of floating image reference bins gt 3 lt type of interpolation in reference image gt 2 lt number of resolution levels for column dimension gt 2 lt number of resolution levels for row dimension gt 2 lt number of resolution levels for slice dimension gt 3 lt resolution factor for column dimension gt 3 lt resolution factor for row dimension gt 2 lt resolution factor for slice dimension gt These are explained below Registration criterion This parameter specifies the information measure that is computed from the joint image intensity histogram and used as the matching criterion The parameter is to be specified as follows H F R the joint entropy of both images H F R the conditional entro
29. imal registration position of the images specified by MIRITimages using the registration parameters specified by MIRITparams and writing the registration results to file MIRIToutput e 2 tracing computes a trace of the registration criterion through parameter space for the images specified by MIRITimages using the trace parameters specified by MIRITparams and writing the trace values to file MIRIToutput e 3 histogram computes the joint intensity histogram of the images specified by MIRITimages at their relative position specified by MIRITparams and writing the histogram image to file MIRIToutput 18 e 4 reformatting reformats the reference image specified by MIRITimages along the grid of the floating image at the relative position of the images specified by MIRITparams and writing the reformatted image to file MIRIToutput If any other number is specified the program exits 5 2 Image data specification MIRITimages The program reads the images from user specified files assuming that the image data are stored contigu ously row by row and slice by slice using 16 bits per voxel short The image data file may have a header and the user can specify the number of bytes to skip from the beginning of the file The dimensions and voxelsizes of the images have to be specified in the file MIRITimages This file has the following contents 2 lt number of images to be registered gt 123 lt reference image orientation gt 256 2
30. inate to the slice number assuming that the data are stored row by row and slice by slice row major order With each of the images is also associated a world coordinate frame in millimeter units with its origin positioned in the center of the image volume defined in image coordinates by Z cs y cz d 1 2 with d dz dy d being the dimensions of the image in directions x y and z respectively The axes of the world coordinate frame are defined with respect to the patient the x axis running from patient s right to patient s left R L the y axis from patients s anterior to patient s posterior A P and the z axis from patient s inferior to patient s superior I S The orientation of each of the image coordinate axes column row and slice direction with respect to the world or patient coordinate frame has to be specified for each of the images This allows to compute image to image transformations of for instance transversal CT images to sagittal MR images Each axis can have 6 possible orientations as defined in table 1 1 The image orientation d oz oy 0 is specified by 3 numbers with values 1 1 2 2 3 or 3 The default orientations on Siemens scanners for CT and axial sagittal and coronal MR images are depicted in figure 1 1 From patient s to patient s right left 1 left right anterior posterior posterior anterior inferior superior superior inferior Table 1 1 Definition of image axis orien
31. initial and at the final position NrOfIterations the number of Powell iterations at each resolution level and overall NrOfEvaluations the number of criterion evaluations at each resolution level and overall TimeInSeconds the CPU time at each resolution level and overall ResolutionFactors the total floating image subsampling factors at each resolution level Transformation the transformation matrices for the floating to reference and for the reference to flaoting image transforms including the affine transformation parameters the image to world transforma tion matrix T 1 the affine world to world transformation matrix A the world to image transfor mation matrix Tos and the overall image to image transformation matrix 719 qz 2 A Ti wi For instance when points are being indicated in the reference image their coordinates in voxels have to be multiplied with the RefToFltImageTransform matrix to find the corresponding points in the floating image while its inverse matrix the FltToRefTransform has to be used to transform points in the floating image into the reference image The image to image matrices take both the differences in resolution and orientation as the difference in relative position of the images into account 28 MIRIT Version 97 08 registration Tue Jul 29 12 06 08 1997 REFERENCE IMAGE RefOrientation 1 2 3 RefDimensions 256 256 180 RefSizes 0 976562 0 976562 1 000000 RefGantry 0 000000
32. ling factors applied from one level to another is specified by the user for each dimension of the floating image If resolution levels 1 1 1 are specified no additional resolution levels will be considered and the complete optimization is done at full resolution If resolution levels 2 2 1 are specified together with resolution factors 2 2 2 2 resolution levels will be used The floating image is first subsampled by a factor of 2 along the column and row direction but not along the slice direction The optimization starts using this subsampled image and proceeds with the full resolution image If resolution levels 3 3 2 and resolution factors 2 2 2 are specified 3 resolution levels will be used with subsampling factors of 4 4 2 2 2 2 and 1 1 1 respectively In the example above 2 resolution levels are used one using subsampling factors 3 3 2 and the other at full floating image resolution Optimization parameters These parameters are used in registration mode only They specify the convergence parameters for the Powell algorithm and for the Brent line minimization routine Re ferring to the example file above and the description of the optimization parameters in sections 3 1 and 3 2 the interpretation of the parameters is as follows 1 lt optimization method should be 1 Powell gt 10 0 lt Brent maximum bracket step gt 1e 4 lt Brent absolute tolerance for convergence gt 1e 3 lt Brent relative tolerance
33. m may run in one of 4 different operation modes as specified by the operation mode parameter in the MIRITinput file described above The functionality of each of these modes is described hereafter 6 1 Registration mode The reference and the floating image are read from file as specified by the parameters in the file MIRITimages The floating and reference images are resampled within the specified ROI The world 23 coordinate frame for the cropped and resampled images is the same as the one assigned to the original images Both the reference and the floating image intensities are rescaled linearly to the range 1 Nbins 1 with Nbins being the number of histogram bins as specified in the MIRITparams file for the reference and floating image respectively The images are initially positioned according to the position parameters specified in the MIRITparams file which are echoed to screen Registration then involves optimizing the selected matching criterion During optimization the user can monitor the behaviour of the program from the output on the screen figure 6 1 At each iteration the criterion is optimized along a set of one dimensional lines specified by Direction These directions are initialized with unit vectors in each of the parameters seperately and are considered in the order as specified by the order parameters in MIRITparams as discussed above Bracket and Brent denote the bracketing and Brent steps of the line minimization alg
34. mage to image transformation matrix mapping image coordinates in image 1 into image coordinates in image 2 T1 is thus given by Tis Ti EAT ua 1 8 T1 incorporates the differences in orientation resolution and relative position of the images 1 and 2 1 5 Implementation of the image to image transformation One of the images to be registered is selected as the floating image from which samples are taken and transformed into the volume of the other image the reference image The image to image transformation is computed from the floating image into the reference image The image to world coordinate transformation matrices T y 1 and T w 2 are constructed for both the floating and the reference image respectively using expression 1 1 These matrices are independent of the actual geometrical transformation between both images and are fixed given the image dimensions voxel sizes and axes orientation For a given set of affine registration parameters tz ty tz Ox Oy Oz Gas Jy Jz S Sy Sz the matrices S G R T and A are constructed using expressions 1 6 1 5 1 4 1 3 and 1 2 The image to image transformation matrix 715 is computed using expression 1 8 The transformation of all points pi in the floating image into points p gt in the reference image using expression 1 7 can be efficiently evaluated using addition operations only Indeed because of the linearity of the transformation the columns of matrix 712 represent the tr
35. n IFR Hp Hg and App are computed from the histogram H using expressions 2 1 2 2 and 2 3 and all other quantities are computed from these MIRIT allows all these quantities to be specified as the registration criterion although only IFR erg UFR and prg are meaningful 12 Chapter 3 Search The optimal set of registration parameters d is defined as k a arg max I d 3 1 with the mutual information of the intensities of corresponding voxels in both images to be registered or any other criterion specified in MIRIT The optimization method that is used for registration is a straightforward implementation of the Powell direction set method as discussed in 3 chapter 10 pp 412 420 For the line minimizations the Brent algorithm is used see 3 chapter 10 pp 402 405 Because the measures pg cpr and Upp need to be maximized while the implementation of the Powell algorithm allows minimization only these measures are internally converted into measures that are minimal when the original measures are maximal The conversions used are log min np ng I 1 e 1 U with np and ng the number of histogram bins used for the floating and the reference image respectively Because 0 lt I lt log min nr ng 0 lt e lt landO U lt 1 these measures obey 0 I lt loga min np ng 0 x 1and 0 U x 1 3 1 The line minimization algorithm The implementation of the lin
36. n acquiring image 1 the gantry angle y can be recovered from Gx S as y atan g The translation matrix T is defined by the 3 translation parameters tz t and t expressed in millimeter representing translation along the x y and z axis respectively 8 lt 1 3 O O O aa ou or oo N The rotation matrix R is the product of 3 matrices Rg Ry and R defined by the 3 rotation parameters Q5 Py and expressed in radians representing rotation around the z y and z axis respectively R Rz Ry Rz 1 0 0 0 R 0 cos Q sin x 0 7 0 sin cos x 0 0 0 0 1 cos Qy 0 sin dy 0 TE 0 1 0 0 i sin dy 0 cos dy 0 0 0 0 1 cos sin 9 0 0 R sin cos 0 0 EG 0 0 1 0 0 0 0 1 such that R is given by denoting cos and sin by cos and sin respectively COS COS COS SIN siny Sinj siny cOos coS sin sinj siny sinz FCcOS COS Siny COSy COS Siny COS SiN SIM cos siny sinz siny COS COS COS 0 0 0 R 1 4 e OOO The skew matrix G is the product of 3 matrices G Gy and G defined by the 3 skew parameters gz gy and gz representing skew in the x y and z direction respectively G Gr Gy G 10g 0 01 0 0 G l0010 00 0 1 1 00 0 g 1 0 0 Gy sqq 0001 1 0 00 0100 E ges E 0 0 0 1 such that G is given by l 9x9z Jx O0 ME E M G 0 os 10 1 5 0 0 0 1 The scale matriz S is defined by the 3 scale parameters sz s an
37. nput The file MIRITinput specifies the image data the user controllable parameters the mode in which the program should operate and the name of the file in which the program writes its output It has the following contents MIRITimages lt filename image data gt MIRITparams lt filename parameters gt 1 lt number Operation mode 1 2 3 or 4 gt MIRIToutput m lt filename output file gt MIRITimages MIRITimages is the name of the file that specifies the floating and reference image data The format of this file is discussed in section 5 2 This file should be accessible from the working directory from which the program is started MIRITparams MIRITparams is the name of the file that specifies all user controllable parameters The format of this file is discussed in section 5 3 This file should be accessible from the working directory from which the program is started MIRIToutput MIRIToutput is the name of the file to which the program writes its output The output format depends on the operation mode in which the program runs For registration the outputfile is in Matlab 4 format If a directory path is specified dirpath filename this should be accessible from the working directory If no directory path is specified the file is written in the working directory Operation mode The operation mode defines the behaviour of the program Currently the following modes are defined e 1 registration computes the opt
38. nsity value at position pg and update the corresponding joint histogram entry while PV interpolation distributes the contribution of this sample over multiple histogram entries defined by its nearest neighbour intensities using the same weights as for TRI interpolation 2 3 Criterion Once the joint image intensity histogram Hg f r has been constructed as described above estimations for the joint and marginal probability density functions of the image intensities in the overlapping volume of images F and R are obtained by normalization of H nr lnmg l N M SO Ha f r f 1 r 1 11 np l r Y p f r f 1 Histogram bins corresponding to value zero in either image are not being taken into account as explained in sections 2 1 and 2 2 Various information theoretic measures are then computed as follows e the marginal entropy of F np l Hy M p f logo p f 2 1 fa the marginal entropy of R nr 1 Hg 3 p r logap r 2 2 r 1 the joint entropy of F and R nrp lnR l Hrr 3 M p f r log2p f r 2 3 f 1 r 1 the conditional entropy of F given R Hygin Hrn Hn the conditional entropy of R given F Hg Hrn Hr the mutual information of F and R Irr Hrp Hr HFR 2 4 the entropy correlation coefficient of F and R TR TEE ae Hr Hg e the uncertainty coefficient of F and R IFR Urr ERE Hp the entropy metric of F and R PFR Hy n Arr Hp
39. on or to decrease the resolution of the floating image to increase speed performance Sampling of the images is specified by the sampling factors fs f fz along each dimension A sampling factor of f for the x dimension means that 1 out of every fy voxels along the image row axis is sampled Sampling is done using trilinear interpolation Each image can be either supersampled f lt 1 or subsampled f gt 1 within the specified ROI Sampling factors of 1 0 in each dimension have to be specified if no sub or supersampling has to be applied 5 3 Parameter specification MIRITparams All user controllable parameters are specified in the file MIRITparams Their interpretation may depend on the operation mode of the program registration tracing histogram or reformatting All parameters have to be specified in all cases although some parameters are not used in some operation modes The contents of the parameter file MIRITparams are 00 0 0 0 0 0 initial position gt 63 40000 parameter order gt O lt 5 lt 10 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 lt nr of steps step sizes gt lt lt 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 12 00 4 256 256 3222332 criterion parameters gt 1 10 0 le 4 1e 3 20 1e 5 20 0 optimization parameters gt Initial position The initial position specifies the initial value for each of the 12 registration parameters in the order Qs Oy bz ta ty tz Sx Sy Sz Qu Jy Gz Q being the rota
40. orithm respectively The initial guess for the bracket is determined by the step size paramaters in the MIRITparams file The bracket values the current guess for the optimum within this bracket and the corresponding function value are displayed at each step Because some criteria needs to be minimized while others have to be maximized and because our implementation of the Powell algorithm only allows minimization those criteria that require maximization are in fact converted to a related measure that is minimal when the true criterion is maximal At the end of each line minimization step the updated position is displayed together with the number of evaluations so far Convergence of each of the Brent steps and of the overall minimization process is determined by the convergence parameters specified in MIRITparams The program terminates by writing the registration result to file which can be interpreted by Matlab 4 provided it has been given the extension m The format of this output file is discussed in section 7 6 2 Tracing mode The images are treated as in registration mode Values of the registration measures are computed along a one dimensional line in the parameter space and written to file which can be interpreted by Matlab 4 provided it has been given the extension m The format of this output file is discussed in section 7 The position corresponding to the first value in the trace is determined from the number of steps and th
41. ples in the overlapping volume of both images the marginal entropy of both the floating and the reference image and their joint entropy and the value of the registration criterion Criterion as specified in MIRITparams From these values all other information theoretic quantities can be computed 32 4 MIRIT Version 97 08 trace Thu Apr 17 18 08 25 1997 REFERENCE IMAGE RefOrientation 1 2 3 RefDimensions 256 256 180 RefSizes 0 976562 0 976562 1 000000 RefGantry 0 000000 RefFile images registration MR img RefHeader 512 RefROI O 255 O 255 O 179 RefSampling 1 000000 1 000000 1 000000 RefRange 1 4094 FLOATING IMAGE FltOrientation 1 2 3 FltDimensions 256 256 100 FltSizes 0 937500 0 937500 1 550000 FltGantry 0 000000 FltFile images registration CT img FltHeader 512 F1tROI O 255 O 255 O 99 FltSampling 1 000000 1 000000 1 000000 FltRange 1 4093 CRITERION Criterion 4 I X Y 4 Interpolation 3 partial volume Criterion 4 RefNrOfBins 256 FloatNrOfBins 256 Interpolation 3 ResolutionLevels 2 2 2 NrOfResolutionLevels 2 ResolutionSteps 3 3 2 TRACE 4 Trace parameters InitialPosition 5 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 1 StepSizes 1 0000000
42. py of the floating image H R F the conditional entropy of the reference image I F R the mutual information of both images U F R the uncertainty coefficient of both images ECC F R the entropy correlation coefficient of both images rho F R the entropy metric of both images other invalid NO oPWNEK i The program exits if a non valid number is specified Histogram size The joint histogram of the overlapping part of both images is computed by bin ning all pairs of floating and reference image intensity values of geometrically corresponding voxels in both images The number of bins to be used for either image is specified by these parameters Typically a value of 256 is used for all images Interpolation The type of interpolation to be used to obtain image intensity values in the reference image is specified by 1 nearest neighbour 2 trilinear 3 partial volume other invalid The program exits if a non valid number is specified 21 Multi resolution parameters These parameters have only effect in registration mode To in crease speed performance MIRIT incorporates a multiresolution strategy for registration see section 3 3 When the optimization starts the criterion is computed from a coarsely sampled floating image while at the end as the optimization convergences the full floating image resolution is considered The number of resolution levels to be considered during optimization and the subsamp
43. s the optimization proceeds A straightforward multiresolution strategy has been incorporated in MIRIT The number of resolution levels and the resolution factor from one level to the next has to be specified for each axis of the floating image The floating image is simply subsampled without prior lowpass filtering in order not to introduce new intensity values in the subsampled images The optimization is started using the subsampled image at the lowest resolution level and the method switches to the next resolution level when convergence is reached at the previous level Experiments have indicated that the robustness of the method is not affected by subsampling the floating image as long as moderate subsampling factors 2 to 4 along each axis are used depending on the resolution of the original floating image Using more than 2 resolution levels a lower resolution one and the full resolution level results generally in slower overall convergence due to the additional evaluations that are needed to reach convergence at the intermediate resolution level When more than 1 resolution level is specified the speed performance can usually be increased by limiting the number of Powell iterations at each level for instance by setting the maximum allowed number of Powell iterations to 5 instead of 20 by default 14 Part II User s manual Chapter 4 Running the program 4 1 Platforms The Mirit software was developed in C on an IBM RS600
44. tations 1 2 The image to world transformation The transformation of homogeneous image coordinates p to homogeneous world coordinates w is defined by the center c of the image in image coordinates its voxelsizes Y vz Vy vz in millimeter and the orientation o of the image axes w r t the world coordinate frame w OxVxC p CT Anterior Axial MR Anterior Sagittal MR Anterior Coronal MR Figure 1 1 Standard orientation of image coordinate frames for different modalities Superior gt Posterior Inferior Superior Posterior Inferior Superior Posterior Inferior Superior Posterior Inferior Origin RAI o 1 2 3 1000 0100 da ERE ME 0 0 0 1 Origin RAS o 1 2 3 10 0 0 01 0 0 m 0 0 1 0 00 0 1 Origin RAS o 2 3 1 0 0 1 0 1 000 Cs 0 1 0 0 0 0 0 1 Origin RPS o 1 3 2 1 0 0 0 0 0 1 0 E 0 1 0 0 0 0 0 1 with C V and O 4 x 4 transformation matrices that represent translation of the origin to the image center C scaling by the voxelsizes Y and permutation of the axes according to 0 respectively 0 0 Cx 1 0 p 0 1 c 0 0 1 S Oe se ooofg E o R o00oo and O O is given by see figure 1 1 ae sign o ifi abs o 1 1 0 otherwise Eventually for CT images a non zero gantry tilt can also be taken into account w OxTxVx C p with I defined by the gantry tilt y me
45. tion parameters t the translation parameters s the scaling factors and g the skew parameters Rotation parameters are specified in degrees translation parameters in millimeter In the example file above the rotation translation and skew parameters are initialized with 0 0 while the scaling factors are set to 1 0 These values are used in each of the operation modes In registration mode the optimization will be started from this initial position In histogram and reformatting mode the images are considered to be at this position In tracing mode the trace is computed with its center point corresponding to this initial position Parameter order The next line specifies the order in which the parameters will be optimized over in the Powell algorithm These values are only used in registration mode In the example above tz is optimized first followed by ty dz 4 and and finally t An entry of 0 means that the corresponding parameter will not be optimied over but will be kept constant and equal to its initial value In the example above the three scaling parameters and the three skewing parameters are unaffected by the optimization and the affine transformation is restricted to a rigid body trans formation only The number of non zero entries determines the number of directions in Powell s direction set method Different non zero entries should start at 1 and be numbered consecutively Registration parameters may be coupled and optimized to
46. tting mode In this mode the reference image is reformatted along the grid of the transformed floating image The reference and floating image are read from file cropped and resampled according to the ROI and sam 24 CRITERION PARAMETERS Criterion 4 I X Y Reference image 256 bins Floating image 256 bins Interpolation 3 partial volume Resolution levels 2 2 2 2 levels Resolution steps 3 3 2 OPTIMISATION PARAMETERS Method 1 Powell Line minimisation parameters max bracket step 1 00e 01 absolute tolerance 1 fractional tolerance max iterations 20 Convergence parameters fractional tolerance max iterations 20 Powell parameters 00e 04 1 00e 03 1 00e 05 order 563124000000 steps 1 00 1 00 1 00 1 00 1 00 1 00 0 10 0 10 0 10 0 10 0 10 0 10 accept new direction 0 REGISTRATION HAS STARTED Initial position 1 0000000000 6 3510195276 6 3510195276 6 3510195276 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 1 0000000000 1 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 Iteration 1 direction 1 0 0000000000 0 0000000000 0 0000000000 1 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 0 0000000000 Bracket 0 0000000000 2 6180340000 2 6180340000 12 3911170327 Brent 2 6180340000 8 6581314112 4 9251459775 8 6581314112 Final position 10
Download Pdf Manuals
Related Search