Efficient "bandmax" function in pyvips?

Question

I'm looking for a way to efficiently find the band that has the maximum value in a large multi-band image.

For example if we had an image that looked like

band 1   band 2   band 3
[1 2 3]  [3 4 5]  [0 1 2]
[3 2 1]  [1 2 3]  [1 0 1]
[4 5 6]  [6 7 5]  [0 0 7]

I'd like to make something like

bandFind = A.bandmax()

Where bandFind would look like

band 1   band 2   band 3
[0 0 0]  [1 1 1]  [0 0 0]
[1 1 0]  [0 0 1]  [0 0 0]
[0 0 0]  [1 1 0]  [0 0 1]

Alternatively if I can get an index image I can pretty easily convert it to what I want.

I wrote a function in python that iterates through the bands, accumulating the max value and the index in a buffer, and the converts it to the kind of map shown above, but it doesn't seem very performant on large images.

  def pickMax(inImg):
    imBandSelect = pyvips.Image.black(inImg.width, inImg.height, bands=1)
    imBandMax = pyvips.Image.black(inImg.width, inImg.height, bands=1)
    for bandIndex, band in enumerate(inImg.bandsplit()):
      runningSelect = band > imBandMax
      imBandMax = runningSelect.ifthenelse(band, imBandMax)
      imBandSelect = runningSelect.ifthenelse(bandIndex, imBandSelect)

    bandList = [ (imBandSelect == bi) / 255.0 for bi in range(inImg.bands) ]
    return bandList[0].bandjoin(bandList[1:])

Update:

Thanks to @jcupitt I tried this version of the code using bandrank:

  def pickMaxUchar(inImg):
    short = (inImg.cast(pyvips.enums.BandFormat.USHORT)) << 8
    index = pyvips.Image.black(short.width, short.height, bands=1).bandjoin_const([b for b in range(1, inImg.bands)]).cast(pyvips.enums.BandFormat.UCHAR)
    combo = short | index
    list = combo.bandsplit()
    ranked = list[0].bandrank(list[1:], index=inImg.bands-1)
    rankedindex = (ranked & 255).cast(pyvips.enums.BandFormat.UCHAR)
    bandList = [ (rankedindex == bi) / 255.0 for bi in range(inImg.bands) ]
    return bandList[0].bandjoin(bandList[1:])

This now assumes the input is a char, where in my original function the input was float. Maybe there's a way to do this without re-casting the data but as I've implemented it there's zero speedup relative to my "naive" code above, so perhaps this just can't be accelerated any further.

Answer 1

You can use bandrank for this. It's like a rank filter, so just ask for the final index.

It returns a one-band image with the selected value in each pixel and you need the index, so the trick is to hide the index in the low bits of the pixel values, and pull it out again from the result.

Something like:

#!/usr/bin/env python

import sys
import pyvips

image = pyvips.Image.new_from_file(sys.argv[1], access='sequential')

# a band index image, ie. each band contains its own index
index = (pyvips.Image.black(image.width, image.height) + [0, 1, 2]).cast("uchar")

# put that into the bottom two bits of the image
image = (image << 2) | index

# now use bandrank to find the max at each pixel
# index=2 means the max of a three band image
bands = image.bandsplit()
mx = bands[0].bandrank(bands[1:], index=2)

# and extract the index value
index = mx & 3

index.write_to_file(sys.argv[2])

With a 30,000 x 30,000 pixel jpg I see:

$ time ./rank.py ~/pics/st-francis.jpg x.v

real    0m5.296s
user    0m51.233s
sys 0m1.501s