diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..b970088
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "delaunator-cpp"]
+	path = delaunator-cpp
+	url = https://github.com/abellgithub/delaunator-cpp.git
diff --git a/delaunator-cpp b/delaunator-cpp
new file mode 160000
index 0000000..6f28799
--- /dev/null
+++ b/delaunator-cpp
@@ -0,0 +1 @@
+Subproject commit 6f2879967bc96a9bcdbacf418e560e9f2e170ace
diff --git a/python/galario_defs.pxd b/python/galario_defs.pxd
index fd06835..c214325 100644
--- a/python/galario_defs.pxd
+++ b/python/galario_defs.pxd
@@ -23,8 +23,10 @@ cdef extern from "galario_py.h" namespace "galario":
     # Main user functions
     void _sample_profile(int nr, void* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis) except +
     void _sample_image(int nx, int ny, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis) except +
+    void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis) except +
     dreal _chi2_profile(int nr, void* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* vis_obs_w) except +
     dreal _chi2_image(int nx, int ny, void* image, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* vis_obs_w) except +
+    dreal _chi2_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* data, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights) except +
     void _sweep(int nr, void* intensity, dreal Rmin, dreal dR, int nxy, dreal dxy, dreal inc, void* image) except +
     void _uv_rotate(dreal PA, dreal dRA, dreal dDec, void* dRArot, void* dDecrot, int nd, void* u, void* v, void* urot, void* vrot) except +
 
diff --git a/python/libcommon.pyx b/python/libcommon.pyx
index 5420888..d9a0f71 100644
--- a/python/libcommon.pyx
+++ b/python/libcommon.pyx
@@ -22,9 +22,13 @@ from cpython cimport PyObject, Py_INCREF
 
 # Numpy must be initialized. When using numpy from C or Cython you must
 # _always_ do that, or you will have segfaults
+import time
 import numpy as np
 np.import_array()
 
+from scipy.interpolate import LinearNDInterpolator
+from scipy.spatial import Voronoi, ConvexHull
+
 include "galario_config.pxi"
 
 cimport galario_defs as cpp
@@ -33,7 +37,7 @@ __all__ = ['arcsec', 'deg', 'cgs_to_Jy', 'pc', 'au',
            '_init', '_cleanup', 'set_v_origin',
            'ngpus', 'use_gpu', 'threads',
            'check_obs', 'check_image_size', 'get_image_size',
-           'sampleImage', 'sampleProfile', 'chi2Image', 'chi2Profile',
+           'sampleImage', 'sampleUnstructuredImage', 'sampleProfile', 'chi2Image', 'chi2Profile',
            'get_coords_meshgrid',
            'sweep', 'uv_rotate', 'interpolate', 'apply_phase_vis', 'reduce_chi2',
            '_fft2d', '_fftshift', '_fftshift_axis0']
@@ -443,6 +447,93 @@ def sampleImage(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v,
     return vis
 
 
+def sampleUnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, 
+                int nxy, dreal dxy, dreal[::1] u, dreal[::1] v,
+                dRA=0., dDec=0., PA=0., check=False, origin='upper'):
+    """
+    Compute the synthetic visibilities of a model image at the specified (u, v) locations.
+
+    The 2D surface brightness in `image` is Fourier transformed and sampled in the
+    (u, v) locations given in the `u` and `v` arrays.
+
+    Typical call signature::
+
+        vis = sampleImage(image, dxy, u, v, dRA=0, dDec=0, PA=0, check=False, origin='upper')
+
+    Parameters
+    ----------
+    x : 1D array_like, float
+        List of x coordinates at which intensities are known.
+        **units**: rad
+    y : 1D array_like, float
+        List of y coordinates at which intensities are known.
+        **units**: rad
+    image : 1D array_like, float
+        Array containing the surface brightness of the model.
+        Assume the x-axis (R.A.) increases from right (West) to left (East)
+        and the y-axis (Dec.) increases from bottom (South) to top (North).
+        `nxy` must be even.
+        **units**: Jy/st
+    nxy : int
+        Number of pixels to use for the interpolated gridded image.
+    dxy : float
+        Size of the image cell in the interpolated image, assumed equal in both x and y direction.
+        **units**: rad
+    u : array_like, float
+        u coordinate of the visibility points where the FT has to be sampled.
+        **units**: wavelength
+    v : array_like, float
+        v coordinate of the visibility points where the FT has to be sampled.
+        The length of v must be equal to the length of u.
+        **units**: wavelength
+    dRA : float, optional
+        R.A. offset w.r.t. the phase center by which the image is translated.
+        If dRA > 0 translate the image towards the left (East). Default is 0.
+        **units**: rad
+    dDec : float, optional
+        Dec. offset w.r.t. the phase center by which the image is translated.
+        If dDec > 0 translate the image towards the top (North). Default is 0.
+        **units**: rad
+    PA : float, optional
+        Position Angle, defined East of North. Default is 0.
+        **units**: rad
+    check : bool, optional
+        If True, check whether `image` and `dxy` satisfy Nyquist criterion for
+        computing the synthetic visibilities in the (u, v) locations provided.
+        Additionally check that the (u, v) points fall in the image to avoid
+        segmentation violations. Default is False since the check might take
+        time. For executions where speed is important, set to False.
+    origin : ['upper' | 'lower'], optional
+        Set the [0,0] pixel index of the matrix in the upper left or lower left corner of the axes.
+        It follows the same convention as in matplotlib `matshow` and `imshow` commands.
+        Declination axis and the matrix y axis are parallel for `origin='lower'`, anti-parallel for `origin='upper'`.
+        The central pixel corresponding to the (RA, Dec) = (0, 0) is always [Nxy/2, Nxy/2].
+        For more details see the Technical Requirements page in the online docs.
+
+    Returns
+    -------
+    vis : array_like, complex
+        Synthetic visibilities sampled in the (u, v) locations given in `u` and `v`.
+        **units**: Jy
+
+    """
+
+    # Now pick back up with what is typically done for regular grids.
+    duv = 1 / (dxy*nxy)
+
+    if check:
+        check_image_size(u, v, nxy, dxy, duv)
+
+    vis = np.zeros(len(u), dtype=complex_dtype)
+    v_origin = set_v_origin(origin)
+    cpp._sample_unstructured_image(<void*>&x[0], <void*>&y[0], nxy, nxy, dxy, len(x), <void*>&image[0], v_origin, dRA, dDec, duv, PA, len(u), <void*>&u[0], <void*>&v[0], <void*>np.PyArray_DATA(vis))
+
+    return vis
+
+
+
+
+
 def sampleProfile(dreal[::1] intensity, Rmin, dR, nxy, dxy, dreal[::1] u, dreal[::1] v,
                   dRA=0., dDec=0., PA=0., inc=0., check=False):
     """
@@ -626,6 +717,110 @@ def chi2Image(dreal[:,::1] image, dxy, dreal[::1] u, dreal[::1] v,
     return cpp._chi2_image(image.shape[0], image.shape[1], <void*>&image[0,0], v_origin, dRA, dDec, duv, PA, len(u), <void*> &u[0],  <void*> &v[0],  <void*>&vis_obs_re[0], <void*>&vis_obs_im[0], <void*>&vis_obs_w[0])
 
 
+def chi2UnstructuredImage(dreal[::1] x, dreal[::1] y, dreal[::1] image, 
+              int nxy, dreal dxy, dreal[::1] u, dreal[::1] v,
+              dreal[::1] vis_obs_re, dreal[::1] vis_obs_im, dreal[::1] vis_obs_w,
+              dRA=0., dDec=0., PA=0., check=False, origin='upper'):
+    """
+    Compute the chi square of a model unstructured image given the observed visibilities.
+
+    The chi square is computed from the observed and synthetic visibilities as:
+
+    .. math::
+
+        \chi^2 = \sum_{j=1}^N w_j * [(Re V_{obs\ j}-Re V_{mod\ j})^2 + (Im V_{obs\ j}-Im V_{mod\ j})^2]
+
+    where :math:`V_{mod}` are the synthetic visibilities, which are computed internally
+    as in :func:`.sampleUnstructuredImage`.
+
+    Typical call signature::
+
+        chi2 = chi2UnstructuredImage(x, y, image, nxy, dxy, u, v, vis_obs_re, vis_obs_im, vis_obs_w,
+                         dRA=0, dDec=0, PA=0, check=False, origin='upper')
+
+    Parameters
+    ----------
+    x : 1D array_like, float
+        List of x coordinates at which intensities are known.
+        **units**: rad
+    y : 1D array_like, float
+        List of y coordinates at which intensities are known.
+        **units**: rad
+    image : 1D array_like, float
+        Array containing the surface brightness of the model.
+        Assume the x-axis (R.A.) increases from right (West) to left (East)
+        and the y-axis (Dec.) increases from bottom (South) to top (North).
+        `nxy` must be even.
+        **units**: Jy/st
+    nxy : int
+        Number of pixels to use for the interpolated gridded image.
+    dxy : float
+        Size of the image cell in the interpolated image, assumed equal in both x and y direction.
+        **units**: rad
+    u : array_like, float
+        u coordinate of the visibility points where the FT has to be sampled.
+        **units**: wavelength
+    v : array_like, float
+        v coordinate of the visibility points where the FT has to be sampled.
+        The length of `v` must be equal to the length of `u`.
+        **units**: wavelength
+    vis_obs_re : array_like, float
+        Real part of the observed visibilities.
+        **units**: Jy
+    vis_obs_im: array_like, float
+        Imaginary part of the observed visibilities.
+        The length of `vis_obs_im` must be equal to the length of `vis_obs_re`.
+        **units**: Jy
+    vis_obs_w: array_like, float
+        Weight associated to the observed visibilities.
+        The length of `vis_obs_w` must be equal to the length of `vis_obs_re`.
+        **units**:
+    dRA : float, optional
+        R.A. offset w.r.t. the phase center by which the image is translated.
+        If dRA > 0 translate the image towards the left (East). Default is 0.
+        **units**: rad
+    dDec : float, optional
+        Dec. offset w.r.t. the phase center by which the image is translated.
+        If dDec > 0 translate the image towards the top (North). Default is 0.
+        **units**: rad
+    PA : float, optional
+        Position Angle, defined East of North. Default is 0.
+        **units**: rad
+    check : bool, optional
+        If True, check whether `image` and `dxy` satisfy Nyquist criterion for
+        computing the synthetic visibilities in the (u, v) locations provided.
+        Additionally check that the (u, v) points fall in the image to avoid
+        segmentation violations. Default is False since the check might take
+        time. For executions where speed is important, set to False.
+    origin : ['upper' | 'lower'], optional
+        Set the [0,0] pixel index of the matrix in the upper left or lower left corner of the axes.
+        It follows the same convention as in matplotlib `matshow` and `imshow` commands.
+        Declination axis and the matrix y axis are parallel for `origin='lower'`, anti-parallel for `origin='upper'`.
+        The central pixel corresponding to the (RA, Dec) = (0, 0) is always [Nxy/2, Nxy/2].
+        For more details see the Technical Requirements page in the online docs.
+
+    Returns
+    -------
+    chi2: float
+        The chi square, not normalized.
+
+    See also
+    --------
+    :func:`.sampleImage`
+
+    """
+    check_obs(vis_obs_re, vis_obs_im, vis_obs_w, u=u, v=v)
+
+    duv = 1 / (dxy*nxy)
+
+    if check:
+        check_image_size(u, v, nxy, dxy, duv)
+
+    v_origin = set_v_origin(origin)
+
+    return cpp._chi2_unstructured_image(<void*>&x[0], <void*>&y[0], nxy, nxy, dxy, len(x), <void*>&image[0], v_origin, dRA, dDec, duv, PA, len(u), <void*>&u[0], <void*>&v[0], <void*>&vis_obs_re[0], <void*>&vis_obs_im[0], <void*>&vis_obs_w[0])
+
+
 def chi2Profile(dreal[::1] intensity, Rmin, dR, nxy, dxy, dreal[::1] u, dreal[::1] v,
                 dreal[::1] vis_obs_re, dreal[::1] vis_obs_im, dreal[::1] vis_obs_w,
                 dRA=0., dDec=0., PA=0., inc=0., check=False):
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 5c16505..953a52b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -71,7 +71,7 @@ OPTION(GALARIO_TIMING "Output timing of selected functions. For testing only. De
 
 foreach(t IN ITEMS galario_single galario)
   target_link_libraries(${t} ${FFTW3_LIBRARIES})
-  target_include_directories(${t} PUBLIC ${FFTW3_INCLUDE_DIRS} ${FFTW3_INCLUDE_DIR_PARALLEL})
+  target_include_directories(${t} PUBLIC ${FFTW3_INCLUDE_DIRS} ${FFTW3_INCLUDE_DIR_PARALLEL} ../delaunator-cpp/include)
   if(GALARIO_TIMING)
     target_compile_definitions(${t} PRIVATE GALARIO_TIMING=1)
   endif()
diff --git a/src/galario.cpp b/src/galario.cpp
index 3d31e88..e35cc62 100644
--- a/src/galario.cpp
+++ b/src/galario.cpp
@@ -19,6 +19,9 @@
 
 #include "galario.h"
 #include "galario_py.h"
+#include <cassert>
+#include <delaunator-header-only.hpp>
+#include <unordered_map>
 
 // full function makes code hard to read
 #define tpb galario::threads()
@@ -1351,6 +1354,320 @@ void sample_h(int nx, int ny, dcomplex* data, const dreal v_origin, dreal dRA, d
 
 #endif
 
+/**
+ * Find the index of the triangle that a point is in using brute force.
+ */
+int find_triangle_bruteforce_h(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy) {
+
+    bool found_triangle = false;
+    int which_triangle = -1;
+
+    // Loop through all the triangles to brute force-find which one a point is in.
+    for (int k = 0; k < d->triangles.size(); k+=3) {
+        int ia = d->triangles[k];
+        double ax = x[ia];
+        double ay = y[ia];
+
+        int ib = d->triangles[k+1];
+        double bx = x[ib];
+        double by = y[ib];
+
+        int ic = d->triangles[k+2];
+        double cx = x[ic];
+        double cy = y[ic];
+
+        double vbx = bx - ax;
+        double vby = by - ay;
+        double vcx = cx - ax;
+        double vcy = cy - ay;
+
+        double det_vv2 = gx*vcy - gy*vcx;
+        double det_v0v2 = ax*vcy - ay*vcx;
+        double det_v1v2 = vbx*vcy - vby*vcx;
+        double det_vv1 = gx*vby - gy*vbx;
+        double det_v0v1 = ax*vby - ay*vbx;
+
+        double a = (det_vv2 - det_v0v2) / det_v1v2;
+        double b = -(det_vv1 - det_v0v1) / det_v1v2;
+
+        // We've found the right triangle, now interpolate.
+        if ((a > 0) & (b > 0) & (a + b < 1)) {
+            which_triangle = k;
+            found_triangle = true;
+            break;
+        }
+    }
+
+    return which_triangle;
+}
+
+/**
+ * Find which triangle a point is in using a directed walk.
+ */
+int find_triangle_directedwalk_h(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double *time) {
+    int which_triangle = -2;
+    int count = 0;
+    dreal eps = 1.0e-3;
+    bool found_triangle = false;
+    while (count < d->triangles.size() / (3*4)) {
+        int ia = d->triangles[start];
+        double ax = x[ia];
+        double ay = y[ia];
+        int ib = d->triangles[start+1];
+        double bx = x[ib];
+        double by = y[ib];
+        int ic = d->triangles[start+2];
+        double cx = x[ic];
+        double cy = y[ic];
+
+        double wa = ((by - cy)*(gx - cx) + (cx - bx)*(gy - cy)) / 
+            ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy));
+        double wb = ((cy - ay)*(gx - cx) + (ax - cx)*(gy - cy)) /
+            ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy));
+        double wc = 1 - wa - wb;
+
+        if (wa < -eps) {
+            start = d->halfedges[start+1];
+        } else if (wb < -eps) {
+            start = d->halfedges[start+2];
+        } else if (wc < -eps) {
+            start = d->halfedges[start+0];
+        } else {
+            which_triangle = start;
+            found_triangle = true;
+        }
+
+        if (start >= 0) {
+            start = start - start%3;
+            *last_good = start;
+        }
+        else
+            which_triangle = start;
+
+        if ((found_triangle) or (which_triangle == -1))
+            break;
+
+        count++;
+    }
+
+    return which_triangle;
+}
+
+/**
+ * First try to find the triangle index using a directed walk, and if that fails switch to brute force.
+ */
+int find_triangle_h(delaunator::Delaunator *d, const dreal *x, const dreal *y, dreal gx, dreal gy, int start, int* last_good, double* time) {
+    int which_triangle = find_triangle_directedwalk_h(d, x, y, gx, gy, start, last_good, time);
+    if (which_triangle == -2) {
+        printf("Switching to brute force \n");
+        which_triangle = find_triangle_bruteforce_h(d, x, y, gx, gy);
+    }
+
+    return which_triangle;
+}
+
+
+/**
+ * Run the Delauney triangulation.
+ */
+delaunator::Delaunator triangulate_h(int ni, const dreal* x, const dreal* y, dreal v_origin) {
+    // Set up the Delauney triangulation.
+
+    std::vector<double> coords;
+
+    dreal xmin = std::numeric_limits<dreal>::max(); dreal xmax = -std::numeric_limits<dreal>::max();
+    dreal ymin = std::numeric_limits<dreal>::max(); dreal ymax = -std::numeric_limits<dreal>::max();
+    for (int i=0; i < ni; i++) {
+        coords.push_back(x[i]);
+        coords.push_back(y[i]);
+
+        if (x[i] > xmax) xmax = x[i];
+        if (x[i] < xmin) xmin = x[i];
+        if (y[i] > ymax) ymax = y[i];
+        if (y[i] < ymin) ymin = y[i];
+    }
+
+    delaunator::Delaunator d(coords);
+
+    return d;
+}
+
+/**
+ * For each triangle, calculate the centroid and which grid cell it falls in.
+ */
+void bin_triangles_h(int nx, int ny, dreal dxy, const dreal *x, const dreal *y, const dreal *realdata, delaunator::Delaunator &d, std::unordered_map<int,dreal> &binned_image, 
+        std::unordered_map<int,dreal> &binned_weights, std::unordered_map<int,int> &npoints, dreal v_origin) {
+    auto tx = static_cast<dreal*>(malloc(sizeof(dreal)*d.triangles.size()/3));
+    auto ty = static_cast<dreal*>(malloc(sizeof(dreal)*d.triangles.size()/3));
+    auto tf = static_cast<dreal*>(malloc(sizeof(dreal)*d.triangles.size()/3));
+    auto ta = static_cast<dreal*>(malloc(sizeof(dreal)*d.triangles.size()/3));
+
+    auto itx = static_cast<int*>(malloc(sizeof(int)*d.triangles.size()/3));
+    auto ity = static_cast<int*>(malloc(sizeof(int)*d.triangles.size()/3));
+
+    dreal gx_max = 0.5*nx*dxy;
+    dreal gy_max = 0.5*ny*dxy*v_origin;
+
+    #pragma omp parallel for
+    for (int i = 0; i < d.triangles.size()/3; i++) {
+        tx[i] = (x[d.triangles[3*i]] + x[d.triangles[3*i+1]] + x[d.triangles[3*i+2]]) / 3.;
+        ty[i] = (y[d.triangles[3*i]] + y[d.triangles[3*i+1]] + y[d.triangles[3*i+2]]) / 3.;
+        tf[i] = (realdata[d.triangles[3*i]] + realdata[d.triangles[3*i+1]] + realdata[d.triangles[3*i+2]]) / 3.;
+        ta[i] = std::fabs((y[d.triangles[3*i+1]] - y[d.triangles[3*i]]) * (x[d.triangles[3*i+2]] - x[d.triangles[3*i+1]]) - 
+                (x[d.triangles[3*i+1]] - x[d.triangles[3*i]]) * (y[d.triangles[3*i+2]] - y[d.triangles[3*i+1]]));
+
+        itx[i] = trunc((tx[i] - gx_max) / (-dxy) + 0.5);
+        ity[i] = trunc((ty[i] - gy_max) / (-dxy*v_origin) + 0.5);
+    }
+
+    for (int i = 0; i < d.triangles.size()/3; i++) {
+        // Note: cant do this in parallel because two threads could access same
+        // grid cell at the same time. Locking made this very slow.
+        if ((itx[i] >= 0) and (itx[i] < nx) and (ity[i] >= 0) and (ity[i] < ny)) {
+            if (npoints.find(ity[i] * nx + itx[i]) == npoints.end()) {
+                npoints[ity[i] * nx + itx[i]] = 1;
+                binned_image[ity[i] * nx + itx[i]] = tf[i] * ta[i];
+                binned_weights[ity[i] * nx + itx[i]] = ta[i];
+            } else {
+                npoints[ity[i] * nx + itx[i]] += 1;
+                binned_image[ity[i] * nx + itx[i]] += tf[i] * ta[i];
+                binned_weights[ity[i] * nx + itx[i]] += ta[i];
+            }
+        }
+    }
+
+    std::unordered_map<int, int>::iterator it = npoints.begin();
+    while (it != npoints.end()) {
+        // Erase any places where npoints = 1
+        if (it->first <= 1)
+            it = npoints.erase(it);
+        else
+            it++;
+    }
+
+    free(tx); free(ty); free(tf); free(ta); free(itx); free(ity);
+}
+
+/**
+ * Do the interpolation onto a single point in a single triangle.
+ */
+double interpolate_on_triangle_h(delaunator::Delaunator &d, int which_triangle, const dreal *x, const dreal *y, const dreal *realdata, dreal gx, dreal gy) {
+    int ia = d.triangles[which_triangle];
+    double ax = x[ia];
+    double ay = y[ia];
+    int ib = d.triangles[which_triangle+1];
+    double bx = x[ib];
+    double by = y[ib];
+    int ic = d.triangles[which_triangle+2];
+    double cx = x[ic];
+    double cy = y[ic];
+
+    double wa = ((by - cy)*(gx - cx) + (cx - bx)*(gy - cy)) / 
+        ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy));
+    double wb = ((cy - ay)*(gx - cx) + (ax - cx)*(gy - cy)) /
+        ((by - cy)*(ax - cx) + (cx - bx)*(ay - cy));
+    double wc = 1 - wa - wb;
+
+    return wa*realdata[ia] + wb*realdata[ib] + wc*realdata[ic];
+}
+
+/**
+ * Interpolate when the triangles are bigger than the grid cells, and use the binned image when triangles are smaller.
+ */
+dreal* interpolate_or_bin_to_image_h(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* realdata, dreal v_origin, 
+        delaunator::Delaunator &d, std::unordered_map<int,dreal> &binned_image, std::unordered_map<int,dreal> &binned_weights, 
+        std::unordered_map<int,int> &npoints) {
+
+    // Get the max and min x and y values from the triangulation.
+    dreal xmin = std::numeric_limits<dreal>::max(); dreal xmax = -std::numeric_limits<dreal>::max();
+    dreal ymin = std::numeric_limits<dreal>::max(); dreal ymax = -std::numeric_limits<dreal>::max();
+    for (int i=0; i < ni; i++) {
+        if (x[i] > xmax) xmax = x[i];
+        if (x[i] < xmin) xmin = x[i];
+        if (y[i] > ymax) ymax = y[i];
+        if (y[i] < ymin) ymin = y[i];
+    }
+
+    // Create an image including the appropriate coordinates.
+    auto gx = static_cast<dreal*>(malloc(sizeof(dreal)*nx));
+    auto gy = static_cast<dreal*>(malloc(sizeof(dreal)*ny));
+    auto image = static_cast<dreal*>(malloc(sizeof(dreal)*nx*ny));
+
+    #pragma omp parallel for
+    for (int i = 0; i < nx; i++)
+        gx[i] = (0.5 - i * 1./nx) * nx * dxy;
+    #pragma omp parallel for
+    for (int i = 0; i < ny; i++)
+        gy[i] = (0.5 - i * 1./ny) * ny * dxy * v_origin;
+
+    #pragma omp parallel
+    {
+    int which_triangle = 0;
+    int last_triangle = 0;
+    int col_start_triangle = -1;
+    double time = 0.;
+
+    // Now loop through the pixels in the image pixels, find the triangle each point is in, and interpolate.
+    #pragma omp for schedule(static)
+    for (int i = 0; i < ny; i++) {
+        if ((i > 0) and (col_start_triangle > -1)) {
+            which_triangle = col_start_triangle;
+            last_triangle = col_start_triangle;
+            col_start_triangle = -1;
+        }
+        for (int j = 0; j < nx; j++) {
+            // Check whether the triangle is out of the triangulation.
+            if ((gx[j] > xmin) and (gx[j] < xmax) and (gy[i] > ymin) and (gy[i] < ymax))
+                // Find which triangle this grid point is in.
+                which_triangle = find_triangle_h(&d, x, y, gx[j], gy[i], which_triangle, &last_triangle, &time);
+            else
+                which_triangle = -1;
+
+            // We've found the right triangle, now interpolate.
+            if (which_triangle > -1) {
+                if (npoints.find(i * nx + j) != npoints.end())
+                    image[i * nx + j] = binned_image[i * nx + j] / binned_weights[i * nx + j] * dxy * dxy;
+                else
+                    image[i * nx + j] = interpolate_on_triangle_h(d, which_triangle, x, y, realdata, gx[j], gy[i])*dxy*dxy;
+
+                if (col_start_triangle == -1)
+                    col_start_triangle = last_triangle;
+            }
+            // If no triangle was found, the point is outside the area with data so set to 0.
+            else {
+                image[i * nx + j] = 0.;
+                which_triangle = last_triangle;
+            }
+        }
+    }
+    }
+
+    // Clean up
+    free(gx); free(gy);
+
+    return image;
+}
+
+/**
+ * Interpolate from an unstructured image onto a regular grid.
+ */
+dreal* unstructured_to_grid_h(int nx, int ny, int ni, dreal dxy, const dreal* x, const dreal* y, const dreal* realdata, dreal v_origin) {
+    // Set up the Delauney triangulation.
+    OPENMPTIME(delaunator::Delaunator d = triangulate_h(ni, x, y, v_origin), "unstructured_to_grid::triangulation");
+
+    // For each triangle, calculate the centroid and which grid cell it falls in.
+    std::unordered_map<int,dreal> binned_image;
+    std::unordered_map<int,dreal> binned_weights;
+    std::unordered_map<int,int> npoints;
+
+    OPENMPTIME(bin_triangles_h(nx, ny, dxy, x, y, realdata, d, binned_image, binned_weights, npoints, v_origin), "unstructured_to_grid::bin_trixels");
+
+    // Interpolate or bin, as appropriate to get to an image.
+    OPENMPTIME(auto image = interpolate_or_bin_to_image_h(nx, ny, ni, dxy, x, y, realdata, v_origin, d, binned_image, binned_weights, npoints), "unstructured_to_grid::generate_gridded_image");
+
+    return image;
+}
+
 
 namespace galario {
 /**
@@ -1396,6 +1713,52 @@ void _sample_image(int nx, int ny, void* data, dreal v_origin, dreal dRA, dreal
     sample_image(nx, ny, static_cast<dreal*>(data), v_origin, dRA, dDec, duv, PA, nd, static_cast<dreal*>(u), static_cast<dreal*>(v), static_cast<dcomplex*>(vis_int));
 }
 
+/**
+ * return result in `vis_int`
+ */
+void sample_unstructured_image(const dreal* realx, const dreal* realy, int nx, int ny, dreal dxy, int ni, const dreal* realdata, dreal v_origin, dreal dRA, dreal dDec, dreal duv,
+                          const dreal PA, int nd, const dreal* u, const dreal* v, dcomplex* vis_int) {
+    CPUTimer t_start;
+
+    // Initialization for uv_idx and interpolate
+    CHECK_INPUT(nx);
+
+/*#ifdef __CUDACC__
+    GPUTimer t_total;
+    CudaMemory<dcomplex> vis_int_d(nd);
+
+    auto data_d = copy_input_d(nx, ny, realdata);
+
+    // do the actual computation
+    sample_d(nx, ny, data_d.ptr, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int_d.ptr);
+
+    // retrieve interpolated values
+    CCheck(cudaDeviceSynchronize());
+
+    GPUTimer t;
+    vis_int_d.Retrieve(vis_int);
+    t.Elapsed("sample_image::vis_int_ D->H");
+
+    t_total.Elapsed("sample_image_tot");
+#else*/
+
+    auto data = unstructured_to_grid_h(nx, ny, ni, dxy, realx, realy, realdata, v_origin);
+
+    CPUTimer t; 
+    auto image = copy_input(nx, ny, data); t.Elapsed("sample_image::copy_input");
+
+    sample_h(nx, ny, image, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int);
+
+    t = CPUTimer(); galario_free(data); galario_free(image); t.Elapsed("sample_image::free_data");
+//#endif
+    t_start.Elapsed("sample_image_tot");
+}
+
+void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_int) {
+    sample_unstructured_image(static_cast<dreal*>(x), static_cast<dreal*>(y), nx, ny, dxy, ni, static_cast<dreal*>(data), v_origin, dRA, dDec, duv, PA, nd, static_cast<dreal*>(u), static_cast<dreal*>(v), static_cast<dcomplex*>(vis_int));
+}
+
+
 
 /**
  * return result in `vis_int`
@@ -1612,6 +1975,65 @@ dreal _chi2_image(int nx, int ny, void* realdata, const dreal v_origin, dreal dR
                  static_cast<dreal*>(weights));
 }
 
+dreal chi2_unstructured_image(const dreal* realx, const dreal* realy, int nx, int ny, dreal dxy, int ni, const dreal* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, const dreal* vis_obs_re, const dreal* vis_obs_im, const dreal* weights) {
+    CPUTimer t_start;
+
+    CHECK_INPUTXY(nx, ny);
+    dreal chi2 = 0;
+#ifdef __CUDACC__
+    GPUTimer t;
+     // ################################
+     // ### ALLOCATION, INITIALIZATION ###
+     // ################################
+
+     /* async memory copy:
+      TODO copy memory asynchronously or create streams to define dependencies
+      use nonzero cudaStream_t
+      kernel<<< blocks, threads, bytes=0, stream =! 0>>>();
+
+      all cufft calls are asynchronous, can specify the stream explicitly (cf. doc)
+      same for cublas
+      draw dependcies on paper: first thing is to do fft while other data is transferred
+
+      While the FFT etc. are calculated, we can copy over the weights and observed values.
+     */
+    // reserve memory for the interpolated values
+    //CudaMemory<dcomplex> vis_int_d(nd);
+    //t.Elapsed("chi2_image::malloc_vis_int");
+
+    // Initialization for comparison and chi square computation
+    /* allocate and copy observational data */
+    /*CudaMemory<dreal> vis_obs_re_d(nd, vis_obs_re);
+    CudaMemory<dreal> vis_obs_im_d(nd, vis_obs_im);
+    CudaMemory<dreal> weights_d(nd, weights);
+    t.Elapsed("chi2_image::copy_observations");
+
+    auto data_d = copy_input_d(nx, ny, realdata);
+
+    sample_d(nx, ny, data_d.ptr, v_origin, dRA, dDec, nd, duv, PA, u, v, vis_int_d.ptr);
+    chi2 = reduce_chi2_d(nd, vis_obs_re_d.ptr, vis_obs_im_d.ptr, vis_int_d.ptr, weights_d.ptr);*/
+#else
+    CPUTimer t;
+
+    auto vis_int = reinterpret_cast<dcomplex*>(FFTW(alloc_complex)(nd)); t.Elapsed("chi2_imag::fftw_alloc");
+    sample_unstructured_image(realx, realy, nx, ny, dxy, ni, realdata, v_origin, dRA, dDec, duv, PA, nd, u, v, vis_int);
+
+    chi2 = reduce_chi2(nd, vis_obs_re, vis_obs_im, vis_int, weights);
+
+    t = CPUTimer(); galario_free(vis_int); t.Elapsed("chi2_imag::free_vis_int");
+#endif
+    t_start.Elapsed("chi2_image_tot");
+    flush_timing();
+
+    return chi2;
+}
+
+dreal _chi2_unstructured_image(void* realx, void* realy, int nx, int ny, dreal dxy, int ni, void* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights) {
+    return chi2_unstructured_image(static_cast<dreal*>(realx), static_cast<dreal*>(realy), nx, ny, dxy, ni, static_cast<dreal*>(realdata), v_origin, dRA, dDec, duv, PA, nd, 
+                 static_cast<dreal*>(u), static_cast<dreal*>(v), static_cast<dreal*>(vis_obs_re), static_cast<dreal*>(vis_obs_im),
+                 static_cast<dreal*>(weights));
+}
+
 dreal chi2_profile(int nr, dreal *const intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA,
                           dreal dDec, dreal duv, dreal PA, int nd, const dreal *u, const dreal *v, const dreal *vis_obs_re,
                           const dreal *vis_obs_im, const dreal *weights) {
diff --git a/src/galario.h b/src/galario.h
index 9a3c4ab..1787edc 100644
--- a/src/galario.h
+++ b/src/galario.h
@@ -27,10 +27,12 @@ namespace galario {
 void sample_profile(int nr, const dreal* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA,
                     dreal dDec, dreal duv, dreal PA, int nd, const dreal *u, const dreal *v, dcomplex *vis_int);
 void sample_image(int nx, int ny, const dreal* image, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, dcomplex* vis_int);
+void sample_unstructured_image(const dreal* x, const dreal *y, int nx, int ny, dreal dxy, int ni, const dreal* image, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, dcomplex* vis_int);
 dreal chi2_profile(int nr, const dreal* intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA,
                    dreal dDec, dreal duv, dreal PA, int nd, const dreal *u, const dreal *v, const dreal *vis_obs_re,
                    const dreal *vis_obs_im, const dreal *weights);
 dreal chi2_image(int nx, int ny, const dreal* image, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, const dreal* vis_obs_re, const dreal* vis_obs_im, const dreal* weights);
+dreal chi2_unstructured_image(const dreal* realx, const dreal* realy, int nx, int ny, dreal dxy, int ni, const dreal* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, const dreal* u, const dreal* v, const dreal* vis_obs_re, const dreal* vis_obs_im, const dreal* weights);
 void sweep(int nr, const dreal* intensity, dreal Rmin, dreal dR, int nxy, dreal dxy, dreal inc, dcomplex *image);
 void uv_rotate(dreal PA, dreal dRA, dreal dDec, dreal* dRArot, dreal* dDecrot, int nd, const dreal* u, const dreal* v, dreal* urot, dreal* vrot);
 
diff --git a/src/galario_py.h b/src/galario_py.h
index 1321fa8..fca8cca 100644
--- a/src/galario_py.h
+++ b/src/galario_py.h
@@ -30,9 +30,11 @@ namespace galario {
 void _sample_profile(int nr, void *intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec,
                      dreal duv, dreal PA, int nd, void *u, void *v, void *vis_int);
 void _sample_image(int nx, int ny, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_int);
+void _sample_unstructured_image(void* x, void* y, int nx, int ny, dreal dxy, int ni, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_int);
 dreal _chi2_profile(int nr, void *intensity, dreal Rmin, dreal dR, dreal dxy, int nxy, dreal inc, dreal dRA, dreal dDec,
                     dreal duv, dreal PA, int nd, void *u, void *v, void *vis_obs_re, void *vis_obs_im, void *weights);
 dreal _chi2_image(int nx, int ny, void* data, dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights);
+dreal _chi2_unstructured_image(void* realx, void* realy, int nx, int ny, dreal dxy, int ni, void* realdata, const dreal v_origin, dreal dRA, dreal dDec, dreal duv, dreal PA, int nd, void* u, void* v, void* vis_obs_re, void* vis_obs_im, void* weights);
 void _sweep(int nr, void *intensity, dreal Rmin, dreal dR, int nxy, dreal dxy, dreal inc, void *image);
 void _uv_rotate(dreal PA, dreal dRA, dreal dDec, void* dRArot, void* dDecrot, int nd, void* u, void* v, void* urot, void* vrot);
 
diff --git a/test/test_centering.py b/test/test_centering.py
new file mode 100755
index 0000000..bf751e4
--- /dev/null
+++ b/test/test_centering.py
@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+
+from galario import double
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import pyDOE
+import numpy
+
+# Make an image.
+
+grid = pyDOE.lhs(2, samples=400)
+
+r = grid[:,0] * 2
+phi = grid[:,1] * 2*numpy.pi
+
+grid = pyDOE.lhs(2, samples=2000)
+
+r = numpy.hstack((r, grid[:,0]*0.04 + 0.98))
+phi = numpy.hstack((phi, grid[:,1]*2*numpy.pi))
+
+x = r * numpy.cos(phi)
+y = r * numpy.sin(phi)
+
+flux = numpy.where(r < 1., 1., 0.)
+
+# Do the Fourier transform with TrIFT
+
+u, v = numpy.meshgrid(numpy.linspace(-3.,3.,100),numpy.linspace(-3.,3.,100))
+
+u = u.reshape((u.size,))
+v = v.reshape((v.size,))
+
+vis = double.sampleUnstructuredImage(x, y, flux, 4096, 0.02, u, v, 0.5, 0.25)
+
+# Now shift the image manually.
+
+x += 0.5
+y += 0.25
+
+vvis = double.sampleUnstructuredImage(x, y, flux, 4096, 0.02, u, v, 0., 0.)
+
+# Plot the image to make sure we did the correct shifting.
+
+triang = tri.Triangulation(x, y)
+
+fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,4))
+
+ax.tripcolor(triang, flux, "ko-")
+ax.triplot(triang, "k.-", linewidth=0.1, markersize=0.1)
+
+ax.set_aspect("equal")
+
+ax.set_xlim(1.6,-1.6)
+ax.set_ylim(-1.6,1.6)
+
+ax.set_xlabel("x", fontsize=14)
+ax.set_ylabel("y", fontsize=14)
+
+ax.tick_params(labelsize=14)
+
+plt.show()
+
+# Finally, plot the visibilities.
+
+fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8,8))
+
+ax[0,0].scatter(u, v, c=vis.real, marker=".")
+ax[0,1].scatter(u, v, c=vis.imag, marker=".")
+
+ax[1,0].scatter(u, v, c=vvis.real, marker=".")
+ax[1,1].scatter(u, v, c=vvis.imag, marker=".")
+
+for a in ax.flatten():
+    a.set_xlabel("u", fontsize=14)
+    a.set_ylabel("v", fontsize=14)
+
+    a.tick_params(labelsize=14)
+
+plt.show()
diff --git a/test/test_circle.py b/test/test_circle.py
new file mode 100755
index 0000000..880c57c
--- /dev/null
+++ b/test/test_circle.py
@@ -0,0 +1,77 @@
+#!/usr/bin/env python3
+
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import scipy.special
+import pyDOE
+import numpy
+import time
+
+from galario import double
+
+# Make an image.
+
+grid = pyDOE.lhs(2, samples=100)
+
+r = grid[:,0] * 2
+phi = grid[:,1] * 2*numpy.pi
+
+grid = pyDOE.lhs(2, samples=2000)
+
+r = numpy.hstack((r, grid[:,0]*0.04 + 0.98))
+phi = numpy.hstack((phi, grid[:,1]*2*numpy.pi))
+
+x = r * numpy.cos(phi)
+y = r * numpy.sin(phi)
+
+flux = numpy.where(r < 1., 1., 0.)
+
+# Plot the image.
+
+triang = tri.Triangulation(x, y)
+
+plt.tripcolor(triang, flux, "ko-")
+plt.triplot(triang, "k.-", linewidth=0.1, markersize=0.1)
+
+plt.axes().set_aspect("equal")
+
+plt.xlim(-1.1,1.1)
+plt.ylim(-1.1,1.1)
+
+plt.xlabel("x", fontsize=14)
+plt.ylabel("y", fontsize=14)
+
+plt.axes().tick_params(labelsize=14)
+
+plt.show()
+
+# Do the Fourier transform with TrIFT
+
+u = numpy.linspace(0.001,10.,1000)
+v = numpy.repeat(0., 1000)
+
+t1 = time.time()
+vis = double.sampleUnstructuredImage(x, y, flux, 4096, 0.02, u, v, 0.25, 0.25)
+t2 = time.time()
+print(t2 - t1)
+
+# Calculate the analytic result.
+
+vis_analytic = scipy.special.jv(1, 2*numpy.pi*u) / u * numpy.exp(2*numpy.pi*\
+        1j*(0.25*u + 0.25*v))
+
+# Finally, plot the visibilities.
+
+plt.plot(u, vis.real, "k.-", label="Unstructured Fourier Transform")
+plt.plot(u, vis_analytic.real, "r-", label="Analytic Solution")
+
+plt.xlabel("u", fontsize=14)
+plt.ylabel("Real Component", fontsize=14)
+
+plt.legend(fontsize=14)
+
+plt.axes().tick_params(labelsize=14)
+
+plt.subplots_adjust(left=0.17, right=0.95, top=0.99)
+
+plt.show()
diff --git a/test/test_orientation.py b/test/test_orientation.py
new file mode 100755
index 0000000..731a028
--- /dev/null
+++ b/test/test_orientation.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+
+from galario import double
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+import pyDOE
+import numpy
+
+# Make an image.
+
+grid = pyDOE.lhs(2, samples=400)
+
+r = grid[:,0] * 2
+phi = grid[:,1] * 2*numpy.pi
+
+grid = pyDOE.lhs(2, samples=2000)
+
+r = numpy.hstack((r, grid[:,0]*0.04 + 0.98))
+phi = numpy.hstack((phi, grid[:,1]*2*numpy.pi))
+
+x = r * numpy.cos(phi) * numpy.cos(numpy.pi/3)
+y = r * numpy.sin(phi)
+
+flux = numpy.where(r < 1., y - y.min(), 0.)
+
+pa = numpy.pi/4
+
+xp = x * numpy.cos(-pa) - y * numpy.sin(-pa)
+yp = x * numpy.sin(-pa) + y * numpy.cos(-pa)
+
+x = xp
+y = yp
+
+# Also make a traditional image to compare with.
+
+xx, yy = numpy.meshgrid(numpy.linspace(15.,-15.,1024, endpoint=False), \
+        numpy.linspace(15.,-15.,1024, endpoint=False))
+
+xp = xx * numpy.cos(pa) - yy * numpy.sin(pa)
+yp = xx * numpy.sin(pa) + yy * numpy.cos(pa)
+
+rr = numpy.sqrt((xp/numpy.cos(numpy.pi/3))**2 + yp**2)
+
+fflux = numpy.where(rr < 1., yp - y.min(), 0.)
+
+# Plot the image.
+
+triang = tri.Triangulation(x, y)
+
+fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8,4))
+
+ax[0].tripcolor(triang, flux, "ko-")
+ax[0].triplot(triang, "k.-", linewidth=0.1, markersize=0.1)
+
+ax[1].imshow(fflux, interpolation="nearest")
+
+for i in range(1):
+    ax[i].set_aspect("equal")
+
+    ax[i].set_xlim(1.1,-1.1)
+    ax[i].set_ylim(-1.1,1.1)
+
+    ax[i].set_xlabel("x", fontsize=14)
+    ax[i].set_ylabel("y", fontsize=14)
+
+    ax[i].tick_params(labelsize=14)
+
+plt.show()
+
+# Do the Fourier transform with TrIFT
+
+u, v = numpy.meshgrid(numpy.linspace(-3.,3.,100),numpy.linspace(-3.,3.,100))
+
+u = u.reshape((u.size,))
+v = v.reshape((v.size,))
+
+vis = double.sampleUnstructuredImage(x, y, flux, 1024*4, 0.025/4, u, v, 0., 0.)
+
+# Do the Fourier transform with GALARIO.
+
+dxy = abs(xx[0,1] - xx[0,0])
+
+vvis = double.sampleImage(fflux, dxy, u, v)
+
+# Finally, plot the visibilities.
+
+fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8,8))
+
+ax[0,0].scatter(u, v, c=vis.real/vis.real.max(), marker=".")
+ax[0,1].scatter(u, v, c=vis.imag/vis.imag.max(), marker=".")
+
+ax[1,0].scatter(u, v, c=vvis.real/vvis.real.max(), marker=".")
+ax[1,1].scatter(u, v, c=vvis.imag/vvis.imag.max(), marker=".")
+
+for a in ax.flatten():
+    a.set_xlabel("u", fontsize=14)
+    a.set_ylabel("v", fontsize=14)
+
+    a.tick_params(labelsize=14)
+
+plt.show()