4

I have a 2D array, a, comprising a set of 100 x,y,z coordinates:

[[ 0.81  0.23  0.52]
 [ 0.63  0.45  0.13]
 ...
 [ 0.51  0.41  0.65]]

I would like to create a 3D binary image, b, with 101 pixels in each of the x,y,z dimensions, of coordinates ranging between 0.00 and 1.00. Pixels at locations defined by a should take on a value of 1, all other pixels should have a value of 0.

I can create an array of zeros of the right shape with b = np.zeros((101,101,101)), but how do I assign coordinate and slice into it to create the ones using a?

Ulrich Stern
  • 10,761
  • 5
  • 55
  • 76
MyCarta
  • 808
  • 2
  • 12
  • 37

2 Answers2

4

You could do something like this -

# Get the XYZ indices
idx = np.round(100 * a).astype(int)

# Initialize o/p array
b = np.zeros((101,101,101))

# Assign into o/p array based on linear index equivalents from indices array
np.put(b,np.ravel_multi_index(idx.T,b.shape),1)

Runtime on the assignment part -

Let's use a bigger grid for timing purposes.

In [82]: # Setup input and get indices array
    ...: a = np.random.randint(0,401,(100000,3))/400.0
    ...: idx = np.round(400 * a).astype(int)
    ...: 

In [83]: b = np.zeros((401,401,401))

In [84]: %timeit b[list(idx.T)] = 1 #@Praveen soln
The slowest run took 42.16 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 6.28 ms per loop

In [85]: b = np.zeros((401,401,401))

In [86]: %timeit np.put(b,np.ravel_multi_index(idx.T,b.shape),1) # From this post
The slowest run took 45.34 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 5.71 ms per loop

In [87]: b = np.zeros((401,401,401))

In [88]: %timeit b[idx[:,0],idx[:,1],idx[:,2]] = 1 #Subscripted indexing
The slowest run took 40.48 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 6.38 ms per loop
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • 1
    `(a * 100).astype(int)` is extremely dangerous! You could get erroneous indices! Try: `a = np.arange(101) / 100`, `b = (100 * a).astype(int)`, `(b == np.arange(101, dtype=int)).all()`, and you'll get `False`! The number 28 shows up twice in `b`, for instance, due to floating point precision issues. – Praveen Aug 30 '16 at 06:10
  • @Praveen Damn, you are right! Rounding would be the only way here. Updated the post. Now, it's similar to yours at the start. Do you want me to delete this post? – Divakar Aug 30 '16 at 06:21
  • I wouldn't ask you to delete it. Your `np.put` method is interesting. I'd be interested in a timing analysis between that and just using `list(a_indices.T)`. – Praveen Aug 30 '16 at 06:23
  • @Praveen Sure, let me add that. – Divakar Aug 30 '16 at 06:24
  • @Praveen Added. Seems the performance figures are comparable. – Divakar Aug 30 '16 at 06:51
  • Thanks! I tried the same thing myself, but the intermediate caching was putting me off... While this is probably a separate question, is there a way to avoid that, that you know of? These arrays seem too big to be in CPU cache as [this](http://stackoverflow.com/a/29765693/525169) answer suggests... – Praveen Aug 30 '16 at 06:54
  • @Praveen Hmm not that I know about off the top of my head. – Divakar Aug 30 '16 at 06:58
4

First, start off by safely rounding your floats to ints. In context, see this question.

a_indices = np.rint(a * 100).astype(int)

Next, assign those indices in b to 1. But be careful to use an ordinary list instead of the array, or else you'll trigger the usage of index arrays. It seems as though performance of this method is comparable to that of alternatives (Thanks @Divakar! :-)

b[list(a_indices.T)] = 1

I created a small example with size 10 instead of 100, and 2 dimensions instead of 3, to illustrate:

>>> a = np.array([[0.8, 0.2], [0.6, 0.4], [0.5, 0.6]])
>>> a_indices = np.rint(a * 10).astype(int)
>>> b = np.zeros((10, 10))
>>> b[list(a_indices.T)] = 1
>>> print(b) 
[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]
Community
  • 1
  • 1
Praveen
  • 6,872
  • 3
  • 43
  • 62