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.
1 Answer
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 4