# [Numpy-discussion] Iterative Matrix Multiplication

Ian Mallett geometrian@gmail....
Fri Mar 5 00:13:31 CST 2010

```Firstly, I want to thank you for all the time and attention you've obviously
put into this code.

On Tue, Mar 2, 2010 at 12:27 AM, Friedrich Romstedt <
friedrichromstedt@gmail.com> wrote:

> The loop I can replace by numpy operations:
>
> >>> v_array
> array([[1, 2, 3],
>       [4, 5, 6],
>       [7, 8, 9]])
> >>> n_array
> array([[ 0.1,  0.2,  0.3],
>       [ 0.4,  0.5,  0.6]])
> >>> f_list
> array([[0, 1, 2, 0],
>       [2, 1, 0, 1]])
>
> Retrieving the v1 vectors:
>
> >>> v1s = v_array[f_list[:, 0]]
> >>> v1s
> array([[1, 2, 3],
>       [7, 8, 9]])
>
> Retrieving the normal vectors:
>
> >>> ns = n_array[f_list[:, 3]]
> >>> ns
> array([[ 0.1,  0.2,  0.3],
>       [ 0.4,  0.5,  0.6]])
>
I don't think you understand quite what I'm looking for here.  Every vec4 in
f_list describes a triangle.  The first, second, and third are indices of
vertices in v_array.  The fourth is an index of n_array.  From your code,
I've learned that I can do this, which is more what I want:
v1s = v_array[f_list[:,0:3]]
ns = n_array[f_list[:,3]]
Obviously, this changes the arrays quite a bit.  With changes, the code
doesn't actually crash until the corners = line.  Do the following lines
that do the dot product and comparison still behave correctly?  They should
be finding, for each triangle, whether or not the associated normal is less
than 90 degrees.  The triangle's edges should then be added to an array.
See end.

> Now how to calculate the pairwise dot product (I suppress the
> difference of v1 to some_point for now):
>
> >>> inner = numpy.inner(ns, v1s)
> >>> inner
> array([[  1.4,   5. ],
>       [  3.2,  12.2]])
>
> This calculates *all* pairwise dot products, we have to select the
> diagonal of this square ndarray:
>
> >>> dotprods = inner[[numpy.arange(0, 2), numpy.arange(0, 2)]]
> >>> dotprods
> array([  1.4,  12.2])
>
> Now we can create a boolean array saying where the dotprod is > 0
> (i.e, angle < 90°), and select those triangles:
>
> >>> select = dotprods > 0
> >>> select
> array([ True,  True], dtype=bool)
> >>> selected = f_list[select]
> >>> selected
> array([[0, 1, 2, 0],
>       [2, 1, 0, 1]])
>
This seems like a clever idea.

> In this case it's the full list.  Now build the triangles corners array:
>
> >>> corners = v_array[selected[:, :3]]
> >>> corners
> array([[[1, 2, 3],
>        [4, 5, 6],
>        [7, 8, 9]],
>
>       [[7, 8, 9],
>        [4, 5, 6],
>        [1, 2, 3]]])
> >>>
>
> This has indices [triangle, vertex number (0, 1, or 2), xyz].
> And compute the edges (I think you can make use of them):
>
> >>> edges_dist = numpy.asarray([corners[:, 1] - corners[:, 0], corners[:,
> 2] - corners[:, 0], corners[:, 2] - corners[:, 1]])
> >>> edges_dist
> array([[[ 3,  3,  3],
>        [-3, -3, -3]],
>
>       [[ 6,  6,  6],
>        [-6, -6, -6]],
>
>       [[ 3,  3,  3],
>        [-3, -3, -3]]])
>
> This has indices [corner number, triangle, xyz].
> I think it's easier to compare then "reversed" edges, because then
> edge[i, j] == -edge[k, l]?
>
> But of course:
>
> >>> edges = numpy.asarray([[corners[:, 0], corners[:, 1]], [corners[:, 1],
> corners[:, 2]], [corners[:, 2], corners[:, 0]]])
> >>> edges
> array([[[[1, 2, 3],
>         [7, 8, 9]],
>
>        [[4, 5, 6],
>         [4, 5, 6]]],
>
>
>       [[[4, 5, 6],
>          [4, 5, 6]],
>
>        [[7, 8, 9],
>          [1, 2, 3]]],
>
>
>       [[[7, 8, 9],
>         [1, 2, 3]],
>
>        [[1, 2, 3],
>         [7, 8, 9]]]])
>
> This has indices [edge number (0, 1, or 2), corner number in edge (0
> or 1), triangle].
> But this may not be what you want (not flattened in triangle number).
> Therefore:
>
> >>> edges0 = numpy.asarray([corners[:, 0], corners[:, 1]])
> >>> edges1 = numpy.asarray([corners[:, 1], corners[:, 2]])
> >>> edges2 = numpy.asarray([corners[:, 2], corners[:, 0]])
> >>> edges0
> array([[[1, 2, 3],
>        [7, 8, 9]],
>
>       [[4, 5, 6],
>        [4, 5, 6]]])
> >>> edges1
> array([[[4, 5, 6],
>         [4, 5, 6]],
>
>       [[7, 8, 9],
>         [1, 2, 3]]])
> >>> edges2
> array([[[7, 8, 9],
>        [1, 2, 3]],
>
>       [[1, 2, 3],
>        [7, 8, 9]]])
>
> >>> edges = numpy.concatenate((edges0, edges1, edges2), axis = 0)
> >>> edges
> array([[[1, 2, 3],
>        [7, 8, 9]],
>
>       [[4, 5, 6],
>        [4, 5, 6]],
>
>       [[4, 5, 6],
>         [4, 5, 6]],
>
>       [[7, 8, 9],
>         [1, 2, 3]],
>
>       [[7, 8, 9],
>        [1, 2, 3]],
>
>        [[1, 2, 3],
>        [7, 8, 9]]])
>
> This should be as intended.
> The indices are [flat edge number, edge side (left or right), xyz].
>
> Now I guess you have to iterate over all pairs of them, don't know a
> numpy accelerated method.  Maybe it's even faster to draw the edges
> twice than to invest O(N_edges ** 2) complexity for comparing?
>
Unfortunately, no.   The whole point of the algorithm is to extrude
back-facing triangles (those with normals facing away from the light)
backward, leaving polygons behind in a light volume.  Although a shadow
volume tutorial can explain this in a more detailed way, basically, for
every triangle, if it is back-facing, add its edges to a list.  Remove the
duplicate edges.  So, the edge between two back-facing
triangles-that-meet-along-an-edge is not kept.  However, if a back-facing
triangle and a front-facing triangle meet, only the back-facing triangle's
edge is present in the list, and so it is not removed.  Thus, only the
border edges between the front-facing triangles and the back facing
triangles remain in the list.  As front-facing triangles face the light and
back-facing triangles don't, a silhouette edge is built up.  When these
edges are extruded, they're extremely useful.  The following image
http://www.gamedev.net/reference/articles/1873/image010.gif shows the
process.  The four triangles are all back facing.  If duplicate edges are
removed, only edges I0-I2, I2-I4, I4-I3, and I3-I0 remain--the silhouette
edges.

Still to do, remove the duplicate edges (actually where a good deal of the
optimization lies too).

So, for every back-facing triangle [v1,v2,v3,n], (where v*n* is a vertex *
index*), the edges [v1,v2], [v2,v3], and [v3,v1] need to be added to a
list.  I.e., f_list needs to be converted into a list of edges in this way.
Then, duplicate edge pairs need to be removed, noting that [v1,v2] and
[v2,v1] are still a pair (in my Python code, I simply sorted the edges
before removing duplicates: [123,4] -> [4,123] and [56,968] -> [56,968]).
The final edge list then should be converted back into actual vertices by
indexing it into v_array (I think I understand how to do this now!): [
[1,2], [4,6], ... ] -> [ [[0.3,1.6,4.5],[9.1,4.7,7.7]],
[[0.4,5.5,8.3],[9.6,8.1,0.3]], ... ]

> It may seem a bit complicated, but I hope this impression is mainly
> because of the many outputs ...
>
> I double-checked everything, *hope* everything is correct.

> So far from me,
> Friedrich
>
Once again, thanks so much,
Ian
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20100304/28bb206e/attachment.html
```