[Insight-users] ITK Image Coordinates / Problem withPhysicalPoint to Index Conversion

Sat Feb 2 06:37:17 EST 2008

Hi Andreas,

it seems to me that the problems you indicate are related to the question: what is the co-ordinate of a pixel?
The way I see it, a pixel has no co-ordinate (but it has an area).
I could rephrase your step1 one as follows:

1. For 
the 
given 
output 
index, calculate the co-ordinate >at the center< of the corresponding output pixel.

Ideally, this should be performed by a helper function: IndexToPhysicalCoordinateAtPixelCenter(). Note that the output of this function depends on whether you have the origin of the image at the bottom left or center of a pixel. For a 'bottom-left' representation, the function IndexToPhysicalCoordinateAtPixelCenter(0) would return 0.5 * imageSpacing, whereas for a 'center' representation, it would be zero.
However, the writer of the resample algorithm will not be concerned with this, because in both cases  IndexToPhysicalCoordinateAtPixelCenter() will return the co-ordinate at the center of the pixel, and the algorithm will perform correctly.

Best regards, and looking forward to your replies too,
Maarten

----- Original Message ----
From: Andreas Keil <andreas.keil at cs.tum.edu>
To: insight-users at itk.org
Sent: Tuesday, January 29, 2008 6:16:31 PM
Subject: RE: [Insight-users] ITK Image Coordinates / Problem withPhysicalPoint to Index Conversion

Dear 
all,

let 
me 
elaborate 
why 
I 
am 
voting 
for 
a 
center-based 
coordinate 
system:

Let's 
consider 
the 
simple 
example 
of 
resampling 
an 
image. 
The 
resample
filter 
does 
the 
following 
for 
every 
output 
pixel:
1) 
Calculate 
the 
physical 
coordinate 
for 
the 
given 
output 
index.
2) 
Transform 
the 
physical 
coordinate 
by 
the 
given 
transform.
3) 
Calculate 
the 
continuous 
index 
for 
the 
transformed 
physical 
point.
4) 
Linearly 
interpolate 
the 
input 
value 
at 
the 
continuous 
index 
and

write 
it 
to 
the 
output.

Now, 
think 
about 
this 
process 
with 
corner-based 
coordinated:
1) 
Would 
work 
fine.
2) 
Transforms 
are 
there 
for 
converting 
physical 
coordinates. 
Giving 
them

corner 
coordinates 
would, 
of 
course, 
work, 
too. 
However, 
the

interpolation 
in 
step 
4 
is 
done 
based 
on 
a 
single 
point 
representing

the 
whole 
pixel. 
In 
general, 
we 
don't 
have 
an 
idea 
on 
how 
the 
box 
of

a 
pixel 
is 
transformed. 
We 
can 
just 
pass 
a 
single 
coordinate 
through 

through 
to 
the 
interpolator. 
Keep 
this 
in 
mind 
for 
now.
3) 
Works 
fine.
4) 
Now 
we 
have 
a 
problem 
with 
a 
corner-based 
physical 
point 
having 
been

transformed: 
The 
interpolator 
has 
to 
decide, 
which 
closest 
pixels 
to

take 
into 
account 
for 
interpolation 
(and 
how 
to 
weight 
them). 
This

would 
be 
a 
problem 
even 
for 
a 
the 
linear 
interpolator. 
How 
should

it 
know, 
which 
4 
pixels 
(in 
the 
2D 
case) 
to 
take? 
The 
problem 
is,

that 
the 
interpolator 
has 
to 
estimate 
the 
location 
of 
the 
backward

warped 
input 
pixel. 
This 
is 
never 
really 
possible 
due 
to 
the 
unknown

transformation. 
But 
the 
guess 
of 
a 
center-based 
rectangle 
is 
better

than 
guessing 
a 
corner-based 
rectangle 
(where 
you 
don't 
know, 
which

of 
the 
corners 
is 
given).

A 
numerical 
example 
(sorry, 
I 
was 
too 
lazy 
to 
paint): 
Let's 
resample

a 
2D 
image 
with 
an 
input 
and 
output 
spacing 
of 
1mm 
x 
1mm. 
We 
want 
to

interpolate 
the 
output 
pixel 
at 
(0mm..1mm,0mm..1mm) 
which 
has 
the

corner 
coordinates 
(0mm,0mm).

A) 
For 
a 
translation 
of 
(0.1mm,0.2mm), 
we 
get 
an 
input 
coordinate

at 
(0.1mm,0.2mm) 
and 
now 
have 
to 
decide, 
which 
image 
range 
of

2x2 
pixels 
to 
use 
for 
interpolation. 
Since 
we 
translated 
a 
bit

in 
positive 
x 
and 
y 
direction, 
we 
should 
use 
the 
pixels 
at

(0mm..2mm,0mm..2mm). 
Alright. 
So 
the 
rule 
should 
be 
to 
round 
the

physical 
coordinate 
and 
this 
yields 
the 
lower 
left 
corner 
of 
our

2x2 
pixels 
interpolation 
region.

B) 
Let's 
try 
a 
different 
transformation: 
A 
180 
deg. 
rotation 
around

(0mm,0mm) 
and 
a 
subsequent 
translation 
of 
(0.1mm,0.2mm). 
We 
again

compute 
(0.1mm,0.2mm) 
as 
input 
coordinate. 
So 
let's 
follow 
the

rule 
and 
use 
again 
(0mm..2mm,0mm..2mm) 
as 
interpolation 
region.

Ups! 
We 
should 
have 
used 
(-1mm..1mm,-1mm..1mm) 
instead! 
(You 
can

see 
this 
by 
also 
transforming 
the 
"upper-right" 
corner 
output

pixel 
which 
is 
(1mm,1mm) 
giving 
(-0.9mm,-0.8mm) 
as 
input 
corner.

Therfore, 
the 
output 
pixel 
lies 
at 
(-0.9mm..0.1mm,-0.8mm..0.2mm).

This 
input 
pixel 
overlaps 
the 
4 
input 
pixels 
in 
the 
region

(-1mm..1mm,-1mm..1mm) 
as 
stated 
above.)

Conclusion: 
The 
transformed 
corner 
coordinate 
in 
reality 
is 
no

longer 
a 
"lower-left" 
corner 
but 
how 
should 
the 
interpolater 
have

guessed 
that? 
The 
only 
remedy 
would 
be 
that 
every 
class 
using 
a

transformation 
would 
not 
only 
transform 
one 
single 
coordinate 
per

pixel 
but 
all 
4 
(in 
2D) 
or 
8 
(in 
3D) 
corners. 
Nobody 
wants 
this.

Guessing 
an 
input 
region 
after 
having 
backward 
warped 
a 
pixel's

center 
is 
much 
easier 
and 
faster. 
(Although 
I 
would 
opt 
for 
taking

circles 
instead 
of 
rectangles 
here, 
but 
that's 
another 
issue.)

Hopefully, 
I 
made 
my 
point 
clear 
that 
a 
pixel's 
center 
coordinate 
is 
the
more 
useful 
representation 
for 
it 
than 
its 
corner 
coordinate. 
I 
am 
pretty
sure, 
the 
original 
developers 
or 
writer 
of 
the 
docs 
had 
this 
point 
in 
mind
when 
the 
wrote 
the 
Software 
Guide. 
And 
since 
the 
implementation 
is 
flawed
anyways, 
and 
since 
a 
lot 
of 
classes 
could 
not 
be 
adopted 
to 
corner-based
coordinates, 
I 
still 
strongly 
vote 
for 
center-based 
ones. 
This 
"physical
point 
of 
view" 
is 
also 
more 
in 
the 
style 
of 
ITK 
than 
the 
"computer 
science
point 
of 
view". 
ITK 
developers 
should 
not 
think 
in 
pixels 
but 
rather 
in
millimeters.

Thanks 
for 
reading. 
Waiting 
for 
your 
replies,
Andreas.

-----Original 
Message-----
From: 
Maarten 
Nieber
Sent: 
Monday, 
January 
28, 
2008 
14:45
To: 
insight-users at itk.org
Subject: 
Re: 
[Insight-users] 
ITK 
Image 
Coordinates 
/ 
Problem
withPhysicalPoint 
to 
Index 
Conversion

Hi 
Andreas,

I 
agree 
with 
you 
that 
this 
issue 
is 
of 
very 
high 
importance!

-- 
quote 
from 
your 
previous 
post 
--
I 
think 
that 
the 
personal 
preference
for 
one 
of 
the 
two 
options 
we 
have 
basically 
depends 
on 
whether 
one 
is
more 
interested 
in 
the 
overall 
dimensions 
of 
a 
dataset 
or 
whether 
one 
has
to 
do 
voxel-wise 
operations 
depending 
on 
the 
physical 
coordinate.
-- 
end 
quote 
--

Can 
you 
explain 
why 
you 
think 
- 
for 
doing 
voxel-wise 
operations 
- 
it 
would
be 
better 
if 
the 
image 
origin 
co-incides 
with 
the 
center 
of 
a 
pixel?

Best 
regards,
Maarten

----- 
Original 
Message 
----
From: 
Andreas 
Keil
To: 
insight-users at itk.org
Sent: 
Tuesday, 
January 
8, 
2008 
10:37:32 
AM
Subject: 
RE: 
[Insight-users] 
ITK 
Image 
Coordinates 
/ 
Problem 
with
PhysicalPoint 
to 
Index 
Conversion

Dear 
all,

I 
would 
like 
to 
try 
a 
restart 
of 
the 
discussion 
related 
to 
the 
definition
of 
physical 
coordinates 
if 
ITK 
images:

As 
mentioned 
in 
my 
earlier 
post 
(see 
below), 
the 
documentation 
and
implementation 
differ 
in 
this 
point. 
The 
two 
proposed 
solutions 
are 
either
fixing 
the 
documentation 
(ITK 
Software 
Guide, 
p. 
40) 
and 
defining 
the
origin 
to 
lie 
in 
the 
"lower-left(-back)" 
corner 
of 
the 
pixel 
with 
index
0/0(/0) 
or 
fixing 
the 
implementation 
of 
the 
respective 
methods
(IndexToPhysicalPoint, 
PhysicalPointToIndex, 
etc.) 
and 
defining 
the 
origin
to 
lie 
in 
the 
center 
of 
this 
pixel.

So 
far, 
only 
Maarten 
and 
I 
were 
discussing 
about 
this, 
with 
him 
favoring
to 
fix 
the 
docs 
and 
me 
favoring 
to 
fix 
the 
implementation. 
Another
argument 
I 
found 
for 
my 
preference 
is 
that 
the 
origin 
would 
not 
have 
to
change 
when 
downsampling 
an 
image. 
I 
think 
that 
the 
personal 
preference
for 
one 
of 
the 
two 
options 
we 
have 
basically 
depends 
on 
whether 
one 
is
more 
interested 
in 
the 
overall 
dimensions 
of 
a 
dataset 
or 
whether 
one 
has
to 
do 
voxel-wise 
operations 
depending 
on 
the 
physical 
coordinate.

--> 
Therefore, 
I 
strongly 
ask 
other 
ITK 
users 
and 
developers 
to 
take 
part
in 
this 
discussion 
so 
that 
we 
can 
make 
a 
decision. 
This 
problem 
is 
really
of 
high 
importance 
since 
it 
is 
related 
to 
the 
core 
of 
the 
lib 
(namely 
the
ImageBase 
class) 
and 
it 
has 
a 
big 
influence 
on 
reconstruction 
problems. 
As
soon 
as 
we 
have 
come 
to 
a 
conclusion, 
I 
would 
file 
a 
bug 
report 
with 
the
proposed 
solution.

Thank 
you,
Andreas.

-----Original 
Message-----
From: 
Andreas 
Keil
Sent: 
Tuesday, 
December 
18, 
2007 
13:30
To: 
insight-users at itk.org
Subject: 
RE: 
[Insight-users] 
ITK 
Image 
Coordinates 
/ 
Problem 
with
PhysicalPoint 
to 
Index 
Conversion

Hi 
Maarten,

thank 
you 
for 
your 
reply 
(which 
other 
list 
subscribers 
may 
find 
below). 
My
initial 
posting 
was 
biased 
towards 
changing 
the 
ITK 
implementation 
rather
than 
the 
documentation 
for 
the 
following 
reasons:

Working 
with 
the 
physical 
coordinates 
of 
a 
voxel 
usually 
means 
that 
one
needs 
a 
single 
coordinate 
representation 
of 
a 
voxel. 
Algorithms 
relying
one 
this 
single 
physical 
coordinate 
would 
usually 
prefer 
the 
center 
of 
a
voxel 
for 
their 
space-related 
computations.

Moreover, 
ITK 
images 
only 
reflect 
the 
spacing 
between 
voxels 
which 
is 
not
necessarily 
equal 
to 
the 
width 
of 
a 
voxel. 
I 
am 
pretty 
sure 
that 
there 
are
cases 
where 
the 
voxel 
width 
is 
smaller 
than 
the 
spacing 
(in 
CT 
for
example). 
In 
this 
case, 
the 
term 
"spacing" 
would 
also 
be 
ambigious 
if 
it
would 
not 
refer 
to 
the 
distance 
between 
adjacent 
voxels' 
centers.

What 
do 
others 
think 
about 
this? 
I 
would 
appreciate 
a 
clearification 
as
the 
definitions 
of 
image 
origin 
and 
spacing 
really 
matter 
for 
my
application 
in 
3D 
reconstruction. 

Thank 
you,
Andreas.

-----Original 
Message-----
From: 
Maarten 
Nieber
Sent: 
Tuesday, 
December 
18, 
2007 
10:11
To: 
Andreas 
Keil
Subject: 
Re: 
[Insight-users] 
ITK 
Image 
Coordinates 
/ 
Problem 
with 
Physical
Point 
to 
Index 
Conversion

Hi 
Andreas,

I 
agree 
with 
you 
that 
according 
to 
the 
software 
guide, 
mapping 
(0.6, 
0.6)
should 
yield 
an 
index 
of 
(0,0).
On 
the 
other 
hand, 
I 
think 
that 
the 
definition 
of 
origin 
in 
the 
software
guide 
is 
not 
intuitive.
It 
says 
"the 
image 
origin 
is 
associated 
with 
the 
co-ordinates 
of 
the 
first
pixel 
in 
the 
image". 
Probably, 
it 
would 
be 
more 
accurate 
to 
say 
that 
the
image 
origin 
is 
associated 
with 
the 
>center< 
of 
the 
first 
pixel 
in 
the
image. 
However, 
by 
using 
this 
definition, 
the 
extent 
of 
an 
image 
with
origin 
(0,0) 
contains 
negative 
co-ordinates. 
I 
would 
find 
it 
more
intuitive 
if 
the 
extent 
of 
such 
an 
image 
is 
something 
like 
(0,0) 
- 
(size1,
size2). 
This 
would 
be 
the 
case 
if 
the 
origin 
of 
the 
image 
is 
associated
with 
the 
bottom-left 
corner 
of 
the 
first 
pixel.

If 
in 
the 
itk 
code, 
the 
origin 
of 
an 
itk 
image 
co-incides 
with 
the 
bottom
left 
corner 
of 
the 
first 
pixel, 
then 
I 
would 
prefer 
to 
change 
the 
software
guide 
and 
not 
the 
code.

Best 
regards
Maarten 
Nieber

----- 
Original 
Message 
----
From: 
Andreas 
Keil
To: 
insight-users at itk.org
Sent: 
Friday, 
December 
14, 
2007 
6:10:25 
PM
Subject: 
[Insight-users] 
ITK 
Image 
Coordinates 
/ 
Problem 
with 
Physical
Point 
to 
Index 
Conversion

Hi,

I 
think 
I 
have 
discovered 
an 
inconsistency 
between 
the 
ITK 
Software 
Guide
(p.40) 
and 
the 
implementation 
of 
itk::Image 
(all 
the 
methods 
taking
physical 
points 
/ 
continuous 
indices 
as 
arguments):

The 
trunctation 
of 
continuous 
index 
coordinates 
to 
integers 
does 
not 
yield
the 
correct 
pixel 
coordinates 
as 
expected 
by 
looking 
at 
the 
definition 
in
the 
software 
guide.

A 
simple 
example 
is:
Image 
spacing: 
1mm
Image 
origin: 
0mm/0mm

The 
pixel 
with 
index 
(1/1) 
would 
(according 
to 
the 
software 
guide) 
have
the 
following 
extents:
0.5mm/0.5mm 
to 
1.5mm/1.5mm

However, 
the 
physical 
point 
0.6mm/0.6mm 
gets 
mapped 
to 
the 
index 
0/0 
by
the 
method 
TransformPhysicalPointToIndex.

The 
solution 
would 
be 
to 
check 
those 
conversion 
methods 
as 
well 
as 
others
(like 
IsInside) 
and 
replace 
integer 
truncations 
with 
rounding.

If 
my 
reasoning 
is 
correct, 
I'll 
file 
a 
bug 
report. 
However, 
I'd 
like 
to
have 
some 
confirmation 
first.

Thanks,
Andreas.

_______________________________________________
Insight-users 
mailing 
list
Insight-users at itk.org
http://www.itk.org/mailman/listinfo/insight-users

      ____________________________________________________________________________________
Never miss a thing.  Make Yahoo your home page. 
http://www.yahoo.com/r/hs
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20080202/eef14321/attachment.htm