github/pubgeo
1
1
Fork 0
mirror of https://github.com/pubgeo/pubgeo.git synced 2026-01-08 21:48:04 -05:00
Code Issues Packages Projects Releases Wiki Activity

Updated build for PDAL 1.8, added intensity image output to shr3d, and refactored shr3d to better support use in PDAL pipelines.

Browse Source
This commit is contained in:
Scott Almes
2019-01-18 15:06:17 -05:00
parent 0005c375d3
commit 394fb6fdf4
15 changed files with 499 additions and 223 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@@ -1,4 +1,5 @@
build
.gitlab-ci.yml
# Prerequisites
*.d

27
CMakeLists.txt
View File

@@ -8,14 +8,21 @@ option(PUBGEO_INSTALL_ALIGN3D "Install align3d executable." OFF)
option(PUBGEO_INSTALL_SHR3D "Install shr3d executable." OFF)
IF(WIN32)
# Default install directory
MESSAGE(WARNING "Advance INSTALL OSGeo4W64. Include: (Command Line Tools) pdal and its dependencies. Potentially requires change to PDALConfig.cmake.")
SET(CMAKE_PREFIX_PATH "C:\\OSGeo4W64\\")
SET(CMAKE_BUILD_TYPE RelWithDebInfo)
# Default install directory
MESSAGE(WARNING "Advance INSTALL OSGeo4W64. Include: (Command Line Tools) pdal and its dependencies. Potentially requires change to PDALConfig.cmake.")
SET(CMAKE_PREFIX_PATH "C:\\OSGeo4W64\\")
# This part is for auto-defines in windows libraries that cause macro errors in our code
add_definitions(-DWIN32_LEAN_AND_MEAN -DNOMINMAX)
SET(BOOST_INCLUDEDIR "C:\\OSGeo4W64\\include\\boost-1_55")
IF(NOT EXISTS ${BOOST_INCLUDEDIR})
message(SEND_ERROR "Boost directory ${BOOST_INCLUDEDIR} does not exist!")
ENDIF()
SET(CMAKE_BUILD_TYPE RelWithDebInfo)
# This part is for auto-defines in windows libraries that cause macro errors in our code
add_definitions(-DWIN32_LEAN_AND_MEAN -DNOMINMAX)
ELSE()
SET(CMAKE_INSTALL_PREFIX "/usr")
ENDIF()
FIND_PACKAGE(GDAL 1.9.0 REQUIRED)# Using Open Source Geo for Windows installer
@@ -28,9 +35,13 @@ else()
endif()
FIND_PACKAGE(PDAL 1.4.0 REQUIRED)
LINK_DIRECTORIES(${PDAL_LIBRARY_DIRS})
IF(WIN32)
# PDAL also has a requirement for winsock dlls
set(PDAL_LIBRARIES ${PDAL_LIBRARIES} ws2_32)
ENDIF()
if (CMAKE_COMPILER_IS_GNUCXX)
SET( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Werror -O3" )
SET( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Werror -O3" )
ENDIF()
# This makes linking/including a little cleaner looking later

22
Dockerfile
View File

@@ -1,21 +1,18 @@
FROM ubuntu:17.10
FROM pdal/pdal:1.8
MAINTAINER JHUAPL <pubgeo@jhuapl.edu>
RUN apt update && apt upgrade -y && apt install -y --fix-missing --no-install-recommends\
build-essential \
ca-certificates \
cmake \
curl \
gdal-bin \
git \
libgdal-dev \
libpdal-dev \
pdal \
cmake \
libgdal-dev \
&& rm -rf /var/lib/apt/lists/*
RUN cd / && git clone https://github.com/pubgeo/pubgeo
# Make a directory to work out of, and change to it
RUN cd / && mkdir /pubgeo
COPY src /pubgeo/src
COPY CMakeLists.txt /pubgeo/
COPY CMakeSettings.json /pubgeo/
WORKDIR /pubgeo/build/
RUN cmake .. && make -j 10 && make install && mv shr3d align3d /usr/local/bin
WORKDIR /
@@ -24,17 +21,14 @@ WORKDIR /
RUN rm -rf /pubgeo
RUN apt purge -y \
build-essential \
libgdal-dev \
libpdal-dev \
cmake \
git \
libgdal-dev \
&& apt autoremove -y
CMD echo "Please run a valid executable:" && \
echo "docker run --rm -v <path to 3D data>:<mount point (MP)> jhuapl/pubgeo shr3d <MP>/<3D file> DH=2 DZ=1 AGL=2 AREA=50" && \
echo "docker run --rm -v <path to point clouds (PC)>:<mount point (MP)> jhuapl/pubgeo align3d <MP>/<reference pc> <MP>/<target pc> gsd=1.0 maxt=10.0"
# Reminders:
# Data coming out will be in <path to 3D data> or <path to point clouds (pc)>
# --rm will auto delete the container when complete (no need to take up disk space)

12
src/align3d/plugin.cpp
View File

@@ -29,16 +29,16 @@
#include "orthoimage.h"
#include "plugin.hpp"
#include <pdal/pdal_macros.hpp>
namespace pdal
{
static PluginInfo const s_info = PluginInfo(
"filters.align3d", "Shareable High Resolution 3D (Almes et al., 2017)",
"http://pdal.io/stages/filters.align3d.html");
static StaticPluginInfo const s_info = {
"filters.align3d",
"Shareable High Resolution 3D (Almes et al., 2017)",
"http://pdal.io/stages/filters.align3d.html"
};
CREATE_SHARED_PLUGIN(1, 0, Align3dFilter, Filter, s_info)
CREATE_SHARED_STAGE(Align3dFilter, s_info)
std::string Align3dFilter::getName() const
{

5
src/common/Image.h
View File

@@ -108,6 +108,11 @@ namespace pubgeo {
}
}
// Test if object does not contain data
bool empty() {
return !data;
}
/**
* Apply a filter to this image, copying off a reference version for use by the filter algorithm
* and overwriting this image with the new results

3
src/common/PointCloud.cpp
View File

@@ -54,8 +54,7 @@ namespace pubgeo {
pdal::BOX3D box;
pv->calculateBounds(box);
pdal::SpatialReference sr = pv->spatialReference();
zone = sr.computeUTMZone(box);
zone = pv->spatialReference().getUTMZone();
// used later to return points
xOff = (int) floor(box.minx);

12
src/common/PointCloud.h
View File

@@ -45,19 +45,25 @@ namespace pubgeo {
inline float x(int i) const {
if (pv != nullptr)
return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::X, i) - xOff);
throw "Point set has not be initialized.";
throw "Point set has not been initialized.";
}
inline float y(int i) const {
if (pv != nullptr)
return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::Y, i) - yOff);
throw "Point set has not be initialized.";
throw "Point set has not been initialized.";
}
inline float z(int i) const {
if (pv != nullptr)
return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::Z, i) - zOff);
throw "Point set has not be initialized.";
throw "Point set has not been initialized.";
}
inline float i(int i) const {
if (pv != nullptr)
return pv->getFieldAs<float>(pdal::Dimension::Id::Intensity, i);
throw "Point set has not been initialized.";
}
private:

28
src/common/orthoimage.h
View File

@@ -40,6 +40,10 @@ namespace pubgeo {
template<> struct make_long<unsigned int> { typedef unsigned long type; };
template<> struct make_long<unsigned long> { typedef unsigned long long type; };
// Define type trait to make signed version of specific types when applicable
template<typename T> struct make_signed { typedef typename std::make_signed<T>::type type; };
template<> struct make_signed<double> { typedef double type; };
template<> struct make_signed<float> { typedef float type; };
//
// Ortho image template class
//
@@ -47,7 +51,7 @@ namespace pubgeo {
class OrthoImage : public Image<TYPE> {
public:
typedef typename std::make_signed<TYPE>::type STYPE; // Define a signed version of TYPE
typedef typename make_signed<TYPE>::type STYPE; // Define a signed version of TYPE
typedef typename make_long<TYPE>::type LTYPE; // Define a long version of TYPE
double easting;
@@ -71,7 +75,7 @@ namespace pubgeo {
~OrthoImage() {}
// Read any GDAL-supported image.
bool read(char *fileName) {
bool read(const char *fileName) {
// Open the image.
GDALAllRegister();
CPLSetConfigOption("GDAL_DATA", ".\\gdaldata");
@@ -344,14 +348,16 @@ namespace pubgeo {
this->offset = minVal;
this->scale = (maxVal - minVal) / maxImageVal;
// Calculate image width and height.
this->width = (unsigned int) ((pset.bounds.xMax - pset.bounds.xMin) / gsdMeters + 1);
this->height = (unsigned int) ((pset.bounds.yMax - pset.bounds.yMin) / gsdMeters + 1);
if (this->empty()) {
// Calculate image width and height.
this->width = (unsigned int) ((pset.bounds.xMax - pset.bounds.xMin) / gsdMeters + 1);
this->height = (unsigned int) ((pset.bounds.yMax - pset.bounds.yMin) / gsdMeters + 1);
// Allocate an ortho image.
this->Allocate(this->width, this->height);
this->easting = pset.bounds.xMin;
this->northing = pset.bounds.yMin;
// Allocate an ortho image.
this->Allocate(this->width, this->height);
this->easting = pset.bounds.xMin;
this->northing = pset.bounds.yMin;
}
this->zone = pset.zone;
this->gsd = gsdMeters;
@@ -514,7 +520,7 @@ namespace pubgeo {
// Conceptionally, this is the same as a median filter, but instead of comparing the
// cell value to the median, we're comparing it against the specified quantile
void quantileFilter(int rad, TYPE dzScaled, float quantile) {
void quantileFilter(int rad, TYPE dzScaled, float quantile, TYPE voidVal = 0, bool skipVoids = true) {
// Filter quantile
quantile = std::max(0.0f,std::min(1.0f,quantile));
@@ -527,7 +533,7 @@ namespace pubgeo {
// Only replace if it differs by more than dz from the median
if (abs(qValue - static_cast<STYPE>(ref)) > dzScaled)
*val = qValue;
}, rad);
}, rad, voidVal, skipVoids);
}
// Apply a minimum filter to an image (erosion)

4
src/shr3d/CMakeLists.txt
View File

@@ -4,6 +4,10 @@ SET(SHR3D_HEADER_FILES
SET(SHR3D_SOURCE_FILES
shr3d.cpp)
IF(WIN32)
include_directories(${BOOST_INCLUDEDIR})
ENDIF()
ADD_LIBRARY(SHR3D_LIB STATIC ${SHR3D_HEADER_FILES} ${SHR3D_SOURCE_FILES})
TARGET_LINK_LIBRARIES(SHR3D_LIB
PUBLIC

53
src/shr3d/main.cpp
View File

@@ -3,6 +3,7 @@
#include <cstdio>
#include <chrono>
#include <regex>
#include "orthoimage.h"
#include "shr3d.h"
@@ -14,8 +15,9 @@ void printArguments() {
printf(" DZ= vertical uncertainty (meters)\n");
printf(" AGL= minimum building height above ground level (meters)\n");
printf("Options:\n");
printf(" AREA= minimum building area (meters)\n");
printf(" AREA= minimum building area (meters)\n");
printf(" EGM96 set this flag to write vertical datum = EGM96\n");
printf(" BOUNDS=MINX,MAXX,MINY,MAXY set to define image bounds\n");
printf("Examples:\n");
printf(" For EO DSM: shr3d dsm.tif DH=5.0 DZ=1.0 AGL=2 AREA=50.0 EGM96\n");
printf(" For lidar DSM: shr3d dsm.tif DH=1.0 DZ=1.0 AGL=2.0 AREA=50.0\n");
@@ -46,6 +48,17 @@ int main(int argc, char **argv) {
if (strstr(argv[i], "AGL=")) { shr3dder.agl_meters = atof(&(argv[i][4])); }
if (strstr(argv[i], "AREA=")) { shr3dder.min_area_meters = atof(&(argv[i][5])); }
if (strstr(argv[i], "EGM96")) { shr3dder.egm96 = true; }
if (strstr(argv[i], "BOUNDS=")) {
std::string values(&(argv[i][7]));
std::regex fmt("^(.+),(.+),(.+),(.+)$");
if (!std::regex_match(values,fmt)) {
printf("Error: Bounds are in an incorrect format.\n");
printArguments();
return -1;
}
std::istringstream str_bnds(std::regex_replace(values,fmt,"([$1,$2],[$3,$4])"));
str_bnds >> shr3dder.bounds;
}
if (strstr(argv[i], "CONVERT")) { convert = true; }
}
if ((shr3dder.dh_meters == 0.0) || (shr3dder.dz_meters == 0.0) || (shr3dder.agl_meters == 0.0)) {
@@ -75,31 +88,15 @@ int main(int argc, char **argv) {
std::string basename(inputFileName);
// Read DSM as SHORT.
// Input can be GeoTIFF or LAS.
shr3d::OrthoImage<unsigned short> dsmImage, minImage;
printf("Reading DSM as SHORT.\n");
// Input can be GeoTIFF or LAS/BPF.
int len = (int) strlen(inputFileName);
char *ext = &inputFileName[len - 3];
printf("File Type = .%s\n", ext);
if (strcmp(ext, "tif") == 0) {
bool ok = dsmImage.read(readFileName);
if (!ok) return -1;
// Copy DSM
minImage = dsmImage;
// Min filter, replacing only points differing by more than the AGL threshold.
minImage.minFilter(4, (unsigned int) (shr3dder.agl_meters / minImage.scale));
// Fill small voids in the DSM.
minImage.fillVoidsPyramid(true, 2);
if (!shr3dder.setDSM(readFileName))
return -1;
} else if ((strcmp(ext, "las") == 0) || (strcmp(ext, "bpf") == 0)) {
// Read a PSET file (e.g., BPF or LAS) & create DSM and MIN images
shr3d::PointCloud pset;
if (!pset.Read(inputFileName)
|| !shr3dder.createDSM(pset, dsmImage)
|| !shr3dder.createMIN(pset, minImage))
if (!shr3dder.setPSET(inputFileName))
return -1;
outputFilenames[shr3d::DSM] = basename + "_DSM.tif";
@@ -109,14 +106,14 @@ int main(int argc, char **argv) {
}
// Convert horizontal and vertical uncertainty values to bin units.
int dh_bins = MAX(1, (int) floor(shr3dder.dh_meters / dsmImage.gsd));
int dh_bins = MAX(1, (int) floor(shr3dder.dh_meters / shr3dder.getDSM().gsd));
printf("DZ_METERS = %f\n", shr3dder.dz_meters);
printf("DH_METERS = %f\n", shr3dder.dh_meters);
printf("DH_BINS = %d\n", dh_bins);
unsigned int dz_short = (unsigned int) (shr3dder.dz_meters / dsmImage.scale);
unsigned int dz_short = (unsigned int) (shr3dder.dz_meters / shr3dder.getDSM().scale);
printf("DZ_SHORT = %d\n", dz_short);
printf("AGL_METERS = %f\n", shr3dder.agl_meters);
unsigned int agl_short = (unsigned int) (shr3dder.agl_meters / dsmImage.scale);
unsigned int agl_short = (unsigned int) (shr3dder.agl_meters / shr3dder.getDSM().scale);
printf("AGL_SHORT = %d\n", agl_short);
printf("AREA_METERS = %f\n", shr3dder.min_area_meters);
@@ -124,17 +121,23 @@ int main(int argc, char **argv) {
#ifdef DEBUG
outputFilenames[shr3d::MIN] = basename + "_MIN.tif";
outputFilenames[shr3d::DSM2] = basename + "_DSM2.tif";
outputFilenames[shr3d::MINAGL] = basename + "_MINAGL.tif";
outputFilenames[shr3d::LABEL] = basename + "_label.tif";
outputFilenames[shr3d::LABELED_BUILDINGS] = basename + "_building_labels.tif";
outputFilenames[shr3d::LABELED_BUILDINGS_3] = basename + "_building_labels_3.tif";
#endif
outputFilenames[shr3d::DTM] = basename + "_DTM.tif";
outputFilenames[shr3d::INTENSITY] = basename + "_INT.tif";
outputFilenames[shr3d::CLASS] = basename + "_class.tif";
outputFilenames[shr3d::BUILDING] = basename + "_buildings.tif";
outputFilenames[shr3d::BUILDING_OUTLINES] = basename + "_buildings.shp";
// Process
shr3dder.process(dsmImage, minImage, outputFilenames);
try {
shr3dder.process(outputFilenames);
} catch (std::exception& e) {
std::cerr << "Error: " << e.what() << std::endl;
}
// Report total elapsed time.
std::chrono::steady_clock::time_point t1 = std::chrono::steady_clock::now();

98
src/shr3d/plugin.cpp
View File

@@ -20,25 +20,42 @@
// SOFTWARE.
// PDAL Plugin of S. Almes, S. Hagstrom, D. Chilcott, H. Goldberg, M. Brown,
// “Open Source Geospatial Tools to Enable Large Scale 3D Scene Modeling,”
// "Open Source Geospatial Tools to Enable Large Scale 3D Scene Modeling,"
// FOSS4G, 2017.
#include <algorithm>
#include <iterator>
#include <boost/algorithm/string/trim.hpp>
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/tokenizer.hpp>
#include "orthoimage.h"
#include "plugin.hpp"
#include "shr3d.h"
#include <pdal/pdal_macros.hpp>
namespace pdal
{
static PluginInfo const s_info = PluginInfo(
"writers.shr3d", "Shareable High Resolution 3D (Almes et al., 2017)",
"http://pdal.io/stages/writers.shr3d.html");
static StaticPluginInfo const s_info = {
"writers.shr3d",
"Shareable High Resolution 3D (Almes et al., 2017)",
"http://pdal.io/stages/writers.shr3d.html"
};
CREATE_SHARED_PLUGIN(1, 0, Shr3dWriter, Writer, s_info)
CREATE_SHARED_STAGE(Shr3dWriter, s_info)
Shr3dWriter::Shr3dWriter() : shr3dder(), m_filename() {
output_types = (boost::format("%d,%d,%d,%d,%d,%d,%d")
% shr3d::DSM
% shr3d::MIN
% shr3d::DTM
% shr3d::INTENSITY
% shr3d::CLASS
% shr3d::BUILDING
% shr3d::BUILDING_OUTLINES).str();
}
std::string Shr3dWriter::getName() const
{
@@ -47,39 +64,70 @@ std::string Shr3dWriter::getName() const
void Shr3dWriter::addArgs(ProgramArgs& args)
{
args.add("filename", "Output filename", m_filename).setPositional();
args.add("filename", "Output filename (if multiple image types are specified, will be used as a basename)", m_filename).setPositional();
args.add("dh", "Horizontal uncertainty", shr3dder.dh_meters, shr3dder.dh_meters);
args.add("dz", "Vertical uncertainty", shr3dder.dz_meters, shr3dder.dz_meters);
args.add("agl", "Minimum building height above ground level", shr3dder.agl_meters, shr3dder.agl_meters);
args.add("area", "Minimum building area", shr3dder.min_area_meters, shr3dder.min_area_meters);
args.add("egm96", "Set vertical datum to EGM96", shr3dder.egm96, shr3dder.egm96);
args.add("image_type", "Image Type (0=DSM,1=MIN,2=DTM,3=CLASS,4=BUILDING)", output, output);
args.add("bounds", "Bounds for cropped image, e.g.: '([xmin, xmax], [ymin, ymax])'", shr3dder.bounds);
args.add("image_type",
(boost::format("Output image type (%d=DSM,%d=MIN,%d=DTM,%d=INTENSITY,%d=CLASS,%d=BUILDING,%d=BUILDING SHAPEFILE)")
% shr3d::DSM % shr3d::MIN % shr3d::DTM % shr3d::INTENSITY % shr3d::CLASS % shr3d::BUILDING % shr3d::BUILDING_OUTLINES).str(),
output_types, output_types);
}
std::vector<shr3d::ImageType> getImageTypeList(std::string str) {
std::vector<shr3d::ImageType> output;
//Remove whitespace from string
boost::algorithm::trim(str);
// Parse and convert
boost::tokenizer<> tok(str);
transform(tok.begin(), tok.end(), std::back_inserter(output),
[](std::string s) -> shr3d::ImageType {
return (shr3d::ImageType) boost::lexical_cast<unsigned int>(s);
});
return output;
}
void Shr3dWriter::write(const PointViewPtr view)
{
if (output < shr3d::DSM || output > shr3d::LABEL)
throw pdal_error("Unknown output image type\n");
// Set output filenames
std::map<shr3d::ImageType,std::string> outputFilenames;
outputFilenames[(shr3d::ImageType) output] = m_filename;
std::vector<shr3d::ImageType> types = getImageTypeList(output_types);
if (types.empty()) {
throw pdal_error("No image type specified");
} else if (types.size() == 1) {
outputFilenames[types.front()] = m_filename;
} else {
for (shr3d::ImageType t : types) {
switch (t) {
case shr3d::DSM: outputFilenames[t] = m_filename + "_DSM.tif"; break;
case shr3d::MIN: outputFilenames[t] = m_filename + "_MIN.tif"; break;
case shr3d::DTM: outputFilenames[t] = m_filename + "_DTM.tif"; break;
case shr3d::INTENSITY: outputFilenames[t] = m_filename + "_INT.tif"; break;
case shr3d::CLASS: outputFilenames[t] = m_filename + "_class.tif"; break;
case shr3d::BUILDING: outputFilenames[t] = m_filename + "_buildings.tif"; break;
case shr3d::BUILDING_OUTLINES: outputFilenames[t] = m_filename + "_buildings.shp"; break;
case shr3d::MINAGL: outputFilenames[t] = m_filename + "_MINAGL.tif"; break;
default: throw pdal_error("Unknown image type");
}
}
}
// Read in point cloud
shr3d::PointCloud pset;
if (!pset.Read(view))
if (!shr3dder.setPSET(view))
throw pdal_error("Error reading from view\n");
// Create DSM
shr3d::OrthoImage<unsigned short> dsmImage;
if (!shr3dder.createDSM(pset, dsmImage))
throw pdal_error("Error creating DSM\n");
// Run process
shr3dder.process(outputFilenames);
// Create MIN
shr3d::OrthoImage<unsigned short> minImage;
if (!shr3dder.createMIN(pset,minImage))
throw pdal_error("Error creating minimum Z image\n");
// Create everything else
shr3dder.process(dsmImage, minImage, outputFilenames);
// Reset for subsequent uses
shr3dder.reset();
}
} // namespace pdal

4
src/shr3d/plugin.hpp
View File

@@ -40,7 +40,7 @@ namespace pdal
class PDAL_DLL Shr3dWriter : public Writer
{
public:
Shr3dWriter() : shr3dder(), m_filename(), output(shr3d::DTM) {}
Shr3dWriter();
static void* create();
static int32_t destroy(void*);
@@ -53,7 +53,7 @@ private:
shr3d::Shr3dder shr3dder;
std::string m_filename;
int output;
std::string output_types;
Shr3dWriter& operator=(const Shr3dWriter&) = delete;
Shr3dWriter(const Shr3dWriter&) = delete;

357
src/shr3d/shr3d.cpp
View File

@@ -12,88 +12,64 @@
#include <climits>
#include <vector>
#include "geo_polygon.h"
using namespace shr3d;
void Shr3dder::process(const OrthoImage<unsigned short> &dsmImage, const OrthoImage<unsigned short> &minImage,
std::map<ImageType,std::string> outputFilenames) {
if (outputFilenames.empty())
return;
if (!outputFilenames[DSM].empty())
dsmImage.write(outputFilenames[DSM].c_str(), true, egm96);
if (!outputFilenames[MIN].empty())
minImage.write(outputFilenames[MIN].c_str(), true, egm96);
ImageType last_output = DSM;
for (const std::pair<ImageType,std::string>& o : outputFilenames)
if (o.first > last_output && o.first <= BUILDING_OUTLINES)
last_output = o.first;
if (last_output <= MIN)
return;
shr3d::OrthoImage<unsigned short> dsm2Image, dtmImage;
shr3d::OrthoImage<unsigned long> labelImage;
createDTM(dsmImage,minImage,dtmImage,dsm2Image,labelImage);
if (!outputFilenames[DSM2].empty())
dsm2Image.write(outputFilenames[DSM2].c_str(), true, egm96);
if (!outputFilenames[DTM].empty())
dtmImage.write(outputFilenames[DTM].c_str(), true, egm96);
if (!outputFilenames[LABEL].empty())
labelImage.write(outputFilenames[LABEL].c_str(), false);
if (last_output <= DTM)
return;
shr3d::OrthoImage<unsigned char> classImage = labelClasses(dsmImage,dtmImage,dsm2Image,labelImage);
if (!outputFilenames[CLASS].empty())
classImage.write(outputFilenames[CLASS].c_str(), false);
if (last_output <= CLASS)
return;
shr3d::OrthoImage<unsigned char> bldgImage = labelBuildings(classImage);
if (!outputFilenames[BUILDING].empty())
bldgImage.write(outputFilenames[BUILDING].c_str(), false);
if (last_output <= BUILDING)
return;
shr3d::OrthoImage<int> bldgLabels(&bldgImage);
bldgLabels.labelConnectedComponentsFrom(&bldgImage);
if (!outputFilenames[LABELED_BUILDINGS].empty())
bldgLabels.write(outputFilenames[LABELED_BUILDINGS].c_str(), false);
shr3d::OrthoImage<int> bldgLabels3;
bldgLabels3.nn_upsample(&bldgLabels,3);
if (!outputFilenames[LABELED_BUILDINGS_3].empty())
bldgLabels3.write(outputFilenames[LABELED_BUILDINGS_3].c_str(), false);
std::map<int,GeoPolygon<double>> bounds = GeoPolygon<double>::traceBoundaries(bldgLabels3);
printf("Traced %lu building outlines.\n",bounds.size());
std::map<int,GeoPolygon<double>> new_bounds;
for (std::pair<int,GeoPolygon<double>> pr : bounds) {
GeoPolygon<double> poly = pr.second.buildingSimplify();
if (poly.ring.empty())
void Shr3dder::process(std::map<ImageType,std::string> outputFilenames) {
for (std::pair<ImageType,std::string> elem : outputFilenames) {
if (elem.second.empty())
continue;
new_bounds[pr.first] = poly;
auto fname = elem.second.c_str();
bool b = true;
switch (elem.first) {
case DSM:
b = getDSM().write(fname, true, egm96); break;
case MIN:
b = getMIN().write(fname, true, egm96); break;
case DTM:
b = getDTM().write(fname, true, egm96); break;
case INTENSITY:
b = getINT().write(fname, false); break;
case CLASS:
b = getCLS().write(fname, false); break;
case BUILDING:
b = getBLDG().write(fname, false); break;
case BUILDING_OUTLINES:
b = GeoPolygon<double>::write(fname, getBLDGPOLY()); break;
case DSM2:
b = getDSM2().write(fname, true, egm96); break;
case MINAGL:
b = getMINAGL().write(fname, true, egm96); break;
case LABEL:
b = getLBL().write(fname, false); break;
case LABELED_BUILDINGS:
b = getBLDGLBL().write(fname, false); break;
case LABELED_BUILDINGS_3:
b = getBLDGLBL3().write(fname, false); break;
default:
throw pdal::pdal_error("Unknown output image type\n");
}
if (!b)
std::cerr << "Failed to write " << elem.second << std::endl;
}
printf("After simplifying, left with %lu building outlines.\n",new_bounds.size());
if (!outputFilenames[BUILDING_OUTLINES].empty())
GeoPolygon<double>::write(outputFilenames[BUILDING_OUTLINES],new_bounds);
}
bool Shr3dder::createDSM(const PointCloud& pset, OrthoImage<unsigned short> &dsmImage) {
void sizeImageFromBox(OrthoImage<unsigned short> &im, pdal::BOX2D box, float dh_meters) {
if (!box.empty()) {
im.easting = box.minx;
im.northing = box.miny;
im.width = (unsigned int) ((box.maxx - box.minx) / dh_meters);
im.height = (unsigned int) ((box.maxy - box.miny) / dh_meters);
im.Allocate(im.width, im.height);
}
}
void Shr3dder::createDSM() {
printf("Creating DSM\n");
if (!pset.numPoints)
throw pdal::pdal_error("Point cloud undefined / empty\n");
sizeImageFromBox(dsmImage, bounds.to2d(), dh_meters);
if (!dsmImage.readFromPointCloud(pset, (float) dh_meters, shr3d::MAX_VALUE))
return false;
throw pdal::pdal_error("Error filling DSM from point cloud\n");
// Median filter, replacing only points differing by more than the AGL threshold.
dsmImage.quantileFilter(1, (unsigned int) (agl_meters / dsmImage.scale), 0.4);
@@ -107,26 +83,35 @@ bool Shr3dder::createDSM(const PointCloud& pset, OrthoImage<unsigned short> &dsm
// Fill small voids in the DSM.
dsmImage.fillVoidsPyramid(true, 2);
return true;
}
bool Shr3dder::createMIN(const PointCloud& pset, OrthoImage<unsigned short> &minImage) {
// Now get the minimum Z values for the DTM.
if (!minImage.readFromPointCloud(pset, (float) dh_meters, shr3d::MIN_VALUE))
return false;
void Shr3dder::createMIN() {
printf("Creating MIN\n");
if (pset.numPoints > 0) {
sizeImageFromBox(minImage, bounds.to2d(), dh_meters);
// Now get the minimum Z values for the DTM.
if (!minImage.readFromPointCloud(pset, (float) dh_meters, shr3d::MIN_VALUE))
throw pdal::pdal_error("Error filling MIN from point cloud\n");
// Median filter, replacing only points differing by more than the AGL threshold.
minImage.quantileFilter(2, (unsigned int) (agl_meters / minImage.scale), 0.33);
// Median filter, replacing only points differing by more than the AGL threshold.
minImage.quantileFilter(2, (unsigned int) (agl_meters / minImage.scale), 0.33);
} else {
// Copy DSM
minImage = getDSM();
// Min filter, replacing only points differing by more than the AGL threshold.
minImage.minFilter(4, (unsigned int) (agl_meters / minImage.scale));
}
// Fill small voids in the DSM.
minImage.fillVoidsPyramid(true, 2);
return true;
}
void Shr3dder::createDTM(const OrthoImage<unsigned short> &dsmImage, const OrthoImage<unsigned short> &minImage,
OrthoImage<unsigned short> &dtmImage, OrthoImage<unsigned short> &dsm2Image, OrthoImage<unsigned long> &labelImage) {
void Shr3dder::createDTM() {
const OrthoImage<unsigned short> &dsmImage = getDSM();
const OrthoImage<unsigned short> &minImage = getMIN();
printf("Creating DTM\n");
// Convert horizontal and vertical uncertainty values to bin units.
int dh_bins = MAX(1, (int) floor(dh_meters / dsmImage.gsd));
@@ -153,7 +138,7 @@ void Shr3dder::createDTM(const OrthoImage<unsigned short> &dsmImage, const Ortho
// Generate label image.
labelImage = OrthoImage<unsigned long>(&dsmImage);
labelImage = OrthoImage<unsigned int>(&dsmImage);
// Allocate a DTM image as SHORT and copy in the Min DSM values.
dtmImage = minImage;
@@ -181,9 +166,133 @@ void Shr3dder::createDTM(const OrthoImage<unsigned short> &dsmImage, const Ortho
dtmImage.fillVoidsPyramid(true, 2);
}
OrthoImage<unsigned char> Shr3dder::labelClasses(
const OrthoImage<unsigned short> &dsmImage, const OrthoImage<unsigned short> &dtmImage,
const OrthoImage<unsigned short> &dsm2Image, const OrthoImage<unsigned long> &labelImage) {
void Shr3dder::createIntensity() {
const OrthoImage<unsigned short> &dsmImage = getDSM();
printf("Creating INT\n");
if (!pset.numPoints)
throw pdal::pdal_error("Point cloud undefined / empty\n");
intImage = OrthoImage<unsigned short>(&dsmImage);
if (pset.numPoints) {
OrthoImage<double> sumImage = OrthoImage<double>(&dsmImage);
OrthoImage<unsigned short> cntImage = OrthoImage<unsigned short>(&dsmImage);
// Determine scalings
double mx = 1.0/dsmImage.gsd;
double bx = mx*(pset.xOff - dsmImage.easting) - 0.5;
double my =-1.0/dsmImage.gsd;
double by = dsmImage.height-1 + my*(pset.yOff - dsmImage.northing) + 0.5;
double mz = 1/dsmImage.scale;
double bz = mz*(pset.zOff - dsmImage.offset);
int dz_short = (unsigned int) (dz_meters / dsmImage.scale);
// Loop over points
double max_i = 0;
for (unsigned long i = 0; i < pset.numPoints; i++) {
int x = mx*pset.x(i) + bx;
int y = my*pset.y(i) + by;
int z = mz*pset.z(i) + bz;
for (int y1 = std::max(y,0); y1 <= std::min(y+1,(int) dsmImage.height-1); ++y1) {
for (int x1 = std::max(x,0); x1 <= std::min(x+1,(int) dsmImage.width-1); ++x1) {
int z0 = dsmImage.data[y1][x1];
if (abs(z-z0) < dz_short) {
sumImage.data[y1][x1] += pset.i(i);
++cntImage.data[y1][x1];
max_i = std::max(max_i, sumImage.data[y1][x1] / cntImage.data[y1][x1]);
}
}
}
}
// Divide & copy
double scale = std::numeric_limits<unsigned short>::max() / max_i;
for (unsigned int y = 0; y < intImage.height; ++y) {
for (unsigned int x = 0; x <= intImage.width; ++x) {
intImage.data[y][x] = (unsigned short) (scale * sumImage.data[y][x] / cntImage.data[y][x]);
}
}
// Fill small voids in the INT.
intImage.fillVoidsPyramid(true, 2);
}
}
void Shr3dder::createMinAGL() {
const OrthoImage<unsigned short> &dsmImage = getDSM();
const OrthoImage<unsigned short> &dtmImage = getDTM();
printf("Creating MinAGL\n");
if (!pset.numPoints)
throw pdal::pdal_error("Point cloud undefined / empty\n");
unsigned int dz_short = (unsigned int) (dz_meters / dsmImage.scale);
unsigned int agl_short = (unsigned int) (agl_meters / dsmImage.scale);
// Create the Minimum AGL Image; image gets the minimum values that are greater than the DTM+AGL
if (pset.numPoints) {
// Initialize size and allocate zeros
minAglImage = OrthoImage<unsigned short>(&dsmImage);
// Determine scalings
double mx = 1.0/minAglImage.gsd;
double bx = mx*(pset.xOff - minAglImage.easting) - 0.5;
double my =-1.0/minAglImage.gsd;
double by = minAglImage.height-1 + my*(pset.yOff - minAglImage.northing) + 0.5;
double mz = 1/minAglImage.scale;
double bz = mz*(pset.zOff - minAglImage.offset);
// Loop over points
for (unsigned long i = 0; i < pset.numPoints; i++) {
int x = mx*pset.x(i) + bx;
int y = my*pset.y(i) + by;
unsigned short z = mz*pset.z(i) + bz;
for (int y1 = std::max(y,0); y1 <= std::min(y+1,(int) minAglImage.height-1); ++y1) {
for (int x1 = std::max(x,0); x1 <= std::min(x+1,(int) minAglImage.width-1); ++x1) {
unsigned short& z0 = minAglImage.data[y1][x1];
unsigned short& z1 = dtmImage.data[y1][x1];
if (z1 && (z > z1+agl_short) && (!z0 || z < z0))
z0 = z;
}
}
}
} else {
// For when the input is not a pointset
minAglImage = dsmImage;
for (unsigned int y=0; y<minAglImage.height; ++y)
for (unsigned int x=0; x<minAglImage.width; ++x)
if (minAglImage.data[y][x] <= dtmImage.data[y][x]+agl_short)
minAglImage.data[y][x] = 0;
}
// Median filter, replacing only points differing by more than the AGL threshold.
//minAglImage.quantileFilter(1, agl_short, 0.4, 0, false);
// Filter quantile
float quantile = 0.5f;
// Apply filter to the image
minAglImage.filter([&](unsigned short* val, const unsigned short& ref, std::vector<unsigned short> &ngbrs) {
// Find quantile
size_t ix = std::min((size_t) floor(quantile * ngbrs.size()),ngbrs.size()-1);
std::partial_sort(ngbrs.begin(), ngbrs.begin() + (ix + 1), ngbrs.end());
short qValue = static_cast<short>(ngbrs[ix]);
// Only replace if it differs by more than dz from the median
if (ref && abs(qValue - static_cast<short>(ref)) > dz_short)
*val = qValue;
},2);
// Fill small voids in the image.
//minAglImage.fillVoidsPyramid(true, 2);
}
void Shr3dder::labelClasses() {
const OrthoImage<unsigned short> &dsmImage = getDSM();
const OrthoImage<unsigned short> &dtmImage = getDTM();
const OrthoImage<unsigned short> &dsm2Image = getDSM2();
const OrthoImage<unsigned int> &labelImage = getLBL();
printf("Creating CLS\n");
// Produce a classification raster image with LAS standard point classes.
// Convert horizontal and vertical uncertainty values to bin units.
@@ -192,7 +301,7 @@ OrthoImage<unsigned char> Shr3dder::labelClasses(
unsigned int dz_short = (unsigned int) (dz_meters / dsmImage.scale);
unsigned int agl_short = (unsigned int) (agl_meters / dsmImage.scale);
OrthoImage<unsigned char> classImage(&dsmImage);
classImage = OrthoImage<unsigned char>(&dsmImage);
for (unsigned int j = 0; j < classImage.height; j++) {
for (unsigned int i = 0; i < classImage.width; i++) {
// Set default as unlabeled.
@@ -237,23 +346,46 @@ OrthoImage<unsigned char> Shr3dder::labelClasses(
// Fill missing labels inside building regions.
fillInsideBuildings(classImage);
return classImage;
}
OrthoImage<unsigned char> Shr3dder::labelBuildings(const OrthoImage<unsigned char> &classImage) {
OrthoImage<unsigned char> bldgImage(classImage);
void Shr3dder::labelBuildings() {
const OrthoImage<unsigned char>& classImage = getCLS();
printf("Labeling buildings\n");
bldgImage = OrthoImage<unsigned char>(classImage);
for (unsigned int j = 0; j < classImage.height; j++) {
for (unsigned int i = 0; i < classImage.width; i++) {
if (classImage.data[j][i] != LAS_BUILDING) bldgImage.data[j][i] = 0;
}
}
return bldgImage;
}
void Shr3dder::createOutlines() {
printf("Creating building outlines\n");
const shr3d::OrthoImage<unsigned char>&bldgImage = getBLDG();
bldgLabels = shr3d::OrthoImage<int>(&bldgImage);
bldgLabels.labelConnectedComponentsFrom(&bldgImage);
bldgLabels3.nn_upsample(&bldgLabels,3);
std::map<int,GeoPolygon<double>> bounds = GeoPolygon<double>::traceBoundaries(bldgLabels3);
printf("Traced %lu building outlines.\n",bounds.size());
bldgOutlines.clear();
for (std::pair<int,GeoPolygon<double>> pr : bounds) {
GeoPolygon<double> poly = pr.second.buildingSimplify();
if (poly.ring.empty())
continue;
bldgOutlines[pr.first] = poly;
}
printf("After simplifying, left with %lu building outlines.\n",bldgOutlines.size());
bldgProcessed = true;
}
// Extend object boundaries to capture points missed around the edges.
void extendObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned long> &labelImage,
void extendObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned int> &labelImage,
int edgeResolution, unsigned int minDistanceShortValue) {
// Loop enough to capture the edge resolution.
for (int k = 0; k < edgeResolution; k++) {
@@ -309,7 +441,7 @@ void extendObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<uns
}
// Label boundaries of objects above ground level.
void labelObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned long> &labelImage,
void labelObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned int> &labelImage,
int edgeResolution, unsigned int minDistanceShortValue) {
// Initialize the labels to LABEL_GROUND.
for (unsigned int j = 0; j < labelImage.height; j++) {
@@ -343,7 +475,7 @@ void labelObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsi
}
}
bool findObjectBoundsInColumn(OrthoImage<unsigned long> &labelImage, ObjectType &obj, unsigned int column, int &min_row, int &max_row) {
bool findObjectBoundsInColumn(OrthoImage<unsigned int> &labelImage, ObjectType &obj, unsigned int column, int &min_row, int &max_row) {
min_row = -1;
max_row = -1;
@@ -373,7 +505,7 @@ bool findObjectBoundsInColumn(OrthoImage<unsigned long> &labelImage, ObjectType
}
// Fill inside the object countour labels if points are above the nearby ground level.
void fillObjectBounds(OrthoImage<unsigned long> &newLabelImage, OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned short> &dsmImage, ObjectType &obj,
void fillObjectBounds(OrthoImage<unsigned int> &newLabelImage, OrthoImage<unsigned int> &labelImage, OrthoImage<unsigned short> &dsmImage, ObjectType &obj,
int edgeResolution) {
unsigned int label = obj.label;
@@ -413,8 +545,8 @@ void fillObjectBounds(OrthoImage<unsigned long> &newLabelImage, OrthoImage<unsig
}
}
// If entire column is labeled, then continue.
if ((startIndex == 0) && (stopIndex == (int) labelImage.width-1)) continue;
// If entire column is labeled, then continue.
if ((startIndex == 0) && (stopIndex == (int) labelImage.width-1)) continue;
// Get max ground level height for this row.
// If the DSM image value is void, then the ground level value is zero, so that's ok.
@@ -443,8 +575,8 @@ void fillObjectBounds(OrthoImage<unsigned long> &newLabelImage, OrthoImage<unsig
if (!findObjectBoundsInColumn(labelImage, obj, i, startIndex, stopIndex))
continue;
// If entire row is labeled, then continue.
if ((startIndex == 0) && (stopIndex == (int) labelImage.height-1)) continue;
// If entire row is labeled, then continue.
if ((startIndex == 0) && (stopIndex == (int) labelImage.height-1)) continue;
// Get max ground level height for this row.
unsigned short groundLevel;
@@ -508,8 +640,8 @@ void fillObjectBounds(OrthoImage<unsigned long> &newLabelImage, OrthoImage<unsig
}
// Add neighboring pixels to an object.
bool addNeighbors(std::vector<PixelType> &neighbors, OrthoImage<unsigned long> &labelImage,
OrthoImage<unsigned short> &dsmImage, ObjectType &obj, unsigned int dzShort) {
bool addNeighbors(std::vector<PixelType> &neighbors, OrthoImage<unsigned int> &labelImage,
const OrthoImage<unsigned short> &dsmImage, ObjectType &obj, unsigned int dzShort) {
// Get neighbors for all pixels in the list.
std::vector<PixelType> newNeighbors;
for (size_t k = 0; k < neighbors.size(); k++) {
@@ -549,7 +681,7 @@ bool addNeighbors(std::vector<PixelType> &neighbors, OrthoImage<unsigned long> &
}
// Group connected labeled pixels into objects.
void groupObjects(OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned short> &dsmImage,
void groupObjects(OrthoImage<unsigned int> &labelImage, const OrthoImage<unsigned short> &dsmImage,
std::vector<ObjectType> &objects, long maxCount, unsigned int dzShort) {
// Sweep from top left to bottom right, assigning object labels.
long maxGroupSize = 0;
@@ -606,7 +738,7 @@ void groupObjects(OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned sho
// Finish label image for display as image overlay in QT Reader.
// Set all labeled values to 1. Leave all unlabeled values LABEL_GROUND.
void finishLabelImage(OrthoImage<unsigned long> &labelImage) {
void finishLabelImage(OrthoImage<unsigned int> &labelImage) {
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
if ((labelImage.data[j][i] != LABEL_GROUND)) {
@@ -617,7 +749,7 @@ void finishLabelImage(OrthoImage<unsigned long> &labelImage) {
}
// Classify ground points, fill the voids, and generate a bare earth terrain model.
void Shr3dder::classifyGround(OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned short> &dsmImage,
void Shr3dder::classifyGround(OrthoImage<unsigned int> &labelImage, OrthoImage<unsigned short> &dsmImage,
OrthoImage<unsigned short> &dtmImage, int dhBins, unsigned int dzShort) {
// Fill voids.
printf("Filling voids...\n");
@@ -650,7 +782,7 @@ void Shr3dder::classifyGround(OrthoImage<unsigned long> &labelImage, OrthoImage<
// Generate object groups and void fill them in the DEM image.
printf("Labeling and removing objects...\n");
OrthoImage<unsigned long> newlabelImage(labelImage);
OrthoImage<unsigned int> newlabelImage(labelImage);
for (size_t i = 0; i < objects.size(); i++) {
fillObjectBounds(newlabelImage,labelImage, dtmImage, objects[i], dhBins);
}
@@ -734,7 +866,7 @@ void Shr3dder::classifyGround(OrthoImage<unsigned long> &labelImage, OrthoImage<
// Classify non-ground points.
void Shr3dder::classifyNonGround(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned short> &dtmImage,
OrthoImage<unsigned long> &labelImage, unsigned int dzShort, unsigned int aglShort,
OrthoImage<unsigned int> &labelImage, unsigned int dzShort, unsigned int aglShort,
float minAreaMeters) {
// Compute minimum number of points based on threshold given for area.
// Note that ISPRS challenges indicate that performance is dramatically better for structures larger than 50m area.
@@ -933,7 +1065,6 @@ addClassNeighbors(std::vector<PixelType> &neighbors, OrthoImage<unsigned char> &
for (size_t k = 0; k < newNeighbors.size(); k++) {
neighbors.push_back(newNeighbors[k]);
}
// neighbors = newNeighbors;
return true;
} else return false;
}

90
src/shr3d/shr3d.h
View File

@@ -10,7 +10,12 @@
#include <cstdio>
#include <cstring>
#include <map>
#include <memory>
#include <string>
#include <pdal/util/Bounds.hpp>
#include "geo_polygon.h"
#include "orthoimage.h"
#define LABEL_OBJECT 1
@@ -64,10 +69,12 @@ namespace shr3d {
DSM,
MIN,
DTM,
INTENSITY,
CLASS,
BUILDING,
BUILDING_OUTLINES,
DSM2,
MINAGL,
LABEL,
LABELED_BUILDINGS,
LABELED_BUILDINGS_3,
@@ -82,32 +89,91 @@ namespace shr3d {
double min_area_meters;
double max_tree_height_meters;
bool egm96;
pdal::Bounds bounds;
// Constructor
Shr3dder() : dh_meters(1), dz_meters(1), agl_meters(2),
min_area_meters(50), max_tree_height_meters(40), egm96(false) {}
min_area_meters(50), max_tree_height_meters(40), egm96(false), bldgProcessed(false) {}
// Function declarations.
void process(const OrthoImage<unsigned short> &dsmImage, const OrthoImage<unsigned short> &minImage,
std::map<ImageType,std::string> outputFilenames);
void process(std::map<ImageType,std::string> outputFilenames);
bool createDSM(const PointCloud& pset, OrthoImage<unsigned short> &dsmImage);
bool setDSM(std::string dsmFile) { return dsmImage.read(dsmFile.c_str()); }
bool setPSET(std::string psetFile) { return pset.Read(psetFile.c_str()); }
bool setPSET(const pdal::PointViewPtr view) { return pset.Read(view); }
void reset() {
dsmImage.Deallocate();
minImage.Deallocate();
dtmImage.Deallocate();
intImage.Deallocate();
classImage.Deallocate();
bldgImage.Deallocate();
dsm2Image.Deallocate();
minAglImage.Deallocate();
labelImage.Deallocate();
bldgLabels.Deallocate();
bldgLabels3.Deallocate();
bldgOutlines.clear();
bldgProcessed = false;
}
bool createMIN(const PointCloud& pset, OrthoImage<unsigned short> &minImage);
#define IMGET(im,func) if(im.empty()) func(); return im
const OrthoImage<unsigned short> &getDSM() {IMGET(dsmImage,createDSM);}
const OrthoImage<unsigned short> &getMIN() {IMGET(minImage,createMIN);}
const OrthoImage<unsigned short> &getINT() {IMGET(intImage,createIntensity);}
const OrthoImage<unsigned short> &getDTM() {IMGET(dtmImage,createDTM);}
const OrthoImage<unsigned short> &getDSM2() {IMGET(dsm2Image,createDTM);}
const OrthoImage<unsigned int> &getLBL() {IMGET(labelImage,createDTM);}
const OrthoImage<unsigned short> &getMINAGL() {IMGET(minAglImage,createMinAGL);}
const OrthoImage<unsigned char> &getCLS() {IMGET(classImage,labelClasses);}
const OrthoImage<unsigned char> &getBLDG() {IMGET(bldgImage,labelBuildings);}
const OrthoImage<int> &getBLDGLBL() {IMGET(bldgLabels,createOutlines);}
const OrthoImage<int> &getBLDGLBL3(){IMGET(bldgLabels3,createOutlines);}
const std::map<int,GeoPolygon<double>>& getBLDGPOLY() {
if (!bldgProcessed)
createOutlines();
return bldgOutlines;
}
#undef IMGET
protected:
PointCloud pset;
OrthoImage<unsigned short> dsmImage;
OrthoImage<unsigned short> minImage;
OrthoImage<unsigned short> dtmImage;
OrthoImage<unsigned short> intImage;
OrthoImage<unsigned char> classImage;
OrthoImage<unsigned char> bldgImage;
OrthoImage<unsigned short> dsm2Image;
OrthoImage<unsigned short> minAglImage;
OrthoImage<unsigned int> labelImage;
OrthoImage<int> bldgLabels;
OrthoImage<int> bldgLabels3;
std::map<int,GeoPolygon<double>> bldgOutlines;
bool bldgProcessed;
void createDSM();
void createDTM(const OrthoImage<unsigned short> &dsmImage, const OrthoImage<unsigned short> &minImage,
OrthoImage<unsigned short> &dtmImage, OrthoImage<unsigned short> &dsm2Image, OrthoImage<unsigned long> &labelImage);
void createMIN();
OrthoImage<unsigned char> labelClasses(const OrthoImage<unsigned short> &dsmImage, const OrthoImage<unsigned short> &dtmImage,
const OrthoImage<unsigned short> &dsm2Image, const OrthoImage<unsigned long> &labelImage);
void createDTM();
OrthoImage<unsigned char> labelBuildings(const OrthoImage<unsigned char> &classImage);
void createMinAGL();
void classifyGround(OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned short> &dsmImage,
void createIntensity();
void labelClasses();
void labelBuildings();
void createOutlines();
void classifyGround(OrthoImage<unsigned int> &labelImage, OrthoImage<unsigned short> &dsmImage,
OrthoImage<unsigned short> &dtmImage, int dhBins, unsigned int dzShort);
void classifyNonGround(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned short> &dtmImage,
OrthoImage<unsigned long> &labelImage, unsigned int dzShort,
OrthoImage<unsigned int> &labelImage, unsigned int dzShort,
unsigned int aglShort,
float minAreaMeters);

6
windowsInstallation.txt
View File

@@ -3,13 +3,15 @@ While only the docker version of align3d and shr3d is fully supported, it is pos
Download and run OSGeo4W64 (64-bit) Installer (http://download.osgeo.org/osgeo4w/osgeo4w-setup-x86_64.exe)
Select advanced installation. Everything else until packages is default.
From Command Line Utilities select:
Gdal, gdal-dev, pdal
Gdal, gdal-dev, pdal, boost
Install packages to meet dependencies (page after package selection).
If you did not install OSGeo4W64 in the standard location (C:\OSGeo4W64), do the following:
In CMakeLists.txt, modify SET(CMAKE_PREFIX_PATH "...") to point to your OSGeo4W64 install directory
Edit ...\OSGeo4W64\lib\pdal\cmake\PDALConfig.cmake, ...\OSGeo4W64\lib\pdal\cmake\PDALTargets.cmake, and ...\OSGeo4W64\lib\pdal\cmake\PDALTargets-relwithdebinfo.cmake so that all references to C:\OSGeo4W64 are replaced with the actual install location.
In CMakeLists.txt, modify SET(BOOST_INCLUDEDIR "C:\\OSGeo4W64\\include\\boost-1_55") as needed to point to the include directory for the version of boost installed with OSGeo4W64.
Download and install visual studio 2017 community (REQUIRES SYSTEM-RESTART)
Install packages:
Desktop development with c++
@@ -24,7 +26,7 @@ Open Visual Studio
Verify shr3d.exe and align3d.exe exist in the pubgeo/build/x64-Debug folder
Open command prompt
Run ...\OSGeo4W64bin\o4w_env.bat; this will configure required environment variables for the session
Run ...\OSGeo4W64\bin\o4w_env.bat; this will configure required environment variables for the session
cd to ...\pubgeo\build\x64-Debug
Run align3d or shr3d:
align3d.exe ...\path\to\reference.las ...\path\to\alignMe.las gsd=1.0 maxt=10.0