WSL/SLF GitLab Repository

ResamplingAlgorithms2D.cc 7.77 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
22
23
24
25
26
27
28
29
30
31
32
#include <meteoio/ResamplingAlgorithms2D.h>
#include <cmath>

using namespace std;

namespace mio {

/**
 * @brief Bilinear spatial data resampling
 */
const Grid2DObject ResamplingAlgorithms2D::BilinearResampling(const Grid2DObject &i_grid, const double &factor)
{
	const double cellsize = i_grid.cellsize/factor;
Mathias Bavay's avatar
Mathias Bavay committed
33
34
	const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.ncols)*factor ));
	const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
35
36
	Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);

37
38
39
40
41
42
43
	Bilinear(o_grid, i_grid); //GridObjects always keep nodata
	return o_grid;
}

const Grid2DObject ResamplingAlgorithms2D::cubicBSplineResampling(const Grid2DObject &i_grid, const double &factor)
{
	const double cellsize = i_grid.cellsize/factor;
Mathias Bavay's avatar
Mathias Bavay committed
44
45
	const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.ncols)*factor ));
	const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
46
47
48
	Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);

	cubicBSpline(o_grid, i_grid); //GridObjects always keep nodata
49
50
51
	return o_grid;
}

52
53
54
const Grid2DObject ResamplingAlgorithms2D::NearestNeighbour(const Grid2DObject &i_grid, const double &factor)
{
	const double cellsize = i_grid.cellsize/factor;
Mathias Bavay's avatar
Mathias Bavay committed
55
56
	const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.ncols)*factor ));
	const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
57
58
59
60
61
62
63
64
65
66
	Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);

	NearestNeighbour(o_grid, i_grid); //GridObjects always keep nodata
	return o_grid;
}

///////////////////////////////////////////////////////////////////////
//Private Methods
///////////////////////////////////////////////////////////////////////
void ResamplingAlgorithms2D::NearestNeighbour(Grid2DObject &o_grid, const Grid2DObject &i_grid)
67
{
Mathias Bavay's avatar
Mathias Bavay committed
68
69
	const size_t org_ncols = i_grid.ncols;
	const size_t org_nrows = i_grid.nrows;
70
71
72
	const double scale_x = (double)o_grid.ncols / (double)org_ncols;
	const double scale_y = (double)o_grid.nrows / (double)org_nrows;

Mathias Bavay's avatar
Mathias Bavay committed
73
74
	for (size_t jj=0; jj<o_grid.nrows; jj++) {
		size_t org_jj = (size_t) Optim::round( (double)jj/scale_y );
75
		if(org_jj>=org_nrows) org_jj=org_nrows-1;
76

Mathias Bavay's avatar
Mathias Bavay committed
77
78
		for (size_t ii=0; ii<o_grid.ncols; ii++) {
			size_t org_ii = (size_t) Optim::round( (double)ii/scale_x );
79
			if(org_ii>=org_ncols) org_ii=org_ncols-1;
80
81
82
83
84
			o_grid(ii,jj) = i_grid(org_ii, org_jj);
		}
	}
}

Mathias Bavay's avatar
Mathias Bavay committed
85
double ResamplingAlgorithms2D::bilinear_pixel(const Grid2DObject &i_grid, const size_t &org_ii, const size_t &org_jj, const size_t &org_ncols, const size_t &org_nrows, const double &x, const double &y)
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
{
	if(org_jj>=(org_nrows-1) || org_ii>=(org_ncols-1)) return i_grid(org_ii, org_jj);

	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;
	if(f_0_0!=IOUtils::nodata) {
		avg_value += f_0_0;
		avg_count++;
	}
	if(f_1_0!=IOUtils::nodata) {
		avg_value += f_1_0;
		avg_count++;
	}
	if(f_0_1!=IOUtils::nodata) {
		avg_value += f_0_1;
		avg_count++;
	}
	if(f_1_1!=IOUtils::nodata) {
		avg_value += f_1_1;
		avg_count++;
	}

	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;

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

	double value = 0.;
	const double avg = avg_value/(double)avg_count;
	if(f_0_0!=IOUtils::nodata) value += f_0_0 * (1.-x)*(1.-y);
	else value += avg * (1.-x)*(1.-y);
	if(f_1_0!=IOUtils::nodata) value += f_1_0 * x*(1.-y);
	else value += avg * x*(1.-y);
	if(f_0_1!=IOUtils::nodata) value += f_0_1 * (1.-x)*y;
	else value += avg * (1.-x)*y;
	if(f_1_1!=IOUtils::nodata) value += f_1_1 *x*y;
	else value += avg *x*y;

	return value;
}

void ResamplingAlgorithms2D::Bilinear(Grid2DObject &o_grid, const Grid2DObject &i_grid)
133
{
Mathias Bavay's avatar
Mathias Bavay committed
134
135
	const size_t org_ncols = i_grid.ncols;
	const size_t org_nrows = i_grid.nrows;
136
137
138
	const double scale_x = (double)o_grid.ncols / (double)org_ncols;
	const double scale_y = (double)o_grid.nrows / (double)org_nrows;

Mathias Bavay's avatar
Mathias Bavay committed
139
	for (size_t jj=0; jj<o_grid.nrows; jj++) {
140
		const double org_y = (double)jj/scale_y;
Mathias Bavay's avatar
Mathias Bavay committed
141
		const size_t org_jj = static_cast<size_t>( org_y );
142
143
		const double y = org_y - (double)org_jj; //normalized y, between 0 and 1

Mathias Bavay's avatar
Mathias Bavay committed
144
		for (size_t ii=0; ii<o_grid.ncols; ii++) {
145
			const double org_x = (double)ii/scale_x;
Mathias Bavay's avatar
Mathias Bavay committed
146
			const size_t org_ii = static_cast<size_t>( org_x );
147
148
			const double x = org_x - (double)org_ii; //normalized x, between 0 and 1

149
			o_grid(ii,jj) = bilinear_pixel(i_grid, org_ii, org_jj, org_ncols, org_nrows, x, y);
150
151
152
153
		}
	}
}

154
155
156
157
158
159
160
161
162
163
164
165
double ResamplingAlgorithms2D::BSpline_weight(const double &x) {
	double R = 0.;
	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.);

	return 1./6.*R;
}

void ResamplingAlgorithms2D::cubicBSpline(Grid2DObject &o_grid, const Grid2DObject &i_grid)
{//see http://paulbourke.net/texture_colour/imageprocess/
Mathias Bavay's avatar
Mathias Bavay committed
166
167
	const size_t org_ncols = i_grid.ncols;
	const size_t org_nrows = i_grid.nrows;
168
169
170
	const double scale_x = (double)o_grid.ncols / (double)org_ncols;
	const double scale_y = (double)o_grid.nrows / (double)org_nrows;

Mathias Bavay's avatar
Mathias Bavay committed
171
	for (size_t jj=0; jj<o_grid.nrows; jj++) {
172
		const double org_y = (double)jj/scale_y;
Mathias Bavay's avatar
Mathias Bavay committed
173
		const size_t org_jj = static_cast<size_t>( org_y );
174
		const double dy = org_y - (double)org_jj; //normalized y, between 0 and 1
175

Mathias Bavay's avatar
Mathias Bavay committed
176
		for (size_t ii=0; ii<o_grid.ncols; ii++) {
177
			const double org_x = (double)ii/scale_x;
Mathias Bavay's avatar
Mathias Bavay committed
178
			const size_t org_ii = static_cast<size_t>( org_x );
179
			const double dx = org_x - (double)org_ii; //normalized x, between 0 and 1
180

181
			double F = 0., max=-std::numeric_limits<double>::max(), min=std::numeric_limits<double>::max();
182
			unsigned int avg_count = 0;
183
184
			for(short n=-1; n<=2; n++) {
				for(short m=-1; m<=2; m++) {
185
					if(((signed)org_ii+m)<0 || ((signed)org_ii+m)>=(signed)org_ncols || ((signed)org_jj+n)<0 || ((signed)org_jj+n)>=(signed)org_nrows) continue;
186
					const double pixel = i_grid(org_ii+m, org_jj+n);
187
188
189
190
191
192
193
					if(pixel!=IOUtils::nodata) {
						F += pixel * BSpline_weight(m-dx) * BSpline_weight(dy-n);
						avg_count++;
						if(pixel>max) max=pixel;
						if(pixel<min) min=pixel;
					}
				}
194
195
			}

196
197
198
199
200
201
202
			if(avg_count==16) { //normal bicubic
				o_grid(ii,jj) = F*(16/avg_count);
				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
			else //not enought data points -> bilinear for this pixel
				o_grid(ii,jj) = bilinear_pixel(i_grid, org_ii, org_jj, org_ncols, org_nrows, dx, dy);
203

204
205
206
		}
	}

207
}
208
209
210

} //namespace