geomalgo aims at prodiving basic geometric 2D and 3D algorithms, for example computing the area of a 2D triangle, or finding intersection points of 3D triangle and a 3D segment.

geomalgo algorithms are primary based on geomalgorithms.

geomalgo provides a convenient Python API, and a Cython API to obtain C performances, typically for computation inside loops.

Installation

Install the Conda package manager, and install geomalgo with:

conda install -c dfroger geomalgo  # -c conda-forge soon

Documentation

Warning

This documentation is being drafted. It should be complete around in a few weeks (around May 2017).

Tutorial

This tutorial introduces geomalgo through simple examples, exploring three differents aspects of its API:

• Python API: Manipulate Python objects in the Python interactive interpreter.
• Cython API: Use Cython compiler to add static typing and manipulate C structures.
• Work on large data sets: Use SoA memory layout and compile to C code to efficiently work large date sets.

These concepts will be demonstrated on a simple example: computing intersection between segments.

Using the Python interface

This tutorial is written in a Jupyter notebook. pylab is imported, and matplotlib figures are inlined in the notebook. geomalgo is imported as ga.

%pylab inline

import geomalgo as ga

Populating the interactive namespace from numpy and matplotlib

%load_ext Cython

%%cython
cimport geomalgo as ga
cdef:
    ga.CSegment2D seg

Creating points and segments

The first geomalgo class (actually, extension type) we manipulate is geomalgo.Point2D. It takes point coordinates as argument. The Python API also accept an optional name, and provides a plot method, usefull to quickly create matplotlib figures.

A = ga.Point2D(2, 1, name='A')
B = ga.Point2D(6, 4, name='B')

C = ga.Point2D(5, 2, name='C')
D = ga.Point2D(2, 3.5, name='D')

for obj in [A, B, C, D]:
    obj.plot()

These four points can be used to create two segments, using geomalgo.Segment2D. Like points, segments have a plot method.

AB = ga.Segment2D(A, B)
CD = ga.Segment2D(C, D)

for obj in [A, B, C, D, AB, CD]:
    obj.plot()

Computing intersection

The next task we want to achieve is to compute the intersection between the two segments.

There is a Segment2D.intersect_segment() method to compute intersection with another segment. It returns two objects, which can be:

• both None if segments do not intersect,
• a point and None if there is one intersection point,
• two points when segments are colinear.

I0, I1 = AB.intersect_segment(CD)
print(I0)
print(I1)

for obj in [A, B, C, D, AB, CD]:
    obj.plot()

I0.plot(name='I', color='red')

Point2D(4.0, 2.5)
None

Segment relative coordinates

Another usefull information is to locate the intersection point not only with cartesian coordiante, but also with coordiante relative to each segment. This is achieve with the return_coords flag.

Note that as there is not second intersection point, we can ignore None value, assigning them to a _ variable.

I0, _, (AB_coord, CD_coord, _, _)  = AB.intersect_segment(CD, return_coords=True)
print(AB_coord, CD_coord)

0.5 0.3333333333333333

AB_coord is \(\frac{1}{2}\), meaning the point is in the middle of the segment, while CD_coord is \(\frac{1}{3}\), meaning the point is at a distance \(\frac{|CD|}{3}\) from the point C.

To change from cartesian coordinates to segment coordiante, use the Point2D.at() method:

print(AB.at(AB_coord))
print(CD.at(CD_coord))

Point2D(4.0, 2.5)
Point2D(4.0, 2.5)

To change from segment coordinate to cartesian coordiantes, use the Point2D.where() method:

print(AB.where(I0))
print(CD.where(I0))

0.5
0.3333333333333333

Conlusion

This concludes the first tutorial part.

The Python interpreter has allowed us to interactively construct points, segments, and compute intersection between them. This a very handy way to perform temporary computation, using Python interpreter as an advanced calculator, or test a idea, write a prototype, etc.

The code can then be written to a Python script to be reused, and to build a library or an application. All the code will be interpreted and dynamically typed. If performance become an issue, we will want to add static typing. This is what will be covered in next section.

Using the Cython interface

geomalgo provides a cython API, making it easy to:

• compile code instead of interpret it,
• add static typing for extension types, C built in, functions, etc,
• compile to full performance C code, using C structures instead of extension types.

We first activate cython compilation in the notebook.

%load_ext Cython
%load_ext wurlitzer

The Cython extension is already loaded. To reload it, use:
  %reload_ext Cython

Compiling the code

Instead of interpreting source code, cython compiles cython code to C code, and the build a native Python module. Learn more about Cython. For example, let’s rewrite the code in Cython. The %%cython magic do (learn more about building cython code):

• compile Cython code to C code,
• compile C code to Python native module,
• execute the module (hence, the Cython code).

%%cython

import geomalgo as ga

A = ga.Point2D(2, 1, name='A')
B = ga.Point2D(6, 4, name='B')
C = ga.Point2D(5, 2, name='C')
D = ga.Point2D(2, 3.5, name='D')
AB = ga.Segment2D(A, B)
CD = ga.Segment2D(C, D)

I0, I1 = AB.intersect_segment(CD)
print(I0)
print(I1)

Point2D(4.0, 2.5)
None

Adding static typing

Once the code is written in Cython, which compile to C, static typing can be added. The above code rewrittes to (notice the cimport instead of import):

%%cython

cimport geomalgo as ga

cdef:
    ga.Point2D A, B, C, D
    ga.Segment2D AB, CD

A = ga.Point2D(2, 1, name='A')
B = ga.Point2D(6, 4, name='B')
C = ga.Point2D(5, 2, name='C')
D = ga.Point2D(2, 3.5, name='D')
AB = ga.Segment2D(A, B)
CD = ga.Segment2D(C, D)

I0, I1 = AB.intersect_segment(CD)
print(I0)
print(I1)

Point2D(4.0, 2.5)
None

Using C structure instead of extension type

Calling extention type methods is as fast as calling C functions on C structure, provided that extension type attributes and methods are all statically typed (cdef). However, extension types can not be stack allocated as efficiently as C structure. GeomAlgo relies on stack allocation for basic types (points, segments, etc). For performance, we need be able to use C structures, instead of extension types.

At the end of this section, we will see that extension are wrapper around these C structures, and methods wrapper around C functions.

The above code is rewritten using CPoint2D, CSegment2D and intersect_segment2d_segment2d().

%%cython

from libc.stdio cimport printf

cimport geomalgo as ga

cdef:
    ga.CPoint2D A, B, C, D, I0, I1
    ga.CSegment2D AB, CD
    double coords[4]
    int n

A.x, A.y = 2, 1
B.x, B.y = 6, 4

C.x, C.y = 5, 2
D.x, D.y = 2, 3.5

ga.segment2d_set(&AB, &A, &B)
ga.segment2d_set(&CD, &C, &D)

n = ga.intersect_segment2d_segment2d(&AB, &CD, &I0, &I1, coords)

printf("Number of intersection: %d\n", n)
printf("Intersection point: (%.2f, %.2f)\n", I0.x, I0.y)

Number of intersection: 1
Intersection point: (4.00, 2.50)

This generate “pure” C code, without reference to Python. It is as efficient as writting in C directly.

Relations between Python and Cython API

A Point2D is a wrapper around a CPoint2D, and coordinates are hold by CPoint2D.

For example:

%%cython

cimport geomalgo as ga

cdef:
    ga.Point2D X
    ga.CPoint2D* C_X

X = ga.Point2D(1, 2)

C_X = X.cpoint2d

The Python API of a Segment2D extension type has two attributes:

• Point2D A (extension type)
• Point2D B (extension type)

The Cython API has access to an additionnal attribute:

• CSegment2D (C structure)

The CSegment2D struture has two members:

• CPoint2D* A (C structure)
• CPoint2D* B (C structure)

CSegment2D and Point2D A and B points to the same CPoint2D* A and B.

For example:

%%cython

cimport geomalgo as ga

cdef:
    ga.Point2D X, Y
    ga.Segment2D XY

X = ga.Point2D(1, 2)
Y = ga.Point2D(3, 4)

XY = ga.Segment2D(X, Y)

cdef:
    ga.CSegment2D* C_XY
    ga.CPoint2D* C_X
    ga.CPoint2D* C_Y

C_XY = &XY.csegment2d

# Points can be accessed via X and Y
C_X = X.cpoint2d

# Or via C_XY
C_Y = C_XY.B

Working with large dataset

To conclude this tutorial, in this last part, the segment intersection computation will be applied to a (potentialy large) collection of segment.

A = ga.Point2D(0, 0, name='A')
B = ga.Point2D(8, 8, name='B')
AB = ga.Segment2D(A, B)

segment_list = [
    ga.Segment2D((2,1), (1,2)),
    ga.Segment2D((3,2), (3,4)),
    ga.Segment2D((2,6), (3,7)),
    ga.Segment2D((5,5), (6,6)),
    ga.Segment2D((4,4), (5,4)),
    ga.Segment2D((6,8), (8,6)),
    ga.Segment2D((5,1), (6,1)),
]

figure(figsize=(8,8))

def plot_segments():
    for obj in [A, B, AB]:
        obj.plot()

    for i, seg in enumerate(segment_list):
        seg.plot(style='r-')
        text(seg.B.x, seg.B.y+0.1, i, color='red')

plot_segments()
axis('scaled')
grid()

For example, we want to compute all intersection of red segments with the blue one, indexed by red segment indices. This should be optimized for a large number of red segments.

Memory layout

Segment coordinates will be stored as SoA. geomalgo provides a SegmentCollection extension type to groups the 4 arrays together.

xa = np.array([seg.A.x for seg in segment_list], dtype='d')
xb = np.array([seg.B.x for seg in segment_list], dtype='d')

ya = np.array([seg.A.y for seg in segment_list], dtype='d')
yb = np.array([seg.B.y for seg in segment_list], dtype='d')

segments = ga.Segment2DCollection(xa, xb, ya, yb)

A CSegment2D and its 2 CPoints2D are declared to iterate on the collection.

The Segment2DCollection.get is used to set PQ for the current segment index.

%%cython

from libc.stdio cimport printf

cimport geomalgo as ga

def iter_segments(ga.Segment2DCollection segments):
    cdef:
        ga.CSegment2D PQ
        ga.CPoint2D P, Q
        int S

    ga.segment2d_set(&PQ, &P, &Q)

    for S in range(segments.size):
        segments.get(S, &PQ)
        printf("(%.0f, %.0f), (%.0f, %.0f)\n", P.x, P.y, Q.x, Q.y)

iter_segments(segments)

(2, 1), (1, 2)
(3, 2), (3, 4)
(2, 6), (3, 7)
(5, 5), (6, 6)
(4, 4), (5, 4)
(6, 8), (8, 6)
(5, 1), (6, 1)

Computing intersection

Once we can iterate the segment collection, intersection can be computed as in the previous part with the intersect_segment2d_segment2d function.

We need to choose where to store the results. Each red segment may intersection the blue one on 0, 1 or 2 points. The total number of intersections points in not knonw in advance. We choose to allocate arrays to store results for the maximal number of intersection points. We also choose to store only the coordinates relative to the segment AB.

%%cython

import numpy as np
cimport geomalgo as ga

def compute_intersections(ga.Segment2D AB, ga.Segment2DCollection segments):
    cdef:
        ga.CSegment2D PQ
        ga.CPoint2D P, Q
        ga.CPoint2D I0, I1
        double coords[4]
        int S, C

        # Count of intersections for each segment. May be 0, 1 or 2.
        int[:] count = np.empty(segments.size, dtype='int32')
        int countmax = 2

        # Intersection coordinate relative to AB.
        double[:,:] coord = np.empty((segments.size, countmax), dtype='d')

    ga.segment2d_set(&PQ, &P, &Q)

    for S in range(segments.size):
        segments.get(S, &PQ)
        count[S] = ga.intersect_segment2d_segment2d(&AB.csegment2d,
                                                    &PQ, &I0, &I1, coords)
        for C in range(count[S]):
            coord[S,C] = coords[C*2]

    return np.asarray(count), np.asarray(coord)

Let’s call this function on our segments, and plot the result.

counts, coords = compute_intersections(AB, segments)

print(counts)
print(coords)

[1 1 0 2 1 1 0]
[[  1.87500000e-001   2.64239983e-316]
 [  3.75000000e-001   0.00000000e+000]
 [  0.00000000e+000   0.00000000e+000]
 [  6.25000000e-001   7.50000000e-001]
 [  5.00000000e-001   0.00000000e+000]
 [  8.75000000e-001   0.00000000e+000]
 [  0.00000000e+000   0.00000000e+000]]

figure(figsize=(8,8))
plot_segments()
grid()
axis('scaled')

for S, (count, coord) in enumerate(zip(counts, coords)):
    for C in range(count):
        coord = coords[S, C]
        P = AB.at(coord)
        plot(P.x, P.y, 'ro')

A high level Python extension type may now be create to wrap the compute_intersections code and its results in a Pythonic API.

Conlusion

In this tutorial, we’ve explored the Python and Cython API of geomalgo, and how to easily move from one to another depending on the needs.

Python API allows interactive processing using the Python scientific stack (numpy, matplotlib, jupyter notebook, etc), while Cython API give use C like performances.

Collection of object are stored as SoA, which can be iterate using a stack allocated element structure: for example, a Segment2DCollection is allocated using 4 arrays, and can be iterated using a CSegment2D structure to be processed.

Points

There are two data structures to represent a point in two-dimensional space, of coordiantes \((x, y)\):

• the C structure CPoint2D, - the Python extension type Point2D, a wrapper around CPoint2D.

CPoint2D C structure

A CPoint2D has just two members, to store its coordinates:

CPoint2D

double x

double y

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A

A.x = 1
A.y = 2

printf("A: (%.1f, %.1f)\n", A.x, A.y)

A: (1, 2)

CPoint2D can be allocated and destroyed with these functions:

CPoint2D* new_point2d()

Allocate a new CPoint2D

del_point2d(CPoint2D* cpoint2d)

Delete a CPoint2D

Note

This functions are rarely used, as C structure in GeomAlgo are usually stack-allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D* A

A = ga.new_point2d()

A.x = 1
A.y = 2

printf("A: (%.1f, %.1f)\n", A.x, A.y)

ga.del_point2d(A)

A: (1, 2)

Point2D Python extension type

A Point2D takes coordinates as arguments, an optional index (can be for example vertice indices in a triangulation), and optional name (used for example by its Point2D.plot method).

class Point2D(x, y, index=0, name=None)

Attributes:

cpoint2d: CPoint2D*

index: int

name: string

A = ga.Point2D(1, 2, name='A')
B = ga.Point2D(4, 6, name='B')
print('A:', A)
A.plot()
B.plot()

A: Point2D(1.0, 2.0)

The wrapped C structure Point2D.cpoint2d is accessible only using Cython.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.Point2D A
    ga.CPoint2D* ptr

A = ga.Point2D(1, 2)

ptr = A.cpoint2d
printf("(%.1f, %.1f)\n", ptr.x, ptr.y)

(1, 2)

Get and set coordiantes

CPoint2D members are accessed directly to get and set coordinates:

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A

A.x, A.y = 1, 2
printf("A: (%.1f, %.1f)\n", A.x, A.y)

A.x, A.y = 4, -5
printf("A: (%.1f, %.1f)\n", A.x, A.y)

A: (1.0, 2.0)
A: (4.0, -5.0)

Point2D as two properties x and y to get/set coordiantes from/to its underlying cpoint2d attribute:

Point2D.x

Point2D.y

A = ga.Point2D(1, 2)
print(A)

A.x, A.y = 4, -5
print(A)

Point2D(1.0, 2.0)
Point2D(4.0, -5.0)

Point operators

Points subtraction

Two points \(B\) and \(A\) can be subtracted to compute the vector \(\mathbf{AB} = B - A\).

void subtract_points2d(CVector2D* AB, const CPoint2D* B, const CPoint2D* A)

Variable AB must be already allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A, B
    ga.CVector2D AB

A.x, A.y = 1, 2
B.x, B.y = 4, 6

ga.subtract_points2d(&AB, &B, &A)

printf("AB: (%.1f, %.1f)\n", AB.x, AB.y)

AB: (3.0, 4.0)

Point2D.__sub__(B, A)

A = ga.Point2D(1, 2)
B = ga.Point2D(4, 6)
AB = B - A
print('AB: ({}, {})'.format(AB.x, AB.y))

AB: (3.0, 4.0)

Point plus vector

A vector \(\mathbf{AB}\) can be added to a point \(A\) to compute the point \(B = A + \alpha \mathbf{AB}\).

void point2d_plus_vector2d(CPoint2D* B, CPoint2D* A, double alpha, CVector2D* AB)

B must be already allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A, B
    ga.CVector2D AB

A.x, A.y = 1, 2
AB.x, AB.y = 1.5, 2

ga.point2d_plus_vector2d(&B, &A, 2., &AB)

printf("B: (%.1f, %.1f)", B.x, B.y)

B: (4.0, 6.0)

Point2D.__add__(A, AB):

A = ga.Point2D(1, 2)
AB = ga.Vector2D(1.5, 2)

B = A + AB*2

print(B)

Point2D(4.0, 6.0)

Points equality

Test whether two points \(A\) and \(B\) are strictly equal.

bint point2d_equal(CPoint2D* A, CPoint2D* B)

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A, B, C

A.x, A.y = 1, 2
B.x, B.y = 1, 2
C.x, C.y = 2, 2

printf("A==B: %d\n", ga.point2d_equal(&A, &B))
printf("A==C: %d\n", ga.point2d_equal(&A, &C))

A==B: 1
A==C: 0

Distance computation

Compute the distance between points \(A\) and \(B\).

double point2d_distance(CPoint2D* A, CPoint2D* B)

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A, B
    double dist

A.x, A.y = 1, 2
B.x, B.y = 4, 6
dist = ga.point2d_distance(&A, &B)
printf("distance: %.1f", dist)

distance: 5.0

distance(A, B)

A = ga.Point2D(1, 2)
B = ga.Point2D(4, 6)
A.distance(B)

5.0

Sometime, the knowledge of the square distance is enough, and for performance, computing the square root can be avoided.

double point2d_square_distance(CPoint2D* A, CPoint2D* B)

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A, B
    double sdist

A.x, A.y = 1, 2
B.x, B.y = 4, 6
sdist = ga.point2d_square_distance(&A, &B)
printf("square distance: %.1f", sdist)

square distance: 25.0

Various

Test if the point \(P\) is left, right, or of an infinite line \((AB)\), through \(A\) to \(B\).

is_left(P, A, B, isclose=math.isclose)

Returns True if \(P\) is left, False if \(P\) is right, and raises a ValueError if \(P\) is on line AB. comparer is a function to test floating point equality.

A = ga.Point2D(1, 1, name='A')
B = ga.Point2D(4, 4, name='B')
P = ga.Point2D(2, 3, name='P')
Q = ga.Point2D(3, 2, name='Q')

AB = ga.Segment2D(A, B)

for obj in [A, B, AB]:
    obj.plot()

for obj in [P, Q]:
    obj.plot(color='red')

print('P is left AB:', P.is_left(A, B))
print('P is left BA:', P.is_left(B, A))
print()
print('Q is left AB:', Q.is_left(A, B))
print('Q is left BA:', Q.is_left(B, A))

P is left AB: True
P is left BA: False

Q is left AB: False
Q is left BA: True

for obj in [A, B, AB]:
    obj.plot()

R = ga.Point2D(2, 2, name='R')
S = ga.Point2D(3, 3.1, name='S')
for obj in [R, S]:
    obj.plot(color='red')

try:
    R.is_left(B, A)
except ValueError as err:
    print(err)

print('S is left AB:', S.is_left(A, B))

try:
    S.is_left(A, B, isclose=lambda x,y: abs(x-y)<0.5)
except ValueError as err:
    print(err)

R in on line (AB)
S is left AB: True
S in on line (AB)

double is_left(CPoint2D* A, CPoint2D* B, CPoint2D* P)

The value returned is:

• Strictly negative if \(P\) is right of the line through \(A\) to \(B\).
• Strictly positive if \(P\) is left of the line through A to B.
• Zero if \(P\) is on the line \((AB)\).

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint2D A, B, P, R

A.x, A.y = 1, 1
B.x, B.y = 4, 4

P.x, P.y = 2, 3
R.x, R.y = 2, 2

printf('P is left AB: %.2f\n', ga.is_left(&A, &B, &P))
printf('P is left BA: %.2f\n', ga.is_left(&B, &A, &P))
printf('R is left AB: %.2f\n', ga.is_left(&A, &B, &R))

P is left AB: 3.00
P is left BA: -3.00
R is left AB: 0.00

Vectors

There are two data structures to represent a vector in two-dimensional space, of components \((x, y)\):

• the C structure CVector2D,
• the Python extension type Vector2D, a wrapper around CVector2D.

CVector2D C structure

A CVector2D has just two members, to store its components:

CVector2D

double x

double y

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u

u.x = 1
u.y = 2

printf('u: (%.1f, %.1f)\n', u.x, u.y)

u: (1.0, 2.0)

CVector2D can be allocated and destroyed with these functions:

CVector2D* new_vector2d()

void del_vector2d(CVector2D* V)

Note

This functions are rarely used, as C structure in GeomAlgo are usually stack-allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D* u

u = ga.new_vector2d()

u.x = 1
u.y = 2

printf('u: (%.1f, %.1f)\n', u.x, u.y)

ga.del_vector2d(u)

u: (1.0, 2.0)

Vector2D Python extension

A Vector2D take components as arguments.

class Vector2D(x, y)

Attributes:

cvector2d: CVector2D*

u = ga.Vector2D(1, 2)
print(u)

<Vector2D(1.0,2.0>

The wrapped C structure Vector2D.cvector2d is accessible only using Cython.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.Vector2D u
    ga.CVector2D* ptr

u = ga.Vector2D(1, 2)

ptr = u.cvector2d
printf("(%.1f, %.1f)\n", ptr.x, ptr.y)

(1.0, 2.0)

Get and set components

CVector2D members are accessed directly to get and set coordinates:

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u

u.x, u.y = 1, 2
printf("u: (%.1f, %.1f)\n", u.x, u.y)

u.x, u.y = 4, -5
printf("u: (%.1f, %.1f)\n", u.x, u.y)

u: (1.0, 2.0)
u: (4.0, -5.0)

Vector2D as two properties x and y to get/set coordinates from/to its underlying cvector2d attribute:

Vector2D.x

Vector2D.y

u = ga.Vector2D(1, 2)
print(u)

u.x, u.y = 4, -5
print(u)

<Vector2D(1.0,2.0>
<Vector2D(4.0,-5.0>

Vector operators

Vectors addition

Compute the vector \(\mathbf{AC} = \mathbf{AB} + \mathbf{BC}\)

Vector2D.__add__(self, other)

AB = ga.Vector2D(0, 1)
BC = ga.Vector2D(1, 1)

AC = AB + BC
print(AC)

<Vector2D(1.0,2.0>

void add_vector2d(CVector2D *AC, CVector2D *AB, CVector2D *BC)

AC must be already allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D AB, BC, AC

AB.x, AB.y = 0, 1
BC.x, BC.y = 1, 1

ga.add_vector2d(&AC, &AB, &BC)
printf("(%.1f, %.1f)\n", AC.x, AC.y)

(1.0, 2.0)

Vectors substraction

Compute the vector \(\mathbf{AB} = \mathbf{AC} - \mathbf{BC}\)

Vector2D.__sub__(self, other)

AC = ga.Vector2D(1, 2)
BC = ga.Vector2D(1, 1)

AB = AC - BC
print(AB)

<Vector2D(0.0,1.0>

void subtract_vector2d(CVector2D *AB, CVector2D *AC, CVector2D *BC)

AB must be already allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D AB, BC, AC

AC.x, AC.y = 1, 2
BC.x, BC.y = 1, 1

ga.subtract_vector2d(&AB, &AC, &BC)
printf("(%.1f, %.1f)\n", AB.x, AB.y)

(0.0, 1.0)

Vector multiplication

Compute the vector \(\mathbf{t} = \alpha \mathbf{u}\)

Point2D.__mul__(self, alpha)

u = ga.Vector2D(1, 2)
t = u*2
print(t)

<Vector2D(2.0,4.0>

void vector2d_times_scalar(CVector2D *t, double alpha, CVector2D *u)

t must be already allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u, t

u.x, u.y = 1, 2
ga.vector2d_times_scalar(&t, 2., &u)
printf("(%.1f, %.1f)", t.x, t.y)

(2.0, 4.0)

Vectors dot product

Compute the dot product \(u_x v_x + u_y v_y\) between vectors \(\mathbf{u}\) and \(\mathbf{v}\).

Point2D.dot(self, other)

u = ga.Vector2D(1, 2)
v = ga.Vector2D(3, 4)
u.dot(v)

11.0

double dot_product2d(CVector2D *u, CVector2D *v)

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u, v

u.x, u.y = 1, 2
v.x, v.y = 3, 4

printf("%.1f\n", ga.dot_product2d(&u, &v))

11.0

Vectors cross product

Compute the cross product \(u_x v_y - u_y v_x\) between vectors \(\mathbf{u}\) and \(\mathbf{v}\).

Point2D.__xor__(self, other)

u = ga.Vector2D(1, 2)
v = ga.Vector2D(3, 4)
u^v

-2.0

double cross_product2d(CVector2D *u, CVector2D *v)

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u, v

u.x, u.y = 1, 2
v.x, v.y = 3, 4

printf("%.1f\n", ga.cross_product2d(&u, &v))

-2.0

Norms

Compute the norm \(|\mathbf{u}| = \sqrt{u_x^2 + u_y^2}\).

Vector2D.norm

u = ga.Vector2D(3, 4)
u.norm

5.0

double compute_norm2d(CVector2D *u)

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u

u.x, u.y = 3, 4
printf("%.1f\n", ga.compute_norm2d(&u))

5.0

Normalize vector \(\mathbf{u}\) so that its norm \(|\mathbf{u}|\) is 1.

Point2D.normalize(self)

u = ga.Vector2D(3, 4)
u.normalize()
print(u)
print(u.norm)

<Vector2D(0.6,0.8>
1.0

void normalize_vector2d(CVector2D *u)

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u

u.x, u.y = 3, 4
ga.normalize_vector2d(&u)
printf("(%.1f, %.1f)\n", u.x, u.y)
printf("%.1f\n", ga.compute_norm2d(&u))

(0.6, 0.8)
1.0

Normals

Compute the unitary (\(|\mathbf{u}=1|\) ) normal \(\mathbf{u}\) of vector \(\mathbf{u}\).

Point2D.normal

(property)

u = ga.Vector2D(0, 1)
print(u.normal)

<Vector2D(1.0,-0.0>

void compute_normal2d(CVector2D *n, CVector2D *u, double norm)

norm argument can be previously obtained with compute_norm2d()

n must be already allocated.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CVector2D u, n
    double norm

u.x, u.y = 0, 1

norm = ga.compute_norm2d(&u)
ga.compute_normal2d(&n, &u, norm)
printf("(%.1f, %.1f)\n", n.x, n.y)

(1.0, -0.0)

Segments

Contruction

class Segment2D(A, B)

Create a segment \([AB]\), where A and B can be either Point2D or tuple (x,y).

Read-write attributes

A: Point2D

B: Point2D

Segment Points. They may be replaced by another points, for example:

AB.A = Point2D(1, 2)

or coordinate may be changed inplace, for example:

AB.A.y = -2

Read only attributes

area: length

Segment length \(|AB|\)

Methods

In either case, the method recompute() must be called after change(s).

recompute(self)

Recompute segment property (length, vector). Must be called when a segment point is changed, or its coordinate is chanded.

Structures

CSegment2D

CPoint2D* A

\(A\) is the first segment point.

CPoint2D* B

\(B\) is the second segment point.

Note

CSegment2D represents a segment by two points. A segment may also be represented by a point and a vector. There is no C structure that represent a segment this way. If needed by a function, a Point2D and a Vector2D must be passed explicitly. See for example segment2d_where.

CSegment2D* new_segment2d()

Allocate a new CSegment2D.

This do not allocate members A, B.

void del_segment2d(CSegment2D* csegment2d)

Delete a CSegment2D.

This do not delete members A, B.

void segment2d_set(CSegment2D* AB, CPoint2D* A, CPoint2D* B)

Set members A and B.

Computational functions

double segment2d_distance_point2d(CSegment2D* AB, CVector2D* u, CPoint2D* P)

Compute distance of a point \(P\) to the segment \([AB]\).

\(\mathbf{u}\) is the vector from \(A\) to \(A\).

double segment2d_square_distance_point2d(CSegment2D* AB, CVector2D* u, CPoint2D* P)

Compute distance of a point \(P\) to the segment \([AB]\).

\(\mathbf{u}\) is the vector from \(A\) to \(A\).

Sometime, the knowledge of the square distance is enough, and for performance, computing the square root can be avoided.

double segment2d_where(CPoint2D* A, CVector2D* AB, CPoint2D* P)

Compute \(\alpha\) such as \(P = A + \alpha \mathbf{AB}\).

This assumes point \(P\) in on line \((AB)\).

void segment2d_middle(CPoint2D* M, CSegment2D* AB)

Compute the point \(M\), middle of the segment \([AB]\).

Triangles

Contruction

class Triangle2D(A, B, C, index=0, force_counterclockwise=False)

Create a triangle \(ABC\), where A, B and C are Point2D.

index is an optional integer index for this triangle.

If force_counterclockwise is True, points A, B and C may be reordered so that triangle is counterclockwise

Read-write attributes

A: Point2D

B: Point2D

C: Point2D

Triangle points. They may be replaced by another points, for example:

ABC.A = Point2D(1, 2)

or coordinate may be changed inplace, for example:

ABC.A.y = -2

In either case, the method recompute() must be called after change(s).

Read only attributes

area: float

Triangle area.

center: Point2D

Triangle center.

counterclockwise: bool

Whether triangle is counterclockwise.

Methods

recompute(self)

Recompute triangle properties (signed area, ...). Must be called when a triangle point is changed, or its coordinate is chanded.

includes_point(self, P, edge_width=0.)

Test if the triangle includes (contains) the point \(P\).

If point \(P\) is on one of the edge of triangle, due to numerical accuracy issues, the test may failed.

To solve this issue, if edge_width_square is not 0, triangle will be considered to include \(P\) if distance between \(P\) and one of triangle edge is less than edge_width (the root square of edge_width_square).

interpolate(self, data: double[3], P: Point2D)

Interpolate data defined on triangle vertices, to a point P inside the triangle.

plot(self, style='b-')

Plot triangle points using matplotlib

Structures

CTriangle2D

CPoint2D* A

\(A\) is the first triangle vertice.

CPoint2D* B

\(B\) is the second triangle vertice.

CPoint2D* C

\(C\) is the third triangle vertice.

CTriangle2D* new_triangle2d()

Allocate a new CTriangle2D.

This do not allocate members A, B, and C.

void del_triangle2d(CTriangle2D* ctri2d)

Delete a CTriangle2D.

This do not delete members A, B, and C.

void triangle2d_set(CTriangle2D* ABC, CPoint2D* A, CPoint2D* B, CPoint2D* C)

Set triangle points \(A\), \(B\), \(C\).

Computational functions

bint triangle2d_includes_point2d(CTriangle2D* ABC, CPoint2D* P, double edge_width_square)

Test if the triangle \(ABC\) includes (contains) the point \(P\).

If point \(P\) is on one of the edge of triangle \(ABC\), due to numerical accuracy issues, the test may failed.

To solve this issue, if edge_width_square is not 0, ABC will be considered to include \(P\) if distance between \(P\) and one of \(ABC\) edge is less than edge_width (the root square of edge_width_square).

int triangle2d_on_edges(CTriangle2D* ABC, CPoint2D* P, double edge_width_square)

Test if point \(P\) is on one of the edge of triangle \(ABC\).

\(P\) will be considered to be on one of the edge of triangle \(ABC\). if distance between \(P\) and one of \(ABC\) edge is less than edge_width (the root square of edge_width_square).

Returns:

• 0 if \(P\) is on edge \([AB]\),
• 1 if \(P\) is on edge \([BC]\),
• 2 if \(P\) is on edge \([CA]\),
• -2 if \(P\) is not on any edge.

double triangle2d_signed_area(CTriangle2D* T)

Compute the signed area of triangle \(T\).

Triangle area is the absolute value of the signed area.

• If signed area is positive, triangle is counterclockwise.
• If signed area is negative, triangle is clockwise.
• If signed area is zero, triangle is degenerated.

double triangle2d_area(CTriangle2D* T)

Compute the area of triangle \(T\).

void triangle2d_center(CTriangle2D* T, CPoint2D* C)

Compute the center \(C\) of triangle \(T\).

Variable C must be already allocated.

void triangle2d_gradx_grady_det(CTriangle2D* tri, double signed_area, double gradx[3], double grady[3], double det[3])

Compute factors for linear interpolation of data defined on triangle vertices \(A\), \(B\) and \(C\) to points included in the triangle.

Relations between Python and Cython API

Triangulation

Introduction

x, y and trivtx

A triangulation is a subdivision of a two-dimensional domain into triangles, defined by:

• its vertices: two 1D arrays x and y for coordiantes,
• its triangles: one 2D array trivtx, giving the 3 vertice indices of each triangle.

x = array([0, 1, 2, 0, 1, 2, 0, 1], dtype='d')
y = array([0, 0, 0, 1, 1, 1, 2, 2], dtype='d')
trivtx = array([[0, 1, 3], [1, 2, 4], [1, 4, 3], [2, 5, 4],
                [3, 4, 6], [4, 7, 6]], dtype='int32')

def plot_with_indices():
    # Plot triangulation.
    triplot(x, y, trivtx)

    # Plot vertices indices.
    for ivert, (xvert, yvert) in enumerate(zip(x, y)):
        text(xvert, yvert, ivert)

    # Plot triangle indices.
    for itri, (v0, v1, v2) in enumerate(trivtx):
        xcenter = (x[v0] + x[v1] + x[v2]) / 3.
        ycenter = (y[v0] + y[v1] + y[v2]) / 3.
        text(xcenter, ycenter, itri, color='red')

    axis('scaled')
    show()

plot_with_indices()

For example, the vertice indices of triangle 2 are:

print(trivtx[2])

[1 4 3]

ga.Triangulation2D

The main purpose of the Triangulation2D Python extension type is to group this arrays together.

TG = ga.Triangulation2D(x, y, trivtx)

Note

Point2D, Segment2D Triangle2D, etc, are objects that exist in a large number in a program. For example, a numerical simulation may create 10 millions of points. These points will be stored in SoA x and y. There is a Point2D extension type, and an associated CPoint2D C structure that can be stack-allocated for high performance computing on these million of points.

For Triangulation2D, there is typically one or two objects in a program. There do not need stack allocation, and as a such there is no need for C structure.

Iterating on triangle

To iterate a triangle, one first need toextract triangle vertice indices from trivtx, then vertice coordinates from x and y.

Triangulation2D.get(Triangulation2D self, int triangle_index, CTriangle2D* triangle)

The Triangulation2D.get() method make it a lot easier. It take as argument the triangle iteration index, and a c:type:CTriangle2D setted with tree c:type:CPoint2D.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

def iter_triangles(ga.Triangulation2D TG):
    cdef:
        int T
        ga.CTriangle2D ABC
        ga.CPoint2D A, B, C

    ga.triangle2d_set(&ABC, &A, &B, &C)

    for T in range(TG.NT):
        TG.get(T, &ABC)
        printf("%d: A(%.0f,%.0f), B(%.0f,%.0f), C(%.0f,%.0f)\n",
               T, A.x, A.y, B.x, B.y, C.x, C.y)

iter_triangles(TG)

0: A(0,0), B(1,0), C(0,1)
1: A(1,0), B(2,0), C(1,1)
2: A(1,0), B(1,1), C(0,1)
3: A(2,0), B(2,1), C(1,1)
4: A(0,1), B(1,1), C(0,2)
5: A(1,1), B(1,2), C(0,2)

Triangulation2D.__getitem__(self, triangle_index: int)

Return a Triangle2D instance

Computational functions

Triangle centers

triplot(x, y, trivtx)

xcenter, ycenter = ga.triangulation.compute_centers(TG)

plot(xcenter, ycenter, 'bo')

axis('scaled')
show()

bounding box

edge lengths

area

Note

Triangulation2D extension type is voluntary keept minimial with 3 attributes and no methods, like above ones could be. The rational behind this is to not impose the application that geomalgo a particular structure, but to let it keep is own: triangular mesh, dual mesh... geomalgo goal is to provide some useful but non-intruise computational primitives.

Boundary and intern edges

A lot of operation on triangulation, for example building a dual mesh, or refining the mesh can be made simple by dealing with boundary edges and intern edges.

build_edges(trivtx, NV, intern_edges_order=None, boundary_edges_order=None)

Parameters:

• trivtx (int[NT,3]) – Triangle vertice indices.
• NV (int) – Number of vertices.

It returns:

a BoundaryEdges object a InternEdges object a EdgeMap object

Boundary edges belongs to 1 triangle, while intern edges belongs to 2 triangles. Their are build with the build_edges().

intern_edges, boundary_edges, edge_map = ga.build_edges(TG.trivtx, TG.NV)

InternEdges has vertices and triangles attributes, giving the 2 vertices indices and the 2 triangle indices for each intern edge. For example:

triplot(x, y, trivtx)
I = 2

IA = intern_edges.vertices[I, 0]
IB = intern_edges.vertices[I, 1]
plot([x[IA], x[IB]],
     [y[IA], y[IB]], 'ro')

T0 = intern_edges.triangles[I, 0]
T1 = intern_edges.triangles[I, 1]
plot([xcenter[T0], xcenter[T1]],
     [ycenter[T0], ycenter[T1]], 'go')

axis('scaled')
show()

BoundaryEdges.size: int

Number of boundary edges in the triangulation.

BoundaryEdges.vertices: int[:,:]

BoundaryEdges.triangle: int[:]

BoundaryEdges.next_boundary_edge: int[:]

BoundaryEdges.edge_map: EdgeMap

BoundaryEdges has vertices and triangle attributes, giving the 2 vertices indices and the triangle index for each boundary edge. For example:

triplot(x, y, trivtx)
I = 2

IA = boundary_edges.vertices[I, 0]
IB = boundary_edges.vertices[I, 1]
plot([x[IA], x[IB]],
     [y[IA], y[IB]], 'ro')

T = boundary_edges.triangle[I]
plot(xcenter[T], ycenter[T], 'go')

axis('scaled')
show()

The third returned object, edge_map, is used to find a intern_edge or boundary_edge given its vertices. It is implemented to be very fast with O(1) complexity.

EdgeMap indexes both InternalEdge and BoundaryEdges. The map_search_edge_idx returns whether the edge is found, and it EdgeMap index. The EdgeMap index is then used to determine whether the edge is a boundary or internal one, and its corresponding index.

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

def find_edge(ga.EdgeMap edge_map, int IA, int IB):
    cdef:
        bint found
        int E, I

    E = edge_map.search_edge_idx(IA, IB, &found)

    if not found:
        printf("No such edge: (%d,%d)\n", IA, IB)
        return

    I = edge_map.idx[E]

    if edge_map.location[E] == ga.INTERN_EDGE:
        printf("Intern edge for (%d,%d): %d\n", IA, IB, I)
    else:
        printf("Boudnary edge for (%d,%d): %d\n", IA, IB, I)

def text_middle(IA, IB, color):
    A = ga.Point2D(x[IA], y[IA])
    B = ga.Point2D(x[IB], y[IB])
    AB = ga.Segment2D(A, B)
    M = AB.compute_middle()
    text(M.x, M.y, I, color=color, size=16)

for I in range(boundary_edges.size):
    IA, IB = boundary_edges.vertices[I]
    text_middle(IA, IB, 'green')

for I in range(intern_edges.size):
    IA, IB = intern_edges.vertices[I]
    text_middle(IA, IB, 'blue')

plot_with_indices()

find_edge(edge_map, 1, 4)
find_edge(edge_map, 7, 6)
find_edge(edge_map, 7, 5)

Intern edge for (1,4): 1
Boudnary edge for (7,6): 7
No such edge: (7,5)

Locate point

Introduction

The TriangulationLocator is used to locate in a triangulation which triangle contains a point. grid defines the grid size used to speed up the triangle searching, see Change grid size. edge_width is used for triange point inclusion tests, see Triangle2D.include_points(). bounds is an int[NT:4] array allocated to hold temporary value, an already allocated buffer may be provided.

class TriangulationLocator(self, TG, grid=None, edge_width=-1, bounds=None)

Attributes:

TG: Triangulation2D

grid: Grid2D

edge_width: double

edge_width_square: double

celltri: int[:]

celltri_idx: int[:]

Cell to triangle map in CSR format.

import triangle

def generate_triangulation():
    A = {'vertices' : array([[0,0], [1,0], [1, 1], [0, 1]])}
    B = triangle.triangulate(A, 'qa0.001')
    return ga.Triangulation2D(B['vertices'][:,0],
                              B['vertices'][:,1],
                              B['triangles'])
TG = generate_triangulation()

xpoints = np.array([0.2, 0.4, 0.6, 0.8])
ypoints = np.array([0.5, 0.5, 0.5, 0.5])

locator = ga.TriangulationLocator(TG)
triangles = locator.search_points(xpoints, ypoints)

fig = figure(figsize=(8,8))
ax = fig.add_subplot(1,1,1)
triplot(TG.x, TG.y, TG.trivtx, color='gray')
plot(xpoints, ypoints, 'ro')
for T in triangles:
    TG[T].plot()
axis('scaled')
show()

Change grid size

A grid of size \(nx \times ny\) (depending of the number of triangles) is used to map each cells to the triangles it contains. This allow the triangle search performance to be \(O(1)\):

• for a given point, the including cell is found (trivial computation)
• only the triangles in the cells are tested for point inclusion.

The counterpart is the memory usage expanse, as celltri and celltri_idx grows.

For example:

cell_numbers = array([[1,  2,  4], [8, 16, 32]])

nrows = cell_numbers.shape[0]
ncols = cell_numbers.shape[1]
fig, axes = subplots(ncols=3, nrows=2, figsize=(16,12))

for irow in range(nrows):
    for icol in range(ncols):
        n = cell_numbers[irow, icol]
        fig.sca(axes[irow, icol])

        triplot(TG.x, TG.y, TG.trivtx, color='gray', lw=0.5)

        grid = ga.triangulation.build_grid(TG, nx=n, ny=n)
        grid.plot(color='black', lw=0.5)
        locator = ga.TriangulationLocator(TG, grid)

        P = ga.Point2D(0.4, 0.4)
        P.plot(color='red')

        cell = locator.grid.find_cell(P)

        for T in locator.cell_to_triangles(cell.ix, cell.iy):
            TG[T].plot(lw=0.5)

        print("nx, ny = {}, celltri size: {}, celltri_idx size: {}".format(
                  n, len(locator.celltri), len(locator.celltri_idx)))

        axis('scaled')
show()

nx, ny = 1, celltri size: 1539, celltri_idx size: 2
nx, ny = 2, celltri size: 1679, celltri_idx size: 5
nx, ny = 4, celltri size: 1945, celltri_idx size: 17
nx, ny = 8, celltri size: 2513, celltri_idx size: 65
nx, ny = 16, celltri size: 3909, celltri_idx size: 257
nx, ny = 32, celltri size: 7549, celltri_idx size: 1025

Interpolator

TriangulationInterpolator interpolates data defined on vertices of a triangulation to given points. NP is the number of points to interpolate to. gradx, grady and det are double[NT, 3] arrays, computed with interpolation.compute_interpolator(), which will be invoked if one the three value is None.

class TriangulationInterpolator(self, TG, locator, NP, gradx=None, grady=None, det=None)

Attributes:

locator: TriangulationLocator

TG: Triangulation2D

NP: int

gradx: double[NT, 3]

grady: double[NT, 3]

det: double[NT, 3]

triangles: int[NP]

factors: double[NT, 3]

Points positions are with with TriangulationInterpolator.set_points() which take a inputs xpoints and ypoints, two double[NP] arrays, and returns the number of points find out of the triangulation.

TriangulationInterpolator.set_points(self, xpoints, ypoints)

Once the points are set, TriangulationInterpolator.interpolate() is called to interpolate vertdata (double[NV]) to pointdata (``double[NP]). Of course, it can be called many times of different verdata values.

TriangulationInterpolator.interpolate(self, vertdata, pointdata)

def f(x,y):
    x, y  = np.asarray(x), np.asarray(y)
    return 3*x**2 - 2*y + 2


def g(x,y):
    x, y  = np.asarray(x), np.asarray(y)
    return -x**2 + 0.5*y


# Plot data.
vertdata = f(TG.x, TG.y)
tripcolor(TG.x, TG.y, TG.trivtx, vertdata)
axis('scaled')

# Plot points to interpolate on.
NP = 15
xline = linspace(0, 1, NP)
yline = full(NP, fill_value=0.5)
plot(xline, yline, 'ko', ms=1)

[<matplotlib.lines.Line2D at 0x7fecc5831f98>]

# Interpolate
locator = ga.TriangulationLocator(TG)
interpolator = ga.TriangulationInterpolator(TG, locator, NP)
interpolator.set_points(xline, yline)
linedata = np.empty(NP)
interpolator.interpolate(vertdata, linedata)

# Plot interpolation
plot(xline, linedata, label='f at y=0.5')

# Interpolate other data at the same points
vertdata = g(TG.x, TG.y)
interpolator.interpolate(vertdata, linedata)
plot(xline, linedata, label='g at y=0.5')

# Change points, and reinterpolate
yline[:] = 0.75
interpolator.set_points(xline, yline)
interpolator.interpolate(vertdata, linedata)
plot(xline, linedata, label='g at y=0.75')

legend()

<matplotlib.legend.Legend at 0x7fecc6155748>

Note

Changing points position imply re-locating triangle and re-computing interpolation factors.

Note

The example interpolate on a line, but there no assumption about it. Points can be unstructured, on a circle, on a cartesian grid, etc.

Triangle 3D

There are two data structures to represent a triangle in three-dimensinal space, of points \((A, B, C)\):

• the C structure CTriangle3D,
• the Python extension type Triangle3D, a wrapper around CTriangle3D.

CTriangle3D C structure

A CTriangle3D has just three members, pointing to triangle vertices:

CTriangle3D

CPoint3D* A

CPoint3D* B

CPoint3D* C

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.CPoint3D A, B, C
    ga.CTriangle3D ABC

A.x, A.y, A.z = 0, 0, 0
B.x, B.y, B.z = 1, 1, 0
C.x, C.y, C.z = 1, 1, 0

ABC.A = &A
ABC.B = &B
ABC.C = &C

printf("ABC.B.x: %.1f\n", ABC.B.x)

ABC.B.x: 1.0

Triangle3D Python extension type

A Triangle3D takes the three vertices A, B and C as arguments (Point3D), and an optional index.

class Triangle3D(A, B, C, index=0)

Attributes:

ctri3d: CTriangle3D

index: int

A = ga.Point3D(0, 0, 0, name='A')
B = ga.Point3D(1, 1, 0, name='B')
C = ga.Point3D(0, 0, 1, name='C')

ABC = ga.Triangle3D(A, B, C)

fig = plt.figure(figsize=(8,8))

for obj in [ABC, A, B, C]:
    obj.plot()

ax = gca(projection='3d')
ax.set_aspect('equal')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

<matplotlib.text.Text at 0x7f129317a358>

the wrapped C structure Triangle3D.ctri3D is accessible only using Cython

%%cython

from libc.stdio cimport printf
cimport geomalgo as ga

cdef:
    ga.Point3D A = ga.Point3D(0, 0, 0)
    ga.Point3D B = ga.Point3D(1, 1, 0)
    ga.Point3D C = ga.Point3D(0, 0, 1)
    ga.Triangle3D ABC = ga.Triangle3D(A, B, C)
    ga.CTriangle3D* ptr

ABC = ga.Triangle3D(A, B, C)

ptr = &ABC.ctri3d
printf('ABC.B.x: %.1f\n', ptr.B.x)

ABC.B.x: 1.0

Intersections with segments

P = ga.Point3D(0.2, 0., 0.4, name='P')
Q = ga.Point3D(0.2, 1., 0.4, name='Q')
PQ = ga.Segment3D(P, Q)

I = ga.intersec3d_triangle_segment(ABC, PQ)
I.name = 'I'

fig = plt.figure(figsize=(8,8))

for obj in [ABC, A, B, C, PQ, P, Q]:
    obj.plot()

I.plot(color='r')

ax = gca(projection='3d')
ax.set_aspect('equal')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

<matplotlib.text.Text at 0x7f1290621550>

Developer’s Guide

Build

conda install python=3.5 cython craftr

Test

conda install nose

Documentation

conda install matplotlib jupyter nbsphinx sphinx_rtd_theme sphinx-gallery
pip install wurlitzer

Release

Update CHANGELOG file.

Bump version number in files:

• conda-recipe/meta.yaml
• setup.py
• geomalgo/__init__.py

Commit and tag:

git commit -m ‘bump to version X.Y.Z’ git tag X.Y.Z git push –tags

Change version number to X.Y.<Z+1>dev in files:

• conda-recipe/meta.yaml
• setup.py
• geomalgo/__init__.py

git commit -m ‘change version to X.Y.<Z+1>dev’

Warning

Changing version number to X.Y.<Z+1>dev is important, because on each commit, conda packages are built on Travis, and uploaded to anaconda.org.

If version number remains to X.Y.Z, development package will erase tagged package on anaconda.org

Development and contact

geomalgo is developed on GitHub, were issues and pull requests can be made. Do not hesitate to send me an email at david.froger@mailoo.org .