OGS
Raster.cpp
Go to the documentation of this file.
1
14#include "Raster.h"
15
16#include <fstream>
17
18// BaseLib
19#include "BaseLib/FileTools.h"
20#include "BaseLib/StringTools.h"
21#include "Triangle.h"
22
23namespace GeoLib
24{
25void Raster::refineRaster(std::size_t scaling)
26{
27 std::vector<double> new_raster_data(_header.n_rows * _header.n_cols *
28 scaling * scaling);
29
30 for (std::size_t row(0); row < _header.n_rows; row++)
31 {
32 for (std::size_t col(0); col < _header.n_cols; col++)
33 {
34 const std::size_t idx(row * _header.n_cols + col);
35 for (std::size_t new_row(row * scaling);
36 new_row < (row + 1) * scaling;
37 new_row++)
38 {
39 const std::size_t idx0(new_row * _header.n_cols * scaling);
40 for (std::size_t new_col(col * scaling);
41 new_col < (col + 1) * scaling;
42 new_col++)
43 {
44 new_raster_data[idx0 + new_col] = _raster_data[idx];
45 }
46 }
47 }
48 }
49
50 std::swap(_raster_data, new_raster_data);
51 _header.cell_size /= scaling;
52 _header.n_cols *= scaling;
53 _header.n_rows *= scaling;
54}
55
57{
58 if (pnt[0] >= _header.origin[0] &&
59 pnt[0] < (_header.origin[0] + (_header.cell_size * _header.n_cols)) &&
60 pnt[1] >= _header.origin[1] &&
61 pnt[1] < (_header.origin[1] + (_header.cell_size * _header.n_rows)))
62 {
63 auto cell_x = static_cast<int>(
64 std::floor((pnt[0] - _header.origin[0]) / _header.cell_size));
65 auto cell_y = static_cast<int>(
66 std::floor((pnt[1] - _header.origin[1]) / _header.cell_size));
67
68 // use raster boundary values if node is outside raster due to rounding
69 // errors or floating point arithmetic
70 cell_x = (cell_x < 0) ? 0
71 : ((cell_x > static_cast<int>(_header.n_cols))
72 ? static_cast<int>(_header.n_cols - 1)
73 : cell_x);
74 cell_y = (cell_y < 0) ? 0
75 : ((cell_y > static_cast<int>(_header.n_rows))
76 ? static_cast<int>(_header.n_rows - 1)
77 : cell_y);
78
79 const std::size_t index = cell_y * _header.n_cols + cell_x;
80 return _raster_data[index];
81 }
82 return _header.no_data;
83}
84
86{
87 // position in raster
88 double const xPos((pnt[0] - _header.origin[0]) / _header.cell_size);
89 double const yPos((pnt[1] - _header.origin[1]) / _header.cell_size);
90 // raster cell index
91 double const xIdx(std::floor(xPos)); // carry out computions in double
92 double const yIdx(std::floor(yPos)); // so not to over- or underflow.
93
94 // weights for bilinear interpolation
95 double const xShift = std::abs((xPos - xIdx) - 0.5);
96 double const yShift = std::abs((yPos - yIdx) - 0.5);
97 Eigen::Vector4d weight = {(1 - xShift) * (1 - yShift),
98 xShift * (1 - yShift), xShift * yShift,
99 (1 - xShift) * yShift};
100
101 // neighbors to include in interpolation
102 int const xShiftIdx = (xPos - xIdx >= 0.5) ? 1 : -1;
103 int const yShiftIdx = (yPos - yIdx >= 0.5) ? 1 : -1;
104 std::array<int, 4> const x_nb = {{0, xShiftIdx, xShiftIdx, 0}};
105 std::array<int, 4> const y_nb = {{0, 0, yShiftIdx, yShiftIdx}};
106
107 // get pixel values
108 Eigen::Vector4d pix_val{};
109 unsigned no_data_count(0);
110 for (unsigned j = 0; j < 4; ++j)
111 {
112 // check if neighbour pixel is still on the raster, otherwise substitute
113 // a no data value. This also allows the cast to unsigned type.
114 if ((xIdx + x_nb[j]) < 0 || (yIdx + y_nb[j]) < 0 ||
115 (xIdx + x_nb[j]) > (_header.n_cols - 1) ||
116 (yIdx + y_nb[j]) > (_header.n_rows - 1))
117 {
118 pix_val[j] = _header.no_data;
119 }
120 else
121 {
122 pix_val[j] = _raster_data[static_cast<std::size_t>(yIdx + y_nb[j]) *
124 static_cast<std::size_t>(xIdx + x_nb[j])];
125 }
126
127 // remove no data values
128 if (std::abs(pix_val[j] - _header.no_data) <
129 std::numeric_limits<double>::epsilon())
130 {
131 weight[j] = 0;
132 no_data_count++;
133 }
134 }
135
136 // adjust weights if necessary
137 if (no_data_count > 0)
138 {
139 if (no_data_count == 4)
140 { // if there is absolutely no data just use the default value
141 return _header.no_data;
142 }
143
144 weight /= weight.sum();
145 }
146
147 // new value
148 return weight.dot(pix_val);
149}
150
152{
153 return !(
154 (pnt[0] < _header.origin[0]) ||
155 (pnt[0] > _header.origin[0] + (_header.n_cols * _header.cell_size)) ||
156 (pnt[1] < _header.origin[1]) ||
157 (pnt[1] > _header.origin[1] + (_header.n_rows * _header.cell_size)));
158}
159
160} // end namespace GeoLib
Filename manipulation routines.
Definition of the GeoLib::Raster class.
Definition of string helper functions.
GeoLib::RasterHeader _header
Definition: Raster.h:150
std::vector< double > _raster_data
Definition: Raster.h:151
void refineRaster(std::size_t scaling)
Definition: Raster.cpp:25
double getValueAtPoint(const MathLib::Point3d &pnt) const
Definition: Raster.cpp:56
bool isPntOnRaster(MathLib::Point3d const &pnt) const
Definition: Raster.cpp:151
double interpolateValueAtPoint(const MathLib::Point3d &pnt) const
Definition: Raster.cpp:85
double cell_size
Definition: Raster.h:33
MathLib::Point3d origin
Definition: Raster.h:32
std::size_t n_cols
Definition: Raster.h:29
std::size_t n_rows
Definition: Raster.h:30