Reference: http://www.flegrand.fr/scidoc/docmml/graphie/geometrie/polygone/polygone.html
(in french)
reference herein : M. de Berg, O. Cheong, M. van Kreveld, M. Overmars,
Computational geometry, (Springer, 2010)
Direct translation of the Frederic Legrand's Python code
Meaning of the different french terms
voisin_sommet (n, i, di)
let L the list [1, ..., n] where n is the number of vertices
This procedure returns the index of the neighbour (voisin) of the vertex (sommet) i when L is rotated by di
equation_droite (P0, P1, M)
Let P0 and P1 two vertices and M an arbitrary point.
Let (P0, P1) the vector originated at P0 and ending at P1 (idem for (P0, M)) and the unitary vector in the Z direction.
This procedure returns (P0, P1) o (P0, M) .
point_dans_triangle (triangle, M) P1, P2]
This procedure returns "true" if point M is within (strictly) the triangle "triangle" and "false" if not.
sommet_distance_maximale (polygone, P0, P1, P2, indices)
Given a polygon (polygone) threes vertices P0, P1 and P2 and a list of indices , this procedure returns
the vertex of the polygon "polygone" which verifies: 1/ this vertex is strictly within
the triangle [P0, P1, P2] and 2/ it is the farthest from side [P1, P2] amid all the vertices that verifies point 1/.
If there is no such vertex the procedure returns NULL.
sommet_gauche (polygone)
This procedure returns the index of the leftmost ("gauche" means "left") vertex in polygon "polygone".
If more than one vertices have the same minimum abscissa value then only the first one is returned.
nouveau_polygone(polygone,i_debut,i_fin)
This procedure creates a new polygon from index i_debut (could be translated by i_first) to i_end (i_last)
trianguler_polygone_recursif(polygone)
This procedure recursively divides a polygon in two parts A and B from its leftmost vertex.
If A (B) is a triangle the list "liste_triangles" (mening "list of triangles") is augmented by A (B);
if not the procedure executes recursively on A and B.
trianguler_polygone(polygone)
This procedure triangulates the polygon "polygon"
hachurer(shapes, hatch_angle, hatch_number, hatch_color)
This procedure generates stes of hatches of different angles, colors and numbers
Limitations:
1/ "polygone" is a simply connected polygon
2/ two different sides S and S', either do not intersect or either have a common vertex
> 
trianguler := module()
export voisin_sommet, equation_droite, interieur_forme, point_dans_triangle, sommet_distance_maximale,
sommet_gauche, nouveau_polygone, trianguler_polygone_recursif, trianguler_polygone, hachurer:
#
voisin_sommet := (n, i, di) > ListTools:Rotate([$1..n], di)[i]:
#
equation_droite := proc(P0, P1, M) (P1[1]P0[1])*(M[2]P0[2])  (P1[2]P0[2])*(M[1]P0[1]) end proc:
#
interieur_forme := proc(forme, M)
local N:
N := numelems(forme);
{ seq( equation_droite(forme[n], forme[piecewise(n=N, 1, n+1)], M) >= 0, n=1..N) }
end proc:
#
point_dans_triangle := proc(triangle, M)
`and`(
is( equation_droite(triangle[1], triangle[2], M) > 0 ),
is( equation_droite(triangle[2], triangle[3], M) > 0 ),
is( equation_droite(triangle[3], triangle[1], M) > 0 )
)
end proc:
#
sommet_distance_maximale := proc(polygone, P0, P1, P2, indices)
local n, distance, j, i, M, d;
n := numelems(polygone):
distance := 0:
j := NULL:
for i from 1 to n do
if `not`(member(i, indices)) then
M := polygone[i];
if point_dans_triangle([P0, P1, P2], M) then
d := abs(equation_droite(P1, P2, M)):
if d > distance then
distance := d:
j := i
end if:
end if:
end if:
end do:
return j:
end proc:
#
sommet_gauche := polygone > sort(polygone, key=(x>x[1]), output=permutation)[1]:
#
nouveau_polygone := proc(polygone, i_debut, i_fin)
local n, p, i:
n := numelems(polygone):
p := NULL:
i := i_debut:
while i <> i_fin do
p := p, polygone[i]:
i := voisin_sommet(n, i, 1)
end do:
p := [p, polygone[i_fin]]
end proc:
#
trianguler_polygone_recursif := proc(polygone)
local n, j0, j1, j2, P0, P1, P2, j, polygone_1, polygone_2:
global liste_triangles:
n := numelems(polygone):
j0 := sommet_gauche(polygone):
j1 := voisin_sommet(n, j0, +1):
j2 := voisin_sommet(n, j0, 1):
P0 := polygone[j0]:
P1 := polygone[j1]:
P2 := polygone[j2]:
j := sommet_distance_maximale(polygone, P0, P1, P2, [j0, j1, j2]):
if `not`(j::posint) then
liste_triangles := liste_triangles, [P0, P1, P2]:
polygone_1 := nouveau_polygone(polygone,j1,j2):
if numelems(polygone_1) = 3 then
liste_triangles := liste_triangles, polygone_1:
else
thisproc(polygone_1)
end if:
else
polygone_1 := nouveau_polygone(polygone, j0, j ):
polygone_2 := nouveau_polygone(polygone, j , j0):
if numelems(polygone_1) = 3 then
liste_triangles := liste_triangles, polygone_1:
else
thisproc(polygone_1)
end if:
if numelems(polygone_2) = 3 then
liste_triangles := liste_triangles, polygone_2:
else
thisproc(polygone_2)
end if:
end if:
return [liste_triangles]:
end proc:
#
trianguler_polygone := proc(polygone)
trianguler_polygone_recursif(polygone):
return liste_triangles:
end proc:
#
hachurer := proc(shapes, hatch_angle::list, hatch_number::list, hatch_color::list)
local A, La, Lp;
local N, P, _sides, L_sides, Xshape, ch, rel, p_rel, n, sol, p_range:
local AllHatches, window, p, _segment:
local NT, ka, N_Hatches, p_range_t, nt, shape, p_hatches, WhatHatches:
#
# internal functions:
#
# La(x, y, alpha, p) is the implicit equation of a straight line of angle alpha relatively
# to the horizontal axis and intercept p
#
# Lp(x, y, P) is the implicit equation of a straight line passing through points P[1] and P[2]
#
# interieur_triangle(triangle, M)
La := (x, y, alpha, p) > cos(alpha)*x  sin(alpha)*y + p;
Lp := proc(x, y, P::list) (xP[1][1])*(P[2][2]P[1][2])  (yP[1][2])*(P[2][1]P[1][1] = 0) end proc;
p_range := [+infinity, infinity]:
NT := numelems(shapes):
AllHatches := NULL:
for ka from 1 to numelems(hatch_angle) do
A := hatch_angle[ka]:
N_Hatches := hatch_number[ka]:
p_range_t := NULL:
_sides := []:
L_sides := []:
rel := []:
for nt from 1 to NT do
shape := shapes[nt]:
# _sides : two points description of the sides of the shape
# L_sides : implicit equations of the straight lines that support the sides
N := [1, 2, 3];
P := [2, 3, 1];
_sides := [ _sides[] , [ seq([shape[n], shape[P[n]]], n in N) ] ];
L_sides := [ L_sides[], Lp~(x, y, _sides[1]) ];
# Inequalities that define the interior of the shape
rel := [ rel[], trianguler:interieur_forme(shape, [x, y]) ];
# Given the orientation of the hatches we search here the extreme values of
# the intercept p for wich a straight line of equation La(x, y, alpha, p)
# cuts the shape.
p_rel := NULL:
for n from 1 to numelems(L_sides[1]) do
sol := solve({La(x, y, A, q), lhs(L_sides[1][n])} union rel[1], [x, y]);
p_rel := p_rel, `if`(sol <> [], [rhs(op(1, %)), rhs(op(3, %))], [+infinity, infinity]);
end do:
p_range_t := p_range_t, evalf(min(op~(1, [p_rel]))..max(op~(2, [p_rel])));
p_range := evalf(min(op~(1, [p_rel]), op(1, p_range))..max(op~(2, [p_rel]), op(2, p_range)));
end do: # end of the loop over triangles
p_range_t := [p_range_t]:
p_hatches := [seq(p_range, (op(2, p_range)op(1, p_range))/N_Hatches)]:
# Building of the hatches
#
# This construction is far from being optimal.
# Here again the main goal was to obtain the hatches with a minimum effort
# if algorithmic development.
window := min(op~(1..shape))..max(op~(1..shape)):
WhatHatches := map(v > map(u > if verify(u, v, 'interval'('closed') ) then u end if, p_hatches), p_range_t):
for nt from 1 to NT do
for p in WhatHatches[nt] do
_segment := []:
for n from 1 to numelems(L_sides[nt]) do
_segment := _segment, evalf( solve({La(x, y, A, p), lhs(L_sides[nt][n])} union rel[nt], [x, y]) );
end do;
map(u > u[], [_segment]);
AllHatches := AllHatches, plot(map(u > rhs~(u), %), color=hatch_color[ka]):
end do:
end do;
end do: # end of loop over hatch angles
plots:display(
PLOT(POLYGONS(polygone, COLOR(RGB, 1$3))),
AllHatches,
scaling=constrained
)
end proc:
end module:

Legrand's example (see reference above)
> 
global liste_triangles:
liste_triangles := NULL:

> 
polygone := [[0,0],[0.5,1],[1.5,0.2],[2,0.5],[2,0],[1.5,1],[0.3,0],[0.5,1]]:
trianguler:trianguler_polygone(polygone):
PLOT(seq(POLYGONS(u, COLOR(RGB, rand()/10^12, rand()/10^12, rand()/10^12)), u in liste_triangles), VIEW(0..2, 2..2))

> 
trianguler:hachurer([liste_triangles], [Pi/4, Pi/4], [40, 40], [red, blue])

> 
F := (P, a, b) > map(p > [p[1]+a, p[2]+b], P):

> 
MOTIF := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]];
motifs := [ seq(seq(F(MOTIF, 0+i*0.075, 0+j*0.075), i=0..26), j=14..13) ]:
plots:display(
plot([polygone[], polygone[1]], color=red, filled),
map(u > plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):
texture := NULL:
rel_motifs := map(u > trianguler:interieur_forme(u, [x, y]), motifs):
for ref in liste_triangles do
ref;
#
# the three lines below are used to define REF counter clockwise
#
g := trianguler:sommet_gauche(ref):
bas := sort(op~(2, ref), output=permutation);
REF := ref[[g, op(map(u > if u<>g then u end if, bas))]];
rel_ref := trianguler:interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
texture_ref := map(u > plots:inequal(rel_ref union u, x=0..2, y=1..1, color=blue, 'nolines'), rel_motifs):
texture := texture, texture_ref:
end do:
plots:display(
plot([polygone[], polygone[1]], color=red, scaling=constrained),
texture
)





> 
MOTIF := [[0, 0], [0.05, 0], [0.05, 0.05], [0, 0.05]];
motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.05, 0.1)+i*0.1, 0+j*0.05), i=0.2..20), j=20..20) ]:
plots:display(
plot([polygone[], polygone[1]], color=red, filled),
map(u > plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):
texture := NULL:
rel_motifs := map(u > trianguler:interieur_forme(u, [x, y]), motifs):
for ref in liste_triangles do
ref;
g := trianguler:sommet_gauche(ref):
bas := sort(op~(2, ref), output=permutation);
REF := ref[[g, op(map(u > if u<>g then u end if, bas))]];
rel_ref := trianguler:interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
texture_ref := map(u > plots:inequal(rel_ref union u, x=0..2, y=1..1, color=blue, 'nolines'), rel_motifs):
texture := texture, texture_ref:
end do:
plots:display(
plot([polygone[], polygone[1]], color=red, scaling=constrained),
texture
)





> 
MOTIF := [[0, 0], [0.4, 0], [0.4, 0.14], [0, 0.14]]:
motifs := [ seq(seq(F(MOTIF, piecewise(j::odd, 0.4, 0.2)+i*0.4, 0+j*0.14), i=1..4), j=8..7) ]:
plots:display(
plot([polygone[], polygone[1]], color=red, filled),
map(u > plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):
palettes := ColorTools:PaletteNames():
ColorTools:GetPalette("HTML"):
couleurs := [SandyBrown, PeachPuff, Peru, Linen, Bisque, Burlywood, Tan, AntiqueWhite, NavajoWhite, BlanchedAlmond, PapayaWhip, Moccasin, Wheat]:
nc := numelems(couleurs):
roll := rand(1..nc):
motifs_nb := numelems(motifs):
motifs_couleur := [ seq(cat("HTML ", couleurs[roll()]), n=1..motifs_nb) ]:
texture := NULL:
rel_motifs := map(u > trianguler:interieur_forme(u, [x, y]), motifs):
for ref in liste_triangles do
ref;
g := trianguler:sommet_gauche(ref):
bas := sort(op~(2, ref), output=permutation);
REF := ref[[g, op(map(u > if u<>g then u end if, bas))]];
rel_ref := trianguler:interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
texture_ref := map(n > plots:inequal(rel_ref union rel_motifs[n], x=0..2, y=1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]):
texture := texture, texture_ref:
end do:
plots:display(
plot([polygone[], polygone[1]], color=red, scaling=constrained),
texture
)



> 
MOTIF := [ seq(0.1*~[cos(Pi/6+Pi/3*i), sin(Pi/6+Pi/3*i)], i=0..5) ]:
motifs := [ seq(seq(F(MOTIF, i*0.2*cos(Pi/6)+piecewise(j::odd, 0, 0.08), j*0.3*sin(Pi/6)), i=0..12), j=6..6) ]:
plots:display(
plot([polygone[], polygone[1]], color=red, filled),
map(u > plot([u[], u[1]], color=blue, filled, scaling=constrained), motifs)
):
motifs_nb := numelems(motifs):
motifs_couleur := [ seq(`if`(n::even, yellow, brown) , n=1..motifs_nb) ]:
texture := NULL:
rel_motifs := map(u > trianguler:interieur_forme(u, [x, y]), motifs):
for ref in liste_triangles do
ref;
g := trianguler:sommet_gauche(ref):
bas := sort(op~(2, ref), output=permutation);
REF := ref[[g, op(map(u > if u<>g then u end if, bas))]];
rel_ref := trianguler:interieur_forme(REF, [x, y]): #print(ref, REF, rel_ref);
texture_ref := map(n > plots:inequal(rel_ref union rel_motifs[n], x=0..2, y=1..1, color=motifs_couleur[n], 'nolines'), [$1..motifs_nb]):
texture := texture, texture_ref:
end do:
plots:display(
plot([polygone[], polygone[1]], color=red, scaling=constrained),
texture
)



