FW: FW: [Insight-users] MultiResMIRegistration Example
Lydia Ng
lng@insightful.com
Thu, 11 Apr 2002 09:58:54 -0700
Another post from David.
-----Original Message-----
From: David Rich [mailto:David.Rich@trw.com]
Sent: Wednesday, April 10, 2002 11:03 PM
To: Lydia Ng
Cc: luis.ibanez@kitware.com
Subject: RE: FW: [Insight-users] MultiResMIRegistration Example
Lydia,
(Will you post this to the mailing list? My mailer would not recognize
the address. Thanks.)
Okay, you have definitely piqued my interest! I first ran the program
with the image registered to an exact copy. I used 100 sample points,
as before. The results were much better this time:
Target filename: regfile2.raw
Big Endian: 0
Image Size: [128, 128, 105]Image Spacing: 2.72 2.72 2.72=20
Reference filename: CopyOfRegfile2.raw
Big Endian: 0
Image Size: [128, 128, 105]Image Spacing: 2.72 2.72 2.72=20
Number of levels: 5
Target shrink factors: 4 4 4=20
Reference shrink factors: 4 4 4=20
Number of iterations: 250 250 250 250 250=20
Learning rates: 0.0001 5e-005 1e-005 5e-006 1e-006=20
Translation scales: 348 348 348 348 348=20
Number of spatial samples: 100
Standard deviation of Target: 0.4
Standard deviation of Reference: 0.4
Registered image filename: regfile3.dat
Big Endian: 0
Dump PGM files: 1
PGM directory: pgmdir
Reading in target.
Reading in reference.
Normalizing the target.
Mean: 30.5417 StdDev: 63.514
Normalizing the reference.
Mean: 30.5417 StdDev: 63.514
Setting up the registrator.
Start the registration.
Final parameters:=20
0.000285657 0.000150937 8.54086e-005 1 0.0254802 -0.000343909
-0.00488092 =20
Transforming the reference.
Writing registered image to regfile3.dat.
Writing PGM files of the target.
Writing PGM files of the reference.
Writing PGM files of the registered image.
You will notice that I have changed both the input and output to include
both the number of samples and the standard deviations for the target
and reference, which indicate the 100 samples and default standard
deviations.
The results imply almost an exact match.
I then used the reference image with slice 0 missing:
Target filename: regfile2.raw
Big Endian: 0
Image Size: [128, 128, 105]Image Spacing: 2.72 2.72 2.72=20
Reference filename: Regfile2-slice0.raw
Big Endian: 0
Image Size: [128, 128, 104]Image Spacing: 2.72 2.72 2.72=20
Number of levels: 5
Target shrink factors: 4 4 4=20
Reference shrink factors: 4 4 4=20
Number of iterations: 250 250 250 250 250=20
Learning rates: 0.0001 5e-005 1e-005 5e-006 1e-006=20
Translation scales: 348 348 348 348 348=20
Number of spatial samples: 100
Standard deviation of Target: 0.4
Standard deviation of Reference: 0.4
Registered image filename: regfile3.dat
Big Endian: 0
Dump PGM files: 1
PGM directory: pgmdir
Reading in target.
Reading in reference.
Normalizing the target.
Mean: 30.5417 StdDev: 63.514
Normalizing the reference.
Mean: 30.8354 StdDev: 63.7476
Setting up the registrator.
Start the registration.
Final parameters:=20
0.000299443 0.000164776 8.04235e-005 1 0.0265114 -0.00140282
-1.36595 =20
Transforming the reference.
Writing registered image to regfile3.dat.
Writing PGM files of the target.
Writing PGM files of the reference.
Writing PGM files of the registered image.
Again, the registration is excellent. The offset in z is half a pixel,
which I am not sure is quite the right answer, but the process did have
to stretch the reduced image set by one slice.
I then registered the image with the first 5 slices removed, and again
the results were excellent:
Target filename: regfile2.raw
Big Endian: 0
Image Size: [128, 128, 105]Image Spacing: 2.72 2.72 2.72=20
Reference filename: Regfile2-slices04.raw
Big Endian: 0
Image Size: [128, 128, 100]Image Spacing: 2.72 2.72 2.72=20
Number of levels: 5
Target shrink factors: 4 4 4=20
Reference shrink factors: 4 4 4=20
Number of iterations: 250 250 250 250 250=20
Learning rates: 0.0001 5e-005 1e-005 5e-006 1e-006=20
Translation scales: 348 348 348 348 348=20
Number of spatial samples: 100
Standard deviation of Target: 0.4
Standard deviation of Reference: 0.4
Registered image filename: regfile3.dat
Big Endian: 0
Dump PGM files: 1
PGM directory: pgmdir
Reading in target.
Reading in reference.
Normalizing the target.
Mean: 30.5417 StdDev: 63.514
Normalizing the reference.
Mean: 32.0688 StdDev: 64.7051
Setting up the registrator.Start the registration.
Final parameters:=20
0.000315852 0.000151921 6.09405e-005 1 0.0267852 -0.00191591
-6.81043 =20
Transforming the reference.
Writing registered image to regfile3.dat.
Writing PGM files of the target.
Writing PGM files of the reference.
Writing PGM files of the registered image.
Again, the offset in z is about half the thickness of the slices
removed.
During the registration, I observed the memory allocation for the
executable. It typically jumped up and down, indicating some memory
cleanup. However, it still ran between 30MB to 50MB, which raises
questions as to whether or not the intent of using heavily templated
code to improve operation efficiency is being successful. However, in
the meantime, I can now say that the code does work effectively (at
least in this case)--even if the target gets sampled repeatedly.
Now, can you help me understand the differences in the operations
between before and now? What is the specific purpose of the translation
scale? Is it the total amount that the image can be moved? If it is
compared with the rotation parameter of 1.0 for complete rotation, is it
necessary to use the full width of the image to get comparable offsets?
And how are the learning rates used in conjunction with the translation
and rotation paameters? Anything else you can tell me to help me
understand the proper application here would be greatly appreciated.=20
Dave
>>> "Lydia Ng" <lng@insightful.com> 04/10/02 06:32PM >>>
David,
> I hope I haven't bored you with all this stuff. If I have=20
> missed something critical, I would be glad to try again. =20
> Please let me know if you have other questions.
Would you try an experiment for me?
Can you up the *translation scale* so that it is
348 348 348 348
And decrease the *learning rate* to something like
1e-4 5e-5 1e-5 5e-6
Change the shrink factors for both image to
4 4 4
and use the default standard deviation, number of samples etc.
Could first try this on registering the same image
and then registering when you have shifted one of the images?
I am just curious as to what results you get?
- Lydia
> -----Original Message-----
> From: David Rich [mailto:David.Rich@trw.com]=20
> Sent: Wednesday, April 10, 2002 4:59 PM
> To: Lydia Ng
> Cc: luis.ibanez@kitware.com; insight-users@public.kitware.com=20
> Subject: RE: FW: [Insight-users] MultiResMIRegistration Example
>=20
>=20
> Lydia,
>=20
> Here is the output file from a run with 100 data points (I=20
> usually redirect the output to a file so I can capture it for=20
> the record). Since everything is essentially here, I have=20
> not included the imput file:
>=20
> Target filename: regfile2.raw
> Big Endian: 0
> Image Size: [128, 128, 105]Image Spacing: 2.72 2.72 2.72=20
>=20
> Reference filename: CopyOfRegfile2.raw
> Big Endian: 0
> Image Size: [128, 128, 105]Image Spacing: 2.72 2.72 2.72=20
>=20
> Number of levels: 5
> Target shrink factors: 4 4 1=20
> Reference shrink factors: 4 4 1=20
> Number of iterations: 250 250 250 250 250=20
> Learning rates: 1 0.1 0.01 0.0001 1e-008=20
> Translation scales: 7 3.5 1.75 0.875 0.438=20
>=20
> Registered image filename: regfile3.dat
> Big Endian: 0
>=20
> Dump PGM files: 1
> PGM directory: pgmdir
>=20
> Reading in target.
> Reading in reference.
> Normalizing the target.
> Mean: 30.5417 StdDev: 63.514
> Normalizing the reference.
> Mean: 30.5417 StdDev: 63.514
> Setting up the registrator.
> Start the registration.
> Final parameters: -0.889207 -0.434599 -0.0411224 0.136909 =20
> -2.05609 6.46628 -8.74688 =20
> =20
> The "CopyOfRegfile2" is the exact copy of "Regfile2." As you=20
> can see, the rotations and translations all look pretty bad. =20
> I could include shots of some of the frames, but the numbers=20
> should be self explanatory.
>=20
> I also tried shifting the reference file by removing the=20
> first slice and again by removing the first five slices. The=20
> data of interest does not start until slice 27, which is what=20
> generated the question Luis responded to about limiting the region.
>=20
> With a single slice removed and using 1000 sample points, the=20
> input and results were good:
>=20
> regfile2.raw
> 0
> 128 128 105
> 2.72 2.72 2.72
> regfile2-slice0.raw
> 0
> 128 128 104
> 2.72 2.72 2.72
> 5
> 4 4 1
> 4 4 1
> 250 250 250 250 250
> 1.0 0.1 0.01 0.0001 1.0e-8
> 7.0 3.5 1.75 0.875 0.438
> regfile3.dat
> 0
> 1
> pgmdir
>=20
> Ran for over an hour, fitting 1000 data points. Results:
>=20
> 0.000822807 0.000224679 2.67971e-5 1 0.143948 =20
> 0.239717 -2.32141
>=20
> However, when I reduced the sample points to 100, the results=20
> were not so good:
>=20
> Mean: 30.5417 StdDev: 63.514
> Normalizing the reference.
> Mean: 30.8354 StdDev: 63.7476
> Setting up the registrator.
> Start the registration.
> Final parameters:=20
> 0.174784 0.281457 -0.247609 -0.910452 -2.86575 3.51592 =20
> -3.99107 =20
>=20
>=20
> With the first 5 slices removed and using 100 sample points,=20
> the results were again pretty poor:
>=20
>=20
> Final parameters:=20
> 0.432085 0.85646 0.266564 -0.0933978 1.0369 3.02536 -5.7002 =20
>=20
>=20
> I changed the translaction scale to pixel size (remembering=20
> that only the first is used), but the results did not=20
> improve, just changed:
>=20
>=20
> Translation scales: 2.72 2.72 1.36 1.36 0.68=20
> Results:
> Mean: 30.5417 StdDev: 63.514
> Normalizing the reference.
> Mean: 32.0688 StdDev: 64.7051
> Setting up the registrator.
> Start the registration.
> Final parameters:=20
> -0.0625984 0.935416 0.346148 0.0354966 0.902614 0.302101 0.252
>=20
> I increased the translation scale to 40 without much success:
>=20
> Final parameters:=20
> -0.0675268 0.164706 -0.661179 0.728803 -51.6737 -68.8085=20
> 224.462 =20
>=20
> Changing the translation scale back to 5.5, about 2 pixel=20
> widths, and increasing the samples to 1000, using the image=20
> with the 5 slices removed, still gave poor results:
>=20
> Final parameters:=20
> 0.594023 -0.0829148 -0.798207 0.0559267 -0.404611 1.3842=20
> -6.94312 =20
>=20
> I hope I haven't bored you with all this stuff. If I have=20
> missed something critical, I would be glad to try again. =20
> Please let me know if you have other questions.
>=20
> Dave
>=20
>=20
>=20
> >>> "Lydia Ng" <lng@insightful.com> 04/10/02 05:05PM >>>
> Luis, David:
>=20
> In terms of memory - I think the multi-resolution image pyramid
> is indeed the culprit.=20
>=20
> Historically, it only computed the images at each level one at
> at time. That was subsequently changed to computing all the level
> following comments of the code review. To support large images I think
> this will need to change back to outputing one level at a time.
>=20
> -however-
>=20
> What I am more concern about is that David is reporting
> that the underlying algorithm is not working.
>=20
> Luis: you are correct when you say that when dealing with
> intra-modality one might be better off using something other
> than MI. But both you and I have seen working demos of when the
> same MI code can successfully register MR to MR with=20
> significant translation and rotation.
>=20
> David:
> [1] How quantitatively bad were the results=20
> when you only used 50 sample points?
> [2] Have you tried tests where you have shifted and/or rotated
> the image by a known amount? What happened?
> [3] Would you mind posting your parameter file?
>=20
> - Lydia=20
>=20
> > -----Original Message-----
> > From: Luis Ibanez [mailto:luis.ibanez@kitware.com]=20
> > Sent: Wednesday, April 10, 2002 3:35 PM
> > To: David Rich
> > Cc: Lydia Ng; insight-users@public.kitware.com=20
> > Subject: Re: FW: [Insight-users] MultiResMIRegistration Example
> >=20
> >=20
> >=20
> > Hi David,
> >=20
> > Thanks for your feedback on you experience with
> > the toolkit.
> >=20
> > ITK is not intended only for education and research
> > purposes. It has been designed to be used for
> > commercial products as well. Particlar attention
> > has been paid to the terms of the license so that
> > commercial use becomes a reality.
> >=20
> > Performance both in memory and computing time are
> > important elements in ITK design. The heavy use of
> > templated code, for example, responds to the need of
> > simplifying the maintenance of the code and improving
> > its performance. Templates allow a lot of code to be
> > inlined and also substitute the run-time mechanism
> > of virtual function overloading with more expedite
> > compile-time instantiations.
> >=20
> > However the toolkit needs to pass through the natural
> > process of code maturation. In this sense, feedback
> > from users is of fundamental importance and is always
> > greately appreciated.
> >=20
> > Your report on memory consumption is quite interesting.
> > We certainly want to track down the reasons for the
> > excess of memory use and come up with a solution for it.
> >=20
> > The fact that you observed the memory growing as the
> > registration ran, leads to suspect that the image pyramids
> > may be implicated. At each stage of the registration
> > only images at the same level of the pyramid are required
> > to be in memory. It is likely that the current algorithm
> > is not releasing the memory used by previous levels of
> > the pyramid. Also, the subsampling is performed by first
> > smoothing the images with Gaussian Blurring filters.
> > These filters allocate memory for internal copies of the
> > images for intermediate computations. It is also likely
> > that this memory is failing to be released. If this is
> > the case it should be relatively easy to fix. We will
> > explore these possibilities.
> >=20
> >=20
> > You also mention an option for increasing performance
> > by computing the metric over a restricted region of
> > interest in the image. This modification has recently
> > being made on the metrics. In order to take advantage
> > of it you just have to provide the metric with a region
> > on the Fixed image. The rest of the Fixed image will
> > be excluded from the computations required for registration.
> >=20
> > The use is like:
> >=20
> > metric->SetFixedImageRegion( region );
> >=20
> >=20
> >=20
> >=20
> > A comment about Mutual Information: this metric is not
> > particularly well suited for images of the same modality,
> > and specially not for registering an image with itself.
> > The reason is that when you produce the joint histogram
> > of an image with itself (assuming they are already
> > perfectly registered) it will result in filling only the
> > diagonal elements of the joint histogram. In consequence
> > only N out of NxN bins will be contributing information to
> > the matching. Mutual Information behaves better when it is
> > used with Multimodal images because in that case the joint
> > histogram have a better distribution of samples (larger
> > joint entropy). The registration will be naturally slow when
> > the images get closer to the optimal position because in
> > that case only the regions of the image with large gradients
> > will contribute to variations in the metric (and hence to the
> > computation of the metric derivatives that drive the gradient
> > descent optimization method). The chances of random points
> > to fall into gradient areas are pretty small in general, so
> > most of the points will be only contributing to the diagonal
> > of the joint histogram. That is, most of them are just
> > reporting that the images seems to be registrered. That may
> > explain why your test required a larger-than-normal population
> > of points to progress and why it fails to register after
> > getting to 1 or 2 pixels of distance from the optimal position.
> > (which seems reasonable to expect for an image with itself).
> > Multimodal images have better chances of producing joint
> > entropy gradients with smaller populations of sampling points.
> > Same modality images will register much better by using a
> > MeanSquare metric for low resolutions and a PatternIntensity
> > metric for high resolution. [I understand that you did this
> > only for the purpouse of testing the code so this is just a
> > comment for the record]
> >=20
> >=20
> > We'll profile the memory consumption of the method and get
> > back to you as soon as we identify the sources.
> >=20
> >=20
> > Thanks a lot for your feedback
> >=20
> >=20
> >=20
> > Luis
> >=20
> >=20
> > =
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
> > >>
> > >>While watching this operate for over an hour, I observed that=20
> > >>the executable started out at about 20MB but increased to=20
> > >>almost 60MB before completion! Such excessive memory usage=20
> > >>becomes intolerable in real-world applications-what if I were=20
> > >>trying to register two 40MB images?=20
> > >>
> > >>So, my first question is whether or not ITK is intended to be=20
> > >>used in commercial applications or is it only designed as a=20
> > >>tool for people who want to experiment with many different=20
> > >>models? The templated classes would certainly facilitate the=20
> > >>latter but the complexity combined with size and time=20
> > >>constraints do not contribute to the former.
> > >>
> > >>In an attempt to better understand the parameters that could=20
> > >>be entered, I continued testing. The initial results in=20
> > >>about 23 minutes of operation were:
> > >>
> > >>Final parameters:=20
> > >>-0.000313974 -0.000876174 -0.000350401 1 -0.431262 =20
> > >>0.368347 -0.0612012 =20
> > >>
> > >>The rotation parameters indicate a close fit although the=20
> > >>translations are disconcerting, considering that an exact=20
> > >>match should be more feasible (Note: the voxel size is 2.72=20
> > >>in all directions).
> > >>
> > >>I changed the default input of 0.4 standard deviation for the=20
> > >>registration to 0.2 and used 100 sample points, and the=20
> > >>results were miserable:
> > >>
> > >>Final parameters:=20
> > >>0.574853 -0.761886 -0.298448 -0.00139123 -19.8732 =20
> > >>28.4372 -79.6355 =20
> > >>
> > >>With a standard deviation of 0.6 and 100 sample points the=20
> > >>results are different but still miserable:
> > >>
> > >>Final parameters:=20
> > >>-0.0671739 0.0841046 -0.994183 0.00385258 -14.183 =20
> > >>-1.85797 -1.66959 =20
> > >>
> > >>The conclusion seems to be that with enough sample points the=20
> > >>basic algorithm may provide satisfactory results. However,=20
> > >>from past experience the time delay becomes critical for=20
> > >>medical personnel. And if the algorithm requires a minimum=20
> > >>of 500 data points, 60MB RAM, and 23 minutes on a 1.6MB image=20
> > >>registered with itself, what would be required for two=20
> > >>less-similar images of larger sizes?
> > >>
> > >>An examination of the source code to try to understand the=20
> > >>parameters and whether or not the memory can be handled more=20
> > >>efficiently again reminded me of the Microsoft ATL wizard. =20
> > >>Only this time it seemed like it was necessary to understand=20
> > >>the complexities of the ATL wizard for the purpose of=20
> > >>creating a specific implementation. And again it occurred to=20
> > >>me that the purpose of ITK seems to be that of creating a=20
> > >>framework for research and experimentation for those who are=20
> > >>not concerned with commercial requirements of speed and=20
> > >>hardware constraints. Am I trying to use this contrary to=20
> > >>the intent of the toolkit?
> > >>
> > >>On the other hand, is it possible that the ITK could be=20
> > >>developed into something more like the ATL wizard? That is,=20
> > >>ITK with the power of the templates built in, could be used=20
> > >>to generate a basic code structure for the specific=20
> > >>implementation requested. With such a code framework=20
> > >>generated, it might be more feasible for the researcher as=20
> > >>well as the commercial user to customize the specific=20
> > >>implementation for speed and memory or to modify the=20
> > >>mathematical calculations to fit specific requirements. =20
> > >>
> > >>At this point I feel little more than frustrated. Although I=20
> > >>would like to be able to provide commercial support to public=20
> > >>code, customizing it for customers but being able to leave=20
> > >>the source with them, I can only recommend that ITK is not=20
> > >>the right framework for commercial applications. It is slow,=20
> > >>cumbersome, and requires significant hardware resources. If=20
> > >>anyone wants to delve into the complexities of the code, they=20
> > >>have to peel away multiple layers of typedefs and templated=20
> > >>classes to try to evaluate implemented constructs and=20
> > >>operations. If the end result in operation, size, and=20
> > >>effectiveness were tremendous, the complexities could be=20
> > >>overlooked. That does not seem to be the case.
> > >>
> > >>I would like to report more positive feedback. Can you help me? =20
> > >>1) For the MultResMIRegistration example, how can I identify=20
> > >>parameters that will be stable for generic image sets? =20
> > >>2) How can I improve the time and memory requirements?
> > >>3) Can I tell the algorithm to skip the blank slices at the=20
> > >>front and/or back of the data of interest?
> > >>4) Or, is this toolkit just not intended to be used in=20
> > >>software used by medical personnel and others whose time is=20
> > critical?
> > >>
> > >>Any suggestions that might help me understand the toolkit=20
> > >>better or how to make it effective would be greatly appreciated.
> > >>
> > >>Dave Rich
> > >>
> > >>
> > >>
> >=20
> >=20
> >=20
>=20
_______________________________________________
Insight-users mailing list
Insight-users@public.kitware.com=20
http://public.kitware.com/mailman/listinfo/insight-users