AtHeartEngineer/pubgeo
1
0
Fork 0
mirror of https://github.com/JHUAPL/pubgeo.git synced 2026-01-08 22:18:11 -05:00
Code Issues Packages Projects Releases Wiki Activity

Pubgeo first release. Enjoy.

Browse Source
This commit is contained in:
Scott Almes
2017-08-16 21:23:28 -04:00
parent f8251901e6
commit ef6faab171
22 changed files with 2475 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
CMakeLists.txt Normal file
View File

@@ -0,0 +1,28 @@
CMAKE_MINIMUM_REQUIRED(VERSION 3.0.2)
PROJECT(pubgeo)
SET(CMAKE_CXX_STANDARD 11)
IF(WIN32)
MESSAGE(ERROR "Windows builds are currently not working due to erroneous linking errors.")
# 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)
ENDIF()
FIND_PACKAGE(GDAL 1.9.0 REQUIRED)# Using Open Source Geo for Windows installer
FIND_PACKAGE(PDAL 1.4.0 REQUIRED)
# This makes linking/including a little cleaner looking later
SET(DAL_LIBS ${GDAL_LIBRARY} ${PDAL_LIBRARIES})
SET(DAL_INCLUDE_DIRS ${GDAL_INCLUDE_DIR} ${PDAL_INCLUDE_DIRS})
MESSAGE(STATUS "${DAL_INCLUDE_DIRS}")
SET(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
ADD_SUBDIRECTORY(src/common)
ADD_SUBDIRECTORY(src/align3d)
ADD_SUBDIRECTORY(src/shr3d)

18
DISTRIBUTION Normal file
View File

@@ -0,0 +1,18 @@
Distribution statement:
Approved for Public Release; Distribution is Unlimited
Please reference the following when reporting any results using this software:
M. Bosch, A. Leichtman, D. Chilcott, H. Goldberg, M. Brown, “Metric Evaluation Pipeline for 3D Modeling of Urban Scenes,” ISPRS Archives, 2017.
S. Almes, S. Hagstrom, D. Chilcott, H. Goldberg, M. Brown, “Open Source Geospatial Tools to Enable Large Scale 3D Scene Modeling,” FOSS4G, 2017.
For more information, please see: http://www.jhuapl.edu/pubgeo.html

32
Dockerfile Normal file
View File

@@ -0,0 +1,32 @@
FROM ubuntu:17.04
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 \
&& rm -rf /var/lib/apt/lists/*
RUN git clone https://github.com/pubgeo/pubgeo
# Make a directory to work out of, and change to it
WORKDIR /build
RUN cmake ../pubgeo && make -j 10
# cleanup
RUN rm -rf /pubgeo
RUN apt purge -y \
build-essential \
libgdal-dev \
libpdal-dev \
cmake \
git
CMD echo "Please run a valid executable:" && \
echo "docker run -v <path to 3D data> shr3d <3D file> DH=2 DZ=1 AGL=2 AREA=50" && \
echo "docker run -v <path to point cloud> align3d <reference point cloud> <target pc> gsd=1.0 maxt=10.0"

22
LICENSE
View File

@@ -1,21 +1,7 @@
MIT License
Copyright 2016 The Johns Hopkins University Applied Physics Laboratory
Copyright (c) 2017 pubgeo
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:
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 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.
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.

2
README Normal file
View File

@@ -0,0 +1,2 @@
To be updated 8/17/2017

28
src/align3d/CMakeLists.txt Normal file
View File

@@ -0,0 +1,28 @@
SET(ALIGN3D_HEADER_FILES
align3d.h)
SET(ALIGN3D_SOURCE_FILES
align3d.cpp)
ADD_LIBRARY(ALIGN3D_LIB STATIC ${ALIGN3D_HEADER_FILES} ${ALIGN3D_SOURCE_FILES})
TARGET_INCLUDE_DIRECTORIES(ALIGN3D_LIB
PUBLIC
PUBGEO_LIB)
TARGET_LINK_LIBRARIES(ALIGN3D_LIB
PUBLIC
PUBGEO_LIB)
### Add some executables to play with.
### Main exe will showcase all features
ADD_EXECUTABLE(align3d main.cpp)
TARGET_LINK_LIBRARIES(align3d
PUBLIC
ALIGN3D_LIB)
if(TESTING)
ADD_EXECUTABLE(align3dTest test-main.cpp)
TARGET_LINK_LIBRARIES(align3dTest
PUBLIC
ALIGN3D_LIB)
endif()

261
src/align3d/align3d.cpp Normal file
View File

@@ -0,0 +1,261 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
#include <stdlib.h>
#include <iostream>
#include <algorithm>
#include "align3d.h"
namespace align3d {
bool computeRMS(float dx, float dy, long numSamples, long maxSamples, std::vector<double> &xlist,
std::vector<double> &ylist, OrthoImage<unsigned short> &referenceDSM,
OrthoImage<unsigned short> &targetDSM, float &medianDZ, float &rms, long &ndx,
float &completeness) {
// Loop on number of samples and compute RMS.
// Just in case there aren't many valid points, don't let this loop forever.
long count = 0;
ndx = 0;
std::vector<float> differences;
while ((count < numSamples) && (ndx < maxSamples)) {
// Get the next point for matching.
double x = xlist[ndx];
double y = ylist[ndx];
ndx++;
// Map the point into the target.
// Skip if this isn't a valid point.
long col = long((x - targetDSM.easting + dx) / targetDSM.gsd + 0.5);
long row = targetDSM.height - 1 - long((y - targetDSM.northing + dy) / targetDSM.gsd + 0.5);
if (col <= 0) continue;
if (row <= 0) continue;
if (col >= targetDSM.width - 1) continue;
if (row >= targetDSM.height - 1) continue;
if (targetDSM.data[row][col] == 0) continue;
float targetZ = float(targetDSM.data[row][col]) * targetDSM.scale + targetDSM.offset;
// Map the point into the reference.
// Skip if this isn't a valid point.
col = long((x - referenceDSM.easting) / referenceDSM.gsd + 0.5);
row = referenceDSM.height - 1 - long((y - referenceDSM.northing) / referenceDSM.gsd + 0.5);
if (col <= 0) continue;
if (row <= 0) continue;
if (col >= referenceDSM.width - 1) continue;
if (row >= referenceDSM.height - 1) continue;
if (referenceDSM.data[row][col] == 0) continue;
float referenceZ = float(referenceDSM.data[row][col]) * referenceDSM.scale + referenceDSM.offset;
// Keep going until we have enough points.
float difference = referenceZ - targetZ;
differences.push_back(difference);
count++;
}
// Skip if not enough sampled points.
if (count < numSamples) return false;
// Compute median Z offset and a robust estimate of the RMS difference.
rms = 0.0;
sort(differences.begin(), differences.end());
medianDZ = differences[count / 2];
for (long k = 0; k < count; k++) {
differences[k] = fabs(differences[k] - medianDZ);
}
sort(differences.begin(), differences.end());
rms = differences[(long) (count * 0.67)];
// Compute the completeness.
long good = 0;
for (long k = 0; k < count; ++k) {
if (differences[k] < 1.0) ++good;
}
completeness = good / (float) numSamples;
return true;
}
// Estimate 3D rigid body transform parameters to align target points with reference.
void EstimateRigidBody(OrthoImage<unsigned short> &referenceDSM, OrthoImage<unsigned short> &targetDSM, float maxt,
AlignBounds &bounds, AlignResult &result) {
float step = MIN(referenceDSM.gsd, targetDSM.gsd);
long numSamples = 10000;
long maxSamples = numSamples * 10;
// Allocate RMS array.
long bins = long(maxt / step * 2) + 1;
float **rmsArray = new float *[bins];
for (long i = 0; i < bins; i++) rmsArray[i] = new float[bins];
// Get random samples.
srand(0);
long seed = 0;
long count = 0;
std::vector<double> xlist;
std::vector<double> ylist;
for (long i = 0; i < maxSamples; i++) {
// Get a random point in the overlap area.
double x = bounds.xmin + ((rand() / (double) RAND_MAX) * bounds.width);
double y = bounds.ymin + ((rand() / (double) RAND_MAX) * bounds.height);
xlist.push_back(x);
ylist.push_back(y);
}
// Start with brute force, but sample points to reduce timeline.
// Loop on X and Y translation within search distance.
float threshold = MAX_FLOAT;
float bestDX = 0.0;
float bestDY = 0.0;
float bestDZ = 0.0;
long besti = 0;
long bestj = 0;
float bestRMS = MAX_FLOAT;
float medianDZ = 0.0;
bestRMS = threshold;
float bestCompleteness = 0.0;
long numSampled = 0;
for (long i = 0; i < bins; i++) {
float dx = -maxt + i * step;
for (long j = 0; j < bins; j++) {
float dy = -maxt + j * step;
rmsArray[i][j] = 0.0;
float rms = 0.0;
float completeness = 0.0;
bool ok = computeRMS(dx, dy, numSamples, maxSamples, xlist, ylist, referenceDSM, targetDSM, medianDZ,
rms, numSampled, completeness);
if (!ok) continue;
rmsArray[i][j] = rms;
if (rms < bestRMS) {
bestCompleteness = completeness;
bestRMS = rms;
bestDX = dx;
bestDY = dy;
bestDZ = medianDZ;
besti = i;
bestj = j;
}
}
}
// Apply quadratic interpolation to localize the peak.
if ((besti > 0) && (besti < bins - 1) && (bestj > 0) && (bestj < bins - 1)) {
float dx = (rmsArray[besti + 1][bestj] - rmsArray[besti - 1][bestj]) / 2.f;
float dy = (rmsArray[besti][bestj + 1] - rmsArray[besti][bestj - 1]) / 2.f;
float dxx = (rmsArray[besti + 1][bestj] + rmsArray[besti - 1][bestj] - 2 * rmsArray[besti][bestj]);
float dyy = (rmsArray[besti][bestj + 1] + rmsArray[besti][bestj - 1] - 2 * rmsArray[besti][bestj]);
float dxy =
(rmsArray[besti + 1][bestj + 1] - rmsArray[besti + 1][bestj - 1] - rmsArray[besti - 1][bestj + 1] +
rmsArray[besti - 1][bestj - 1]) / 4.f;
float det = dxx * dyy - dxy * dxy;
if (det != 0.0) {
float ix = besti - (dyy * dx - dxy * dy) / det;
float iy = bestj - (dxx * dy - dxy * dx) / det;
bestDX = -maxt + ix * step;
bestDY = -maxt + iy * step;
}
}
// Deallocate RMS array.
for (long i = 0; i < bins; i++) delete[]rmsArray[i];
delete[]rmsArray;
// Update the result and return.
result.rms = bestRMS;
result.tx = -bestDX;
result.ty = -bestDY;
result.tz = bestDZ;
printf("Percent less than 1m Z difference = %6.2f%%\n", bestCompleteness * 100.0);
printf("X offset = %f m\n", result.tx);
printf("Y offset = %f m\n", result.ty);
printf("Z offset = %f m\n", result.tz);
printf("Z RMS = %f m\n", result.rms);
}
// Align target file to match reference file.
bool AlignTarget2Reference(char *referenceFileName, char *targetFileName, AlignParameters params) {
// Read the reference LAS file as a DSM.
// Fill small voids.
// Remove points along edges which are difficult to match.
printf("Reading reference point cloud: %s\n", referenceFileName);
OrthoImage<unsigned short> referenceDSM;
bool ok = referenceDSM.readFromPointCloud(referenceFileName, params.gsd, MAX_VALUE);
if (!ok) {
printf("Failed to read %s\n", referenceFileName);
return false;
}
referenceDSM.fillVoidsPyramid(true, 2);
printf("Filtering reference point cloud.\n");
referenceDSM.edgeFilter((long) (params.maxdz / referenceDSM.scale));
// Read the target LAS file as a DSM.
// Fill small voids.
// Remove points along edges which are difficult to match.
printf("Reading target point cloud: %s\n", targetFileName);
OrthoImage<unsigned short> targetDSM;
ok = targetDSM.readFromPointCloud(targetFileName, params.gsd, MAX_VALUE);
if (!ok) {
printf("Failed to read %s\n", targetFileName);
return false;
}
targetDSM.fillVoidsPyramid(true, 2);
printf("Filtering target point cloud.\n");
targetDSM.edgeFilter((long) (params.maxdz / targetDSM.scale));
// Get overlapping bounds.
AlignBounds bounds;
bounds.xmin = MAX(referenceDSM.easting, targetDSM.easting);
bounds.ymin = MAX(referenceDSM.northing, targetDSM.northing);
bounds.xmax = MIN(referenceDSM.easting + (referenceDSM.width * referenceDSM.gsd),
targetDSM.easting + (targetDSM.width * targetDSM.gsd));
bounds.ymax = 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;
printf("Overlap = %d m x %d m = %f km\n", (long) bounds.width, (long) bounds.height, overlap_km);
if (overlap_km == 0.0) return false;
// Estimate rigid body transform to align target points to reference.
AlignResult result;
printf("Estimating rigid body transformation.\n");
EstimateRigidBody(referenceDSM, targetDSM, params.maxt, bounds, result);
// Write offsets text file.
printf("Writing offsets text file.\n");
int len = strlen(targetFileName);
char outFileName[1024];
sprintf(outFileName, "%s", targetFileName);
sprintf(&outFileName[len - 4], "_offsets.txt");
FILE *fptr = fopen(outFileName, "w");
if (fptr) {
fprintf(fptr, "X Offset Y Offset Z Offset Z RMS\n");
fprintf(fptr, "%08.3f %08.3f %08.3f %08.3f\n", result.tx, result.ty, result.tz, result.rms);
fclose(fptr);
} else {
printf("Failed to write %s\n", outFileName);
return false;
}
// Write aligned TIF file.
printf("Writing aligned TIF file.\n");
sprintf(&outFileName[len - 4], "_aligned.tif");
targetDSM.offset += result.tz;
targetDSM.easting += result.tx;
targetDSM.northing += result.ty;
ok = targetDSM.write(outFileName, true);
if (!ok) {
printf("Failed to write %s\n", outFileName);
return false;
}
// Write aligned BPF file.
// For now, this requires an extra read of the point cloud file.
printf("Writing aligned LAS file.\n");
sprintf(&outFileName[len - 4], "_aligned.las");
ok = PointCloud::TransformPointCloud(targetFileName, outFileName, result.tx, result.ty, result.tz);
if (!ok) {
printf("Failed to write %s\n", outFileName);
return false;
}
}
}

51
src/align3d/align3d.h Normal file
View File

@@ -0,0 +1,51 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
#ifndef PUBGEO_ALIGN3D_H
#define PUBGEO_ALIGN3D_H
#include <vector>
#include "orthoimage.h"
#include "Image.h"
namespace align3d {
using namespace pubgeo;
typedef struct {
float tx;
float ty;
float tz;
float rms;
}
AlignResult;
typedef struct {
float gsd;
float maxdz;
float maxt;
}
AlignParameters;
typedef struct {
double xmin;
double xmax;
double ymin;
double ymax;
double width;
double height;
}
AlignBounds;
bool computeRMS(float dx, float dy, long numSamples, long maxSamples, std::vector<double> &xlist,
std::vector<double> &ylist, OrthoImage<unsigned short> &referenceDSM,
OrthoImage<unsigned short> &targetDSM, float &medianDZ, float &rms, long &ndx, float &completeness);
// Estimate 3D rigid body transform parameters to align target points with reference.
void EstimateRigidBody(OrthoImage<unsigned short> &referenceDSM, OrthoImage<unsigned short> &targetDSM, float maxt,
AlignBounds &bounds, AlignResult &result);
// Align target file to match reference file.
bool AlignTarget2Reference(char *referenceFileName, char *targetFileName, AlignParameters params);
}
#endif // PUBGEO_ALIGN3D_H

72
src/align3d/main.cpp Normal file
View File

@@ -0,0 +1,72 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
#include "align3d.h"
void printArguments();
// Main program for simple 3d alignment.
int main(int argc, char **argv) {
// If no parameters, then print command line arguments.
if (argc < 4) {
printf("\nNumber of arguments = %d\n", argc);
for (int i = 1; i < argc; i++) {
printf(" ARG[%d] = %s\n", i, argv[i]);
}
printArguments();
return -1;
}
// Parse command line parameters.
char referenceFileName[1024];
char targetFileName[1024];
sprintf(referenceFileName, argv[1]);
sprintf(targetFileName, argv[2]);
align3d::AlignParameters params;
params.gsd = 1.0;
params.maxt = 10.0;
params.maxdz = 0.0;
for (int i = 2; i < argc; i++) {
if (strstr(argv[i], "maxdz=")) { params.maxdz = (float) atof(&(argv[i][6])); }
if (strstr(argv[i], "gsd=")) { params.gsd = (float) atof(&(argv[i][4])); }
if (strstr(argv[i], "maxt=")) { params.maxt = (float) atof(&(argv[i][5])); }
}
// Default MAXDZ = GSD x 2 to ensure reliable performance on steep slopes.
if (params.maxdz == 0.0) params.maxdz = params.gsd * 2.0;
printf("Selected Parameters:\n");
printf(" ref = %s\n", referenceFileName);
printf(" tgt = %s\n", targetFileName);
printf(" gsd = %f\n", params.gsd);
printf(" maxdz = %f\n", params.maxdz);
printf(" maxt = %f\n", params.maxt);
// Initialize the timer.
time_t t0;
time(&t0);
try {
// Align the target point cloud to the reference.
AlignTarget2Reference(referenceFileName, targetFileName, params);
} catch (char const *err) {
std::cerr << "Reporting error: " << err << std::endl;
}
// Report total elapsed time.
time_t t1;
time(&t1);
printf("Total time elapsed = %f seconds\n", double(t1 - t0));
return 0;
}
// Print command line arguments.
void printArguments() {
printf("Command Line: align-3d <reference> <target> <parameters>\n");
printf("Parameters:\n");
printf(" maxdz= Max local Z difference (meters) for matching\n");
printf(" gsd= Ground Sample Distance (GSD) for gridding (meters)\n");
printf(" maxt= Maximum XYZ translation in search (meters); default = 10.0\n");
printf("Examples:\n");
printf(" align-3d ref.las tgt.las maxt=10.0 gsd=0.5 maxdz=0.5 \n\n");
}

14
src/align3d/test-main.cpp Normal file
View File

@@ -0,0 +1,14 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
#include <PointCloud.h>
using namespace pubgeo;
int main() {
if (!PointCloud::TransformPointCloud("someFile.las", "someOut.las", 3.0, 4.0, 5.0)) {
return 1;
}
return 0;
}

18
src/common/CMakeLists.txt Normal file
View File

@@ -0,0 +1,18 @@
SET(PUBGEO_HEADER_FILES
util.h
Image.h
orthoimage.h
PointCloud.h)
SET(PUBGEO_SOURCE_FILES
PointCloud.cpp
)
ADD_LIBRARY(PUBGEO_LIB STATIC ${PUBGEO_HEADER_FILES} ${PUBGEO_SOURCE_FILES})
TARGET_INCLUDE_DIRECTORIES(PUBGEO_LIB
PUBLIC
${DAL_INCLUDE_DIRS}
${CMAKE_CURRENT_SOURCE_DIR})
TARGET_LINK_LIBRARIES(PUBGEO_LIB
PUBLIC
${DAL_LIBS})

46
src/common/Image.h Normal file
View File

@@ -0,0 +1,46 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
// Created by almess1 on 4/21/17.
//
#ifndef PUBGEO_IMAGE_H
#define PUBGEO_IMAGE_H
#include <cstring>
namespace pubgeo {
template<class TYPE>
class Image {
public:
unsigned int width;
unsigned int height;
unsigned int bands;
float scale;
float offset;
TYPE **data; // Order is HEIGHT (outside), WIDTH, BAND (inside).
// Deallocate memory.
void Deallocate() {
if (!data) return;
for (unsigned int y = 0; y < height; y++) delete[]data[y];
delete[]data;
data = nullptr;
}
// Allocate memory.
void Allocate(unsigned int numColumns, unsigned int numRows, unsigned int numBands = 1) {
if (data) Deallocate();
width = numColumns;
height = numRows;
bands = numBands;
data = new TYPE *[height];
for (unsigned int y = 0; y < height; y++) {
data[y] = new TYPE[width * bands];
std::memset(data[y], 0, width * bands * sizeof(TYPE));
}
}
};
}
#endif //PUBGEO_IMAGE_H

91
src/common/PointCloud.cpp Normal file
View File

@@ -0,0 +1,91 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
#include "PointCloud.h"
// Pipeline needs to read in point cloud file of any type, and read it in. ideally in meters
static std::string PDAL_PIPELINE_OPEN_ENGINE = R"({ "pipeline": [ ")";
static std::string PDAL_PIPELINE_OPEN_CABOOSE = R"("] } )";
std::string buildPipelineStr(const char *fileName) {
return PDAL_PIPELINE_OPEN_ENGINE + fileName + PDAL_PIPELINE_OPEN_CABOOSE;
}
namespace pubgeo {
PointCloud::PointCloud() : executor(nullptr), pv(nullptr), zone(0), xOff(0), yOff(0), zOff(0), numPoints(0) {
bounds = MinMaxXYZ{0, 0, 0, 0, 0, 0};
}
PointCloud::~PointCloud() {
CleanupPdalPointers();
}
bool PointCloud::Read(const char *fileName) {
try {
CleanupPdalPointers();
executor = new pdal::PipelineExecutor(buildPipelineStr(fileName));
executor->execute();
const pdal::PointViewSet &pvs = executor->getManagerConst().views();
pv = pvs.begin()->get();
if (pvs.size() > 1) {
std::cerr << "[PUBGEO::PointCloud::READ] File contains additional unread sets." << std::endl;
}
numPoints = pv->size();
if (numPoints < 1) {
std::cerr << "[PUBGEO::PointCloud::READ] No points found in file." << std::endl;
return false;
}
pdal::BOX3D box;
pv->calculateBounds(box);
pdal::SpatialReference sr = pv->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;
}
catch (pdal::pdal_error &pe) {
std::cerr << pe.what() << std::endl;
return false;
}
}
bool PointCloud::TransformPointCloud(const char *inputFileName, const char *outputFileName,
float translateX = 0, float translateY = 0, float translateZ = 0) {
std::ostringstream pipeline;
pipeline << "{\n\t\"pipeline\":[\n\t\t\"" << inputFileName
<< "\",\n\t\t{\n\t\t\t\"type\":\"filters.transformation\",\n"
<< "\t\t\t\"matrix\":\""
<< "1 0 0 " << translateX << " "
<< "0 1 0 " << translateY << " "
<< "0 0 1 " << translateZ << " "
<< "0 0 0 1\"\n\t\t},\n\t\t{"
<< "\n\t\t\t\"filename\":\"" << outputFileName << "\"\n\t\t}\n\t]\n}";
try {
const std::string pipe = pipeline.str();
pdal::PipelineExecutor executor(pipe);
executor.execute();
} catch (pdal::pdal_error &pe) {
std::cerr << pe.what() << std::endl;
return false;
}
return true;
}
void PointCloud::CleanupPdalPointers() {
if (pv != nullptr) {
pv = nullptr;
}
if (executor != nullptr) {
delete executor;
executor = nullptr;
}
}
};

69
src/common/PointCloud.h Normal file
View File

@@ -0,0 +1,69 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
//
// Created by almess1 on 4/20/17.
// This file was created in order to adapt the SHR3D and ALIGN3D to open source IO.
//
#ifndef PUBGEO_NOT_POINT_SETS_H
#define PUBGEO_NOT_POINT_SETS_H
#include <pdal/PointView.hpp>
#include <pdal/PipelineExecutor.hpp>
namespace pubgeo {
struct MinMaxXYZ {
double xMin;
double xMax;
double yMin;
double yMax;
double zMin;
double zMax;
};
class PointCloud {
public:
PointCloud();
~PointCloud();
static bool TransformPointCloud(const char *inputFileName, const char *outputFileName,
float translateX, float translateY, float translateZ);
bool Read(const char *fileName);
MinMaxXYZ bounds;
int zone;
unsigned long numPoints;
int xOff;
int yOff;
int zOff;
inline float x(int i) {
if (pv != nullptr)
return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::X, i) - xOff);
throw "Point set has not be initialized.";
}
inline float y(int i) {
if (pv != nullptr)
return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::Y, i) - yOff);
throw "Point set has not be initialized.";
}
inline float z(int i) {
if (pv != nullptr)
return (float) (pv->getFieldAs<double>(pdal::Dimension::Id::Z, i) - zOff);
throw "Point set has not be initialized.";
}
private:
pdal::PipelineExecutor *executor;
pdal::PointView *pv;
void CleanupPdalPointers();
};
}
#endif //PUBGEO_NOT_POINT_SETS_H

549
src/common/orthoimage.h Normal file
View File

@@ -0,0 +1,549 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
// orthoimage.h
//
#ifndef PUBGEO_ORTHO_IMAGE_H
#define PUBGEO_ORTHO_IMAGE_H
#include <stdio.h>
#include <math.h>
#include <vector>
#include <typeinfo>
#include <cstring>
#include <algorithm>
#ifdef WIN32
#include "gdal_priv.h"
#include "ogr_spatialref.h"
#else
#include <gdal/gdal_priv.h>
#include "gdal/ogr_spatialref.h"
#endif
#include "util.h"
#include "PointCloud.h"
#include "Image.h"
namespace pubgeo {
typedef enum {
MIN_VALUE, MAX_VALUE
} MIN_MAX_TYPE;
//
// Ortho image template class
//
template<class TYPE>
class OrthoImage : public Image<TYPE> {
public:
double easting;
double northing;
int zone;
float gsd;
// Constructor
OrthoImage() {
memset(this, 0, sizeof(OrthoImage));
this->scale = 1.0f;
}
// Destructor
~OrthoImage() {
this->Deallocate();
}
// Read any GDAL-supported image.
bool read(char *fileName) {
// Open the image.
GDALAllRegister();
CPLSetConfigOption("GDAL_DATA", ".\\gdaldata");
GDALDataset *poDataset = (GDALDataset *) GDALOpen(fileName, GA_ReadOnly);
if (poDataset == NULL) {
printf("Error opening %s.\n", fileName);
return false;
}
// Get geospatial metadata.
double adfGeoTransform[6];
printf("Driver: %s/%s\n", poDataset->GetDriver()->GetDescription(),
poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));
this->width = (unsigned int) poDataset->GetRasterXSize();
this->height = (unsigned int) poDataset->GetRasterYSize();
this->bands = (unsigned int) poDataset->GetRasterCount();
printf("Width = %d\nHeight = %d\nBands = %d\n", this->width, this->height, this->bands);
if (poDataset->GetProjectionRef() != NULL) printf("Projection is %s.\n", poDataset->GetProjectionRef());
if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
printf("GeoTransform = (%.6f,%.6f,%.6f,%.6f,%.6f,%.6f)\n",
adfGeoTransform[0], adfGeoTransform[1], adfGeoTransform[2], adfGeoTransform[3],
adfGeoTransform[4],
adfGeoTransform[5]);
}
double xscale = adfGeoTransform[1];
double yscale = -adfGeoTransform[5];
this->gsd = (float) ((xscale + yscale) / 2.0);
this->easting = adfGeoTransform[0] - adfGeoTransform[2] * xscale;
this->northing = adfGeoTransform[3] - (this->height) * yscale;
OGRSpatialReference myOGRS(poDataset->GetProjectionRef());
int north = 1;
this->zone = myOGRS.GetUTMZone(&north);
if (north == 0) {
this->zone *= -1;
}
printf("UTM Easting = %lf\n", this->easting);
printf("UTM Northing = %lf\n", this->northing);
printf("UTM Zone = %d\n", this->zone);
printf("GSD = %f\n", this->gsd);
// Get band information.
GDALRasterBand *poBand = poDataset->GetRasterBand(1);
if (poBand == NULL) {
printf("Error opening first band.");
return false;
}
GDALDataType BandDataType = poBand->GetRasterDataType();
int bytesPerPixel = GDALGetDataTypeSize(BandDataType) / 8;
unsigned char *poBandBlock = (unsigned char *) CPLMalloc(bytesPerPixel * this->width);
int ok = 0;
double noData = poBand->GetNoDataValue(&ok);
if (!ok) {
// Set noData only for floating point images.
if (strcmp(typeid(TYPE).name(), "float") == 0)
noData = -10000.0;
else
noData = 0;
}
// Get scale and offset values.
if (strcmp(typeid(TYPE).name(), "float") == 0) {
// Do not scale if floating point values.
this->scale = 1.0;
this->offset = 0.0;
} else {
float minVal = MAX_FLOAT;
float maxVal = MAX_FLOAT;
for (unsigned int i = 0; i < this->bands; i++) {
poBand = poDataset->GetRasterBand(i + 1);
if (poBand == NULL) {
printf("Error opening band %d", i + 1);
return false;
}
int bGotMin, bGotMax;
double adfMinMax[2];
adfMinMax[0] = poBand->GetMinimum(&bGotMin);
adfMinMax[1] = poBand->GetMaximum(&bGotMax);
if (!(bGotMin && bGotMax)) GDALComputeRasterMinMax((GDALRasterBandH) poBand, false, adfMinMax);
if (minVal > adfMinMax[0]) minVal = (float) adfMinMax[0];
if (maxVal < adfMinMax[1]) maxVal = (float) adfMinMax[1];
}
minVal -= 1; // Reserve zero for noData value
maxVal += 1;
float maxImageVal = (float) (pow(2.0, int(sizeof(TYPE) * 8)) - 1);
this->offset = minVal;
this->scale = (maxVal - minVal) / maxImageVal;
}
printf("Offset = %f\n", this->offset);
printf("Scale = %f\n", this->scale);
// Allocate memory for the image.
this->Allocate(this->width, this->height, this->bands);
// Read the image, one band at a time.
for (unsigned int irow = 0; irow < this->height; irow++) {
for (unsigned int ib = 0; ib < this->bands; ib++) {
// Read the next row.
poBand = poDataset->GetRasterBand(ib + 1);
poBand->RasterIO(GF_Read, 0, irow, this->width, 1, poBandBlock, this->width, 1, BandDataType, 0, 0);
switch (BandDataType) {
case GDT_Byte: {
for (unsigned int icols = 0; icols < this->width; icols++) {
long k = icols * this->bands;
if (poBandBlock[icols] == (unsigned char) noData)
this->data[irow][k + ib] = 0;
else
this->data[irow][k + ib] = (TYPE) ((poBandBlock[icols] - this->offset) /
this->scale);
}
break;
}
case GDT_UInt16: {
unsigned int *poBandBlockShort = (unsigned int *) poBandBlock;
for (unsigned int icols = 0; icols < this->width; icols++) {
long k = icols * this->bands;
if (poBandBlockShort[icols] == (unsigned int) noData)
this->data[irow][k + ib] = 0;
else
this->data[irow][k + ib] = (TYPE) ((poBandBlockShort[icols] - this->offset) /
this->scale);
}
break;
}
case GDT_Float32: {
float *poBandBlockFloat = (float *) poBandBlock;
for (unsigned int icols = 0; icols < this->width; icols++) {
long k = icols * this->bands;
if (poBandBlockFloat[icols] == (float) noData)
this->data[irow][k + ib] = 0;
else
this->data[irow][k + ib] = (TYPE) ((poBandBlockFloat[icols] - this->offset) /
this->scale);
}
break;
}
default: {
// TODO: shouldn't get here, relay error
}
}
}
}
// Deallocate memory.
CPLFree(poBandBlock);
delete poDataset;
return true;
}
// Write GEOTIFF image using GDAL.
bool write(char *fileName, bool convertToFloat = false, bool egm96 = false) {
GDALAllRegister();
const char *pszFormat = "GTiff";
GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if (poDriver == NULL) return false;
char **papszMetadata = poDriver->GetMetadata();
if (!CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, false)) return false;
// Get the GDAL data type to match the image data type.
GDALDataType theBandDataType = (GDALDataType) sizeof(TYPE);
if (strcmp(typeid(TYPE).name(), "float") == 0) theBandDataType = GDT_Float32;
// If converting this image to FLOAT, then update the output format type.
if (convertToFloat) theBandDataType = GDT_Float32;
// Write geospatial metadata.
char **papszOptions = nullptr;
GDALDataset *poDstDS = poDriver->Create(fileName, this->width, this->height, this->bands, theBandDataType,
papszOptions);
double adfGeoTransform[6] = {this->easting, this->gsd, 0, this->northing + this->height * this->gsd, 0,
-1 * this->gsd};
poDstDS->SetGeoTransform(adfGeoTransform);
OGRSpatialReference oSRS;
char *pszSRS_WKT = NULL;
oSRS.SetProjCS("UTM (WGS84)");
oSRS.SetUTM(abs(this->zone), (this->zone > 0));
int theZone = 255;
int theZoneIsNorthern = 255;
theZone = oSRS.GetUTMZone(&theZoneIsNorthern);
oSRS.SetWellKnownGeogCS("WGS84");
if (!egm96)
oSRS.SetVertCS("WGS 84", "World Geodetic System 1984");
else
oSRS.SetVertCS("EGM96 Geoid", "EGM96 geoid");
oSRS.SetExtension("VERT_DATUM", "PROJ4_GRIDS", "g2009conus.gtx");
oSRS.exportToWkt(&pszSRS_WKT);
poDstDS->SetProjection(pszSRS_WKT);
CPLFree(pszSRS_WKT);
// Write the image.
if (convertToFloat) {
// Write the image as FLOAT and convert values.
float *raster = new float[this->width * this->bands];
float noData = -10000.0;
for (unsigned int i = 1; i <= this->bands; i++) {
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++) {
if (this->data[j][k] == 0)
raster[k] = noData;
else
raster[k] = (float(this->data[j][k]) * this->scale) + this->offset;
}
poBand->RasterIO(GF_Write, 0, j, this->width, 1, raster, this->width, 1, theBandDataType,
sizeof(float) * this->bands,
this->width * this->bands * sizeof(float));
}
}
delete[]raster;
} else {
// Write the image without conversion.
for (unsigned int i = 1; i <= this->bands; i++) {
GDALRasterBand *poBand = poDstDS->GetRasterBand(i);
for (unsigned int j = 0; j < this->height; j++) {
poBand->RasterIO(GF_Write, 0, j, this->width, 1, this->data[j], this->width, 1, theBandDataType,
sizeof(TYPE) * this->bands,
this->width * this->bands * sizeof(TYPE));
}
}
}
GDALClose((GDALDatasetH) poDstDS);
return true;
}
// Read image from point cloud file.
bool readFromPointCloud(char *fileName, float gsdMeters, MIN_MAX_TYPE mode = MIN_VALUE) {
// Read a PSET file (e.g., BPF or LAS).
PointCloud pset;
bool ok = pset.Read(fileName);
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 (int i = 0; i < pset.numPoints; i++) {
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);
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);
if ((x < 0) || (x > this->width - 1)) continue;
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;
}
}
return true;
}
// Count voids in an image.
// Note that voids are always labeled zero in this data structure.
long countVoids() {
long count = 0;
for (unsigned int j = 0; j < this->height; j++) {
for (unsigned int i = 0; i < this->width; i++) {
if (this->data[j][i] == 0) {
count++;
}
}
}
return count;
}
// Fill any voids in the image using a simple multigrid scheme.
// Note that voids are always labeled zero.
// MaxLevel by default is the maximum value of int
void fillVoidsPyramid(bool noSmoothing, unsigned int maxLevel = MAX_INT) {
// Check for voids.
long count = countVoids();
if (count == 0) return;
// Create image pyramid.
std::vector<OrthoImage<TYPE> *> pyramid;
pyramid.push_back(this);
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;
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);
// 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++) {
if (pyramid[level]->data[jj][ii] != 0) {
z += pyramid[level]->data[jj][ii];
ct++;
}
}
}
if (ct != 0) {
z = z / ct;
newImagePtr->data[j][i] = (TYPE) z;
}
}
}
pyramid.push_back(newImagePtr);
level++;
count = pyramid[level]->countVoids();
}
// 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++) {
// 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);
if (noSmoothing) {
// Just use the closest pixel from above.
pyramid[k]->data[j][i] = pyramid[k + 1]->data[j2][i2];
} else {
// Average neighboring pixels from above.
float z = 0;
int ct = 0;
for (int jj = MAX(0, j2 - 1);
jj <= MIN(j2 + 1, pyramid[k + 1]->height - 1); jj++) {
for (int ii = MAX(0, i2 - 1);
ii <= MIN(i2 + 1, pyramid[k + 1]->width - 1); ii++) {
z += pyramid[k + 1]->data[jj][ii];
ct++;
}
}
z = z / ct;
pyramid[k]->data[j][i] = (TYPE) z;
}
}
}
}
}
// Deallocate memory for all but the input DSM.
for (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++) {
// 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);
// 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++) {
if (this->data[jj][ii] != 0) {
values.push_back(this->data[jj][ii]);
}
}
}
// Find the median value.
unsigned long num = values.size();
if (num > 0) {
std::sort(values.begin(), values.end());
unsigned short medianValue = values[num / 2];
if (fabs(float(medianValue) - float(this->data[j][i])) > dzShort)
this->data[j][i] = medianValue;
}
}
}
}
// Apply a minimum filter to an image.
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++) {
tempImage.data[j][i] = this->data[j][i];
}
}
for (int j = 0; j < this->height; j++) {
for (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);
// Check all neighbors.
for (int jj = j1; jj <= j2; jj++) {
for (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];
}
}
}
}
}
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
this->data[j][i] = tempImage.data[j][i];
}
}
}
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++) {
tempImage.data[j][i] = this->data[j][i];
}
}
for (int j = 0; j < this->height; j++) {
for (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);
// If there's an edge with any neighbor, then remove.
for (int jj = j1; jj <= j2; jj++) {
for (int ii = i1; ii <= i2; ii++) {
if (fabs(float(this->data[jj][ii] - this->data[j][i])) > dzShort) {
tempImage.data[j][i] = 0;
}
}
}
}
}
for (int j = 0; j < this->height; j++) {
for (int i = 0; i < this->width; i++) {
this->data[j][i] = tempImage.data[j][i];
}
}
}
};
}
#endif //PUBGEO_ORTHO_IMAGE_H

27
src/common/util.h Normal file
View File

@@ -0,0 +1,27 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
//
// Created by almess1 on 4/20/17.
//
#ifndef PUBGEO_UTIL_H
#define PUBGEO_UTIL_H
#include <limits>
namespace pubgeo {
template<typename T>
inline T &Max(const T &a, const T &b) {
return a < b ? b : a;
}
template<typename T>
inline T &Min(const T &a, const T &b) {
return a > b ? b : a;
}
const int MAX_INT = std::numeric_limits<int>::max();
const float MAX_FLOAT = std::numeric_limits<float>::max();
}
#endif //PUBGEO_UTIL_H

25
src/shr3d/CMakeLists.txt Normal file
View File

@@ -0,0 +1,25 @@
SET(SHR3D_HEADER_FILES
shr3d.h)
SET(SHR3D_SOURCE_FILES
shr3d.cpp)
ADD_LIBRARY(SHR3D_LIB STATIC ${SHR3D_HEADER_FILES} ${SHR3D_SOURCE_FILES})
TARGET_LINK_LIBRARIES(SHR3D_LIB
PUBLIC
PUBGEO_LIB)
### Add some executables to play with.
### Main exe will showcase all features
ADD_EXECUTABLE(shr3d main.cpp)
TARGET_LINK_LIBRARIES(shr3d
PUBLIC
SHR3D_LIB)
IF(TESTING)
## This executable will be used to test IO
ADD_EXECUTABLE(shr3dTest test-main.cpp)
TARGET_LINK_LIBRARIES(shr3dTest
PUBLIC
SHR3D_LIB)
ENDIF()

265
src/shr3d/main.cpp Normal file
View File

@@ -0,0 +1,265 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
#include <cstdio>
#include <ctime>
#include "orthoimage.h"
#include "shr3d.h"
// Print command line arguments.
void printArguments() {
printf("Command line arguments: <Input File (LAS|TIF)> <Options>\n");
printf("Options:\n");
printf(" DH= horizontal uncertainty (meters)\n");
printf(" DZ= vertical uncertainty (meters)\n");
printf(" AGL= minimum building height above ground level (meters)\n");
printf(" AREA= minimum building area (meters)\n");
printf(" EGM96 set this flag to write vertical datum = EGM96\n");
printf(" FOPEN set this flag for ground detection under dense foliage\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");
printf(" For lidar LAS: shr3d pts.las DH=1.0 DZ=1.0 AGL=2.0 AREA=50.0\n");
}
// Main program for bare earth classification.
int main(int argc, char **argv) {
// If no parameters, then print command line arguments.
if (argc < 4) {
printf("Number of arguments = %d\n", argc);
for (int i = 0; i < argc; i++) {
printf("ARG[%d] = %s\n", i, argv[i]);
}
printArguments();
return -1;
}
// Get command line arguments and confirm they are valid.
double gsd_meters = 0.0;
double dh_meters = 0.0;
double dz_meters = 0.0;
double agl_meters = 0.0;
double min_area_meters = 50.0;
bool egm96 = false;
bool convert = false;
char inputFileName[1024];
sprintf(inputFileName, argv[1]);
for (int i = 2; i < argc; i++) {
if (strstr(argv[i], "DH=")) { dh_meters = atof(&(argv[i][3])); }
if (strstr(argv[i], "DZ=")) { dz_meters = atof(&(argv[i][3])); }
if (strstr(argv[i], "AGL=")) { agl_meters = atof(&(argv[i][4])); }
if (strstr(argv[i], "AREA=")) { min_area_meters = atof(&(argv[i][5])); }
if (strstr(argv[i], "EGM96")) { egm96 = true; }
if (strstr(argv[i], "CONVERT")) { convert = true; }
}
if ((dh_meters == 0.0) || (dz_meters == 0.0) || (agl_meters == 0.0)) {
printf("DH_METERS = %f\n", dh_meters);
printf("DZ_METERS = %f\n", dz_meters);
printf("AGL_METERS = %f\n", agl_meters);
printf("Error: All three values must be nonzero.\n");
printArguments();
return -1;
}
// Initialize the timer.
time_t t0;
time(&t0);
// If specified, then convert to GDAL TIFF.
char readFileName[1024];
strcpy(readFileName, inputFileName);
if (convert) {
char cmd[4096];
sprintf(readFileName, "temp.tif\0");
sprintf(cmd, ".\\gdal\\gdal_translate %s temp.tif\n", inputFileName);
system(cmd);
}
// Read DSM as SHORT.
// Input can be GeoTIFF or LAS.
shr3d::OrthoImage<unsigned short> dsmImage;
shr3d::OrthoImage<unsigned short> minImage;
printf("Reading DSM as SHORT.\n");
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;
} else if ((strcmp(ext, "las") == 0) || (strcmp(ext, "bpf") == 0)) {
// First get the max Z values for the DSM.
bool ok = dsmImage.readFromPointCloud(inputFileName, (float) dh_meters, shr3d::MAX_VALUE);
if (!ok) return -1;
// Median filter, replacing only points differing by more than the AGL threshold.
dsmImage.medianFilter(1, (unsigned int) (agl_meters / dsmImage.scale));
// Fill small voids in the DSM.
dsmImage.fillVoidsPyramid(true, 2);
// Write the DSM image as FLOAT.
char dsmOutFileName[1024];
sprintf(dsmOutFileName, "%s_DSM.tif\0", inputFileName);
dsmImage.write(dsmOutFileName, true);
// Now get the minimum Z values for the DTM.
ok = minImage.readFromPointCloud(inputFileName, (float) dh_meters, shr3d::MIN_VALUE);
if (!ok) return -1;
// Median filter, replacing only points differing by more than the AGL threshold.
minImage.medianFilter(1, (unsigned int) (agl_meters / minImage.scale));
// Fill small voids in the DSM.
minImage.fillVoidsPyramid(true, 2);
// Write the MIN image as FLOAT.
char minOutFileName[1024];
sprintf(minOutFileName, "%s_MIN.tif\0", inputFileName);
minImage.write(minOutFileName, true);
// 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++) {
float minValue = (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 (((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 = 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++) {
float diff = (float) dsmImage.data[jj][ii] - minValue;
if (diff < threshold) found = true;
}
}
if (!found) dsmImage.data[j][i] = 0;
}
}
}
// Write the DSM2 image as FLOAT.
char dsm2OutFileName[1024];
sprintf(dsm2OutFileName, "%s_DSM2.tif\0", inputFileName);
dsmImage.write(dsm2OutFileName, true);
} else {
printf("Error: Unrecognized file type.");
return -1;
}
// Convert horizontal and vertical uncertainty values to bin units.
int dh_bins = MAX(1, (int) floor(dh_meters / dsmImage.gsd));
printf("DZ_METERS = %f\n", dz_meters);
printf("DH_METERS = %f\n", dh_meters);
printf("DH_BINS = %d\n", dh_bins);
unsigned int dz_short = (unsigned int) (dz_meters / dsmImage.scale);
printf("DZ_SHORT = %d\n", dz_short);
printf("AGL_METERS = %f\n", agl_meters);
unsigned int agl_short = (unsigned int) (agl_meters / dsmImage.scale);
printf("AGL_SHORT = %d\n", agl_short);
printf("AREA_METERS = %f\n", min_area_meters);
// 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 (int j = 0; j < dsmImage.height; j++) {
for (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.
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++) {
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, (unsigned int) (dz_meters / dsmImage.scale));
// Refine the object label image and export building outlines.
shr3d::Shr3dder::classifyNonGround(dsmImage, dtmImage, labelImage, dz_short, agl_short, (float) min_area_meters);
// Fill small voids in the DTM after all processing is complete.
dtmImage.fillVoidsPyramid(true, 2);
// Write the DTM image as FLOAT.
char dtmOutFileName[1024];
sprintf(dtmOutFileName, "%s_DTM.tif\0", inputFileName);
dtmImage.write(dtmOutFileName, true, egm96);
// Produce a classification raster image with LAS standard point classes.
shr3d::OrthoImage<unsigned char> classImage;
classImage.Allocate(labelImage.width, labelImage.height);
classImage.easting = dsmImage.easting;
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++) {
// Set default as unlabeled.
classImage.data[j][i] = LAS_UNCLASSIFIED;
// Label trees.
if ((dsmImage.data[j][i] == 0) ||
(fabs((float) dsmImage.data[j][i] - (float) dtmImage.data[j][i]) > agl_short))
classImage.data[j][i] = LAS_TREE;
// Label buildings.
if (labelImage.data[j][i] == 1) classImage.data[j][i] = LAS_BUILDING;
// Label ground.
if (fabs((float) dsmImage.data[j][i] - (float) dtmImage.data[j][i]) < dz_short)
classImage.data[j][i] = LAS_GROUND;
}
}
// Fill missing labels inside building regions.
shr3d::Shr3dder::fillInsideBuildings(classImage);
// Write the classification image.
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++) {
if (classImage.data[j][i] != LAS_BUILDING) classImage.data[j][i] = 0;
}
}
sprintf(classOutFileName, "%s_buildings.tif\0", inputFileName);
classImage.write(classOutFileName, false, egm96);
// Report total elapsed time.
time_t t1;
time(&t1);
printf("Total time elapsed = %f seconds\n", double(t1 - t0));
}

727
src/shr3d/shr3d.cpp Normal file
View File

@@ -0,0 +1,727 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
// shr3d.cpp
//
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <climits>
#include <vector>
#include "shr3d.h"
using namespace shr3d;
// Extend object boundaries to capture points missed around the edges.
void extendObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned long> &labelImage,
int edgeResolution, unsigned int minDistanceShortValue) {
// Loop enough to capture the edge resolution.
for (unsigned int k = 0; k < edgeResolution; k++) {
// First, label any close neighbor LABEL_TEMP.
for (unsigned int j = 1; j < labelImage.height - 1; j++) {
for (unsigned int i = 1; i < labelImage.width - 1; i++) {
// For any labeled point, check all neighbors.
if (labelImage.data[j][i] == 1) {
for (unsigned int jj = j - 1; jj <= j + 1; jj++) {
for (unsigned int ii = i - 1; ii <= i + 1; ii++) {
if (labelImage.data[jj][ii] == 1) continue;
if (((float) dsmImage.data[j][i] - (float) dsmImage.data[jj][ii]) <
minDistanceShortValue / 2.0) {
labelImage.data[jj][ii] = LABEL_TEMP;
}
}
}
}
}
}
// Then label any high points labeled LABEL_TEMP as an object of interest.
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
if (labelImage.data[j][i] == LABEL_TEMP) {
// Check to make sure this point is also higher than one of its neighbors.
unsigned int j1 = MAX(0, j - 1);
unsigned int j2 = MIN(j + 1, labelImage.height - 1);
unsigned int i1 = MAX(0, i - 1);
unsigned int i2 = MIN(i + 1, labelImage.width - 1);
for (unsigned int jj = j1; jj <= j2; jj++) {
for (unsigned int ii = i1; ii <= i2; ii++) {
if (((float) dsmImage.data[j][i] - (float) dsmImage.data[jj][ii]) >
minDistanceShortValue / 2.0) {
labelImage.data[j][i] = 1;
}
}
}
}
}
}
}
// Reset any temporary values.
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
if (labelImage.data[j][i] == LABEL_TEMP) {
labelImage.data[j][i] = LABEL_GROUND;
}
}
}
}
// Label boundaries of objects above ground level.
void labelObjectBoundaries(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned long> &labelImage,
int edgeResolution, unsigned int minDistanceShortValue) {
// Initialize the labels to LABEL_GROUND.
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
labelImage.data[j][i] = LABEL_GROUND;
}
}
// Mark the label image with object boundaries.
int radius = 2;
float threshold = (float) minDistanceShortValue;
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
// Look for Z steps greater than a threshold.
float value = (float) dsmImage.data[j][i];
// Interestingly, this works about as well as checking every step.
for (int dj = -edgeResolution; dj <= edgeResolution; dj += edgeResolution) {
for (int di = -edgeResolution; di <= edgeResolution; di += edgeResolution) {
int j2 = MIN(MAX(0, j + dj), dsmImage.height - 1);
int i2 = MIN(MAX(0, i + di), dsmImage.width - 1);
if (dsmImage.data[j2][i2] != 0) {
// Remove local slope to avoid tagging rough terrain.
int j3 = MIN(MAX(0, j + dj * 2), dsmImage.height - 1);
int i3 = MIN(MAX(0, i + di * 2), dsmImage.width - 1);
float myGradient = (float) dsmImage.data[j][i] - (float) dsmImage.data[j2][i2];
float neighborGradient = (float) dsmImage.data[j2][i2] - (float) dsmImage.data[j3][i3];
float distance = (myGradient - neighborGradient);
if (distance > threshold) labelImage.data[j][i] = 1;
}
}
}
}
}
}
// Fill inside the object countour labels if points are above the nearby ground level.
void fillObjectBounds(OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned short> &dsmImage, ObjectType &obj,
int edgeResolution, unsigned int dzShort) {
unsigned int label = obj.label;
// Loop on rows, filling in labels.
for (unsigned int j = MAX(0, obj.ymin - 1); j <= MIN(obj.ymax + 1, labelImage.height - 1); j++) {
// Get start index.
int startIndex = -1;
for (unsigned int i = MAX(0, obj.xmin - 1); i <= MIN(obj.xmax + 1, labelImage.width - 1); i++) {
if (labelImage.data[j][i] == label) {
startIndex = i;
break;
}
}
// If no labels in this row, then continue.
if (startIndex == -1) continue;
// Get stop index.
int stopIndex = -1;
for (unsigned int i = MIN(obj.xmax + 1, labelImage.width - 1); i >= MAX(0, obj.xmin - 1); i--) {
if (labelImage.data[j][i] == label) {
stopIndex = i;
break;
}
}
// If entire column is labeled, then continue.
if ((startIndex == 0) && (stopIndex == 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.
unsigned short groundLevel;
if (startIndex == 0)
groundLevel = dsmImage.data[j][stopIndex + 1];
else if (stopIndex == labelImage.width - 1)
groundLevel = dsmImage.data[j][startIndex - 1];
else
groundLevel = MAX(dsmImage.data[j][startIndex - 1], dsmImage.data[j][stopIndex + 1]);
// Fill in the label for any point in between that has height above ground level.
for (unsigned int i = startIndex; i <= stopIndex; i++) {
if (dsmImage.data[j][i] > groundLevel) {
if (labelImage.data[j][i] != label) labelImage.data[j][i] = LABEL_IN_ONE;
} else {
if (labelImage.data[j][i] == label) labelImage.data[j][i] = LABEL_GROUND;
}
}
}
// Loop on columns, filling in labels.
for (unsigned int i = MAX(0, obj.xmin - 1); i <= MIN(obj.xmax + 1, labelImage.width - 1); i++) {
// Get start index.
int startIndex = -1;
for (unsigned int j = MAX(0, obj.ymin - 1); j <= MIN(obj.ymax + 1, labelImage.height - 1); j++) {
if (labelImage.data[j][i] == label) {
startIndex = j;
break;
}
}
// If no labels in this row, then continue.
if (startIndex == -1) continue;
// Get stop index.
int stopIndex = -1;
for (unsigned int j = MIN(obj.ymax + 1, labelImage.height - 1); j >= MAX(0, obj.ymin - 1); j--) {
if (labelImage.data[j][i] == label) {
stopIndex = j;
break;
}
}
// If entire row is labeled, then continue.
if ((startIndex == 0) && (stopIndex == labelImage.height-1)) continue;
// Get max ground level height for this row.
unsigned short groundLevel;
if (startIndex == 0)
groundLevel = dsmImage.data[stopIndex + 1][i];
else if (stopIndex == labelImage.height - 1)
groundLevel = dsmImage.data[startIndex - 1][i];
else
groundLevel = MAX(dsmImage.data[startIndex - 1][i], dsmImage.data[stopIndex + 1][i]);
// Fill in the label for any point in between that has height above ground level.
// This time make sure both the horizontal and vertical check pass and set to LABEL_ACCEPTED.
for (unsigned int j = startIndex; j <= stopIndex; j++) {
if (dsmImage.data[j][i] > groundLevel) {
if ((labelImage.data[j][i] == label) || (labelImage.data[j][i] == LABEL_IN_ONE)) {
labelImage.data[j][i] = LABEL_ACCEPTED;
}
}
}
}
// Erode the labels with a kernel size based on edge resolution.
int rad = edgeResolution;
for (unsigned int j = MAX(0, obj.ymin - 1); j <= MIN(obj.ymax + 1, labelImage.height - 1); j++) {
for (unsigned int i = MAX(0, obj.xmin - 1); i <= MIN(obj.xmax + 1, labelImage.width - 1); i++) {
if (labelImage.data[j][i] == LABEL_ACCEPTED) {
int i1 = MAX(i - rad, 0);
int i2 = MIN(i + rad, labelImage.width - 1);
int j1 = MAX(j - rad, 0);
int j2 = MIN(j + rad, labelImage.height - 1);
for (unsigned int jj = j1; jj <= j2; jj++) {
for (unsigned int ii = i1; ii <= i2; ii++) {
if (labelImage.data[jj][ii] != LABEL_ACCEPTED) labelImage.data[jj][ii] = LABEL_TEMP;
}
}
}
}
}
// Update object bounds to include erosion.
obj.xmin = MAX(0, obj.xmin - edgeResolution - 1);
obj.ymin = MAX(0, obj.ymin - edgeResolution - 1);
obj.xmax = MIN(obj.xmax + edgeResolution + 1, labelImage.width - 1);
obj.ymax = MIN(obj.ymax + edgeResolution + 1, labelImage.height - 1);
for (unsigned int j = obj.ymin; j <= obj.ymax; j++) {
for (unsigned int i = obj.xmin; i <= obj.xmax; i++) {
if (labelImage.data[j][i] == LABEL_TEMP) {
labelImage.data[j][i] = LABEL_ACCEPTED;
}
}
}
// Finish up the labels.
for (unsigned int j = MAX(0, obj.ymin - 1); j <= MIN(obj.ymax + 1, labelImage.height - 1); j++) {
for (unsigned int i = MAX(0, obj.xmin - 1); i <= MIN(obj.xmax + 1, labelImage.width - 1); i++) {
if (labelImage.data[j][i] == label) labelImage.data[j][i] = LABEL_GROUND;
if (labelImage.data[j][i] == LABEL_ACCEPTED) labelImage.data[j][i] = LABEL_OBJECT;
}
}
}
// 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) {
// 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;
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++) {
// 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;
// Skip if height is too different.
if (fabs((float) dsmImage.data[jj][ii] - (float) dsmImage.data[j][i]) > (float) dzShort) continue;
// Update the label.
labelImage.data[jj][ii] = labelImage.data[j][i];
// Add to the new list.
PixelType pixel = {ii, jj};
newNeighbors.push_back(pixel);
// Update object bounds.
obj.xmin = MIN(obj.xmin, ii);
obj.xmax = MAX(obj.xmax, ii);
obj.ymin = MIN(obj.ymin, jj);
obj.ymax = MAX(obj.ymax, jj);
obj.count++;
}
}
}
// Replace the neighbors list with the new one. No need to merge them.
if (newNeighbors.size() > 0) {
neighbors = newNeighbors;
return true;
} else return false;
}
// Group connected labeled pixels into objects.
void groupObjects(OrthoImage<unsigned long> &labelImage, 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;
unsigned int label = 1;
int bigCount = 0;
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
// Skip unlabeled pixels.
if (labelImage.data[j][i] == LABEL_GROUND) continue;
// Skip already labeled pixels.
if (labelImage.data[j][i] > 1) continue;
// Create a new label.
label++;
// Create a new object structure.
ObjectType obj;
// Initialize the object and neighbors structures.
obj.label = label;
obj.xmin = i;
obj.ymin = j;
obj.xmax = i;
obj.ymax = j;
obj.count = 1;
std::vector<PixelType> neighbors;
PixelType pixel = {i, j};
neighbors.push_back(pixel);
labelImage.data[j][i] = label;
// Get points in new group.
bool keepSearching = true;
while (keepSearching) {
keepSearching = addNeighbors(neighbors, labelImage, dsmImage, obj, dzShort);
// This is very quick but not especially smart.
// Try to find a reasonably quick way to split the regions more sensibly.
// That said, this also seems to work well enough.
if (obj.count > maxCount) {
bigCount++;
break;
}
}
// Add this object to the list.
objects.push_back(obj);
maxGroupSize = MAX(maxGroupSize, obj.count);
}
}
printf("Max group size = %ld\n", maxGroupSize);
printf("Number of cropped groups = %d\n", bigCount);
}
// 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) {
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)) {
labelImage.data[j][i] = 1;
}
}
}
}
// Classify ground points, fill the voids, and generate a bare earth terrain model.
void Shr3dder::classifyGround(OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned short> &dsmImage,
OrthoImage<unsigned short> &dtmImage, int dhBins, unsigned int dzShort) {
// Fill voids.
printf("Filling voids...\n");
dtmImage.fillVoidsPyramid(true);
// Allocate a binary label image to indicate voids to be filled.
// The long integer label image has unique labels for objects detected in each iteration.
OrthoImage<unsigned char> voidImage;
voidImage.Allocate(labelImage.width, labelImage.height);
// Iteratively label and remove objects from the DEM.
// Each new iteration removes debris not identified by the previous iteration.
int numIterations = 5;
long maxCount = 10000 / (dsmImage.gsd * dsmImage.gsd); // max count is in meters
for (unsigned int k = 0; k < numIterations; k++) {
printf("Iteration #%d\n", k + 1);
// Label the object boundaries.
printf("Labeling object boundaries...\n");
labelObjectBoundaries(dtmImage, labelImage, dhBins, dzShort);
// Extend labels for object boundaries..
printf("Extending object boundaries...\n");
extendObjectBoundaries(dtmImage, labelImage, dhBins, dzShort);
// Group the objects.
printf("Grouping objects...\n");
std::vector<ObjectType> objects;
groupObjects(labelImage, dtmImage, objects, maxCount, dzShort);
printf("Number of objects = %ld\n", objects.size());
// Generate object groups and void fill them in the DEM image.
printf("Labeling and removing objects...\n");
for (long i = 0; i < objects.size(); i++) {
fillObjectBounds(labelImage, dtmImage, objects[i], dhBins, dzShort);
}
// Update the label image values for easy viewing.
printf("Finishing label image for display...\n");
finishLabelImage(labelImage);
// Update the void image.
printf("Updating void image...\n");
for (unsigned int j = 0; j < voidImage.height; j++) {
for (unsigned int i = 0; i < voidImage.width; i++) {
if (labelImage.data[j][i] == 1) voidImage.data[j][i] = 1;
}
}
// Fill voids.
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
if (voidImage.data[j][i] == 1) dtmImage.data[j][i] = 0;
}
}
bool noSmoothing = true;
if (k == numIterations - 1) noSmoothing = false;
printf("Filling voids with noSmoothing = %d\n", noSmoothing);
dtmImage.fillVoidsPyramid(noSmoothing);
}
// If any DTM points are above the DSM, then restore the DSM values.
for (unsigned int j = 0; j < dtmImage.height; j++) {
for (unsigned int i = 0; i < dtmImage.width; i++) {
if (dtmImage.data[j][i] >= dsmImage.data[j][i]) {
dtmImage.data[j][i] = dsmImage.data[j][i];
labelImage.data[j][i] = LABEL_GROUND;
voidImage.data[j][i] = 0;
}
}
}
// Remove any leftover single point spikes.
printf("Removing spikes...\n");
for (unsigned int j = 0; j < dtmImage.height; j++) {
for (unsigned int i = 0; i < dtmImage.width; i++) {
float minDiff = FLT_MAX;
for (int jj = -1; jj <= 1; jj++) {
int j2 = MAX(0, MIN(dsmImage.height - 1, j + jj));
for (int ii = -1; ii <= 1; ii++) {
if ((ii == 0) && (jj == 0)) continue;
int i2 = MAX(0, MIN(dsmImage.width - 1, i + ii));
float diff = MAX(0, (float) dtmImage.data[j][i] - (float) dtmImage.data[j2][i2]);
minDiff = MIN(minDiff, diff);
}
}
if (minDiff > dzShort / 2.0) {
labelImage.data[j][i] = 1;
voidImage.data[j][i] = 1;
dtmImage.data[j][i] = 0.0;
}
}
}
// Fill voids.
printf("Filling voids...\n");
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
if (voidImage.data[j][i] == 1) dtmImage.data[j][i] = 0;
}
}
dtmImage.fillVoidsPyramid(false);
// Mark all voids.
printf("Marking voids in label image after all iterations are complete...\n");
for (unsigned int j = 0; j < voidImage.height; j++) {
for (unsigned int i = 0; i < voidImage.width; i++) {
if (voidImage.data[j][i] == 1) labelImage.data[j][i] = 1;
else labelImage.data[j][i] = LABEL_GROUND;
}
}
}
// 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,
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.
int minPointCount = int(minAreaMeters / (dsmImage.gsd * dsmImage.gsd));
printf("Min points for removing small objects = %d\n", minPointCount);
// Apply AGL threshold to individual points to reduce clutter.
// Note that ground level clutter tends to be less than 2m AGL.
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) {
if (dsmImage.data[j][i] == 0) labelImage.data[j][i] = LABEL_GROUND;
else if (((float) dsmImage.data[j][i] - (float) dtmImage.data[j][i]) < (float) aglShort)
labelImage.data[j][i] = LABEL_GROUND;
}
}
}
// Group the labeled objects based on height similarity.
{
std::vector<ObjectType> objects;
groupObjects(labelImage, dsmImage, objects, LONG_MAX, dzShort / 2);
// Accept or reject each object independently based on boundary values and flatness.
int numRejected = 0;
for (unsigned int k = 0; k < objects.size(); k++) {
bool reject = false;
// Reject if boundary gradients with neighbors labeled ground are too small.
float meanGradient = 0.0;
int count = 0;
for (int j = objects[k].ymin; j <= objects[k].ymax; j++) {
for (unsigned int i = objects[k].xmin; i <= objects[k].xmax; i++) {
if (labelImage.data[j][i] == objects[k].label) {
for (int jj = -1; jj <= 1; jj++) {
unsigned int j2 = MAX(0, MIN(dsmImage.height - 1, j + jj));
for (int ii = -1; ii <= 1; ii++) {
unsigned int i2 = MAX(0, MIN(dsmImage.width - 1, i + ii));
//if (labelImage.data[j2][i2] != objects[k].label)
if (labelImage.data[j2][i2] == LABEL_GROUND) {
unsigned int j3 = MAX(0, MIN(dsmImage.height - 1, j + jj * 2));
unsigned int i3 = MAX(0, MIN(dsmImage.width - 1, i + ii * 2));
// These assume I'm higher than my neighbors.
float myGradient = MAX(0, ((float) dsmImage.data[j][i] -
(float) dsmImage.data[j2][i2]));
float neighborGradient = MAX(0, ((float) dsmImage.data[j2][i2] -
(float) dsmImage.data[j3][i3]));
meanGradient += MAX(0, (myGradient - neighborGradient));
count++;
}
}
}
}
}
}
meanGradient /= count;
if ((meanGradient < dzShort / 2.0) && (meanGradient != 0.0)) reject = true;
// If rejected, then label each point.
if (reject) {
numRejected++;
for (unsigned int j = objects[k].ymin; j <= objects[k].ymax; j++) {
for (unsigned int i = objects[k].xmin; i <= objects[k].xmax; i++) {
if (labelImage.data[j][i] == objects[k].label) {
labelImage.data[j][i] = LABEL_GROUND;
}
}
}
}
}
}
// Erode and then dilate labels to remove narrow objects.
{
OrthoImage<unsigned long> tempImage;
tempImage.Allocate(labelImage.width, labelImage.height);
for (unsigned int j = 0; j < labelImage.height; j++) {
for (unsigned int i = 0; i < labelImage.width; i++) {
tempImage.data[j][i] = labelImage.data[j][i];
}
}
for (int j = 0; j < labelImage.height; j++) {
for (int i = 0; i < labelImage.width; i++) {
if (labelImage.data[j][i] != LABEL_GROUND) {
// Define neighbor bounds;
int i1 = MAX(0, i - 1);
int i2 = MIN(i + 1, labelImage.width - 1);
int j1 = MAX(0, j - 1);
int j2 = MIN(j + 1, labelImage.height - 1);
// Unlabel any point with an unlabeled neighbor.
for (int jj = j1; jj <= j2; jj++) {
for (int ii = i1; ii <= i2; ii++) {
if (labelImage.data[jj][ii] == LABEL_GROUND) {
tempImage.data[j][i] = LABEL_GROUND;
}
}
}
}
}
}
for (int j = 0; j < labelImage.height; j++) {
for (int i = 0; i < labelImage.width; i++) {
if (labelImage.data[j][i] != LABEL_GROUND) {
// Define neighbor bounds;
int i1 = MAX(0, i - 1);
int i2 = MIN(i + 1, labelImage.width - 1);
int j1 = MAX(0, j - 1);
int j2 = MIN(j + 1, labelImage.height - 1);
// Unlabel any point with no labeled neighbors after erosion.
bool found = false;
for (int jj = j1; jj <= j2; jj++) {
for (int ii = i1; ii <= i2; ii++) {
if (tempImage.data[jj][ii] != LABEL_GROUND) found = true;
}
}
if (!found) labelImage.data[j][i] = LABEL_GROUND;
}
}
}
}
// Reset all non-ground labels to one.
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) labelImage.data[j][i] = 1;
}
}
// Group labeled points to remove small objects.
// Do not split groups based on height or point count.
{
printf("Grouping to remove small objects...\n");
std::vector<ObjectType> objects;
groupObjects(labelImage, dsmImage, objects, LONG_MAX, INT_MAX);
int numRejected = 0;
for (unsigned int k = 0; k < objects.size(); k++) {
bool reject = false;
// Reject if too small.
if (objects[k].count < minPointCount) reject = true;
// If rejected, then label each point.
if (reject) {
numRejected++;
for (unsigned int j = objects[k].ymin; j <= objects[k].ymax; j++) {
for (unsigned int i = objects[k].xmin; i <= objects[k].xmax; i++) {
if (labelImage.data[j][i] == objects[k].label) {
labelImage.data[j][i] = LABEL_GROUND;
}
}
}
}
}
printf("Number of small objects rejected = %d\n", numRejected);
}
// Reset all non-ground labels to one.
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) labelImage.data[j][i] = 1;
}
}
}
// Add neighboring pixels to a void group.
bool
addClassNeighbors(std::vector<PixelType> &neighbors, OrthoImage<unsigned char> &classImage,
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++) {
// If already labeled, then skip.
if (labeled.data[jj][ii] == 1) continue;
// If not the same label, then skip.
if (classImage.data[jj][ii] != label) continue;
// Update the label.
labeled.data[jj][ii] = 1;
// Add to the new list.
PixelType pixel = {ii, jj};
newNeighbors.push_back(pixel);
}
}
}
// Update the neighbors list.
if (newNeighbors.size() > 0) {
for (int k = 0; k < newNeighbors.size(); k++) {
neighbors.push_back(newNeighbors[k]);
}
// neighbors = newNeighbors;
return true;
} else return false;
}
// Fill in any pixels labeled tree that fall entirely within a labeled building group.
void Shr3dder::fillInsideBuildings(OrthoImage<unsigned char> &classImage) {
int numFilled = 0;
OrthoImage<unsigned char> labeled;
labeled.Allocate(classImage.width, classImage.height);
for (unsigned int j = 0; j < classImage.height; j++) {
for (unsigned int i = 0; i < classImage.width; i++) {
bool consider = ((labeled.data[j][i] == 0) && (classImage.data[j][i] == LAS_TREE));
if (!consider) continue;
// Get all pixels in this contiguous group.
std::vector<PixelType> neighbors;
PixelType pixel = {i, j};
neighbors.push_back(pixel);
bool keepSearching = true;
unsigned int label = (unsigned int) classImage.data[j][i];
labeled.data[j][i] = 1;
while (keepSearching) {
keepSearching = addClassNeighbors(neighbors, classImage, labeled, label);
}
// Check all pixels for neighboring labels.
bool inside = true;
for (int k = 0; k < neighbors.size(); k++) {
long i1 = neighbors[k].i;
long j1 = neighbors[k].j;
for (long jj = MAX(0, j1 - 1); jj <= MIN(j1 + 1, classImage.height - 1); jj++) {
for (long ii = MAX(0, i1 - 1); ii <= MIN(i1 + 1, classImage.width - 1); ii++) {
if ((labeled.data[jj][ii] == 0) && (classImage.data[jj][ii] != LAS_BUILDING)) {
inside = false;
}
}
}
}
// If the region is completely inside a building region, then fill it.
if (inside) {
for (int k = 0; k < neighbors.size(); k++) {
numFilled++;
classImage.data[neighbors[k].j][neighbors[k].i] = LAS_BUILDING;
}
}
}
}
printf("Removed %d tree pixels inside building label groups.\n", numFilled);
}

76
src/shr3d/shr3d.h Normal file
View File

@@ -0,0 +1,76 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
// shr3d.h
//
#ifndef PUBGEO_SHR3D_H
#define PUBGEO_SHR3D_H
#include <cstdio>
#include <cstring>
#include "orthoimage.h"
#define LABEL_OBJECT 1
#define LABEL_GROUND UINT_MAX
#define LABEL_IN_ONE UINT_MAX-1
#define LABEL_ACCEPTED UINT_MAX-2
#define LABEL_TEMP UINT_MAX-3
#define LABEL_BOUNDARY UINT_MAX-5
/*
* Las Classifications list
0 Created, never classified
1 Unclassified1
2 Ground
3 Low Vegetation
4 Medium Vegetation
5 High Vegetation
6 Building
7 Low Point (noise)
8 Model Key-point (mass point)
9 Water
10 Reserved for ASPRS Definition
11 Reserved for ASPRS Definition
12 Overlap Points2
13-31 Reserved for ASPRS Definition
*/
#define LAS_UNCLASSIFIED (1*40)
#define LAS_GROUND (2*40)
#define LAS_TREE (5*40)
#define LAS_BUILDING (6*40)
// For now, multiply these by 40 so we can view the images easily.
namespace shr3d {
using namespace pubgeo;
typedef struct {
unsigned int i;
unsigned int j;
} PixelType;
typedef struct {
unsigned int label;
long int xmin;
long int xmax;
long int ymin;
long int ymax;
long int count;
} ObjectType;
class Shr3dder {
public:
// Function declarations.
static void classifyGround(OrthoImage<unsigned long> &labelImage, OrthoImage<unsigned short> &dsmImage,
OrthoImage<unsigned short> &dtmImage, int dhBins, unsigned int dzShort);
static void classifyNonGround(OrthoImage<unsigned short> &dsmImage, OrthoImage<unsigned short> &dtmImage,
OrthoImage<unsigned long> &labelImage, unsigned int dzShort,
unsigned int aglShort,
float minAreaMeters);
static void fillInsideBuildings(OrthoImage<unsigned char> &classImage);
};
}
#endif // PUBGEO_SHR3D_H

36
src/shr3d/test-main.cpp Normal file
View File

@@ -0,0 +1,36 @@
// Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
// Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
#include <iostream>
#include <random>
#include "PointCloud.h"
int main(int argc, char *argv[]) {
const char *fileName = "test.las";
switch (argc) {
case 2:
fileName = argv[1];
case 1: {
pubgeo::PointCloud pointCloud;
std::cerr << "Processing file: " << fileName << std::endl;
if (pointCloud.Read(fileName)) {
// Grab a truly random c++ 11 number
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_int_distribution<> dis(0, pointCloud.numPoints);
long int index = dis(gen);
std::cerr << pointCloud.numPoints << " points in file." << std::endl
<< "Here is a sample point from index " << index << ":"
<< pointCloud.x(index) << ", " << pointCloud.y(index) << ", "
<< pointCloud.z(index) << std::endl;
}
break;
}
default:
std::cerr << "Invalid number of arguments" << std::endl;
}
#ifdef WIN32
system("pause");
#endif
}

36
windowsInstallation.txt Normal file
View File

@@ -0,0 +1,36 @@
These instructions are not complete and do not currently result in a working build. It is close, and subsets of the project have worked with PDAL. The final project could not be entirely integrated.
Download and run OSGeo4W64 Installer.
Select advanced installation. Everything else until packages is default.
From Command Line Utilities select:
Gdal, gdal-dev, pdal
Install packages to meet dependencies. (page after package selection)
Download and install visual studio 2017 community (REQUIRES SYSTEM-RESTART)
Install packages:
.Net desktop development (possibly unnecessary)
Desktop development with c++
** include c++/cli support, std lib mods, mfc and atl, win8.1 sdk **
https://stackoverflow.com/questions/42701019/problems-generating-solution-for-vs-2017-with-cmake
This combination of packages generates a compatible visual studio environment. Some packages may be unnecessary (mfc +atl, or .net desktop env)
Load project (folder with CmakeLists.txt) AS 64 BIT!! VS will probably default you to X86
Visual studio should start processing the cmake file
cmake will probably fail:
Missing c/cxx compiler
see above package installation, you probably missed a package that was necessary in 2017 vs installation
pdal config fails "file or directory x referenced by variable _foo does no exist ! (typo from source)
edit C:\OSGeo4W64\lib\pdal\cmake\PDALConfig.cmake
replace lines 31+ with: (Changing the OSGEO4W64_INSTALLATION_LOCATION respective to your installation directory) This is due PDAL being compiled on a different machine with a different folder structure
set(OSGEO4W64_INSTALLATION_LOCATION "C:/OSGeo4W64/")
set(PDAL_INCLUDE_DIRS "${OSGEO4W64_INSTALLATION_LOCATION}include")
set(PDAL_LIBRARY_DIRS "${OSGEO4W64_INSTALLATION_LOCATION}lib")
include("${CMAKE_CURRENT_LIST_DIR}/PDALTargets.cmake")
set(PDAL_LIBRARIES ${PDAL_LIBRARY_DIRS}/pdalcpp.lib ${PDAL_LIBRARY_DIRS}/pdal_util.lib)
check_required_components(PDAL)
These instructions will generate a linker error missing what appear to be network conversions for ints/shorts/etc referenced in PDALs BEExtractor.
Possible solutions (Scott's best guesses)
Perhaps un-defining 'WIN32_LEAN_AND_MEAN'? It was used to remove macros that were breaking the build. IE MAX(..) and MIN(..)
PDAL build for windows on OSGeo package improperly built/configured (Partly already true
Build PDAL on local machine and