Added odm_dem2points module, better skirts algo

pull/889/head
Piero Toffanin 2018-07-06 16:41:56 -04:00
rodzic aa8eeaa26c
commit 7ec33ffb09
8 zmienionych plików z 770 dodań i 59 usunięć

Wyświetl plik

@ -8,6 +8,7 @@ add_subdirectory(odm_extract_utm)
add_subdirectory(odm_georef)
add_subdirectory(odm_orthophoto)
add_subdirectory(odm_cleanmesh)
add_subdirectory(odm_dem2points)
if (ODM_BUILD_SLAM)
add_subdirectory(odm_slam)
endif ()

Wyświetl plik

@ -58,7 +58,7 @@ int main(int argc, char **argv) {
logWriter.verbose = Verbose.set;
logWriter.outputFile = "odm_cleanmesh_log.txt";
// logWriter.outputFile = "odm_cleanmesh_log.txt";
logArgs(params, logWriter);
vtkSmartPointer<vtkPLYReader> reader =

Wyświetl plik

@ -0,0 +1,19 @@
project(odm_dem2points)
cmake_minimum_required(VERSION 2.8)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR})
set (CMAKE_CXX_STANDARD 11)
find_package(GDAL REQUIRED)
include_directories(${GDAL_INCLUDE_DIR})
# Add compiler options.
add_definitions(-Wall -Wextra)
# Add source directory
aux_source_directory("./src" SRC_LIST)
# Add exectuteable
add_executable(${PROJECT_NAME} ${SRC_LIST})
target_link_libraries(${PROJECT_NAME} ${GDAL_LIBRARY})

Wyświetl plik

@ -0,0 +1,106 @@
/*
Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer. Redistributions in binary form must reproduce
the above copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the distribution.
Neither the name of the Johns Hopkins University nor the names of its contributors
may be used to endorse or promote products derived from this software without specific
prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
*/
#ifndef CMD_LINE_PARSER_INCLUDED
#define CMD_LINE_PARSER_INCLUDED
#include <stdarg.h>
#include <cstring>
#include <cstdlib>
#include <string>
#include <vector>
#ifdef WIN32
int strcasecmp( const char* c1 , const char* c2 );
#endif // WIN32
class cmdLineReadable
{
public:
bool set;
char *name;
cmdLineReadable( const char *name );
virtual ~cmdLineReadable( void );
virtual int read( char** argv , int argc );
virtual void writeValue( char* str ) const;
};
template< class Type > void cmdLineWriteValue( Type t , char* str );
template< class Type > void cmdLineCleanUp( Type* t );
template< class Type > Type cmdLineInitialize( void );
template< class Type > Type cmdLineCopy( Type t );
template< class Type > Type cmdLineStringToType( const char* str );
template< class Type >
class cmdLineParameter : public cmdLineReadable
{
public:
Type value;
cmdLineParameter( const char *name );
cmdLineParameter( const char *name , Type v );
~cmdLineParameter( void );
int read( char** argv , int argc );
void writeValue( char* str ) const;
bool expectsArg( void ) const { return true; }
};
template< class Type , int Dim >
class cmdLineParameterArray : public cmdLineReadable
{
public:
Type values[Dim];
cmdLineParameterArray( const char *name, const Type* v=NULL );
~cmdLineParameterArray( void );
int read( char** argv , int argc );
void writeValue( char* str ) const;
bool expectsArg( void ) const { return true; }
};
template< class Type >
class cmdLineParameters : public cmdLineReadable
{
public:
int count;
Type *values;
cmdLineParameters( const char* name );
~cmdLineParameters( void );
int read( char** argv , int argc );
void writeValue( char* str ) const;
bool expectsArg( void ) const { return true; }
};
void cmdLineParse( int argc , char **argv, cmdLineReadable** params );
char* FileExtension( char* fileName );
char* LocalFileName( char* fileName );
char* DirectoryName( char* fileName );
char* GetFileExtension( const char* fileName );
char* GetLocalFileName( const char* fileName );
char** ReadWords( const char* fileName , int& cnt );
#include "CmdLineParser.inl"
#endif // CMD_LINE_PARSER_INCLUDED

Wyświetl plik

@ -0,0 +1,300 @@
/* -*- C++ -*-
Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer. Redistributions in binary form must reproduce
the above copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the distribution.
Neither the name of the Johns Hopkins University nor the names of its contributors
may be used to endorse or promote products derived from this software without specific
prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
*/
#include <cassert>
#include <string.h>
#if defined( WIN32 ) || defined( _WIN64 )
inline int strcasecmp( const char* c1 , const char* c2 ){ return _stricmp( c1 , c2 ); }
#endif // WIN32 || _WIN64
template< > void cmdLineCleanUp< int >( int* t ){ *t = 0; }
template< > void cmdLineCleanUp< float >( float* t ){ *t = 0; }
template< > void cmdLineCleanUp< double >( double* t ){ *t = 0; }
template< > void cmdLineCleanUp< char* >( char** t ){ if( *t ) free( *t ) ; *t = NULL; }
template< > int cmdLineInitialize< int >( void ){ return 0; }
template< > float cmdLineInitialize< float >( void ){ return 0.f; }
template< > double cmdLineInitialize< double >( void ){ return 0.; }
template< > char* cmdLineInitialize< char* >( void ){ return NULL; }
template< > void cmdLineWriteValue< int >( int t , char* str ){ sprintf( str , "%d" , t ); }
template< > void cmdLineWriteValue< float >( float t , char* str ){ sprintf( str , "%f" , t ); }
template< > void cmdLineWriteValue< double >( double t , char* str ){ sprintf( str , "%f" , t ); }
template< > void cmdLineWriteValue< char* >( char* t , char* str ){ if( t ) sprintf( str , "%s" , t ) ; else str[0]=0; }
template< > int cmdLineCopy( int t ){ return t; }
template< > float cmdLineCopy( float t ){ return t; }
template< > double cmdLineCopy( double t ){ return t; }
#if defined( WIN32 ) || defined( _WIN64 )
template< > char* cmdLineCopy( char* t ){ return _strdup( t ); }
#else // !WIN32 && !_WIN64
template< > char* cmdLineCopy( char* t ){ return strdup( t ); }
#endif // WIN32 || _WIN64
template< > int cmdLineStringToType( const char* str ){ return atoi( str ); }
template< > float cmdLineStringToType( const char* str ){ return float( atof( str ) ); }
template< > double cmdLineStringToType( const char* str ){ return double( atof( str ) ); }
#if defined( WIN32 ) || defined( _WIN64 )
template< > char* cmdLineStringToType( const char* str ){ return _strdup( str ); }
#else // !WIN32 && !_WIN64
template< > char* cmdLineStringToType( const char* str ){ return strdup( str ); }
#endif // WIN32 || _WIN64
/////////////////////
// cmdLineReadable //
/////////////////////
#if defined( WIN32 ) || defined( _WIN64 )
inline cmdLineReadable::cmdLineReadable( const char *name ) : set(false) { this->name = _strdup( name ); }
#else // !WIN32 && !_WIN64
inline cmdLineReadable::cmdLineReadable( const char *name ) : set(false) { this->name = strdup( name ); }
#endif // WIN32 || _WIN64
inline cmdLineReadable::~cmdLineReadable( void ){ if( name ) free( name ) ; name = NULL; }
inline int cmdLineReadable::read( char** , int ){ set = true ; return 0; }
inline void cmdLineReadable::writeValue( char* str ) const { str[0] = 0; }
//////////////////////
// cmdLineParameter //
//////////////////////
template< class Type > cmdLineParameter< Type >::~cmdLineParameter( void ) { cmdLineCleanUp( &value ); }
template< class Type > cmdLineParameter< Type >::cmdLineParameter( const char *name ) : cmdLineReadable( name ){ value = cmdLineInitialize< Type >(); }
template< class Type > cmdLineParameter< Type >::cmdLineParameter( const char *name , Type v ) : cmdLineReadable( name ){ value = cmdLineCopy< Type >( v ); }
template< class Type >
int cmdLineParameter< Type >::read( char** argv , int argc )
{
if( argc>0 )
{
cmdLineCleanUp< Type >( &value ) , value = cmdLineStringToType< Type >( argv[0] );
set = true;
return 1;
}
else return 0;
}
template< class Type >
void cmdLineParameter< Type >::writeValue( char* str ) const { cmdLineWriteValue< Type >( value , str ); }
///////////////////////////
// cmdLineParameterArray //
///////////////////////////
template< class Type , int Dim >
cmdLineParameterArray< Type , Dim >::cmdLineParameterArray( const char *name , const Type* v ) : cmdLineReadable( name )
{
if( v ) for( int i=0 ; i<Dim ; i++ ) values[i] = cmdLineCopy< Type >( v[i] );
else for( int i=0 ; i<Dim ; i++ ) values[i] = cmdLineInitialize< Type >();
}
template< class Type , int Dim >
cmdLineParameterArray< Type , Dim >::~cmdLineParameterArray( void ){ for( int i=0 ; i<Dim ; i++ ) cmdLineCleanUp< Type >( values+i ); }
template< class Type , int Dim >
int cmdLineParameterArray< Type , Dim >::read( char** argv , int argc )
{
if( argc>=Dim )
{
for( int i=0 ; i<Dim ; i++ ) cmdLineCleanUp< Type >( values+i ) , values[i] = cmdLineStringToType< Type >( argv[i] );
set = true;
return Dim;
}
else return 0;
}
template< class Type , int Dim >
void cmdLineParameterArray< Type , Dim >::writeValue( char* str ) const
{
char* temp=str;
for( int i=0 ; i<Dim ; i++ )
{
cmdLineWriteValue< Type >( values[i] , temp );
temp = str+strlen( str );
}
}
///////////////////////
// cmdLineParameters //
///////////////////////
template< class Type >
cmdLineParameters< Type >::cmdLineParameters( const char* name ) : cmdLineReadable( name ) , values(NULL) , count(0) { }
template< class Type >
cmdLineParameters< Type >::~cmdLineParameters( void )
{
if( values ) delete[] values;
values = NULL;
count = 0;
}
template< class Type >
int cmdLineParameters< Type >::read( char** argv , int argc )
{
if( values ) delete[] values;
values = NULL;
if( argc>0 )
{
count = atoi(argv[0]);
if( count <= 0 || argc <= count ) return 1;
values = new Type[count];
if( !values ) return 0;
for( int i=0 ; i<count ; i++ ) values[i] = cmdLineStringToType< Type >( argv[i+1] );
set = true;
return count+1;
}
else return 0;
}
template< class Type >
void cmdLineParameters< Type >::writeValue( char* str ) const
{
char* temp=str;
for( int i=0 ; i<count ; i++ )
{
cmdLineWriteValue< Type >( values[i] , temp );
temp = str+strlen( str );
}
}
inline char* FileExtension( char* fileName )
{
char* temp = fileName;
for( unsigned int i=0 ; i<strlen(fileName) ; i++ ) if( fileName[i]=='.' ) temp = &fileName[i+1];
return temp;
}
inline char* GetFileExtension( const char* fileName )
{
char* fileNameCopy;
char* ext=NULL;
char* temp;
fileNameCopy=new char[strlen(fileName)+1];
assert(fileNameCopy);
strcpy(fileNameCopy,fileName);
temp=strtok(fileNameCopy,".");
while(temp!=NULL)
{
if(ext!=NULL){delete[] ext;}
ext=new char[strlen(temp)+1];
assert(ext);
strcpy(ext,temp);
temp=strtok(NULL,".");
}
delete[] fileNameCopy;
return ext;
}
inline char* GetLocalFileName( const char* fileName )
{
char* fileNameCopy;
char* name=NULL;
char* temp;
fileNameCopy=new char[strlen(fileName)+1];
assert(fileNameCopy);
strcpy(fileNameCopy,fileName);
temp=strtok(fileNameCopy,"\\");
while(temp!=NULL){
if(name!=NULL){delete[] name;}
name=new char[strlen(temp)+1];
assert(name);
strcpy(name,temp);
temp=strtok(NULL,"\\");
}
delete[] fileNameCopy;
return name;
}
inline char* LocalFileName( char* fileName )
{
char* temp = fileName;
for( int i=0 ; i<(int)strlen(fileName) ; i++ ) if( fileName[i] =='\\' ) temp = &fileName[i+1];
return temp;
}
inline char* DirectoryName( char* fileName )
{
for( int i=int( strlen(fileName) )-1 ; i>=0 ; i-- )
if( fileName[i] =='\\' )
{
fileName[i] = 0;
return fileName;
}
fileName[0] = 0;
return fileName;
}
inline void cmdLineParse( int argc , char **argv , cmdLineReadable** params )
{
while( argc>0 )
{
if( argv[0][0]=='-' )
{
cmdLineReadable* readable=NULL;
for( int i=0 ; params[i]!=NULL && readable==NULL ; i++ ) if( !strcasecmp( params[i]->name , argv[0]+1 ) ) readable = params[i];
if( readable )
{
int j = readable->read( argv+1 , argc-1 );
argv += j , argc -= j;
}
else
{
fprintf( stderr , "[WARNING] Invalid option: %s\n" , argv[0] );
for( int i=0 ; params[i]!=NULL ; i++ ) printf( "\t-%s\n" , params[i]->name );
}
}
else fprintf( stderr , "[WARNING] Parameter name should be of the form -<name>: %s\n" , argv[0] );
++argv , --argc;
}
}
inline char** ReadWords(const char* fileName,int& cnt)
{
char** names;
char temp[500];
FILE* fp;
fp=fopen(fileName,"r");
if(!fp){return NULL;}
cnt=0;
while(fscanf(fp," %s ",temp)==1){cnt++;}
fclose(fp);
names=new char*[cnt];
if(!names){return NULL;}
fp=fopen(fileName,"r");
if(!fp){
delete[] names;
cnt=0;
return NULL;
}
cnt=0;
while(fscanf(fp," %s ",temp)==1){
names[cnt]=new char[strlen(temp)+1];
if(!names){
for(int j=0;j<cnt;j++){delete[] names[j];}
delete[] names;
cnt=0;
fclose(fp);
return NULL;
}
strcpy(names[cnt],temp);
cnt++;
}
fclose(fp);
return names;
}

Wyświetl plik

@ -0,0 +1,33 @@
#include <cstdio>
#include <cstdarg>
#include "CmdLineParser.h"
struct Logger{
bool verbose;
const char* outputFile;
Logger(){
this->verbose = false;
this->outputFile = NULL;
}
void operator() ( const char* format , ... )
{
if( outputFile )
{
FILE* fp = fopen( outputFile , "a" );
va_list args;
va_start( args , format );
vfprintf( fp , format , args );
fclose( fp );
va_end( args );
}
if( verbose )
{
va_list args;
va_start( args , format );
vprintf( format , args );
va_end( args );
}
}
};

Wyświetl plik

@ -0,0 +1,234 @@
#include <iostream>
#include <string>
#include <fstream>
#include "CmdLineParser.h"
#include "Logger.h"
#include "gdal_priv.h"
#include "cpl_conv.h" // for CPLMalloc()
Logger logWriter;
typedef struct Point2D{
double x;
double y;
Point2D(): x(0.0), y(0.0){}
Point2D(double x, double y): x(x), y(y){}
} Point2D;
typedef struct BoundingBox{
Point2D max;
Point2D min;
BoundingBox(): max(Point2D()), min(Point2D()){}
BoundingBox(Point2D min, Point2D max): max(max), min(min){}
} BoundingBox;
struct PlyPoint{
float x;
float y;
float z;
float nx;
float ny;
float nz;
} p;
size_t psize = sizeof(float) * 6;
float *rasterData;
int arr_width, arr_height;
#define IS_BIG_ENDIAN (*(uint16_t *)"\0\xff" < 0x100)
cmdLineParameter< char* >
InputFile( "inputFile" ) ,
OutputFile( "outputFile" );
cmdLineParameter< float >
SkirtHeightThreshold( "skirtHeightThreshold" ),
SkirtIncrements("skirtIncrements"),
SkirtHeightCap( "skirtHeightCap");
cmdLineReadable
Verbose( "verbose" );
cmdLineReadable* params[] = {
&InputFile , &OutputFile , &SkirtHeightThreshold, &SkirtIncrements, &SkirtHeightCap, &Verbose ,
NULL
};
void help(char *ex){
std::cout << "Usage: " << ex << std::endl
<< "\t -" << InputFile.name << " <input DSM raster>" << std::endl
<< "\t -" << OutputFile.name << " <output PLY points>" << std::endl
<< "\t [-" << SkirtHeightThreshold.name << " <Height threshold between cells that triggers the creation of a skirt>]" << std::endl
<< "\t [-" << SkirtIncrements.name << " <Skirt height increments when adding a new skirt>]" << std::endl
<< "\t [-" << SkirtHeightCap.name << " <Height cap that blocks the creation of a skirt>]" << std::endl
<< "\t [-" << Verbose.name << "]" << std::endl;
exit(EXIT_FAILURE);
}
void logArgs(cmdLineReadable* params[], Logger& logWriter){
logWriter("Running with parameters:\n");
char str[1024];
for( int i=0 ; params[i] ; i++ ){
if( params[i]->set ){
params[i]->writeValue( str );
if( strlen( str ) ) logWriter( "\t--%s %s\n" , params[i]->name , str );
else logWriter( "\t--%s\n" , params[i]->name );
}
}
}
Point2D geoLoc(float xloc, float yloc, double *affine){
return Point2D(affine[0] + xloc*affine[1] + yloc*affine[2], affine[3] + xloc*affine[4] + yloc*affine[5]);
}
BoundingBox getExtent(GDALDataset *dataset){
double affine[6];
dataset->GetGeoTransform(affine);
return BoundingBox(geoLoc(0, dataset->GetRasterYSize(), affine), geoLoc(dataset->GetRasterXSize(), 0, affine));
}
int countSkirts(int nx, int ny){
float neighbor_z = rasterData[ny * arr_width + nx];
float current_z = p.z;
int result = 0;
if (current_z - neighbor_z > SkirtHeightThreshold.value && current_z - neighbor_z < SkirtHeightCap.value){
while (current_z > neighbor_z){
current_z -= SkirtIncrements.value;
result++;
}
result++;
}
return result;
}
void writeSkirts(std::ofstream &f, int nx, int ny){
float neighbor_z = rasterData[ny * arr_width + nx];
float z_val = p.z;
if (p.z - neighbor_z > SkirtHeightThreshold.value && p.z - neighbor_z < SkirtHeightCap.value){
while (p.z > neighbor_z){
p.z -= SkirtIncrements.value;
f.write(reinterpret_cast<char*>(&p), psize);
}
p.z = neighbor_z;
f.write(reinterpret_cast<char*>(&p), psize);
}
// Restore original p.z value
p.z = z_val;
}
int main(int argc, char **argv) {
cmdLineParse( argc-1 , &argv[1] , params );
if( !InputFile.set || !OutputFile.set ) help(argv[0]);
if (!SkirtHeightThreshold.set) SkirtHeightThreshold.value = 1.5f;
if (!SkirtIncrements.set) SkirtIncrements.value = 0.1f;
if (!SkirtHeightCap.set) SkirtHeightCap.value = 100.0f;
logWriter.verbose = Verbose.set;
// logWriter.outputFile = "odm_dem2points_log.txt";
logArgs(params, logWriter);
GDALDataset *dataset;
GDALAllRegister();
dataset = (GDALDataset *) GDALOpen( InputFile.value, GA_ReadOnly );
if( dataset != NULL )
{
arr_width = dataset->GetRasterXSize();
arr_height = dataset->GetRasterYSize();
logWriter("Raster Size is %dx%d\n", arr_width, arr_height);
BoundingBox extent = getExtent(dataset);
logWriter("Extent is (%f, %f), (%f, %f)\n", extent.min.x, extent.max.x, extent.min.y, extent.max.y);
float ext_width = extent.max.x - extent.min.x;
float ext_height = extent.max.y - extent.min.y;
int vertex_count = (arr_height - 4) * (arr_width - 4);
int skirt_points = 0;
GDALRasterBand *band = dataset->GetRasterBand(1);
rasterData = new float[arr_width * arr_height];
if (band->RasterIO( GF_Read, 0, 0, arr_width, arr_height,
rasterData, arr_width, arr_height, GDT_Float32, 0, 0 ) == CE_Failure){
std::cerr << "Cannot access raster data" << std::endl;
exit(1);
}
// On the first pass we calculate the number of skirts
// on the second pass we sample and write the points
logWriter("Calculating skirts... ");
for (int y = 2; y < arr_height - 2; y++){
for (int x = 2; x < arr_width - 2; x++){
p.z = rasterData[y * arr_width + x];
skirt_points += countSkirts(x, y + 1);
skirt_points += countSkirts(x, y - 1);
skirt_points += countSkirts(x + 1, y);
skirt_points += countSkirts(x - 1, y);
}
}
logWriter("%d vertices will be added\n", skirt_points);
logWriter("Total vertices: %d\n", (skirt_points + vertex_count));
logWriter("Sampling and writing to file...");
// Starting writing ply file
std::ofstream f (OutputFile.value);
f << "ply" << std::endl;
if (IS_BIG_ENDIAN){
f << "format binary_big_endian 1.0" << std::endl;
}else{
f << "format binary_little_endian 1.0" << std::endl;
}
f << "element vertex " << (vertex_count + skirt_points) << std::endl
<< "property float x" << std::endl
<< "property float y" << std::endl
<< "property float z" << std::endl
<< "property float nx" << std::endl
<< "property float ny" << std::endl
<< "property float nz" << std::endl
<< "end_header" << std::endl;
// Normals are always pointing up
p.nx = 0;
p.ny = 0;
p.nz = 1;
for (int y = 2; y < arr_height - 2; y++){
for (int x = 2; x < arr_width - 2; x++){
p.z = rasterData[y * arr_width + x];
p.x = extent.min.x + (static_cast<float>(x) / static_cast<float>(arr_width)) * ext_width;
p.y = extent.max.y - (static_cast<float>(y) / static_cast<float>(arr_height)) * ext_height;
f.write(reinterpret_cast<char*>(&p), psize);
writeSkirts(f, x, y + 1);
writeSkirts(f, x, y - 1);
writeSkirts(f, x + 1, y);
writeSkirts(f, x - 1, y);
}
}
logWriter(" done!\n");
f.close();
GDALClose(dataset);
}else{
std::cerr << "Cannot open " << InputFile.value << std::endl;
}
}

Wyświetl plik

@ -1,6 +1,5 @@
from __future__ import absolute_import
import os, shutil, sys, struct, random
from io import BytesIO
from gippy import GeoImage
from opendm.dem import commands
from opendm import system
@ -36,7 +35,7 @@ def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples=
max_workers=max_workers
)
dsm_points = dem_to_points(os.path.join(tmp_directory, 'mesh_dsm.tif'), os.path.join(tmp_directory, 'dsm_points.ply'))
dsm_points = dem_to_points(os.path.join(tmp_directory, 'mesh_dsm.tif'), os.path.join(tmp_directory, 'dsm_points.ply'), verbose)
mesh = screened_poisson_reconstruction(dsm_points, outMesh, depth=depth,
samples=samples,
maxVertexCount=maxVertexCount,
@ -49,68 +48,87 @@ def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples=
return mesh
def dem_to_points(inGeotiff, outPointCloud):
def dem_to_points(inGeotiff, outPointCloud, verbose=False):
log.ODM_INFO('Sampling points from DSM: %s' % inGeotiff)
image = GeoImage.open([inGeotiff], bandnames=['z'], nodata=-9999)
arr = image['z'].read_raw()
resolution = max(abs(image.resolution().x()), abs(image.resolution().y()))
kwargs = {
'bin': context.odm_modules_path,
'outfile': outPointCloud,
'infile': inGeotiff,
'verbose': '-verbose' if verbose else ''
}
log.ODM_INFO('Writing points...')
mem_file = BytesIO()
xmin, xmax, ymin, ymax = image.extent().x0(), image.extent().x1(), image.extent().y0(), image.extent().y1()
ext_width, ext_height = xmax - xmin, ymax - ymin
arr_height, arr_width = arr.shape
vertex_count = (arr_height - 4) * (arr_width - 4)
skirt_points = 0
skirt_height_threshold = 1 # meter
skirt_increments = 0.1
for x in range(2, arr_width - 2):
for y in range(2, arr_height - 2):
z = arr[y][x]
tx = xmin + (float(x) / float(arr_width)) * ext_width
ty = ymax - (float(y) / float(arr_height)) * ext_height
mem_file.write(struct.pack('ffffff', tx, ty, z, 0, 0, 1))
# Skirting
for (nx, ny) in ((x, y + 1), (x, y - 1), (x + 1, y), (x - 1, y)):
current_z = z
neighbor_z = arr[ny][nx]
if current_z - neighbor_z > skirt_height_threshold:
while current_z > neighbor_z:
current_z -= skirt_increments
mem_file.write(struct.pack('ffffff', tx, ty, current_z, 0, 0, 1))
skirt_points += 1
mem_file.write(struct.pack('ffffff', tx, ty, neighbor_z, 0, 0, 1))
skirt_points += 1
with open(outPointCloud, "wb") as f:
f.write("ply\n")
f.write("format binary_%s_endian 1.0\n" % sys.byteorder)
f.write("element vertex %s\n" % (vertex_count + skirt_points))
f.write("property float x\n")
f.write("property float y\n")
f.write("property float z\n")
f.write("property float nx\n")
f.write("property float ny\n")
f.write("property float nz\n")
f.write("end_header\n")
f.write(mem_file.getvalue())
mem_file.close()
log.ODM_INFO("Points count: %s (%s samples, %s skirts)", vertex_count + skirt_points, vertex_count, skirt_points)
log.ODM_INFO('Wrote points to: %s' % outPointCloud)
system.run('{bin}/odm_dem2points -inputFile {infile} '
'-outputFile {outfile} '
'-skirtHeightThreshold 1.5 '
'-skirtIncrements 0.2 '
'-skirtHeightCap 100 '
' {verbose} '.format(**kwargs))
return outPointCloud
# Old Python implementation of dem_to_points
# def dem_to_points(inGeotiff, outPointCloud):
# log.ODM_INFO('Sampling points from DSM: %s' % inGeotiff)
# image = GeoImage.open([inGeotiff], bandnames=['z'], nodata=-9999)
# arr = image['z'].read_raw()
# mem_file = BytesIO()
# xmin, xmax, ymin, ymax = image.extent().x0(), image.extent().x1(), image.extent().y0(), image.extent().y1()
# ext_width, ext_height = xmax - xmin, ymax - ymin
# arr_height, arr_width = arr.shape
# vertex_count = (arr_height - 4) * (arr_width - 4)
# skirt_points = 0
# skirt_height_threshold = 1 # meter
# skirt_increments = 0.1
# for x in range(2, arr_width - 2):
# for y in range(2, arr_height - 2):
# z = arr[y][x]
# tx = xmin + (float(x) / float(arr_width)) * ext_width
# ty = ymax - (float(y) / float(arr_height)) * ext_height
# mem_file.write(struct.pack('ffffff', tx, ty, z, 0, 0, 1))
# # Skirting
# for (nx, ny) in ((x, y + 1), (x, y - 1), (x + 1, y), (x - 1, y)):
# current_z = z
# neighbor_z = arr[ny][nx]
# if current_z - neighbor_z > skirt_height_threshold:
# while current_z > neighbor_z:
# current_z -= skirt_increments
# mem_file.write(struct.pack('ffffff', tx, ty, current_z, 0, 0, 1))
# skirt_points += 1
# mem_file.write(struct.pack('ffffff', tx, ty, neighbor_z, 0, 0, 1))
# skirt_points += 1
# log.ODM_INFO("Points count: %s (%s samples, %s skirts)", vertex_count + skirt_points, vertex_count, skirt_points)
# log.ODM_INFO('Writing points...')
# mem_file.seek(0)
# with open(outPointCloud, "wb") as f:
# f.write("ply\n")
# f.write("format binary_%s_endian 1.0\n" % sys.byteorder)
# f.write("element vertex %s\n" % (vertex_count + skirt_points))
# f.write("property float x\n")
# f.write("property float y\n")
# f.write("property float z\n")
# f.write("property float nx\n")
# f.write("property float ny\n")
# f.write("property float nz\n")
# f.write("end_header\n")
# shutil.copyfileobj(mem_file, f)
# mem_file.close()
# log.ODM_INFO('Wrote points to: %s' % outPointCloud)
# return outPointCloud
def screened_poisson_reconstruction(inPointCloud, outMesh, depth = 8, samples = 1, maxVertexCount=100000, pointWeight=4, threads=context.num_cores, verbose=False):