computational geometry n.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
Computational Geometry PowerPoint Presentation
Download Presentation
Computational Geometry

Loading in 2 Seconds...

play fullscreen
1 / 76

Computational Geometry - PowerPoint PPT Presentation


  • 137 Views
  • Uploaded on

Computational Geometry. Piyush Kumar (Lecture 2: NN Search). Welcome to CIS5930. Our First Problem . Nearest neighbor searching Applications? Pattern Classification Graphics Data Compression Document retrieval Statistics Machine Learning …. Similarity Measure.

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
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
computational geometry

Computational Geometry

Piyush Kumar

(Lecture 2: NN Search)

Welcome to CIS5930

our first problem
Our First Problem
  • Nearest neighbor searching
    • Applications?
      • Pattern Classification
      • Graphics
      • Data Compression
      • Document retrieval
      • Statistics
      • Machine Learning
similarity measure
Similarity Measure
  • In terms of Euclidean distance

(4,5)

(2,3)

similarity measure2
Similarity measure
  • Other similarity Measures
the dimension
The dimension
  • Lets assume that our points are in one dimensional space. ( d = 1 ). We will generalize to higher dimension ( When d = some constant ).
question

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

question1

Nearest neighbor search

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
Answers?
  • Fixed NN Search
    • O( n2 ) ?
    • O(nlogn + k) ?
    • O(n + k) ?
  • NN Search
    • O( n ) ?
    • O( log n) ?
answers1
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?
nn searching balanced binary tree
NN Searching : Balanced binary tree

O( log n )

q

nn(q)

Points in S

k nearest neighbor search
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
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 implementation1
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
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 sorting1
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
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.
solution using bucketing

b= -2

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 bucketing1
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
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
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
Observation
  • Since each pair gets counted twice :
running time1
Running Time
  • Depends on the number of distance computations D.

Total Running Time = O(n+k)

higher dimensions
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
Introduction: Geometry Basics
  • Geometric Systems
    • Vector Space
    • Affine Geometry
    • Euclidean Geometry
      • AG + Inner Products = Euclidean Geometry
vector space
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 space1
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
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 geometry1
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 geometry2
Affine Geometry

Affine Combination

Convex Combination

affine combinations
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
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 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
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
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
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
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
Object Oriented Programming
  • Identify functional units in your design
  • Write classes to implement these functional units
  • Separate functionality for code-reuse.
class membership
Class membership
  • Public
  • Private
    • Always : Keep member variables private
    • This ensures that the class knows when the variable changes
  • Protected
inheritance
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
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
Generic Programming
  • How do we implement a linked list with a general type inside?
    • void pointers?
    • Using macros?
    • Using Inheritance?
templates1
Templates
  • Function Templates
  • Class Templates
  • Template templates *
  • Full Template specialization
  • Partial template specialization
metaprogramming
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
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)

templates2
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
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 templates1
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
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
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 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
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
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
More Templates

template<unsigned u>

class MyClass {

enum { X = u };

};

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

template metaprograms
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 metaprograms1
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 metaprograms2
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
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
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
Traits

Template< typename T >

struct average_traits{

typedef T T_average;

};

Template<>

struct average_traits<int>{

typedef float T;

}

traits1
Traits

average_type(T) = T

average_type(int) = float

traits2
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
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
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
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
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
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
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
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
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.
convexity

p

q

p

convex

non-convex

q

Convexity

A set S is convex if for any pair of points p,q S we have pq  S.

convex hulls

p6

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.