System and method for geodesic image matching using edge points interpolation
View Patent ↗A method for deformable registration of 2 digital images includes providing a pair of digital images, including a fixed image and a moving image, extracting a set of edge images from each image of the pair of images, each edge set being extracted at a different resolution, selecting a pair of edge images with a lowest resolution, determining a mapping from edge points of the fixed image to edge points of moving image using a geodesic thin plate spline interpolation, applying the mapping to a next higher resolution edge point image of the moving image, selecting a pair of edge images at a next higher resolution, where a moving edge image is the moving edge image to which the mapping has been applied, repeating the steps at a next higher resolution for all edge images in the set of edge images, and applying the mapping to an entire moving image.
1. A method for deformable registration of 2 digital images, comprising the steps of:
providing a pair of edge images, said pair including a fixed edge image and a moving edge image, each said edge image comprising a set of edge points extracted from a corresponding 2-dimensional digital image;
initializing moving points q i from fixed points p i , wherein said fixed points are selected from said fixed edge image;
initializing a set of N 2-dimensional curbs ζ i (t), iε[1,N], on a domain of said edge images parameterized by tε[0,T] associated with paths from said fixed points p i to moving points q i wherein
ζ
i
(
t
)
=
p
i
+
t
T
(
q
i
-
p
i
)
;
updating said curbs for parameters tε[2,T−1] using a gradient descent of an energy associated with paths from said landmark points p i to landmark points q i ;
calculating from said curbs for parameters tε[2,T−1] a set of coefficients a, b, α that define a diffeomorphic function g defined on said image domain wherein
g
(
p
i
)
=
~
q
i
,
∀
i
∈
[
1
,
N
]
,
and
g
(
x
)
=
ax
+
b
+
∑
i
=
1
N
α
i
f
p
i
(
x
)
,
wherein a is a 2×2 matrix, b is a 2D vector, α is a 2×N matrix, and ƒ p (x) are a set of basis functions;
calculating g from said coefficients;
updating said curbs for t=T from said diffeomorphic function g and a gradient descent of said constraint energy calculated for points transformed by g; and
repeating said steps of updating said curbs for parameters tε[2,T−1], calculating said coefficients a, b, α, calculating g, and updating said curbs for t=T until g converges,
wherein g is a deformation field that maps said fixed edge points to said moving edge points, and said curbs for t=T represent the image of the initial points as transformed by g.
2. The method of claim 1 , wherein each corresponding digital image comprises a plurality of intensities associated with a 2-dimensional grid of points, and further comprising extracting a set of edge images from said corresponding pair of digital images, each edge set being extracted at a different resolution;
selecting a pair of edge images with a lowest resolution;
applying said deformation field mapping to a next higher resolution edge point image of said moving image;
selecting a pair of edge images at a next higher resolution, wherein a moving edge image is the moving edge image to which said mapping has been applied;
calculating a deformation field and set of curbs for each set of next higher resolution edge images; and
applying said mapping to an entire moving image.
3. The method of claim 1 , wherein the q i are initialized as
q
i
=
p
i
-
λ
∂
C
F
∂
p
i
(
p
i
,
…
,
p
N
)
,
wherein
C
F
(
P
I
)
=
∑
x
,
y
(
S
F
(
x
,
y
)
-
S
M
(
x
,
y
)
)
2
is said constraint energy, S F and S M are structural intensities for a fixed image and a moving image, respectively, wherein
S
I
=
∑
i
∈
P
I
v
i
I
K
*
δ
(
P
i
I
)
,
where δ represents the Dirac distribution, K represents a Gaussian kernel, points P I ={P 1 I , . . . , P i I , . . . , P n I } are points defining the edges of image I, wherein I is either the moving or the fixed image, v I ={v 1 I , . . . , v i I , . . . , v n I } are their associated values, and γ is a predefined constant.
4. The method of claim 3 , wherein points in the moving set are constrained to be in a set
P
^
=
{
k
∈
N
|
K
*
δ
(
P
k
M
)
⋂
⋃
j
(
K
*
δ
(
P
j
M
)
)
≠
{
∅
}
}
.
5. The method of claim 1 , wherein updating said curbs for parameters tε[2,T−1] using a gradient descent comprises calculating
ζ
i
(
t
)
←
ζ
i
(
t
)
-
β
∂
E
∂
ζ
i
(
t
)
wherein
E
(
α
,
ζ
)
=
∑
i
,
j
=
1
N
∑
t
=
0
T
〈
α
i
(
t
)
,
α
j
(
t
)
〉
U
(
ζ
i
(
t
)
-
ζ
j
(
t
)
)
+
λ
∑
i
=
1
N
∑
t
=
0
T
-
1
Z
i
(
t
)
2
,
wherein
U
(
r
)
=
r
log
r
,
Z
i
(
t
)
=
D
t
(
ζ
i
)
-
a
(
t
)
ζ
i
(
t
)
-
b
(
t
)
-
∑
k
=
1
N
α
k
U
(
ζ
k
(
t
)
-
ζ
i
(
t
)
)
,
wherein for any time dependent function u, D(u)=T×(u(t+1)−u(t)),tε[0,T−1], and β is a predefined constant.
6. The method of claim 1 , wherein calculating said coefficients a, b, α comprises:
forming a 3×N matrix
P
=
[
1
ζ
1
x
(
t
-
1
)
ζ
1
y
(
t
-
1
)
1
ζ
2
x
(
t
-
1
)
ζ
2
y
(
t
-
1
)
…
…
…
1
ζ
N
x
(
t
-
1
)
ζ
N
y
(
t
-
1
)
]
;
forming an N×N matrix K of U(r i,j ) as
K
=
[
0
U
(
r
1
,
2
)
…
U
(
r
1
,
N
)
U
(
r
2
,
1
)
0
…
U
(
r
2
,
N
)
…
…
…
…
U
(
r
N
,
1
)
U
(
r
N
,
2
)
…
0
]
,
and matrix K λ =K+1/λK wherein λ is a smoothness parameter greater than 0, r i,j =|p i −p j | is a distance between the points i and j, and U(r)=r log(r),
forming a 3×3 matrix of zeros O;
forming a (N+3)×(N+3) matrix
L
=
[
K
λ
P
t
P
O
]
;
calculating L −1 Q wherein
Q
=
t
(
ζ
1
x
(
t
)
,
…
,
ζ
N
x
(
t
)
,
000
ζ
1
y
(
t
)
,
…
,
ζ
N
y
(
t
)
,
000
)
,
wherein t (α|a b)=L −1 Q is an (N+3)×2 matrix, and wherein | represents concatenation.
7. The method of claim 1 , wherein calculating g comprises initializing g to an identity mapping, and evaluating
g
(
x
)
←
a
(
t
)
g
(
x
)
+
b
(
t
)
+
∑
i
=
1
N
α
i
(
t
)
f
ζ
i
(
t
)
(
g
(
x
)
)
for all t=0, . . . , T.
8. The method of claim 1 , wherein updating said curbs for t=T comprises calculating
ζ
(
T
)
←
g
(
ζ
(
T
)
)
-
γ
∂
C
F
∂
g
(
ζ
(
T
)
)
(
g
(
ζ
(
T
)
)
)
,
C
F
(
P
I
)
=
∑
x
,
y
(
S
F
(
x
,
y
)
-
S
M
(
x
,
y
)
)
2
is said constraint energy, S F and S M are structural intensities for a fixed image and a moving image, respectively, wherein
S
I
=
∑
i
∈
P
I
v
i
I
K
*
δ
(
P
i
I
)
,
where δ represents the Dirac distribution, K represents a Gaussian kernel, points P I ={P 1 I , . . . , P i I , . . . , P n I } are points defining the edges of image I, wherein I is either the moving or the fixed image, v I ={v 1 I , . . . , v i I , . . . , v n I } are their associated values, g(ζ(T)) represent the moving points, and γ is a predefined constant.
9. A non-transitory program storage device readable by a computer, tangibly embodying a program of instructions executable by the computer to perform the method steps for deformable registration of 2 digital images, said method comprising the steps of:
providing a pair of edge images, said pair including a fixed edge image and a moving edge image, each said edge image comprising a set of edge points extracted from a corresponding 2-dimensional digital image;
initializing moving points q i from fixed points p i , wherein said fixed points are selected from said fixed edge image;
initializing a set of N 2-dimensional curbs ζ i (t), iε[1,N], on a domain of said edge images parameterized by tε[0,T] associated with paths from said fixed points p i to moving points q i wherein
ζ
i
(
t
)
=
p
i
+
t
T
(
q
i
-
p
i
)
;
updating said curbs for parameters tε[2,T−1] using a gradient descent of an energy associated with paths from said landmark points p i to landmark points q i ;
calculating from said curbs for parameters tε[2,T−1] a set of coefficients a, b, α that define a diffeomorphic function g defined on said image domain wherein
g
(
p
i
)
=
~
q
i
,
∀
i
∈
[
1
,
N
]
,
and
g
(
x
)
=
ax
+
b
+
∑
i
=
1
N
α
i
f
p
i
(
x
)
,
wherein a is a 2×2 matrix, b is a 2D vector, α is a 2×N matrix, and ƒ p (x) are a set of basis functions;
calculating g from said coefficients;
updating said curbs for t=T from said diffeomorphic function g and a gradient descent of said constraint energy calculated for points transformed by g; and
repeating said steps of updating said curbs for parameters tε[2,T−1], calculating said coefficients a, b, α, calculating g, and updating said curbs for t=T until g converges,
wherein g is a deformation field that maps said fixed edge points to said moving edge points, and said curbs for t=T represent the image of the initial points as transformed by g.
10. The computer readable program storage device of claim 9 , wherein each corresponding digital image comprises a plurality of intensities associated with a 2-dimensional grid of points, and further comprising extracting a set of edge images from said corresponding pair of digital images, each edge set being extracted at a different resolution;
selecting a pair of edge images with a lowest resolution;
applying said deformation field mapping to a next higher resolution edge point image of said moving image;
selecting a pair of edge images at a next higher resolution, wherein a moving edge image is the moving edge image to which said mapping has been applied;
calculating a deformation field and set of curbs for each set of next higher resolution edge images; and
applying said mapping to an entire moving image.
11. The computer readable program storage device of claim 10 , the method further comprising reducing a number of points in each of said edge images, prior to determining said mapping.
12. The computer readable program storage device of claim 11 , wherein reducing a number of points comprises following a curve defined by a set of edge points, keeping one point every α, wherein α is a predetermined distance interval, and keeping those points wherein an angle made by three points in β points is greater than θ, wherein β is a predetermined number of points and θ is a predetermined angle.
13. The computer readable program storage device of claim 10 , the method further comprising using a Gaussian pyramid to generate images at progressively lower resolutions to extract edge sets at different resolutions.
14. The computer readable program storage device of claim 9 , wherein the q i are initialized as
q
i
=
p
i
-
γ
∂
C
F
∂
p
i
(
p
i
,
…
,
p
N
)
,
wherein
C
F
(
P
I
)
=
∑
x
,
y
(
S
F
(
x
,
y
)
-
S
M
(
x
,
y
)
)
2
is said constraint energy, S F and S M are structural intensities for a fixed image and a moving image, respectively, wherein
S
I
=
∑
i
∈
P
I
v
i
I
K
*
δ
(
P
i
I
)
,
where δ represents the Dirac distribution, K represents a Gaussian kernel, points P I ={P 1 I , . . . , P i I , . . . , P n I } are points defining the edges of image I, wherein I is either the moving or the fixed image, v I ={v 1 I , . . . , v i I , . . . , v n I } are their associated values, and γ is a predefined constant.
15. The computer readable program storage device of claim 14 , wherein points in the moving set are constrained to be in a set
P
^
=
{
k
∈
N
|
K
*
δ
(
P
k
M
)
⋂
⋃
j
(
K
*
δ
(
P
j
M
)
)
≠
{
∅
}
}
.
16. The computer readable program storage device of claim 9 , wherein updating said curbs for parameters tε[2,T−1] using a gradient descent comprises calculating
ζ
i
(
t
)
←
ζ
i
(
t
)
-
β
∂
E
∂
ζ
i
(
t
)
wherein
E
(
α
,
ζ
)
=
∑
i
,
j
=
1
N
∑
t
=
0
T
〈
α
i
(
t
)
,
α
j
(
t
)
〉
U
(
ζ
i
(
t
)
-
ζ
j
(
t
)
)
+
λ
∑
i
=
1
N
∑
t
=
0
T
-
1
Z
i
(
t
)
2
,
wherein
U
(
r
)
=
r
log
r
,
Z
i
(
t
)
=
D
t
(
ζ
i
)
-
a
(
t
)
ζ
i
(
t
)
-
b
(
t
)
-
∑
k
=
1
N
α
k
U
(
ζ
k
(
t
)
-
ζ
i
(
t
)
)
,
wherein for any time dependent function u, D t (u)=T×(u(t+1)−u(t)), tε[0,T−1], and β is a predefined constant.
17. The computer readable program storage device of claim 9 , wherein calculating said coefficients a, b, α comprises:
forming a 3×N matrix
P
=
[
1
ζ
1
x
(
t
-
1
)
ζ
1
y
(
t
-
1
)
1
ζ
2
x
(
t
-
1
)
ζ
2
y
(
t
-
1
)
…
…
…
1
ζ
N
x
(
t
-
1
)
ζ
N
y
(
t
-
1
)
]
;
forming an N×N matrix K of U(r i,j ) as
K
=
[
0
U
(
r
1
,
2
)
…
U
(
r
1
,
N
)
U
(
r
2
,
1
)
0
…
U
(
r
2
,
N
)
…
…
…
…
U
(
r
N
,
1
)
U
(
r
N
,
2
)
…
0
]
,
matrix K λ =K+1/λK wherein λ is a smoothness parameter greater than 0, r i,j =|p i −p j | is a distance between the points i and j, and U(r)=r log(r),
forming a 3×3 matrix of zeros O;
forming a (N+3)×(N+3) matrix
L
=
[
K
λ
P
t
P
O
]
;
calculating L −1 Q wherein
Q
=
t
(
ζ
1
x
(
t
)
,
…
,
ζ
N
x
(
t
)
,
000
ζ
1
y
(
t
)
,
…
,
ζ
N
y
(
t
)
,
000
)
,
wherein t (α|a b)=L −1 is an (N+3)×2 matrix and wherein | represents concatenation.
18. The computer readable program storage device of claim 9 , wherein calculating g comprises initializing g to an identity mapping, and evaluating
g
(
x
)
←
a
(
t
)
g
(
x
)
+
b
(
t
)
+
∑
i
=
1
N
α
i
(
t
)
f
ζ
i
(
t
)
(
g
(
x
)
)
for
all
t
=
0
,
…
,
T
.
19. The computer readable program storage device of claim 9 , wherein updating said curbs for t=T comprises calculating
ζ
(
T
)
←
g
(
ζ
(
T
)
)
-
γ
∂
C
F
∂
g
(
ζ
(
T
)
)
(
g
(
ζ
(
T
)
)
)
,
wherein
C
F
(
P
I
)
=
∑
x
,
y
(
S
F
(
x
,
y
)
-
S
M
(
x
,
y
)
)
2
said constraint energy, S F and S M are structural intensities for a fixed image and a moving image, respectively, wherein
S
I
=
∑
i
∈
P
I
v
i
I
K
*
δ
(
P
i
I
)
,
where δ represents the Dirac distribution, K represents a Gaussian kernel, points P I ={P 1 I , . . . , P i I , . . . , P n I } are points defining the edges of image I, wherein I is either the moving or the fixed image, v I ={v 1 I , . . . , v i I , . . . , v n I } are their associated values, g(ζ(T)) represent the moving points, and γ is a predefined constant.