Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[DE] Added support for world files to georeference raster files #2460

Merged
merged 3 commits into from Apr 10, 2019
Merged
Changes from all commits
Commits
File filter...
Filter file types
Jump to…
Jump to file or symbol
Failed to load files and symbols.

Always

Just for now

@@ -20,31 +20,31 @@

#include <QFileInfo>

#include <vtkBMPReader.h>
#include <vtkImageData.h>
#include <vtkImageImport.h>
#include <vtkImageReader2.h>
#include <vtkPNGReader.h>
#include <vtkJPEGReader.h>
#include <vtkBMPReader.h>
#include <vtkPNGReader.h>
#include <vtkTIFFReader.h>

#ifdef GEOTIFF_FOUND
#include "geo_tiffp.h"
#include "xtiffio.h"
#endif

#include <memory>
#include <logog/include/logog.hpp>
#include <memory>

#include "Applications/FileIO/AsciiRasterInterface.h"
#include "BaseLib/StringTools.h"
#include "GeoLib/Raster.h"


vtkImageAlgorithm* VtkRaster::loadImage(const std::string &fileName,
double& x0, double& y0, double& delta)
{
QFileInfo fileInfo(QString::fromStdString(fileName));


std::unique_ptr<GeoLib::Raster> raster(nullptr);
if (fileInfo.suffix().toLower() == "asc")
raster.reset(FileIO::AsciiRasterInterface::getRasterFromASCFile(fileName));
@@ -61,8 +61,9 @@ vtkImageAlgorithm* VtkRaster::loadImage(const std::string &fileName,
(void)x0;
(void)y0;
(void)delta;
ERR("VtkRaster::loadImage(): Tiff file format not supported in this version!");
return nullptr;
ERR("VtkRaster::loadImage(): GeoTiff file format not supported in this "
"version! Trying to parse as Tiff-file.");
return loadImageFromFile(fileName);
#endif
}

@@ -102,16 +103,23 @@ vtkImageImport* VtkRaster::loadImageFromArray(double const*const data_array, Geo
}

#ifdef GEOTIFF_FOUND
vtkImageImport* VtkRaster::loadImageFromTIFF(const std::string &fileName,
double &x0, double &y0,
double &cellsize)
vtkImageAlgorithm* VtkRaster::loadImageFromTIFF(const std::string& fileName,
double& x0, double& y0,
double& cellsize)
{
TIFF* tiff = XTIFFOpen(fileName.c_str(), "r");

if (tiff)
{
GTIF* geoTiff = GTIFNew(tiff);

int version[3];
int count (0);
GTIFDirectoryInfo(geoTiff, version, &count);

if (count == 0)
WARN("VtkRaster::loadImageFromTIFF - file is not georeferenced.");

if (geoTiff)
{
int imgWidth = 0, imgHeight = 0, nImages = 0, pntCount = 0;
@@ -244,22 +252,101 @@ vtkImageReader2* VtkRaster::loadImageFromFile(const std::string &fileName)

if (fi.suffix().toLower() == "png")
image = vtkPNGReader::New();
else if ((fi.suffix().toLower() == "jpg") || (fi.suffix().toLower() == "jpeg"))
else if ((fi.suffix().toLower() == "tif") ||
(fi.suffix().toLower() == "tiff"))
image = vtkTIFFReader::New();
else if ((fi.suffix().toLower() == "jpg") ||
(fi.suffix().toLower() == "jpeg"))
image = vtkJPEGReader::New();
else if (fi.suffix().toLower() == "bmp")
image = vtkBMPReader::New();
else
{
ERR("VtkRaster::readImageFromFile(): File format not supported, please convert to BMP, JPG or PNG.");
ERR("VtkRaster::readImageFromFile(): File format not supported, please "
"convert to BMP, JPG, PNG or TIFF.");
return nullptr;
}

image->SetFileName(fileName.c_str());
image->GetOutput()->AllocateScalars(VTK_FLOAT, 1);
image->Update();
readWorldFile(fileName, image);
return image;
}

std::string VtkRaster::findWorldFile(std::string const& filename)
{
std::string const no_ext = BaseLib::dropFileExtension(filename);

std::vector<std::string> const supported_extensions =
{ ".pgw", ".pngw", ".pgwx",
".jgw", ".jpgw", ".jgwx",
".tfw", ".tifw", ".tfwx",
".bpw", ".bmpw", ".bpwx",
".wld" };

for (auto& ext : supported_extensions)
{
if (BaseLib::IsFileExisting(no_ext + ext))
return no_ext + ext;
}

// no world file found
return {};
}

bool VtkRaster::readWorldFile(std::string const& filename,
vtkImageReader2* image)
{
std::string const world_file = findWorldFile(filename);
if (world_file.empty())
{
WARN("No world file found. Image is not georeferenced.");
return false;
}

std::ifstream in(world_file.c_str());
if (!in.is_open())
{
ERR("VtkRaster::readWorldFile(): Could not open file %s.",
filename.c_str());
return false;
}

std::string line("");
// x-scaling
if (!std::getline(in, line))
return false;
double const delta_x = BaseLib::str2number<double>(line);
// 2x rotation (ignored)
if (!(std::getline(in, line) && std::getline(in, line)))
return false;
// negative y-scaling
if (!std::getline(in, line))
return false;
double const delta_y = BaseLib::str2number<double>(line);
if (delta_x != -delta_y)
WARN("Anisotropic pixel size detected (%f vs %f). An isotropic "
"spacing of %f is assumed, be aware results may be wrong.",
delta_x, delta_y, delta_x)
// x-translation
if (!std::getline(in, line))
return false;
double const x0 = BaseLib::str2number<double>(line);
// y-translation
if (!std::getline(in, line))
return false;
double const y0 = BaseLib::str2number<double>(line);

int extent[3];
image->GetOutput()->GetDimensions(extent);
image->SetDataSpacing(delta_x, delta_x, delta_x);
// for GIS the origin is on the lower left, for VTK it is on the upper left
double const vtk_y0 = y0 + (extent[1] * delta_y);
image->SetDataOrigin(x0, vtk_y0, 0);
return true;
}

void VtkRaster::uint32toRGBA(const unsigned int s, int* p)
{
p[3] = s / (256 * 256 * 256);
@@ -58,8 +58,9 @@ class VtkRaster
* \param delta The size of each pixel in the image which is needed for correctly displaying the data.
* \return A vtkImageImport-object (derived from vtkImageAlgorithm).
*/
static vtkImageImport* loadImageFromTIFF(const std::string &fileName,
double& x0, double& y0, double& delta);
static vtkImageAlgorithm* loadImageFromTIFF(const std::string& fileName,
double& x0, double& y0,
double& delta);
#endif

/**
@@ -69,6 +70,20 @@ class VtkRaster
*/
static vtkImageReader2* loadImageFromFile(const std::string &fileName);

/**
* Tries to find a world file for the image given by the filename.
* World files can have a number of extensions depending on the programme
* used to write the image and this method just cycles through the
* possibilities, returning the first match it finds.
*/
static std::string findWorldFile(const std::string& filename);

/**
* Tries to find and load the world file associated with a
* BMP/JPG/PNG-file and create a RasterHeader for the image.
*/
static bool readWorldFile(std::string const& filename, vtkImageReader2* image);

/// Converts an uint32-number into a quadruple representing RGBA-colours for a pixel.
static void uint32toRGBA(const unsigned int s, int* p);
};
@@ -0,0 +1,6 @@
60
0
0
-60
5000
200
Git LFS file not shown
@@ -0,0 +1,89 @@
/**
* \copyright
* Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*
*/

#include <memory>

#include "gtest/gtest.h"

#include <vtkImageAlgorithm.h>
#include <vtkImageData.h>
#include <vtkPNGReader.h>
#include <vtkSmartPointer.h>

#include "Applications/DataExplorer/VtkVis/VtkRaster.h"
#include "Applications/FileIO/AsciiRasterInterface.h"
#include "BaseLib/BuildInfo.h"
#include "GeoLib/Raster.h"

TEST(TestVtkRaster, TestPNGReader)
{
std::string name = BaseLib::BuildInfo::data_path + "/FileIO/testraster.png";
double x0, y0, delta;
vtkSmartPointer<vtkImageAlgorithm> img =
VtkRaster::loadImage(name, x0, y0, delta);
img->Update();
EXPECT_TRUE(img != nullptr);
EXPECT_TRUE(img->GetOutput() != nullptr);
double val[3];
img->GetOutput()->GetSpacing(val);
for (std::size_t i = 0; i < 3; ++i)
{
EXPECT_NEAR(60, val[i], std::numeric_limits<double>::epsilon());
}
img->GetOutput()->GetOrigin(val);
EXPECT_NEAR(5000, val[0], std::numeric_limits<double>::epsilon());
EXPECT_NEAR(-400, val[1], std::numeric_limits<double>::epsilon());
EXPECT_NEAR(0, val[2], std::numeric_limits<double>::epsilon());
int extent[3];
img->GetOutput()->GetDimensions(extent);
EXPECT_EQ(12, extent[0]);
EXPECT_EQ(10, extent[1]);
EXPECT_EQ(1, extent[2]);
}

TEST(TestVtkRaster, TestASCReader)
{
std::string name =
BaseLib::BuildInfo::data_path + "/MeshGeoToolsLib/Hamburg/00-raster.asc";
double x0, y0, delta;
vtkSmartPointer<vtkImageAlgorithm> img =
VtkRaster::loadImage(name, x0, y0, delta);
img->Update();
EXPECT_TRUE(img != nullptr);
EXPECT_TRUE(img->GetOutput() != nullptr);
std::unique_ptr<GeoLib::Raster> raster(FileIO::AsciiRasterInterface::getRasterFromASCFile(name));
EXPECT_EQ(298, raster->getHeader().n_cols);
EXPECT_EQ(205, raster->getHeader().n_rows);
EXPECT_NEAR(25, raster->getHeader().cell_size,
std::numeric_limits<double>::epsilon());

double val[3];
img->GetOutput()->GetSpacing(val);
for (std::size_t i = 0; i < 3; ++i)
{
EXPECT_NEAR(raster->getHeader().cell_size, val[i],
std::numeric_limits<double>::epsilon());
}

double const half_pix = raster->getHeader().cell_size / 2.0;
img->GetOutput()->GetOrigin(val);
EXPECT_NEAR(raster->getHeader().origin[0], val[0] - half_pix,
std::numeric_limits<double>::epsilon());
EXPECT_NEAR(raster->getHeader().origin[1], val[1] - half_pix,
std::numeric_limits<double>::epsilon());
EXPECT_NEAR(raster->getHeader().origin[2], val[2],
std::numeric_limits<double>::epsilon());

int extent[3];
img->GetOutput()->GetDimensions(extent);
EXPECT_EQ(raster->getHeader().n_cols, extent[0]);
EXPECT_EQ(raster->getHeader().n_rows, extent[1]);
EXPECT_EQ(raster->getHeader().n_depth, extent[2]);
}
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.