python3-reproject 0.3.1-4 / usr / lib / python3 / dist-packages / reproject / spherical_intersect / overlapAreaPP.c

/* Methods to compute pixel overlap areas in the plane.
 *
 * Originally developed in 2003 / 2004 by John Good.
 */

#include <math.h>

double computeOverlapPP(double *ix, double *iy, double minX, double maxX,
                        double minY, double maxY, double pixelArea);

double polyArea(int npts, double *nx, double *ny);

int rectClip(int n, double *x, double *y, double *nx, double *ny, double minX,
             double minY, double maxX, double maxY);
int lineClip(int n, double *x, double *y, double *nx, double *ny, double val,
             int dir);
int inPlane(double test, double divider, int direction);
int ptInPoly(double x, double y, int n, double *xp, double *yp);

// Global variables
double tmpX0[100];
double tmpX1[100];
double tmpY0[100];
double tmpY1[100];

/*
 * Sets up the polygons, runs the overlap
 * computation, and returns the area of overlap.
 * This version works in pixel space rather than
 * on the celestial sphere.
 */
double computeOverlapPP(double *ix, double *iy, double minX, double maxX,
                        double minY, double maxY, double pixelArea) {
  int npts;
  double area;

  double nx[100];
  double ny[100];

  double xp[4], yp[4];

  // Clip the input pixel polygon with the output pixel range

  npts = rectClip(4, ix, iy, nx, ny, minX, minY, maxX, maxY);

  // If no points, it may mean that the output is
  // completely contained in the input

  if (npts < 3) {
    xp[0] = minX;
    yp[0] = minY;
    xp[1] = maxX;
    yp[1] = minY;
    xp[2] = maxX;
    yp[2] = maxY;
    xp[3] = minX;
    yp[3] = maxY;

    if (ptInPoly(ix[0], iy[0], 4, xp, yp)) {
      area = pixelArea;
      return area;
    }

    return 0.;
  }

  area = polyArea(npts, nx, ny) * pixelArea;

  return area;
}

int rectClip(int n, double *x, double *y, double *nx, double *ny, double minX,
             double minY, double maxX, double maxY) {
  int nCurr;

  nCurr = lineClip(n, x, y, tmpX0, tmpY0, minX, 1);

  if (nCurr > 0) {
    nCurr = lineClip(nCurr, tmpX0, tmpY0, tmpX1, tmpY1, maxX, 0);

    if (nCurr > 0) {
      nCurr = lineClip(nCurr, tmpY1, tmpX1, tmpY0, tmpX0, minY, 1);

      if (nCurr > 0) {
        nCurr = lineClip(nCurr, tmpY0, tmpX0, ny, nx, maxY, 0);
      }
    }
  }

  return nCurr;
}

int lineClip(int n, double *x, double *y, double *nx, double *ny, double val,
             int dir) {
  int i;
  int nout;
  int last;

  double ycross;

  nout = 0;
  last = inPlane(x[n - 1], val, dir);

  for (i = 0; i < n; ++i) {
    if (last) {
      if (inPlane(x[i], val, dir)) {
        // Both endpoints in, just add the new point

        nx[nout] = x[i];
        ny[nout] = y[i];

        ++nout;
      } else {
        // Moved out of the clip region, add the point we moved out

        if (i == 0)
          ycross = y[n - 1]
              + (y[0] - y[n - 1]) * (val - x[n - 1]) / (x[0] - x[n - 1]);
        else
          ycross = y[i - 1]
              + (y[i] - y[i - 1]) * (val - x[i - 1]) / (x[i] - x[i - 1]);

        nx[nout] = val;
        ny[nout] = ycross;

        ++nout;

        last = 0;
      }
    } else {
      if (inPlane(x[i], val, dir)) {
        // Moved into the clip region.
        // Add the point we moved in, and the end point.

        if (i == 0)
          ycross = y[n - 1]
              + (y[0] - y[n - 1]) * (val - x[n - 1]) / (x[i] - x[n - 1]);
        else
          ycross = y[i - 1]
              + (y[i] - y[i - 1]) * (val - x[i - 1]) / (x[i] - x[i - 1]);

        nx[nout] = val;
        ny[nout] = ycross;

        ++nout;

        nx[nout] = x[i];
        ny[nout] = y[i];

        ++nout;

        last = 1;
      } else {
        // Segment entirely clipped.
      }
    }
  }

  return nout;
}

int inPlane(double test, double divider, int direction) {
  if (direction)
    return test >= divider;
  else
    return test <= divider;
}

double polyArea(int npts, double *nx, double *ny) {
  int i, inext;
  double area;

  area = 0.;

  for (i = 0; i < npts; ++i) {
    inext = (i + 1) % npts;

    area += nx[i] * ny[inext] - nx[inext] * ny[i];
  }

  area = fabs(area) / 2;

  return area;
}

int ptInPoly(double x, double y, int n, double *xp, double *yp) {
  int i, inext, count;
  double t;

  count = 0;

  for (i = 0; i < n; ++i) {
    inext = (i + 1) % n;

    if (((yp[i] <= y) && (yp[inext] > y))
        || ((yp[i] > y) && (yp[inext] <= y))) {
      t = (y - yp[i]) / (yp[inext] - yp[i]);

      if (x < xp[i] + t * (xp[inext] - xp[i]))
        ++count;
    }
  }

  return (count & 1);
}