Custom operations

The Hypervision SDK provides facilities for running custom code inside its pipeline system. This means that you can take advantage of the same lazy calculations and efficient streaming, that the built-in operations use.

In order to ease the use of these custom operations, the SDK provides a number of helper functions, including @qtec_hv_sdk.util.operation and qtec_hv_sdk.util.predictor().

Important

The helper functions are currently only available in the Python extension. They may be added to the C library later if there is a need for high-level integrations.

A custom operation can be created by using the qtec_hv_sdk.Image.ufunc() method like so:

Warning

The C API is very verbose due to closures not being natively supported. The API will likely change in the future to something more suited to use from C.

def my_function(meta: hs.FrameMeta, plane: np.ndarray) -> np.ndarray:
    print(meta.seq)
    return plane.mean(axis=1, keepdims=True)

# Apply `my_function` to each plane of the image.
res = img.ufunc(my_function)

# The resulting image uses the same lazy pipeline as the built-in operations.
res.array_plane(200, hs.bands)

# When resolving the entire image, the function is applied on a plane-by-plane basis.
out = res.to_numpy()

#include <assert.h>
#include "hv-sdk.h"
#include <stdio.h>

typedef struct {
    int ref_count;
    int count;

    hv_array_t *array;
} callback_env;

hv_array_t *callback(callback_env *env_ptr, hv_array_t *plane) {
    printf("Plane no: %d\n", env_ptr->count++);
    return plane;
}

void release_fn(callback_env *env_ptr) {
    printf("release %d", env_ptr->ref_count);
    env_ptr->ref_count--;
    if (env_ptr->ref_count == 0) {
        hv_array_free(&env_ptr->array);
    }
}

void retain_fn(callback_env *env_ptr) {
    env_ptr->ref_count++;
    printf("retain %d", env_ptr->ref_count);
}

int main() {
    // Loads PAM file due to the extension
    const char *path = "resources/docs/ex1.pam";
    hv_hsi_file_t *file;
    assert(!hv_hsi_file_open(path, &file));

    // Convert the file to an image (transfers ownership - file is NULL after this)
    hv_image_t *img = hv_hsi_file_to_image(&file);

    hv_shape_meta_t *shape_meta;
    assert(!hv_image_shape(img, &shape_meta));

    slice_ref_size_t sh = {
        .len = 3,
        .ptr = (size_t[3]){
            *hv_shape_meta_lines(shape_meta),
            *hv_shape_meta_samples(shape_meta),
            *hv_shape_meta_bands(shape_meta)
        }
    };

    callback_env env = {
        .ref_count = 0,
        .array = hv_array_zeros(HV_DTYPE_U8, sh)
    };
    ArcDynFn1_hv_array_ptr_hv_array_ptr_t closure = {
        .env_ptr = &env,
        .retain = (void (*)(void *)) &retain_fn,
        .release = (void (*)(void *)) &release_fn,
        .call = (hv_array_t* (*)(void *, hv_array_t *)) &callback
    };

    hv_image_t *res = hv_image_op(img, closure);

    hv_array_t *plane;
    hv_image_array_plane(res, 10, HV_AXIS_BANDS, &plane);

    hv_image_free(&res);
    hv_array_free(&plane);
}

This method does not account for different interleaves and always operates on the inner two dimensions of the image.

General operations

The SDK provides a decorator qtec_hv_sdk.util.operation(), that ensures a desired interleave and wraps a function so it can be used directly with qtec_hv_sdk.Image:

from qtec_hv_sdk.util import operation

@operation(hs.bip)
def my_function(meta: hs.FrameMeta, plane: np.ndarray) -> np.ndarray:
    return plane.mean(axis=1, keepdims=True)

# The call now ensures that the interleave is consistent.
res = my_function(img)

Advanced operations

In many cases, qtec_hv_sdk.Image.ufunc() is able to correctly infer how to propagate slicing to earlier operations. In the example above, it will correctly infer that the bands dimension will be flattened and that, even if we only request a single band, e.g. with res.array_plane(0, hs.bands), the operation needs to request all bands from its source.

But this does not always work! A great example is a per-plane SNV normalization operation. Let us investigate by looking at a naive implementation:

@operation(hs.bip)
def snv(meta: hs.FrameMeta, plane: np.ndarray) -> np.ndarray:
    mean = plane.mean(axis=1, keepdims=True)
    std = plane.mean(axis=1, keepdims=True)

    return (plane - mean) / std

res = snv(img)

If we call res.to_numpy(), the operation works as intended, but if we call res.array_plane(500, hs.bands), the operation only gets a Nx1 plane as its input. The reason is that the output is no longer flat but depends on the slicing provided and as a result, the slice ratio inference breaks down.

The solution is to provide a slice_transform argument to the decorator that maps output slices to input slices. This may seem a bit abstract, but the solution for SNV should make it clearer:

import numpy as np
import qtec_hv_sdk as hs
from qtec_hv_sdk.util import operation

# The transformation function has to map a slice received from a downstream operation to a slice to
# send to the upstream (parent) operation.
#
# In this case we use BIP, which means that the slice has the form (samples, bands).
# The transform then simply replaces the bands part with None to always get all bands from the parent.
def full_bands(plane_slice: tuple[slice, slice]) -> tuple[slice, slice]:
    samples, _bands = plane_slice
    return (samples, slice(None))

# The only change here is to add the slice_transform argument
@operation(hs.bip, slice_transform=full_bands)
def snv(meta: hs.FrameMeta, plane: np.ndarray) -> np.ndarray:
    mean = plane.mean(axis=1, keepdims=True)
    std = plane.mean(axis=1, keepdims=True)

    return (plane - mean) / std

res = snv(img)

# Examples

# Get single band.
# The output slice is [:, 500] (we skip the lines since the operation works on one line at a time)
# The mapping function turns this slice into [:, :] and thus gets all bands so SNV works correctly.
a = res[:, :, 500].to_numpy()

# Arbitrary crop
# The output slice is [100:150, 200:300] and maps to input [100:150, :]
b = res[250, 100:150, 200:300]

# You can also test the mapping function directly to get a feel for how it works
assert full_bands((slice(100, 150), slice(200, 300))) == (slice(100, 150), slice(None))

Machine learning models

The SDK also provides a helper function qtec_hv_sdk.util.predictor() for adapting certain scikit-learn like ML models to qtec_hv_sdk.Image s. It support models that use individual spatial pixels as input.

from qtec_hv_sdk.util import predictor
from sklearn.linear_model import LinearRegression

img = hs.open(<path>)

model = LinearRegression()
model.fit(...) # Training step

hsimage_predictor = predictor(model)

res = hsimage_predictor(img)

# Again, the operation is lazy and is only applied when necessary.
out = res.to_numpy()

# This means that operations like slicing will reduce the number of computations.
out_small = res[50:60, :, :].to_numpy()