- By
**azure** - Follow User

- 137 Views
- Uploaded on

Download Presentation
## PowerPoint Slideshow about 'Computational Geometry' - azure

**An Image/Link below is provided (as is) to download presentation**

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript

Our First Problem

- Nearest neighbor searching
- Applications?
- Pattern Classification
- Graphics
- Data Compression
- Document retrieval
- Statistics
- Machine Learning
- …

Similarity Measure

Similar?

Similarity measure

- Other similarity Measures

The dimension

- Lets assume that our points are in one dimensional space. ( d = 1 ). We will generalize to higher dimension ( When d = some constant ).

Fixed radius near neighbor problem

Question- Given a set of points S on the real line, preprocess them to answer the following question :
- Find all the pair of points (p,q) such that distance of (p,q) < r .

q

Points in S

Question

- Given a set of points S on the real line, preprocess them to answer the following query :
- Given a query point q, find the neighbor of q which is closest in S.

q

nn(q)

Points in S

Answers?

- Fixed NN Search
- O( n2 ) ?
- O(nlogn + k) ?
- O(n + k) ?
- NN Search
- O( n ) ?
- O( log n) ?

Answers?

- NN Search
- O( n ) ? Brute Force [ Trivial ]
- O( log n) ? Binary Search Tree
- Fixed NN Search
- O( n2 ) ? Brute Force
- O(nlogn + k) ? Sorting
- O(n + k) ? Hashing?

K-nearest neighbor search

- Problem: Given a set of points P on the real line and a query point q, find the k-nearest neighbors of q in P.
- O(nlogn) Trivial bruteforce
- Do you see how?
- Thought Problem: How do we do this in O(n) time?
- (Hint: Median finding works in O(n) time).

Brute Force implementation

- min = +
- for( i = 0; i < n; ++i ) for( j = 0; j < n; ++j )
- (D = dist (pj, pi )) < r ) ? ( cout << (pj, pi ) ) : ();
- float dist ( point p, point q )
- sum = 0;
- for ( i = 0; i < d; ++i ) sum += (pi – qi)2
- return sqrt(sum)

What can we speed up here?

Brute Force implementation

- min = +
- for( i = 0; i < n; ++i ) for( j = 0; j < n; ++j )
- (D = dist (pj, pi )) < r ) ? ( cout << (pj, pi ) ) : ();
- float dist ( point p, point q )
- sum = 0;
- for ( i = 0; i < d; ++i ) sum += (pi – qi)2
- return sqrt(sum)

How do we speed this up?

Fixed NN Search: By Sorting

- Once we sort the points on the real line, we can just go left and right to identify the pairs we need.
- Each pair is visited at most twice, so the asymptotics do not change.

Total work done after sorting

- ki denotes the pairs generated when visiting pi
- With this approach, we need at least Ω(nlogn) (for sorting).

Fixed radius near neighbor searching

- How do we avoid sorting?
- How do we get a running time of O(n+k) ?
- How do we use hashing to simplify our problem.

b=0

{

r

0

Solution using bucketing- Put points in infinite number of buckets (Assume an infinite array B)

Interval b is [ br, (b+1)r ]

x lies in b = floor (x/r)

Solution using bucketing

- Only n buckets of B might get occupied at most.
- How do we convert this infinite array into a finite one: Use hashing
- In O(1) time, we can determine which bucket a point falls in.
- In O(1) expected time, we can look the bucket up

in the hash table

- Total time for bucketing is expected O(n)
- The total running time can be made O(n) with high probability using multiple hash functions ( essentially using more than one hash function and choosing one at run time to fool the adversary ).

The Algorithm

- Store all the points in buckets of size r
- In a hash table [ Total complexity = O(n) ]
- For each point x
- b = floor(x/r); Retrieve buckets b, b+1
- Output all pairs (x,y) such that y is either in bucket b or b+1 and x < y and ||xy|| < r

x

0

Running Time

- Let nb denote the number of points in bucket b of the input pointset P.
- Define

Note that there are nb2 pairs in bucket b alone that are within distance r of each other.

Observation

- Since each pair gets counted twice :

Higher Dimensions

- Send (x,y) (floor(x/r),floor(y/r))
- Apply hash with two arguments
- Running time still O(n+k)
- Running time increases exponentially with dimension

Introduction: Geometry Basics

- Geometric Systems
- Vector Space
- Affine Geometry
- Euclidean Geometry
- AG + Inner Products = Euclidean Geometry

Vector Space

- Scalar ( + , * ) = Number Types
- Usual example is Real Numbers R.
- Let V be a set with two operations
- + : V x V V
- * : F x V V
- Here F is the set of Scalars

Vector Space

- If (V , +, * ) follows the following properties, its called a vector space :
- (A1) u + (v + w) = (u + v) + w for all u,v,w in V.
- (A2) u + v = v + u for all u,v in V.
- (A3) there is unique 0 in V such that 0 + u = u for all u in V.
- (A4) for every u in V, there is unique -u in V such that u + -u = 0.
- (S1) r(su) = (rs)u for every r,s in R and every u in V.
- (S2) (r +s)u = ru + su for every r,s in R and every u in V.
- (S3) r(u + v) = ru + rv for every r in R and every u,v in V.
- (S4) 1u = u for every u in V.

Note: Vectors are closed under linear combinations

A basis is a set of n linearly independent vectors that span V.

Affine Geometry

- Geometry of vectors
- Not involving any notion of length or angle.
- Consists of
- A set of scalars
- Say Real numbers
- A set of points
- Position specification
- A set of free vectors.
- Direction specification

Affine Geometry

- Legal operations
- Point - Point = Vector
- Point +/- Vector =Point
- Vector +/- Vector = Vector
- Scalar * Vector = Vector
- Affine Combination
- ∑Scalar * Points = Point
- such that ∑Scalar = 1
- Note that scalars can range from –Infinity to +Infinity

Affine Combinations

- [ Affine Span or Affine Closure ] The set of all affine combinations of three points generates a plane.
- [ Convex Closure ] The set of all convex combinations of three points generates all points inside a triangle.

Euclidean Geometry

- One more element added
- Inner Products
- Maps two vectors into a scalar
- A way to `multiply’ two vectors

Example of Inner Products

- Example : Dot products
- (u.v) = u0v0+ u1v1+…+ ud-1vd-1
- Where u,v are d-dimensional vectors.
- u = (u0,u1,…, ud-1); v = (v0,v1,…, vd-1)
- Length of |u| = sqrt(u.u) (Distance from origin)
- Normalization to unit length : u/|u|
- Distance between points |p - q|
- Angle (u’,v’) = cos-1(u’.v’)
- where u’=u/|u| and v’=v/|v|

Dot products

- (u.v) = (+/-)|u|(projection of v on u).
- u is perpendicular to v (u,v) = 0
- u.(v+w) = u.v + u.w
- If u.u not equal to zero then u.u > 0
- positive definite

Some proofs using Dot Products

- Cauchy Schwarz Inequality
- (u.v) <= |u||v|
- Homework.
- Hint: For any real number x
- (u+xv).(u+xv) >= 0
- Triangle Inequality
- |u+v|<=|u|+|v|
- Hint: expand |u+v|2 and use Cauchy Schwarz.

v

u

Next Lecture

- Orientation and Convex hulls

Due Wed

Homework: (Programming)

Play with dpoint.hpp

and example.hpp (Implement orientation in 2D).

Implement your own dvector.hpp, dsegment.hpp using metaprogramming

Make sure you understand how things work.

Due Monday

Homework: (Theory)

Cauchy Schwarz

Triangle Inequality

dot n cross prod.

Reading Assignment:

Page 1-10, Notes of Dr. Mount

Page 1-2 and Section 1.3 of the text book.

Sources for this lecture: Dr. David Mount’s Notes.

WWW.

Crash course on C++“dpoint.hpp”

- Namespaces
- Solve the problem of classes/variables with same name
- namespace _cg {code}
- Outside the namespace use _cg::dpoint
- Code within _cg should refer to code outside _cg explicitly. E.g. std::cout instead of cout.

Object Oriented Programming

- Identify functional units in your design
- Write classes to implement these functional units
- Separate functionality for code-reuse.

Class membership

- Public
- Private
- Always : Keep member variables private
- This ensures that the class knows when the variable changes
- Protected

Inheritance

- ‘is a’ relationship is public inheritance
- Class SuperDuperBoss : public Boss
- Polymorphism : Refer an object thru a reference or pointer of the type of a parent class of the object
- SuperDuperBoss JB;
- Boss *b = &JB;
- Virtual functions

Templates

- Are C macros on Steroids
- Give you the power to parametrize
- Compile time computation
- Performance

“The art of programming programsthat read, transform, or write other programs.”

- François-René Rideau

Generic Programming

- How do we implement a linked list with a general type inside?
- void pointers?
- Using macros?
- Using Inheritance?

Templates

- Function Templates
- Class Templates
- Template templates *
- Full Template specialization
- Partial template specialization

Metaprogramming

- Programs that manipulate other programs or themselves
- Can be executed at compile time or runtime.
- Template metaprograms are executed at compile time.

Good old C

- C code
- Double square(double x) { return x*x; }
- sqare(3.14) Computed at compile time
- #define square(x) ((x)*(x))
- Static double sqrarg;

#define SQR(a) (sqrarg=(a), sqrarg*sqrarg)

Templates

- Help us to write code without being tied to particular type.
- Question: How do you swap two elements of any type? How do you return the square of any type?

Function Templates

- C++
- template< typename T >

inline T square ( T x ) { return x*x; }

- A specialization is instantiated if needed :
- square<double>(3.14)
- Template arguments maybe deduced from the function arguments square(3.14)
- MyType m; … ; square(m); expands to square<MyType>(m)

Operator * must be overloaded for MyType

Function Templates

- template<typename T>

void swap( T& a, T& b ){

T tmp(a); // cc required

a = b; // ao required

b = tmp;

}

Mytype x = 1111;

Mytype y = 100101;

swap(x,y); swap<Mytype>(…) is instantiated

- Note reliance on T’s concepts (properties):
- In above, T must be copyable and assignable
- Compile-time enforcement (concept-checking) techniques available

Function Template Specialization

- template<>

void swap<myclass>( myclass& a, myclass& b){

a = b = 0;

}

Custom version of a template for a specific class

Class Templates

- template<typename NumType, unsigned D> class dpoint{

public:

NumType x[D];

};

A simple 3-dimensional point.

dpoint<float,3> point_in_3d;

point_in_3d.x[0] = 5.0;

point_in_3d.x[1] = 1.0;

point_in_3d.x[2] = 2.0;

Note the explicit instantiation

Class Templates: Unlike function templates

- Class template arguments can’t be deduced in the same way, so usually written explicitly:

dpoint <double, 2> c; // 2d point

- Class template parameters may have default values:

template< class T, class U = int >

class MyCls { … };

- Class templates may be partially specialized:

template< class U >

class MyCls< bool, U > { … };

Using Specializations

- First declare (and/or define) the general case:

template< class T >

class C { /* handle most types this way */ };

- Then provide either or both kinds of special cases as desired:

template< class T >

class C< T * > { /* handle pointers specially */ };

template<>// note: fully specialized

class C< int * > { /* treat int pointers thusly */ };

- Compiler will select the most specialized applicable class template

Template Template Parameters

- A class template can be a template argument
- Example:

template< template<class> class Bag >

class C { // …

Bag< float > b;

};

- Or even:

template<class E, template<class> class Bag >

class C { // …

Bag< E > b;

};

More Templates

template<unsigned u>

class MyClass {

enum { X = u };

};

Cout << MyClass<2>::X << endl;

Template Metaprograms

- Factorials at compile time

template<int N> class Factorial {

public:

enum { value = N * Factorial<N-1>::value };

};

class Factorial<1> {

public:

enum { value = 1 };

};

Int w = Factorial<10>::value;

Template Metaprograms

- Metaprogramming using C++ can be used to implement a turning machine. (since it can be used to do conditional and loop constructs).

Template Metaprograms

- template< typename NumType, unsigned D, unsigned I > struct origin
- {
- static inline void eval( dpoint<NumType,D>& p )
- {
- p[I] = 0.0;
- origin< NumType, D, I-1 >::eval( p );
- }
- };
- // Partial Template Specialization
- template <typename NumType, unsigned D> struct origin<NumType, D, 0>
- {
- static inline void eval( dpoint<NumType,D>& p )
- {
- p[0] = 0.0;
- }
- };

const int D = 3;

inline void move2origin() { origin<NumType, D, D-1>::eval(*this); };

Food for thought

- You can implement
- IF
- WHILE
- FOR

Using metaprogramming. And then use them in your code that needs to run at compile time

For the more challenged: Implement computation of determinant

/ orientation / volume using metaprogramming. (Extra credit)

Traits Technique

- Operate on “types” instead of data
- How do you implement a “mean” class without specifying the types.
- For double arrays it should output double
- For integers it should return a float
- For complex numbers it should return a complex number

Traits

Template< typename T >

struct average_traits{

typedef T T_average;

};

Template<>

struct average_traits<int>{

typedef float T;

}

Traits

template<typename T>

typename average_traits<T>::T_average

average(T* array, int N){

typename avearage_traits<T>::T_average

result = sum(array,N);

return result/N;

}

Sources

- Template Metaprogramming by Todd Veldhuizen
- C++ Meta<Programming> Concepts and results by Walter E. Brown
- C++ For Game programmers by Noel Llopis
- C++ Primer by Lippman and Lajoie

GNU MP

- GMP is a free library for arbitrary precision arithmetic, operating on signed integers, rational numbers, and floating point numbers.
- Get yourself acquinted with :
- Mpz : Integer arithmetic
- mpq : Rationals
- Mpf : Floats

Orientation primitive for points

+ve in this case

Orientation

y

+

The sign (+ or –) of cross

product depends in which

half plane (relative to p1)

lies p2

Intuition for 2D and 3D using

the right hand rule?

p1

–

(0,0)

x

Homeworks

- Prove that for two vectors v1 = (v1x,v1y) and v2 = (v2x,v2y)
- v1x*v2x + v1y *v2y = |v1||v2| cos Ө
- v1y*v2x - v2y *v1x = |v1||v2| sin Ө
- What do these quantities mean geometrically? Note that one changes sign when Ө is negative.

Problems with fixed point arithmetic

: 3.141592653589793238462643383....

float pi = 2 * asin(1); printf("%.35f\n", pi);

Outputs

3.14159274101257324218750000000000000

float typically are 32 bits and can deal with 7 digits after decimal

Exaggerating the error

- Suppose our computer computes in base 10, but only keeps the two most significant digits. We call this fixed precision computation
- Let us take three points:
- p = (-94,0)
- q = (92,68)
- r = (400,180)
- Triangle pqr is oriented clockwise, hence area < 0.

Determinant

- True value: (92+94)*180-(400+94)*68=186*180-494*68=33480-33592 =-112
- Fixed precision computation : (92+94)*180-(400+94)*68=190*180-490*68=34000-33000 =+1000

Floating point computation

- Sometimes produces wrong results
- Exact computation is possible (And hence we will use GMP : An easy way out)
- Exact computations are slow and can be avoided.
- One simple method : Use floating point to estimate whether you really need Exactness, only then invoke exact computation.

p9

p5

p12

p7

p4

p11

p1

p8

p2

p0

Convex hulls- The smallest convex body that contains a set of points
- The intersection of all convex sets that contain S.

Download Presentation

Connecting to Server..