Surface flattening - from local to global Emil Saucan EE - - PowerPoint PPT Presentation

Surface flattening - from local to global Emil Saucan

EE Department, Technion Joint work with

Eli Appleboim and Yehoshua Y. Zeevi.

Stony Brook

September 11, 2008.

Flattening of folded surfaces S2 ⊂ R3, (i.e. “2D” represen- tation by flattening of “3D” object scans) for representa- tion and analysis of medical images is a fundamental step, required in medical volumetric imaging, for example in:

• fMRI scans of the brain cortex, e.g. to better observe and

follow up developments of neural activity within the folds.

• CT virtual colonoscopy, for determining the presence of

colon pathologies such as small polyps.

1

But, because of the extensive folding of the colon and cor- tex, the obtained images are distorted and this renders them largely ineffectual for Medical (Sci- entific) purposes, since one cannot perform any relevant measurements on them. Therefore, for the diagnosis to be accurate, it is essential that the geometric distortion (in terms of distances, area, angles) be minimal.

2

One would hope for maps that at least preserve angles precisely, i.e. conformal mappings. There exists an extensive literature on conformal repre- sentation.∗ In fact mainly produce only quasi-conformal mappings – as approximations to conformal mappings – without any computation of the dilatation (definition to follow!) How ever, one can admit his/hers limitations and (1) Settle for the measurement of approximate length (area, angles, etc.) hence (2) Try and estimate the measurement error or lengths dis- tortion.

∗Hurdal-Bowers-Stephenson, Haker-Tannenbaum, Gu-Yau, etc.

3

Formally these are provided by the following definitions: Definition 1 Let D ⊂ R3 be a domain. A homeomorphism f : D → R3 is called a quasi-isometry (or a bi-lipschitz mapping), if there exists 1 ≤ C < ∞, such that 1 C|p1 − p2| ≤ |f(p1) − f(p2)| < C|p1 − p2| , for all p1, p2 ∈ D ; where | · | denotes (Euclidean) distance. The number C(f) = min{C | f is a quasi − isometry} is called the minimal distortion of f (in D). Remark 2 In fact, one can define quasi-isometries between any two metric spaces (X, d) and (Y, ρ). Remark 3 For the case of a surface in R3, distances are the induced intrinsic distances on the surface.

4

Remark 4 The requirement to almost preserve distances is to strong, and not always necessary. For many imple- mentations, such as Human Brain Mapping – where shape and relative areas are important suffices to almost preserve angles, i.e quasi-conformal map- pings. Remark 5 f quasi-isometry ⇒ f quasi-conformal; but f quasi-conformal f quasi-isometry.

5

We are now introducing the technical definition of quasi- conformality, but first we have to introduce some prelimi- nary notations and definitions: Let D ⊂ R3 be a domain, let f : D ∼ → f(D) be a homeomor- phism and let p ∈ D. We make the following notations: L(p) = Lf(p) = lim sup

x→p

|f(x) − f(p)| |x − p| ; l(p) = lf(p) = lim inf

x→p

|f(x) − f(p)| |x − p| ; J(p) = Jf(p) = lim sup

r→0

V ol

• f(B3(p, r))
• V ol
• B3(p, r)
• .

6

Definition 6 L(p), l(p) are called the maximal, respective minimal stretching, of f. Remark 7 If f is differentiable then Jf(p) = |Jacobian(f)(p)|. Definition 8 Let f be such that Jf > 0. Then the maximal and minimal dilatation of f are defined as follows: KI(f) = KI(f, p) =

• sup

p∈D

J(p) l3(p) ; KO(f) = KO(f, p) =

• sup

p∈D

L3(p) J(p) . Definition 9 The number K(f) = max{KI(f), KO(f)} is called the dilatation of f.

7

Definition 10 f : D

∼

→ f(D) is called quasi-conformal iff K(f) is bounded on D. The connections between the different dilatations are given by: Lemma 11 KO ≤ K2

I and KI ≤ K2 O .

Moreover, we also have the following relationship between distortion and dilatations: Lemma 12 KI(f) ≤ (C(f))2 and KO(f) ≤ (C(f))2 .

8

We also get the following geometric interpretations of quasi- conformal mappings:

• take infinitesimal balls into infinitesimal ellipsoids,
• take almost balls into almost ellipsoids
• distort infinitesimal spheres by a constant factor

9

Since we are interested in images of 3-dim objects on 2-dim “screens”, i.e. projections∗, we are conducted to ask the following questions:

∗and since human vision is essentially and intrinsically 2-dimensional

10

Question 1 When is orthogonal projection a quasi-isomorphism (quasi-conformal mapping)? Question 2 And if it is, what are its distortion and dilata- tion? The answer to these questions is to be found in the classical work of F. Gehring and Y. V¨ ais¨ al¨ a (1965). It applies for surfaces that do not “deviate” (or “turn”) “to fast” i.e. that satisfy:

11

The Geometric (Gehring) Condition: We say S ⊂ R3 satisfies the Gehring Condition if for any p ∈ S there exists

• np such that for any ε > 0, there exists Up ≃ D2, such that

for any q1, q2 ∈ Up the acute angle ∡(q1q2, np) ≥ α, where: 0 < inf

p∈S αp < π

2 − ε .

S

p n

_ α α

~ _D

q q

1 2

2

Up

p

Example 13 Any surface in S ⊂ R3 that admits a well- defined continuous turning tangent plane at any point p ∈ S satisfies the geometric condition.

12

Then for any x ∈ Up there is a unique representation of the following form: x = qx + u n ; where qx lies on the plane through p which is orthogonal to

• n and u ∈ R. Define:

Pr(x) = qx .

S

p n x

T (S)

p

q

_

x

p

U

p

Remark 14 n need not be the normal vector to S at p.

13

Moreover, we can compute bounds for C(f) and K(f), for f ≡ Pr. We get: C(f) ≤ cot α + 1 ; (1) K(f) ≤

1

2(cot α)2 + 4

1

2 + 1

2 cot α

3

2 ≤ (cot α + 1) 3 2 . (2) 14

An algorithm for triangulated surfaces∗ is readily produced from the results above:

• Let Np stand for the normal vector to the surface at a

point p on the surface.

• Choose a triangle ∆, of the triangulation. There are two

possibilities to chose ∆: one is in a random manner and the

• ther is based on curvature considerations. Trivially project

∆ onto itself. †

• Suppose ∆′ is a neighbor of ∆ having edges e1, e2, e3,

where e1 is the edge common to both ∆ and ∆′. We will call ∆′ Gehring compatible w.r.t ∆, if the maximal angle between e2 or e3 and N∆ (the normal vector to ∆), is greater then a predefined measure suited to the desired predefined maximal allowed distortion, i.e. max {ϕ1, ϕ2} ≥ α, where ϕ1 = ∡(e2, N∆), ϕ2 = ∡(e3, N∆).

∗that represent the basic tool in Computer Graphics †A variant of this algorithm will be discussed.

15

• Project ∆′ orthogonally onto the plane included in ∆ and

insert it to the patch of ∆, iff it is Gehring compatible w.r.t ∆.

e e e 1

2 3

N

∆

∆ ∆'

N

∆'

• If by this time all triangles where added to the patch

we have completed constructing the mapping. Otherwise, chose a new triangle that has not been projected yet, to be the starting triangle of a new patch.

16

†Variant of the Basic Algorithm:

• Project the faces adjacent to the vertex v on the plane

TPv through v, orthogonal to the mean normal: Nv = 1 k

k

• 1

Ni ; where Ni is the normal to the face Fi adjacent to v.

v Nv _ N N N N N N _ _ _ _ _ _

1 2 3 i j k

F F F F F

1 2 3 i k

F

j

TP

v

17

• Choose starting vertex using Gauss Curvature K:

For triangulated (PL) surfaces we define∗ K at every vertex as the defect of the sum of angles surrounding it: K(v) = δ(v) = 2π −

• i

αi . Remark 15 The curvature based method is better fitted for:

• Low curvature (“almost flat”) surfaces;
• High α.

∗following Descartes and,

in more recent times, Hilbert–Cohn- Vossen, P´

• lya, Banchoff,...

18

We present some experimental results∗, both on synthetic surfaces and on data obtained from actual CT scans (of the Human Brain Cortex and Colon):

∗LNCS 4241 ECCV - CVAMIA 2006, LNCS 4040 IWCIA 2006, EU-

SIPCO 2006

19

A lower curvature produces a larger patch (with more tri- angles)...

In The Skull Model the resolution is of 60,339 triangles. Here α = 10◦ and the dilatation is 1.1763.

20

...than when flattening regions of higher curvature:

Here α = 6◦ and the dilatation is 1.1051.

21

It is also evident in the development of the Human Brain Cortex:

22

23

Remark 16 Note that non-simply connected patches may be obtained.

24

25

To summarize:

• The above given algorithm is local.

Indeed, in a sense the (proposed) algorithm gives a measure of “globality” of this intrinsically local process.

• Our algorithm is best suited for highly folded surfaces,

because of its intrinsic locality, on the one hand, and com- putational simplicity, on the other. It is highly sensitive to high curvature (but not to mesh quality!) – thus being excellently suited for discovering tumors (e.g. in the colon, where they are associated to high curvature)

• However, on “quasi-developable” surfaces (i.e. surfaces

that are almost cylindrical or conical) the algorithm be- haves similar to other algorithms, with practically identical results).

26

• The theory and algorithm guarantee minimal (and com-

putable!) metric, angular and area distortion.

• Relatively simple – yet correct(!), robust and computa-

tionally efficient, since it does not require computational derivatives.

• Holds in any dimension.

27

Towards Globality

Gluing Different Patches to obtain a Global Flattening: The need for gluing patches together into a global picture is well known in Radiography as “pantomograph” This is done very approximatively and with no control of the dilatation.

28

We have applied a “naive” (but with dilatation control) gluing process to the triangulated surface obtained from 3 slices of human colon scan:

(a) (b)

29

The reason for these “cuts” and “holes” resides in the fact that (evidently) one can have two neighbouring patches, with markedly different dilatations/distorsions, which re- sults in different lengths for the common boundary edges. Therefore, “cuts” and “holes” appear when applying a “naive” gluing – as the colon flattening example shows. The discontinuities appear at the common boundary of two patches obtained from regions with very different curvature. Indeed, the “back part” seems close enough to be half of a cylinder (and thus developable)...

30

...but in fact it is highly folded:

31

32

−6.5 −6 −5.5 −5 −4.5 14 14.5 15 15.5 −1.5 −1 −0.5 0.5 1 1.5 −7 −6 −5 −4 −3 −2 −1 14 15 16 17 18 19 −2 2 4 6 −7 −6 −5 −4 −3 −2 −1 1 14 16 18 20 −2 2 4 6 −8 −6 −4 −2 2 14 16 18 20 22 −4 −2 2 4 6 −8 −6 −4 −2 2 10 15 20 25 −3 −2 −1 1 2 3 4 5 −8 −6 −4 −2 2 10 15 20 25 −3 −2 −1 1 2 3 4 5

33

−15 −10 −5 5 10 20 30 −3 −2 −1 1 2 3 4 5 −15 −10 −5 5 5 10 15 20 25 −5 5 −20 −10 10 20 30 40 50 20 40 −5 5 −20 −10 10 20 30 40 50 20 40 −5 5 −20 20 40 60 5 10 15 20 25 30 −5 5 −20 −10 10 20 30 40 50 5 10 15 20 25 30 −5 5

34

Conformal Modulus

• For Colon Unfolding: Use global conformal invariants, e.g.

Ring/ Modulus, which, in the our particular setting can be taken to be: :∗ M(C) = 2πminc(length(c))2 Area(C) ; or M(C) = 2π Area(C) minγ(length(γ))2 . where c and γ are types of curves depicted in the figure below:

C

c γ

∗One can employ also the Quadrilateral Modulus.

35

In this case we used quasi-geodesics and we obtain the following estimate: 15.04 ≤ M(Q) ≤ 21.26 .

36

The quasi-geodesics were computed using an adaptation of the algorithm of Polthier and Schmies for the computation

• f geodesic curvature for curves on polyhedral surfaces:

κg(γ, v) = 2π θ

θ

2 − β

• where β represents a choice (i.e. the smallest) of the angle
• fγ at p where δ represents the discrete straightest geodesic

(at p)

δ γ β θ/2 θ/2−β

p

37

In this discrete approach to geodesic curvature, κg(e) ≡ 0, for any edge e. To include the geodesic curvature of the edges one can apply the the following formula: κg(e) = ∠( ¯ Nv1, ¯ Nv2) where ¯ Nv1, ¯ Nv2 represent the computed normals at the vertices that define e.

v v v v

N N e _ _

1 2 1 2

38

v

k

e v v v

1 n m

δ

δ

β θ/2

e e v

1

vk v

m

vn e

1

δ

e

1

e

1

vk v

m

vn v

1

v

1

γ

39

Another Application: Images

In recent years it became common amongst the image pro- cessing community, to consider images as Riemannian man- ifolds embedded in higher dimensional spaces the embed- ding manifold usually being Rn. For example, a gray scale image is a surface in R3, whereas a color image is a surface embedded in R5, each color channel representing a coordi-

• nate. In both cases the intensity, either gray scale or color,

is considered as a function of the two spatial coordinates (x, y). D be a domain in Rn−1 and let u : D → R1, u ∈ C1 de- note the grey scale (or luminosity) function. Define: S = {f(x) | x ∈ D}, where f(x) = x + u(x)en, x ∈ D, and where en = (0, . . . , 0) ∈ Rn.

40

Then KI(f) = KO(f) = K(f) = 1 2n

• a +
• a2 + 4

n

; where a = sup

x∈D

|u′(x)| . In particular, for the cases of practical interest in (medical) imaging, that is n = 2, n = 3 the maximal dilatation of f is K(f) = 1

4

• a +
• a2 + 4

2

, and K(f) = 1

8

• a +
• a2 + 4

3

, respectively. Note that f is qc iff a < ∞. It also follows immediately that K(f) is minimal for a = 0, that is for u ≡ const. This produces the expected – yet trivial – result that the dilata- tion of f, hence the that of the “photographic surface” is minimal for gray shading.

41

Natural (gray-scale) image (left) and its associated gray-scale surfaces (right): “Lena”. Here the maximal dilatation is K = 1.9753.

42

Natural (gray-scale) image (left) and its associated gray-scale surface (right): a brain scan image. Here the maximal dilatation is K = 2.3682.

43

Further applications

• Edge detection

44

• Face recognition

However, due to the high sensitivity to edges (high curva- ture) and the fact that K is a supremum, it can not be applied to whole faces.

45

However, it functions very well when applied to:

• Eyes
• Mouth

46

Tested on:

• ≈ 40 people-images
• Under different illumination conditions
• With or without glasses

47