9.8 A Database of Shapes

Until now, we have used the CAD tool just as a visualizer of the shapes that our programs generate. We will now see that a CAD tool is not only a drawing program, but also a database of geometric figures. In fact, every time we draw something, the CAD tool records the created graphic entity in its database.

There are several ways to access the created entities. One of the simplest ways for a CAD tool user is to use the mouse to click on the graphic entity to which he wants to access. Another way, more useful for those who want to program, will be by calling Julia’s functions that return, as results, the existing geometrical entities. There are several of those functions at our disposal but, for now, let us look at one of the simplest: the function all_shapes.

The all_shapes function does not receive any arguments, and it returns an array of all existing geometric entities in the CAD tool. Naturally, if there are no entities, the function returns an empty array.

The following interaction demonstrates the behavior of this function:

> all_shapes()

Shape[]

> circle(xy(1, 2), 3)

Circle(...)

> sphere(xyz(1, 2, 3), 4)

Sphere(...)

> all_shapes()

2-element Vector{Shape}:

Circle(...)

Sphere(...)

As we can see from the previous interaction, the function all_shapes creates representations of the geometric entities existent in the CAD tool without having any idea of how these entities got there. When necessary (and possible), these forms can be identified by the recognizers is_point, is_circle, is_line, is_closed_line, is_spline, is_closed_spline, is_surface, among others.

Given a geometric entity, we may be interested in knowing its properties. The access to these properties depends much on what the CAD tool is able to provide. For example, for a circle, we can know its center through the function circle_center and its radius through the function circle_radius. For a point, the function point_position returns its position and for a polygonal line, the function line_vertices produces an array with the positions of the vertices of that line. However, for a generic solid, we have no access to any property.

In summary, we have:

point_position(point(p)) \(=\) p

circle_center(circle(p, r)) \(=\) p

circle_center(circle(p, r)) \(=\) r

line_vertices(line(p1, p2, ..., pn)) \(=\) [p1, p2, ..., pn]

Although the above equivalences show the relations between the constructors of geometric entities (point, circle, etc.) and the selectors of those entities (point_position, circle_center, etc.), it is important to acknowledge that it is possible to use the selectors with entities that were not created from Khepri, but directly in the CAD tool. This can be particularly useful when we want to write programs that, instead of generating geometry, only process existing geometry.

9.8.1 City Models from City Plans

We will now exemplify the use of these operations in solving a real-world problem: the creation of three-dimensional models of buildings from two-dimensional plans provided by town hall services. These plans are characterized by representing each building with its delimiting polygon, containing, in its interior, a point at the building’s height. This figure shows a fragment of one of such plans, where the points are displayed by icons large enough to make them visible.

Plan supplied by town hall services. Each building is represented by its delimiting polygon and a point at the building’s height. Due to errors in plans, some buildings may not have any height point, while others may have more than one height point. There may also be height points that are not associated with any building.

Unfortunately, it is not uncommon to find plans with errors and, the example shown in this figure illustrates precisely this: a careful observation will reveal buildings that do not possess any height point, as well as height points that are not associated with any building. Furthermore, it is possible to find buildings with more than one height point. These are problems that we have to take into account when we think about creating programs that manipulate these plans.

In order to create a three-dimensional model from one of these plans we can, for each polygon, create a region that we will extrude up to the height indicated by the corresponding point. To do so, we have to be able to identify, for each polygon, which point that is (or what points those are in case of there is more than one).

Detecting if a point is contained inside a polygon is a classic problem of computational geometry, for which there are various solutions. One of the simplest is to draw a ray from the point and count the number of intersections that this ray has with the edges of the polygon, as illustrated in this figure with several points. If the number of intersections is zero, the point is outside the polygon; if the number is one, the point is inside the polygon; if the number is two, the ray entered and exited the polygon, and therefore the point is outside the polygon; and if the number is three, the ray exited, entered, and re-exited the polygon, and therefore the point is inside the polygon. Since each ray entry in the polygon implies its exit, it becomes evident that if the number of intersections of the ray with the polygon is even then the point is outside the polygon, and if the number of intersections is odd then the point is inside the polygon.

The use of a ray to determine if a point is contained inside a polygon.

The implementation of this algorithm is relatively simple: given a point \(P=(P_x,P_y)\) and an array \(V_s\) with the vertices of the polygon \(V_0,V_1,\ldots{},V_n\), we need to create a ray with its origin in \(P\) and whose destination is any point \(Q=(Q_x,Q_y)\) far enough away from the polygon. To simplify, we will arbitrate that this ray is horizontal, which allows the destination point \(Q\) to have the same ordinate as \(P\), i.e., \(Q_y=P_y\). To ensure that \(Q\) is far enough away from the polygon, we can calculate the greatest abscissa of the vertices of the polygon \(\max({V_0}_x,{V_1}_x,\ldots,{V_n}_x)\) and add a small distance, for example, one unit.

The calculation of the greatest abscissa of the vertices is simple when using higher-order functions: we can start by mapping the array of vertices onto the array of its abscissas, and then operate a reduction of this second array with the max function. This reasoning can be translated into Julia through the following function:

Next, we have to determine the intersections between the segment defined by points \(P-Q\) and the segments that make up the polygon. Those segments are defined by the vertices \(V_0-V_1\), \(V_1-V_2\), ..., \(V_{n-1}-V_n\) and \(V_n-V_0\), which can be obtained through a mapping along two arrays, one with the points \(V_0,V_1,\ldots{},V_n\) and the other with the points \(V_1,\ldots{},V_n,V_0\). Thus, our function takes the following format:

The function that we map along the two arrays of vertices should produce, for each two consecutive vertices \(V_i\) and \(V_j\), the intersection of the segment \(P-Q\) with the segment \(V_i-V_j\).

To determine the intersection of two line segments, we can think as follows: given the points \(P_0\) and \(P_1\) that delimit a line segment, every point \(P_u\) between \(P_0\) and \(P_1\) is determined by the equation \[P_u = P_0 + u(P_1 - P_0), 0\leq u\leq 1\]

If the line segment \(P_0-P_1\) intersects another line segment delimited by the points \(P_2\) and \(P_3\) and defined by \[P_v = P_2 + v(P_3 - P_2), 0\leq v\leq 1\]

The intersection point between the segments \(P_0-P_1\) and \(P_2-P_3\) equals \(P_u\) and \(P_v\) such that: \[P_0 + u(P_1 - P_0) = P_2 + v(P_3 - P_2)\]

In the two-dimensional case, we have that \(P_i=(x_i,y_i)\), thus the previous equation unfolds in the two equations:

\[x_0 + u(x_1 - x_0) = x_2 + v(x_3 - x_2)\]

\[y_0 + u(y_1 - y_0) = y_2 + v(y_3 - y_2)\]

Solving the equation system, we obtain:

\[u=\frac{(x_3-x_2)(y_0-y_2)-(y_3-y_2)(x_0-x_2)}{(y_3-y_2)(x_1-x_0)-(x_3-x_2)(y_1-y_0)}\] \[v=\frac{(x_1-x_0)(y_0-y_2)-(y_1-y_0)(x_0-x_2)}{(y_3-y_2)(x_1-x_0)-(x_3-x_2)(y_1-y_0)}\]

Of course, for the divisions to be performed, the denominator \((y_3-y_2)(x_1-x_0)-(x_3-x_2)(y_1-y_0)\) has to be different from zero. If that is not the case, it is because the lines are parallel. On the other hand, the intersection can occur outside the line segments, thus we have to check if the obtained values ​​\(u\) and \(v\) obey the conditions \(0\leq u\leq 1\) and \(0\leq v\leq 1\). When this happens, the exact point of intersection can be calculated by replacing \(u\) in the equation \(P_u = P_0 + u(P_1 - P_0)\).

This leads us to the following definition for a function that calculates the intersection:

intersection_point(p0, p1, p2, p3) =

let denom = ((p3.y - p2.y)*(p1.x - p0.x) - (p3.x - p2.x)*(p1.y - p0.y))

denom == 0 ?

nothing :

let u = ((p3.x - p2.x)*(p0.y - p2.y) - (p3.y - p2.y)*(p0.x - p2.x))/denom

v = ((p1.x - p0.x)*(p0.y - p2.y) - (p1.y - p0.y)*(p0.x - p2.x))/denom

0 <= u <= 1 && 0 <= v <= 1 ?

xy(p0.x + u*(p1.x - p0.x), p0.y + u*(p1.y - p0.y)) :

nothing

end

end

Using this function, we can determine all the intersections between the segment \(P-Q\) and the edges of a polygon:

The mapping result will be, for each pair of consecutive vertices \(V_i-V_j\), the intersection point with the line segment \(P-Q\) or nothing when there is no intersection. Since we only want to know the intersections, we can now filter this array, keeping only the elements that are points, i.e., those that are not nothing:

Lastly, to know how many intersections occur, we just need to measure the length of the resulting array and check if the result is odd:

is_point_in_polygon(p, vs) =

let q = xy(reduce(max, map(cx, vs)) + 1, p.y)

rs = filter(e -> e != nothing,

map((vi, vj) -> intersection_point(p, q, vi, vj),

vs,

[vs[2:end]..., vs[1]]))

isodd(length(rs))

end

With the function is_point_in_polygon it is now possible to establish the association between each point and each polygon displayed in the plans provided by town hall services. However, as we mentioned initially, we need to be careful with the fact that it is also possible that a given polygon does not have an associated point or has more than one associated point. One way of treating all these situations identically will be to compute, for each polygon, the array of points it contains. Ideally, this array should have only one element, but if there are errors in a plan, it can have no elements or more than one element. In any case, the function always returns an array, so it will never produce an error.

To create this array of points for a given polygon, we must test each of the plan’s points to see if it belongs to the polygon. This can be obtained by filtering an array of points, keeping only those that are contained in the polygon. Thus, admitting that pts is the plan’s array of points and vs are the vertices of a particular polygon, we can define:

points_inside_polygon(pts, vs) =

filter(pt -> is_point_in_polygon(xyz(pt.x, pt.y, vs[1].z), vs),

pts)

Finally, given an array of points representing the coordinates of the top of the buildings and an array of polygons representing the building’s base perimeter (each implemented by the array of its vertices), we will create a three-dimensional representation of these buildings by determining, for each polygon, the number of points it contains and act accordingly:

• If that number is zero, we do not know what the height of the building is. This represents a plan error that we might not want to fix or that we want to fix by simply arbitrating a height. For now, we will choose to do nothing;

• If that number is one, then this point is located on top of the building and its \(z\) coordinate gives us the corresponding height;

• If that number is greater than one, then there are multiple applicable points, which can correspond to several possible situations:

  • If one of the points in question also belongs to another polygon, then it is because there are polygons inside other polygons or polygons that intersect each other. The first case may correspond to a building that varies in height (for example, having an engine room at the top). The second case is, typically, an error in the plan;

  • If the points all exclusively belong to the same polygon, then this can represent buildings with tops that are not horizontal planes.

The correct treatment of these latter cases is relatively complex. As we only want to create a prismatic approximation to the shape of the building, we will employ a simple approach: we will use the highest coordinate found as the building’s height. The following function implements these behaviors:

create_buildings(points, polygons) =

for polygon in polygons

let pts = points_inside_polygon(points, polygon),

n_pts = length(pts)

if n_pts == 0

nothing

elseif n_pts == 1

create_building(polygon, cz(pts[1]))

else

create_building(polygon, reduce(max, map(cz, pts)))

end

end

end

To create a building from the array of vertices of its perimeter and its height we will simply create a polygonal region at the base of the building that we then extrude to the top:

create_building(vs, z) =

z > 0 ?

extrusion(surface_polygon(vs), z) :

nothing

So far we solved the problem only from a geometrical point of view, representing the plan’s points by its coordinates and the polygons by the coordinates of its vertices. However, as we know, what the CAD tool provides are geometrical entities that contain the information we need. Now, given an array of entities, we can select those that correspond to points by filtering the entities that satisfy the predicate is_point. Then, from these entities, we can obtain their coordinates by applying a mapping process with the function point_position. This means that we can define the function that obtains the coordinates of all points that exist in an array of entities:

points(shapes) =

map(point_position,

filter(is_point,

shapes))

Likewise, given an array of entities, we can select the array of polygons’ vertices through a filtering process followed by a mapping process, as follows:

polygons(shapes) =

map(line_vertices,

filter(is_polygon,

shapes))

We should note that, in the previous definition, we are selecting both the closed polygonal lines and the opened ones. This is for dealing with the possibility of some polygonal lines being visually closed, even if they are not from the point of view of the geometrical operation that was used to create them.

Finally, we are in conditions of defining a function that, from the set of entities of a plan, creates the corresponding three-dimensional representation:

buildings_from_shapes(shapes) =

create_buildings(points(shapes), polygons(shapes))

Naturally, the set of entities to be processed can be selectively produced or, alternatively, be selected from a given plan. In the latter case, we only need to do:

buildings_from_shapes(all_shapes())

As an example, we show in this figure the result of evaluating the previous expression for the plan presented earlier in this figure.

Three-dimensional modeling of the buildings present in the plan shown in this figure.

9.8.2 Exercises 47

9.8.2.1 Question 175

The function is_point_in_polygon is not as efficient as it could be, since it performs a mapping followed by a filtering only to count the number of true elements in the final array. In practice, this combination of operations is no more than the counting of the number of times a binary predicate is satisfied along successive elements of two arrays. Define the operation how_many that implements this process. For example, consider:

> how_many(>, [1, 5, 3, 4, 6, 2], [1, 6, 2, 5, 4, 3])

2

9.8.2.2 Question 176

Re-implement the function is_point_in_polygon so that it uses the function how_many.

9.8.2.3 Question 177

Another possible approach to calculate the height of a building whose corresponding polygon contains multiple points is to use the average of the \(z\) coordinates of these points. Implement this approach.