Fast histogram

classic Classic list List threaded Threaded
13 messages Options
Reply | Threaded
Open this post in threaded view
|

Fast histogram

Zachary Pincus-2
Hi folks,

I'm working on a live-video display for some microscope control tools  
I'm building. For this, I need a  fast histogram function to work on  
large-ish images (1000x2000 or so) at video rate, with cycles left  
over for more interesting calculations (like autofocus).

Now, numpy.histogram is a bit slower than I'd like, probably because  
it's pretty general (and of course cf. the recent discussion about its  
speed). I just need even bins within a set range. This is easy enough  
to do with a C-extension, or perhaps even cython, but before I go  
there, I was wondering if there's a numpy function that can help.

Here's what I have in mind:

def histogram(arr, bins, range):
   min, max = range
   indices = numpy.clip(((arr.astype(float) - min) * bins / (max -  
min)).astype(int), 0, bins-1)
   histogram = numpy.zeros(bins, numpy.uint32)
   for i in indices:
     histogram[i] += 1

Now, clearly, the last loop is what needs speeding up. Are there any  
numpy functions that can do this kind of operation? Also perhaps  
unnecessarily slow is the conversion of 'arr' to a float -- I do this  
to avoid overflow issues with integer arrays.

Any advice? Should I go ahead and write this up in C (easy enough), or  
can I do this in numpy? Probably the indices-computation line I'll  
speed up with numexpr, if I use a pure-numpy approach.

Thanks,
Zach
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Charles R Harris


On Thu, Apr 17, 2008 at 10:02 AM, Zachary Pincus <[hidden email]> wrote:
Hi folks,

I'm working on a live-video display for some microscope control tools
I'm building. For this, I need a  fast histogram function to work on
large-ish images (1000x2000 or so) at video rate, with cycles left
over for more interesting calculations (like autofocus).

Now, numpy.histogram is a bit slower than I'd like, probably because
it's pretty general (and of course cf. the recent discussion about its
speed). I just need even bins within a set range. This is easy enough
to do with a C-extension, or perhaps even cython, but before I go
there, I was wondering if there's a numpy function that can help.

Here's what I have in mind:

def histogram(arr, bins, range):
  min, max = range
  indices = numpy.clip(((arr.astype(float) - min) * bins / (max -
min)).astype(int), 0, bins-1)
  histogram = numpy.zeros(bins, numpy.uint32)
  for i in indices:
    histogram[i] += 1

Now, clearly, the last loop is what needs speeding up. Are there any
numpy functions that can do this kind of operation? Also perhaps
unnecessarily slow is the conversion of 'arr' to a float -- I do this
to avoid overflow issues with integer arrays.

How about a combination of sort, followed by searchsorted right/left using the bin boundaries as keys? The difference of the two resulting vectors is the bin value. Something like:

In [1]: data = arange(100)

In [2]: bins = [0,10,50,70,100]

In [3]: lind = data.searchsorted(bins)

In [4]: print lind[1:] - lind[:-1]
[10 40 20 30]

This won't be as fast as a c implementation, but at least avoids the loop.

Chuck



_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Sebastian Haase
On Thu, Apr 17, 2008 at 6:18 PM, Charles R Harris
<[hidden email]> wrote:

>
>
> On Thu, Apr 17, 2008 at 10:02 AM, Zachary Pincus <[hidden email]>
> wrote:
> > Hi folks,
> >
> > I'm working on a live-video display for some microscope control tools
> > I'm building. For this, I need a  fast histogram function to work on
> > large-ish images (1000x2000 or so) at video rate, with cycles left
> > over for more interesting calculations (like autofocus).
> >
> > Now, numpy.histogram is a bit slower than I'd like, probably because
> > it's pretty general (and of course cf. the recent discussion about its
> > speed). I just need even bins within a set range. This is easy enough
> > to do with a C-extension, or perhaps even cython, but before I go
> > there, I was wondering if there's a numpy function that can help.
> >
> > Here's what I have in mind:
> >
> > def histogram(arr, bins, range):
> >   min, max = range
> >   indices = numpy.clip(((arr.astype(float) - min) * bins / (max -
> > min)).astype(int), 0, bins-1)
> >   histogram = numpy.zeros(bins, numpy.uint32)
> >   for i in indices:
> >     histogram[i] += 1
> >
> > Now, clearly, the last loop is what needs speeding up. Are there any
> > numpy functions that can do this kind of operation? Also perhaps
> > unnecessarily slow is the conversion of 'arr' to a float -- I do this
> > to avoid overflow issues with integer arrays.
> >
>
> How about a combination of sort, followed by searchsorted right/left using
> the bin boundaries as keys? The difference of the two resulting vectors is
> the bin value. Something like:
>
> In [1]: data = arange(100)
>
> In [2]: bins = [0,10,50,70,100]
>
>  In [3]: lind = data.searchsorted(bins)
>
> In [4]: print lind[1:] - lind[:-1]
>  [10 40 20 30]
>
> This won't be as fast as a c implementation, but at least avoids the loop.
>
> Chuck

How many bits per pixel does your camera actually generate !?
If its for example a 12 bit camera, you could just fill in directly
into 4096 preallocated bins.
You would not need any sorting !!
That's what I did for a 16 bit camera -- but I wrote it in C and I had
4 cameras at 30 Hz.

-Sebastian Haase
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Zachary Pincus-2
In reply to this post by Charles R Harris
Hi,

> How about a combination of sort, followed by searchsorted right/left  
> using the bin boundaries as keys? The difference of the two  
> resulting vectors is the bin value. Something like:
>
> In [1]: data = arange(100)
>
> In [2]: bins = [0,10,50,70,100]
>
> In [3]: lind = data.searchsorted(bins)
>
> In [4]: print lind[1:] - lind[:-1]
> [10 40 20 30]
>
> This won't be as fast as a c implementation, but at least avoids the  
> loop.

This is, more or less, what the current numpy.histogram does, no? I  
was hoping to avoid the O(n log n) sorting, because the image arrays  
are pretty big, and numpy.histogram doesn't get close to video rate  
for me...

Perhaps, though, some of the slow-down from numpy.histogram is from  
other overhead, and not the sorting. I'll try this, but I think I'll  
probably just have to write the c loop...

Zach
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Zachary Pincus-2
In reply to this post by Sebastian Haase
Hi, and thanks for the suggestion!

> How many bits per pixel does your camera actually generate !?
> If its for example a 12 bit camera, you could just fill in directly
> into 4096 preallocated bins.
> You would not need any sorting !!
> That's what I did for a 16 bit camera -- but I wrote it in C and I had
> 4 cameras at 30 Hz.

That approach avoids the bin-index calculation line:
indices = numpy.clip(((array.astype(float) - min) * bins / (max -  
min)).astype(int), 0, bins-1)

But even if indices = array, one still needs to do something like:
for index in indices: histogram[index] += 1

Which is slow in python and fast in C.

I'm guessing that there's no utility function in numpy that does a  
loop like this? If so, that would be handy, but if now, I think I need  
to dig out the numpy book and write a little extension...

Zach
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Sebastian Haase
On Thu, Apr 17, 2008 at 6:54 PM, Zachary Pincus <[hidden email]> wrote:

> Hi, and thanks for the suggestion!
>
>
>  > How many bits per pixel does your camera actually generate !?
>  > If its for example a 12 bit camera, you could just fill in directly
>  > into 4096 preallocated bins.
>  > You would not need any sorting !!
>  > That's what I did for a 16 bit camera -- but I wrote it in C and I had
>  > 4 cameras at 30 Hz.
>
>  That approach avoids the bin-index calculation line:
>  indices = numpy.clip(((array.astype(float) - min) * bins / (max -
>
> min)).astype(int), 0, bins-1)
>
>  But even if indices = array, one still needs to do something like:
>  for index in indices: histogram[index] += 1
>
>  Which is slow in python and fast in C.
>
>  I'm guessing that there's no utility function in numpy that does a
>  loop like this? If so, that would be handy, but if now, I think I need
>  to dig out the numpy book and write a little extension...
>
>
I thought of a broadcasting approach... what are the chances that a simple

bins[:] = 0
bins[ img.flat ] += 1

works !?  Or something along those lines ....
-Sebastian
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Zachary Pincus-2
>> But even if indices = array, one still needs to do something like:
>> for index in indices: histogram[index] += 1
>>
>> Which is slow in python and fast in C.
>>
>>
> I thought of a broadcasting approach... what are the chances that a  
> simple
>
> bins[:] = 0
> bins[ img.flat ] += 1

That doesn't work cumulatively, it seems. Some bins get a value of  
'1', but it never gets incremented past that.

It was worth a try, though... in some cases in-place updating seems to  
work this way (which is usually a source of confusion on this list,  
not the desired result!)

Zach
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Jae-Joon Lee
In reply to this post by Zachary Pincus-2
>  But even if indices = array, one still needs to do something like:
>  for index in indices: histogram[index] += 1
>
>  Which is slow in python and fast in C.
>
>  I'm guessing that there's no utility function in numpy that does a
>  loop like this? If so, that would be handy, but if now, I think I need
>  to dig out the numpy book and write a little extension...
>

numpy.bincount?

Docstring:
    bincount(x,weights=None)

    Return the number of occurrences of each value in x.

    x must be a list of non-negative integers.  The output, b[i],
    represents the number of times that i is found in x.  If weights
    is specified, every occurrence of i at a position p contributes
    weights[p] instead of 1.

    See also: histogram, digitize, unique.

-JJ
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Peter Creasey-2
In reply to this post by Zachary Pincus-2
On 17/04/2008, Zachary Pincus wrote:
>  But even if indices = array, one still needs to do something like:
>  for index in indices: histogram[index] += 1
>
>  Which is slow in python and fast in C.
>

I haven't tried this, but if you want the sum in C you could do

for x in unique(indices):
    histogram[x] = (indices==x).sum()

Of course, this just replaces an O(N log N) algorithm by an O(N * M)
(M is the number of bins), which is only going to be worth for very
small M.

Peter
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Fwd: Fast histogram

pec27
On 17/04/2008, Zachary Pincus wrote:
 >  But even if indices = array, one still needs to do something like:
 >  for index in indices: histogram[index] += 1
 >
 >  Which is slow in python and fast in C.
 >

 I haven't tried this, but if you want the sum in C you could do

 for x in unique(indices):
    histogram[x] = (indices==x).sum()

 Of course, this just replaces an O(N log N) algorithm by an O(N * M)
 (M is the number of bins), which is only going to be worth it for very
 small M.


 Peter
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Andrew Straw
In reply to this post by Zachary Pincus-2
Hi Zach,

I have a similar loop which I wrote using scipy.weave. This was my first
foray into weave, and I had to dig through the intermediate C sources to
find the macros that did the indexing in the way I make use of here, but
this snipped may get you started. There are 2 functions, which each do
the same thing, but one is in python, the other is in C. Note that this
is for a 3D histogram -- presumably you could remove B and C from this
example.

I'm sure there are better (more documented) ways to do this using weave
-- but I had this code written, it works, and it appears it may be
useful to you... (Sorry it's not documented, however.)

-Andrew

Zachary Pincus wrote:

> Hi,
>
>  
>> How about a combination of sort, followed by searchsorted right/left  
>> using the bin boundaries as keys? The difference of the two  
>> resulting vectors is the bin value. Something like:
>>
>> In [1]: data = arange(100)
>>
>> In [2]: bins = [0,10,50,70,100]
>>
>> In [3]: lind = data.searchsorted(bins)
>>
>> In [4]: print lind[1:] - lind[:-1]
>> [10 40 20 30]
>>
>> This won't be as fast as a c implementation, but at least avoids the  
>> loop.
>>    
>
> This is, more or less, what the current numpy.histogram does, no? I  
> was hoping to avoid the O(n log n) sorting, because the image arrays  
> are pretty big, and numpy.histogram doesn't get close to video rate  
> for me...
>
> Perhaps, though, some of the slow-down from numpy.histogram is from  
> other overhead, and not the sorting. I'll try this, but I think I'll  
> probably just have to write the c loop...
>
> Zach
> _______________________________________________
> Numpy-discussion mailing list
> [hidden email]
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>  

from scipy import weave
import numpy

def increment_slow(final_array_shape,  A,B,C):
    counts = numpy.zeros(final_array_shape, dtype=numpy.uint64)
    for i in range(len(A)):
        a=A[i]; b=B[i]; c=C[i]
        counts[a,b,c] += 1
    return counts

def increment_fast(final_array_shape,  idxa,idxb,idxc):
    counts = numpy.zeros(final_array_shape, dtype=numpy.uint64)
    assert len(idxa.shape)==1
    assert len(idxa)==len(idxb)
    assert len(idxa)==len(idxc)

    code = r"""
    for (int i=0; i<Nidxa[0]; i++) {
        COUNTS3( IDXA1(i), IDXB1(i), IDXC1(i) )++;
    }
    """
    weave.inline( code, ['counts', 'idxa','idxb','idxc'])
    return counts

#increment = increment_slow
increment = increment_fast

_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Zachary Pincus-2
In reply to this post by Jae-Joon Lee
Hello,

>> But even if indices = array, one still needs to do something like:
>> for index in indices: histogram[index] += 1
> numpy.bincount?

That is indeed what I was looking for! I knew I'd seen such a function.

However, the speed is a bit disappointing. I guess the sorting isn't  
too much of a penalty:

def histogram(array, bins, range):
   min, max = range
   indices = numpy.clip(((array.astype(float) - min) * bins / (max -  
min)).astype(int), 0, bins-1).flat
   return numpy.bincount(indices)

import numexpr
def histogram_numexpr(array, bins, range):
   min, max = range
   min = float(min)
   max = float(max)
   indices = numexpr.evaluate('(array - min) * bins / (max - min)')
   indices = numpy.clip(indices.astype(int), 0, bins-1).flat
   return numpy.bincount(indices)

 >>> arr.shape
(1300, 1030)

 >>> timeit numpy.histogram(arr, 12, [0, 5000])
10 loops, best of 3: 99.9 ms per loop

 >>>  timeit histogram(arr, 12, [0, 5000])
10 loops, best of 3: 127 ms per loop

 >>> timeit histogram_numexpr(arr, 12, [0, 5000])
10 loops, best of 3: 109 ms per loop

 >>>  timeit numpy.histogram(arr, 5000, [0, 5000])
10 loops, best of 3: 111 ms per loop

 >>>  timeit histogram(arr, 5000, [0, 5000])
10 loops, best of 3: 127 ms per loop

 >>> timeit histogram_numexpr(arr, 5000, [0, 5000])
10 loops, best of 3: 108 ms per loop

So, they're all quite close, and it seems that numpy.histogram is the  
definite winner. Huh. I guess I will have to go to C or maybe weave to  
get up to video-rate, unless folks can suggest some further  
optimizations...

Zach
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: Fast histogram

Zachary Pincus-2
Combining Sebastian and Jae-Joon's suggestions, I have something that  
might work:

 >>> timeit numpy.bincount(array.flat)
10 loops, best of 3: 28.2 ms per loop

This is close enough to video-rate... And I can then combine bins as  
needed to get a particular bin count/range after the fact.

Thanks, everyone,

Zach
_______________________________________________
Numpy-discussion mailing list
[hidden email]
http://projects.scipy.org/mailman/listinfo/numpy-discussion