Eurographics 2 0 1 3 Projective Geom etry, Duality and Precision of - - PowerPoint PPT Presentation

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 1

Projective Geom etry, Duality and Precision of Com putation in Com puter Graphics, Visualization and Gam es

Tutorial

Vaclav Skala

University of West Bohemia, Plzen, Czech Republic VSB-Technical University, Ostrava, Czech Republic http://www.VaclavSkala.eu

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 2

Plzen ( Pilsen) City

Plzen is an old city [first records of Plzen castle 976] city of culture, industry, and brewery. City, where today’s beer fermentation process was invented that is why today’s beers are called Pilsner - world wide

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 3

Ostrava City

Ostrava is

• an industrial city of coal mining &

iron melting

• 3rd largest city

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 4

University of W est Bohem ia 1 7 5 3 0 students + 9 8 7 PhD students Com puter Science and Engineering Mathem atics (+ Geomatics) Physics Cybernetics Mechanics (Computational)

• Over 5 0 % of income from research and application projects
• New research center EU project - investment 64 mil. EUR
• 2nd in the ranking of Czech technical / informatics faculties 2009, 2012

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 5

“Real science” in the XXI century

Courtesy of Czech Film, Barrandov

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 6

An overview

• Precision and robustness
• Euclidean space and projective extension
• Principle of duality and its applications
• Geometric computation in the projective space
• Intersection of two planes in E3 with additional constrains
• Barycentric coordinates and intersections
• Interpolation and intersection algorithms
• Implementation aspects and GPU
• Conclusion and summary

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 7

Num erical system s

• Binary system is used nearly exclusively
• Octal & hexadecimal representation is used
• If we would be direct descendants of tetrapods –

we would have a great advantage – “simple counting in hexadecimal system”

Nam e Base Digits E m in E m ax BI NARY B 16 Half 2 10+1 −14 15 B 32 Single 2 23+1 −126 127 B 64 Double 2 52+1 −1022 1023 B 128 Quad 2 112+1 −16382 16383 DECI MAL D 32 10 7 −95 96 D 64 10 16 −383 384 D 128 10 34 −6143 6144

IEEE 758-2008 standard

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 8

Mathem atically perfect algorithm s fail due to instability Main issues

• stability, robustness of algorithms
• acceptable speed
• linear speedup – results depends on HW, CPU …. parameters !

Num erical stability

• limited precision of float / double
• tests A ? B with floats

if A = B then ….. else ….. ; if A = 0 then ….. else …. should be forbidden in programming languages

• division operation should be removed or postponed to the last moment if

possible - “blue screens”, system resets

E

T T d w { } /

Eurographics 2

Typical e

• inters
• point
• detec

throu Typical p double x w hile ( x { if (x x += } /* if p

2013

exam ple section o lies on a Ax + ction if a ugh problem x = -1; d x < +1) == p) Co p; p = 0.1 th

Eu

Vaclav Sk

es of ins

• f 2 lines

a line in E By + C line inte m double p

• nsole.O

hen no o

urogr

kala http

stability in E3 E2 or a p = 0 or rsects a = ….; ut.WriteL

• utput,

raphic

p://www.Vacla

plane in E r Ax + By polygon, Line(” ** if p = 0.

cs 2 0

avSkala.eu 2

E3 y + Cz + , touches ** ”) 25 then

1 3

2013-05-03

D = 0 s a vertex expecte x or pass ed output ses t */

9

E

D P i I i p ?

Eurographics 2

Delauna Point insi n regular It can be s very se position c ?? ROBU

2013

y triang de of a c r rectang seen tha ensitive t change USTNESS

Eu

Vaclav Sk

gulation circle give gular grid at the DT to a point S ??

urogr

kala http

& Voron en by thr d. T & t

raphic

p://www.Vacla

noi diag ree point

cs 2 0

avSkala.eu 2

gram ts – prob

1 3

2013-05-03

lems with meshin ng points V

10

s VD

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 11

Courtesy of http://mathworld.wolfram.com/ ConfocalEllipsoidalCoordinates.html

Vectors and Points in Geom etry

• Vectors – movable, no fixed position
• Points – no size, position fixed in the

GIVEN coordinate system Coordinate system s:

• Cartesian – left / right handed

right handed system is used

• Polar
• Spherical, Cylindrical
• and many others, e.g. Confocal Ellipsoidal Coordinates

http://mathworld.wolfram.com/ConfocalEllipsoidalCoordinates.html

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 12

Floating point

• Not all numbers are represented

correctly

• Logarithmic arithmetic
• Continuous fractions
• Interval arithmetic

Num erically NOT valid identities due to limited precision

• 1
• x y

4 1 1 3 2 5 3 …

x + y = [a + c, b + d] x = [ a , b ] x - y = [a - d, b - c] y = [ c , d ] x × y = [min(ac, ad, bc, bd), max(ac, ad, bc, bd)] x / y = [min(a/c, a/d, b/c, b/d), max(a/c, a/d, b/c, b/d)] if y ≠ 0

3; 7,15,1,292,1,1,1,2,1,3,1 …

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 13

Statem ents like if <float> = <float> then …. or if <float> ≠ <float> then …. should not be allow ed in program m ing languages Quadratic equation - m ore reliable results usually solved as , √ 4 2 If 4 then 4 /2

• ⁄

The discriminant should be computed with a twice precision Vieta’s form ula

• ⁄

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 14

Function value com putation at 77617, 33096 , 333.75 11 121 2 5.5 /2 6.33835 10 single precision 1,1726039400532 double precision 1,1726039400531786318588349045201838 extended precision The correct result is “somewhere” in the interval of 0,82739605994682136814116509547981629, 0,82739605994682136814116509547981629 Exact solution , 2 2 54767 66192

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 15

Sum m ation is one of often used com putations. 10 0.999990701675415

• r

10 1.000053524971008

• The result should be only one.

The correctness in sum m ation is very im portant in pow er series com putations. !!!!ORDER of sum m ation 1 14.357357

• 1

14.392651

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 16

Recursion Towers of Hanoi Ackermann function The value of the function grows very fast as 4,4 2 2 MOVE (A, C, n); { MOVE (A, B, n-1); MOVE (A, C, 1); MOVE (B, C, n-1) } # MOVE (from, to, number) #

, 1 1,1 0 0 1, , 1 0 0

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 17

Num erical com putations Hilbert’s Matrix I nversion

• 1 1 1

1 2 1

• 1,0E-13

1,0E-11 1,0E-09 1,0E-07 1,0E-05 1,0E-03 1,0E-01 1,0E+01 1,0E+03 1,0E+05 1,0E+07 1,0E+09 1,0E+11 1,0E+13 1,0E+15 1,0E+17 1,0E+19 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Error

Order of the Hilbert matrix

ε εp ξ

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 18

Mathem atical “form s” There are several “form s”: I m plicit , , 0 or 0 or (system of equations) There is no orientation, e.g.

• if 0 is a iso-curve there is no hint how to find another point of

this curve, resp. a line segment approximating the curve => tracing algorithms

• if 0 is a iso-surface there is no hint how to find another point of

this surface => iso-surface extraction algorithms Param etrical , Points of a curve are “ORDERED” according to a parameter , resp. , Explicit , [actually 2 ½ D] For the given value , resp. , we get function value

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 19

I m plicit form

• Is used for separation - for detection if a point is inside or outside, e.g. a

half-plane or a circle etc.

• There is always a question how to compute of 0 as there are

several solutions in general

• Complexity of computations × precision of computation

Compiler optimization is DANGEROUS in general - numerical precision

• 1
• 1
• 1
• 1

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 20

Another exam ple 1 cos

• 0 0,5

Computed values are w rong in an interval close to zero. In the interval , the function values are ZERO instead of 0.5!!! 1 cos

• Interval 10 , 10

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

• 10
• 5

5 10 (1-cos(x))/x^ 2 x

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.0000000

• 0.0000000
• 0.0000000
• 0.0000000
• 0.0000000

0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000 (1-cos(x))/ x^ 2 x

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 21

Exam ples – w hat happened? There are famous examples of numerical disasters. When reading the original reports and followed comments and details one must be really surprised how simple errors occur and should be worried what could happen in complex problems solution. Let us shortly explore some “traditional” cases.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 22

Courtesy CNN

Explosion of Ariane 5 An Ariane 5 rocket was launched by the European Space Agency (ESA) on June 4,

• 1996. The development cost over

$7 billion. The rocket exploded after lift-off in about 40 sec. Destroyed rocket and cargo were valued at $500 million. The cause of a failure was a software error in inertial reference system. From the CNN article: “The internal SRI [ Inertial Reference System] software exception was caused during execution of a data conversion from 64-bit floating point to 16-bit signed integer value. The floating point number which was converted had a value greater than what could be represented by a 16-bit signed integer.” The conversion from the floating point to the integer representation is very dangerous as it is not reported by an exception and stored value represents an existing number.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 23

Courtesy of GAO report

Patriot Missile Failure

The system was originally designed in mid- 1960 for a short and flexible operation. The system was actually running for more than 100 hours) and for intercepting cruise missiles running at MACH 2 speed and was used to intercept the Scud missile running at MACH 5. The computation of intercepting and hitting was based on time counting with 24 bits integers with the clock of 1/10 and speed computation in floats. The clock setting to 1/10 was a critical issue and not acceptable even for application in sport activities at that time. Unfortunately 1/10 = 1/24+1/25+1/28+1/29+1/212+.... and therefore the error on 24 bits is about 0.000000095 and in 100 hours the error is 0.34. As the Scud flies at MACH 5, the error was actually 687[m] and the missile was out of the “range gate” area. As a result of the fault assumptions, incorrect software design and irresponsible attitude of the army officials, 28 Americans were killed and over 100 other people injured in the Iraq’s Scud missile attack in Dhahran, Saudi Arabia on February 25, 1991 according to the GAO report.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 24

Courtesy of SINTEF

Sleipner offshore platform sinking

Another well known example is the Sleipner

• ffshore platform sinking. The top deck is

about 57 000 tons, drilling and support equipments weight about 40 000 tons and the deck provides an accommodation for about 200 people. The Sleipner platform structure was “optimized” using finite element system and the shear stresses were underestimated nearly by 50%. It led to serious cracks in the structure and leakage that the pumps were unable to cope with. The sinking of the platform estimated cost is about $700 million.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 25

We have presented some basic facts on numerical precision and examples of some disasters. Many engineering problems are somehow connected with geometry and geometrical computations with respecting physical phenomena etc. The majority of computations are made in the Euclidean space representation and with the Cartesian coordinate system. In the following we will show how the non-Euclidean representation, actually its projective extension, and the principle of duality can be used to solve some problems in a simple, robust and elegant ways.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 26

x y w w=1 x X Y (a) p P

2

E

2

ρ a b c c=1 D(p) D( ) ρ A B (b) D(P )

2

D(E )

2

Projective Space , , : Conversion: /, / 0 If 0 then represents “an ideal point” - a point in infinity, i.e. it is a directional vector. The Euclidean space E2 is represented as a plane 1. Equivalent “mathematical” notation often used: : , generally for En : , … , i.e. homogeneous coordinate is the first

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 27

Points and vectors

• Vectors are “freely m ovable” – not having a fixed position

, : 0

• Points are not “freely m ovable” – they are fixed to an origin of the

current coordinate system , : and , : usually in textbooks 1 A vector in the Euclidean coordinate system – CORRECT ,

, ,

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 28

Horrible “construction”! DO NOT USE I T – I T I S TOTALLY W RONG , : as 1 , : 1 1 , : 0 What happen if due to numerical representation? , : , :

• and 0

This is considered as a point !!! This was presented as “How a vector” is constructed in the projective space in a textbook!! W RONG, W RONG, W RONG This construction has been found in SW!!

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 29

A Euclidean vector given by two points expressed in

• the Euclidean coordinates
• ,
• the homogeneous coordinates

, : We use homogeneous coordinates to represent a denominator of a fraction This is the CORRECT SOLUTI ON, but what is the interpretation? A “difference” of coordinates of tw o points is a vector in the m athem atical m eaning and is a “scaling” factor actually Actually the division operation is postponed and not performed immediately. A vector in projective notation , , : 0

• where: means projectively equivalent

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 30

A Euclidean vector in the projective representation (if the vector length matters) , , : 0

• A vector in the projective space is given by coordinates , , as

, : [=>Linear interpolation with a non-linear monotonic interpolation] We have to strictly distinguish meaning of one dimensional array [vector], i.e. if we are working with:

• points, i.e. a data structure represent point coordinates, or
• vectors, i.e. a data structure represent a vector in the mathematical

meaning VECTORS x POI NTS

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 31

x y w w=1 x X Y (a) p P

2

E

2

ρ a b c c=1 D(p) D( ) ρ A B (b) D(P )

2

D(E )

2

Duality For simplicity, let us consider a line defined as: We can multiply it by 0 and we get: As and 0 i.e. , : , : , : A line is actually a plane in the projective space (point 0,0: 0excluded)

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 32

Duality From the mathematical notation we cannot distinguish whether is a line and is a point or vice versa in the case of . It means that

• a point and a line are dual in the case of , and
• a point and a plane are dual in the case of .

The principle of duality in states that: Any theorem in E2 remains true when we interchange the words “point” and “line”, “lie on” and “pass through”, “join” and “intersection”, “collinear” and “concurrent” and so on. Similarly for the E3 case. Once the theorem has been established, the dual theorem is obtained as described above. This helps a lot to solve some geometrical problems.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 33

Exam ples of dual objects and operators Prim itive Dual prim itive

• Point

Line Line Point

• Point

Plane Plane Point Operator Dual operator Join Intersect Intersect Join Com putational sequence for a problem is the sam e as for a dual problem .

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 34

Definition The cross product of the two vectors , : and , : is defined as:

• where: 1,0: 0 0,1: 0 0,0: 1
• r as

1 1 1 1 1 1 2 2 2 Please, note that homogeneous coordinates are used.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 35

I ntersection of tw o lines Let two lines and are given by , : and , : We have to solve a system of linear equations i.e.

• and

* Then well known formula is used

• But what if is small? What is ?

Usually a sequence like det. . is used.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 36

Theorem Let two lines and be given. Then the coordinates of an intersection point , which is defined by those two lines, are determined as the cross product of homogeneous coefficients of those lines as , : Proof We are actually looking for a solution to the following equations:

0 0

where: , : Note * usually a line is in its implicit form as instead of 0, or in the

explicit form as .

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 37

A Line given by tw o points Given two points and and we want to compute a line given by those two points, i.e. we need to compute 3 values , , from two values , . One parametric set of solutions In a matrix form

• 1
• 1
• How to solve it?

Select 1 ? W hat happen if a line passing the origin?

• r 1 ? or 1 or similarly? NO, NO, NO!

BUT HOW ?

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 38

• 1
• 1
• Additional condition 1?
• 1
• 1

1 1

• 1
• Another approach

We know that a line is dual to a point in and vice versa. Due to the duality principle in : <= duality => <= w hy different? =>

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 39

Theorem Let two points and be given in the projective space. Then the coefficients of the line, which is defined by those two points, are determined as the cross product of their homogeneous coordinates

, :

Proof Let the line be defined in homogeneous coordinates as We are actually looking for a solution to the following equations:

0 0

where: , : Note that represents a “distance” from the origin of the coordinate system.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 40

It means that any point that lies on the line must satisfy both the equation above and the equation 0 in other words the vector is defined as

• We can write
• Note that the cross product and the dot product are the instructions in

Cg/HLSL on GPU.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 41

Evaluating the determinant

• we get the line coefficients of the line as:
• Note:

1.A line 0 is a one parametric set of coefficients , : From two values and we have to compute 3 values, coefficients , and 2.For 1 we get the standard cross product formula and the cross product defines the line, i.e. 1 2 where: , :

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 42

DUALI TY APPLI CATI ON In the projective space points and lines are dual. Due to duality we can directly intersection of two lines as

• , :

If the lines are parallel or close to parallel, the homogeneous coordinate 0 and users have to take a decision – so there is no sequence in the code like det. . …in the procedure. Generally computation can continue even if 0 if projective space is used.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 43

Com putation in Projective Space

• Extended cross product definition
• A plane is determined as a cross

product of three given points Due to the duality

• An intersection point of three

planes is determined as a cross product of three given planes.

• Computation of generalized cross product is equivalent to a solution of a

linear system of equations = > no division operation!

• Using the cross product we can continue with symbolic operations which

could not be made if solution of is used.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 44

We have seen that computation of

• an intersection of two lines is given as
• a line given by two points is given as

Different schem e BUT Those problem s are DUAL. W hy algorithm s should be different??

Cross product is equivalent to a solution of both linear system s of equations, i.e. and ! No division operations!

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 45

DI STANCE Geometry is strongly connected with distances and their measurement. Geometry education deals strictly with the Euclidean geometry, where the distance is measured as ∆ ∆ , resp. ∆ ∆ ∆ . This concept is convenient for a solution of basic geometric problems, but in many cases it results into quite complicated formula. There are severe questions of stability and robustness in many cases. The m ain objection against the projective representation is that there is no m etric.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 46

The distance of two points can be easily computed as / where: Also a distance of a point from a line in E2 can be computed as

• √

where: , : , : The extension to E3/P3 is simple and the distance of a point from a plane in E3 can be computed as

• √

where: , , : , , : .

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 47

In many cases we do not need actually a distance, e.g. for a decision which

• bject is closer, and distance2 can be used instead, i.e. for the E2 case
• where: , : : and the normal vector is not normalized.

If we are comparing distances of points from the given line p we can use “pseudo-distance” for comparisons _

• Similarly for a plane in the case of E3
• and

_

• where: , , : :

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 48

Geom etric transform ations w ith points (note ⁄ , , 0 ): Translation by a vector , , : , i.e. / , / , 0: In the Euclidean space:

• 1
• 1
• 1
• /

/ 1 1

• In the projective space:
• /

/ 1

• ⁄

⁄

• 1
• 1
• and det . For 1 we get a standard formula
• 1

1

• 1
• 1
• 1

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 49

Rotation by an angle ,

• ,
• , :

In the Euclidean space:

• 1
• 1
• 1
• In the projective space:
• /

/ 1

• 1
• 1
• as by definition, det

As we are working in the projective space, it does not matter

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 50

Scaling by a factor ,

• ,
• , :
• 1
• det

It is necessary to note that the determinant of a transformation matrix , i.e. matrices , , , is det 1 in general, but as the formulation is in the projective space, there is no need to “normalize” transformations to det 1 even for rotation. It can be seen that if the parameters of a geometric transformation are given in the homogeneous coordinates, no division operation is needed at all.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 51

Transform ation of lines and planes E2 E3 Dual problem In graphical applications position of points are changed by an interaction, i.e. The question is how coefficients of a line, resp. a plane are changed if the points are transformed without a need to be recomputed from the definition. It can be proved that

• r

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 52

Transform ation of lines and planes As the computation is made in the projective space we can write , : for lines in E2

• r

, , : for planes in E3 THI S SI MPLI FI ES COMPUTATI ONS Transformation matrices for lines, resp. for planes are DI FFERENT from transformations for points! Note that a normal vector of a line is actually a co-vector, i.e. an oriented “surface”.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 53

Transform ation of lines and planes Transformation about a general axis in E3 / P3 Usually used transformation ( is translation):

• 1
• √
• √
• √
• √

1

• , , is an axis directional vector. This is unstable if √ 0 and

not precise if or vice versa. That is generally computationally complex and unstable as a user has to select which axis is to be used for a rotation ϕ Rotation axis

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 54

Transform ation of lines and planes Transformation about an axis in the Euclidean space E3 1 . 1 where: . is a matrix. In the Euclidean space E3 the vector has to be normalized The matrix is defined as:

• in our case

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 55

Projection Standard projection

• 1
• How to determine

coordinate? 1 1 1

• 1 1
• 1
• “Reverse” depth com putation, e.g. for correct intensity computation

if perspective projection is used. xP xQ xI zP zQ zI x z P I Q 1

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 56

1,90 2,10 2,30 2,50 2,70 2,90 3,10 0,00 2,00 4,00 6,00 X - Y 0,60 1,60 2,60 3,60 4,60 5,60 0,0 0,2 0,4 0,6 0,8 1,0

X(t) Y(t)

Com putation in Projective Space Linear interpolation Linear parameterization: ∞, ∞ Non-linear monotonous parameterization: ∞, ∞

• we can interpolate using homogeneous

coordinates without “normalization” to Ek !!

• homogeneous coordinate 0

I n m any algorithm s, w e need “m onotonous” param eterization, only

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 57

Com putation in Projective Space Spherical interpolation , , sin1 Ω sin Ω sinΩ sin Ω Instability occurs if Ω . Mathematically formula is correct; in practice the code is generally incorrect! [

• ]

, , , , sin1 Ω sinΩ sin Ω

• sin1 Ω sinΩ : sin Ω projective scalar used

Homogeneous coordinates = > better numerical stability & division operation can be postponed Homogeneous coordinate

Courtesy of Wikipedia

Courtesy of wikipedia

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 58

Com putation in Projective Space , , sin1 Ω sin Ω sinΩ sin Ω sin1 Ω sinΩ : sin ΩT What is a result in the Euclidean space of a spherical interpolation with non-linear parameterization, i.e. 1? sin1 Ω sin Ω sinΩ sin Ω sin1 Ω sin Ω sinΩ sin Ω sin1 Ω sin Ω sinΩ sin Ω sin1 Ω sin Ω sinΩ sin Ω sin1 Ω sinΩ sin1 Ω sinΩ If represented as “projective scalar” value sin1 Ω sinΩ : sin1 Ω sinΩ => Better numerical stability

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 59

I nfluence of hom ogeneous values 1 & 0 (Fun1: 1) Can be used for:

• relative position

comparison

• interpolation if a

“linear” parameterization is required 1

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 60

Com putation in Projective Space Intersection line – plane (Cyrus-Beck clipping algorithm) Computation in projective space! , , : ∞, ∞ 0 0 , , : , , : Euclidean solution:

• Projective solution:

projective scalar: : and : if 0 then Test: if then ...modification: 0 if then .....

• An intersection of a plane with a line in E2 / E3 can be computed efficiently,

no division operation, but com parison operations m ust be m odified!

• Cyrus-Beck line clipping algorithm 10-25% faster

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 61

Line Clipping – convex polygon procedure CLIP_Line ( x A , x B ); /* x A=[xA,yA:wA]T x B=[xB,yB:wB]T */ begin /* p=[a,b:c]T given - NO STEP 1 */ {1} p : = x A x x B; /* p: ax+by+c = 0 */ {2} for k:=0 to N-1 do /*x k=[xk,yk,wk]T*/ {3} if pTx k ≥ 0 then ck:=1 else ck:=0; {4} if c=[0…0]T or c=[1…1]T then EXI T; {5} i:= TAB1[c]; j:= TAB2[c]; {6} x A := p x e i ; x B := p x e j ; {7} DRAW (x A; x B ) * e i – i-th edge */ end /* CLIP_Line */ /* c identifies an edge intersected */ TOO COMPLEX? NO SI MPLE, ROBUST and FAST

• Skala,V.: A new approach to line and line segment clipping in homogeneous coordinates, The Visual Computer,

SpringerVol.21, No.11, pp.905-914, 2005

Line clipping algorithms in E2

− Cohen-Sutherland − Liang-Barsky − Hodgman − Skala – modification of

Clip_L for line segments

x0 xA x2 x3 e3 e0 e1 e2 x1 xB F(x)<0 F(x)>0 s

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 62

For a rectangular normalized window x A := p x e i For the edge 1, i.e. 1 0 , :

• 1

1 , :

• ,

: 1

• , 1: 1

Actually expression for , resp. for is given by the window edge. No m ultiplication or division operations A simple modification if a line is given parametrically (in the Euclidean or projective space) as Simple modification for non-convex polygon but it requires intersections sorting =>

x0 xA x2 x3 e3 e0 e1 e2 x1 xB F(x)<0 F(x)>0 s

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 63

/* Additional code for Line Segment clipping colored */ function CODE (x); /* Cohen Sutherland window coding */ begin c:= [0000]; /* land / lor – bitwise operations and / or */ if x<xmin then c:=[1000] else if x>xmax then c:= [0100]; if y<ymin then c:=c lor [1001] else if y>ymax then c:=c lor [0010]; CODE := c end /*CODE */; procedure Clip_LineSegment (x A , x B); /* short x long line segments */ begin cA := CODE (x A); cB := CODE (x B); if (cA lor cB) = 0 then { output (x A; x B ); EXI T } /* inside */ if (cA land cB) ≠ 0 then EXI T; /* outside */ p := x A x x B; /* ax+by+c = 0; p = [a,b:c]T */ for k:=0 to 3 do /* x k=[xk,yk: wk]T c=[ c3, c2, c1, c0 ]T */ if pTx k ≥ 0 then ck:=1 else ck:=0; if c = [0000]T or c = [1111]T then EXI T; i:= TAB1[c]; j:= TAB2[c]; /* original code x A := p x e i ; x B := p x e j ; DRAW (x A; x B ) */

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 64

/* additional code for Line Segment clipping */ if cA ≠ 0 and cB ≠ 0 then { x A := p x e i ; x B:= p x e j } /* two intersections */ else /* only one intersection point */ if cA = 0 then /* x B is outside */ { if cB land MASK[c] ≠ 0 then x B := p x e i else x B := p x e j } else if cB =0 then /* x A is outside */ { if cA land MASK[c] ≠ 0 then x A := p x e i else x A := p x e j };

• utput (x A , x B )

end /* Clip_LineSegment */

Algorithm can be extended to a convex polygon clipping & modified for parametric lines/rays as well For a convex -sided convex polygon the table can be generated synthetically. Modification for non-convex polygon possible.

c

c3 c2 c1 c0 TAB1 TAB2 MASK None None None 1 1 3 0100 2 1 1 0100 3 1 1 1 3 0010 4 1 1 2 0010 5 1 1 N/A N/A N / A 6 1 1 2 0100 7 1 1 1 2 3 1000 8 1 2 3 1000 9 1 1 2 0100 10 1 1 N/A N/A N / A 11 1 1 1 1 2 0010 12 1 1 1 3 0010 13 1 1 1 1 0100 14 1 1 1 3 0100 15 1 1 1 1 None None None

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 65

P1 x1 x x3 P3 x2 P2

Com putation in Projective Space - Barycentric coordinates Let us consider a triangle with vertices , , A position of any point ∈ can be expressed as

• Additional condition:

1 0 1 1, … , 3

• 1

1 1

• 1
• 1, . . . ,3
• 1

1 1 1

• Linear system has to be solved

If points are given as , , : and 1 then have to be “normalized” to 1, i.e. 4 * 3 = 12 division operations are used.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 66

Com putation in Projective Space

• 0 0 1, … , 3

Rewriting

• 1

1 1 1

• , , ,

, , ,

, , ,

1, 1, 1,1 A solution of the linear system of equations (LSE) is equivalent to generalized cross product:

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 67

Com putation in Projective Space For 0 each column can be multiplied by the relevant

• where again

, , , , , , , , , , , , Barycentric coordinates (Euclidean): 0 : , , 1 0 : , , 1 0 : , , 1 => Tests – point in triangle, point in tetrahedron - all in hom ogeneous coordinates => new entities: projective scalar, projective vector

(Skala,V.: Barycentric coordinates computation in homogeneous coordinates, Computers&Graphics, 2008)

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 68

Com putation in Projective Space Area of a triangle Volume of a tetrahedron 1 2

.

⁄ 1 6

.

⁄ As the principle of duality is valid, one could ask: What is a “dual” value to a computation of the area if the triangle is given by three lines in the “normalized” form, e.g.

. instead of three points? .

• 1
• . 2

⁄ / It can be seen that / , where: is the length of the line segment on and is a radius of the circumscribing circle. => value can be used as criterion for a quality triangular meshes.

• Skala,V.: Geometric Computation, Duality and Projective Space, IW-LGK workshop proceedings, ISBN 978-3-

86780-244-4, pp.105-111, Dresden University of Technology, 2011

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 69

Com putation in Projective Space Line in E3 as Tw o Plane I ntersection Standard formula in the Euclidean space , , :

: , , : :

Line given as an intersection of two planes , ,

• The form ula is quite “horrible” one and for students not acceptable

as it is too com plex and they do not see from the form ula com es from .

E

C S L I

Eurographics 2

Com puta Space Line in E I ntersec

• Stand

Euclid

• Plück

a line point is dua a pla an in plane a dua But c

• Proje

2013

ation in E3 as Tw ction dard form dean spa ker’s coor e is given

• s. DUALI

al to ne and v tersectio es can be al problem computat ective for

Eu

Vaclav Sk

Project w o Plane mula in th ace rdinates n by two ITY – a p vice versa

• n of two

e comput m tionally e mulation

urogr

kala http

ive e he point a => ted as expensiv n and sim

raphic

p://www.Vacla

ve compu mple com z

cs 2 0

avSkala.eu 2

utation putation

1 3

2013-05-03 70

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 71

Com putation in Projective Space Line in E3 as Tw o Plane I ntersection , , : , , : normal vectors are , , , , Directional vector of a line

• f two planes and is given as

“starting” point ??? A plane passing the origin with a normal vector , , , : 0 The point is defined as Sim ple form ula for m atrix-vector architectures like GPU and parallel

• processing. Compare the standard and projective formulas

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 72

Com putation in Projective Space – the nearest point Find the nearest point on an intersection of two planes to the given point ξ Simple solution:

• Translate planes and so the ξ

is in the origin

• Compute intersection of two planes

i.s. and

• Translate using

Again – an elegant solution, sim ple form ula supporting m atrix-vector architectures like GPU and parallel processing Solution DETAILS next

Known solution using Lagrange multipliers

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ ⋅ = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡

2 2 1 1 2 2 2 1 1 1 2 1 2 1 2 1

2 2 2 2 2 2 n p n p p p p p p p n n n n n n n n n n n n

z y x z y x z y x z y x z z y y x x

• μ

λ

Krumm,J.: Intersection of Two Planes, Microsoft Research

ξ

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 73

The closest point to an intersection of tw o planes In some applications we need to find a closest point on a line given as an intersection

• f two planes. We want to

find a point , the closest point to the given point , which lies on an intersection

• f two planes

: and :

This problem was recently solved by using Lagrange multipliers and an optimization approach leading to a solution of a system of linear equations with 5 equations.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 74

Solution in the projective space 1.Translate the given point , , : 1 to the origin – matrix 2.Compute parameters of the given planes and after the transformation as

and ,

3.Compute the intersection of those two planes

and

• 4.Transform the point to the original coordinate system using

transformation

: 0

• 1
• 1
• 1
• 1
• 1

1 1

• 1
• 1
• 1
• 1
• 1
• It is simple, easy to implement on GPU.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 75

Curves and surfaces Rational Bézier curve – Euclidean , , ∑

• ∑
• 0 1
• 1

1st derivative

• quite complicated

Projective , , :

• How simple !

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 76

Com putation in Projective Space Disadvantages

• Careful handling with formula in the projective space
• “Oriented” projective space is to be used, i.e. 0;

HW could support it – simple solution

• Exponents of homogeneous vectors can overflow – special handling
• exponents should be normalized; HW could support it

however not supported by the current hardware

• P_Lib – library for computation in the projective space - uses SW

solution for normalization on GPU (C# and C++), GPU support

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 77

Com putation in Projective Space Advantages

• “Infinity” is well represented
• No division operation is needed, a division operation can be hidden to the

homogeneous coordinate

• Many mathematical formula are simpler and elegant
• One code sequence solves primary and dual problems
• Supports matrix–vector operations in hardware – like GPU etc.
• Numerical computation can be faster
• Precision is nearly 2-times higher in mantissa and exponent -
• More robust and stable solutions can be achieved
• System of linear equations can be solved directly without division
• peration, if exponent normalization is provided

E

I

Eurographics 2

I m plem e

• GPU

vecto

• Nativ

expon

• Progr

2013

entation (Graphic

• r operat

ve arithm nent “no rammabl

Eu

Vaclav Sk

n aspect cal Proces ion – esp metic oper rmalizati e HW – p

urogr

kala http

ts and G ssing Uni pecially f rations w ion” parallel p

raphic

p://www.Vacla

PU it) -optim

• r [x,y,z

with homo processin

cs 2 0

avSkala.eu 2

mized for :w]T

• geneou

g

1 3

2013-05-03

matrix-v s coordin vector, v nates – w ector- without

78

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 79

I m plem entation aspects and GPU 4D cross product can be implemented in Cg/HLSL on GPU (not optimal implementation) as: float4 cross_4D(float4 x1, float4 x2, float4 x3) { float4 a; # simple formula # a.x=dot(x1.yzw, cross(x2.yzw, x3.yzw)); a.y=-dot(x1.xzw, cross(x2.xzw, x3.xzw)); a.z=dot(x1.xyw, cross(x2.xyw, x3.xyw)); a.w=-dot(x1.xyz, cross(x2.xyz, x3.xyz)); return a; } # more compact formula available #

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 80

Appendix Data processing - main field in computer science Data processing itself can be split to two main areas:

• processing of textual data

limited interval of values, unlimited dimensionality [a char as one dimension - Methionylthreonylthreonylglutaminylarginyl..isoleucine 189,819 chars] No interpolation is defined

• processing of num erical data

unlimited interval of values, limited dimensionality – usually 2 or 3 Interpolation can be used

Textual Graphical Dim ∞ 2, 3 Interval 0-255(ASCII) (-∞, ∞)

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 81

18 1 10 100 1 201 401 601 801 1001 Number of buckets Bucket length

Hash functions

• usually used for textual data processing
• prime numbers and modulo
• perations are used

Usual form 3 5 7 multiplication int * float needed

271 1 1 10 100 1000 1 2001 4001 6001 8001 10001 12001 Number of buckets Bucket length

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 82

173878

1 10 100 1000 10000 100000 1000000 1 10 100

Number of clusters Cluster length

data set: A4_unterbau1.stl

If the hash function is constructed as where , , are “irrational” numbers and 2 better distribution is obtained => much faster processing. is a “magic” constant which maps the expression for , to 0, 2 1

• Hradek,J., Skala,V.: Hash Function and Triangular Mesh Reconstruction, Vol.29, No.6., pp.741-751,

Computers&Geosciences, Pergamon Press, ISSN 0098-3004, 2003

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 83

1 10 100 1000 10000 0,5

q

English dictionary (start) Textual processing The has function is constructed as

• „irrational“ 0 1

2 1 Both geometrical and textual hash function design have the same approach coefficients are “irrational” and no division operation is needed. Some differences for Czech, Hebrew, English, German, Arabic, … languages and “chemical” words.

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 84

Sum m ary and conclusion We have got within this course an understanding of:

• projective representation use for geometric transformations with points,

lines and planes

• principle of duality and typical examples of dual problems, influence to

computational complexity

• intersection computation of two planes in E3, dual Plücker coordinates

and simple projective solution

• geometric problems solution with additional constrains
• intersection computations and interpolation algorithms directly in the

projective space

• barycentric coordinates computation on GPU
• avoiding or postponing division operations in computations

Projective space representation supports matrix-vector architectures like GPU – faster, robust and easy to implement algorithms achieved

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 85

References

• Skala,V.: Projective Geometry and Duality for Graphics, Games and Visualization - Course

SIGGRAPH Asia 2012, Singapore, ISBN 978-1-4503-1757-3, 2012

• Skala,V.: Barycentric Coordinates Computation in Homogeneous Coordinates, Computers &

Graphics, Elsevier, ISSN 0097-8493, Vol. 32, No.1, pp.120-127, 2008

• Skala,V.: Intersection Computation in Projective Space using Homogeneous Coordinates, Int.

Journal of Image and Graphics, ISSN 0219-4678, Vol.7, No.4, pp.615-628, 2008

• Skala,V.: Length, Area and Volume Computation in Homogeneous Coordinates, Int. Journal of

Image and Graphics, Vol.6., No.4, pp.625-639, ISSN 0219-4678, 2006

• Skala,V., Kaiser,J., Ondracka,V.: Library for Computation in the Projective Space, 6th Int.Conf.

Aplimat, Bratislava, ISBN 978-969562-4-1, pp. 125-130, 2007

• Skala,V.: GPU Computation in Projective Space Using Homogeneous Coordinates , Game

Programming GEMS 6 (Ed.Dickheiser,M.), pp.137-147, Charles River Media, 2006

• Skala,V.: A new approach to line and line segment clipping in homogeneous coordinates, The

Visual Computer, Vol.21, No.11, pp.905-914, Springer Verlag, 2005 Generally: “Publications w ith on-line DRAFTs” via http:/ / w w w .VaclavSkala.eu References related

• Agoston,M.K.: Computer Graphics and Geometric Modeling - Mathematics, Springer, 2005
• Ammeral,L: Programming Principles in Computer Graphics, John Wiley, 1986
• Miller,J.R.: The Mathematics of Graphical Transformations: Vector Geometric and Coordinate-

Based Approaches, DesignLab, 1997

• Yamaguchi,F.: Computer-Aided Geometric Design: A Totally Four-Dimensional Approach,

Springer, 2002

Eurographics 2 0 1 3

Eurographics 2013 Vaclav Skala http://www.VaclavSkala.eu 2013-05-03 86

?? Questions ??

Contact: Vaclav Skala c/o University of West Bohemia CZ 306 14 Plzen, Czech Republic http:/ / w w w .VaclavSkala.eu skala@kiv.zcu.cz subj. EG 2013 http:/ / Geom etryAlgebra.zcu.cz – Geometric Algebra

I NVI TATI ON

2 1 st W SCG conference on Com puter Graphics, Visualization and Com puter Vision 2 0 1 3 Plzen [ Pilsen] close to Prague, June 2 4 -2 7 , 2 0 1 3 http:/ / w w w .w scg.eu http:/ / w w w .w scg.cz

Supported by the Ministry of Education of the Czech Republic, projects No.LH12181, LG13047