1 /* Robert Hijmans, June 2011, July 2016
2 // Based on public-domain code by Darel Rex Finley, 2007
3 // http://alienryderflex.com/polygon_fill/
4
5 */
6 #include <Rcpp.h>
7 using namespace Rcpp;
8 using namespace std;
9 #include <vector>
10 #include "spat.h"
11
rasterize_polygon(std::vector<double> r,double value,std::vector<double> pX,std::vector<double> pY,unsigned nrows,unsigned ncols,double xmin,double ymax,double rx,double ry)12 std::vector<double> rasterize_polygon(std::vector<double> r, double value, std::vector<double> pX, std::vector<double> pY, unsigned nrows, unsigned ncols, double xmin, double ymax, double rx, double ry) {
13
14 unsigned n = pX.size();
15 std::vector<unsigned> nCol(n);
16
17 for (size_t row=0; row<nrows; row++) {
18
19 double y = ymax - (row+0.5) * ry;
20 // Get nodes.
21 unsigned nodes = 0;
22 size_t j = n-1;
23 for (size_t i=0; i<n; i++) {
24 if (((pY[i] < y) && (pY[j] >= y)) || ((pY[j] < y) && (pY[i] >= y))) {
25 double nds = ((pX[i] - xmin + (y-pY[i])/(pY[j]-pY[i]) * (pX[j]-pX[i])) + 0.5 * rx ) / rx;
26 nds = nds < 0 ? 0 : nds;
27 nds = nds > ncols ? ncols : nds;
28 nCol[nodes] = (unsigned) nds;
29 nodes++;
30 }
31 j = i;
32 }
33
34 std::sort(nCol.begin(), nCol.begin()+nodes);
35 unsigned ncell = ncols * row;
36
37 // Fill the cells between node pairs.
38 for (size_t i=0; i < nodes; i+=2) {
39 if (nCol[i+1] > 0 && nCol[i] < ncols) {
40 //if (nCol[i] >= ncols || nCol[i+1] <= 0) break;
41 for (size_t col = nCol[i]; col < nCol[i+1]; col++) {
42 r[col + ncell] = value;
43 }
44 }
45 }
46 }
47 return(r);
48 }
49
50
51
52
rasterize(unsigned nrow,unsigned ncol,std::vector<double> extent,std::vector<double> values,double background)53 std::vector<double> SpPolygons::rasterize(unsigned nrow, unsigned ncol, std::vector<double> extent, std::vector<double> values, double background) {
54
55 unsigned n = size();
56
57 std::vector<double> v(nrow*ncol, background);
58
59 double resx = (extent[1] - extent[0]) / ncol;
60 double resy = (extent[3] - extent[2]) / nrow;
61
62 for (size_t j = 0; j < n; j++) {
63
64 SpPoly poly = getPoly(j);
65 double value = values[j];
66 unsigned np = poly.size();
67 for (size_t k = 0; k < np; k++) {
68 SpPolyPart part = poly.getPart(k);
69
70 if (part.hasHoles()) {
71 std::vector<double> vv = rasterize_polygon(v, value, part.x, part.y, nrow, ncol, extent[0], extent[3], resx, resy);
72 for (size_t h=0; h < part.nHoles(); h++) {
73 vv = rasterize_polygon(vv, background, part.xHole[h], part.yHole[h], nrow, ncol, extent[0], extent[3], resx, resy);
74 }
75 for (size_t q=0; q < vv.size(); q++) {
76 if ((vv[q] != background) && (!std::isnan(vv[q]))) {
77 //if (vv[q] != background) {
78 v[q] = vv[q];
79 }
80 }
81 } else {
82 v = rasterize_polygon(v, value, part.x, part.y, nrow, ncol, extent[0], extent[3], resx, resy);
83 }
84 }
85 }
86
87 return(v);
88
89 }
90
91
92