diff --git a/BaseLib/quicksort.h b/BaseLib/quicksort.h
index b887b8c10fa04c957e990b34aae478fed9b5c8fe..bccebb1fb645ed7e3057c4c1b9c520e2a510010d 100644
--- a/BaseLib/quicksort.h
+++ b/BaseLib/quicksort.h
@@ -16,17 +16,17 @@
#define QUICKSORT_H_
// STL
-#include <cstddef>
#include <algorithm>
+#include <cstddef>
-namespace BaseLib {
-
+namespace BaseLib
+{
template <class T>
unsigned partition_(T* array, unsigned beg, unsigned end)
{
- unsigned i = beg+1;
- unsigned j = end-1;
- T m = array[beg];
+ unsigned i = beg + 1;
+ unsigned j = end - 1;
+ T m = array[beg];
for (;;) {
while ((i<end) && (array[i] < m)) i++;
@@ -36,8 +36,8 @@ unsigned partition_(T* array, unsigned beg, unsigned end)
std::swap(array[i], array[j]);
}
- std::swap(array[beg], array[j]);
- return j;
+ std::swap(array[beg], array[j]);
+ return j;
}
template <class T>
@@ -62,7 +62,7 @@ void quickSort(T* array, unsigned beg, unsigned end)
* @return
*/
template <typename T1, typename T2>
-std::size_t partition_(T1* array, std::size_t beg, std::size_t end, T2 *second_array)
+std::size_t partition_(T1 *array, std::size_t beg, std::size_t end, T2 *second_array)
{
std::size_t i = beg + 1;
std::size_t j = end - 1;
@@ -74,7 +74,8 @@ std::size_t partition_(T1* array, std::size_t beg, std::size_t end, T2 *second_a
while ((j > beg) && !(array[j] <= m))
j--;
- if (i >= j) break;
+ if (i >= j)
+ break;
std::swap(array[i], array[j]);
std::swap(second_array[i], second_array[j]);
@@ -100,7 +101,7 @@ void quicksort(T1* array, std::size_t beg, std::size_t end, T2* second_array)
if (beg < end) {
std::size_t p = partition_(array, beg, end, second_array);
quicksort(array, beg, p, second_array);
- quicksort(array, p+1, end, second_array);
+ quicksort(array, p + 1, end, second_array);
}
}
@@ -109,21 +110,25 @@ void quicksort(T1* array, std::size_t beg, std::size_t end, T2* second_array)
// STL
#include <vector>
-namespace BaseLib {
-
-template <typename T>
-class Quicksort {
+namespace BaseLib
+{
+template <typename T1, typename T2 = std::size_t>
+class Quicksort
+{
public:
- Quicksort (std::vector<T>& array, std::size_t beg, std::size_t end, std::vector<std::size_t>& perm)
+ Quicksort (std::vector<T1>& array, std::size_t beg, std::size_t end, std::vector<T2>& perm)
{
quicksort (array, beg, end, perm);
}
private:
- std::size_t partition_(std::vector<T>& array, std::size_t beg, std::size_t end, std::vector<std::size_t>& perm)
+ std::size_t partition_(std::vector<T1>& array,
+ std::size_t beg,
+ std::size_t end,
+ std::vector<T2>& perm)
{
std::size_t i = beg + 1;
std::size_t j = end - 1;
- T m = array[beg];
+ T1 m = array[beg];
for (;;) {
while ((i < end) && (array[i] <= m))
@@ -142,36 +147,44 @@ private:
return j;
}
- void quicksort(std::vector<T>& array, std::size_t beg, std::size_t end, std::vector<std::size_t>& perm)
+ void quicksort(std::vector<T1>& array,
+ std::size_t beg,
+ std::size_t end,
+ std::vector<T2>& perm)
{
if (beg < end) {
std::size_t p = partition_(array, beg, end, perm);
quicksort(array, beg, p, perm);
- quicksort(array, p+1, end, perm);
+ quicksort(array, p + 1, end, perm);
}
}
};
// specialization for pointer types
-template <typename T>
-class Quicksort <T *> {
+template <typename T1, typename T2>
+class Quicksort <T1*, T2>
+{
public:
- Quicksort (std::vector<T*>& array, std::size_t beg, std::size_t end, std::vector<std::size_t>& perm)
+ Quicksort (std::vector<T1*>& array, std::size_t beg, std::size_t end, std::vector<T2>& perm)
{
quicksort (array, beg, end, perm);
}
- Quicksort (std::vector<std::size_t>& perm, std::size_t beg, std::size_t end, std::vector<T*>& array)
+ Quicksort (std::vector<std::size_t>& perm, std::size_t beg, std::size_t end,
+ std::vector<T1*>& array)
{
quicksort (perm, beg, end, array);
}
private:
- std::size_t partition_(std::vector<T*>& array, std::size_t beg, std::size_t end, std::vector<std::size_t>& perm)
+ std::size_t partition_(std::vector<T1*>& array,
+ std::size_t beg,
+ std::size_t end,
+ std::vector<T2>& perm)
{
std::size_t i = beg + 1;
std::size_t j = end - 1;
- T* m = array[beg];
+ T1* m = array[beg];
for (;;) {
while ((i < end) && (*array[i] <= *m))
@@ -190,16 +203,22 @@ private:
return j;
}
- void quicksort(std::vector<T*>& array, std::size_t beg, std::size_t end, std::vector<std::size_t>& perm)
+ void quicksort(std::vector<T1*>& array,
+ std::size_t beg,
+ std::size_t end,
+ std::vector<T2>& perm)
{
if (beg < end) {
std::size_t p = partition_(array, beg, end, perm);
quicksort(array, beg, p, perm);
- quicksort(array, p+1, end, perm);
+ quicksort(array, p + 1, end, perm);
}
}
- std::size_t partition_(std::vector<std::size_t> &perm, std::size_t beg, std::size_t end, std::vector<T*>& array)
+ std::size_t partition_(std::vector<std::size_t> &perm,
+ std::size_t beg,
+ std::size_t end,
+ std::vector<T1*>& array)
{
std::size_t i = beg + 1;
std::size_t j = end - 1;
@@ -222,16 +241,18 @@ private:
return j;
}
- void quicksort(std::vector<std::size_t>& perm, std::size_t beg, std::size_t end, std::vector<T*>& array)
+ void quicksort(std::vector<std::size_t>& perm,
+ std::size_t beg,
+ std::size_t end,
+ std::vector<T1*>& array)
{
if (beg < end) {
std::size_t p = partition_(perm, beg, end, array);
quicksort(perm, beg, p, array);
- quicksort(perm, p+1, end, array);
+ quicksort(perm, p + 1, end, array);
}
}
};
-
} // end namespace BaseLib
#endif /* QUICKSORT_H_ */
diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp
index d260dc1aa77a21bfabfb7e2fc7048797f5a31bc9..67ab69db48459971998af6d0ed4c7a263f918f66 100644
--- a/GeoLib/PointVec.cpp
+++ b/GeoLib/PointVec.cpp
@@ -30,8 +30,9 @@ namespace GeoLib
{
PointVec::PointVec (const std::string& name, std::vector<Point*>* points,
std::map<std::string, size_t>* name_id_map, PointType type, double rel_eps) :
- TemplateVec<Point> (name, points, name_id_map),
- _type(type), _sqr_shortest_dist (std::numeric_limits<double>::max()), _aabb(points->begin(), points->end())
+ TemplateVec<Point> (name, points, name_id_map),
+ _type(type), _sqr_shortest_dist (std::numeric_limits<double>::max()),
+ _aabb(points->begin(), points->end())
{
assert (_data_vec);
size_t number_of_all_input_pnts (_data_vec->size());
@@ -174,7 +175,7 @@ void PointVec::makePntsUnique (std::vector<GeoLib::Point*>* pnt_vec,
pnt_id_map[perm[k + 1]] = pnt_id_map[perm[k]];
// reverse permutation
- BaseLib::Quicksort<GeoLib::Point*> (perm, 0, n_pnts_in_file, *pnt_vec);
+ BaseLib::Quicksort<std::size_t, GeoLib::Point*> (perm, 0, n_pnts_in_file, *pnt_vec);
// remove the second, third, ... occurrence from vector
for (size_t k(0); k < n_pnts_in_file; k++)
diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
index 8e75dcdbf0cce305ebc1363dd331a23638df9819..4099dda882f48a229cbc92349ff1ba570d0a62ae 100644
--- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
+++ b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
@@ -11,57 +11,72 @@
* http://www.opengeosys.org/project/license
*
*/
+#include <cmath>
+#include <limits>
#include "PiecewiseLinearInterpolation.h"
-#include <iostream>
+// BaseLib
+#include "quicksort.h"
namespace MathLib
{
PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(const std::vector<double>& supporting_points,
- const std::vector<double>& values_at_supp_pnts)
- : _supporting_points (supporting_points), _values_at_supp_pnts (values_at_supp_pnts)
-{}
+ const std::vector<double>& values_at_supp_pnts,
+ bool supp_pnts_sorted) :
+ _supp_pnts(supporting_points), _values_at_supp_pnts(values_at_supp_pnts)
+{
+ if (!supp_pnts_sorted) {
+ BaseLib::Quicksort<double, double>(_supp_pnts, static_cast<std::size_t> (0),
+ _supp_pnts.size(), _values_at_supp_pnts);
+ }
+}
PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(const std::vector<double>& supporting_points,
- const std::vector<double>& values_at_supp_pnts,
- const std::vector<double>& points_to_interpolate,
- std::vector<double>& values_at_interpol_pnts)
- : _supporting_points (supporting_points), _values_at_supp_pnts (values_at_supp_pnts)
+ const std::vector<double>& values_at_supp_pnts,
+ const std::vector<double>& points_to_interpolate,
+ std::vector<double>& values_at_interpol_pnts,
+ bool supp_pnts_sorted) :
+ _supp_pnts(supporting_points), _values_at_supp_pnts(values_at_supp_pnts)
{
-// std::cout << "PiecewiseLinearInterpolation::PiecewiseLinearInterpolation support_points, values_at_supp_pnts: " << std::endl;
-// for (size_t k(0); k<supporting_points.size(); k++) {
-// std::cout << supporting_points[k] << " " << values_at_supp_pnts[k] << std::endl;
-// }
-// std::cout << std::endl;
- values_at_interpol_pnts.clear();
- for (size_t k(0); k < points_to_interpolate.size(); k++)
- values_at_interpol_pnts.push_back (this->getValue (points_to_interpolate[k]));
+ if (!supp_pnts_sorted) {
+ BaseLib::Quicksort<double, double>(_supp_pnts, static_cast<std::size_t> (0),
+ _supp_pnts.size(),
+ _values_at_supp_pnts);
+ }
+
+ values_at_interpol_pnts.resize(points_to_interpolate.size());
+ std::transform(points_to_interpolate.begin(), points_to_interpolate.end(),
+ values_at_interpol_pnts.begin(),
+ [&](double const& p) { return this->getValue(p); } );
}
PiecewiseLinearInterpolation::~PiecewiseLinearInterpolation()
{}
-double PiecewiseLinearInterpolation::getValue ( double pnt_to_interpolate )
+double PiecewiseLinearInterpolation::getValue(double pnt_to_interpolate) const
{
// search interval that has the point inside
- size_t interval_idx (std::numeric_limits<size_t>::max());
- for (size_t k(1);
- k < _supporting_points.size() && interval_idx == std::numeric_limits<size_t>::max();
- k++)
- if (_supporting_points[k - 1] <= pnt_to_interpolate && pnt_to_interpolate <=
- _supporting_points[k])
+ std::size_t interval_idx(std::numeric_limits<std::size_t>::max());
+ for (std::size_t k(1); k < _supp_pnts.size()
+ && interval_idx == std::numeric_limits<std::size_t>::max(); k++)
+ if (_supp_pnts[k - 1] <= pnt_to_interpolate
+ && pnt_to_interpolate <= _supp_pnts[k])
interval_idx = k - 1;
+ if (interval_idx == std::numeric_limits<std::size_t>::max()) {
+ const double dist_first(fabs(_supp_pnts[0] - pnt_to_interpolate));
+ const double dist_last(fabs(_supp_pnts[_supp_pnts.size() - 1] - pnt_to_interpolate));
+ if (dist_first < dist_last)
+ interval_idx = 0;
+ else
+ interval_idx = _supp_pnts.size() - 2;
+ }
+
// compute linear interpolation polynom: y = m * x + n
- long double m (
- (_values_at_supp_pnts[interval_idx +
- 1] -
- _values_at_supp_pnts[interval_idx]) /
- (_supporting_points[interval_idx + 1] - _supporting_points[interval_idx]));
-// double m ((_values_at_supp_pnts[interval_idx] - _values_at_supp_pnts[interval_idx+1]) / (_supporting_points[interval_idx] - _supporting_points[interval_idx+1]));
-// double n (_values_at_supp_pnts[interval_idx+1] - m * _supporting_points[interval_idx+1]);
- long double n (_values_at_supp_pnts[interval_idx] - m * _supporting_points[interval_idx]);
+ long double m((_values_at_supp_pnts[interval_idx + 1] - _values_at_supp_pnts[interval_idx])
+ / (_supp_pnts[interval_idx + 1] - _supp_pnts[interval_idx]));
+ long double n(_values_at_supp_pnts[interval_idx] - m * _supp_pnts[interval_idx]);
return m * pnt_to_interpolate + n;
}
diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h
index 975f3582fc6a8ee705a75e8a3c48079013bb486e..b5596622d5602efe6428cd442be6155d6bd6fa9d 100644
--- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h
+++ b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h
@@ -15,27 +15,71 @@
#ifndef PIECEWISELINEARINTERPOLATION_H_
#define PIECEWISELINEARINTERPOLATION_H_
-#include <limits>
#include <vector>
namespace MathLib
{
+/**
+ * This class implements a one dimensional piecewise linear interpolation algorithm.
+ */
class PiecewiseLinearInterpolation
{
public:
+ /**
+ * The constructor copies the entries of the vector of supporting points
+ * \f$(x_0, x_1, \dots, x_n)\f$ and the entries of the vector of values at
+ * the supporting points \f$(y_0, y_1, \dots, y_n)\f$ where \f$n\f$
+ * is the number of entries of the vector. The number of supporting
+ * points must be the same like the number of values at the supporting
+ * points. It is assumed that \f$x_j\f$ corresponds to
+ * \f$y_j\f$ for all \f$j \in [0, n]\f$.
+ *
+ * It is not assumed that the supporting points are sorted, i.e.
+ * \f$x_0 < x_1 < \dots < x_n\f$. It is assumed, that the supporting points
+ * are pairwise different. The user can set the flag supp_pnts_sorted to
+ * true, if the supporting points are sorted. This will save some setup
+ * time.
+ * @param supporting_points vector of supporting points
+ * @param values_at_supp_pnts vector of values at the supporting points
+ * @param supp_pnts_sorted false (default), if it is sure the supporting points are sorted
+ * one can set the switch to true
+ */
PiecewiseLinearInterpolation(const std::vector<double>& supporting_points,
- const std::vector<double>& values_at_supp_pnts);
+ const std::vector<double>& values_at_supp_pnts,
+ bool supp_pnts_sorted = false);
+ /**
+ * The same requirements on the input vectors supporting_points and values_at_supp_pnts
+ * have to apply as for the previous constructor.
+ * @param supporting_points vector of supporting points
+ * @param values_at_supp_pnts vector of values at the supporting points
+ * @param points_to_interpolate The points should be located within the range
+ * \f$[x_{\min}, x_{\max}]\f$, where \f$x_{\min} = \min_{1 \le j \le n} x_j\f$ and
+ * \f$x_{\max} = \max_{1 \le j \le n} x_j\f$. Points outside of this interval are extrapolated.
+ * @param values_at_interpol_pnts The interpolated values corresponding to the entries
+ * of the vector points_to_interpolate.
+ * @param supp_pnts_sorted Is set to false per default. If it is sure that the supporting
+ * points are sorted one can set the switch to true.
+ */
PiecewiseLinearInterpolation(const std::vector<double>& supporting_points,
- const std::vector<double>& values_at_supp_pnts,
- const std::vector<double>& points_to_interpolate,
- std::vector<double>& values_at_interpol_pnts);
+ const std::vector<double>& values_at_supp_pnts,
+ const std::vector<double>& points_to_interpolate,
+ std::vector<double>& values_at_interpol_pnts,
+ bool supp_pnts_sorted = false);
virtual ~PiecewiseLinearInterpolation();
- double getValue ( double pnt_to_interpolate );
+ /**
+ * \brief Calculates the interpolation value.
+ * @param pnt_to_interpolate The point should be located within the range
+ * \f$[x_{\min}, x_{\max}]\f$, where \f$x_{\min} = \min_{1 \le j \le n} x_j\f$ and
+ * \f$x_{\max} = \max_{1 \le j \le n} x_j\f$. Points outside of this interval are
+ * extrapolated.
+ * @return The interpolated value.
+ */
+ double getValue(double pnt_to_interpolate) const;
private:
- const std::vector<double>& _supporting_points;
- const std::vector<double>& _values_at_supp_pnts;
+ std::vector<double> _supp_pnts;
+ std::vector<double> _values_at_supp_pnts;
};
} // end namespace MathLib
diff --git a/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp b/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1772e8860b5b75bf86593964562fba8639427043
--- /dev/null
+++ b/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp
@@ -0,0 +1,90 @@
+/**
+ * @file TestPiecewiseLinearInterpolation.cpp
+ * @author Thomas Fischer
+ * @date Feb 12, 2013
+ * @brief
+ *
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+// stl
+#include <limits>
+
+// google test
+#include "gtest/gtest.h"
+
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+
+TEST(MathLibInterpolationAlgorithms, PiecewiseLinearInterpolation)
+{
+ const std::size_t size(1000);
+ std::vector<double> supp_pnts, values;
+ for (std::size_t k(0); k<size; ++k) {
+ supp_pnts.push_back(static_cast<double>(k));
+ if (k % 2 == 0) {
+ values.push_back(0.0);
+ } else {
+ values.push_back(1.0);
+ }
+ }
+
+ MathLib::PiecewiseLinearInterpolation interpolation(supp_pnts, values);
+ // Interpolation
+ for (std::size_t k(0); k<size-1; ++k) {
+ ASSERT_NEAR(0.5, interpolation.getValue(k+0.5), std::numeric_limits<double>::epsilon());
+ }
+
+ for (std::size_t k(0); k<size-1; ++k) {
+ if (k % 2 == 0) {
+ ASSERT_NEAR(0.25, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.75, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
+ } else {
+ ASSERT_NEAR(0.75, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.25, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
+ }
+ }
+
+ // Extrapolation
+ ASSERT_NEAR(-0.5, interpolation.getValue(-0.5), std::numeric_limits<double>::epsilon());
+ // Extrapolation
+ ASSERT_NEAR(1.5, interpolation.getValue(size-0.5), std::numeric_limits<double>::epsilon());
+}
+
+TEST(MathLibInterpolationAlgorithms, PiecewiseLinearInterpolationSupportPntsInReverseOrder)
+{
+ const std::size_t size(1000);
+ std::vector<double> supp_pnts, values;
+ for (std::size_t k(0); k<size; ++k) {
+ supp_pnts.push_back(static_cast<double>(size-1-k));
+ if (k % 2 == 0) {
+ values.push_back(1.0);
+ } else {
+ values.push_back(0.0);
+ }
+ }
+
+ MathLib::PiecewiseLinearInterpolation interpolation(supp_pnts, values);
+ // Interpolation
+ for (std::size_t k(0); k<size-1; ++k) {
+ ASSERT_NEAR(0.5, interpolation.getValue(k+0.5), std::numeric_limits<double>::epsilon());
+ }
+
+ for (std::size_t k(0); k<size-1; ++k) {
+ if (k % 2 == 0) {
+ ASSERT_NEAR(0.25, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.75, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
+ } else {
+ ASSERT_NEAR(0.75, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0.25, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
+ }
+ }
+
+ // Extrapolation
+ ASSERT_NEAR(-0.5, interpolation.getValue(-0.5), std::numeric_limits<double>::epsilon());
+ // Extrapolation
+ ASSERT_NEAR(1.5, interpolation.getValue(size-0.5), std::numeric_limits<double>::epsilon());
+}
diff --git a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
index 32d372b31059079a4c254cf9881c088664b21acc..d4eb23402531d0e3e32e555a9297b1f8b667b1bc 100644
--- a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
+++ b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
@@ -17,14 +17,14 @@
// ThirdParty/logog
#include "logog/include/logog.hpp"
// BaseLib
+#include "LogogSimpleFormatter.h"
#include "StringTools.h"
#include "quicksort.h"
-#include "LogogSimpleFormatter.h"
// FileIO/Legacy
+#include "FileTools.h"
#include "MeshIO.h"
#include "readMeshFromFile.h"
-#include "FileTools.h"
// GeoLib
#include "Raster.h"
@@ -33,11 +33,11 @@
#include "VtkMeshConverter.h"
// MeshLib
-#include "Mesh.h"
+#include "ConvertRasterToMesh.h"
#include "Elements/Element.h"
-#include "MshEnums.h"
+#include "Mesh.h"
#include "Mesh2MeshPropertyInterpolation.h"
-#include "ConvertRasterToMesh.h"
+#include "MshEnums.h"
int main (int argc, char* argv[])
{
@@ -46,24 +46,58 @@ int main (int argc, char* argv[])
BaseLib::LogogSimpleFormatter *custom_format (new BaseLib::LogogSimpleFormatter);
logog_cout->SetFormatter(*custom_format);
- TCLAP::CmdLine cmd("Generates properties for mesh elements of an input mesh deploying a ASC raster file", ' ', "0.1");
-
- TCLAP::ValueArg<std::string> out_mesh_arg("o", "out-mesh", "the mesh is stored to a file of this name", false, "", "filename for mesh output");
+ TCLAP::CmdLine cmd(
+ "Generates properties for mesh elements of an input mesh deploying a ASC raster file",
+ ' ',
+ "0.1");
+
+ TCLAP::ValueArg<std::string> out_mesh_arg("o",
+ "out-mesh",
+ "the mesh is stored to a file of this name",
+ false,
+ "",
+ "filename for mesh output");
cmd.add( out_mesh_arg );
- TCLAP::ValueArg<bool> refinement_raster_output_arg("","output-refined-raster", "write refined raster to a new ASC file", false, false, "0");
+ TCLAP::ValueArg<bool> refinement_raster_output_arg("",
+ "output-refined-raster",
+ "write refined raster to a new ASC file",
+ false,
+ false,
+ "0");
cmd.add( refinement_raster_output_arg );
- TCLAP::ValueArg<unsigned> refinement_arg("r","refine","refinement factor that raises the resolution of the raster data", false, 1, "factor (default = 1)");
+ TCLAP::ValueArg<unsigned> refinement_arg(
+ "r",
+ "refine",
+ "refinement factor that raises the resolution of the raster data",
+ false,
+ 1,
+ "factor (default = 1)");
cmd.add( refinement_arg );
- TCLAP::ValueArg<std::string> mapping_arg("","mapping-name","file name of mapping", true, "", "file name");
+ TCLAP::ValueArg<std::string> mapping_arg("",
+ "mapping-name",
+ "file name of mapping",
+ true,
+ "",
+ "file name");
cmd.add( mapping_arg );
- TCLAP::ValueArg<std::string> raster_arg("","raster-file","the name of the ASC raster file", true, "", "file name");
+ TCLAP::ValueArg<std::string> raster_arg("",
+ "raster-file",
+ "the name of the ASC raster file",
+ true,
+ "",
+ "file name");
cmd.add( raster_arg );
- TCLAP::ValueArg<std::string> mesh_arg("m", "mesh", "the mesh is read from this file", true, "test.msh", "file name");
+ TCLAP::ValueArg<std::string> mesh_arg("m",
+ "mesh",
+ "the mesh is read from this file",
+ true,
+ "test.msh",
+ "file name");
cmd.add( mesh_arg );
cmd.parse( argc, argv );
@@ -88,7 +122,7 @@ int main (int argc, char* argv[])
// put raster data in a std::vector
GeoLib::Raster::const_iterator raster_it(raster->begin());
unsigned n_cols(raster->getNCols()), n_rows(raster->getNRows());
- std::vector<double> src_properties(n_cols*n_rows);
+ std::vector<double> src_properties(n_cols * n_rows);
for (unsigned row(0); row<n_rows; row++) {
for (unsigned col(0); col<n_cols; col++) {
src_properties[row*n_cols+col] = *raster_it;
@@ -98,9 +132,8 @@ int main (int argc, char* argv[])
{
double src_mean_value(src_properties[0]);
- for (size_t k(1); k<n_cols*n_rows; k++) {
+ for (size_t k(1); k < n_cols * n_rows; k++)
src_mean_value += src_properties[k];
- }
src_mean_value /= n_cols*n_rows;
std::cout << "mean value of source: " << src_mean_value << std::endl;
@@ -108,16 +141,17 @@ int main (int argc, char* argv[])
for (size_t k(1); k<n_cols*n_rows; k++) {
src_varianz += MathLib::fastpow(src_properties[k] - src_mean_value, 2);
}
- src_varianz /= n_cols*n_rows;
+ src_varianz /= n_cols * n_rows;
std::cout << "variance of source: " << src_varianz << std::endl;
}
MeshLib::Mesh* src_mesh(MeshLib::ConvertRasterToMesh(*raster, MshElemType::QUAD,
MeshLib::UseIntensityAs::MATERIAL).execute());
- std::vector<size_t> src_perm(n_cols*n_rows);
- for (size_t k(0); k<n_cols*n_rows; k++) src_perm[k] = k;
- BaseLib::Quicksort<double>(src_properties, 0, n_cols*n_rows, src_perm);
+ std::vector<size_t> src_perm(n_cols * n_rows);
+ for (size_t k(0); k < n_cols * n_rows; k++)
+ src_perm[k] = k;
+ BaseLib::Quicksort<double, std::size_t>(src_properties, 0, n_cols * n_rows, src_perm);
// compress the property data structure
const size_t mat_map_size(src_properties.size());
@@ -187,7 +221,7 @@ int main (int argc, char* argv[])
if (! out_mesh_arg.getValue().empty()) {
std::vector<size_t> dest_perm(n_dest_mesh_elements);
for (size_t k(0); k<n_dest_mesh_elements; k++) dest_perm[k] = k;
- BaseLib::Quicksort<double>(dest_properties, 0, n_dest_mesh_elements, dest_perm);
+ BaseLib::Quicksort<double, std::size_t>(dest_properties, 0, n_dest_mesh_elements, dest_perm);
// reset materials in destination mesh
for (size_t k(0); k<n_dest_mesh_elements; k++) {