[Numpy-discussion] Iterative Matrix Multiplication
Hi, I have a list of vec3 lists (e.g. [[1,2,3],[4,5,6],[7,8,9],...]). To every single one of the vec3 sublists, I am currently applying transformations. I need to optimize this with numpy. To get proper results, as far as I can tell, the vec3 lists must be expressed as vec4s: [[1,2,3],[4,5,6],[7,8,9],...] -> [[1,2,3,1],[4,5,6,1],[7,8,9,1],...]. Each of these needs to be multiplied by either a translation matrix, or a rotation and translation matrix. I don't really know how to approach the problem . . . Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
On Sun, Feb 28, 2010 at 6:31 PM, Charles R Harris wrote: > As I understand it, you want *different* matrices applied to each vector? Nope--I need the same matrix applied to each vector. Because 3D translation matrices must, if I understand correctly be 4x4, the vectors must first be changed to length 4 (adding a 1 for the last term). Then, the matrices would be applied. Then, the resulting n*4 array would be converted back into a n*3 array (simply by dropping the last term). Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
On Sun, Feb 28, 2010 at 6:54 PM, Charles R Harris wrote: > Why not just add a vector to get translation? There is no need to go the > homogeneous form. Or you can just leave the vectors at length 4 and use a > slice to access the first three components. That way you can leave the ones > in place. > Oh . . . duh :D Excellent--and a 3D rotation matrix is 3x3--so the list can remain n*3. Now the question is how to apply a rotation matrix to the array of vec3? > Chuck > Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
Excellent--this setup works perfectly! In the areas I was concentrating on, the the speed increased an order of magnitude. However, the overall speed seems to have dropped. I believe this may be because the heavy indexing that follows on the result is slower in numpy. Is this a correct analysis? What would be amazing is if I could port the other half of the algorithm to numpy too . . . The algorithm is for light volume extrusion of geometry. Pseudocode of the rest of the algorithm: 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 If that needs clarification, let me know. The goal with this is obviously to remove as many looping operations as possible. I think the slowdown may be in lines 8, 9, and 16, as these are the lines that index into v_array or n_array. In line 9, the normal "n_array[nindex]" is tested against any vector from a vertex (in this case "v1") to a constant point (here, the light's position) see if it is less than 90deg--that is, that the triangle's front is towards the light. I thought this might be a particularly good candidate for a boolean array? The edge pair removal is something that I have worked fairly extensively on. By testing and removing pairs as they are added (line 11), a bit more performance can be gained. I have put it here in 12-13 because I'm not sure this can be done in numpy. My current algorithm works using Python's standard sets, and using the "edge" as a key. If an edge is already in the set of edges (in the same or reversed form), then the edge is removed from the set. I may be able to do something with lines 14-17 myself--but I'm not sure. If my code above can't be simplified using numpy, is there a way to efficiently convert numpy arrays back to standard Python lists? As mentioned before, I'm certainly no pro with numpy--but I am learning! Thanks! Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
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 < friedrichromst...@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 i
Re: [Numpy-discussion] Iterative Matrix Multiplication
Cool--this works perfectly now :-) Unfortunately, it's actually slower :P Most of the slowest part is in the removing doubles section. Some of the costliest calls: #takes 0.04 seconds inner = np.inner(ns, v1s - some_point) #0.0840001106262 sum_1 = sum.reshape((len(sum), 1)).repeat(len(sum), axis = 1) #0.032923706 sum_2 = sum.reshape((1, len(sum))).repeat(len(sum), axis = 0) #0.026504089 comparison_sum = (sum_1 == sum_2) #0.0909998416901 diff_1 = diff.reshape((len(diff), 1)).repeat(len(diff), axis = 1) #0.0340001583099 diff_2 = diff.reshape((1, len(diff))).repeat(len(diff), axis = 0) #0.026504089 comparison_diff = (diff_1 == diff_2) #0.023019073 same_edges = comparison_sum * comparison_diff #0.12848502 doublet_count = same_edges.sum(axis = 0) Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
Hi, On Sat, Mar 6, 2010 at 1:20 AM, Friedrich Romstedt < friedrichromst...@gmail.com> wrote: > Hmm. Let's see ... Can you tell me how I can test the time calls in a > script take? I have no idea. > I've been doing: t1 = time.time() #code t2 = time.time() print "[code]:",t2-t1 > > #0.0840001106262 > > sum_1 = sum.reshape((len(sum), 1)).repeat(len(sum), axis = 1) > > > > #0.032923706 > > sum_2 = sum.reshape((1, len(sum))).repeat(len(sum), axis = 0) > > > > #0.026504089 > > comparison_sum = (sum_1 == sum_2) > > We can leave out the repeat() calls and leave only the reshape() calls > there. Numpy will substitute dimi == 1 dimensions with stride == 0, > i.e., it will effectively repeat those dimension, just as we did it > explicitly. > Wow! Drops to immeasurably quick :D > > #0.0909998416901 > > diff_1 = diff.reshape((len(diff), 1)).repeat(len(diff), axis = 1) > > > > #0.0340001583099 > > diff_2 = diff.reshape((1, len(diff))).repeat(len(diff), axis = 0) > > > > #0.026504089 > > comparison_diff = (diff_1 == diff_2) > > Same here. Delete the repeat() calls, but not the reshape() calls. > Once again, drops to immeasurably quick :D > > #0.023019073 > > same_edges = comparison_sum * comparison_diff > > Hmm, maybe use numpy.logical_and(comparison_sum, comparison_diff)? I > don't know, but I guess it is in some way optimised for such things. > It's marginally faster. > > #0.12848502 > > doublet_count = same_edges.sum(axis = 0) > > Maybe try axis = 1 instead. I wonder why this is so slow. Or maybe > it's because he does the conversion to ints on-the-fly, so maybe try > same_edges.astype(numpy.int8).sum(axis = 0). > Actually, it's marginally slower :S > Hope this gives some improvement. I attach the modified version. > > Ah, one thing to mention, have you not accidentally timed also the > printout functions? They should be pretty slow. > Nope--I've been timing as above. > Friedrich > On Sat, Mar 6, 2010 at 1:42 AM, Friedrich Romstedt < friedrichromst...@gmail.com> wrote: > 2010/3/5 Ian Mallett : > > #takes 0.04 seconds > > inner = np.inner(ns, v1s - some_point) > > Ok, I don't know why I was able to overlook this: > > dotprod = (ns * (v1s - some_point)).sum(axis = 1) > Much faster :D So, the most costly lines: comparison_sum = (sum_1 == sum_2) #0.024 sec comparison_diff = (diff_1 == diff_2) #0.024 sec same_edges = np.logical_and(comparison_sum, comparison_diff) #0.029 sec doublet_count = same_edges.sum(axis = 0) #0.147 Thanks again, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Iterative Matrix Multiplication
On Sat, Mar 6, 2010 at 12:03 PM, Friedrich Romstedt <
friedrichromst...@gmail.com> wrote:
> At the moment, I can do nothing about that. Seems that we have
> reached the limit. Anyhow, is it now faster than your Python list
> implementation, and if yes, how much? How large was your gain by
> using numpy means at all? I'm just curious.
>
Unfortunately, the pure Python implementation is actually an order of
magnitude faster. The fastest solution right now is to use numpy for the
transformations, then convert it back into a list (.tolist()) and use Python
for the rest.
Here's the actual Python code.
def glLibInternal_edges(object,lightpos):
edge_set = set([])
edges = {}
for sublist in xrange(object.number_of_lists): #There's only one sublist
here
face_data = object.light_volume_face_data[sublist]
for indices in face_data: #v1,v2,v3,n
normal = object.transformed_normals[sublist][indices[3]]
v1,v2,v3 = [ object.transformed_vertices[sublist][indices[i]]
for i in xrange(3) ]
if
abs_angle_between_rad(normal,vec_subt(v1,lightpos))
Thanks,
Ian
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Why is the shape of a singleton array the empty tuple?
>>> x = numpy.array(3) >>> x array(3) >>> x.shape () >>> y = numpy.array([3]) >>> y array([3]) >>> y.shape (1,) Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Why is the shape of a singleton array the empty tuple?
On Sat, Mar 6, 2010 at 9:46 PM, David Goldsmith wrote: > Thanks, Ian. I already figured out how to make it not so, but I still want > to understand the design reasoning behind it being so in the first place > (thus the use of the question "why (is it so)," not "how (to make it > different)"). > Well, I can't help you with that. I would also ask why this design even exists? Equating an array with a single number doesn't make sense to me. Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] matrix operation
>>> import numpy >>> A = numpy.array([[2,3],[10,12]]) >>> B = numpy.array([[1,4],[9,13]]) >>> C = numpy.array([A,B]) >>> numpy.min(C,0) array([[ 1, 3], [ 9, 12]]) Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] get information/position of value from matrix
On Wed, Apr 7, 2010 at 7:40 AM, ioannis syntychakis wrote: > Hallo Everybody, > > I am new in this mail list and python. But I am working at something and I > need your help. > > I have a very big matrix. What I want is to search in that matrix for > values above the (for example:) 150. If there are values above the 150, I > also want to get their position and value. > I think: your_matrix = # indices = np.where(matrix>150.0) values = your_matrix[indices] should work. > Is this possible? > Thanks in advance!! > > ps. could you also tell me how i can show the whole matrix? because now i > see: > > [[ 0. 0. 0. ..., 87. 4. 4.] > [ 0. 0. 0. ..., 89. 6. 4.] > [ 0. 1. 0. ..., 87. 4. 3.] > ..., > [ 0. 0. 0. ..., 87. 4. 4.] > [ 0. 3. 0. ..., 87. 4. 4.] > [ 0. 0. 4. ..., 87. 4. 4.]] > set_printoptions(precision=None, threshold=None, edgeitems=None, linewidth=None, suppress=None, nanstr=None, infstr=None) Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Physics Engine
Hi, I'm working on various projects, I think it would be great to have a built-in physics engine to my rendering library. I would rather not rely on a large dedicated Physics library (PyOGRE, PyODE, etc.). I thought it might be interesting to try to make a simple soft-body physics engine with NumPy. To this end, over the past few hours, I've constructed this per-vertex physics object class, which plugs directly into my rendering object's class. So far, it implements global forces, collision detection with a horizontal plane, and length tensors: class glLibCPUSoftPhysMesh: def __init__(self,object): self.object = object self.original_vertices = np.array(object.raw_vertices) #object.raw_vertices is the unique list of vertices self.vertices = np.array(object.raw_vertices) self.vertex_indices = np.array(object.indexed_vertices[0]) #object.indexed_vertices[0] is a list of indices into object.raw_vertices self.vertex_speeds = np.zeros((self.vertices.shape[0],3)) self.delta_pos = np.zeros((self.vertices.shape[0],3)) v1s = self.original_vertices[self.vertex_indices[0::3]] v2s = self.original_vertices[self.vertex_indices[1::3]] v3s = self.original_vertices[self.vertex_indices[2::3]] self.target_lengths1 = np.sum((v1s-v2s)**2.0,axis=1)**0.5 self.target_lengths2 = np.sum((v2s-v3s)**2.0,axis=1)**0.5 self.target_lengths3 = np.sum((v3s-v1s)**2.0,axis=1)**0.5 self.dampening = 0.999 def add_force(self,force): self.vertex_speeds += np.array(force) def add_position(self,pos): self.delta_pos += np.array(pos) def move(self): #Length tensors v1s = self.vertices[self.vertex_indices[0::3]] v2s = self.vertices[self.vertex_indices[1::3]] v3s = self.vertices[self.vertex_indices[2::3]] side1 = v2s - v1s side2 = v3s - v2s side3 = v1s - v3s delta_side1 = np.transpose( ((np.sum(side1**2.0,axis=1)**0.5)-self.target_lengths1,)*3 ) * side1 delta_side2 = np.transpose( ((np.sum(side2**2.0,axis=1)**0.5)-self.target_lengths2,)*3 ) * side2 delta_side3 = np.transpose( ((np.sum(side3**2.0,axis=1)**0.5)-self.target_lengths3,)*3 ) * side3 tensor = 0.01 self.vertex_speeds[self.vertex_indices[0::3]] += tensor*delta_side1 self.vertex_speeds[self.vertex_indices[1::3]] -= tensor*delta_side1 self.vertex_speeds[self.vertex_indices[1::3]] += tensor*delta_side2 self.vertex_speeds[self.vertex_indices[2::3]] -= tensor*delta_side2 self.vertex_speeds[self.vertex_indices[2::3]] += tensor*delta_side3 self.vertex_speeds[self.vertex_indices[0::3]] -= tensor*delta_side3 self.delta_pos += self.vertex_speeds self.vertex_speeds *= self.dampening self.vertices += self.delta_pos self.delta_pos.fill(0.0) def collision_detect(self): indices = np.where(self.vertices[:,1]<0.0001) self.vertices[indices,1] = 0.0001 self.vertex_speeds[indices] *= -0.1 def get_data(self): return [self.vertices[self.vertex_indices].tolist()] The code works fairly well for its complexity (and runs very quickly). I believe it also needs angle tensors, which would be more complicated. There are a few strange errors I can't account for (like the object very slowly drifts, as it rests on the ground). At the moment, I'm wondering whether a physics engine has already been done with NumPy before, and if not, whether anyone has any interest in making one. Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Smoothing a Grid
Hi, So I have a 4D grid (x,y,z,vec3), where the vec3 represents a position. What I want to do is to move vec3 elements of the grid based on surrounding vec3 elements so that the grid's values overall are more orthogonal. For example, consider the following 3 x 3 x 3 x vec3 grid: -1.0,-1.0,-1.0],[ 0.0,-1.0,-1.0],[ 1.0,-1.0,-1.0]], [[-1.0,-1.0, 0.0],[ 0.0,-1.0, 0.0],[ 1.0,-1.0, 0.0]], [[-1.0,-1.0, 1.0],[ 0.0,-1.0, 1.0],[ 1.0,-1.0, 1.0]]], [[[-1.0, 0.0,-1.0],[ 0.0, 0.0,-1.0],[ 1.0, 0.0,-1.0]], [[-1.0, 0.0, 0.0],[ 0.2, 0.0, 0.0],[ 1.0, 0.0, 0.0]], [[-1.0, 0.0, 1.0],[ 0.0, 0.0, 1.0],[ 1.0, 0.0, 1.0]]], [[[-1.0, 1.0,-1.0],[ 0.0, 1.0,-1.0],[ 1.0, 1.0,-1.0]], [[-1.0, 1.0, 0.0],[ 0.0, 1.0, 0.0],[ 1.0, 1.0, 0.0]], [[-1.0, 1.0, 1.0],[ 0.0, 1.0, 1.0],[ 1.0, 1.0, 1.0 Notice that, with the exception of the middle value, [0.2,0.0,0.0], the vec3 elements in the grid define a 3 x 3 x 3 grid. The middle value should thus be moved perhaps 10% of the distance to the target location [0.0,0.0,0.0], which can be inferred from the surrounding vec3 elements. After this movement, the value is [0.18,0.00,0.00]. To be more precise, the target position of each vec3 is determined from the surrounding values. Every vec3 moves towards being a certain x_distance from adjacent vec3s in the x direction. Every vec3 moves towards being a certain y_distance from adjacent vec3s in the y direction. Every vec3 moves towards being a certain z_distance from adjacent vec3s in the z direction. Each vec3 also moves towards a certain diagonal_distance from all vec3s along a major diagonal. I am trying to accomplish this. Any ideas? Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] list to array is slow
On Thu, Jul 22, 2010 at 2:09 PM, marco cammarata wrote: > To convert the list into an array takes about 5 sec ... > Not too familiar with typical speeds, but at a guess, perhaps because it must convert 61.4 million (640*480*200) values? Just to *count* that high with xrange takes 1.6 seconds for my computer. I imagine doing further operations (like converting into a numpy array) on other Python objects (like lists) might be similarly slow. Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Finding Unique Pixel Values
Hi, So, I'm working on a radiosity renderer, and it's basically finished. I'm now trying to optimize it. Currently, by far the most computationally expensive operation is visibility testing, where pixels are counted by the type of patch that was drawn on them. Here's my current code, which I'm fairly sure can be optimized: line 1) pixel_data = np.cast['int16'](255.0*glReadPixels(0,0,patch_view.width,patch_view.height,GL_RGB,GL_FLOAT)).reshape((-1,3)) line 2) pixel_data[:,1] *= 255 line 3) pixel_data[:,0] *= (255*255) line 4) pixel_data = np.sum(pixel_data,1) line 5) for pixel in pixel_data: line 6)if pixel != 0: line 7)try:self.visibility[poly1num][pixel-1] line 8)except: self.visibility[poly1num][pixel-1] = 0 line 9)self.visibility[poly1num][pixel-1] += 1 The basic idea: Convert unique an array of RGB triplets into and array of unique numbers. The triplet [0,0,1] maps to 1 and the triplet [255,255,255] maps to 16581375. The triplet [0,0,0], mapping to 0, is always discarded. The resulting array of integers in the range [1,16581375] are the (index+1)'s of a list, self.visibility[poly1num]. For each integer "i" in the array, the value in self.visibility[poly1num][i-1] must be incremented by 1. By line: line 1: convert the read pixels from the range [0.0,1.0] (as floats) to the range [0,255] (as integers). Change the shape from (width,height,3) to (width*height,3). lines 2, 3, and 4: the pixels are an index to patches. The green channel counts 255 times as much. The red, 255 times as much as the green. This way, 255^3 possible pixel types can be had. Sum the red, green, and blue channels to arrive at a single number, which is the patch index in the range [0,255^3]. line 5: for each index in the data line 6: I'm reserving the index [0,0,0] for no patch, so if the index is that color, no patch was drawn. lines 7, 8, and 9: increment self.visibility[poly1num][pixel-1] by 1. I have a feeling that this code is likely the bottleneck. It's taking about an hour to execute this process 8720 times when width and height are both 512. I have a feeling my code can be greatly improved in many ways--if not all over, certainly in lines 5 through 9. Help? Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Finding Unique Pixel Values
Hi again, I've condensed the problem down a lot, because I both presented it in an overcomplicated way, and did not explain it particularly well. Condensed problem: a = np.zeros(num_patches) b = np.array(...) #created, and is size 512^512 = 262,144 #Each value in "b" is an index into "a". #For each occurrence of a given index in "b", increment the corresponding value in "a" by +1. #How should I do that? Simpler, and more importantly, better explained problem. Help? Thanks again, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Finding Unique Pixel Values
On Thu, Jul 22, 2010 at 10:05 PM, Charles R Harris < charlesr.har...@gmail.com> wrote: > Is that what you want, or do you just want to know how many unique indices > there are? As to encoding the RGB, unless there is a existing program your > best bet is probably to use a dot product, i.e., if pix_data is (n,3) uint8, > then > > pix_index = np.dot(pix_data, array([1, 2**8, 2**16], dtype=uint32) > > then check out the documentation for np.unique. > I like the dot product idea for the indices. To the second, actually, I need to increment the number of times the index is there. For example, if b=[1,5,6,6,6,9], then a[6-1] would have to be incremented by +3 = +1+1+1. I tried simply a[b-1]+=1, but it seems to only increment values once, even if there are more than one occurrence of the index in "b". What to do about that? Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Finding Unique Pixel Values
On Fri, Jul 23, 2010 at 12:13 AM, Jon Wright wrote: > Ian Mallett wrote: > > > > To the second, actually, I need to increment the number of times the > > index is there. For example, if b=[1,5,6,6,6,9], then a[6-1] would have > > to be incremented by +3 = +1+1+1. I tried simply a[b-1]+=1, but it > > seems to only increment values once, even if there are more than one > > occurrence of the index in "b". What to do about that? > > Is this what you mean? > > >>> numpy.bincount( [1,5,6,6,6,9] ) > array([0, 1, 0, 0, 0, 1, 3, 0, 0, 1]) > Well, I've implemented it thusly: pixel_data = np.cast['int16'](255.0*glReadPixels(0,0,patch_view.width,patch_view.height,GL_RGB,GL_FLOAT)).reshape((-1,3)) pixel_data = np.dot(pixel_data,np.array([65025,255,1])) count = np.bincount(pixel_data)[1:] self.visibility[poly1num][0:count.shape[0]] += count And it works! Thanks everyone! Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Arrays of Python Values
Hi,
So working on the radiosity renderer:
http://a.imageshack.us/img186/2479/image2f.png.
The code now runs fast enough to generate the data required to draw that.
Now, I need to optimize the radiosity calculation, so that it will converge
in a reasonable amount of time. Right now, the individual blocks
("patches") are stored as instances of Python classes in a list.
I'd heard that NumPy supports arrays of Python objects, so, I simply made an
array out of the list, which seemed ok. Unfortunately, there is a custom
sorting operation that sorted the list of by an attribute of each class:
self.patches.sort( lambda x,y:cmp(x.residual_radiance,y.residual_radiance),
reverse=True )
Because I've never used arrays of Python objects (and Googling didn't turn
up any examples), I'm stuck on how to sort the corresponding array in NumPy
in the same way.
Of course, perhaps I'm just trying something that's absolutely impossible,
or there's an obviously better way. I get the feeling that having no Python
objects in the NumPy array would speed things up even more, but I couldn't
figure out how I'd handle the different attributes (or specifically, how to
keep them together during a sort).
What're my options?
Thanks again for the extremely valuable help,
Ian
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Control format for array of integers
How about numpy.set_printoptions(suppress=True)? See http://docs.scipy.org/doc/numpy/reference/generated/numpy.set_printoptions.html . HTH, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Sine Distribution on a 2D Array
Hi, So I have a square 2D array, and I want to fill the array with sine values. The values need to be generated by their coordinates within the array. The center of the array should be treated as the angle 90 degrees. Each of the four edges should be 0 degrees. The corners, therefore, ought to be -sqrt(2)*90 degrees. The angle is equal to (distance_from_center/(dimension_of_array/2))*90 degrees. Then take the sine of this angle. To describe another way, if the array is treated like a height-field, a single mound of the sine wave should just fit inside the array. Right now, I'm having trouble because I don't know how to operate on an array's values based on the index of the values themselves. Help? Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Sine Distribution on a 2D Array
Hi, After much deliberation, I found a passable solution: distances = np.abs(np.arange(0,resolution,1)+0.5-(resolution/2.0)) x_gradient = np.tile(distances,(resolution,1)) y_gradient = np.copy(x_gradient) y_gradient = np.swapaxes(y_gradient,0,1) distances_to_center = np.hypot(x_gradient,y_gradient) angles = np.radians( (distances_to_center/((resolution/2.0)-0.5)) * 90.0 ) lambert_weights = np.cos(angles) This seems like it could be improved. Ian PS, I meant "Sinusoidal" in the title. I carried that over into my description, which should have been "cosine", instead of "sine". ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Arrays of Python Values
Hi, I've converted all of the code to use record arrays, for a 10-fold speed boost. Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] two dimensional array of sets
Hi, You can use np.mgrid to construct a grid of coordinates. From there, you can make your new array based on these coordinates however you like. Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] summing over more than one axis
Hi, Couldn't you do it with several sum steps? E.g.: result = array.sum(axis=1).sum(axis=2) Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] point line distance
Hi, A background in linear algebra helps. I just came up with this method (which, because I thought of it 5 seconds ago, I don't know if it works): Line p1, p2 Point v costheta = normalize(p2-p1) dot normalize(v-p1) dist = length(v-p1)*sin(acos(costheta) Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Numpy Positional Array
Hello, I'm trying to make an array of size n*n*2. It should be of the form: [[[0,0],[1,0],[2,0],[3,0],[4,0], ... ,[n,0]], [[0,1],[1,1],[2,1],[3,1],[4,1], ... ,[n,1]], [[0,2],[1,2],[2,2],[3,2],[4,2], ... ,[n,2]], [[0,3],[1,3],[2,3],[3,3],[4,3], ... ,[n,3]], [[0,4],[1,4],[2,4],[3,4],[4,4], ... ,[n,4]], ... ... ... ... ... ... ... [[0,n],[1,n],[2,n],[3,n],[4,n], ... ,[n,n]]] Each vec2 represents the x,y position of the vec in the array itself. I'm completely stuck on implementing this setup in numpy. Any pointers? Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Numpy Positional Array
On Tue, Mar 31, 2009 at 5:32 PM, Robert Kern wrote: > How do you want to fill in the array? If you are typing it in > literally into your code, you would do basically the above, without > the ...'s, and wrap it in numpy.array(...). I know that, but in some cases, n will be quite large, perhaps 1000 on a side. I'm trying to generate an array of that form in numpy entirely for speed and aesthetic reasons. Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Numpy Positional Array
The array follows a pattern: each array of length 2 represents the x,y index of that array within the larger array. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Numpy Positional Array
Thanks! ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Numpy Positional Array
Same. Thanks, too. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Specially Constructed Arrays
Hello, I want to make an array of size sqrt(n) by sqrt(n) by 3, filled with special values. The values range from 0.0 to 3.0, starting with 0.0 at one corner and ending at 3.0 in the opposite, increasing going row by row. The value is to be encoded in each color. Because this is somewhat abstract, here's a small example (n=25), generated using the attached code (it also multiplies the number by 255 to obtain a RGB color and not messy floats) to show the concept. The real version should be done by NumPy. This is where I need help; I have no idea how to even approach the problem. [[ 0, 0, 0],[ 32, 0, 0],[ 64, 0, 0],[ 96, 0, 0],[128, 0, 0], [159, 0, 0],[191, 0, 0],[223, 0, 0],[255, 0, 0],[255, 32, 0], [255, 64, 0],[255, 96, 0],[255,128, 0],[255,159, 0],[255,191, 0], [255,223, 0],[255,255, 0],[255,255, 32],[255,255, 64],[255,255, 96], [255,255,128],[255,255,159],[255,255,191],[255,255,223],[255,255,255]] Arrays like this need to be generated quite quickly, so the per-pixel method I presented will not work. How should I do it with NumPy? Thanks, Ian def clamp(num,low,high): if numhigh: return high return num def clamp_vec3(vec3,low,high): return [clamp(vec3[0],low,high), clamp(vec3[1],low,high), clamp(vec3[2],low,high)] def rndint(num): return int(round(num)) def rndint_vec3(vec3): return [rndint(vec3[0]), rndint(vec3[1]), rndint(vec3[2])] def mult_vec3(sc,vec3): return [sc*vec3[0], sc*vec3[1], sc*vec3[2]] n = 25 intensity = 0.0 step = 3.0/(n-1) surf = [] for x in xrange(n): color = clamp_vec3([intensity,intensity-1.0,intensity-2.0],0.0,1.0) color = mult_vec3(255.0,color) color = rndint_vec3(color) surf.append(color) intensity += step surf = str(surf) surf2 = "" for char in surf: if char != " ": surf2 += char print surf2 ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Another Array
Hello, With the help of this list, over the past two days, I have implemented a GPU particle system (picture here: http://img7.imageshack.us/img7/7589/image1olu.png). Every particle is updated entirely on the GPU, with texture data (arrays) being updated iteratively between two framebuffer objects. Here, 65,536 particles are drawn at ~50 fps. Currently, jitter is introduced into the system by means of a "random texture", which should store random unit 3D vectors encoded into the RGB channels (vectors point in all directions, positive and negative, so the vecotr is actually 0.5 long and added to 0.5). The random texture is constructed using Python's random number generation, which is noticeably slower. I also didn't find an efficient way of creating "random" vectors. All this contributes to the pattern artifact at the bottom of the screenshot. I'd like to use NumPy. The array should be constructed as n*n*3, with every vec3 being a normalized random vector with each component in the range [-1.0,1.0]. I could use numpy.random.random(), but this would not give the appropriate vectors. Every 3D vector should have the same probability of being chosen. I don't know how to do that in NumPy--that's where I need help. Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Another Array
This works; thank you. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Another Array
It gives a perfect parabolic shape that looks very nice, but somewhat unrealistic. I'd like to scale the unit vectors by a random length (which can just be a uniform distribution). I tried scaling the unit vector n*n*3 array by a random n*n array, but that didn't work, obviously. Help? ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Another Array
On Thu, Apr 9, 2009 at 11:46 PM, Robert Kern wrote: > Parabolic? They should be spherical. The particle system in the last screenshot was affected by gravity. In the absence of gravity, the results should be spherical, yes. All the vectors are a unit length, which produces a perfectly smooth surface (unrealistic for such an effect). > No, it's not obvious. Exactly what code did you try? What results did > you get? What results were you expecting? It crashed. I have this code: vecs = Numeric.random.standard_normal(size=(self.size[0],self.size[1],3)) magnitudes = Numeric.sqrt((vecs*vecs).sum(axis=-1)) uvecs = vecs / magnitudes[...,Numeric.newaxis] randlen = Numeric.random.random((self.size[0],self.size[1])) randuvecs = uvecs*randlen #It crashes here with a dimension mismatch rgb = ((randvecs+1.0)/2.0)*255.0 I also tried randlen = Numeric.random.random((self.size[0],self.size[1],3)), but this does not scale each of the vector's components equally, producing artifacts again. Each needs to be scaled by the same random value for it to make sense. > Let's take a step back. What kind of distribution are you trying to > achieve? You asked for uniformly distributed unit vectors. Now you are > asking for something else, but I'm really not sure what. What standard > are you comparing against when you say that the unit vectors look > "unrealistic"? The vectors are used to "jitter" each particle's initial speed, so that the particles go in different directions instead of moving all as one. Using the unit vector causes the particles to make the smooth parabolic shape. The jitter vectors much then be of a random length, so that the particles go in all different directions at all different speeds, instead of just all in different directions. Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Another Array
This seems to work: vecs = Numeric.random.standard_normal(size=(self.size[0],self.size[1],3)) magnitudes = Numeric.sqrt((vecs*vecs).sum(axis=-1)) uvecs = vecs / magnitudes[...,Numeric.newaxis] randlen = Numeric.random.random((self.size[0],self.size[1])) randuvecs = uvecs*randlen[...,Numeric.newaxis] rgb = ((randuvecs+1.0)/2.0)*255.0 (I have "import numpy as Numeric" for other reasons, that's why there's Numeric there). Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Distance Formula on an Array
Hi, I have an array sized n*3. Each three-component is a 3D position. Given another 3D position, how is the distance between it and every three-component in the array found with NumPy? So, for example, if the array is: [[0,0,0],[0,1,0],[0,0,3]] And the position is: [0,4,0] I need this array out: [4,3,5] (Just a simple Pythagorean Distance Formula) Ideas? Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
On Sat, Apr 25, 2009 at 12:57 PM, Charles R Harris < charlesr.har...@gmail.com> wrote: > In [3]: vec = array([[0,0,0],[0,1,0],[0,0,3]]) > > In [4]: pos = array([0,4,0]) > > In [5]: sqrt(((vec - pos)**2).sum(1)) > Out[5]: array([ 4., 3., 5.]) > > Chuck > On Sat, Apr 25, 2009 at 1:00 PM, wrote: > if scipy is permitted: > > >>> a = np.array([[0,0,0],[0,1,0],[0,0,3]]) > >>> scipy.spatial.distance_matrix(a, [[0,4,0]]) > array([[ 4.], > [ 3.], > [ 5.]]) > >>> scipy.spatial.minkowski_distance(a, [0,4,0]) > array([ 4., 3., 5.]) > > Josef Thanks you two. I'm going to guess SciPy might be faster (?), but unfortunately it's not going to be available. Thanks, though. Problem solved. Thanks again, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Oops, one more thing. In reference to: vec = array([[0,0,0],[0,1,0],[0,0,3]]) pos = array([0,4,0]) sqrt(((vec - pos)**2).sum(1)) -> array([ 4., 3., 5.]) Can I make "vec" an array of class instances? I tried: class c: def __init__(self): self.position = [0,0,0] vec = array([c(),c(),c()]) pos = array([0,4,0]) sqrt(((vec.position - pos)**2).sum(1)) Which doesn't work. I'm not familiar with class objects in arrays--how should they be referenced? Thanks again, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Well, if it will kill performance, I'm afraid I can't do that. Thanks though. I think it's working now. Now that I have the 1D array of distances, I need the indices of those distances that are less than a number "d". what should I do to do that? Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
It would be: numpy.where(array___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Yeah, I ended up finding the [0] bit at the end through trial and error. I actually do need the indices, though. I'm having a strange new problem though. numpy.array([1,2,3]) is returning a sequence??? I'm really confused. Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Hmmm, I played around with some other code, and it's working right now--not sure what I did... ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Yes, this is pretty much what I'm doing. Right now, I'm having weird troubles with the objects themselves; the objects should and do terminate after a certain time, yet for some reason they're still being drawn. User error, I'm sure. Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
On Sat, Apr 25, 2009 at 1:26 PM, Ian Mallett wrote: > I'm going to guess SciPy might be faster (?), but unfortunately it's not > going to be available. Thanks, though. > ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
The problem is that the object moves too much between frames. A reasonable bounding sphere is 1 for this purpose, but the objects move at least 3. So, given the two arrays, one holding the objects' positions and the other their previous positions, how can I find if, at some point between, the objects collided with the bounding sphere? Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
It was in an error code somewhere. I fixed the problem by messing around with it. I tried the following: a = numpy.array([1, 2, 3]) print a and it gave: [1, 2, 3] instead of: array([1, 2, 3]) Then there were errors about it being a sequence instead of an array somewhere else. Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Using http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html, I came up with: x0 = numpy.array(#point to collide with) x1 = #objects' positions x2 = #objects' previous positions numerator = numpy.sqrt((numpy.cross((x0-x1),(x0-x2))**2).sum(1)) denominator = numpy.sqrt(((x2-x1)**2).sum(1)) dist = numerator/denominator collided_indices = numpy.where(dist<1.0)[0] #1.0==radius Unfortunately, in some cases, it returns indices when it shouldn't--when the objects clearly don't penetrate the sphere. I also get the feeling this will give false positives if the line containing the vector intersects the sphere, but the vector itself does not, though I don't know. Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Thanks, but I don't want to make SciPy a dependency. NumPy is ok though. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Distance Formula on an Array
Hey, this looks cool! I may use it in the future. The problem has already been solved, though, and I don't think changing it is necessary. I'd also like to keep the dependencies (even packaged ones) to a minimum. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Elementary Array Switching
Hi, I have a NumPy array. The array is 3D, n x n x 3. I'm trying to flip the first element of the last dimension with the last. I tried: temp = myarray[:,:,0].copy() myarray[:,:,0] = myarray[:,:,2].copy() myarray[:,:,2] = temp del temp But it doesn't work as expected. I'm definitely not very good at NumPy, so I get the feeling it's something silly I'm doing. What should that code look like? Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Elementary Array Switching
It seems to be working now--I think my problem is elsewhere. Sorry... ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Interleaved Arrays and
Hi, So I'm trying to get a certain sort of 3D terrain working in PyOpenGL. The idea is to get vertex buffer objects to draw a simple 2D plane comprised of many flat polygons, and use a vertex shader to deform that with a heightmap and map that on a sphere. I've managed to do this with a grid (simple points), making the vertex buffer object: threedimensionalgrid = dstack(mgrid[0:size,0:size,0:1])/float(size-1) twodimensionalgrid = threedimensionalgrid.reshape(self.size_squared,3) floattwodimensionalgrid = array(twodimensionalgrid,"f") self.vertex_vbo = vbo.VBO(floattwodimensionalgrid) However, landscapes tend to be, um, solid :D So, the landscape needs to be drawn as quads or triangles. Strips of triangles will be most effective, and the data must be specified to vbo.VBO() in a certain way: n = #blah testlist = [] for x in xrange(n): for y in xrange(n): testlist.append([x,y]) testlist.append([x+1,y]) If "testlist" is an array (i.e., I could go: "array(testlist)"), it works nicely. However, my Python method is certainly improveable with numpy. I suspect the best way is interleaving the arrays [x,y->yn] and [x+1,y->yn] ntimes, but I couldn't figure out how to do that... Help? Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Interleaved Arrays and
Most excellent solutions, thanks! ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Blurring an Array
Hello, I'm working on a program that will draw me a metallic 3D texture. I successfully made a Perlin noise implementation and found that when the result is blurred in one direction, the result actually looks somewhat like brushed aluminum. The plan is to do this for every n*m*3 layer (2D texture) in the 3D texture. My solution to this anisotropic blurring looks like: soften = layerarray.copy() total = 1 for bluriter in xrange(1,5,1): soften[:,bluriter:] += layerarray[:,:-bluriter] soften[:,:-bluriter] += layerarray[:,bluriter:] total += 2 soften /= total Where layerarray is a n*m*3 array, and soften is the array that will be converted into an image with the other 2D images and saved. This code successfully blurs the array in the y-direction. However, it does not do so the way I would like. The blur is accomplished by a simple shift, making the arrays not line up. This leaves space at the edges. When the final soften array is divided by total, those areas are especially dark. Visually, this is unacceptable, and leads to banding, which is especially evident if the texture is repeated. As you can see in this image, which shows about 6 repeats of the texture, http://img13.imageshack.us/img13/5789/image1wgq.png, the dark edges are annoying. I made sure, of course, that the Perlin noise implementation is tileable. The solution to my problem is to make the shifted array wrap around so that its overhang fills in the hole the shift caused it to leave behind. For example, to simulate shifting the texture 8 units up with wrap, the first 8 rows should be removed from the top and added to the bottom. Likewise for columns if the blur goes in that direction. I already tried a couple of times at this, and it's not working. I need a way to blur soften by column and by row. Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Blurring an Array
Sounds like it would work, but unfortunately numpy was one of my dependency constraints. I should have mentioned that. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Blurring an Array
This works perfectly! Is there likewise a similar call for Numeric? ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] numpy/numexpr performance (particle simulation)
As an off-topic solution, there's always the GPU to do the the particle updating. With half decent optimization, I've gotten over a million particles in *real-time*. You could presumably run several of these at the same time to get as many particles as you want. Downside would be ease-of-implementation... ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Even Sphere Volume
Hello, I'm trying to get a cloud particle system working. As part of determining the depth of every point (to get the light that can pass through), it makes sense that the volume should be of even density. The volume is a sphere. Currently, I'm using: vecs = numpy.random.standard_normal(size=(size,size,3)) magnitudes = numpy.sqrt((vecs*vecs).sum(axis=-1)) uvecs = vecs / magnitudes[...,numpy.newaxis] randlen = numpy.random.random((size,size)) randpoints = uvecs*randlen[...,numpy.newaxis] The particles' positions are set to randpoints. However, this creates only an even distribution of *directions *on the sphere. The particles are too closely centered around the center. What I want is a certain number of particles arranged evenly throughout the sphere. How do I do that? Thanks, Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Even Sphere Volume
Presently, it's being rendered using point sprites/VBOs. It's supposed to be for a game I'm working on, but it's a good effect to have in the toolbox too :D ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Even Sphere Volume
@Stéfan: I thought of the first method. Let's hear the second approach. @Gökhan: Yes. Toolbox is my metaphor for being able to do effects in OpenGL. Point sprites are images mapped onto vertices, VBOs are *v*ertex *b *uffer *o*bjects, that make stuff faster. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Even Sphere Volume
Thanks for the example, Stéfan! I'm trying to work this into a position texture, and the arrays need interleaving. I tried a couple times. Here's what I have: az = numpy.random.uniform(0, numpy.pi * 2, size*size) el = numpy.random.uniform(0, numpy.pi, size*size) r = numpy.random.uniform(size*size) # Transform random variable r = r ** (1/3.) * 1 el = numpy.arccos(1 - 2*el) x = r * numpy.cos(az) * numpy.sin(el); y = r * numpy.sin(az) * numpy.sin(el); z = r * numpy.cos(el) randvecs = numpy.array([x,y,z]).reshape((size,size,3)) rgb = ((randvecs+1.0)/2.0)*255.0 surface = pygame.surfarray.make_surface(rgb) I think the arrays are being concatenated. Hence, visually, it doesn't look right... Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Even Sphere Volume
That didn't fix it. I messed around some more, but I couldn't get the spherical coordinates working. I decided to rework my first method. By raising the radius to the one third power, like for the other method, basically the same thing is accomplished. It's working now, thanks. :D vecs = numpy.random.standard_normal((size,size,3)) magnitudes = numpy.sqrt((vecs*vecs).sum(axis=-1)) vecs = vecs / magnitudes[...,numpy.newaxis] randlen = numpy.random.random((size,size)) randlen = randlen ** (1.0/3.0) randpoints = vecs*randlen[...,numpy.newaxis] rgb = ((randpoints+1.0)/2.0)*255.0 Ian ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Even Sphere Volume
These are all great algorithms, thanks for the help! ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] Optimizing a pure Python Workaround
Hello, I have some code that makes vertex buffer object terrain. Because the setup for this (a series of triangle strips) is a bit clunky, I just implemented the tricky parts in Python. The code works, but it's slow. How should I go about optimizing it? Thanks, Ian size = [#int_something,#int_something_else] fsize = map(float,size) array = [] for x in xrange(size[0]): for y in xrange(size[1]): array.append(np.array([ x /fsize[0],y/fsize[1],0.0],"f")) array.append(np.array([(x+1)/fsize[0],y/fsize[1],0.0],"f")) array = np.array(array) #"array" is the final product. I need this final product generated more quickly. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Optimizing a pure Python Workaround
Hi, Sorry, I've been away in Oregon... The result isn't quite the same. The arrays must be in the range [0,1], so I just have it divide x3 and y. I also have it add 1 to size[1], as I realized that was also necessary for that behavior: x = np.arange(size[0]) x2 = np.column_stack([x,x+1]).reshape([-1,1]) x3 = np.array(x2.repeat(size[1]+1,axis=1).flatten(),"f") y = np.array(np.arange(size[1]+1).repeat(size[0]*2),"f") array = np.zeros([len(y),3]) array[:,0] = x3/size[0] array[:,1] = y/size[1] array = np.array(array,"f") When size is [3,2], the result from this code is: [[ 0. 0. 0.] [ 0. 0. 0.] [ 0. 0. 0.] [ 0.3334 0. 0.] [ 0.3334 0. 0.] [ 0.3334 0. 0.] [ 0.3334 0.5 0.] [ 0.3334 0.5 0.] [ 0.3334 0.5 0.] [ 0.6669 0.5 0.] [ 0.6669 0.5 0.] [ 0.6669 0.5 0.] [ 0.6669 1. 0.] [ 0.6669 1. 0.] [ 0.6669 1. 0.] [ 1. 1. 0.] [ 1. 1. 0.] [ 1. 1. 0.]] The correct output is: [[ 0. 0. 0.] [ 0.3334 0. 0.] [ 0. 0.5 0.] [ 0.3334 0.5 0.] [ 0. 1. 0.] [ 0.3334 1. 0.] [ 0.3334 0. 0.] [ 0.6669 0. 0.] [ 0.3334 0.5 0.] [ 0.6669 0.5 0.] [ 0.3334 1. 0.] [ 0.6669 1. 0.] [ 0.6669 0. 0.] [ 1. 0. 0.] [ 0.6669 0.5 0.] [ 1. 0.5 0.] [ 0.6669 1. 0.] [ 1. 1. 0.]] Thanks, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Is this a bug?
Awww, it's fun to be foolish on Fridays! ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] GPU Numpy
On Wed, Aug 5, 2009 at 11:34 AM, Charles R Harris wrote: > It could be you could slip in a small mod that would do what you want. I'll help, if you want. I'm good with GPUs, and I'd appreciate the numerical power it would afford. > The main problems with using GPUs were that CUDA was only available for > nvidia video cards and there didn't seem to be any hope for a CUDA version > of LAPACK. You don't have to use CUDA, although it would make it easier. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Recarray comparison and byte order
On Sat, Oct 31, 2009 at 9:38 AM, Matthew Brett wrote: > c = a.byteswap().newbyteorder() > c == a > In the last two lines, a variable "c" is assigned to a modified "a". The next line tests (==) to see if "c" is the same as (==) the unmodified "a". It isn't, because "c" is the modified "a". Hence, "False". Do you mean: c = a instead of: c == a ...? HTH, Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Fitting a curve on a log-normal distributed data
Theory wise: -Do a linear regression on your data. -Apply a logrithmic transform to your data's dependent variable, and do another linear regression. -Apply a logrithmic transform to your data's independent variable, and do another linear regression. -Take the best regression (highest r^2 value) and execute a back transform. Then, to get your desired extrapolation, simply substitute in the size for the independent variable to get the expected value. If, however, you're looking for how to implement this in NumPy or SciPy, I can't really help :-P Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Fitting a curve on a log-normal distributed data
Hello, My analysis shows that the exponential regression gives the best result (r^2=87%)--power regression gives worse results (r^2=77%). Untransformed data gives r^2=76%. I don't think you want lognorm. If I'm not mistaken, that fits the data to a log(normal distribution random variable). So, take the logarithm (to any base) of all the "conc" values. Then do a linear regression on those values versus "sizes". Try (semi-pseudocode): slope, intercept, p, error = scipy.stats.linregress(sizes,log(conc)) Ian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
