github/pubgeo
1
1
Fork 0
mirror of https://github.com/pubgeo/pubgeo.git synced 2026-01-06 20:53:51 -05:00
Code Issues Packages Projects Releases Wiki Activity

Add PDAL plugin

Browse Source 
SHR3D is available as writers.shr3d and can be appended to any PDAL pipeline.
PDAL will not be altered to infer writers.shr3d (many writers could write a
TIF), so users will have to specify writers.shr3d explicitly.

ALIGN3D is available as filters.align3d and can be inserted into any PDAL
pipeline. ALIGN3D will take the first incoming PointView as the fixed or
reference PointView, and will register all subsequent PointViews to this
reference view.

Minor changes: Ignore build directory for in-source builds. Add clang-format
config file (the plugin code adheres to this).

PubGeo common code updated to read directly from PDAL PointViews. Common code
also updated to silence a bunch of warnings.
This commit is contained in:
Bradley J Chambers
2017-09-18 12:16:23 -04:00
parent f38f4bfb1f
commit fef97f6275
13 changed files with 718 additions and 86 deletions
Download Patch File Download Diff File Expand all files Collapse all files

97
.clang-format Normal file
View File

@@ -0,0 +1,97 @@
---
Language: Cpp
# BasedOnStyle: LLVM
AccessModifierOffset: -4
AlignAfterOpenBracket: Align
AlignConsecutiveAssignments: false
AlignConsecutiveDeclarations: false
AlignEscapedNewlinesLeft: false
AlignOperands: true
AlignTrailingComments: true
AllowAllParametersOfDeclarationOnNextLine: true
AllowShortBlocksOnASingleLine: false
AllowShortCaseLabelsOnASingleLine: false
AllowShortFunctionsOnASingleLine: Empty
AllowShortIfStatementsOnASingleLine: false
AllowShortLoopsOnASingleLine: false
AlwaysBreakAfterDefinitionReturnType: None
AlwaysBreakAfterReturnType: None
AlwaysBreakBeforeMultilineStrings: false
AlwaysBreakTemplateDeclarations: false
BinPackArguments: true
BinPackParameters: true
BraceWrapping:
AfterClass: false
AfterControlStatement: false
AfterEnum: false
AfterFunction: false
AfterNamespace: false
AfterObjCDeclaration: false
AfterStruct: false
AfterUnion: false
BeforeCatch: false
BeforeElse: false
IndentBraces: false
BreakBeforeBinaryOperators: None
BreakBeforeBraces: Allman
BreakBeforeTernaryOperators: true
BreakConstructorInitializersBeforeComma: false
BreakAfterJavaFieldAnnotations: false
BreakStringLiterals: true
ColumnLimit: 80
CommentPragmas: '^ IWYU pragma:'
BreakBeforeInheritanceComma: false
ConstructorInitializerAllOnOneLineOrOnePerLine: false
ConstructorInitializerIndentWidth: 4
ContinuationIndentWidth: 4
Cpp11BracedListStyle: true
DerivePointerAlignment: false
DisableFormat: false
ExperimentalAutoDetectBinPacking: false
FixNamespaceComments: true
ForEachMacros: [ foreach, Q_FOREACH, BOOST_FOREACH ]
IncludeCategories:
- Regex: '^"(llvm|llvm-c|clang|clang-c)/'
Priority: 2
- Regex: '^(<|"(gtest|isl|json)/)'
Priority: 3
- Regex: '.*'
Priority: 1
IncludeIsMainRegex: '$'
IndentCaseLabels: false
IndentWidth: 4
IndentWrappedFunctionNames: false
JavaScriptQuotes: Leave
JavaScriptWrapImports: true
KeepEmptyLinesAtTheStartOfBlocks: true
MacroBlockBegin: ''
MacroBlockEnd: ''
MaxEmptyLinesToKeep: 1
NamespaceIndentation: None
ObjCBlockIndentWidth: 2
ObjCSpaceAfterProperty: false
ObjCSpaceBeforeProtocolList: true
PenaltyBreakBeforeFirstCallParameter: 19
PenaltyBreakComment: 300
PenaltyBreakFirstLessLess: 120
PenaltyBreakString: 1000
PenaltyExcessCharacter: 1000000
PenaltyReturnTypeOnItsOwnLine: 60
PointerAlignment: Left
ReflowComments: true
SortIncludes: true
SpaceAfterCStyleCast: false
SpaceAfterTemplateKeyword: true
SpaceBeforeAssignmentOperators: true
SpaceBeforeParens: ControlStatements
SpaceInEmptyParentheses: false
SpacesBeforeTrailingComments: 1
SpacesInAngles: false
SpacesInContainerLiterals: true
SpacesInCStyleCastParentheses: false
SpacesInParentheses: false
SpacesInSquareBrackets: false
Standard: Cpp11
TabWidth: 4
UseTab: Never
...

2
.gitignore vendored
View File

@@ -1,3 +1,5 @@
build
# Prerequisites
*.d

16
src/align3d/CMakeLists.txt
View File

@@ -25,4 +25,18 @@ if(TESTING)
TARGET_LINK_LIBRARIES(align3dTest
PUBLIC
ALIGN3D_LIB)
endif()
endif()
IF(WIN32)
SET(ALIGN3D_FILTER_NAME libpdal_plugin_filter_align3d)
ELSE(WIN32)
SET(ALIGN3D_FILTER_NAME pdal_plugin_filter_align3d)
ENDIF(WIN32)
ADD_LIBRARY(${ALIGN3D_FILTER_NAME} SHARED plugin.cpp)
LINK_DIRECTORIES(${PDAL_LIBRARY_DIRS})
TARGET_LINK_LIBRARIES(${ALIGN3D_FILTER_NAME} ${PDAL_LIBRARIES} ALIGN3D_LIB)
SET_TARGET_PROPERTIES(${ALIGN3D_FILTER_NAME} PROPERTIES SOVERSION "0.1.0")
INSTALL(TARGETS ${ALIGN3D_FILTER_NAME}
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)

140
src/align3d/plugin.cpp Normal file
View File

@@ -0,0 +1,140 @@
// Copyright (c) 2017, Bradley J Chambers (brad.chambers@gmail.com)
// All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// 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,”
// FOSS4G, 2017.
#include <algorithm>
#include "align3d.h"
#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");
CREATE_SHARED_PLUGIN(1, 0, Align3dFilter, Filter, s_info)
std::string Align3dFilter::getName() const
{
return s_info.name;
}
void Align3dFilter::addArgs(ProgramArgs& args)
{
args.add("gsd", "Ground Sample Distance (GSD) for gridding", m_gsd, 1.0);
args.add("maxt", "Maximum XYZ translation in search", m_maxt, 10.0);
args.add("maxdz", "Max local Z difference for matching", m_maxdz, 0.0);
}
PointViewSet Align3dFilter::run(PointViewPtr view)
{
PointViewSet viewSet;
if (!view->size())
return viewSet;
if (!m_fixed)
{
m_fixed = view;
return viewSet;
}
// Default MAXDZ = GSD x 2 to ensure reliable performance on steep
// slopes.
if (m_maxdz == 0.0)
m_maxdz = m_gsd * 2.0;
// Align the target point cloud to the reference.
// AlignTarget2Reference(referenceFileName, targetFileName, params);
// Read the reference LAS file as a DSM.
// Fill small voids.
// Remove points along edges which are difficult to match.
log()->get(LogLevel::Debug) << "Reading reference point cloud\n";
pubgeo::OrthoImage<unsigned short> referenceDSM;
if (!referenceDSM.readFromPointView(m_fixed, m_gsd, pubgeo::MAX_VALUE))
throw pdal_error("Error reading reference PointView\n");
referenceDSM.fillVoidsPyramid(true, 2);
log()->get(LogLevel::Debug) << "Filtering reference point cloud.\n";
referenceDSM.edgeFilter((long)(m_maxdz / referenceDSM.scale));
// Read the target LAS file as a DSM.
// Fill small voids.
// Remove points along edges which are difficult to match.
log()->get(LogLevel::Debug) << "Reading target point cloud\n";
pubgeo::OrthoImage<unsigned short> targetDSM;
if (!targetDSM.readFromPointView(view, m_gsd, pubgeo::MAX_VALUE))
throw pdal_error("Error reading target PointView\n");
targetDSM.fillVoidsPyramid(true, 2);
log()->get(LogLevel::Debug) << "Filtering target point cloud.\n";
targetDSM.edgeFilter((long)(m_maxdz / targetDSM.scale));
// Get overlapping bounds.
align3d::AlignBounds bounds;
bounds.xmin = std::max(referenceDSM.easting, targetDSM.easting);
bounds.ymin = std::max(referenceDSM.northing, targetDSM.northing);
bounds.xmax =
std::min(referenceDSM.easting + (referenceDSM.width * referenceDSM.gsd),
targetDSM.easting + (targetDSM.width * targetDSM.gsd));
bounds.ymax = std::min(
referenceDSM.northing + (referenceDSM.height * referenceDSM.gsd),
targetDSM.northing + (targetDSM.height * targetDSM.gsd));
bounds.width = bounds.xmax - bounds.xmin;
bounds.height = bounds.ymax - bounds.ymin;
double overlap_km = bounds.width / 1000.0 * bounds.height / 1000.0;
log()->get(LogLevel::Debug)
<< "Overlap = " << bounds.width << " m x " << bounds.height
<< " m = " << overlap_km << " km\n";
if (overlap_km == 0.0)
throw pdal_error("PointViews do not overlap");
// Estimate rigid body transform to align target points to reference.
align3d::AlignResult result;
log()->get(LogLevel::Debug) << "Estimating rigid body transformation.\n";
EstimateRigidBody(referenceDSM, targetDSM, m_maxt, bounds, result);
for (PointId i = 0; i < view->size(); ++i)
{
double x = view->getFieldAs<double>(Dimension::Id::X, i);
double y = view->getFieldAs<double>(Dimension::Id::Y, i);
double z = view->getFieldAs<double>(Dimension::Id::Z, i);
view->setField(Dimension::Id::X, i, x + result.tx);
view->setField(Dimension::Id::Y, i, y + result.ty);
view->setField(Dimension::Id::Z, i, z + result.tz);
}
viewSet.insert(view);
MetadataNode root = getMetadata();
root.add("x_offset", result.tx);
root.add("y_offset", result.ty);
root.add("z_offset", result.tz);
root.add("z_rms", result.rms);
return viewSet;
}
} // namespace pdal

61
src/align3d/plugin.hpp Normal file
View File

@@ -0,0 +1,61 @@
// Copyright (c) 2017, Bradley J Chambers (brad.chambers@gmail.com)
// All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// 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,”
// FOSS4G, 2017.
#pragma once
#include <cstdint>
#include <string>
#include <pdal/Filter.hpp>
#include <pdal/pdal_export.hpp>
#include <pdal/util/ProgramArgs.hpp>
namespace pdal
{
class PDAL_DLL Align3dFilter : public Filter
{
public:
Align3dFilter() : Filter(), m_fixed(nullptr)
{
}
static void* create();
static int32_t destroy(void*);
std::string getName() const;
private:
virtual void addArgs(ProgramArgs& args);
virtual PointViewSet run(PointViewPtr view);
double m_gsd;
double m_maxt;
double m_maxdz;
PointViewPtr m_fixed;
Align3dFilter& operator=(const Align3dFilter&) = delete;
Align3dFilter(const Align3dFilter&) = delete;
};
}

22
src/common/PointCloud.cpp
View File

@@ -56,6 +56,28 @@ namespace pubgeo {
}
}
bool PointCloud::Read(pdal::PointViewPtr view) {
numPoints = view->size();
if (numPoints < 1) {
std::cerr << "[PUBGEO::PointCloud::READ] No points found in file." << std::endl;
return false;
}
pdal::BOX3D box;
view->calculateBounds(box);
pdal::SpatialReference sr = view->spatialReference();
zone = sr.computeUTMZone(box);
// used later to return points
xOff = (int) floor(box.minx);
yOff = (int) floor(box.miny);
zOff = (int) floor(box.minz);
bounds = {box.minx, box.maxx, box.miny, box.maxy, box.minz, box.maxz};
return true;
}
bool PointCloud::TransformPointCloud(const char *inputFileName, const char *outputFileName,
float translateX = 0, float translateY = 0, float translateZ = 0) {
std::ostringstream pipeline;

1
src/common/PointCloud.h
View File

@@ -33,6 +33,7 @@ namespace pubgeo {
float translateX, float translateY, float translateZ);
bool Read(const char *fileName);
bool Read(pdal::PointViewPtr view);
MinMaxXYZ bounds;
int zone;

166
src/common/orthoimage.h
View File

@@ -258,7 +258,7 @@ namespace pubgeo {
GDALRasterBand *poBand = poDstDS->GetRasterBand(i);
poBand->SetNoDataValue(noData);
for (unsigned int j = 0; j < this->height; j++) {
for (int k = 0; k < this->width * this->bands; k++) {
for (unsigned int k = 0; k < this->width * this->bands; k++) {
if (this->data[j][k] == 0)
raster[k] = noData;
else
@@ -312,19 +312,19 @@ namespace pubgeo {
// Copy points into the ortho image.
if (mode == MIN_VALUE) {
for (int i = 0; i < pset.numPoints; i++) {
int x = int((pset.x(i) + pset.xOff - easting) / gsd + 0.5);
for (unsigned long i = 0; i < pset.numPoints; i++) {
unsigned int x = int((pset.x(i) + pset.xOff - easting) / gsd + 0.5);
if ((x < 0) || (x > this->width - 1)) continue;
int y = this->height - 1 - int((pset.y(i) + pset.yOff - northing) / gsd + 0.5);
unsigned int y = this->height - 1 - int((pset.y(i) + pset.yOff - northing) / gsd + 0.5);
if ((y < 0) || (y > this->height - 1)) continue;
TYPE z = TYPE((pset.z(i) + pset.zOff - this->offset) / this->scale);
if ((this->data[y][x] == 0) || (z < this->data[y][x])) this->data[y][x] = z;
}
} else if (mode == MAX_VALUE) {
for (int i = 0; i < pset.numPoints; i++) {
int x = int((pset.x(i) + pset.xOff - easting) / gsd + 0.5);
for (unsigned long i = 0; i < pset.numPoints; i++) {
unsigned int x = int((pset.x(i) + pset.xOff - easting) / gsd + 0.5);
if ((x < 0) || (x > this->width - 1)) continue;
int y = this->height - 1 - int((pset.y(i) + pset.yOff - northing) / gsd + 0.5);
unsigned int y = this->height - 1 - int((pset.y(i) + pset.yOff - northing) / gsd + 0.5);
if ((y < 0) || (y > this->height - 1)) continue;
TYPE z = TYPE((pset.z(i) + pset.zOff - this->offset) / this->scale);
if ((this->data[y][x] == 0) || (z > this->data[y][x])) this->data[y][x] = z;
@@ -333,6 +333,60 @@ namespace pubgeo {
return true;
}
// Read image from PDAL PointView.
bool readFromPointView(pdal::PointViewPtr view, float gsdMeters, MIN_MAX_TYPE mode = MIN_VALUE) {
// Read a PSET file (e.g., BPF or LAS).
PointCloud pset;
bool ok = pset.Read(view);
if (!ok) return false;
// Calculate scale and offset for conversion to TYPE.
float minVal = pset.bounds.zMin - 1; // Reserve zero for noData value
float maxVal = pset.bounds.zMax + 1;
float maxImageVal = (float) (pow(2.0, int(sizeof(TYPE) * 8)) - 1);
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);
// 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;
// Copy points into the ortho image.
if (mode == MIN_VALUE) {
for (unsigned long i = 0; i < pset.numPoints; i++) {
double dx = view->getFieldAs<double>(pdal::Dimension::Id::X, i);
double dy = view->getFieldAs<double>(pdal::Dimension::Id::Y, i);
double dz = view->getFieldAs<double>(pdal::Dimension::Id::Z, i);
unsigned int x = int((dx - easting) / gsd + 0.5);
if ((x < 0) || (x > this->width - 1)) continue;
unsigned int y = this->height - 1 - int((dy - northing) / gsd + 0.5);
if ((y < 0) || (y > this->height - 1)) continue;
TYPE z = TYPE((dz - this->offset) / this->scale);
if ((this->data[y][x] == 0) || (z < this->data[y][x])) this->data[y][x] = z;
}
} else if (mode == MAX_VALUE) {
for (unsigned long i = 0; i < pset.numPoints; i++) {
double dx = view->getFieldAs<double>(pdal::Dimension::Id::X, i);
double dy = view->getFieldAs<double>(pdal::Dimension::Id::Y, i);
double dz = view->getFieldAs<double>(pdal::Dimension::Id::Z, i);
unsigned int x = int((dx - easting) / gsd + 0.5);
if ((x < 0) || (x > this->width - 1)) continue;
unsigned int y = this->height - 1 - int((dy - northing) / gsd + 0.5);
if ((y < 0) || (y > this->height - 1)) continue;
TYPE z = TYPE((dz - this->offset) / this->scale);
if ((this->data[y][x] == 0) || (z > this->data[y][x])) this->data[y][x] = z;
}
}
return true;
}
// Count voids in an image.
// Note that voids are always labeled zero in this data structure.
long countVoids() {
@@ -358,27 +412,27 @@ namespace pubgeo {
// Create image pyramid.
std::vector<OrthoImage<TYPE> *> pyramid;
pyramid.push_back(this);
int level = 0;
unsigned int level = 0;
while ((count > 0) && (level < maxLevel)) {
// Create next level.
OrthoImage<TYPE> *newImagePtr = new OrthoImage<TYPE>;
int nextWidth = pyramid[level]->width / 2;
int nextHeight = pyramid[level]->height / 2;
unsigned int nextWidth = pyramid[level]->width / 2;
unsigned int nextHeight = pyramid[level]->height / 2;
newImagePtr->Allocate(nextWidth, nextHeight, 1);
// Fill in non-void values from level below building up the pyramid with a simple running average.
for (int j = 0; j < nextHeight; j++) {
for (int i = 0; i < nextWidth; i++) {
int j2 = MIN(MAX(0, j * 2 + 1), pyramid[level]->height - 1);
int i2 = MIN(MAX(0, i * 2 + 1), pyramid[level]->width - 1);
for (unsigned int j = 0; j < nextHeight; j++) {
for (unsigned int i = 0; i < nextWidth; i++) {
unsigned int j2 = MIN(MAX(0, j * 2 + 1), pyramid[level]->height - 1);
unsigned int i2 = MIN(MAX(0, i * 2 + 1), pyramid[level]->width - 1);
// Average neighboring pixels from below.
float z = 0;
int ct = 0;
std::vector<TYPE> neighbors;
for (int jj = MAX(0, j2 - 1); jj <= MIN(j2 + 1, pyramid[level]->height - 1); jj++) {
for (int ii = MAX(0, i2 - 1); ii <= MIN(i2 + 1, pyramid[level]->width - 1); ii++) {
for (unsigned int jj = MAX(0, j2 - 1); jj <= MIN(j2 + 1, pyramid[level]->height - 1); jj++) {
for (unsigned int ii = MAX(0, i2 - 1); ii <= MIN(i2 + 1, pyramid[level]->width - 1); ii++) {
if (pyramid[level]->data[jj][ii] != 0) {
z += pyramid[level]->data[jj][ii];
ct++;
@@ -399,12 +453,12 @@ namespace pubgeo {
// Void fill down the pyramid.
for (int k = level - 1; k >= 0; k--) {
for (int j = 0; j < pyramid[k]->height; j++) {
for (int i = 0; i < pyramid[k]->width; i++) {
for (unsigned int j = 0; j < pyramid[k]->height; j++) {
for (unsigned int i = 0; i < pyramid[k]->width; i++) {
// Fill this pixel if it is currently void.
if (pyramid[k]->data[j][i] == 0) {
int j2 = MIN(MAX(0, j / 2), pyramid[k + 1]->height - 1);
int i2 = MIN(MAX(0, i / 2), pyramid[k + 1]->width - 1);
unsigned int j2 = MIN(MAX(0, j / 2), pyramid[k + 1]->height - 1);
unsigned int i2 = MIN(MAX(0, i / 2), pyramid[k + 1]->width - 1);
if (noSmoothing) {
// Just use the closest pixel from above.
@@ -413,9 +467,9 @@ namespace pubgeo {
// Average neighboring pixels from above.
float z = 0;
int ct = 0;
for (int jj = MAX(0, j2 - 1);
for (unsigned int jj = MAX(0, j2 - 1);
jj <= MIN(j2 + 1, pyramid[k + 1]->height - 1); jj++) {
for (int ii = MAX(0, i2 - 1);
for (unsigned int ii = MAX(0, i2 - 1);
ii <= MIN(i2 + 1, pyramid[k + 1]->width - 1); ii++) {
z += pyramid[k + 1]->data[jj][ii];
ct++;
@@ -430,28 +484,28 @@ namespace pubgeo {
}
// Deallocate memory for all but the input DSM.
for (int i = 1; i <= level; i++) {
for (unsigned int i = 1; i <= level; i++) {
pyramid[i]->Deallocate();
}
}
// Apply a median filter to an image.
void medianFilter(int rad, unsigned int dzShort) {
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
// Skip if void.
if (this->data[j][i] == 0) continue;
// Define bounds;
int i1 = MAX(0, i - rad);
int i2 = MIN(i + rad, this->width - 1);
int j1 = MAX(0, j - rad);
int j2 = MIN(j + rad, this->height - 1);
unsigned int i1 = MAX(0, i - rad);
unsigned int i2 = MIN(i + rad, this->width - 1);
unsigned int j1 = MAX(0, j - rad);
unsigned int j2 = MIN(j + rad, this->height - 1);
// Add valid values to the list.
std::vector<unsigned short> values;
for (int jj = j1; jj <= j2; jj++) {
for (int ii = i1; ii <= i2; ii++) {
for (unsigned int jj = j1; jj <= j2; jj++) {
for (unsigned int ii = i1; ii <= i2; ii++) {
if (this->data[jj][ii] != 0) {
values.push_back(this->data[jj][ii]);
}
@@ -474,25 +528,25 @@ namespace pubgeo {
void minFilter(int rad) {
OrthoImage<TYPE> tempImage;
tempImage.Allocate(this->width, this->height);
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
tempImage.data[j][i] = this->data[j][i];
}
}
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
// Skip if void.
if (this->data[j][i] == 0) continue;
// Define bounds;
int i1 = MAX(0, i - rad);
int i2 = MIN(i + rad, this->width - 1);
int j1 = MAX(0, j - rad);
int j2 = MIN(j + rad, this->height - 1);
unsigned int i1 = MAX(0, i - rad);
unsigned int i2 = MIN(i + rad, this->width - 1);
unsigned int j1 = MAX(0, j - rad);
unsigned int j2 = MIN(j + rad, this->height - 1);
// Check all neighbors.
for (int jj = j1; jj <= j2; jj++) {
for (int ii = i1; ii <= i2; ii++) {
for (unsigned int jj = j1; jj <= j2; jj++) {
for (unsigned int ii = i1; ii <= i2; ii++) {
if (this->data[jj][ii] != 0) {
if (this->data[jj][ii] < tempImage.data[j][i])
tempImage.data[j][i] = this->data[jj][ii];
@@ -501,8 +555,8 @@ namespace pubgeo {
}
}
}
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
this->data[j][i] = tempImage.data[j][i];
}
}
@@ -511,25 +565,25 @@ namespace pubgeo {
void edgeFilter(int dzShort) {
OrthoImage<TYPE> tempImage;
tempImage.Allocate(this->width, this->height);
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
tempImage.data[j][i] = this->data[j][i];
}
}
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
// Skip if void.
if (this->data[j][i] == 0) continue;
// Define bounds;
int i1 = MAX(0, i - 1);
int i2 = MIN(i + 1, this->width - 1);
int j1 = MAX(0, j - 1);
int j2 = MIN(j + 1, this->height - 1);
unsigned int i1 = MAX(0, i - 1);
unsigned int i2 = MIN(i + 1, this->width - 1);
unsigned int j1 = MAX(0, j - 1);
unsigned int j2 = MIN(j + 1, this->height - 1);
// If there's an edge with any neighbor, then remove.
for (int jj = j1; jj <= j2; jj++) {
for (int ii = i1; ii <= i2; ii++) {
for (unsigned int jj = j1; jj <= j2; jj++) {
for (unsigned int ii = i1; ii <= i2; ii++) {
if (fabs(float(this->data[jj][ii] - this->data[j][i])) > dzShort) {
tempImage.data[j][i] = 0;
}
@@ -537,8 +591,8 @@ namespace pubgeo {
}
}
}
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
this->data[j][i] = tempImage.data[j][i];
}
}
@@ -546,4 +600,4 @@ namespace pubgeo {
};
}
#endif //PUBGEO_ORTHO_IMAGE_H
#endif //PUBGEO_ORTHO_IMAGE_H

14
src/shr3d/CMakeLists.txt
View File

@@ -23,3 +23,17 @@ IF(TESTING)
PUBLIC
SHR3D_LIB)
ENDIF()
IF(WIN32)
SET(SHR3D_WRITER_NAME libpdal_plugin_writer_shr3d)
ELSE(WIN32)
SET(SHR3D_WRITER_NAME pdal_plugin_writer_shr3d)
ENDIF(WIN32)
ADD_LIBRARY(${SHR3D_WRITER_NAME} SHARED plugin.cpp)
LINK_DIRECTORIES(${PDAL_LIBRARY_DIRS})
TARGET_LINK_LIBRARIES(${SHR3D_WRITER_NAME} ${PDAL_LIBRARIES} SHR3D_LIB)
SET_TARGET_PROPERTIES(${SHR3D_WRITER_NAME} PROPERTIES SOVERSION "0.1.0")
INSTALL(TARGETS ${SHR3D_WRITER_NAME}
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)

35
src/shr3d/main.cpp
View File

@@ -120,8 +120,8 @@ int main(int argc, char **argv) {
minImage.write(minOutFileName, true);
#endif
// Find many of the trees by comparing MIN and MAX. Set their values to void.
for (int j = 0; j < dsmImage.height; j++) {
for (int i = 0; i < dsmImage.width; i++) {
for (unsigned int j = 0; j < dsmImage.height; j++) {
for (unsigned int i = 0; i < dsmImage.width; i++) {
float minValue = (float) minImage.data[j][i];
// This check is to avoid penalizing spurious returns under very tall buildings.
@@ -130,12 +130,12 @@ int main(int argc, char **argv) {
// If this is in the trees, then set to void.
bool found = false;
double threshold = dz_meters / dsmImage.scale;
int i1 = MAX(0, i - 1);
int i2 = MIN(i + 1, dsmImage.width - 1);
int j1 = MAX(0, j - 1);
int j2 = MIN(j + 1, dsmImage.height - 1);
for (int jj = j1; jj <= j2; jj++) {
for (int ii = i1; ii <= i2; ii++) {
unsigned int i1 = MAX(0, i - 1);
unsigned int i2 = MIN(i + 1, dsmImage.width - 1);
unsigned int j1 = MAX(0, j - 1);
unsigned int j2 = MIN(j + 1, dsmImage.height - 1);
for (unsigned int jj = j1; jj <= j2; jj++) {
for (unsigned int ii = i1; ii <= i2; ii++) {
float diff = (float) dsmImage.data[jj][ii] - minValue;
if (diff < threshold) found = true;
}
@@ -186,8 +186,8 @@ int main(int argc, char **argv) {
dtmImage.scale = dsmImage.scale;
dtmImage.offset = dsmImage.offset;
dtmImage.bands = dsmImage.bands;
for (int j = 0; j < dsmImage.height; j++) {
for (int i = 0; i < dsmImage.width; i++) {
for (unsigned int j = 0; j < dsmImage.height; j++) {
for (unsigned int i = 0; i < dsmImage.width; i++) {
dtmImage.data[j][i] = minImage.data[j][i];
}
}
@@ -197,8 +197,8 @@ int main(int argc, char **argv) {
// For DSM voids, also set DTM value to void.
printf("Setting DTM values to VOID where DSM is VOID...\n");
for (int j = 0; j < dsmImage.height; j++) {
for (int i = 0; i < dsmImage.width; i++) {
for (unsigned int j = 0; j < dsmImage.height; j++) {
for (unsigned int i = 0; i < dsmImage.width; i++) {
if (dsmImage.data[j][i] == 0) dtmImage.data[j][i] = 0;
}
}
@@ -224,8 +224,8 @@ int main(int argc, char **argv) {
classImage.northing = dsmImage.northing;
classImage.zone = dsmImage.zone;
classImage.gsd = dsmImage.gsd;
for (int j = 0; j < classImage.height; j++) {
for (int i = 0; i < classImage.width; i++) {
for (unsigned int j = 0; j < classImage.height; j++) {
for (unsigned int i = 0; i < classImage.width; i++) {
// Set default as unlabeled.
classImage.data[j][i] = LAS_UNCLASSIFIED;
@@ -250,8 +250,8 @@ int main(int argc, char **argv) {
char classOutFileName[1024];
sprintf(classOutFileName, "%s_class.tif\0", inputFileName);
classImage.write(classOutFileName, false, egm96);
for (int j = 0; j < classImage.height; j++) {
for (int i = 0; i < classImage.width; i++) {
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) classImage.data[j][i] = 0;
}
}
@@ -263,6 +263,3 @@ int main(int argc, char **argv) {
time(&t1);
printf("Total time elapsed = %f seconds\n", double(t1 - t0));
}

167
src/shr3d/plugin.cpp Normal file
View File

@@ -0,0 +1,167 @@
// Copyright (c) 2017, Bradley J Chambers (brad.chambers@gmail.com)
// All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// 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,”
// FOSS4G, 2017.
#include <algorithm>
#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");
CREATE_SHARED_PLUGIN(1, 0, Shr3dWriter, Writer, s_info)
std::string Shr3dWriter::getName() const
{
return s_info.name;
}
void Shr3dWriter::addArgs(ProgramArgs& args)
{
args.add("filename", "Output filename", m_filename).setPositional();
args.add("dh", "Horizontal uncertainty", m_dh, 0.5);
args.add("dz", "Vertical uncertainty", m_dz, 0.5);
args.add("agl", "Minimum building height above ground level", m_agl, 2.0);
args.add("area", "Minimum building area", m_area, 50.0);
args.add("egm96", "Set vertical datum to EGM96", m_egm96, false);
}
void Shr3dWriter::write(const PointViewPtr view)
{
shr3d::OrthoImage<unsigned short> dsmImage;
if (!dsmImage.readFromPointView(view, m_dh, shr3d::MAX_VALUE))
throw pdal_error("Error createing DSM\n");
dsmImage.medianFilter(1, static_cast<unsigned int>(m_agl / dsmImage.scale));
dsmImage.fillVoidsPyramid(true, 2);
shr3d::OrthoImage<unsigned short> minImage;
if (!minImage.readFromPointView(view, m_dh, shr3d::MIN_VALUE))
throw pdal_error("Error creating minimum Z image\n");
minImage.medianFilter(1, static_cast<unsigned int>(m_agl / minImage.scale));
minImage.fillVoidsPyramid(true, 2);
for (unsigned int j = 0; j < dsmImage.height; ++j)
{
for (unsigned int i = 0; i < dsmImage.width; ++i)
{
float minValue = static_cast<float>(minImage.data[j][i]);
// This check is to avoid penalizing spurious returns under very
// tall buildings. CAUTION: This is a hack to address an observed
// lidar sensor issue and may not generalize well.
if ((static_cast<float>(dsmImage.data[j][i]) - minValue) <
(40.0 / minImage.scale))
{
// If this is in the trees, then set to void.
bool found = false;
double threshold = m_dz / dsmImage.scale;
unsigned int i1 = std::max(0u, i - 1u);
unsigned int i2 = std::min(i + 1u, dsmImage.width - 1u);
unsigned int j1 = std::max(0u, j - 1u);
unsigned int j2 = std::min(j + 1u, dsmImage.height - 1u);
for (unsigned int jj = j1; jj <= j2; ++jj)
{
for (unsigned int ii = i1; ii <= i2; ++ii)
{
float diff = static_cast<float>(dsmImage.data[jj][ii]) -
minValue;
if (diff < threshold)
found = true;
}
}
if (!found)
dsmImage.data[j][i] = 0;
}
}
}
// Convert horizontal and vertical uncertainty values to bin units.
int dh_bins = std::max(1, static_cast<int>(floor(m_dh / dsmImage.gsd)));
unsigned int dz_short = static_cast<unsigned int>(m_dz / dsmImage.scale);
unsigned int agl_short = static_cast<unsigned int>(m_agl / dsmImage.scale);
// Generate label image.
shr3d::OrthoImage<unsigned long> labelImage;
labelImage.Allocate(dsmImage.width, dsmImage.height);
labelImage.easting = dsmImage.easting;
labelImage.northing = dsmImage.northing;
labelImage.zone = dsmImage.zone;
labelImage.gsd = dsmImage.gsd;
// Allocate a DTM image as SHORT and copy in the DSM values.
shr3d::OrthoImage<unsigned short> dtmImage;
dtmImage.Allocate(dsmImage.width, dsmImage.height);
dtmImage.easting = dsmImage.easting;
dtmImage.northing = dsmImage.northing;
dtmImage.zone = dsmImage.zone;
dtmImage.gsd = dsmImage.gsd;
dtmImage.scale = dsmImage.scale;
dtmImage.offset = dsmImage.offset;
dtmImage.bands = dsmImage.bands;
for (unsigned int j = 0; j < dsmImage.height; ++j)
{
for (unsigned int i = 0; i < dsmImage.width; ++i)
{
dtmImage.data[j][i] = minImage.data[j][i];
}
}
// Classify ground points.
shr3d::Shr3dder::classifyGround(labelImage, dsmImage, dtmImage, dh_bins,
dz_short);
// For DSM voids, also set DTM value to void.
for (unsigned int j = 0; j < dsmImage.height; ++j)
{
for (unsigned int i = 0; i < dsmImage.width; ++i)
{
if (dsmImage.data[j][i] == 0)
dtmImage.data[j][i] = 0;
}
}
// Median filter, replacing only points differing by more than the DZ
// threshold.
dtmImage.medianFilter(1, static_cast<unsigned int>(m_dz / dsmImage.scale));
// Refine the object label image and export building outlines.
shr3d::Shr3dder::classifyNonGround(dsmImage, dtmImage, labelImage, dz_short,
agl_short, static_cast<float>(m_area));
// Fill small voids in the DTM after all processing is complete.
dtmImage.fillVoidsPyramid(true, 2);
// Write the DTM image as FLOAT.
dtmImage.write(const_cast<char*>(m_filename.c_str()), true, m_egm96);
}
} // namespace pdal

63
src/shr3d/plugin.hpp Normal file
View File

@@ -0,0 +1,63 @@
// Copyright (c) 2017, Bradley J Chambers (brad.chambers@gmail.com)
// All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// 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,”
// FOSS4G, 2017.
#pragma once
#include <cstdint>
#include <string>
#include <pdal/Writer.hpp>
#include <pdal/pdal_export.hpp>
#include <pdal/util/ProgramArgs.hpp>
namespace pdal
{
class PDAL_DLL Shr3dWriter : public Writer
{
public:
Shr3dWriter()
{
}
static void* create();
static int32_t destroy(void*);
std::string getName() const;
private:
virtual void addArgs(ProgramArgs& args);
virtual void write(const PointViewPtr view);
std::string m_filename;
double m_dh;
double m_dz;
double m_agl;
double m_area;
bool m_egm96;
Shr3dWriter& operator=(const Shr3dWriter&) = delete;
Shr3dWriter(const Shr3dWriter&) = delete;
};
}

20
src/shr3d/shr3d.cpp
View File

@@ -250,12 +250,12 @@ bool addNeighbors(std::vector<PixelType> &neighbors, OrthoImage<unsigned long> &
OrthoImage<unsigned short> &dsmImage, ObjectType &obj, unsigned int dzShort) {
// Get neighbors for all pixels in the list.
std::vector<PixelType> newNeighbors;
for (long k = 0; k < neighbors.size(); k++) {
long i = neighbors[k].i;
long j = neighbors[k].j;
for (size_t k = 0; k < neighbors.size(); k++) {
unsigned int i = neighbors[k].i;
unsigned int j = neighbors[k].j;
float value = (float) dsmImage.data[j][i];
for (long jj = MAX(0, j - 1); jj <= MIN(j + 1, labelImage.height - 1); jj++) {
for (long ii = MAX(0, i - 1); ii <= MIN(i + 1, labelImage.width - 1); ii++) {
for (unsigned int jj = MAX(0, j - 1); jj <= MIN(j + 1, labelImage.height - 1); jj++) {
for (unsigned int ii = MAX(0, i - 1); ii <= MIN(i + 1, labelImage.width - 1); ii++) {
// Skip if pixel is already labeled or if it is LABEL_GROUND.
// Note that non-ground labels are initialized with 1.
if (labelImage.data[jj][ii] > 1) continue;
@@ -646,11 +646,11 @@ addClassNeighbors(std::vector<PixelType> &neighbors, OrthoImage<unsigned char> &
OrthoImage<unsigned char> &labeled, unsigned int label) {
// Get neighbors for all pixels in the list.
std::vector<PixelType> newNeighbors;
for (long k = 0; k < neighbors.size(); k++) {
long i = neighbors[k].i;
long j = neighbors[k].j;
for (long jj = MAX(0, j - 1); jj <= MIN(j + 1, classImage.height - 1); jj++) {
for (long ii = MAX(0, i - 1); ii <= MIN(i + 1, classImage.width - 1); ii++) {
for (size_t k = 0; k < neighbors.size(); k++) {
unsigned int i = neighbors[k].i;
unsigned int j = neighbors[k].j;
for (unsigned int jj = MAX(0, j - 1); jj <= MIN(j + 1, classImage.height - 1); jj++) {
for (unsigned int ii = MAX(0, i - 1); ii <= MIN(i + 1, classImage.width - 1); ii++) {
// If already labeled, then skip.
if (labeled.data[jj][ii] == 1) continue;