Almost Vanishing Polynomials and an Application to Hough Transforms - - PowerPoint PPT Presentation

Almost Vanishing Polynomials and an Application to Hough Transforms

Maria-Laura Torrente

joint work with M.C. Beltrametti

Dipartimento di Matematica Universit` a di Genova

RICAM, Linz November 2013

• M. Torrente

AVP and Applications to HT RICAM, Linz 1 / 29

Introduction

Let’s start from the application.

• M. Torrente

AVP and Applications to HT RICAM, Linz 2 / 29

Introduction

Let’s start from the application. In the analysis of digital images, e.g. medical and astronomical images, the problem of automated recognition of special curves is very important.

• M. Torrente

AVP and Applications to HT RICAM, Linz 2 / 29

The Hough Transform

The main tool is based on the Hough Transform technique. HT is a technique mainly used in image analysis, computer vision, and digital image processing. The purpose of HT is to identify, in a given image, (approximate) instances of a certain class of shapes.

• M. Torrente

AVP and Applications to HT RICAM, Linz 3 / 29

The Hough Transform

The main tool is based on the Hough Transform technique. HT is a technique mainly used in image analysis, computer vision, and digital image processing. The purpose of HT is to identify, in a given image, (approximate) instances of a certain class of shapes. Originally (1962, Hough) HT was concerned with identification of lines in images; later on (1972 Duda & Hart, 1981 Ballard) HT was extended to identify circles and ellipses; many refinements have been investigated since then. HT exploits the duality between image space and parameter space; result is achieved through a voting procedure in the parameter space.

• M. Torrente

AVP and Applications to HT RICAM, Linz 3 / 29

Detecting Aligned Points

Suppose we want to detect aligned points in a given picture.

• M. Torrente

AVP and Applications to HT RICAM, Linz 4 / 29

Detecting Aligned Points

Suppose we want to detect aligned points in a given picture. Represent a straight line as y = ax + b (not the best representation!).

• M. Torrente

AVP and Applications to HT RICAM, Linz 4 / 29

Detecting Aligned Points

Suppose we want to detect aligned points in a given picture. Represent a straight line as y = ax + b (not the best representation!). Consider points in the picture (practically pixels, small cells).

• M. Torrente

AVP and Applications to HT RICAM, Linz 4 / 29

Detecting Aligned Points

Suppose we want to detect aligned points in a given picture. Represent a straight line as y = ax + b (not the best representation!). Consider points in the picture (practically pixels, small cells). Let p1 = (x1, y1) be a point of the picture; a straight line containing it satisfies y1 = ax1 + b.

• M. Torrente

AVP and Applications to HT RICAM, Linz 4 / 29

Detecting Aligned Points

Suppose we want to detect aligned points in a given picture. Represent a straight line as y = ax + b (not the best representation!). Consider points in the picture (practically pixels, small cells). Let p1 = (x1, y1) be a point of the picture; a straight line containing it satisfies y1 = ax1 + b. MAIN IDEA: move to the parameter space , so y1 = ax1 + b is a straight line in this space. Repeat this process for every point p2, p3, . . . in the picture. Let (A, B) the intersection point of many lines, it means that the corresponding points in the picture lie on y = Ax + B!!!

• M. Torrente

AVP and Applications to HT RICAM, Linz 4 / 29

Aligned Points

Image space

✲ ✻

x y O

r ❅ ❅ ❅ ❅

p1 = (x1, y1)

• M. Torrente

AVP and Applications to HT RICAM, Linz 5 / 29

Aligned Points

Image space

✲ ✻

x y O

r ❅ ❅ ❅ ❅

p1 = (x1, y1)

Parameter space

✲ ✻

a b O

❍❍❍❍❍❍

b = −x1a + y1

• M. Torrente

AVP and Applications to HT RICAM, Linz 5 / 29

Aligned Points

Image space

✲ ✻

x y O

r r r ❅ ❅ ❅ ❅

p1 = (x1, y1) p2 p3

Parameter space

✲ ✻

a b O

❍❍❍❍❍❍

b = −x1a + y1

• M. Torrente

AVP and Applications to HT RICAM, Linz 6 / 29

Aligned Points

Image space

✲ ✻

x y O

r r r ❅ ❅ ❅ ❅

p1 = (x1, y1) p2 p3

Parameter space

✲ ✻

a b O

r ❍❍❍❍❍❍

• ✟✟✟✟✟✟

b = −x1a + y1 b = −x2a + y2 b = −x3a + y3 B A

• M. Torrente

AVP and Applications to HT RICAM, Linz 7 / 29

Detecting Aligned Points

Suppose we want to detect aligned points in a given picture. Represent a straight line as y = ax + b (not the best representation!). Consider points in the picture (that is pixels, small cells). Let p1 = (x1, y1) be a point of the picture; a straight line containing it satisfies y1 = ax1 + b. MAIN IDEA: move to the parameter space, so y1 = ax1 + b1 is a straight line in this space. Repeat this process for every point p2, p3, . . . in the picture. Let (A, B) be the intersection point of many lines, it means that the corresponding points in the picture lie on y = Ax + B!!!

• M. Torrente

AVP and Applications to HT RICAM, Linz 8 / 29

Aligned Points

Image space

✲ ✻

x y O

r r r ❅ ❅ ❅ ❅

p1 p2 p3 y = Ax + B

Parameter space

✲ ✻

a b O

r ❍❍❍❍❍❍

• ✟✟✟✟✟✟

b = −x1a + y1 b = −x2a + y2 b = −x3a + y3 B A

• M. Torrente

AVP and Applications to HT RICAM, Linz 9 / 29

The Hough Transform (Formalism)

Our interest is to detect more complicated curves.

• M. Torrente

AVP and Applications to HT RICAM, Linz 10 / 29

The Hough Transform (Formalism)

Our interest is to detect more complicated curves. Let: n, t be positive integers; K = R; x = (x1, . . . , xn) and a = (a1, . . . , at);

• M. Torrente

AVP and Applications to HT RICAM, Linz 10 / 29

The Hough Transform (Formalism)

Our interest is to detect more complicated curves. Let: n, t be positive integers; K = R; x = (x1, . . . , xn) and a = (a1, . . . , at); (Beltrametti, Robbiano 2012) F(x, a) ∈ K[x1, . . . , xn, a1, . . . , at] = K[x, a] such that for each α = (α1, . . . , αt) ∈ At(K) (parameter space) and for each p = (p1, . . . , pn) ∈ An(K) (image space) we have: Hα : F(x, α) := fα(x) = 0

• irreduc. hypersurface degree d

Γp : F(p, a) := Γp(a) = 0 Hough transform of the point p

• M. Torrente

AVP and Applications to HT RICAM, Linz 10 / 29

The Hough Transform (Formalism)

Our interest is to detect more complicated curves. Let: n, t be positive integers; K = R; x = (x1, . . . , xn) and a = (a1, . . . , at); (Beltrametti, Robbiano 2012) F(x, a) ∈ K[x1, . . . , xn, a1, . . . , at] = K[x, a] such that for each α = (α1, . . . , αt) ∈ At(K) (parameter space) and for each p = (p1, . . . , pn) ∈ An(K) (image space) we have: Hα : F(x, α) := fα(x) = 0

• irreduc. hypersurface degree d

Γp : F(p, a) := Γp(a) = 0 Hough transform of the point p Proposition (Regularity Property): the following conditions are equivalent: a) for any Hα, Hα′, we have Hα = Hα′ ⇒ α = α′; b) for any Hα, we have

p∈Hα Γp(a) = {α}.

• M. Torrente

AVP and Applications to HT RICAM, Linz 10 / 29

Detection procedure

Consider the case n = 2 (detection of curves in images). Let F = {Hα} be a suitable (irreducible, with fixed degree...) family of curves. Assume that Regularity Property (RP) holds.

• M. Torrente

AVP and Applications to HT RICAM, Linz 11 / 29

Detection procedure

Consider the case n = 2 (detection of curves in images). Let F = {Hα} be a suitable (irreducible, with fixed degree...) family of curves. Assume that Regularity Property (RP) holds. Recognition Algorithm

1 Choose a set X = {p1, . . . , pν} of points of interest in A2(R). 2 In At(R) find the (unique) intersection of the HT corresponding to

the points pi, that is compute {α} =

i=1,...,ν Γpi(a).

3 Return the parameter α, and the curve Hα uniquely determined by α.

• M. Torrente

AVP and Applications to HT RICAM, Linz 11 / 29

Detection procedure

Consider the case n = 2 (detection of curves in images). Let F = {Hα} be a suitable (irreducible, with fixed degree...) family of curves. Assume that Regularity Property (RP) holds. Recognition Algorithm

1 Choose a set X = {p1, . . . , pν} of points of interest in A2(R). 2 In At(R) find the (unique) intersection of the HT corresponding to

the points pi, that is compute {α} =

i=1,...,ν Γpi(a).

3 Return the parameter α, and the curve Hα uniquely determined by α.

Note that, in practice, {α} =

i=1,...,ν Γpi(a) is NOT CORRECT since:

RP applies to infinite intersection;

• M. Torrente

AVP and Applications to HT RICAM, Linz 11 / 29

Detection procedure

Consider the case n = 2 (detection of curves in images). Let F = {Hα} be a suitable (irreducible, with fixed degree...) family of curves. Assume that Regularity Property (RP) holds. Recognition Algorithm

1 Choose a set X = {p1, . . . , pν} of points of interest in A2(R). 2 In At(R) find the (unique) intersection of the HT corresponding to

the points pi, that is compute {α} =

i=1,...,ν Γpi(a).

3 Return the parameter α, and the curve Hα uniquely determined by α.

Note that, in practice, {α} =

i=1,...,ν Γpi(a) is NOT CORRECT since:

RP applies to infinite intersection; in general

i=1,...,ν Γpi(a) = ∅!!!!

• M. Torrente

AVP and Applications to HT RICAM, Linz 11 / 29

Approximate computation of

i=1,...,ν Γpi(a)

In practice, we want to compute an approximation of

i=1,...,ν Γpi(a), that

is {α} ≈

i=1,...,ν Γpi(a). We proceed as follows:

• M. Torrente

AVP and Applications to HT RICAM, Linz 12 / 29

Approximate computation of

i=1,...,ν Γpi(a)

In practice, we want to compute an approximation of

i=1,...,ν Γpi(a), that

is {α} ≈

i=1,...,ν Γpi(a). We proceed as follows:

1 Consider a bounded region T = [a1, b1] × . . . × [at, bt] ⊂ At(R)

and a discretization step d = (d1, . . . , dt).

2 Construct d-discretization of T ⇔ multi-grid of size J1 × . . . × Jt.

Construct the corresponding multi-matrix A (accumulator function)

• f size J1 × . . . × Jt. Initially A is zero.

3 For each i = 1, . . . , ν and each j = (j1, . . . , jt) assign to A(j)

A(j) = A(j) + 1 if Γpi(a) ∩ C(j) = ∅ A(j) if Γpi(a) ∩ C(j) = ∅ where C(j) denotes the j-th cell of the discretization of T .

4 Find cell C(j∗) such that A(j∗) = maxj A(j); return its center α.

• M. Torrente

AVP and Applications to HT RICAM, Linz 12 / 29

The problem

Problem: Determine (possibly computationally simple) conditions to decide if Γpi(a) ∩ C(j) = ∅;

• M. Torrente

AVP and Applications to HT RICAM, Linz 13 / 29

The problem

Problem: Determine (possibly computationally simple) conditions to decide if Γpi(a) ∩ C(j) = ∅; Notation: let P = R[x] = R[x1, . . . , xn]; let f = f (x) be a polynomial of P; assume f = 0 is a hypersurface;

• M. Torrente

AVP and Applications to HT RICAM, Linz 13 / 29

The problem

Problem: Determine (possibly computationally simple) conditions to decide if Γpi(a) ∩ C(j) = ∅; Notation: let P = R[x] = R[x1, . . . , xn]; let f = f (x) be a polynomial of P; assume f = 0 is a hypersurface; let p = (p1, . . . , pn) be a point of Rn; let ε1, . . . , εn be positive real numbers and ε = (ε1, . . . , εn);

• M. Torrente

AVP and Applications to HT RICAM, Linz 13 / 29

The problem

Problem: Determine (possibly computationally simple) conditions to decide if Γpi(a) ∩ C(j) = ∅; Notation: let P = R[x] = R[x1, . . . , xn]; let f = f (x) be a polynomial of P; assume f = 0 is a hypersurface; let p = (p1, . . . , pn) be a point of Rn; let ε1, . . . , εn be positive real numbers and ε = (ε1, . . . , εn); let B(p) be the (∞, ε)-unit ball centered at p, that is: B(p) = {q ∈ Rn : (q − p)t∞,E ≤ 1} where E = diag(1/ε1, . . . , 1/εn) and v∞,E = Ev∞ =

n

max

i=1 |(Ev)i|,

where v ∈ Rn. B(p) represents the generic cell C(j) of the discretization.

• M. Torrente

AVP and Applications to HT RICAM, Linz 13 / 29

The evaluation of f at p

Obviously, |f (p)| gives useful informations about the crossing.

• M. Torrente

AVP and Applications to HT RICAM, Linz 14 / 29

The evaluation of f at p

Obviously, |f (p)| gives useful informations about the crossing. Nevertheless, the value |f (p)| is not enough since:

|f (p)| small {f = 0} ∩ B(p) = ∅; |f (p)| large {f = 0} ∩ B(p) = ∅.

• M. Torrente

AVP and Applications to HT RICAM, Linz 14 / 29

The evaluation of f at p

Obviously, |f (p)| gives useful informations about the crossing. Nevertheless, the value |f (p)| is not enough since:

|f (p)| small {f = 0} ∩ B(p) = ∅; |f (p)| large {f = 0} ∩ B(p) = ∅.

The crossing criteria will depend on:

the tolerance ε (equivalently the discretization step); the local differential geometry of z = f (x1, . . . , xn) in An+1(R).

• M. Torrente

AVP and Applications to HT RICAM, Linz 14 / 29

Example 1

Let f (x, y) = x2 +

1 100y2 − 1 100 and p = (0, 2).

We have: |f (p)| = 0.03 small dist(p, {f = 0}) = 1

• 2,4
• 2
• 1,6
• 1,2
• 0,8
• 0,4

0,4 0,8 1,2 1,6 2 2,4

• 1,2
• 0,8
• 0,4

0,4 0,8 1,2 1,6 2 2,4

p p

Figure : x2 +

1 100y 2 − 1 100 = 0 and

point p = (0, 2)

• M. Torrente

AVP and Applications to HT RICAM, Linz 15 / 29

Example 1

Let f (x, y) = x2 +

1 100y2 − 1 100 and p = (0, 2).

We have: |f (p)| = 0.03 small dist(p, {f = 0}) = 1

• 2,4
• 2
• 1,6
• 1,2
• 0,8
• 0,4

0,4 0,8 1,2 1,6 2 2,4

• 1,2
• 0,8
• 0,4

0,4 0,8 1,2 1,6 2 2,4

p p

Figure : x2 +

1 100y 2 − 1 100 = 0 and

point p = (0, 2) Figure : z = x2 +

1 100y 2 − 1 100 = 0

and point p′ = (0, 2, 0)

• M. Torrente

AVP and Applications to HT RICAM, Linz 15 / 29

Example 2

Let f (x, y) = y − 10x2 and p = (1.1, 10). We have: |f (p)| = 2.1 big dist(p, {f = 0}) ≈ 0.1

• 7,5
• 5
• 2,5

2,5 5 7,5 10 2,5 5 7,5 10

p

Figure : y − 10x2 = 0 and point p = (1.1, 10)

• M. Torrente

AVP and Applications to HT RICAM, Linz 16 / 29

Example 2

Let f (x, y) = y − 10x2 and p = (1.1, 10). We have: |f (p)| = 2.1 big dist(p, {f = 0}) ≈ 0.1

• 7,5
• 5
• 2,5

2,5 5 7,5 10 2,5 5 7,5 10

p

Figure : y − 10x2 = 0 and point p = (1.1, 10) Figure : z = y − 10x2 = 0 and point p′ = (1.1, 10, 0)

• M. Torrente

AVP and Applications to HT RICAM, Linz 16 / 29

Necessary crossing cell conditions

• Prop. 1. Notation as above.

Let H = maxx∈B(p) Hf (x)∞ and εmax = ε∞. If |f (p)| > Jacf (p)1εmax + H 2 ε2

max := B1

then the hypersurface of equation f = 0 does not cross B(p).

• M. Torrente

AVP and Applications to HT RICAM, Linz 17 / 29

Necessary crossing cell conditions

• Prop. 1. Notation as above.

Let H = maxx∈B(p) Hf (x)∞ and εmax = ε∞. If |f (p)| > Jacf (p)1εmax + H 2 ε2

max := B1

then the hypersurface of equation f = 0 does not cross B(p).

• Prop. 2. Notation as above. If

|f (p)| > Jacf (p)1εmax + 1 2Hf (p)∞ε2

max := B′ 1

then the hypersurface of equation f = 0 does not cross B(p) (neglecting contributions of order O(ε3

max)).

• M. Torrente

AVP and Applications to HT RICAM, Linz 17 / 29

Sufficient crossing cell conditions

• Prop. 3. Notation as above.

Suppose that Jacf (p) is not the zero vector. Let 0 < R < min{εmin, Jacf (p)1

H

}, J = supx∈B(p,R) Jac†

f (x)∞ and

c = max{2, √n}. If |f (p)| < 2R J(c + √nHJR) := B2 then the hypersurface of equation f = 0 crosses B(p).

• M. Torrente

AVP and Applications to HT RICAM, Linz 18 / 29

Sufficient crossing cell conditions

• Prop. 3. Notation as above.

Suppose that Jacf (p) is not the zero vector. Let 0 < R < min{εmin, Jacf (p)1

H

}, J = supx∈B(p,R) Jac†

f (x)∞ and

c = max{2, √n}. If |f (p)| < 2R J(c + √nHJR) := B2 then the hypersurface of equation f = 0 crosses B(p).

• Prop. 4. Assumptions as in Prop. 3. Let 0 < R < min{εmin, Jacf (p)1

n2Hf (p)∞ }

and Θ = Jac†

f (p)∞ + n2(1 + 2√n) Hf (p)∞ Jacf (p)2

1 R. If

|f (p)| < 2R Θ(c + n5/2Hf (p)∞ΘR) := B′

2

then the hypersurface of equation f = 0 crosses B(p) (neglecting contributions of order O(R3)).

• M. Torrente

AVP and Applications to HT RICAM, Linz 18 / 29

Comparison among the bounds

• Prop. 5. Notation as above. Let R be a real number s.t.

0 < R < min

• εmin, Jacf (p)1

H , Jacf (p)1 n2Hf (p)∞

• Then:

B′

2 + O(R3) ≤ B2 < B′ 1 ≤ B1

• M. Torrente

AVP and Applications to HT RICAM, Linz 19 / 29

Comparison among the bounds

• Prop. 5. Notation as above. Let R be a real number s.t.

0 < R < min

• εmin, Jacf (p)1

H , Jacf (p)1 n2Hf (p)∞

• Then:

B′

2 + O(R3) ≤ B2 < B′ 1 ≤ B1

We observe that:

• Prop. 5. yields B2 < B1, so it may happen that |f (p)| ∈ (B1, B2).

In this cases, using the previous results, nothing can be concluded regarding the crossing problem. Because of the local nature of the results, a more accurate analysis, by considering smaller subregions of B(p), may overcome the problem.

• M. Torrente

AVP and Applications to HT RICAM, Linz 19 / 29

Algorithms

The CROSSING CELL algorithm Input: polynomial f , point p s.t. Jacf (p) = 0, tolerance ε = (ε1, . . . , εn). Output: an element of {0, 1, ξ}.

1 Compute |f (p)|, and the bounds B1 and B2. 2 If |f (p)| > B1 return 0;

if |f (p)| > B2 return 1; else return ξ.

• M. Torrente

AVP and Applications to HT RICAM, Linz 20 / 29

Algorithms II

The CROSSING AREA algorithm Input: polynomial f , region T = [a1, b1] × . . . × [an, bn], discretization step d = (d1, . . . , dn). Output: multi-matrix A with values in {0, 1, ξ}.

1 Construct the discretization of T and the multi-matrix A. 2 For each multi-index j set A(j) = CROSSING CELL(f , xj, d

2 ).

3 Return A.

• M. Torrente

AVP and Applications to HT RICAM, Linz 21 / 29

Algorithms III

The RECOGNITION algorithm Input: regular family F of irreducible curves with same degree, set X = {p1, . . . , pν} of points of interest, region T = [a1, b1] × . . . × [an, bn], discretization step d = (d1, . . . , dn). Output: point x∗.

1 For each i = 1, . . . , ν, set Ai =CROSSING AREA(Γpi(F), T , d). 2 Compute A =

i=1,...,ν Ai.

3 Compute j∗ = argmaxjA(j); return the point x∗ = x(j∗).

• M. Torrente

AVP and Applications to HT RICAM, Linz 22 / 29

Example I

In this (toy) example we aim to detect the external profile of vertebral column using a rational cubic curve (Conchoid of Sl¨ use).

• 1
• 0,75
• 0,5
• 0,25

0,25 0,5 0,75 1 1,25

• 0,75
• 0,5
• 0,25

0,25 0,5 0,75

Figure : Points of interest

• M. Torrente

AVP and Applications to HT RICAM, Linz 23 / 29

Example II

The Conchoid of Sl¨ use of parameters a, b has the following equation: Ca,b : a(x − a)(x2 + y2) = b2x2 For each p = (xp, yp) in the image space, the HT is a conic in the parameter space A2

a,b(R) whose equation is:

Γp(a, b) : (x2

p + y2 p)a2 + x2 pb2 − xp(x2 p + y2 p)a = 0

Γp(a, b) are represented in the following figure.

• M. Torrente

AVP and Applications to HT RICAM, Linz 24 / 29

Example III

0,025 0,05 0,075 0,1 0,125 0,2 0,225 0,25 0,275 0,3

• M. Torrente

AVP and Applications to HT RICAM, Linz 25 / 29

Example IV

We consider the region T = [0.03, 0.5] × [0.05, 0.4] ⊂ A2

(a,b)(R), and

a discretization step d = 2ε = (0.008, 0.008). Accumulator matrix A ∈ Mat47×45 with values in {0, 1} The maximum value of A is 17 and corresponds to the cell centered in (A, B) = (0.078, 0.266)

• . . .

7 12 14 8 2 1 2 . . . . . . 8 14 11 5 3 2 . . . . . . 2 10 17 10 1 . . . . . . 2 2 8 15 9 2 . . . . . . 2 2 2 2 5 13 11 3 . . .

• M. Torrente

AVP and Applications to HT RICAM, Linz 26 / 29

Example V

The blue cell is the one which realizes the maximum for A.

0,025 0,05 0,075 0,1 0,125 0,2 0,225 0,25 0,275 0,3

• M. Torrente

AVP and Applications to HT RICAM, Linz 27 / 29

Example VI

Finally, the detected curve is:

• 1
• 0,75
• 0,5
• 0,25

0,25 0,5 0,75 1 1,25

• 0,75
• 0,5
• 0,25

0,25 0,5 0,75

• M. Torrente

AVP and Applications to HT RICAM, Linz 28 / 29

References

[1] M. C. Beltrametti, A. M. Massone and M. Piana, Hough transform of special classes of curves, SIAM J. Imaging Sci. 6(1), (2013), 391–412. [2] M. C. Beltrametti and L. Robbiano, An algebraic approach to Hough transforms, Journal of Algebra 371 (2012), 669–681. [3] CoCoATeam, CoCo A: a system for doing Computations in Commutative Algebra, available at http://cocoa.dima.unige.it. [4] R. O. Duda and P. E. Hart, Use of the Hough transformation to detect lines and curves in pictures, Comm. ACM, 15, vol. 1 (1972), 11–15. [5] C. Fassino and M. Torrente, Simple Varieties for Limited Precision Points, Theoret.

• Comput. Sci. 479 (2013), 174–186.

[6] P. V. C. Hough, Method and means for recognizing complex patterns, US Patent 3069654, December 18, 1962. [7] A. M. Massone, A. Perasso, C. Campi and M. C. Beltrametti, Profile detection in medical and astronomical imaging by mean of the Hough transform of special classes of curves, preprint, 2013.

Thank you!

• M. Torrente

AVP and Applications to HT RICAM, Linz 29 / 29