WSL/SLF GitLab Repository

ResamplingAlgorithms2D.cc 10.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/***********************************************************************************/
/*  Copyright 2011 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <meteoio/IOUtils.h>
19
#include <meteoio/MathOptim.h>
20
21
#include <meteoio/ResamplingAlgorithms2D.h>
#include <cmath>
22
#include <sstream>
23
#include <algorithm>
24
25
26
27
28

using namespace std;

namespace mio {

29
const Grid2DObject ResamplingAlgorithms2D::NearestNeighbour(const Grid2DObject &i_grid, const double &factor)
30
{
31
32
33
34
35
	if (factor<=0) {
		ostringstream ss;
		ss << "Rescaling factor " << factor << " is invalid!";
		throw InvalidArgumentException(ss.str(), AT);
	}
36
	const double cellsize = i_grid.cellsize/factor;
37
38
	const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor ));
	const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor ));
39
40
	Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);

41
	NearestNeighbour(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
42
43
44
	return o_grid;
}

45
46
47
48
/**
 * @brief Bilinear spatial data resampling
 */
const Grid2DObject ResamplingAlgorithms2D::BilinearResampling(const Grid2DObject &i_grid, const double &factor)
49
{
50
51
52
53
54
	if (factor<=0) {
		ostringstream ss;
		ss << "Rescaling factor " << factor << " is invalid!";
		throw InvalidArgumentException(ss.str(), AT);
	}
55
	const double cellsize = i_grid.cellsize/factor;
56
57
	const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor ));
	const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor ));
58
59
	Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);

60
	Bilinear(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
61
62
63
	return o_grid;
}

64
const Grid2DObject ResamplingAlgorithms2D::cubicBSplineResampling(const Grid2DObject &i_grid, const double &factor)
65
{
66
67
68
69
70
	if (factor<=0) {
		ostringstream ss;
		ss << "Rescaling factor " << factor << " is invalid!";
		throw InvalidArgumentException(ss.str(), AT);
	}
71
	const double cellsize = i_grid.cellsize/factor;
72
73
	const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor ));
	const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor ));
74
75
	Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);

76
	cubicBSpline(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
77
78
79
80
	return o_grid;
}

const Array2D<double> ResamplingAlgorithms2D::NearestNeighbour(const Array2D<double> &i_grid, const double &factor_x, const double &factor_y) {
81
82
83
84
85
	if (factor_x<=0 || factor_y<=0) {
		ostringstream ss;
		ss << "Rescaling factors (" << factor_x << "," << factor_y << ") are invalid!";
		throw InvalidArgumentException(ss.str(), AT);
	}
86
87
88
89
90
91
92
93
94
	const size_t nx = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor_x ));
	const size_t ny = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor_y ));
	Array2D<double> o_grid(nx, ny);

	NearestNeighbour(o_grid, i_grid);
	return o_grid;
}

const Array2D<double> ResamplingAlgorithms2D::BilinearResampling(const Array2D<double> &i_grid, const double &factor_x, const double &factor_y) {
95
96
97
98
99
	if (factor_x<=0 || factor_y<=0) {
		ostringstream ss;
		ss << "Rescaling factors (" << factor_x << "," << factor_y << ") are invalid!";
		throw InvalidArgumentException(ss.str(), AT);
	}
100
101
102
103
104
105
106
107
	const size_t nx = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor_x ));
	const size_t ny = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor_y ));
	Array2D<double> o_grid(nx, ny);
	Bilinear(o_grid, i_grid);
	return o_grid;
}

const Array2D<double> ResamplingAlgorithms2D::cubicBSplineResampling(const Array2D<double> &i_grid, const double &factor_x, const double &factor_y) {
108
109
110
111
112
	if (factor_x<=0 || factor_y<=0) {
		ostringstream ss;
		ss << "Rescaling factors (" << factor_x << "," << factor_y << ") are invalid!";
		throw InvalidArgumentException(ss.str(), AT);
	}
113
114
115
116
117
	const size_t nx = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor_x ));
	const size_t ny = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor_y ));
	Array2D<double> o_grid(nx, ny);

	cubicBSpline(o_grid, i_grid);
118
119
120
121
122
123
	return o_grid;
}

///////////////////////////////////////////////////////////////////////
//Private Methods
///////////////////////////////////////////////////////////////////////
124
void ResamplingAlgorithms2D::NearestNeighbour(Array2D<double> &o_grid, const Array2D<double> &i_grid)
125
{
126
127
128
129
	const size_t org_nx = i_grid.getNx(), org_ny = i_grid.getNy();
	const size_t dest_nx = o_grid.getNx(), dest_ny = o_grid.getNy();
	const double scale_x = (double)dest_nx / (double)org_nx;
	const double scale_y = (double)dest_ny / (double)org_ny;
130

131
	for (size_t jj=0; jj<dest_ny; jj++) {
132
		const size_t org_jj = std::min( (size_t) Optim::floor( (double)jj/scale_y ) , org_ny-1 );
133

134
		for (size_t ii=0; ii<dest_nx; ii++) {
135
			const size_t org_ii = std::min( (size_t) Optim::floor( (double)ii/scale_x ) , org_nx-1 );
136
137
138
139
140
			o_grid(ii,jj) = i_grid(org_ii, org_jj);
		}
	}
}

141
double ResamplingAlgorithms2D::bilinear_pixel(const Array2D<double> &i_grid, const size_t &org_ii, const size_t &org_jj, const size_t &org_nx, const size_t &org_ny, const double &x, const double &y)
142
{
143
	if (org_jj>=(org_ny-1) || org_ii>=(org_nx-1)) return i_grid(org_ii, org_jj);
144
145
146
147
148
149
150
151

	const double f_0_0 = i_grid(org_ii, org_jj);
	const double f_1_0 = i_grid(org_ii+1, org_jj);
	const double f_0_1 = i_grid(org_ii, org_jj+1);
	const double f_1_1 = i_grid(org_ii+1, org_jj+1);

	double avg_value = 0.;
	unsigned int avg_count = 0;
152
	if (f_0_0!=IOUtils::nodata) {
153
154
155
		avg_value += f_0_0;
		avg_count++;
	}
156
	if (f_1_0!=IOUtils::nodata) {
157
158
159
		avg_value += f_1_0;
		avg_count++;
	}
160
	if (f_0_1!=IOUtils::nodata) {
161
162
163
		avg_value += f_0_1;
		avg_count++;
	}
164
	if (f_1_1!=IOUtils::nodata) {
165
166
167
168
		avg_value += f_1_1;
		avg_count++;
	}

169
	if (avg_count==4) return f_0_0 * (1.-x)*(1.-y) + f_1_0 * x*(1.-y) + f_0_1 * (1.-x)*y + f_1_1 *x*y;
170
171

	//special cases: less than two neighbours or three neighbours
172
	if (avg_count<=2) return IOUtils::nodata;
173
174
175

	double value = 0.;
	const double avg = avg_value/(double)avg_count;
176
	if (f_0_0!=IOUtils::nodata) value += f_0_0 * (1.-x)*(1.-y);
177
	else value += avg * (1.-x)*(1.-y);
178
	if (f_1_0!=IOUtils::nodata) value += f_1_0 * x*(1.-y);
179
	else value += avg * x*(1.-y);
180
	if (f_0_1!=IOUtils::nodata) value += f_0_1 * (1.-x)*y;
181
	else value += avg * (1.-x)*y;
182
	if (f_1_1!=IOUtils::nodata) value += f_1_1 *x*y;
183
184
185
186
187
	else value += avg *x*y;

	return value;
}

188
void ResamplingAlgorithms2D::Bilinear(Array2D<double> &o_grid, const Array2D<double> &i_grid)
189
{
190
191
192
193
	const size_t org_nx = i_grid.getNx(), org_ny = i_grid.getNy();
	const size_t dest_nx = o_grid.getNx(), dest_ny = o_grid.getNy();
	const double scale_x = (double)dest_nx / (double)org_nx;
	const double scale_y = (double)dest_ny / (double)org_ny;
194

195
	for (size_t jj=0; jj<dest_ny; jj++) {
196
		const double org_y = (double)jj/scale_y;
Mathias Bavay's avatar
Mathias Bavay committed
197
		const size_t org_jj = static_cast<size_t>( org_y );
198
199
		const double y = org_y - (double)org_jj; //normalized y, between 0 and 1

200
		for (size_t ii=0; ii<dest_nx; ii++) {
201
			const double org_x = (double)ii/scale_x;
Mathias Bavay's avatar
Mathias Bavay committed
202
			const size_t org_ii = static_cast<size_t>( org_x );
203
204
			const double x = org_x - (double)org_ii; //normalized x, between 0 and 1

205
			o_grid(ii,jj) = bilinear_pixel(i_grid, org_ii, org_jj, org_nx, org_ny, x, y);
206
207
208
209
		}
	}
}

210
211
double ResamplingAlgorithms2D::BSpline_weight(const double &x)
{
212
	double R = 0.;
213
214
215
216
	if ((x+2.)>0.) R += Optim::pow3(x+2.);
	if ((x+1.)>0.) R += -4.*Optim::pow3(x+1.);
	if ((x)>0.) R += 6.*Optim::pow3(x);
	if ((x-1.)>0.) R += -4.*Optim::pow3(x-1.);
217
218
219
220

	return 1./6.*R;
}

221
void ResamplingAlgorithms2D::cubicBSpline(Array2D<double> &o_grid, const Array2D<double> &i_grid)
222
{//see http://paulbourke.net/texture_colour/imageprocess/
223
224
225
226
	const size_t org_nx = i_grid.getNx(), org_ny = i_grid.getNy();
	const size_t dest_nx = o_grid.getNx(), dest_ny = o_grid.getNy();
	const double scale_x = (double)dest_nx / (double)org_nx;
	const double scale_y = (double)dest_ny / (double)org_ny;
227

228
	for (size_t jj=0; jj<dest_ny; jj++) {
229
		const double org_y = (double)jj/scale_y;
Mathias Bavay's avatar
Mathias Bavay committed
230
		const size_t org_jj = static_cast<size_t>( org_y );
231
		const double dy = org_y - (double)org_jj; //normalized y, between 0 and 1
232

233
		for (size_t ii=0; ii<dest_nx; ii++) {
234
			const double org_x = (double)ii/scale_x;
Mathias Bavay's avatar
Mathias Bavay committed
235
			const size_t org_ii = static_cast<size_t>( org_x );
236
			const double dx = org_x - (double)org_ii; //normalized x, between 0 and 1
237

238
			double F = 0., max=-std::numeric_limits<double>::max(), min=std::numeric_limits<double>::max();
239
			unsigned int avg_count = 0;
240
241
			for (short n=-1; n<=2; n++) {
				for (short m=-1; m<=2; m++) {
242
					if (((signed)org_ii+m)<0 || ((signed)org_ii+m)>=(signed)org_nx || ((signed)org_jj+n)<0 || ((signed)org_jj+n)>=(signed)org_ny) continue;
243
					const double pixel = i_grid(org_ii+m, org_jj+n);
244
					if (pixel!=IOUtils::nodata) {
245
246
						F += pixel * BSpline_weight(m-dx) * BSpline_weight(dy-n);
						avg_count++;
247
248
						if (pixel>max) max=pixel;
						if (pixel<min) min=pixel;
249
250
					}
				}
251
252
			}

253
			if (avg_count==16) { //normal bicubic
254
				o_grid(ii,jj) = F;
255
256
257
				if (o_grid(ii,jj)>max) o_grid(ii,jj)=max; //try to limit overshoot
				else if (o_grid(ii,jj)<min) o_grid(ii,jj)=min; //try to limit overshoot
			} else if (avg_count==0) o_grid(ii,jj) = IOUtils::nodata; //nodata-> nodata
258
			else //not enought data points -> bilinear for this pixel
259
				o_grid(ii,jj) = bilinear_pixel(i_grid, org_ii, org_jj, org_nx, org_ny, dx, dy);
260
261
		}
	}
262
}
263
264
265

} //namespace