All right, first thing first. The main part of the cure is to construct the admissibility table that, for each pair of directions, lists all directions that may form a vertex with the given two. Surprisingly, it is very short. I'll keep the directions in the 0,$\pm 1$ format like you (not normalized). I assume that you keep an array of these directions p of length 26. The only thing I want you to make sure about the order is that p[2k] and p[2k+1] are opposite for each k=0,1,...,12. This will spare a few operations in the edge tracing algorithm.
First, let us exclude all degenerate pairs. That is pretty easy. The pair (A,B) is degenerate iff A+B is a multiple of some other direction in the list. In practical terms, it means that all non-zero entries of $A+B$ are equal in absolute value. The first thing
is to create the table of pairs of directions and precompute some information for them.
To make it working for both algorithms, I suggest that for each pair (i,j) of indices, you create the structure Q[i,j] with the following fields
Boolean nondegeneracy; true if the pair is non-degenerate, false if it is degenerate.
You should get 288 non-degenerate pairs if the order of i,j matters (144 with i<j)
vector normal; normal=$p[i]\times p[j]$ is the direction orthogonal to both p[i] and p[j]
iprojector, jprojector; those are just found as $v=p[i]\times\operatorname{normal}$; iprojector$=v/(v\cdot p[i])$ and similarly for jprojector. This will allow you to find a point X with given scalar products $a=p[i]\cdot X$ and $b=p[j]\cdot X$ as a times iprojector plus b times jprojector (6 multiplications and 3 additions)
integer array goodindices. It will contain all indices k such that p[i],p[j],p[k] may make a vertex. To construct it, you take your pair a=p[i],b=p[j] of directions and include k into goodindices if the following two conditions hold:
1)the determinant det=det(a,b,c) is not 0
- for every index m different from i,j,k at least one of the determinants det(a,b,d), det(a,d,c), and det(d,b,c) where d=p[m] has sign opposite to the sine of det, i.e., if, say, det > 0, one of them is strictly less than 0 (since all coordinates are integer, compare to -0.5). This checks that there is no direction in the cone formed by a,b,c.
Funnily enough, the total number of indices in goodindices is never greater than 6.
Also, let Max[i] be the maximal scalar product of your object points with p[i]
(I prefer to keep all 26 directions and Max array to keeping 13 and Min and Max arrays just to avoid any casework anywhere)
Now, when the table is ready, it is very easy to run the N^3 algorithm. For each pair i,j with i < j do the following.
If Q[i,j].nondegeneracy=false, just skip the pair altogether
Otherwise put X=M[i]Q[i,j].iprojector+M[j]Q[i,j].jprojector.
Determine tmin as the minimum over all indices k in Q[i,j].goodindices such that a=p[k]$\cdot$Q[i,j].normal < 0 of the ratios (M[k]-X$\cdot$p[k])/a.
Determine tmax as the maximum over all indices k in Q[i,j].goodindices such that a=p[k]$\cdot$Q[i,j].normal > 0 of the ratios (M[k]-X$\cdot$p[k])/a.
You should optimize here updating tmin and tmax when looking at each k and stopping immediately if tmax gets less than or equal to tmin. This cycle over k is the most time-consuming part when defining the edge, so how you write it determines the real speed of this algorithm. Optimize as much as possible here to make sure that you do not do any extraction from array twice, etc. Each little extra operation counts!
If tmin < tmax, compute the points X+tmin Q[i,j].normal and X+tmax Q[i,j].normal and draw the edge between them (or store it for later drawing).
That's all. It should run fast enough to be acceptable already and it is robust in the sense that you do not need to perturb your maxes, check for paralellness anywhere, etc., etc.
The other algorithm is still about 4 times faster but I should check that it has no stability issues before posting it. Try this first and tell me the running time.
