natneighbors.cpp

#include "delaunay_utils.h"
#include "natneighbors.h"

#include <stack>
#include <list>
#include <vector>
#include <set>
#include <math.h>
#include <iostream>
#include <iterator>

using namespace std;

NaturalNeighbors::NaturalNeighbors(int npoints, int ntriangles, double *x, double *y,
    double *centers, int *nodes, int *neighbors)
{
    this->npoints = npoints;
    this->ntriangles = ntriangles;
    this->x = x;
    this->y = y;
    this->centers = centers;
    this->nodes = nodes;
    this->neighbors = neighbors;

    this->radii2 = new double[ntriangles];
    for (int i=0; i<ntriangles; i++) {
        double x2 = x[INDEX3(nodes,i,0)] - INDEX2(centers,i,0);
        x2 = x2*x2;
        double y2 = y[INDEX3(nodes,i,0)] - INDEX2(centers,i,1);
        y2 = y2*y2;
        this->radii2[i] = x2 + y2;
    }

}

NaturalNeighbors::~NaturalNeighbors()
{
    delete[] this->radii2;
}

int NaturalNeighbors::find_containing_triangle(double targetx, double targety, int start_triangle)
{
    int final_triangle;
    final_triangle = walking_triangles(start_triangle, targetx, targety, 
        x, y, nodes, neighbors);
    return final_triangle;
}

double NaturalNeighbors::interpolate_one(double *z, double targetx, double targety,
    double defvalue, int &start_triangle)
{
    int t = find_containing_triangle(targetx, targety, start_triangle);
    if (t == -1) return defvalue;

    start_triangle = t;
    vector<int> circumtri;

    circumtri.push_back(t);
    stack<int> stackA;
    stack<int> stackB;
    int tnew, i;

    for (i=0; i<3; i++) {
        tnew = INDEX3(this->neighbors, t, i);
        if (tnew != -1) {
            stackA.push(tnew);
            stackB.push(t);
        }
    }
    while (!stackA.empty()) {
        tnew = stackA.top();
        stackA.pop();
        t = stackB.top();
        stackB.pop();
        double d2 = (SQ(targetx - INDEX2(this->centers,tnew,0))
                   + SQ(targety - INDEX2(this->centers,tnew,1)));
        if ((this->radii2[tnew]-d2) > TOLERANCE_EPS) {
            // tnew is a circumtriangle of the target
            circumtri.push_back(tnew);
            for (i=0; i<3; i++) {
                int ti = INDEX3(this->neighbors, tnew, i);
                if ((ti != -1) && (ti != t)) {
                    stackA.push(ti);
                    stackB.push(tnew);
                }
            }
        }
    }

    vector<int>::iterator it;
    double f = 0.0;
    double A = 0.0;
    double tA=0.0, yA=0.0, cA=0.0; // Kahan summation temps for A
    double tf=0.0, yf=0.0, cf=0.0; // Kahan summation temps for f

    vector<int> edge;
    bool onedge = false;
    bool onhull = false;

    for (it = circumtri.begin(); it != circumtri.end(); it++) {
        int t = *it;
        double vx = INDEX2(this->centers, t, 0);
        double vy = INDEX2(this->centers, t, 1);
        vector<double> c(6);
        for (int i=0; i<3; i++) {
            int j = EDGE0(i);
            int k = EDGE1(i);

            if (!circumcenter(
                    this->x[INDEX3(this->nodes, t, j)],
                    this->y[INDEX3(this->nodes, t, j)],
                    this->x[INDEX3(this->nodes, t, k)],
                    this->y[INDEX3(this->nodes, t, k)],
                    targetx, targety,
                    INDEX2(c, i, 0), INDEX2(c, i, 1))) {

                // bail out with the appropriate values if we're actually on a 
                // node
                if ((fabs(targetx - this->x[INDEX3(this->nodes, t, j)]) < TOLERANCE_EPS)
                 && (fabs(targety - this->y[INDEX3(this->nodes, t, j)]) < TOLERANCE_EPS)) {
                    return z[INDEX3(this->nodes, t, j)];
                } else if ((fabs(targetx - this->x[INDEX3(this->nodes, t, k)]) < TOLERANCE_EPS)
                        && (fabs(targety - this->y[INDEX3(this->nodes, t, k)]) < TOLERANCE_EPS)) {
                    return z[INDEX3(this->nodes, t, k)];
                } else if (!onedge) {
                    onedge = true;
                    edge.push_back(INDEX3(this->nodes, t, j));
                    edge.push_back(INDEX3(this->nodes, t, k));
                    onhull = (INDEX3(neighbors, t, i) == -1);
                }
            }
        }
        for (int i=0; i<3; i++) {
            int j = EDGE0(i);
            int k = EDGE1(i);
            int q = INDEX3(this->nodes, t, i);
            double ati = 0.0;

            if (!onedge || ((edge[0] != q) && edge[1] != q)) {
                ati = signed_area(vx, vy, 
                                  INDEX2(c, j, 0), INDEX2(c, j, 1),
                                  INDEX2(c, k, 0), INDEX2(c, k, 1));


                yA = ati - cA;
                tA = A + yA;
                cA = (tA - A) - yA;
                A = tA;

                yf = ati*z[q] - cf;
                tf = f + yf;
                cf = (tf - f) - yf;
                f = tf;
            }
        }
    }

    // If we're on an edge, then the scheme of adding up triangles as above
    // doesn't work so well. We'll take care of these two nodes here.
    if (onedge) {

        // If we're on the convex hull, then the other nodes don't actually 
        // contribute anything, just the nodes for the edge we're on. The 
        // Voronoi "polygons" are infinite in extent.
        if (onhull) {
            double a = (hypot(targetx-x[edge[0]], targety-y[edge[0]]) / 
                        hypot(x[edge[1]]-x[edge[0]], y[edge[1]]-y[edge[0]]));
            return (1-a) * z[edge[0]] + a*z[edge[1]];
        }

        set<int> T(circumtri.begin(), circumtri.end());
        vector<int> newedges0; // the two nodes that edge[0] still connect to
        vector<int> newedges1; // the two nodes that edge[1] still connect to
        set<int> alltri0; // all of the circumtriangle edge[0] participates in
        set<int> alltri1; // all of the circumtriangle edge[1] participates in
        for (it = circumtri.begin(); it != circumtri.end(); it++) {
            for (int i=0; i<3; i++) {
                int ti = INDEX3(this->neighbors, *it, i);
                int j = EDGE0(i);
                int k = EDGE1(i);
                int q0 = INDEX3(this->nodes, *it, j);
                int q1 = INDEX3(this->nodes, *it, k);

                if ((q0 == edge[0]) || (q1 == edge[0])) alltri0.insert(*it);
                if ((q0 == edge[1]) || (q1 == edge[1])) alltri1.insert(*it);

                if (!T.count(ti)) {
                    // neighbor is not in the set of circumtriangles
                    if (q0 == edge[0]) newedges0.push_back(q1);
                    if (q1 == edge[0]) newedges0.push_back(q0);
                    if (q0 == edge[1]) newedges1.push_back(q1);
                    if (q1 == edge[1]) newedges1.push_back(q0);
                }
            }
        }

        set<int>::iterator sit;

        double cx, cy;
        ConvexPolygon poly0;
        ConvexPolygon poly1;

        if (edge[1] != newedges0[0]) {
            circumcenter(this->x[edge[0]], this->y[edge[0]],
                         this->x[newedges0[0]], this->y[newedges0[0]],
                         targetx, targety,
                         cx, cy);
            poly0.push(cx, cy);
        }
        if (edge[1] != newedges0[1]) {
            circumcenter(this->x[edge[0]], this->y[edge[0]],
                         this->x[newedges0[1]], this->y[newedges0[1]],
                         targetx, targety,
                         cx, cy);
            poly0.push(cx, cy);
        }

        if (edge[0] != newedges1[0]) {
            circumcenter(this->x[edge[1]], this->y[edge[1]],
                         this->x[newedges1[0]], this->y[newedges1[0]],
                         targetx, targety,
                         cx, cy);
            poly1.push(cx, cy);
        }
        if (edge[0] != newedges1[1]) {
            circumcenter(this->x[edge[1]], this->y[edge[1]],
                         this->x[newedges1[1]], this->y[newedges1[1]],
                         targetx, targety,
                         cx, cy);
            poly1.push(cx, cy);
        }

        for (sit = alltri0.begin(); sit != alltri0.end(); sit++) {
            poly0.push(INDEX2(this->centers, *sit, 0),
                       INDEX2(this->centers, *sit, 1)); 
        }
        for (sit = alltri1.begin(); sit != alltri1.end(); sit++) {
            poly1.push(INDEX2(this->centers, *sit, 0),
                       INDEX2(this->centers, *sit, 1)); 
        }

        double a0 = poly0.area();
        double a1 = poly1.area();


        f += a0*z[edge[0]];
        A += a0;
        f += a1*z[edge[1]];
        A += a1;

        // Anticlimactic, isn't it?
    }

    f /= A;
    return f;
}

void NaturalNeighbors::interpolate_grid(double *z, 
    double x0, double x1, int xsteps,
    double y0, double y1, int ysteps,
    double *output,
    double defvalue, int start_triangle)
{
    int i, ix, iy, rowtri, coltri, tri;
    double dx, dy, targetx, targety;

    dx = (x1 - x0) / (xsteps-1);
    dy = (y1 - y0) / (ysteps-1);

    rowtri = 0;
    i = 0;
    for (iy=0; iy<ysteps; iy++) {
        targety = y0 + dy*iy;
        rowtri = find_containing_triangle(x0, targety, rowtri);
        tri = rowtri;
        for (ix=0; ix<xsteps; ix++) {
            targetx = x0 + dx*ix;
            coltri = tri;
            INDEXN(output, xsteps, iy, ix) = interpolate_one(z, targetx, targety,
                defvalue, coltri);
            if (coltri != -1) tri = coltri;
        }
    }
}

void NaturalNeighbors::interpolate_unstructured(double *z, int size, 
    double *intx, double *inty, double *output, double defvalue)
{
    int i, tri1, tri2;

    tri1 = 0;
    tri2 = 0;
    for (i=0; i<size; i++) {
        tri2 = tri1;
        output[i] = interpolate_one(z, intx[i], inty[i], defvalue, tri2);
        if (tri2 != -1) tri1 = tri2;
    }
}