interp3d.pyx

import numpy as np
cimport numpy as np
from libc.stdlib cimport free, malloc

########### from utils.h

ctypedef struct stack:
    int top
    int slots
    void** arr

ctypedef struct arrayList:
    int   num_slots
    int   num_elements
    void** arr

cdef struct _listNode:
    void *data
    _listNode *next
    _listNode *prev

ctypedef _listNode listNode

ctypedef struct linkedList:
    listNode *head
    listNode *last
    int nelem
    stack   *deadNodes

############ from delaunay.h

ctypedef struct vertex:
    ### This is the location of this point.
    double v[3]

    ### This is the point index in the point list.
    ### We only store this so that it is convenient for using the
    ### tet-mesh function in Matlab.
    ### We can remove it for when the code is actually used.
    int    index

    ### These are the values for our vector field at this location.
    double data[3]

    ### We use this for caching the voronoi volume of this point.
    ### it will provide a good speed-up!
    double voronoiVolume

cdef struct _simplex:
    ### The verticies of this simplex.
    vertex  *p[4]
    ### The neighbouring simlpicies of this simplex.
    ### These are ordered in accordance with our 'get face' routine:
    ### so that the i'th face is shared with the i'th neighbour.
    _simplex *s[4]
    ### This is the node in our auxillary list structure that holds this simplex.
    ### It's far from an elegant solution: but it is fast and space-efficient.
    listNode *node

ctypedef _simplex simplex

ctypedef struct neighbourUpdate:
    stack  *ptrs
    stack  *old

ctypedef struct mesh:
    ### a linked list of all the simplicies.
    linkedList    *tets

    ### The simplex which contains all of the points.
    ### its verticies contain no data values.
    simplex *super
    vertex   superVerticies[4]

    ### Memory pool.
    stack   *deadSimplicies
    stack   *deadVoronoiCells

    ### We modify these when a point is inserted/removed.
    arrayList       *conflicts
    arrayList       *updates
    neighbourUpdate *neighbourUpdates

    ### Keep count of the number of degenerecies we find in the mesh,
    ### so that we can spot errors, and be aware of particularly degenerate data.
    int coplanar_degenerecies
    int cospherical_degenerecies



cdef extern from "3d_interp_v1/natural.h":
    void          interpolate3_3(double x, double y, double z, double *u, double *v, double *w, mesh *m)
    vertex       *loadPoints(char *filename, int *n)
    vertex       *initPoints(double *x, double *y, double *z, double *u, double *v, double *w, int n)
    void          writePointsToFile(vertex *ps, int n)
    # void        lastNaturalNeighbours(vertex *v, mesh *m, arrayList *neighbours, arrayList *neighbourSimplicies)

cdef extern from "3d_interp_v1/delaunay.h":
    mesh         *newMesh()
    void          freeMesh(mesh *m)
    void          buildMesh(vertex *ps, int n, mesh *m)

DT = np.float
ctypedef np.float_t DT_t

########### begin python-available function

cdef class interp3d:

    cdef int n
    cdef vertex* ps
    cdef mesh* workingMesh

    def __init__( self, np.ndarray[DT_t, ndim=1] x, np.ndarray[DT_t, ndim=1] y, np.ndarray[DT_t, ndim=1] z,    \
                         np.ndarray[DT_t, ndim=1] u, np.ndarray[DT_t, ndim=1] v, np.ndarray[DT_t, ndim=1] w):

        # {
        assert x.shape == y.shape == z.shape == u.shape == v.shape == w.shape
        self.n = x.size

        cdef double *xC = <double *> malloc(self.n * sizeof(double))
        cdef double *yC = <double *> malloc(self.n * sizeof(double))
        cdef double *zC = <double *> malloc(self.n * sizeof(double))
        cdef double *uC = <double *> malloc(self.n * sizeof(double))
        cdef double *vC = <double *> malloc(self.n * sizeof(double))
        cdef double *wC = <double *> malloc(self.n * sizeof(double))


        cdef unsigned int i
        for i in range(self.n):
            xC[i] = x[i]
        for i in range(self.n):
            yC[i] = y[i]
        for i in range(self.n):
            zC[i] = z[i]
        for i in range(self.n):
            uC[i] = u[i]
        for i in range(self.n):
            vC[i] = v[i]
        for i in range(self.n):
            wC[i] = u[i]


        self.ps = initPoints(xC, yC, zC, uC, vC, wC, self.n)
        self.workingMesh = newMesh()
        buildMesh(self.ps,self.n, self.workingMesh)
        # }

    cpdef np.ndarray[DT_t, ndim=1] interpolatePoint(self, double x, double y, double z):

        # {
        cdef double uOut
        cdef double vOut
        cdef double wOut

        interpolate3_3(x, y, z, &uOut, &vOut, &wOut, self.workingMesh)

        cdef np.ndarray[DT_t, ndim=1] returnArray
        returnArray = np.array([ uOut, vOut, wOut ])
        return returnArray
        # }

    def clearData(self):

        # {
        freeMesh(self.workingMesh)
        free(self.ps)
        free(&self.n)
        # }