Computing 3D Medial Axis for Chamfer Distances - Semantic Scholar
Laboratoire d’Informatique de Marseille
Computing 3D Medial Axis for Chamfer Distances Eric Remy and Edouard Thiel {Eric.Remy,Edouard.Thiel}@lim.univ-mrs.fr
http://www.lim.univ-mrs.fr/~eremy
Workspace B Discrete square or cubic grid E (Z2 , Z3 ) B Associated to integer values ⇒ Digital images and volumes
Discrete Distances • Frequently used tools in image analysis (distance map, medial axis, skeleton) • Norm l1 = d6 (p, q) = |xp − xq | + |yp − yq | + |zp − zq |
(city block distance)
• Norm l∞ = d26 (p, q) = max(|xp − xq |, |yp − yq |, |zp − zq |) (chessboard distance) ⇒ Chamfer distances dC
Chamfer Mask B Weighted neighbourhood: MC =
(
vi
z }| { M i = (xi , y i , z i , w i ) , 1 ≤ i ≤ m
)
B Set of points, visible from the origin: gcd(x, y, z) = 1 1 3 B By symmetry, the generator MgC ⊆ 48 Z gives MC 1 3 B 48 Z = (x, y, z) ∈ Z3 : 0 ≤ z ≤ y ≤ x
11
11
7
5
5 11
7 11
5
11 → − vb 7 11 MgC → va 5 − 7 11
11
z
− → vc
y
O x
a (1, 0, 0) b (1, 1, 0)
1 3 48 Z
c (1, 1, 1) − → vb − → va
d (2, 1, 0) e (2, 1, 1) f (2, 2, 1)
Path 5
11
5 7 5 A
11 11 5
B 7
W (P◦ ) = 33 W (P◦ ) = 23 W (P◦ ) = 22
11
B Sequence of vectors from MC : P = n1 v1 + . . . + nm vm , ni ≥ 0 Pm B Cost of the path: W (P) = i=1 ni wi
Chamfer Distance dC (A, B) B Minimum of the cost of all finite length paths from A to B : dC (A, B) = min{W (PA→B )} B dC : E → N B Enable to aproximate the Euclidean distance dE B Proof in [Ver91], that this function is a distance: dC (A, B) ≥ 0 ; dC (A, B) = 0 ⇐⇒ A = B
positive defined
dC (A, B) = dC (B, A) dC (A, B) ≤ dC (A, C) + dC (C, B) −−→ −−→ dC (λ.AB) = |λ|.dC (AB)
symmetric triangular inequality
⇒ We shown how to obtain homogeneity in [Rem00a, Rem00b] ⇒ The chamfer balls are convex and symmetric.
homogeneity
Distance Map B Distance Transform (DT): To set each point to its distance to the background
1 1 1 0 0
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
0 1 1 1 0
d3,4 7→
3 3 3 0 0
3 6 4 3 3
3 6 7 6 3
3 6 9 6 3
3 6 8 6 3
3 4 6 4 3
0 3 3 3 0
B Reverse Distance Transform (RDT): To recover the object from some points and their distance to the background → Two passes on the image [Ros66]
Medial Axis (MA) Definition 1 (Maximal Ball) A maximal ball is not included in any single other ball in the shape. Definition 2 (Medial Axis) The Medial Axis of a shape is the set of centres of maximal balls in the shape.
Medial Axis using DT B Reverse Ball: Bd−1 (p, r) = { q : r − d(p, q) > 0 } B How to detect locally which p ∈ MA, using DT image?
3 3 3 0 0
3 6 4 3 3
3 6 7 6 3
3 6 9 6 3
3 6 8 6 3
3 4 6 4 3
0 3 3 3 0
1 1 4 2 5 1 4 1
2 5 8 5 2
1 4 1 5 2 4 1 1
Medial Axis Extraction 1. Local Maxima Criterion (LMC) [Ros66]: → p ∈ MA ⇐⇒ ∀M ∈ M , DT [p + − v ] < DT [p] + w i
C
i
i
B Only valid for l1 and l∞ 2. Equivalent Balls [Arc88]: B To lower some values on the DT, before LMC → − p+− vi → vi R
r
B Limited to da,b in 2D and da,b,c in 3D p
3. Look-Up Tables [Bor91]: p ∈ MA ⇐⇒ ∀Mi ∈ MLut B More general method
R
z }| { → − → − g , DT [p + vi ] < Lut[ vi ][DT [p]] | {z } r
B Suppose that a local test in the neighbourhood MLut is sufficient → B How to compute for any given r and − vi g the smallest → radius R = Lut[− v g ][r] which forbids p ∈ MA ? i
“Na¨ıve” Method B Balls are convex and symmetric 1 3 Z ⇒ restriction to 48 −1 Bd (q, R+ ) → B For a given r and − v g , what is the smallest 1 3 Z 48
Bd−1 (p, r) − → vg q
p
r
radius R which forbids p ∈ MA ? B Method: To decrease R+ while Bd−1 (p, r) ⊆ Bd−1 (q, R+ ) B Problems:
R
R+
1. All the balls must be generated 2. One scan for each overlapping test
Balls and Reverse Balls 6 4 3 4 6 3 0 3 6 4 3 4 6 (a)
1 3 4 3 1 4 7 4 1 3 4 3 1 (b)
Bd (p, 6)
Bd−1 (p, 6 + 1)
B Ball: B Reverse Ball: B Lemma:
Bd (p, r)
=
{ q : d(p, q) ≤ r }
(a)
Bd−1 (p, r)
=
{ q : r − d(p, q) > 0 }
(b)
Bd (p, r)
=
Bd−1 (p, r + 1) .
Proposed Method
1 3 Z 48
− → vg
B Translate the balls to the origin
p2
B Search the maximum of dC (O, p2 ) p1
− → vg r
O
0 0 0 0 0 0 0
0 0 0 0 0 6 4
0 0 0 0 12 9 8
0 0 0 18 15 13 12
R
0 0 24 21 18 17 16
0 30 27 24 22 21 20
36 33 30 27 26 25 24
B CT g : Distance to the origin → computed a single time in a single pass only with MgC Example: Search for Bd (12) = Bd−1 (12 + 1): → → For − v g = (1, 0), → Lut[− v g ][12 + 1] = max{15 + 1, 16 + 1} → → For − v g = (1, 1), → Lut[− v g ][12 + 1] = max{17 + 1, 18 + 1} B a single pass on
1 3 48 Z
→ for each − vg
Inclusions Correction B Partial ordering → → • r1 < r2 but Lut[− v g ][r1 ] 6≤ Lut[− v g ][r2 ] B For any distance, ∀r1 , r2
r1 < r2 ⇒ Bd (r1 ) ⊆ Bd (r2 )
B Total ordering → → • Thus Lut[− v g ][r2 ] at least equals Lut[− v g ][r1 ] → B a single pass for each column Lut[− v g ][ ] → ⇒ Final Lut array (where each column Lut[− v g ][ ] may be computed independently)
MgLut Computation (1/2) : Verification B Problem raised in [Thi94] B By definition, the medial axis of any ball Bd (O, r) equals its sole centre O.
21 14 18 7 11 16 0 5 10 15
(a)
7 (b) 14 11 21 18 14 28 25 20 15
28 25 22 21 20
35 32 29 27 26 25
5 7 5 10 5 10 5
B Hypothesis: MgLut is sufficient for any ball of radius r ≤ rKnown B Problem: Ensure that MgLut is sufficient until r = rTarget B Proposed Method: For each r ∈ [rKnown , rTarget ] do 1. By thresholding of CT g find Bd (O, r) ∩
1 3 48 Z g
(a)
2. Use a restricted and faster DT to obtain DT (b) 3. Verify that ∀p 6= O , p 6∈ MA else..
MgLut Conputation (2/2) : Choice of the direction to add B The current neighbourhood MgLut (around p) is not sufficient to find that Bd (O, r) forbids p ∈ MA
Bd −1 (O, r)
B But Bd (O, r) exists and to find it one must look in −→ −→ direction pO = −Op B Thus we extend the mask −→ MgLut ← MgLut ∪ (Op, CT g [p])
p
−→ B And compute the corresponding Lut[Op][ ]
O r
B Then we resume the verification process from p ⇒ Starting from MgLut = ∅ and rKnown = 0, we compute a sufficient mask until any given rTarget
r
a b c f
k
d 3e j
Example : d11,16,19,j=45
11 12 17 20 36 50 28 91 46
Results verified until r = 1200
19 28 33 36 52 65 44 106 62
16 23 28 31 46 61 39 102 57 22 31 36 39 55 69 46 110 65
MgC =
MgLut
a = (1, 0, 0, 11) b = (1, 1, 0, 16) c = (1, 1, 1, 19) j = (3, 2, 1, 45)
a = (1, 0, 0, 11) k = (3, 2, 2, 49) b = (1, 1, 0, 16) d = (2, 1, 0, 27) = c = (1, 1, 1, 19) 3e = (6, 3, 3, 90) f = (2, 2, 1, 35) j = (3, 2, 1, 45)
27 34 39 42 57 72 50 113 68 30 39 44 46 62 76 55 117 73 32 42 46 50 65 80 57 121 76 33
81
121
35 45 50 53 68 83 61 124 79 38 46 52 55 71 84 62 125 81
41 50 55 58 74 88 66 129 84 43 53 57 61 76 91 68 132 87 44 45 . . . . . .
62 78 91 . . .
79 . . . . . . . . .
178
212
194
228
132 . . .
. . .
. . .
Conclusion B We introduce the more general MgLut (6= MgC ) which removes the problem found in [Thi94] B General method usable on any discrete distance B Particularly effective for discrete norms [Rem00a] (use of symmetries, size of MgLut is bounded, etc.) B Computation of both MgLut and Lut B Automatic and systematic proof of their validity
References [Arc88]
C. Arcelli and G. Sanniti di Baja. Finding local maxima in a pseudo-euclidean distance transform. Computer Vision, Graphics and Image Processing, 43:361–367, 1988.
[Bor91]
G. Borgefors, I. Ragnemalm, and G. Sanniti di Baja. The Euclidean Distance Transform : finding the local maxima and reconstructing the shape. In 7 th Scandinavian Conf. on Image Analysis, volume 2, pages 974–981, Aalborg, Denmark, 1991.
[Rem00a] E. Remy and E. Thiel. Optimizing 3D chamfer masks with distance constraints. In 7 th IWCIA, Int. Workshop on Combinatorial Image Analysis, pages 39–56, Caen, July 2000. [Rem00b] E. Remy and E. Thiel. Structures dans les sph` eres de chanfrein. In 12 e`me RFIA, congr` es Reconnaissance des Formes et I.A, volume 1, pages 483–492, Paris, Fev 2000. [Ros66]
A. Rosenfeld and J.L. Pfaltz. Sequential operations in digital picture processing. Journal of ACM, 13(4):471–494, 1966.
[Thi94]
E. Thiel. Les distances de chanfrein en analyse d’images : fondements et applications. PhD thesis, UJF, Grenoble, Sept 1994. http://www.lim.univ-mrs.fr/~thiel/these .
[Ver91]
J.H. Verwer. Distance transforms: metrics, algorithms and applications. PhD thesis, Technische Universiteit, Delft, 1991.
Computing 3D Medial Axis for Chamfer Distances Eric Remy and Edouard Thiel {Eric.Remy,Edouard.Thiel}@lim.univ-mrs.fr
http://www.lim.univ-mrs.fr/~eremy
Workspace B Discrete square or cubic grid E (Z2 , Z3 ) B Associated to integer values ⇒ Digital images and volumes
Discrete Distances • Frequently used tools in image analysis (distance map, medial axis, skeleton) • Norm l1 = d6 (p, q) = |xp − xq | + |yp − yq | + |zp − zq |
(city block distance)
• Norm l∞ = d26 (p, q) = max(|xp − xq |, |yp − yq |, |zp − zq |) (chessboard distance) ⇒ Chamfer distances dC
Chamfer Mask B Weighted neighbourhood: MC =
(
vi
z }| { M i = (xi , y i , z i , w i ) , 1 ≤ i ≤ m
)
B Set of points, visible from the origin: gcd(x, y, z) = 1 1 3 B By symmetry, the generator MgC ⊆ 48 Z gives MC 1 3 B 48 Z = (x, y, z) ∈ Z3 : 0 ≤ z ≤ y ≤ x
11
11
7
5
5 11
7 11
5
11 → − vb 7 11 MgC → va 5 − 7 11
11
z
− → vc
y
O x
a (1, 0, 0) b (1, 1, 0)
1 3 48 Z
c (1, 1, 1) − → vb − → va
d (2, 1, 0) e (2, 1, 1) f (2, 2, 1)
Path 5
11
5 7 5 A
11 11 5
B 7
W (P◦ ) = 33 W (P◦ ) = 23 W (P◦ ) = 22
11
B Sequence of vectors from MC : P = n1 v1 + . . . + nm vm , ni ≥ 0 Pm B Cost of the path: W (P) = i=1 ni wi
Chamfer Distance dC (A, B) B Minimum of the cost of all finite length paths from A to B : dC (A, B) = min{W (PA→B )} B dC : E → N B Enable to aproximate the Euclidean distance dE B Proof in [Ver91], that this function is a distance: dC (A, B) ≥ 0 ; dC (A, B) = 0 ⇐⇒ A = B
positive defined
dC (A, B) = dC (B, A) dC (A, B) ≤ dC (A, C) + dC (C, B) −−→ −−→ dC (λ.AB) = |λ|.dC (AB)
symmetric triangular inequality
⇒ We shown how to obtain homogeneity in [Rem00a, Rem00b] ⇒ The chamfer balls are convex and symmetric.
homogeneity
Distance Map B Distance Transform (DT): To set each point to its distance to the background
1 1 1 0 0
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
0 1 1 1 0
d3,4 7→
3 3 3 0 0
3 6 4 3 3
3 6 7 6 3
3 6 9 6 3
3 6 8 6 3
3 4 6 4 3
0 3 3 3 0
B Reverse Distance Transform (RDT): To recover the object from some points and their distance to the background → Two passes on the image [Ros66]
Medial Axis (MA) Definition 1 (Maximal Ball) A maximal ball is not included in any single other ball in the shape. Definition 2 (Medial Axis) The Medial Axis of a shape is the set of centres of maximal balls in the shape.
Medial Axis using DT B Reverse Ball: Bd−1 (p, r) = { q : r − d(p, q) > 0 } B How to detect locally which p ∈ MA, using DT image?
3 3 3 0 0
3 6 4 3 3
3 6 7 6 3
3 6 9 6 3
3 6 8 6 3
3 4 6 4 3
0 3 3 3 0
1 1 4 2 5 1 4 1
2 5 8 5 2
1 4 1 5 2 4 1 1
Medial Axis Extraction 1. Local Maxima Criterion (LMC) [Ros66]: → p ∈ MA ⇐⇒ ∀M ∈ M , DT [p + − v ] < DT [p] + w i
C
i
i
B Only valid for l1 and l∞ 2. Equivalent Balls [Arc88]: B To lower some values on the DT, before LMC → − p+− vi → vi R
r
B Limited to da,b in 2D and da,b,c in 3D p
3. Look-Up Tables [Bor91]: p ∈ MA ⇐⇒ ∀Mi ∈ MLut B More general method
R
z }| { → − → − g , DT [p + vi ] < Lut[ vi ][DT [p]] | {z } r
B Suppose that a local test in the neighbourhood MLut is sufficient → B How to compute for any given r and − vi g the smallest → radius R = Lut[− v g ][r] which forbids p ∈ MA ? i
“Na¨ıve” Method B Balls are convex and symmetric 1 3 Z ⇒ restriction to 48 −1 Bd (q, R+ ) → B For a given r and − v g , what is the smallest 1 3 Z 48
Bd−1 (p, r) − → vg q
p
r
radius R which forbids p ∈ MA ? B Method: To decrease R+ while Bd−1 (p, r) ⊆ Bd−1 (q, R+ ) B Problems:
R
R+
1. All the balls must be generated 2. One scan for each overlapping test
Balls and Reverse Balls 6 4 3 4 6 3 0 3 6 4 3 4 6 (a)
1 3 4 3 1 4 7 4 1 3 4 3 1 (b)
Bd (p, 6)
Bd−1 (p, 6 + 1)
B Ball: B Reverse Ball: B Lemma:
Bd (p, r)
=
{ q : d(p, q) ≤ r }
(a)
Bd−1 (p, r)
=
{ q : r − d(p, q) > 0 }
(b)
Bd (p, r)
=
Bd−1 (p, r + 1) .
Proposed Method
1 3 Z 48
− → vg
B Translate the balls to the origin
p2
B Search the maximum of dC (O, p2 ) p1
− → vg r
O
0 0 0 0 0 0 0
0 0 0 0 0 6 4
0 0 0 0 12 9 8
0 0 0 18 15 13 12
R
0 0 24 21 18 17 16
0 30 27 24 22 21 20
36 33 30 27 26 25 24
B CT g : Distance to the origin → computed a single time in a single pass only with MgC Example: Search for Bd (12) = Bd−1 (12 + 1): → → For − v g = (1, 0), → Lut[− v g ][12 + 1] = max{15 + 1, 16 + 1} → → For − v g = (1, 1), → Lut[− v g ][12 + 1] = max{17 + 1, 18 + 1} B a single pass on
1 3 48 Z
→ for each − vg
Inclusions Correction B Partial ordering → → • r1 < r2 but Lut[− v g ][r1 ] 6≤ Lut[− v g ][r2 ] B For any distance, ∀r1 , r2
r1 < r2 ⇒ Bd (r1 ) ⊆ Bd (r2 )
B Total ordering → → • Thus Lut[− v g ][r2 ] at least equals Lut[− v g ][r1 ] → B a single pass for each column Lut[− v g ][ ] → ⇒ Final Lut array (where each column Lut[− v g ][ ] may be computed independently)
MgLut Computation (1/2) : Verification B Problem raised in [Thi94] B By definition, the medial axis of any ball Bd (O, r) equals its sole centre O.
21 14 18 7 11 16 0 5 10 15
(a)
7 (b) 14 11 21 18 14 28 25 20 15
28 25 22 21 20
35 32 29 27 26 25
5 7 5 10 5 10 5
B Hypothesis: MgLut is sufficient for any ball of radius r ≤ rKnown B Problem: Ensure that MgLut is sufficient until r = rTarget B Proposed Method: For each r ∈ [rKnown , rTarget ] do 1. By thresholding of CT g find Bd (O, r) ∩
1 3 48 Z g
(a)
2. Use a restricted and faster DT to obtain DT (b) 3. Verify that ∀p 6= O , p 6∈ MA else..
MgLut Conputation (2/2) : Choice of the direction to add B The current neighbourhood MgLut (around p) is not sufficient to find that Bd (O, r) forbids p ∈ MA
Bd −1 (O, r)
B But Bd (O, r) exists and to find it one must look in −→ −→ direction pO = −Op B Thus we extend the mask −→ MgLut ← MgLut ∪ (Op, CT g [p])
p
−→ B And compute the corresponding Lut[Op][ ]
O r
B Then we resume the verification process from p ⇒ Starting from MgLut = ∅ and rKnown = 0, we compute a sufficient mask until any given rTarget
r
a b c f
k
d 3e j
Example : d11,16,19,j=45
11 12 17 20 36 50 28 91 46
Results verified until r = 1200
19 28 33 36 52 65 44 106 62
16 23 28 31 46 61 39 102 57 22 31 36 39 55 69 46 110 65
MgC =
MgLut
a = (1, 0, 0, 11) b = (1, 1, 0, 16) c = (1, 1, 1, 19) j = (3, 2, 1, 45)
a = (1, 0, 0, 11) k = (3, 2, 2, 49) b = (1, 1, 0, 16) d = (2, 1, 0, 27) = c = (1, 1, 1, 19) 3e = (6, 3, 3, 90) f = (2, 2, 1, 35) j = (3, 2, 1, 45)
27 34 39 42 57 72 50 113 68 30 39 44 46 62 76 55 117 73 32 42 46 50 65 80 57 121 76 33
81
121
35 45 50 53 68 83 61 124 79 38 46 52 55 71 84 62 125 81
41 50 55 58 74 88 66 129 84 43 53 57 61 76 91 68 132 87 44 45 . . . . . .
62 78 91 . . .
79 . . . . . . . . .
178
212
194
228
132 . . .
. . .
. . .
Conclusion B We introduce the more general MgLut (6= MgC ) which removes the problem found in [Thi94] B General method usable on any discrete distance B Particularly effective for discrete norms [Rem00a] (use of symmetries, size of MgLut is bounded, etc.) B Computation of both MgLut and Lut B Automatic and systematic proof of their validity
References [Arc88]
C. Arcelli and G. Sanniti di Baja. Finding local maxima in a pseudo-euclidean distance transform. Computer Vision, Graphics and Image Processing, 43:361–367, 1988.
[Bor91]
G. Borgefors, I. Ragnemalm, and G. Sanniti di Baja. The Euclidean Distance Transform : finding the local maxima and reconstructing the shape. In 7 th Scandinavian Conf. on Image Analysis, volume 2, pages 974–981, Aalborg, Denmark, 1991.
[Rem00a] E. Remy and E. Thiel. Optimizing 3D chamfer masks with distance constraints. In 7 th IWCIA, Int. Workshop on Combinatorial Image Analysis, pages 39–56, Caen, July 2000. [Rem00b] E. Remy and E. Thiel. Structures dans les sph` eres de chanfrein. In 12 e`me RFIA, congr` es Reconnaissance des Formes et I.A, volume 1, pages 483–492, Paris, Fev 2000. [Ros66]
A. Rosenfeld and J.L. Pfaltz. Sequential operations in digital picture processing. Journal of ACM, 13(4):471–494, 1966.
[Thi94]
E. Thiel. Les distances de chanfrein en analyse d’images : fondements et applications. PhD thesis, UJF, Grenoble, Sept 1994. http://www.lim.univ-mrs.fr/~thiel/these .
[Ver91]
J.H. Verwer. Distance transforms: metrics, algorithms and applications. PhD thesis, Technische Universiteit, Delft, 1991.
