How to perform logical operation and logical indexing using VIPS in Python?

Question

I've had following codes that use Python and OpenCV. Briefly, I have a stack of image taken at different focal depth. The codes pick out pixels at every (x,y) position that has the largest Laplacian of Guassian response among all focal depth(z), thus creating a focus-stacked image. Function get_fmap creates a 2d array where each pixel will contains the number of the focal plane having the largest log response. In the following codes, lines that are commented out are my current VIPS implementation. They don't look compatible within the function definition because it's only partial solution.

# from gi.repository import Vips

def get_log_kernel(siz, std):
    x = y = np.linspace(-siz, siz, 2*siz+1)
    x, y = np.meshgrid(x, y)
    arg = -(x**2 + y**2) / (2*std**2)
    h = np.exp(arg)
    h[h < sys.float_info.epsilon * h.max()] = 0
    h = h/h.sum() if h.sum() != 0 else h
    h1 = h*(x**2 + y**2 - 2*std**2) / (std**4)
    return h1 - h1.mean()

def get_fmap(img):    # img is a 3-d numpy array.
    log_response = np.zeros_like(img[:, :, 0], dtype='single')
    fmap = np.zeros_like(img[:, :, 0], dtype='uint8')
    log_kernel = get_log_kernel(11, 2)
    # kernel = get_log_kernel(11, 2)
    # kernel = [list(row) for row in kernel]
    # kernel = Vips.Image.new_from_array(kernel)
    # img = Vips.new_from_file("testimg.tif")
    for ii in range(img.shape[2]):           
        # img_filtered = img.conv(kernel)
        img_filtered = cv2.filter2D(img[:, :, ii].astype('single'), -1, log_kernel)
        index = img_filtered > log_response
        log_response[index] = img_filtered[index]
        fmap[index] = ii
    return fmap

and then fmap will be used to pick out pixels from different focal planes to create a focus-stacked image

This is done on an extremely large image, and I feel VIPS might do a better job than OpenCV on this. However, the official documentation provides rather scant information on its Python binding. From the information I can find on the internet, I'm only able to make image convolution work ( which, in my case, is an order of magnitude faster than OpenCV.). I'm wondering how to implement this in VIPS, especially these lines?

log_response = np.zeros_like(img[:, :, 0], dtype = 'single')

index = img_filtered > log_response

log_response[index] = im_filtered[index]

fmap[index] = ii

Answer 1

log_response and fmap are initialized as 3D arrays in the question code, whereas the question text states that the output, fmap is a 2D array. So, I am assuming that log_response and fmap are to be initialized as 2D arrays with their shapes same as each image. Thus, the edits would be -

log_response = np.zeros_like(img[:,:,0], dtype='single')
fmap = np.zeros_like(img[:,:,0], dtype='uint8')

Now, back to the theme of the question, you are performing 2D filtering on each image one-by-one and getting the maximum index of filtered output across all stacked images. In case, you didn't know as per the documentation of cv2.filter2D, it could also be used on a multi-dimensional array giving us a multi-dimensional array as output. Then, getting the maximum index across all images is as simple as .argmax(2). Thus, the implementation must be extremely efficient and would be simply -

fmap = cv2.filter2D(img,-1,log_kernel).argmax(2)

Answer 2

After consulting the Python VIPS manual and some trial-and-error, I've come up with my own answer. My numpy and OpenCV implementation in question can be translated into VIPS like this:

import pyvips

img = []
for ii in range(num_z_levels):
    img.append(pyvips.Image.new_from_file("testimg_z" + str(ii) + ".tif")

def get_fmap(img)
    log_kernel = get_log_kernel(11,2)  # get_log_kernel is my own function, which generates a 2-d numpy array.
    log_kernel = [list(row) for row in log_kernel]  # pyvips.Image.new_from_array takes 1-d list array.
    log_kernel = pyvips.Image.new_from_array(log_kernel)  # Turn the kernel into Vips array so it can be used by Vips.
    log_response = img[0].conv(log_kernel)

    for ii in range(len(img)):
        img_filtered = img[ii+1].conv(log_kernel)
        log_response = (img_filtered > log_response).ifthenelse(img_filtered, log_response)
        fmap = (img_filtered > log_response).ifthenelse(ii+1, 0)

Logical indexing is achieved through ifthenelse method :

result_img = (test_condition).ifthenelse(value_if_true, value_if_false)

The syntax is rather flexible. The test condition can be a comparison between two images of the same size or between an image and a value, e.g. img1 > img2 or img > 5. Like wise, value_if_true can be a single value or a Vips image.