I should give some justification for the complexity of the definition of "cross-section" given in the terminology section. Why, for instance, is the definition phrased in terms of D(P) and not dim(P)? Why is D(P) an important number? When I was first thinking about taking cross-sections of arbitrary polytopes, it seemed that the hierarchical nature of polytopes directly implied a recursive solution. That is, the cross-section of a polytope should be determined from the cross-sections of its faces. At this point I did not have the exact definition for "cross-section" that I have now, only some vague idea about "cutting" through objects, but I knew it had to do with intersections with R^k for some k. So the question was, for what polytopes P and integers k could the the intersection of P and R^k be determined by the intersection of R^k and the faces of P? Here's an example where this does not work:
There is an intersection of the square and R^1 (the point at the origin), but none of the square's (1-dimensional) faces intersected. On the other hand, the intersection with R^2 can be determined by the intersection of the faces:
There is an intersection between two of the edges and R^2, and drawing a line in between the face intersections gives the correct intersection of the square and R^2.
Another example of this, though not readily visualizable, is the square lying in the z-w plane bounded by the points (0, 0, 1, 1), (0, 0, 1, -1), (0, 0, -1, 1), and (0, 0, -1, -1). The intersection of this square with R^2 (the x-y plane) is a single point at the origin, but R^2 does not touch any of the edges of this square.
So exactly what characterizes a situation where the cross-section can or can not be computed recursively? Or in other words, when is the cross-section geometrically related to the cross-section of the faces, and when is it independant of the faces' cross-section? In the two cases above when the cross-section can not be computed recursively, the dimension of the space that the polytope lies in is more than one dimension greater than the dimension of subspace R^k we are intersecting with. The square in the y-z plane lies in R^3, and we are intersecting it with R^1. The square in the z-w plane lies in R^4, and we are intersecting it with R^2. But "dimension of the space that the polytope lies in" needs to be formalized a bit. Any point (x, y) in R^2 can be trivially identified with the "same" point (x, y, 0) in R^3. So, in some sense, a polygon that lies in the x-y plane could be said to also "lie" in R^3. This is why D(P) is a useful number: it expresses unambiguously the dimension of the space occupied by the polytope.
I don't mean to make it sound like I started my search for a cross-section algorithm with an a priori decision that it be recursive. It's just that I had no idea how to efficiently compute cross-sections if it wasn't going to be based on some sort of recursion. How, for instance, would one find the intersection of the line and square shown above, where the line passes through the square but through none of the edges? It could be done with some linear algebra. One could determine the interection of the 2-flat (spanned by the square) with the 1-flat (the x axis) and then determine if that intersection lay inside the square. In general, one would have to determine the intersection of the n-flat (spanned by some polytope) with a k-flat, and then determine how much of that intersection lay inside the polytope. But I don't know how to do this- given a point location and a polytope, I don't (yet) know how to determine if the point is inside or outside the polytope. Also, I think this sort of solution would rely on a clipping algorithm in n dimensions, something else I know nothing about. On the other hand, a recursive solution seemed intuitively more elegant, simply because there would be far less number crunching. The only time coordinate information would be relevent was at the base case- computing the cross-section of the polytope's edges or vertices. After that, the only necesary information would be the polytope's face lattice.
So, can I prove that if k = D(P) - 1, and if P intersects R^k, then R^k intersects some of the faces of P? The following proof has a few holes which I need to patch but it should be somewhat convincing. Let us restrict our consideration to R^(k + 1), because everything is happening in this space and not in any lower dimensional one. R^k splits R^(k+1) into two open half-spaces: those points with their last coordinate greater than 0, and those with their last coordinate less than 0. Suppose P intersects with R^k, and pick some point p in this intersection. If p is on the boundary of P, then we are done, because being on the boundary of P means being on a face of P, and since p is in R^k, then we have an intersection of one of the faces and R^k. So suppose p is in the interior of P, which means that some of the boundary lies below R^k, and some lies above R^k:
And suppose, for the sake of contradiction, that none of the boundary intersects R^k. If there is no intersection of R^k and the boundary of P, then we have formed a separation of the boundary into two parts: the intersection of the boundary with the open half-space below R^k, and the intersection of the boundary with the open half-space above R^k. But this is impossible, because the boundary of a polytope is connected, and can have no separation. Therefore there must be an intersection of R^k and the boundary of P.
Furthermore, if k is any less than D(P) - 1, then we do not have the implication: if R^k intersects P then R^k intersects boundary of P. This was demonstrated in the example of the square lying in the y-z plane being intersected with R^1 (the x axis). And if k was any bigger, then R^k would contain the polytope, and intersecting it with R^k would recursively produce a copy of the original polytope. So there is exactly one value of k such that R^k "cuts" the polytope in some sense and such that the nature of this cut can be determined recursively: k = D(P) - 1. This is why I defined cross-section to be the intersection of P with a (D(P) - 1)-flat.
It would seem that D should be stored in the object representation, so that when we wanted to do a cross-section, we would know exactly what value of k to use so that the intersection of the object with R^k is a cross-section of the object. But there is another quantity which is a little more relevant to the purpose of Peek. Suppose the user is moving a standard 3-d cube around in R^4, and watching the resulting intersections with R^3 change. This is an interesting thing to do since such intersections appear as different orientations (in R^3) of the familiar 2-d cross-sections of a cube. If the cube C is being rotated in R^4, then D(C) is in general equal to 4, so the intersections with R^3 are in general cross-sections. But of course, the user might happen to orient the cube so that it lies perfectly in R^3, so D(C) would be 3, and the intersection with R^3 would not be a cross-section. But at this point the user wouldn't want the true cross-section (the intersection with R^2) because this result would be inconsistent with all previous ones. So even though D(C) = 3, we still want to show the user the cube's intersection with R^3, since for continuity's sake, the cube should be considered to occupy R^4. So instead of storing D in the object representation, we store another quantity that represents the user's current idea about what dimension the object occupies. The name for this psychological version of D is space_d, after the variable name space_d that appears in the C code. Nearly all of the time D and space_d agree, but for the user, space_d is the relevant value. The algorithms below are called "cross_edge" and "cross_polytope" and so on, as if to suggest that they are always taking cross-sections, and most of the time, they are, but they have been designed so that if they are intersecting the object with R^k in order to "cut" the object, and the object already lies inside R^k, then the result they return is simply a copy of the original object.
Since space_d pertains to the whole object, it is stored in the topmost level of the object (the piece level) in the first piece, as described in the the section on enhancements to the data structure. For some object o which has p as its first piece, and one of p's surfaces is a polytope P, I'll refer to space_d for this object as either space_d(o), space_d(p), or space_d(P), whichever is most appropriate for the context, since they are all the same.
cross_vert(v, k) {
if coord(v)[k] < 0 {
spare(v) <-- -1;
result(v) <-- NULL;
}
else if coord(v)[k] > 0 {
spare(v) <-- 1;
result(v) <-- NULL;
}
else {
spare(v) <-- 0;
create new part r;
copy all info from v into r;
result(v) <-- r;
}
done(v) <-- true;
}
The effect of this function is that vertices remember if they
intersected with R^k- if they did, then their result pointer
points to their duplicate, if not, then spare tells how the
point failed to intersect. spare is just that- a spare
variable that parts have to hold extra information as the need arises.
In this case, it is saying whether the vertex represented by v
lay above (spare = 1), below (spare = -1), or inside (spare = 0) R^k.
This was determined by looking at coord(v)[k]. Remember
that coord(v) is a zero-based coordinate array, so
the coord(v)[k] is the coordinate in the dimension that
is being eliminated in this cross-section. For instance, if
space_d = 3, then <coord(v)[2] is the Z
coordinate. In the cases that a vertex does lie in R^k, then it is
necessary to create a new part, and set result to this, as
opposed to setting result(v) <-- v, because the
cross-section has got to be a new object, independent of the the
original. The original may be disposed of as soon as this
cross-section is done.
Now we create algorithm for the cross-section of an edge. Let E be an edge on the polytope P, let e be the part in the data structure representing E. k is as before.
cross_edge(e, k) {
v1 <-- part representing first endpoint of edge E
v2 <-- part representing second endpoint of edge E
if not done(v1)
cross_vert(v1, k);
if not done(v2)
cross-vert(v2, k);
if (spare(v1)) * (spare(v2)) > 0 {
(both endpoints lay on same side of R^k, there is no intersection
of the edge)
result(e) <-- NULL;
}
else if (spare(v1)) * (spare(v2)) = 0 {
(one or both endpoints lay inside R^k)
if spare(v1) = 0 and spare(v2) = 0 {
(both endpoints lie inside R^k, hence whole edge lies in R^k)
create new part r;
copy all info from e into r;
create image list below r to point to result(v1)
and result(v2);
result(e) <-- r;
}
else if spare(v1) = 0 {
(only first endpoint lies in R^k, edge's intersection is
that endpoint)
result(e) <-- result(v1);
}
else {
(only second endpoint lies in R^k, edge's intersection is
that endpoint)
result(e) <-- result(v2);
}
}
else {
(endpoints lie on different sides of R^k, intersection is
some interior point of edge)
create new part r;
copy color info from e to r
calculate location of intersection of edge and R^k,
put this info in coord(r);
dim(r) <-- 0;
result(e) <-- r;
}
done(e) <-- true;
}
Note how the product of endpoint's spare values was used to determine the nature of the intersection of the edge.
The base case of calculating the cross-section for the edge also illustrates something that will happen again in the general case discussed below. Sometimes the result of an operation on a polytope is simply the result of the operation on one of the faces. In these cases, there is not a new part made to represent the result for the polytope, rather, result is set to the result for the face. For the edges undergoing cross-section, this happens when only one of the endpoints intersects R^k. Then, result for the edge is set to result for that endpoint and no new part is created. Otherwise, a new part is created to represent either the whole edge or an interior point of the edge which intersected R^k.
Also note that the base case is the only time we deal with coordinate information. From now on, the only necessary information (besides the results of cross-sections on lower dimensional faces) is the pattern of links set up by the image lists.
Based on consideration of cases like these, and more complicated higher dimensional cases, I've come up with what I think is a correct solution to this problem. In the following, "intxion" is a boolean flag that is true if any of the faces reported intersections, and the array "report" holds information about how many faces reported intersections of which dimension: report[i] is the number of faces that reported intersections of dimension i. Let n = dim(P). So if a face reports a result of dimension n-1 then we know the result of the cross-section on that face was a copy of the whole face. Rather than give a full account of how the following algorithm came to be, I'll just label all the possible cases and the provide illustrations, hoping that those sufficiently motivate the complexity of the algorithm:
cross_polytope(p, k) {
for each face f of p {
if not done(f)
if dim(f) > 1
cross_polytope(f, k)
else
cross_edge(f, k)
if result(f) != NULL {
intxion <-- true;
increment report[dim(result(f))];
}
}
if not intxion {
(Case A)
(no face had a result, so polytope has no result)
result(p) <-- NULL;
}
else if report[n-1] = number of faces {
(Case B)
(every face is contained by R^k, hence polytope is contained
by R^k)
create new part r;
copy info from p into r;
create image list under r pointing to each result(f)
result(p) <-- r;
}
else if report[n-1] = 1 {
(Case C)
(one face, call it f1, is contained by R^k, report only
that face as the cross-section of the whole polytope P)
result(p) <-- result(f1)
}
(***if report[n-1] is nonzero but is not 1 or the number
of faces, then there is something wrong with the polytope***)
else if report[n-2] > 0 {
(gather unique (n-2) dimensional results under an image
list and count the unique results in num_uniq)
for each face f such that dim(result(f)) = n-2 {
if result(f) has not already been added to image
list of results {
extend image list to point to result(f);
increment num_uniq;
}
}
if num_uniq > n {
(Case D)
(the unique (n-2) dimensional results constitute the boundary
of the (n-1) dimensional cross-section for this polytope)
create new part r;
color(r) <-- color(p);
dim(r) <-- n-1;
make r point down to image list of unique results
created above;
result(p) <-- r;
}
else if num_uniq = 1 {
(Case E)
(one face, call it f1, was among those that reported the
single (n-2) dimensional result, hence that is the result
for the polytope);
result(p) <-- result(f1);
}
(***if num_uniq is between 1 and n then something is wrong
with the polytope***)
}
else {
(Case F)
(we descend through remaining dimensions until we find the
highest dimension face cross-section result, and then pick an
arbitrary result of that dimension because there should only
be one unique result)
d <-- n-3;
while report[d] = 0
decrement d;
result(p) <-- result(any face with d-dimensional result)
}
done(p) <-- true;
}
Case A: no intersection at all:
Case B: report[n-1] = number of faces: (n = 2)
Case C: report[n-1] = 1: (n = 2)
If more than 1 but not all the faces lay in R^k then there must be some non-coplanarity in the polytope- the coordinate information must be corrupted as in the following example where n = 2 and report[n-1] = 2:
Case D: number of unique (n-2) dimensional results > n: (n = 3)
This is the sort of situation where the naive idea of simply stringing together all the face's results would have worked.
Case E: number of unique (n-2) dimensional results = 1: (n = 3)
If the number of unique (n-2) dimensional results is between 1 and n, then there is a non-coplanarity in one of the faces of the polytope, as in this example where n = 3 and report[n-1] = 2:
Case F: single unique (n-3) or lower dimensional result: (n = 3)
Now, we still haven't processed the top two layers of the data structure: the piece layer and the surface layer. These are simple to deal with. Each surface is a polytope, so call cross_polytope() on each surface. For each of the surfaces that have an intersection, create an image list to point to each one, and be sure to maintain the order of the surfaces when creating the image list. If an object has two surfaces, an outer and an inner, then the outer surface is listed first. So if both surfaces report intersections, then the intersection of the outer surface is the outer surface of the intersection, and likewise for the inner surfaces, so the intersection of the outer surface should be the first surface listed. For the piece layer, simply string together any piece results that are non-NULL; the order doesn't matter.
Note in the pseudo-code below that for an object o, the value space_d(o) - 1 is passed as the second parameter in all calls to cross_polytope(). This means that cross_polytope() creates polytopes that do not occupy dimension space_d(o), which in turn implies that the new object o2 can't occupy dimension space_d(o) either. That is, D(o2) must be less than or equal to space_d(o) - 1. Thus we are justified in setting space_d(o2) to space_d(o) - 1. It also further justifies the naming of these functions "cross_object", "cross_polytope", and so on: even if applying cross_object(o) to an object o does not actually take a cross-section of o and thereby decremented the value of D(o), we have decremented the value for space_d(o).
cross_object(o) {
for each piece p in object o {
for each surface s in piece p {
cross_polytope(s, space_d(o) - 1);
}
if any of the surfaces had an intersection {
create an image list to point to each result(s)
which respects the order the surfaces are listed in;
create new piece to point down to this image list;
}
create new object o2 by linking together any new pieces
just created;
space_d(o2) <-- space_d(o) - 1;
return(o2);
}