Section 01:
Ok, probably more info than you want to know, but here is the full scoop:
I work as a medical physicist and as part of the job we have to verify that the radiation that we are giving to the patient is what we thought it was. How do we do this?
We plan the patient on a 3D planning system based off of a CT image set. We use the CT numbers (Hounsfield numbers) to determine the electron density in any given region, and then can determine how the radiation will be affected (attenuated) and converted into dose. When we have a plan that we like and the physician is happy with, we then take all of the parameters from the plan (beam angles, shapes…) and deliver that to a phantom. { a phantom is a device for which the makeup and geometry is well know, often it is made of solid water (a plastic with radiological properties within the range of treatment energies, 4MV ~ 20MV photons and 6MeV~22MeV electrons, nearly identical to that of water)}.
The phantom also has a CT study and we can get a 3D plan similar to the patient. In my case the phantom comes in 2 pieces, a 5cm thick piece on top, and a 5cm thick piece on the bottom. The top piece has 7 holes in the same orientation as shown in my post about the film. We put a piece of XV film between the 2 pieces of the phantom, mark the film with punctures at the holes and deliver the treatment. We can now compare the dose delivered to the film to the calculated does at the same plane as the film. To do this we have to register the film with the calculated plan. We will always align the phantom to the treatment machine, but the film has ~ .5-1cm of play, so we do not have consistency between films, so we can not assume that the points are always at the same location, but we can be fairly certain that they are in the same region. We can process up to 50 of these films a day (currently this number continues to increase).
The minimum number of points to have a good registration is 4.
What I have done so far is to write a macro for imageJ to find points that meet a certain criteria within the region of interest, and then write to a file the matching points.
ImageJ out puts it in the following format:
1 X Y
2 X Y
3 X Y
.
.
.
N X Y
1 X Y
2 X Y
.
.
M X Y
1 -1 -1
1 X Y
.
.
The first column describes the points, every time there is a new “1” it is in a new ROI. I will get you the ROI order with the actual data. The second and third column is self explanatory. If there is a “-1” in it, it means that no points were contained in that ROI.
I wrote a very rough/simple script to import the data (it is ugly, I was only going for expedience in writing, I always clean things up after they work…and I am not good enough to write good code first time through) Here it is if you want it:
my $i = 0;
my $toggle = -1;
my @p = qw|ROI1 ROI2 ROI3 ROI4 ROI5 ROI6 ROI7|;
my @points,%ROI;
open (POINTS, "<points.txt");
while (<POINTS>) {
chomp $_;
my $a,$b,$c;
($a,$b,$c) = split(/\t/,$_);
if ($a == 1) {
$toggle *= -1;
if ($toggle < 0) {
$i++;
$toggle = 1;
}
}
if ($b == -1) {
$i++;
$toggle = -1;
} else {
push @{$p[$i]}, [$b,$c];
}
}
Also I wrote a script to filter the points and remove any points that can not be valid (i.e. less then 4 adequate distances). This I did not test, do not have actual data yet, will have the data today or tomorrow.
# this is sent the roi info and the distances from that point to all o
+f the other points from the template
#$index pulls in the roi info, @deltas then pulls in the distances.
#if there is a match from a point in the region in question and the re
+gion it is testing against, it stops testing against that region and
+increments its success, if it has 4 successes, it stops checking that
+ point. If at the end it does not have 4 successes it splices that p
+oint….. again untested.
sub filter1
{
my $i = $j = $k = 0;
my $index = shift;
my @deltas = @_;
foreach my $testreg (@{$p[$index]}) {
my ($testX,$testY) = @$testreg;
foreach my $region (@p) {
do {
if ($i == $index) {$i++}
else {
$k = 0;
foreach my $points(@{$region}) {
do {
my ($X,$Y) = @$points;
my $distance = sqrt(abs(($testX - $X)**2 +
+ ($testY - $Y)**2));
if (($deltas[$i] - $precision) <= $distanc
+e) && ($deltas[$i] + $precision) >= $distance)){
$k++;
$j++;
}
} unless ($k == 1);
}
}
} unless ($j == 4);
}
if ($j < 4) {
splice (@{$p[$index]},$i,1);
} else {
$i++;
}
}
}
I think that after this first filter, the total number of points should be much more manageable.
Let me know if you have any questions.
Here is the data:
1 599 377
2 829 410
3 643 465
4 522 511
5 626 552
6 762 575
7 577 597
8 699 618
1 1474 13
2 1507 14
3 1882 15
4 1912 15
5 1704 16
6 1573 17
7 1685 21
8 1709 21
9 1795 90
10 1668 107
11 1732 131
12 1819 144
13 1755 159
14 1692 163
15 1936 213
16 1708 218
17 1808 227
18 1704 221
19 1747 229
20 1744 293
1 2673 619
2 2765 620
3 2839 638
4 2729 659
5 2779 688
6 2681 703
7 2744 729
8 2821 738
9 2807 788
10 2740 808
11 2560 822
12 2570 821
1 3308 1800
2 3371 1828
3 3280 1854
4 3412 1856
5 3274 1934
6 3360 1986
7 3308 1998
8 3366 2046
1 2235 2522
2 2408 2541
3 2335 2547
4 2267 2580
5 2383 2602
6 2296 2630
7 2214 2627
8 2426 2652
9 2425 2657
10 2228 2678
11 2364 2684
12 2283 2726
13 2358 2764
1 1768 2505
2 1696 2527
3 1800 2567
4 1650 2586
5 1733 2624
6 1835 2632
7 1621 2675
8 1689 2677
9 1650 2730
10 1770 2738
1 29 1746
2 30 1764
3 78 1781
4 30 1783
5 30 1806
6 320 1897
7 231 1944
8 31 1938
9 30 1968
10 376 1965
11 206 2013
12 260 2068
13 193 2096
14 401 2097
15 350 2133
16 32 2166
17 257 2187
18 32 2258
____________________
1 585 600
1 1715 280
1 2745 815
1 3240 2090
1 2280 2770
1 1715 2770
1 190 2090
The first set are from a film that I just created with a lot of noise, and 2 regions do not contain any true points. The second set are from the template, and the order for the roi's are the same and follow:
----------------------------------
| |
| 2 |
| |
| 1 |
| 3 |
| |
| |
| 7 4 |
| |
| |
| 6 5 |
| |
----------------------------------
Film has a lattace of silver (a molecule containing silver) When the film is exposed to light/x-rays (photons), the silver which has a high probability of interaction becomes ionized. Then when the fixer is used the ionized portion remains and the rest is washed away. This is why in an xray, bone looks white. Compared to other tissue the photons are heavily attenuated by the bone, so most of the silver is washed away. When you damage film, i.e. punture or even just by pressure, it damages the lattace, and when the film is developed, even if it wasn't exposed, the damage can be seen.
As far as the template goes...Think of 2 blocks of styrofoam with dimensions 25cm x 25cm x 5cm, if you drilled holes in one block such that went completely though the 5cm thickness, and the holes were grouped in a known pattern, then when you placed the film between the 2 blocks, you can gain access to the film through the holes. So if you cause damage to the film at the locatoion of the holes, then you will always be able to maintain the alignment of the film to the blocks at the time of the exposure.
Does this help?
Ok...You are very close on many things. I think that you have the problem figure and understood. The geometric consistency IS there, but it is not a perfect 1:1. The jig is fine. The "plunger fits very snug into the holes, BUT there is no way to determine how the damage to the film will be represented (shape) when the film is processed/scanned. I have seen the full range of:
______ __________
/ \ / -- /
| -- | to | -- |
| -- | | /
\_______/ \____/
where the outside ring represents the darkening of the film due to the damage, and the in side marks represent the actual puncture. The resolution that we scan at is 89 microns per pixel, and the average darking per puncture is ~ 3mm across => 33 pixels. The digital processing that I do on the film after scanning finds the points, then it asigns the X,Y value to that point by way of center of mass. So in the first exampl above we would have a near perfect match, but we could be off by 10-15 pixels in the second example.
After getting some prelim. results I checked them against the film. here are the real points +/- 2 pixels:
ROI X Y
1 464 642
2 157 1756
3 687 2779
4
5 2631 2296
6 2625 1732
7
Notice that the X/Y are flipped from the input data, that is because the data given after points are found versus the software where the analysis is done call the columns/rows differently. So reverse the order to find the closest pair which would be:
ROI 1
3 643 465
ROI 2
13 1755 159
ROI 3
5 2779 688
ROI 4
None
ROI 5
6 2296 2630
ROI 6
5 1733 2624
ROI 7
None
If we look at distance beween points:
True Point Distances for ROI1.
ROI2
Distance = 1153.33429672407
ROI3
Distance = 2147.60913576004
ROI5
Distance = 2723.90051213329
ROI6
Distance = 2418.54935860321
Template Point Distances for ROI1.
ROI2
Distance = 1174.4360348695
ROI3
Distance = 2170.67385850569
ROI5
Distance = 2753.52955313721
ROI6
Distance = 2446.58946290545
True Point Distances for ROI2.
ROI1
Distance = 1153.33429672407
ROI3
Distance = 1152.56973758641
ROI5
Distance = 2529.5299958688
ROI6
Distance = 2465.09817248725
Template Point Distances for ROI2.
ROI1
Distance = 1174.4360348695
ROI3
Distance = 1160.65714145048
ROI5
Distance = 2553.29688833868
ROI6
Distance = 2490
True Point Distances for ROI3.
ROI1
Distance = 2147.60913576004
ROI2
Distance = 1152.56973758641
ROI5
Distance = 2001.16291190897
ROI6
Distance = 2200.50266984614
Template Point Distances for ROI3.
ROI1
Distance = 2170.67385850569
ROI2
Distance = 1160.65714145048
ROI5
Distance = 2009.53974830059
ROI6
Distance = 2209.73414690546
True Point Distances for ROI5.
ROI1
Distance = 2723.90051213329
ROI2
Distance = 2529.5299958688
ROI3
Distance = 2001.16291190897
ROI6
Distance = 563.031970673069
Template Point Distances for ROI5.
ROI1
Distance = 2753.52955313721
ROI2
Distance = 2553.29688833868
ROI3
Distance = 2009.53974830059
ROI6
Distance = 565
True Point Distances for ROI6.
ROI1
Distance = 2418.54935860321
ROI2
Distance = 2465.09817248725
ROI3
Distance = 2200.50266984614
ROI5
Distance = 563.031970673069
Template Point Distances for ROI6.
ROI1
Distance = 2446.58946290545
ROI2
Distance = 2490
ROI3
Distance = 2209.73414690546
ROI5
Distance = 565
We can see fairly good agreement, but it does appear to be +/- 30
Well I am going to work further, and will start to implement some of your code for testing....
