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 auto* new_raster_data(
28 new double[_header.n_rows * _header.n_cols * 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 delete[] new_raster_data;
56}
57
59{
60 delete[] _raster_data;
61}
62
64{
65 if (pnt[0] >= _header.origin[0] &&
66 pnt[0] < (_header.origin[0] + (_header.cell_size * _header.n_cols)) &&
67 pnt[1] >= _header.origin[1] &&
68 pnt[1] < (_header.origin[1] + (_header.cell_size * _header.n_rows)))
69 {
70 auto cell_x = static_cast<int>(
71 std::floor((pnt[0] - _header.origin[0]) / _header.cell_size));
72 auto cell_y = static_cast<int>(
73 std::floor((pnt[1] - _header.origin[1]) / _header.cell_size));
74
75 // use raster boundary values if node is outside raster due to rounding
76 // errors or floating point arithmetic
77 cell_x = (cell_x < 0) ? 0
78 : ((cell_x > static_cast<int>(_header.n_cols))
79 ? static_cast<int>(_header.n_cols - 1)
80 : cell_x);
81 cell_y = (cell_y < 0) ? 0
82 : ((cell_y > static_cast<int>(_header.n_rows))
83 ? static_cast<int>(_header.n_rows - 1)
84 : cell_y);
85
86 const std::size_t index = cell_y * _header.n_cols + cell_x;
87 return _raster_data[index];
88 }
89 return _header.no_data;
90}
91
93{
94 // position in raster
95 double const xPos((pnt[0] - _header.origin[0]) / _header.cell_size);
96 double const yPos((pnt[1] - _header.origin[1]) / _header.cell_size);
97 // raster cell index
98 double const xIdx(std::floor(xPos)); // carry out computions in double
99 double const yIdx(std::floor(yPos)); // so not to over- or underflow.
100
101 // weights for bilinear interpolation
102 double const xShift = std::abs((xPos - xIdx) - 0.5);
103 double const yShift = std::abs((yPos - yIdx) - 0.5);
104 Eigen::Vector4d weight = {(1 - xShift) * (1 - yShift),
105 xShift * (1 - yShift), xShift * yShift,
106 (1 - xShift) * yShift};
107
108 // neighbors to include in interpolation
109 int const xShiftIdx = (xPos - xIdx >= 0.5) ? 1 : -1;
110 int const yShiftIdx = (yPos - yIdx >= 0.5) ? 1 : -1;
111 std::array<int, 4> const x_nb = {{0, xShiftIdx, xShiftIdx, 0}};
112 std::array<int, 4> const y_nb = {{0, 0, yShiftIdx, yShiftIdx}};
113
114 // get pixel values
115 Eigen::Vector4d pix_val{};
116 unsigned no_data_count(0);
117 for (unsigned j = 0; j < 4; ++j)
118 {
119 // check if neighbour pixel is still on the raster, otherwise substitute
120 // a no data value. This also allows the cast to unsigned type.
121 if ((xIdx + x_nb[j]) < 0 || (yIdx + y_nb[j]) < 0 ||
122 (xIdx + x_nb[j]) > (_header.n_cols - 1) ||
123 (yIdx + y_nb[j]) > (_header.n_rows - 1))
124 {
125 pix_val[j] = _header.no_data;
126 }
127 else
128 {
129 pix_val[j] = _raster_data[static_cast<std::size_t>(yIdx + y_nb[j]) *
131 static_cast<std::size_t>(xIdx + x_nb[j])];
132 }
133
134 // remove no data values
135 if (std::abs(pix_val[j] - _header.no_data) <
136 std::numeric_limits<double>::epsilon())
137 {
138 weight[j] = 0;
139 no_data_count++;
140 }
141 }
142
143 // adjust weights if necessary
144 if (no_data_count > 0)
145 {
146 if (no_data_count == 4)
147 { // if there is absolutely no data just use the default value
148 return _header.no_data;
149 }
150
151 weight /= weight.sum();
152 }
153
154 // new value
155 return weight.dot(pix_val);
156}
157
159{
160 return !(
161 (pnt[0] < _header.origin[0]) ||
162 (pnt[0] > _header.origin[0] + (_header.n_cols * _header.cell_size)) ||
163 (pnt[1] < _header.origin[1]) ||
164 (pnt[1] > _header.origin[1] + (_header.n_rows * _header.cell_size)));
165}
166
167} // end namespace GeoLib
Filename manipulation routines.
Definition of the GeoLib::Raster class.
Definition of string helper functions.
GeoLib::RasterHeader _header
Definition: Raster.h:137
void refineRaster(std::size_t scaling)
Definition: Raster.cpp:25
double * _raster_data
Definition: Raster.h:138
double getValueAtPoint(const MathLib::Point3d &pnt) const
Definition: Raster.cpp:63
bool isPntOnRaster(MathLib::Point3d const &pnt) const
Checks if the given point is located within the (x,y)-extension of the raster.
Definition: Raster.cpp:158
double interpolateValueAtPoint(const MathLib::Point3d &pnt) const
Interpolates the elevation of the given point based on the 8-neighbourhood of the raster cell it is l...
Definition: Raster.cpp:92
double cell_size
Definition: Raster.h:32
MathLib::Point3d origin
Definition: Raster.h:31
std::size_t n_cols
Definition: Raster.h:28
std::size_t n_rows
Definition: Raster.h:29