9

I have a 2D numpy array that represents a monochrome image from a CCD that has been binned 3x3 (that is, each value in the array represents 9 pixels (3x3) on the physical CCD).

I want to rescale it to match the original CCD layout (so I can easily overlay it with a non-binned image from the same CCD).

I saw Resampling a numpy array representing an image, but that doesn't seem to do what I want.

Suppose I have an array g:

import numpy as np
import scipy.ndimage

 g = np.array([[0, 1, 2],
               [3, 4, 5],
               [6, 7, 8]])

When I try to scale it by a factor of 2:

o = scipy.ndimage.zoom(g, 2, order=0)

I get exactly what I expect - each value is now 2x2 identical values:

array([[0, 0, 1, 1, 2, 2],
       [0, 0, 1, 1, 2, 2],
       [3, 3, 4, 4, 5, 5],
       [3, 3, 4, 4, 5, 5],
       [6, 6, 7, 7, 8, 8],
       [6, 6, 7, 7, 8, 8]])

But when I try to scale by a factor of 3, I get this:

o = scipy.ndimage.zoom(g, 3, order=0)

Gives me:

array([[0, 0, 1, 1, 1, 1, 2, 2, 2],
       [0, 0, 1, 1, 1, 1, 2, 2, 2],
       [3, 3, 4, 4, 4, 4, 5, 5, 5],
       [3, 3, 4, 4, 4, 4, 5, 5, 5],
       [3, 3, 4, 4, 4, 4, 5, 5, 5],
       [3, 3, 4, 4, 4, 4, 5, 5, 5],
       [6, 6, 7, 7, 7, 7, 8, 8, 8],
       [6, 6, 7, 7, 7, 7, 8, 8, 8],
       [6, 6, 7, 7, 7, 7, 8, 8, 8]])

I wanted each value in the original array to become a set of 3x3 values...that's not what I get.

How can I do it? (And why do I get this unintuitive result?)

Community
  • 1
  • 1
nerdfever.com
  • 1,652
  • 1
  • 20
  • 41
  • 1
    Why the -1? Was this a bad question? Too many tags? As a SO newbie, I'd appreciate some feedback about what I'm doing wrong. – nerdfever.com Sep 05 '14 at 01:06
  • My guess is that someone saw this, knew the answer, assumed you hadn't bothered to read the docs on `zoom`, and therefore downvoted you as a help vampire. Personally, I can understand why you still might be confused after reading the `zoom` docs, but I guess you could have explained better? I don't know. – abarnert Sep 05 '14 at 22:11
  • seems to be a duplicate of http://stackoverflow.com/questions/7525214/how-to-scale-a-numpy-array?rq=1 – Mikhail V Mar 30 '15 at 17:38

2 Answers2

11

You can use np.kron:

In [16]: g
Out[16]: 
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

In [17]: np.kron(g, np.ones((3,3), dtype=int))
Out[17]: 
array([[0, 0, 0, 1, 1, 1, 2, 2, 2],
       [0, 0, 0, 1, 1, 1, 2, 2, 2],
       [0, 0, 0, 1, 1, 1, 2, 2, 2],
       [3, 3, 3, 4, 4, 4, 5, 5, 5],
       [3, 3, 3, 4, 4, 4, 5, 5, 5],
       [3, 3, 3, 4, 4, 4, 5, 5, 5],
       [6, 6, 6, 7, 7, 7, 8, 8, 8],
       [6, 6, 6, 7, 7, 7, 8, 8, 8],
       [6, 6, 6, 7, 7, 7, 8, 8, 8]])

The output of zoom(g, 3, order=0) is a bit surprising. Consider the first row: [0, 0, 1, 1, 1, 1, 2, 2, 2]. Why are there four 1s?

When order=0 zoom (in effect) computes np.linspace(0, 2, 9), which looks like

In [80]: np.linspace(0, 2, 9)
Out[80]: array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ,  1.25,  1.5 ,  1.75,  2.  ])

and then rounds the values. If you use np.round(), you get:

In [71]: np.round(np.linspace(0, 2, 9)).astype(int)
Out[71]: array([0, 0, 0, 1, 1, 1, 2, 2, 2])

Note that np.round(0.5) gives 0, but np.round(1.5) gives 2. np.round() uses the "round half to even" tie-breaking rule. Apparently the rounding done in the zoom code uses the "round half down" rule: it rounds 0.5 to 0 and 1.5 to 1, as in the following

In [81]: [int(round(x)) for x in np.linspace(0, 2, 9)]
Out[81]: [0, 0, 1, 1, 1, 1, 2, 2, 2]

and that's why there are four 1s in there.

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • I know I'm not supposed to post "thanks" comments, but thanks for adding the explanation of why zoom() doesn't do what I expected. Now it makes sense. I'd give you another +1 if I could... – nerdfever.com Sep 05 '14 at 15:12
2

And why do I get this unintuitive result?

Because zoom is a spline interpolation function. In other words, it draws a cubic spline from the midpoint of that 1 to the midpoint of that 0, and the values in between get the values of the spline at the appropriate location.

If you want nearest, linear or quadratic interpolation instead of cubic, you can use the order=0 or order=1 or order=2 argument. But if you don't want interpolation at all—which you don't—don't use an interpolation function. This is like asking why using [int(i*2.3) for i in range(10)] to get even numbers from 0 to 20 gives you some odd numbers. It's not a function to get even numbers from 0 to 20, so it doesn't do that, but it does exactly what you asked it to.


How can I do it?

Again, if you want non-interpolated scaling, don't use an interpolation function. The simplest way is probably to use np.kron, to Kroenecker-multiply your array with np.ones((scale, scale)).

abarnert
  • 354,177
  • 51
  • 601
  • 671
  • Hm. If "nearest is order=0" (per http://stackoverflow.com/questions/13242382/resampling-a-numpy-array-representing-an-image) then the result still seems unintuitive. But I've got my answer. – nerdfever.com Sep 05 '14 at 00:51
  • @nerdfever.com: Yes, nearest is order=0. And the nearest value in the original array to the 3rd value in the expanded array is the 2nd, not the 1st. I was going to explain this in more detail, but I just noticed that Warren Weckesser has already done a great job. Too bad I already upvoted his original answer, because his edited version deserves another one. :) – abarnert Sep 05 '14 at 17:37
  • @abanert: Really? I suppose you must be correct or they wouldn't do it that way, but looking at it one-dimensionally we have [ n1 n2 n3 ] mapping to [m1 m2 m3 m4 m5 m6 m7 m8 m9]. I'd assume (maybe wrongly) that n1 maps to m2 and n2 maps to m5. If so, then m3 is closer to n1 than n2, and should get n1's value. Unless you weigh "nearest" in both dimensions somehow. But I thought "nearest" meant literally "nearest", not "nearest after weighting". [Perhaps I thought wrong.] Even in 2 dimensions, m3 is ADJACENT to m2 and not adjacent to any other mapped value. – nerdfever.com Sep 05 '14 at 19:49
  • @nerdfever.com: Read Warren Weckesser's answer; it explains things better than I can in a comment. – abarnert Sep 05 '14 at 21:54