Image Operations

The Hypervision SDK provides a number of built-in operations that work as building blocks for more complex behaviour.

Slicing

The hsi.HSImage.slice() operation (or simply the index operator [] in Python), can be used to construct a pipeline that only performs computations on the selected data.

How slicing works

Slicing is special, because it propagates back through operations as far as possible, which means that it can limit the amount of data being operated on after the previous operations have been constructed. For example, slicing a pipeline that loads data from a file will limit the amount of data being loaded from disk.

This feature is used internally in the library to implement the functionality of methods like hsi.HSImage.array_plane() .

A few examples:

img[:, :, :] # The full image
img[:20, :20, :] # The first 20 lines and samples
img[:, :, ::4] # Take only every 4th band

# Demonstration of how lines/samples/bands map to the index dimensions
lines = slice(start=0, end=10)
samples = slice(end=None)
bands = slice(end=None, step=2)
img[lines, samples, bands]

// The slice descriptor uses BIP order: (lines, samples, bands)
hv_slice_desc_t desc_full = {
    .lines = hv_slice_elem_new_range(0, -1, 1),
    .samples = hv_slice_elem_new_range(0, -1, 1),
    .bands = hv_slice_elem_new_range(0, -1, 1),
};
hv_hs_image_t *full = hv_hs_image_slice(img, desc_full);

// First 20 lines and samples, all bands
hv_slice_desc_t desc_roi = {
    .lines = hv_slice_elem_new_range(0, 20, 1),
    .samples = hv_slice_elem_new_range(0, 20, 1),
    .bands = hv_slice_elem_new_range(0, -1, 1),
};
hv_hs_image_t *roi = hv_hs_image_slice(img, desc_roi);

// Take only every 4th band
hv_slice_desc_t desc_stride = {
    .lines = hv_slice_elem_new_range(0, -1, 1),
    .samples = hv_slice_elem_new_range(0, -1, 1),
    .bands = hv_slice_elem_new_range(0, -1, 4),
};
hv_hs_image_t *sliced_img = hv_hs_image_slice(img, desc_stride);

Note for Python users

The slicing syntax is a bit more restrictive than NumPy. You always need exactly three arguments. The reason for this is that a hsi.HSImage is always three-dimensional.

The slice indices are always in BIP format (lines, samples, bands) as shown in the example above. This can cause confusion when used together with the hsi.HSImage.to_numpy() method.

img = hsi.HSImage.from_numpy(np.ones((10, 10, 10)), interleave=hsi.bil)

# Note how the index order differs between the two types. This is because `img` is in BIL format but the index
# method always accepts BIP format.
assert img[:, :, :5].to_numpy().shape == (10, 5, 10)

Elementwise operations

You can use common arithmetic elementwise operations using the regular operators +, -, *, /. Currently, only these four are supported. The operations work on both hsi.HSImage objects and scalars (with both left and right scalar operations).

# Subtract two bands
diff = img[:, :, 500] - img[:, :, 810]

# Scale image from Hypervision 1700 to (0,1) (image is assumed to be float)
scaled = img / 254.

// Create slices
hv_slice_desc_t b500 = {
    .lines = hv_slice_elem_new_range(0, -1, 1),
    .samples = hv_slice_elem_new_range(0, -1, 1),
    .bands = hv_slice_elem_new_index(500),
};
hv_slice_desc_t b810 = {
    .lines = hv_slice_elem_new_range(0, -1, 1),
    .samples = hv_slice_elem_new_range(0, -1, 1),
    .bands = hv_slice_elem_new_index(810),
};

// Slice image into two bands and subtract them
hv_hs_image_t *band500 = hv_hs_image_slice(img, b500);
hv_hs_image_t *band810 = hv_hs_image_slice(img, b810);
hv_hs_image_t *diff = hv_hs_image_sub(band500, band810);

// Scale the image to (0,1) (image is assumed to be float32_t)
hv_hs_image_t *scaled = hv_hs_image_div_scalar_left_f32(img, 254.0);

Important

The scalar operations do not automatically convert the type. Instead, the output type is always determined by the hsi.HSImage operand.

You can convert the data types using hsi.HSImage.ensure_dtype() .

Reduction operations

The library supports a limited set of reduction operations for calculating the mean, standard deviation, variance, and sum. The operations always work on a single axis only. Because the images are always three-dimensional, the operations reduce the specified axis to length 1, similar to how the keepdims argument works in NumPy.

img.mean_axis(hsi.bands)  # Mean spatial image

# Workaround for performing multiple reductions. The result is the mean spectrum.
img.mean_axis(hsi.lines).mean_axis(hsi.samples)

// Mean over bands (spatial mean image)
hv_hs_image_t *mean_spatial = hv_hs_image_mean_axis(img, HV_AXIS_BANDS);

// Chained reductions to get mean spectrum
hv_hs_image_t *mean_over_lines = hv_hs_image_mean_axis(mean_spatial, HV_AXIS_LINES);
hv_hs_image_t *mean_spectrum = hv_hs_image_mean_axis(mean_over_lines, HV_AXIS_SAMPLES);

Matrix multiplication

Matrix multiplication is essential for many different kinds of processes, including projections and certain machine learning models. The Hypervision SDK provides a generic method hsi.HSImage.dot() which implements both a dot product and vector-matrix multiplication along a specified hsi.Axis .

Let us start by looking at how to use the hsi.HSImage.dot() method in practice. The dropdown below provides additional details on the implementation of the method.

Matrix multiplication calculation details

Let us look at how the calculations are performed. Given a hyperspectral image \(H \in \mathbb{R}^{L\times S\times B}\) and a vector operand \(\mathbf{u} \in \mathbb{R}^{B}\), the dot product is defined as:

\[\begin{split}H_{l,s,1} &= \mathbf{u} \begin{bmatrix}H_{l,s,1}\\\vdots\\ H_{l,s,B}\end{bmatrix}\end{split}\]

For the other two axes, let \(\mathbf{v} \in \mathbb{R}^{L}\) and \(\mathbf{w} \in \mathbb{R}^{S}\), then:

\[\begin{split}H_{1,s,b} &= \mathbf{v} \begin{bmatrix}H_{1,s,b}\\\vdots\\ H_{L,s,b}\end{bmatrix}\\ H_{l,1,b} &= \mathbf{w} \begin{bmatrix}H_{l,1,b}\\\vdots\\ H_{l,S,b}\end{bmatrix}\end{split}\]

For vector-matrix multiplication, let \(U \in \mathbb{R}^{N\times B}\) be the operand. The operation is then defined as:

\[\begin{split}\begin{bmatrix}H_{1,s,b}\\\vdots\\ H_{N,s,b}\end{bmatrix} = U \begin{bmatrix}H_{1,s,b}\\\vdots\\ H_{L,s,b}\end{bmatrix}\end{split}\]

For the other two axes, let \(V \in \mathbb{R}^{N\times L}\) and \(W \in \mathbb{R}^{N\times S}\), then:

\[\begin{split}\begin{bmatrix}H_{l,1,b}\\\vdots\\ H_{l,N,b}\end{bmatrix} &= V \begin{bmatrix}H_{l,1,b}\\\vdots\\ H_{l,S,b}\end{bmatrix} \\ \begin{bmatrix}H_{l,s,1}\\\vdots\\ H_{l,s,N}\end{bmatrix} &= W \begin{bmatrix}H_{l,s,1}\\\vdots\\ H_{1,s,B}\end{bmatrix}\end{split}\]

Note

The implementation uses matrix multiplication as much as possible for efficiency. For example, A vector-matrix multiply along the bands axis in an BIP image is implemented as a matrix-matrix multiplication performed on each image plane, with each plane consisting of bands as its rows and samples as its columns.

# Create a new vector with element 2 and 200 set to 0.5
operand = np.zeros(920)
operand[2] = 0.5
operand[200] = 0.5

# Calculate the dot product. The result is a mix between band 2 and 200.
res = img.dot(operand, hsi.Bands)

# Create a new matrix operand
operand = np.zeros((2, 920))
operand[0, 2] = 0.5
operand[0, 200] = 0.5
operand[1, 25] = 1

# Calculate the matrix-vector product. The result is an image with 2 output bands.
# The first band is equal to the dot-product above. The second band is simply band 25 from the input
# image.
res = img.dot(operand, hsi.bands)

// Create a new vector with element 2 and 200 set to 0.5
float32 operand_array[920];
operand_array[2] = 0.5;
operand_array[200] = 0.5;
slice_ref_float32_t operand_slice = {
    .len = 920,
    .ptr = &operand_array
};

slice_ref_size_t shape = {
    .len = 1,
    .ptr = (size_t[1]){920},
}

hv_array_t operand;
hv_array_from_f32(shape, operand_slice, &operand);

// Calculate the dot product. The result is a mix between band 2 and 200.
hv_hs_image_t* res = hv_hs_image_dot(img, operand, HV_AXIS_BANDS);

// Create a new matrix operand
float32 operand_array2[2*920];
operand_array2[2] = 0.5;
operand_array2[200] = 0.5;
operand_array2[1*920 + 25] = 0.5;

slice_ref_float32_t operand_slice2 = {
    .len = 920*2,
    .ptr = &operand_array2
};

slice_ref_size_t shape2 = {
    .len = 2,
    .ptr = (size_t[1]){2, 920},
}

hv_array_t operand2;
hv_array_from_f32(shape2, operand_slice2, &operand2);

// Calculate the matrix-vector product. The result is an image with 2 output bands.
// The first band is equal to the dot-product above. The second band is simply band 25 from the input
// image.
hv_hs_image_t* res2 = hv_hs_image_dot(img, operand2, HV_AXIS_BANDS);

Other operations

The SDK provides a number of useful operations that

# Apply mean-binning (beware that the dtype will not be changed automatically)
# Use `ensure_dtype()` first if the input type is not big enough to hold the binning results
condensed = img.binning(8, hsi.bands)

# Clip image to the provided range
clipped = img.clip(0, 1)

# Convert NaN values to the provided value
non_nan = img.nan_to_num(0)

# Type conversions
new_type = img.as_dtype(hsi.float32) # This one always does the convertion, even if the type is already the same
ensured_type = img.ensure_dtype(hsi.float32) # This one only converts if the types are different

// Apply mean-binning (beware that the dtype will not be changed automatically)
// Use `ensure_dtype()` first if the input type is not big enough to hold the binning results
hv_hs_image_t* condensed = hv_hs_image_binning(img, 8, HV_AXIS_BANDS);

// Clip image to the provided range
hv_hs_image_t* clipped = hv_hs_image_clip(img, 0.0f, 1.0f);

// Convert NaN values to the provided value
hv_hs_image_t* non_nan = hv_hs_image_nan_to_num(img, 0.0);

// Type conversions

// This one always does the convertion, even if the type is already the same
hv_hs_image_t* new_type = hv_hs_image_as_dtype(img, HV_DTYPE_F32);

// This one only converts if the types are different
hv_hs_image_t* ensured_type;
hv_hs_image_ensure_dtype(img, HV_DTYPE_F32, &ensured_type);