[Numpy-discussion] Iterative Matrix Multiplication

[Friedrich Romstedt]

I'm learning too by answering your questions.

2010/3/2 Ian Mallett <geometrian at gmail.com>:
> 1  v_array #array of beautifully transformed vertices from last step
> 2  n_array #array of beautifully transformed normals from last step
> 3  f_list  #list of vec4s, where every vec4 is three vertex indices and a
> 4          #normal index.  These indices together describe a triangle--
> 5          #three vertex indices into "v_array", and one normal from
> "n_array".
> 6  edges = []
> 7  for v1index,v2index,v3index,nindex in f_list:
> 8      v1,v2,v3 = [v_array[i] for i in [vi1index,v2index,v3index]]
> 9      if angle_between(n_array[nindex],v1-a_constant_point)<90deg:
> 10         for edge in [[v1index,v2index],[v2index,v3index],[v3index,v1index]]:
> 11             #add "edge" to "edges"
> 12 #remove pairs of identical edges (also note that things like
> 13 #[831,326] are the same as [326,831])
> 14 edges2 = []
> 15 for edge in edges:
> 16     edges2.append(v_array[edge[0]],v_array[edge[1]])
> 17 return edges2

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]])

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]])

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?

> I may be able to do something with lines 14-17 myself--but I'm not sure.
Ok, than that's fine, and let us all know about your solution.

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