Add functions for nearest point and linear angle interpolation

fmidev · Feb 7, 2024 · d04edf5 · d04edf5
1 parent 22a938c
commit d04edf5
Show file tree

Hide file tree

Showing 2 changed files with 124 additions and 4 deletions.
diff --git a/himan-lib/include/numerical_functions.h b/himan-lib/include/numerical_functions.h
@@ -191,7 +191,7 @@ std::vector<T> Arange(T start, T stop, T step)
 
 	ret[0] = start;
 
-	std::generate(ret.begin() + 1, ret.end(), [ v = ret[0], &step ]() mutable { return v += step; });
+	std::generate(ret.begin() + 1, ret.end(), [v = ret[0], &step]() mutable { return v += step; });
 	return ret;
 }
 
@@ -235,7 +235,7 @@ std::vector<std::vector<T>> Linspace(const std::vector<T>& start, const std::vec
 
 		ret[i].resize(static_cast<size_t>(length));
 		ret[i][0] = start[i];
-		std::generate(ret[i].begin() + 1, ret[i].end(), [ v = ret[i][0], &step ]() mutable { return v += step; });
+		std::generate(ret[i].begin() + 1, ret[i].end(), [v = ret[i][0], &step]() mutable { return v += step; });
 	}
 
 	return ret;
@@ -278,10 +278,30 @@ namespace interpolation
  * Basic interpolation functions.
  */
 
+template <typename Type>
+CUDA_HOST CUDA_DEVICE inline Type NearestPoint(Type factor, Type Y1, Type Y2)
+{
+	if (factor < 0.5)
+	{
+		return Y1;
+	}
+	else
+	{
+		return Y2;
+	}
+}
+
+template <typename Type>
+CUDA_HOST CUDA_DEVICE inline Type NearestPoint(Type X, Type X1, Type X2, Type Y1, Type Y2)
+{
+	const Type factor = static_cast<Type>((X - X1) / (X2 - X1));
+	return NearestPoint<Type>(factor, Y1, Y2);
+}
+
 template <typename Type>
 CUDA_HOST CUDA_DEVICE inline Type Linear(Type factor, Type Y1, Type Y2)
 {
-	return static_cast<Type> (std::fma(factor, Y2, std::fma(-factor, Y1, Y1)));
+	return static_cast<Type>(std::fma(factor, Y2, std::fma(-factor, Y1, Y1)));
 }
 
 template <typename Type>
@@ -292,7 +312,7 @@ CUDA_HOST CUDA_DEVICE inline Type Linear(Type X, Type X1, Type X2, Type Y1, Type
 		return Y1;
 	}
 
-	const Type factor = static_cast<Type> ((X - X1) / (X2 - X1));
+	const Type factor = static_cast<Type>((X - X1) / (X2 - X1));
 	return Linear<Type>(factor, Y1, Y2);
 }
 
@@ -309,6 +329,30 @@ CUDA_HOST CUDA_DEVICE inline Type BiLinear(Type dx, Type dy, Type a, Type b, Typ
 	return (1 - dx) * (1 - dy) * c + dx * (1 - dy) * d + (1 - dx) * dy * a + dx * dy * b;
 }
 
+/**
+ * @brief Linear interpolation of angles, can deal with modulo 360
+ */
+
+template <typename Type>
+CUDA_HOST CUDA_DEVICE inline Type LinearAngle(Type factor, Type Y1, Type Y2)
+{
+	ASSERT(factor >= 0 && factor <= 1);
+
+	// cast to double so that can use fmod()
+	double y1 = static_cast<double>(Y1);
+	double y2 = static_cast<double>(Y2);
+
+	const double shortest_angle = fmod(fmod((y2 - y1), 360.) + 540., 360.) - 180.;
+
+	return static_cast<Type>(fmod(y1 + shortest_angle * factor + 360., 360.));
+}
+
+template <typename Type>
+CUDA_HOST CUDA_DEVICE inline Type LinearAngle(Type X, Type X1, Type X2, Type Y1, Type Y2)
+{
+	const Type factor = static_cast<Type>((X - X1) / (X2 - X1));
+	return LinearAngle<Type>(factor, Y1, Y2);
+}
 }  // namespace interpolation
 
 }  // namespace numerical_functions

diff --git a/himan-plugins/source/luatool.cpp b/himan-plugins/source/luatool.cpp
@@ -1179,6 +1179,79 @@ matrix<T> ProbLimitEq2D(const matrix<T>& A, const matrix<T>& B, T limit)
 }
 }  // namespace luabind_workaround
 
+namespace numerical_functions_wrapper
+{
+template <typename T>
+object Linear(const object& X, const object& X1, const object& X2, const object& Y1, const object& Y2)
+{
+	auto x = TableToVector<T>(X);
+	auto x1 = TableToVector<T>(X1);
+	auto x2 = TableToVector<T>(X2);
+	auto y1 = TableToVector<T>(Y1);
+	auto y2 = TableToVector<T>(Y2);
+
+	if (y1.size() != y2.size())
+	{
+		throw std::runtime_error("luatool::Linear Y1 and Y2 must be of same size");
+	}
+
+	std::vector<T> interp(y1.size());
+	for (size_t i = 0; i < y1.size(); i++)
+	{
+		interp[i] = numerical_functions::interpolation::Linear<T>(x[i], x1[i], x2[i], y1[i], y2[i]);
+	}
+
+	return VectorToTable<T>(interp);
+}
+
+template <typename T>
+object LinearAngle(const object& X, const object& X1, const object& X2, const object& Y1, const object& Y2)
+{
+	auto x = TableToVector<T>(X);
+	auto x1 = TableToVector<T>(X1);
+	auto x2 = TableToVector<T>(X2);
+	auto y1 = TableToVector<T>(Y1);
+	auto y2 = TableToVector<T>(Y2);
+
+	if (y1.size() != y2.size())
+	{
+		throw std::runtime_error("luatool::LinearAngle Y1 and Y2 must be of same size");
+	}
+
+	std::vector<T> interp(y1.size());
+	for (size_t i = 0; i < y1.size(); i++)
+	{
+		interp[i] = numerical_functions::interpolation::LinearAngle<T>(x[i], x1[i], x2[i], y1[i], y2[i]);
+	}
+
+	return VectorToTable<T>(interp);
+}
+
+template <typename T>
+object NearestPoint(const object& X, const object& X1, const object& X2, const object& Y1, const object& Y2)
+{
+	auto x = TableToVector<T>(X);
+	auto x1 = TableToVector<T>(X1);
+	auto x2 = TableToVector<T>(X2);
+	auto y1 = TableToVector<T>(Y1);
+	auto y2 = TableToVector<T>(Y2);
+
+	if (y1.size() != y2.size())
+	{
+		throw std::runtime_error("luatool::NearestPoint Y1 and Y2 must be of same size");
+	}
+
+	std::vector<T> interp(y1.size());
+	for (size_t i = 0; i < y1.size(); i++)
+	{
+		interp[i] = numerical_functions::interpolation::NearestPoint<T>(x[i], x1[i], x2[i], y1[i], y2[i]);
+	}
+
+	return VectorToTable<T>(interp);
+}
+
+}  // namespace numerical_functions_wrapper
+
 namespace params_wrapper
 {
 params create(const object& table)
@@ -1639,6 +1712,9 @@ void BindLib(lua_State* L)
 	          def("Max2D", &numerical_functions::Max2D<float>),
 	          def("Min2D", &numerical_functions::Min2D<double>),
 	          def("Min2D", &numerical_functions::Min2D<float>),
+		  def("Linear", &numerical_functions_wrapper::Linear<double>),
+		  def("LinearAngle", &numerical_functions_wrapper::LinearAngle<double>),
+		  def("NearestPoint", &numerical_functions_wrapper::NearestPoint<double>),
                   def("ProbLimitGt2D", &luabind_workaround::ProbLimitGt2D<double>),
                   def("ProbLimitGt2D", &luabind_workaround::ProbLimitGt2D<float>),
                   def("ProbLimitGe2D", &luabind_workaround::ProbLimitGe2D<double>),