From 1c2f40fadc5c667156d79925967722ff8e915704 Mon Sep 17 00:00:00 2001
From: rinkk <karsten.rink@ufz.de>
Date: Fri, 29 Jan 2016 15:56:23 +0100
Subject: [PATCH] introduced RasterHeader container and removed lots of
parameters
---
.../DataView/NetCdfConfigureDialog.cpp | 12 +--
.../DataExplorer/VtkVis/VtkRaster.cpp | 26 +++----
Applications/DataExplorer/VtkVis/VtkRaster.h | 9 +--
FileIO/AsciiRasterInterface.cpp | 68 ++++++++---------
FileIO/AsciiRasterInterface.h | 11 +--
GeoLib/Raster.cpp | 73 +++++++++---------
GeoLib/Raster.h | 62 ++++++---------
MeshLib/MeshGenerators/VtkMeshConverter.cpp | 75 ++++++++-----------
MeshLib/MeshGenerators/VtkMeshConverter.h | 20 ++---
9 files changed, 151 insertions(+), 205 deletions(-)
diff --git a/Applications/DataExplorer/DataView/NetCdfConfigureDialog.cpp b/Applications/DataExplorer/DataView/NetCdfConfigureDialog.cpp
index 490f89bf42c..7866dfc23bf 100644
--- a/Applications/DataExplorer/DataView/NetCdfConfigureDialog.cpp
+++ b/Applications/DataExplorer/DataView/NetCdfConfigureDialog.cpp
@@ -3,6 +3,8 @@
#include "NetCdfConfigureDialog.h"
+#include "MathLib/Point3d.h"
+#include "GeoLib/Raster.h"
#include "MeshGenerators/VtkMeshConverter.h"
#include "MeshLib/MeshEnums.h"
#include "VtkGeoImageSource.h"
@@ -352,13 +354,14 @@ void NetCdfConfigureDialog::createDataObject()
double origin_x = (originLon < lastLon) ? originLon : lastLon;
double origin_y = (originLat < lastLat) ? originLat : lastLat;
- double originNetCdf[3] = {origin_x, origin_y, 0};
+ MathLib::Point3d origin({{origin_x, origin_y, 0}});
double resolution = (doubleSpinBoxResolution->value());
if (originLat > lastLat) // reverse lines in vertical direction if the original file has its origin in the northwest corner
this->reverseNorthSouth(data_array, sizeLon, sizeLat);
+ GeoLib::RasterHeader header = {sizeLon,sizeLat,origin,resolution,-9999};
if (this->radioMesh->isChecked())
{
MeshLib::MeshElemType meshElemType = MeshLib::MeshElemType::QUAD;
@@ -375,14 +378,13 @@ void NetCdfConfigureDialog::createDataObject()
}else{
useIntensity = MeshLib::UseIntensityAs::DATAVECTOR;
}
-
- _currentMesh = MeshLib::VtkMeshConverter::convertImgToMesh(data_array,originNetCdf,sizeLon,sizeLat,resolution,meshElemType,useIntensity);
+ _currentMesh = MeshLib::VtkMeshConverter::convertImgToMesh(data_array, header, meshElemType, useIntensity);
}
else
{
- vtkImageImport* image = VtkRaster::loadImageFromArray(data_array, originNetCdf[0], originNetCdf[1], sizeLon, sizeLat, resolution, -9999.0);
+ vtkImageImport* image = VtkRaster::loadImageFromArray(data_array, header);
_currentRaster = VtkGeoImageSource::New();
- _currentRaster->setImage(image, QString::fromStdString(this->getName()), originNetCdf[0], originNetCdf[1], resolution);
+ _currentRaster->setImage(image, QString::fromStdString(this->getName()), origin[0], origin[1], resolution);
}
delete[] length;
diff --git a/Applications/DataExplorer/VtkVis/VtkRaster.cpp b/Applications/DataExplorer/VtkVis/VtkRaster.cpp
index dea54ab1f38..c04b1644b61 100644
--- a/Applications/DataExplorer/VtkVis/VtkRaster.cpp
+++ b/Applications/DataExplorer/VtkVis/VtkRaster.cpp
@@ -54,15 +54,8 @@ vtkImageAlgorithm* VtkRaster::loadImage(const std::string &fileName,
{
raster = FileIO::AsciiRasterInterface::getRasterFromSurferFile(fileName);
}
- if (raster) {
- x0 = raster->getOrigin()[0];
- y0 = raster->getOrigin()[1];
- delta = raster->getRasterPixelSize();
- double const*const data (raster->begin());
- return VtkRaster::loadImageFromArray(data, x0, y0,
- raster->getNCols(), raster->getNRows(), delta,
- raster->getNoDataValue());
- }
+ if (raster)
+ return VtkRaster::loadImageFromArray(raster->begin(), raster->getHeader());
else if ((fileInfo.suffix().toLower() == "tif") || (fileInfo.suffix().toLower() == "tiff"))
{
#ifdef GEOTIFF_FOUND
@@ -75,16 +68,15 @@ vtkImageAlgorithm* VtkRaster::loadImage(const std::string &fileName,
else
return loadImageFromFile(fileName);
}
-
-vtkImageImport* VtkRaster::loadImageFromArray(double const*const data_array, double x0, double y0, std::size_t width, std::size_t height, double delta, double noData)
+vtkImageImport* VtkRaster::loadImageFromArray(double const*const data_array, GeoLib::RasterHeader header)
{
- const unsigned length = height*width;
+ const unsigned length = header.n_rows*header.n_cols;
float* data = new float[length*2];
float max_val = *std::max_element(data_array, data_array+length);
for (unsigned j=0; j<length; ++j)
{
data[j*2] = static_cast<float>(data_array[j]);
- if (fabs(data[j*2]-noData) < std::numeric_limits<double>::epsilon())
+ if (fabs(data[j*2]-header.no_data) < std::numeric_limits<double>::epsilon())
{
data[j*2] = max_val;
data[j*2+1] = 0;
@@ -94,10 +86,10 @@ vtkImageImport* VtkRaster::loadImageFromArray(double const*const data_array, dou
}
vtkImageImport* image = vtkImageImport::New();
- image->SetDataSpacing(delta, delta, delta);
- image->SetDataOrigin(x0+(delta/2.0), y0+(delta/2.0), 0); // translate whole mesh by half a pixel in x and y
- image->SetWholeExtent(0, width-1, 0, height-1, 0, 0);
- image->SetDataExtent(0, width-1, 0, height-1, 0, 0);
+ image->SetDataSpacing(header.cell_size, header.cell_size, header.cell_size);
+ image->SetDataOrigin(header.origin[0]+(header.cell_size/2.0), header.origin[1]+(header.cell_size/2.0), 0); // translate whole mesh by half a pixel in x and y
+ image->SetWholeExtent(0, header.n_rows-1, 0, header.n_cols-1, 0, 0);
+ image->SetDataExtent(0, header.n_rows-1, 0, header.n_cols-1, 0, 0);
image->SetDataExtentToWholeExtent();
image->SetDataScalarTypeToFloat();
image->SetNumberOfScalarComponents(2);
diff --git a/Applications/DataExplorer/VtkVis/VtkRaster.h b/Applications/DataExplorer/VtkVis/VtkRaster.h
index f5ba39c9e03..ae32757a130 100644
--- a/Applications/DataExplorer/VtkVis/VtkRaster.h
+++ b/Applications/DataExplorer/VtkVis/VtkRaster.h
@@ -15,6 +15,7 @@
#define VTKRASTER_H
#include <string>
+#include "GeoLib/Raster.h"
class vtkImageAlgorithm;
class vtkImageImport;
@@ -32,13 +33,7 @@ public:
/**
* \brief Returns a VtkImageAlgorithm from an array of pixel values and some image meta data.
*/
- static vtkImageImport* loadImageFromArray(double const*const data_array,
- double x0,
- double y0,
- std::size_t width,
- std::size_t height,
- double delta,
- double no_data = -9999);
+ static vtkImageImport* loadImageFromArray(double const*const data_array, GeoLib::RasterHeader header);
/**
* \brief Loads an image- or raster-file into an vtkImageAlgorithm-Object.
*
diff --git a/FileIO/AsciiRasterInterface.cpp b/FileIO/AsciiRasterInterface.cpp
index 8ef729016c7..bfa3f081163 100644
--- a/FileIO/AsciiRasterInterface.cpp
+++ b/FileIO/AsciiRasterInterface.cpp
@@ -43,24 +43,21 @@ GeoLib::Raster* AsciiRasterInterface::getRasterFromASCFile(std::string const& fn
}
// header information
- std::size_t n_cols(0), n_rows(0);
- double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), no_data_val(-9999);
-
- if (readASCHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, no_data_val)) {
- double* values = new double[n_cols*n_rows];
+ GeoLib::RasterHeader header;
+ if (readASCHeader(in, header)) {
+ double* values = new double[header.n_cols*header.n_rows];
std::string s;
// read the data into the double-array
- for (std::size_t j(0); j < n_rows; ++j) {
- const std::size_t idx ((n_rows - j - 1) * n_cols);
- for (std::size_t i(0); i < n_cols; ++i) {
+ for (std::size_t j(0); j < header.n_rows; ++j) {
+ const std::size_t idx ((header.n_rows - j - 1) * header.n_cols);
+ for (std::size_t i(0); i < header.n_cols; ++i) {
in >> s;
values[idx+i] = strtod(BaseLib::replaceString(",", ".", s).c_str(),0);
}
}
in.close();
- GeoLib::Raster *raster(new GeoLib::Raster(n_cols, n_rows, xllcorner, yllcorner,
- cell_size, values, values+n_cols*n_rows, no_data_val));
+ GeoLib::Raster *raster(new GeoLib::Raster(header, values, values+header.n_cols*header.n_rows));
delete [] values;
return raster;
} else {
@@ -69,45 +66,45 @@ GeoLib::Raster* AsciiRasterInterface::getRasterFromASCFile(std::string const& fn
}
}
-bool AsciiRasterInterface::readASCHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size, double &no_data_val)
+bool AsciiRasterInterface::readASCHeader(std::ifstream &in, GeoLib::RasterHeader &header)
{
std::string tag, value;
in >> tag;
if (tag.compare("ncols") == 0) {
in >> value;
- n_cols = atoi(value.c_str());
+ header.n_cols = atoi(value.c_str());
} else return false;
in >> tag;
if (tag.compare("nrows") == 0) {
in >> value;
- n_rows = atoi(value.c_str());
+ header.n_rows = atoi(value.c_str());
} else return false;
in >> tag;
if (tag.compare("xllcorner") == 0 || tag.compare("xllcenter") == 0) {
in >> value;
- xllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ header.origin[0] = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
in >> tag;
if (tag.compare("yllcorner") == 0 || tag.compare("yllcenter") == 0) {
in >> value;
- yllcorner = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ header.origin[1] = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
+ header.origin[2] = 0;
in >> tag;
if (tag.compare("cellsize") == 0) {
in >> value;
- cell_size = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ header.cell_size = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
in >> tag;
if (tag.compare("NODATA_value") == 0 || tag.compare("nodata_value") == 0) {
in >> value;
- no_data_val = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
+ header.no_data = strtod(BaseLib::replaceString(",", ".", value).c_str(), 0);
} else return false;
return true;
@@ -123,19 +120,19 @@ GeoLib::Raster* AsciiRasterInterface::getRasterFromSurferFile(std::string const&
}
// header information
- std::size_t n_cols(0), n_rows(0);
- double xllcorner(0.0), yllcorner(0.0), cell_size(0.0), min(0.0), max(0.0);
+ GeoLib::RasterHeader header;
+ double min(0.0), max(0.0);
- if (readSurferHeader(in, n_cols, n_rows, xllcorner, yllcorner, cell_size, min, max))
+ if (readSurferHeader(in, header, min, max))
{
const double no_data_val (min-1);
- double* values = new double[n_cols*n_rows];
+ double* values = new double[header.n_cols*header.n_rows];
std::string s;
// read the data into the double-array
- for (std::size_t j(0); j < n_rows; ++j)
+ for (std::size_t j(0); j < header.n_rows; ++j)
{
- const std::size_t idx (j * n_cols);
- for (std::size_t i(0); i < n_cols; ++i)
+ const std::size_t idx (j * header.n_cols);
+ for (std::size_t i(0); i < header.n_cols; ++i)
{
in >> s;
const double val (strtod(BaseLib::replaceString(",", ".", s).c_str(),0));
@@ -143,8 +140,7 @@ GeoLib::Raster* AsciiRasterInterface::getRasterFromSurferFile(std::string const&
}
}
in.close();
- GeoLib::Raster *raster(new GeoLib::Raster(n_cols, n_rows, xllcorner, yllcorner, cell_size,
- values, values+n_cols*n_rows, no_data_val));
+ GeoLib::Raster *raster(new GeoLib::Raster(header, values, values+header.n_cols*header.n_rows));
delete [] values;
return raster;
} else {
@@ -153,8 +149,8 @@ GeoLib::Raster* AsciiRasterInterface::getRasterFromSurferFile(std::string const&
}
}
-bool AsciiRasterInterface::readSurferHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size, double &min, double &max)
+bool AsciiRasterInterface::readSurferHeader(std::ifstream &in, GeoLib::RasterHeader &header,
+ double &min, double &max)
{
std::string tag;
@@ -167,21 +163,23 @@ bool AsciiRasterInterface::readSurferHeader(std::ifstream &in, std::size_t &n_co
}
else
{
- in >> n_cols >> n_rows;
+ in >> header.n_cols >> header.n_rows;
in >> min >> max;
- xllcorner = min;
- cell_size = (max-min)/static_cast<double>(n_cols);
+ header.origin[0] = min;
+ header.cell_size = (max-min)/static_cast<double>(header.n_cols);
in >> min >> max;
- yllcorner = min;
+ header.origin[1] = min;
+ header.origin[2] = 0;
- if (ceil((max-min)/static_cast<double>(n_rows)) == ceil(cell_size))
- cell_size = ceil(cell_size);
+ if (ceil((max-min)/static_cast<double>(header.n_rows)) == ceil(header.cell_size))
+ header.cell_size = ceil(header.cell_size);
else
{
ERR("Error in readSurferHeader() - Anisotropic cellsize detected.");
return 0;
}
+ header.no_data = min-1;
in >> min >> max;
}
diff --git a/FileIO/AsciiRasterInterface.h b/FileIO/AsciiRasterInterface.h
index 6cd7bc1a413..f4898c3c1c3 100644
--- a/FileIO/AsciiRasterInterface.h
+++ b/FileIO/AsciiRasterInterface.h
@@ -16,9 +16,7 @@
#include <fstream>
-namespace GeoLib {
- class Raster;
-}
+#include "GeoLib/Raster.h"
namespace FileIO {
@@ -44,13 +42,10 @@ public:
private:
/// Reads the header of a Esri asc-file.
- static bool readASCHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size,
- double &no_data_val);
+ static bool readASCHeader(std::ifstream &in, GeoLib::RasterHeader &header);
/// Reads the header of a Surfer grd-file.
- static bool readSurferHeader(std::ifstream &in, std::size_t &n_cols, std::size_t &n_rows,
- double &xllcorner, double &yllcorner, double &cell_size,
+ static bool readSurferHeader(std::ifstream &in, GeoLib::RasterHeader &header,
double &min, double &max);
};
diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp
index 85b2c247ada..e5a2f079462 100644
--- a/GeoLib/Raster.cpp
+++ b/GeoLib/Raster.cpp
@@ -28,13 +28,13 @@ namespace GeoLib {
void Raster::refineRaster(std::size_t scaling)
{
- double *new_raster_data(new double[_n_rows*_n_cols*scaling*scaling]);
+ double *new_raster_data(new double[_header.n_rows*_header.n_cols*scaling*scaling]);
- for (std::size_t row(0); row<_n_rows; row++) {
- for (std::size_t col(0); col<_n_cols; col++) {
- const std::size_t idx(row*_n_cols+col);
+ for (std::size_t row(0); row<_header.n_rows; row++) {
+ for (std::size_t col(0); col<_header.n_cols; col++) {
+ const std::size_t idx(row*_header.n_cols+col);
for (std::size_t new_row(row*scaling); new_row<(row+1)*scaling; new_row++) {
- const std::size_t idx0(new_row*_n_cols*scaling);
+ const std::size_t idx0(new_row*_header.n_cols*scaling);
for (std::size_t new_col(col*scaling); new_col<(col+1)*scaling; new_col++) {
new_raster_data[idx0+new_col] = _raster_data[idx];
}
@@ -43,9 +43,9 @@ void Raster::refineRaster(std::size_t scaling)
}
std::swap(_raster_data, new_raster_data);
- _cell_size /= scaling;
- _n_cols *= scaling;
- _n_rows *= scaling;
+ _header.cell_size /= scaling;
+ _header.n_cols *= scaling;
+ _header.n_rows *= scaling;
delete [] new_raster_data;
}
@@ -57,17 +57,17 @@ Raster::~Raster()
void Raster::setCellSize(double cell_size)
{
- _cell_size = cell_size;
+ _header.cell_size = cell_size;
}
void Raster::setNoDataVal (double no_data_val)
{
- _no_data_val = no_data_val;
+ _header.no_data = no_data_val;
}
-GeoLib::Point const& Raster::getOrigin() const
+GeoLib::Point const Raster::getOrigin() const
{
- return _ll_pnt;
+ return GeoLib::Point(_header.origin,0);
}
Raster* Raster::getRasterFromSurface(Surface const& sfc, double cell_size, double no_data_val)
@@ -100,42 +100,43 @@ Raster* Raster::getRasterFromSurface(Surface const& sfc, double cell_size, doubl
}
}
- return new Raster(n_cols, n_rows, ll[0], ll[1], cell_size, z_vals, z_vals+n_cols*n_rows ,-9999);
+ RasterHeader header = {std::size_t(n_cols), std::size_t(n_rows), MathLib::Point3d(ll), cell_size, static_cast<double>(-9999)};
+ return new Raster(header, z_vals, z_vals+n_cols*n_rows);
}
double Raster::getValueAtPoint(const MathLib::Point3d &pnt) const
{
- if (pnt[0]>=_ll_pnt[0] && pnt[0]<(_ll_pnt[0]+(_cell_size*_n_cols)) &&
- pnt[1]>=_ll_pnt[1] && pnt[1]<(_ll_pnt[1]+(_cell_size*_n_rows)))
+ if (pnt[0]>=_header.origin[0] && pnt[0]<(_header.origin[0]+(_header.cell_size*_header.n_cols)) &&
+ pnt[1]>=_header.origin[1] && pnt[1]<(_header.origin[1]+(_header.cell_size*_header.n_rows)))
{
- int cell_x = static_cast<int>(std::floor((pnt[0] - _ll_pnt[0])/_cell_size));
- int cell_y = static_cast<int>(std::floor((pnt[1] - _ll_pnt[1])/_cell_size));
+ int cell_x = static_cast<int>(std::floor((pnt[0] - _header.origin[0])/_header.cell_size));
+ int cell_y = static_cast<int>(std::floor((pnt[1] - _header.origin[1])/_header.cell_size));
// use raster boundary values if node is outside raster due to rounding errors or floating point arithmetic
- cell_x = (cell_x < 0) ? 0 : ((cell_x > static_cast<int>(_n_cols)) ?
- static_cast<int>(_n_cols-1) : cell_x);
- cell_y = (cell_y < 0) ? 0 : ((cell_y > static_cast<int>(_n_rows)) ?
- static_cast<int>(_n_rows-1) : cell_y);
+ cell_x = (cell_x < 0) ? 0 : ((cell_x > static_cast<int>(_header.n_cols)) ?
+ static_cast<int>(_header.n_cols-1) : cell_x);
+ cell_y = (cell_y < 0) ? 0 : ((cell_y > static_cast<int>(_header.n_rows)) ?
+ static_cast<int>(_header.n_rows-1) : cell_y);
- const std::size_t index = cell_y*_n_cols+cell_x;
+ const std::size_t index = cell_y*_header.n_cols+cell_x;
return _raster_data[index];
}
- return _no_data_val;
+ return _header.no_data;
}
double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
{
// position in raster
- double const xPos ((pnt[0] - _ll_pnt[0]) / _cell_size);
- double const yPos ((pnt[1] - _ll_pnt[1]) / _cell_size);
+ double const xPos ((pnt[0] - _header.origin[0]) / _header.cell_size);
+ double const yPos ((pnt[1] - _header.origin[1]) / _header.cell_size);
// raster cell index
double const xIdx (std::floor(xPos)); //carry out computions in double
double const yIdx (std::floor(yPos)); // so not to over- or underflow.
// weights for bilinear interpolation
- double const half_delta = 0.5*_cell_size;
- double const xShift = std::fabs(xPos-(xIdx+half_delta)) / _cell_size;
- double const yShift = std::fabs(yPos-(yIdx+half_delta)) / _cell_size;
+ double const half_delta = 0.5*_header.cell_size;
+ double const xShift = std::fabs(xPos-(xIdx+half_delta)) / _header.cell_size;
+ double const yShift = std::fabs(yPos-(yIdx+half_delta)) / _header.cell_size;
std::array<double,4> weight = {{ (1-xShift)*(1-yShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift }};
// neighbors to include in interpolation
@@ -153,16 +154,16 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
// a no data value. This also allows the cast to unsigned type.
if ( (xIdx + x_nb[j]) < 0 ||
(yIdx + y_nb[j]) < 0 ||
- (xIdx + x_nb[j]) > (_n_cols-1) ||
- (yIdx + y_nb[j]) > (_n_rows-1) )
- pix_val[j] = _no_data_val;
+ (xIdx + x_nb[j]) > (_header.n_cols-1) ||
+ (yIdx + y_nb[j]) > (_header.n_rows-1) )
+ pix_val[j] = _header.no_data;
else
pix_val[j] = _raster_data[
- static_cast<std::size_t>(yIdx + y_nb[j]) * _n_cols +
+ static_cast<std::size_t>(yIdx + y_nb[j]) * _header.n_cols +
static_cast<std::size_t>(xIdx + x_nb[j])];
// remove no data values
- if (std::fabs(pix_val[j] - _no_data_val) < std::numeric_limits<double>::epsilon())
+ if (std::fabs(pix_val[j] - _header.no_data) < std::numeric_limits<double>::epsilon())
{
weight[j] = 0;
no_data_count++;
@@ -173,7 +174,7 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
if (no_data_count > 0)
{
if (no_data_count == 4) // if there is absolutely no data just use the default value
- return _no_data_val;
+ return _header.no_data;
const double norm = 1.0 / (weight[0]+weight[1]+weight[2]+weight[3]);
std::for_each(weight.begin(), weight.end(), [&norm](double &val){val*=norm;});
@@ -185,8 +186,8 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
bool Raster::isPntOnRaster(MathLib::Point3d const& pnt) const
{
- if ((pnt[0]<_ll_pnt[0]) || (pnt[0]>_ll_pnt[0]+(_n_cols*_cell_size)) ||
- (pnt[1]<_ll_pnt[1]) || (pnt[1]>_ll_pnt[1]+(_n_rows*_cell_size)))
+ if ((pnt[0]<_header.origin[0]) || (pnt[0]>_header.origin[0]+(_header.n_cols*_header.cell_size)) ||
+ (pnt[1]<_header.origin[1]) || (pnt[1]>_header.origin[1]+(_header.n_rows*_header.cell_size)))
return false;
return true;
}
diff --git a/GeoLib/Raster.h b/GeoLib/Raster.h
index abb05c32a69..78183ee72f5 100644
--- a/GeoLib/Raster.h
+++ b/GeoLib/Raster.h
@@ -18,6 +18,16 @@
namespace GeoLib {
+/// Contains the relevant information when handling with geoscientific raster data
+struct RasterHeader
+{
+ std::size_t n_cols; // width
+ std::size_t n_rows; // height
+ MathLib::Point3d origin; // lower left corner
+ double cell_size; // edge length of each pixel
+ double no_data; // no data value
+};
+
/**
* @brief Class Raster is used for managing raster data.
*
@@ -36,20 +46,13 @@ public:
* to be handed over via input iterators. Deploying iterators has the
* advantage that the use of the class is independent from the input
* container.
- * @param n_cols number of columns
- * @param n_rows number of rows
- * @param xllcorner the \f$x\f$ coordinate of lower left point
- * @param yllcorner the \f$y\f$ coordinate of lower left point
- * @param cell_size the size of a raster pixel
+ * @param header meta-information about the raster (height, width, etc.)
* @param begin input iterator pointing to the first element of the data
* @param end input iterator pointing to the last element of the data, end have to be reachable from begin
- * @param no_data_val value for raster pixels that have not valid values
*/
template<typename InputIterator>
- Raster(std::size_t n_cols, std::size_t n_rows, double xllcorner, double yllcorner,
- double cell_size, InputIterator begin, InputIterator end, double no_data_val = -9999) :
- _n_cols(n_cols), _n_rows(n_rows), _ll_pnt(xllcorner, yllcorner, 0.0), _cell_size(cell_size),
- _no_data_val(no_data_val), _raster_data(new double[n_cols*n_rows])
+ Raster(RasterHeader header, InputIterator begin, InputIterator end) :
+ _header(header), _raster_data(new double[_header.n_cols*_header.n_rows])
{
iterator raster_it(_raster_data);
for (InputIterator it(begin); it != end; ++it) {
@@ -58,32 +61,35 @@ public:
}
}
+ /// Returns the complete header information
+ RasterHeader const& getHeader() const { return _header; }
+
/**
* get the number of columns for the raster
*/
- std::size_t getNCols() const { return _n_cols; }
+ std::size_t getNCols() const { return _header.n_cols; }
/**
* get the number of rows for the raster
*/
- std::size_t getNRows() const { return _n_rows; }
+ std::size_t getNRows() const { return _header.n_rows; }
/**
* get the distance between raster pixels
*/
- double getRasterPixelSize() const { return _cell_size; }
+ double getRasterPixelSize() const { return _header.cell_size; }
/**
* get the origin of lower left raster cell
* @return the origin of the raster
*/
- GeoLib::Point const& getOrigin() const;
+ GeoLib::Point const getOrigin() const;
/**
* Refine the raster using scaling as a refinement parameter.
*/
void refineRaster(std::size_t scaling);
- double getNoDataValue() const { return _no_data_val; }
+ double getNoDataValue() const { return _header.no_data; }
/**
* Constant iterator that is pointing to the first raster pixel value.
@@ -94,7 +100,7 @@ public:
* Constant iterator that is pointing to the last raster pixel value.
* @return constant iterator
*/
- const_iterator end() const { return _raster_data + _n_rows*_n_cols; }
+ const_iterator end() const { return _raster_data + _header.n_rows*_header.n_cols; }
/**
* Returns the raster value at the position of the given point.
@@ -116,29 +122,7 @@ private:
void setCellSize(double cell_size);
void setNoDataVal (double no_data_val);
- /**
- * number of columns of the raster
- */
- std::size_t _n_cols;
- /**
- * number of rows of the raster
- */
- std::size_t _n_rows;
- /**
- * lower left point of the raster
- */
- GeoLib::Point _ll_pnt;
- /**
- * cell size - each cell is quadratic
- */
- double _cell_size;
- /**
- * value for cells there is no data for
- */
- double _no_data_val;
- /**
- * raw raster data
- */
+ GeoLib::RasterHeader _header;
double* _raster_data;
};
diff --git a/MeshLib/MeshGenerators/VtkMeshConverter.cpp b/MeshLib/MeshGenerators/VtkMeshConverter.cpp
index 18e876aba01..f7a6f55b086 100644
--- a/MeshLib/MeshGenerators/VtkMeshConverter.cpp
+++ b/MeshLib/MeshGenerators/VtkMeshConverter.cpp
@@ -71,18 +71,17 @@ MeshLib::Mesh* VtkMeshConverter::convertImgToMesh(vtkImageData* img,
return nullptr;
}
- const std::size_t imgHeight = dims[0];
- const std::size_t imgWidth = dims[1];
- const std::size_t incHeight = imgHeight+1;
- const std::size_t incWidth = imgWidth+1;
+ MathLib::Point3d orig ({{origin[0], origin[1], origin[2]}});
+ GeoLib::RasterHeader header = {dims[1], dims[0], orig, scalingFactor, -9999};
+ const std::size_t incHeight = header.n_rows+1;
+ const std::size_t incWidth = header.n_cols+1;
std::vector<double> pix_val (incHeight * incWidth, std::numeric_limits<double>::max());
- std::vector<bool> pix_vis (imgHeight * imgWidth, false);
+ std::vector<bool> pix_vis (header.n_rows * header.n_cols, false);
- for (std::size_t i = 0; i < imgWidth; i++)
- {
- for (std::size_t j = 0; j < imgHeight; j++)
+ for (std::size_t i = 0; i < header.n_cols; i++)
+ for (std::size_t j = 0; j < header.n_rows; j++)
{
- std::size_t const img_idx = i*imgHeight + j;
+ std::size_t const img_idx = i*header.n_rows + j;
std::size_t const fld_idx = i*incHeight + j;
// colour of current pixel
@@ -92,7 +91,7 @@ MeshLib::Mesh* VtkMeshConverter::convertImgToMesh(vtkImageData* img,
if (!visible)
continue;
- double const value = (nTuple < 3) ?
+ double const value = (nTuple < 3) ?
colour[0] : // grey (+ alpha)
(0.3 * colour[0] + 0.6 * colour[1] + 0.1 * colour[2]); // rgb(a)
pix_vis[img_idx] = true;
@@ -101,17 +100,13 @@ MeshLib::Mesh* VtkMeshConverter::convertImgToMesh(vtkImageData* img,
pix_val[fld_idx+incHeight] = value;
pix_val[fld_idx+incHeight+1] = value;
}
- }
-
- return constructMesh(pix_val, pix_vis, origin, imgHeight, imgWidth, scalingFactor, elem_type, intensity_type);
+
+ return constructMesh(pix_val, pix_vis, header, elem_type, intensity_type);
}
MeshLib::Mesh* VtkMeshConverter::convertImgToMesh(
const double* img,
- const double origin[3],
- const std::size_t imgHeight,
- const std::size_t imgWidth,
- const double &scalingFactor,
+ GeoLib::RasterHeader const& header,
MeshElemType elem_type,
UseIntensityAs intensity_type)
{
@@ -121,15 +116,15 @@ MeshLib::Mesh* VtkMeshConverter::convertImgToMesh(
return nullptr;
}
- std::size_t const incHeight (imgHeight+1);
- std::size_t const incWidth (imgWidth+1);
+ std::size_t const incHeight (header.n_rows+1);
+ std::size_t const incWidth (header.n_cols+1);
std::vector<double> pix_val (incHeight * incWidth, std::numeric_limits<double>::max());
std::vector<bool> pix_vis (incHeight * incWidth, false);
- for (std::size_t i = 0; i < imgWidth; i++)
- for (std::size_t j = 0; j < imgHeight; j++)
+ for (std::size_t i = 0; i < header.n_cols; i++)
+ for (std::size_t j = 0; j < header.n_rows; j++)
{
- std::size_t const img_idx = i*imgHeight + j;
+ std::size_t const img_idx = i*header.n_rows + j;
std::size_t const fld_idx = i*incHeight + j;
if (img[img_idx] == -9999)
continue;
@@ -141,26 +136,23 @@ MeshLib::Mesh* VtkMeshConverter::convertImgToMesh(
pix_val[fld_idx+incHeight+1] = img[img_idx];
}
- return constructMesh(pix_val, pix_vis, origin, imgHeight, imgWidth, scalingFactor, elem_type, intensity_type);
+ return constructMesh(pix_val, pix_vis, header, elem_type, intensity_type);
}
MeshLib::Mesh* VtkMeshConverter::constructMesh(
std::vector<double> const& pix_val,
std::vector<bool> const& pix_vis,
- double const origin[3],
- std::size_t const imgHeight,
- std::size_t const imgWidth,
- double const scalingFactor,
+ GeoLib::RasterHeader const& header,
MeshLib::MeshElemType elem_type,
MeshLib::UseIntensityAs intensity_type)
{
- std::vector<int> node_idx_map ((imgHeight+1) * (imgWidth+1), -1);
+ std::vector<int> node_idx_map ((header.n_rows+1) * (header.n_cols+1), -1);
bool const use_elevation (intensity_type == MeshLib::UseIntensityAs::ELEVATION);
- std::vector<MeshLib::Node*> nodes (createNodeVector(pix_val, node_idx_map, imgHeight, imgWidth, origin, scalingFactor, use_elevation));
+ std::vector<MeshLib::Node*> nodes (createNodeVector(pix_val, node_idx_map, header, use_elevation));
if (nodes.empty())
return nullptr;
- std::vector<MeshLib::Element*> elements (createElementVector(pix_val, pix_vis, nodes, node_idx_map, imgHeight, imgWidth, elem_type));
+ std::vector<MeshLib::Element*> elements (createElementVector(pix_val, pix_vis, nodes, node_idx_map, header.n_rows, header.n_cols, elem_type));
if (elements.empty())
return nullptr;
@@ -169,13 +161,13 @@ MeshLib::Mesh* VtkMeshConverter::constructMesh(
{
boost::optional< MeshLib::PropertyVector<int>& > prop_vec =
properties.createNewPropertyVector<int>("MaterialIDs", MeshLib::MeshItemType::Cell, 1);
- fillPropertyVector<int>(*prop_vec, pix_val, pix_vis, imgHeight, imgWidth, elem_type);
- }
+ fillPropertyVector<int>(*prop_vec, pix_val, pix_vis, header.n_rows, header.n_cols, elem_type);
+ }
else if (intensity_type == MeshLib::UseIntensityAs::DATAVECTOR)
{
boost::optional< MeshLib::PropertyVector<double>& > prop_vec =
properties.createNewPropertyVector<double>("Colour", MeshLib::MeshItemType::Cell, 1);
- fillPropertyVector<double>(*prop_vec, pix_val, pix_vis, imgHeight, imgWidth, elem_type);
+ fillPropertyVector<double>(*prop_vec, pix_val, pix_vis, header.n_rows, header.n_cols, elem_type);
}
return new MeshLib::Mesh("RasterDataMesh", nodes, elements, properties);
@@ -184,25 +176,22 @@ MeshLib::Mesh* VtkMeshConverter::constructMesh(
std::vector<MeshLib::Node*> VtkMeshConverter::createNodeVector(
std::vector<double> const& elevation,
std::vector<int> & node_idx_map,
- std::size_t const imgHeight,
- std::size_t const imgWidth,
- double const origin[3],
- double const scalingFactor,
+ GeoLib::RasterHeader const& header,
bool use_elevation)
{
std::size_t node_idx_count(0);
- double const x_offset(origin[0] - scalingFactor/2.0);
- double const y_offset(origin[1] - scalingFactor/2.0);
+ double const x_offset(header.origin[0] - header.cell_size/2.0);
+ double const y_offset(header.origin[1] - header.cell_size/2.0);
std::vector<MeshLib::Node*> nodes;
- for (std::size_t i = 0; i < (imgWidth+1); i++)
- for (std::size_t j = 0; j < (imgHeight+1); j++)
+ for (std::size_t i = 0; i < (header.n_cols+1); i++)
+ for (std::size_t j = 0; j < (header.n_rows+1); j++)
{
- std::size_t const index = i * (imgHeight+1) + j;
+ std::size_t const index = i * (header.n_rows+1) + j;
if (elevation[index] == std::numeric_limits<double>::max())
continue;
double const zValue = (use_elevation) ? elevation[index] : 0;
- MeshLib::Node* node (new MeshLib::Node(x_offset + (scalingFactor * j), y_offset + (scalingFactor * i), zValue));
+ MeshLib::Node* node (new MeshLib::Node(x_offset + (header.cell_size * j), y_offset + (header.cell_size * i), zValue));
nodes.push_back(node);
node_idx_map[index] = node_idx_count;
node_idx_count++;
diff --git a/MeshLib/MeshGenerators/VtkMeshConverter.h b/MeshLib/MeshGenerators/VtkMeshConverter.h
index c8db6ed8b08..b195daade6e 100644
--- a/MeshLib/MeshGenerators/VtkMeshConverter.h
+++ b/MeshLib/MeshGenerators/VtkMeshConverter.h
@@ -22,6 +22,7 @@
#include "logog/include/logog.hpp"
+#include "GeoLib/Raster.h"
#include "MeshLib/Node.h"
#include "MeshLib/Location.h"
#include "MeshLib/MeshEnums.h"
@@ -61,10 +62,7 @@ public:
* \param intensity_type defines how image intensities are interpreted
*/
static MeshLib::Mesh* convertImgToMesh(const double* img,
- const double origin[3],
- const std::size_t imgHeight,
- const std::size_t imgWidth,
- const double &scalingFactor,
+ GeoLib::RasterHeader const& header,
MeshElemType elem_type,
UseIntensityAs intensity_type);
@@ -73,27 +71,19 @@ public:
std::string const& mesh_name = "vtkUnstructuredGrid");
private:
- /// Constructs mesh components based on image data.
static MeshLib::Mesh* constructMesh(
std::vector<double> const& pix_val,
std::vector<bool> const& pix_vis,
- double const origin[3],
- std::size_t const imgHeight,
- std::size_t const imgWidth,
- double const scalingFactor,
+ GeoLib::RasterHeader const& header,
MeshLib::MeshElemType elem_type,
MeshLib::UseIntensityAs intensity_type);
static void convertScalarArrays(vtkUnstructuredGrid &grid, MeshLib::Mesh &mesh);
- /// Creates a mesh node vector based on image data
- static std::vector<MeshLib::Node*> createNodeVector(
+ static std::vector<MeshLib::Node*> VtkMeshConverter::createNodeVector(
std::vector<double> const& elevation,
std::vector<int> &node_idx_map,
- std::size_t const incHeight,
- std::size_t const incWidth,
- double const origin[3],
- double const scalingFactor,
+ GeoLib::RasterHeader const& header,
bool use_elevation);
/// Creates a mesh element vector based on image data
--
GitLab