Partial bugfix for TIRLFile object ID redundancy problem.

src/tirl/cmodules/finv2d_src.h

+/*
+   _______   _____   _____    _
+  |__   __| |_   _| |  __ \  | |
+     | |      | |   | |__) | | |
+     | |      | |   |  _  /  | |
+     | |     _| |_  | | \ \  | |____
+     |_|    |_____| |_|  \_\ |______|
+
+ Copyright (C) 2018-2020 University of Oxford
+ Part of the FMRIB Software Library (FSL)
+ Author: Istvan N. Huszar
+ Copyright (C) 2018-2020 University of Oxford */
+
+/*  CCOPYRIGHT  */
+
+
+#ifndef FINV_FINV2D_H
+#define FINV_FINV2D_H
+
+#include <lapacke.h>
+
+// Basic math
+
+double Dot2D(double v[], double u[]);
+int Difference2D(double v[], double u[], double* Result);
+int DifferenceInt2D(int v[], int u[], int* Result);
+int Unravel2DIndex(int Index, int* Shape, int* MultiIndex);
+int Ravel2DIndex(int* MultiIndex, int* Shape);
+lapack_int InvertMat2D(double* A, int N);
+
+// 2D functions
+
+int IsInTriangle(double* Point, double** Simplex);
+
+int StepTriangle(int** Simplex, int Vertex, int* Shape);
+
+int InitialiseTriangle(int ** Simplex, int* Shape);
+
+int CalculateInverse2D(double** Locations, int* NearestNeighbour, int N,
+    int* InputShape, double** Coordinates, double** FwdMappedCoordinates,
+    double** InvWarp, int Iterations);
+
+int CalculateLocalAffine2D(double** Locations, int* NearestNeighbour,
+    int N, int* InputShape, double** Coordinates,
+    double** FwdMappedCoordinates, double** Jacobian, int Iterations);
+
+int CalculateTxJacobian2D(double** Locations, int* NearestNeighbour,
+    int N, int* InputShape, double** Coordinates,
+    double** Jacobian, int** VertexIndex, int Iterations);
+
+int FillHoles2D(double** InvWarp, int* Shape, int MaxIter);
+
+
+#endif //FINV_FINV2D_H

src/tirl/cmodules/finv3d_src.h

+/*
+   _______   _____   _____    _
+  |__   __| |_   _| |  __ \  | |
+     | |      | |   | |__) | | |
+     | |      | |   |  _  /  | |
+     | |     _| |_  | | \ \  | |____
+     |_|    |_____| |_|  \_\ |______|
+
+ Copyright (C) 2018-2020 University of Oxford
+ Part of the FMRIB Software Library (FSL)
+ Author: Istvan N. Huszar
+ Copyright (C) 2018-2020 University of Oxford */
+
+/*  CCOPYRIGHT  */
+
+
+#ifndef FINV_FINV3D_H
+#define FINV_FINV3D_H
+
+#include <lapacke.h>
+
+// Basic math
+
+double Dot3D(double v[], double u[]);
+int Cross(double v[], double u[], double* Result);
+int Difference3D(double v[], double u[], double* Result);
+int DifferenceInt3D(int v[], int u[], int* Result);
+int Unravel3DIndex(int Index, int* Shape, int* MultiIndex);
+int _lin(int* Shape, int i, int j, int k);
+int Ravel3DIndex(int* MultiIndex, int* Shape);
+lapack_int InvertMat3D(double* A, int N);
+
+// 3D functions
+
+int IsInTetrahedron(double* Point, double** Simplex);
+
+int StepTetrahedron(int** Simplex, int Vertex, int* Shape);
+
+int InitialiseTetrahedron(int ** Simplex, int* Shape);
+
+int CalculateInverse3D(double** Locations, int* NearestNeighbour, int N,
+    int* InputShape, double** Coordinates, double** FwdMappedCoordinates,
+    double** InvWarp, int Iterations);
+
+int CalculateLocalAffine3D(double** Locations,
+    int* NearestNeighbour, int N, int* InputShape,
+    double** Coordinates, double** FwdMappedCoordinates,
+    double** LocalAffine, int Iterations);
+
+int CalculateTxJacobian3D(double** Locations, int* NearestNeighbour,
+    int N, int* InputShape, double** Coordinates,
+    double** Jacobian, int** VertexIndex, int Iterations);
+
+
+
+int FillHoles3D(double** InvWarp, int* Shape, int MaxIter);
+
+
+#endif //FINV_FINV3D_H

src/tirl/cmodules/fslinterpolator.pyx

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+#   _______   _____   _____    _
+#  |__   __| |_   _| |  __ \  | |
+#     | |      | |   | |__) | | |
+#     | |      | |   |  _  /  | |
+#     | |     _| |_  | | \ \  | |____
+#     |_|    |_____| |_|  \_\ |______|
+#
+# Copyright (C) 2018-2020 University of Oxford
+# Part of the FMRIB Software Library (FSL)
+# Author: Istvan N. Huszar
+
+
+# SHBASECOPYRIGHT
+
+
+# DEVELOPMENT NOTES
+
+# PyCharm doesn't fully understand the .pyx syntax, therefore marks
+# certain variables as undefined in this script.
+
+
+# DEPENDENCIES
+
+import numpy as np
+
+
+# IMPLEMENTATION
+
+# Load the C++ function that performs the interpolation over the entire range
+# of coordinates. Note that F is the template argument (typename).
+cdef extern from "fslinterpolator_src.hpp" namespace "FSLINTERPOLATOR":
+    void interpolate_all[F](F * data, size_t * shape, float * coordinates,
+                            size_t n, F * out);
+
+# Compound type definitions are not always syntactically valid, so we assign
+# new names for these.
+ctypedef unsigned char ubyte
+ctypedef unsigned short ushort
+ctypedef unsigned int uint
+ctypedef unsigned long ulong
+ctypedef long double ldouble
+
+
+def fslinterpolator(data, coordinates):
+    """
+    Forefront function that chooses the appropriate implementation of the
+    interpolation algorithm based on the input data type.
+
+    The FSLInterpolator works with 32-bit floating-point coordinates and
+    interpolates data of the following types:
+
+    Integer types:
+        int8, int16, int32, int64
+        uint8, uint16, uint32, uint64
+    Floating-point types:
+        float32, float64, float128
+    Complex types:
+        complex64, complex128, complex256
+
+    The data type is preserved between the input and the output.
+
+    :param data: input data array with known values (any shape)
+    :type data: np.ndarray
+    :param coordinates:
+        (n_points, n_dimensions) table of the sampling coordinates. The number
+        of dimensions must match the number of dimensions of the input array.
+    :type coordinates: np.ndarray
+
+    :returns: (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    coordinates = np.asarray(coordinates, dtype=np.float32)
+    dtype = np.dtype(data.dtype).name
+    function_bindings = {
+        "int8": fslintp_int8,
+        "int16": fslintp_int16,
+        "int32": fslintp_int32,
+        "int64": fslintp_int64,
+        "uint8": fslintp_uint8,
+        "uint16": fslintp_uint16,
+        "uint32": fslintp_uint32,
+        "uint64": fslintp_uint64,
+        "float32": fslintp_float32,
+        "float64": fslintp_float64,
+        "float128": fslintp_float128,
+        "complex64": fslintp_complex64,
+        "complex128": fslintp_complex128,
+        "complex256": fslintp_complex256
+    }
+    func = function_bindings.get(dtype, None)
+    if func is None:
+        raise TypeError(f"FSLInterpolator cannot handle "
+                        f"the following data type: {dtype}.")
+    return func(data, coordinates)
+
+
+# The functions below implement the Python wrapper of the C++ interpolation
+# algorithm for various data types. Note that these cannot be shortened by a
+# conditional iteration, because array definitions must be explicitly declared
+# to be recognised by the compiler.
+
+def fslintp_float32(data, coordinates):
+    """
+    FSLInterpolator: Python binding for 32-bit floating-point input data.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef float[:, :, ::1] d = np.ascontiguousarray(data, np.float32)
+    cdef float[::1] out = np.zeros(n, dtype=np.float32, order="C")
+
+    interpolate_all[float](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_float64(data, coordinates):
+    """
+    FSLInterpolator: Python binding for 64-bit floating-point input data.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef double[:, :, ::1] d = np.ascontiguousarray(data, np.float64)
+    cdef double[::1] out = np.zeros(n, dtype=np.float64, order="C")
+
+    interpolate_all[double](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_float128(data, coordinates):
+    """
+    FSLInterpolator: Python binding for 128-bit floating-point input data.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef long double[:, :, ::1] d = np.ascontiguousarray(data, np.float128)
+    cdef long double[::1] out = np.zeros(n, dtype=np.float128, order="C")
+
+    interpolate_all[ldouble](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_uint8(data, coordinates):
+    """
+    FSLInterpolator: Python binding for unsigned 8-bit integer input data. This
+    is the most common 2D image data type (per colour channel). As an unsigned
+    type, it is vulnerable to overflow on subtraction (negative values roll
+    over to become excessively large positive numbers).
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef ubyte[:, :, ::1] d = np.ascontiguousarray(data, np.uint8)
+    cdef ubyte[::1] out = np.zeros(n, dtype=np.uint8, order="C")
+
+    interpolate_all[ubyte](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_uint16(data, coordinates):
+    """
+    FSLInterpolator: Python binding for unsigned 16-bit integer input data.
+    This format is commonly used by SPM (Statistical Parametric Mapping, UCL).
+    As an unsigned type, it usually represents renormalised data, and is
+    vulnerable to overflow on subtraction (negative values roll over to become
+    excessively large positive numbers).
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef unsigned short[:, :, ::1] d = np.ascontiguousarray(data, np.uint16)
+    cdef unsigned short[::1] out = np.zeros(n, dtype=np.uint16, order="C")
+
+    interpolate_all[ushort](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_uint32(data, coordinates):
+    """
+    FSLInterpolator: Python binding for unsigned 32-bit integer input data.
+    As an unsigned type, it usually represents renormalised data, and is
+    vulnerable to overflow on subtraction (negative values roll over to become
+    excessively large positive numbers).
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef uint[:, :, ::1] d = np.ascontiguousarray(data, np.uint32)
+    cdef uint[::1] out = np.zeros(n, dtype=np.uint32, order="C")
+
+    interpolate_all[uint](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_uint64(data, coordinates):
+    """
+    FSLInterpolator: Python binding for unsigned 64-bit integer input data.
+    As an unsigned type, it usually represents renormalised data, and is
+    vulnerable to overflow on subtraction (negative values roll over to become
+    excessively large positive numbers).
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef ulong[:, :, ::1] d = np.ascontiguousarray(data, np.uint64)
+    cdef ulong[::1] out = np.zeros(n, dtype=np.uint64, order="C")
+
+    interpolate_all[ulong](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_int8(data, coordinates):
+    """
+    FSLInterpolator: Python binding for signed 8-bit integer input data.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef char[:, :, ::1] d = np.ascontiguousarray(data, np.int8)
+    cdef char[::1] out = np.zeros(n, dtype=np.int8, order="C")
+
+    interpolate_all[char](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_int16(data, coordinates):
+    """
+    FSLInterpolator: Python binding for signed 16-bit integer input data.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef short[:, :, ::1] d = np.ascontiguousarray(data, np.int16)
+    cdef short[::1] out = np.zeros(n, dtype=np.int16, order="C")
+
+    interpolate_all[short](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_int32(data, coordinates):
+    """
+    FSLInterpolator: Python binding for signed 32-bit integer input data.
+    This is one of the most common integer data types.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef int[:, :, ::1] d = np.ascontiguousarray(data, np.int32)
+    cdef int[::1] out = np.zeros(n, dtype=np.int32, order="C")
+
+    interpolate_all[int](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_int64(data, coordinates):
+    """
+    FSLInterpolator: Python binding for signed 64-bit integer input data.
+    This is one of the most common integer data types.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef long[:, :, ::1] d = np.ascontiguousarray(data, np.int64)
+    cdef long[::1] out = np.zeros(n, dtype=np.int64, order="C")
+
+    interpolate_all[long](&d[0, 0, 0], &shape[0], &c[0, 0], n, &out[0])
+    return out
+
+
+def fslintp_complex64(data, coordinates):
+    """
+    FSLInterpolator: Python binding for 64-bit complex input data.
+    This is the default complex data type in the Siemens raw-format MRI data.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef float[:, :, ::1] d_real = np.ascontiguousarray(data.real, np.float32)
+    cdef float[:, :, ::1] d_imag = np.ascontiguousarray(data.imag, np.float32)
+    cdef float[::1] o_real = np.zeros(n, dtype=np.float32, order="C")
+    cdef float[::1] o_imag = np.zeros(n, dtype=np.float32, order="C")
+
+    # Real part
+    interpolate_all[float](&d_real[0, 0, 0], &shape[0], &c[0, 0], n,
+                           &o_real[0])
+    # Imaginary part
+    interpolate_all[float](&d_imag[0, 0, 0], &shape[0], &c[0, 0], n,
+                           &o_imag[0])
+    # Combine and return
+    out = np.asarray(o_real) + 1j * np.asarray(o_imag)
+    return out.astype(np.complex64)
+
+
+def fslintp_complex128(data, coordinates):
+    """
+    FSLInterpolator: Python binding for 128-bit complex input data. This data
+    type provides double precision in both the real and the imaginary
+    components, therefore it is most suitable for manipulating complex data.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef float[:, :, ::1] d_real = np.ascontiguousarray(data.real, np.float64)
+    cdef float[:, :, ::1] d_imag = np.ascontiguousarray(data.imag, np.float64)
+    cdef float[::1] o_real = np.zeros(n, dtype=np.float64, order="C")
+    cdef float[::1] o_imag = np.zeros(n, dtype=np.float64, order="C")
+
+    # Real part
+    interpolate_all[float](&d_real[0, 0, 0], &shape[0], &c[0, 0], n,
+                           &o_real[0])
+    # Imaginary part
+    interpolate_all[float](&d_imag[0, 0, 0], &shape[0], &c[0, 0], n,
+                           &o_imag[0])
+    # Combine and return
+    out = np.asarray(o_real) + 1j * np.asarray(o_imag)
+    return out.astype(np.complex128)
+
+
+def fslintp_complex256(data, coordinates):
+    """
+    FSLInterpolator: Python binding for 256-bit complex input data. This is
+    normally too large for practical applications.
+
+    :param data: array of input values
+    :type data: np.ndarray
+    :param coordinates: (n_points, n_dimensions) table of sampling coordinates
+    :type coordinates: np.ndarray
+
+    :returns (n_points,) interpolated values
+    :rtype: np.ndarray
+
+    """
+    cdef int n = coordinates.shape[0]
+    cdef size_t [::1] shape = np.asarray(data.shape, dtype=np.uintp)
+    cdef float[:, ::1] c = np.ascontiguousarray(coordinates, np.float32)
+    cdef float[:, :, ::1] d_real = np.ascontiguousarray(data.real, np.float128)
+    cdef float[:, :, ::1] d_imag = np.ascontiguousarray(data.imag, np.float128)
+    cdef float[::1] o_real = np.zeros(n, dtype=np.float128, order="C")
+    cdef float[::1] o_imag = np.zeros(n, dtype=np.float128, order="C")
+
+    # Real part
+    interpolate_all[float](&d_real[0, 0, 0], &shape[0], &c[0, 0], n,
+                           &o_real[0])
+    # Imaginary part
+    interpolate_all[float](&d_imag[0, 0, 0], &shape[0], &c[0, 0], n,
+                           &o_imag[0])
+    # Combine and return
+    out = np.asarray(o_real) + 1j * np.asarray(o_imag)
+    return out.astype(np.complex256)

src/tirl/cmodules/fslinterpolator_src.cpp

+/*