I'm learning too by answering your questions.
2010/3/2 Ian Mallett <geometr...@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
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
