kopia lustrzana https://github.com/OpenDroneMap/ODM
rodzic
6390a2a7fd
commit
8958f4ff86
|
@ -5,17 +5,15 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR})
|
|||
|
||||
set (CMAKE_CXX_STANDARD 11)
|
||||
find_package(GDAL REQUIRED)
|
||||
find_package(VTK REQUIRED)
|
||||
include_directories(${GDAL_INCLUDE_DIR})
|
||||
include(${VTK_USE_FILE})
|
||||
|
||||
# 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} ${VTK_LIBRARIES})
|
||||
|
||||
target_link_libraries(${PROJECT_NAME} ${GDAL_LIBRARY})
|
||||
|
|
Plik diff jest za duży
Load Diff
|
@ -4,12 +4,7 @@
|
|||
#include <cmath>
|
||||
#include "CmdLineParser.h"
|
||||
#include "Logger.h"
|
||||
#include <vtkSmartPointer.h>
|
||||
#include <vtkGreedyTerrainDecimation.h>
|
||||
#include <vtkTransform.h>
|
||||
#include <vtkTransformFilter.h>
|
||||
#include <vtkPLYWriter.h>
|
||||
#include <vtkImageData.h>
|
||||
#include "Simplify.h"
|
||||
|
||||
#include "gdal_priv.h"
|
||||
#include "cpl_conv.h" // for CPLMalloc()
|
||||
|
@ -32,13 +27,32 @@ typedef struct BoundingBox{
|
|||
BoundingBox(Point2D min, Point2D max): max(max), min(min){}
|
||||
} BoundingBox;
|
||||
|
||||
struct PlyPoint{
|
||||
float x;
|
||||
float y;
|
||||
float z;
|
||||
} p;
|
||||
size_t psize = sizeof(float) * 3;
|
||||
|
||||
struct PlyFace{
|
||||
uint32_t p1;
|
||||
uint32_t p2;
|
||||
uint32_t p3;
|
||||
} face;
|
||||
size_t fsize = sizeof(uint32_t) * 3;
|
||||
|
||||
float *rasterData;
|
||||
int arr_width, arr_height;
|
||||
|
||||
#define IS_BIG_ENDIAN (*(uint16_t *)"\0\xff" < 0x100)
|
||||
|
||||
cmdLineParameter< char* >
|
||||
InputFile( "inputFile" ) ,
|
||||
OutputFile( "outputFile" );
|
||||
cmdLineParameter< int >
|
||||
MaxVertexCount( "maxVertexCount" );
|
||||
cmdLineReadable
|
||||
Verbose( "verbose" );
|
||||
Verbose( "verbose" );
|
||||
|
||||
cmdLineReadable* params[] = {
|
||||
&InputFile , &OutputFile , &MaxVertexCount , &Verbose ,
|
||||
|
@ -85,74 +99,147 @@ int main(int argc, char **argv) {
|
|||
logWriter.verbose = Verbose.set;
|
||||
logArgs(params, logWriter);
|
||||
|
||||
|
||||
GDALDataset *dataset;
|
||||
GDALAllRegister();
|
||||
dataset = (GDALDataset *) GDALOpen( InputFile.value, GA_ReadOnly );
|
||||
if( dataset != NULL )
|
||||
{
|
||||
int width = dataset->GetRasterXSize();
|
||||
int height = dataset->GetRasterYSize();
|
||||
arr_width = dataset->GetRasterXSize();
|
||||
arr_height = dataset->GetRasterYSize();
|
||||
|
||||
logWriter("Raster Size is %dx%d\n", width, height);
|
||||
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);
|
||||
|
||||
double extentX = extent.max.x - extent.min.x;
|
||||
double extentY = extent.max.y - extent.min.y;
|
||||
float ext_width = extent.max.x - extent.min.x;
|
||||
float ext_height = extent.max.y - extent.min.y;
|
||||
|
||||
vtkSmartPointer<vtkImageData> image =
|
||||
vtkSmartPointer<vtkImageData>::New();
|
||||
image->SetDimensions(width, height, 1);
|
||||
image->AllocateScalars(VTK_FLOAT, 1);
|
||||
unsigned long long int vertex_count = static_cast<unsigned long long int>(arr_height) *
|
||||
static_cast<unsigned long long int>(arr_width);
|
||||
|
||||
// TODO:
|
||||
// 2. Figure out bad_alloc in cmparks dataset
|
||||
|
||||
logWriter("Reading raster...\n");
|
||||
|
||||
// If the DSM is really large, we only sample a subset of it
|
||||
// to remain within INT_MAX vertices. This does not happen often,
|
||||
// but it's a safeguard to make sure we'll get an output and not
|
||||
// overflow.
|
||||
int stride = 1;
|
||||
while (vertex_count > INT_MAX){
|
||||
stride += 1;
|
||||
vertex_count = static_cast<int>(std::ceil((arr_height / static_cast<double>(stride))) *
|
||||
std::ceil((arr_width / static_cast<double>(stride))));
|
||||
}
|
||||
|
||||
if (stride != 1){
|
||||
logWriter("Warning: DSM is large, will sample using stride value of %d\n", stride);
|
||||
}
|
||||
logWriter("Total vertices before simplification: %llu\n", vertex_count);
|
||||
|
||||
GDALRasterBand *band = dataset->GetRasterBand(1);
|
||||
float *rasterData = new float[width];
|
||||
rasterData = new float[arr_width];
|
||||
|
||||
for (int y = 0; y < height; y++){
|
||||
if (band->RasterIO( GF_Read, 0, height - y - 1, width, 1,
|
||||
rasterData, width, 1, GDT_Float32, 0, 0 ) == CE_Failure){
|
||||
for (int y = 0; y < arr_height; y += stride){
|
||||
if (band->RasterIO( GF_Read, 0, y, arr_width, 1,
|
||||
rasterData, arr_width, 1, GDT_Float32, 0, 0 ) == CE_Failure){
|
||||
std::cerr << "Cannot access raster data" << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
for (int x = 0; x < width; x++){
|
||||
(static_cast<float*>(image->GetScalarPointer(x,y,0)))[0] = rasterData[x];
|
||||
for (int x = 0; x < arr_width; x += stride){
|
||||
Simplify::Vertex v;
|
||||
v.p.x = extent.min.x + (static_cast<float>(x) / static_cast<float>(arr_width)) * ext_width;
|
||||
v.p.y = extent.max.y - (static_cast<float>(y) / static_cast<float>(arr_height)) * ext_height;
|
||||
v.p.z = rasterData[x];
|
||||
|
||||
Simplify::vertices.push_back(v);
|
||||
}
|
||||
}
|
||||
|
||||
GDALClose(dataset);
|
||||
delete[] rasterData;
|
||||
GDALClose(dataset);
|
||||
|
||||
vtkSmartPointer<vtkGreedyTerrainDecimation> terrain =
|
||||
vtkSmartPointer<vtkGreedyTerrainDecimation>::New();
|
||||
terrain->SetErrorMeasureToNumberOfTriangles();
|
||||
terrain->SetNumberOfTriangles(MaxVertexCount.value * 2); // Approximate
|
||||
terrain->SetInputData(image);
|
||||
terrain->BoundaryVertexDeletionOn();
|
||||
unsigned int cols = static_cast<unsigned int>(std::ceil((arr_width / static_cast<double>(stride))));
|
||||
unsigned int rows = static_cast<unsigned int>(std::ceil((arr_height / static_cast<double>(stride))));
|
||||
|
||||
logWriter("OK\nTransform... ");
|
||||
vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
|
||||
transform->Scale(extentX / width, extentY / height, 1);
|
||||
transform->Translate(extent.min.x * (width / extentX), extent.min.y * (height / extentY), 0);
|
||||
for (unsigned int y = 0; y < rows - 1; y++){
|
||||
for (unsigned int x = 0; x < cols - 1; x++){
|
||||
Simplify::Triangle t1;
|
||||
t1.v[0] = cols * (y + 1) + x;
|
||||
t1.v[1] = cols * y + x + 1;
|
||||
t1.v[2] = cols * y + x;
|
||||
t1.attr = 0;
|
||||
t1.material = -1;
|
||||
|
||||
vtkSmartPointer<vtkTransformFilter> transformFilter =
|
||||
vtkSmartPointer<vtkTransformFilter>::New();
|
||||
transformFilter->SetInputConnection(terrain->GetOutputPort());
|
||||
transformFilter->SetTransform(transform);
|
||||
Simplify::triangles.push_back(t1);
|
||||
|
||||
logWriter("OK\nSaving mesh to file... ");
|
||||
Simplify::Triangle t2;
|
||||
t2.v[0] = cols * (y + 1) + x;
|
||||
t2.v[1] = cols * (y + 1) + x + 1;
|
||||
t2.v[2] = cols * y + x + 1;
|
||||
t2.attr = 0;
|
||||
t2.material = -1;
|
||||
|
||||
vtkSmartPointer<vtkPLYWriter> plyWriter =
|
||||
vtkSmartPointer<vtkPLYWriter>::New();
|
||||
plyWriter->SetFileName(OutputFile.value);
|
||||
plyWriter->SetInputConnection(transformFilter->GetOutputPort());
|
||||
plyWriter->SetFileTypeToBinary();
|
||||
plyWriter->Write();
|
||||
Simplify::triangles.push_back(t2);
|
||||
}
|
||||
}
|
||||
|
||||
logWriter("OK\n");
|
||||
double agressiveness = 7.0;
|
||||
int target_count = std::min(MaxVertexCount.value * 2, static_cast<int>(Simplify::triangles.size()));
|
||||
|
||||
logWriter("Sampled %d faces, target is %d\n", static_cast<int>(Simplify::triangles.size()), target_count);
|
||||
logWriter("Simplifying...\n");
|
||||
|
||||
unsigned long start_size = Simplify::triangles.size();
|
||||
Simplify::simplify_mesh(target_count, agressiveness, true);
|
||||
if ( Simplify::triangles.size() >= start_size) {
|
||||
std::cerr << "Unable to reduce mesh.\n";
|
||||
exit(1);
|
||||
}
|
||||
|
||||
logWriter("Writing to file...");
|
||||
|
||||
// Start 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 " << Simplify::vertices.size() << std::endl
|
||||
<< "property float x" << std::endl
|
||||
<< "property float y" << std::endl
|
||||
<< "property float z" << std::endl
|
||||
<< "element face " << Simplify::triangles.size() << std::endl
|
||||
<< "property list uint8 uint32 vertex_indices" << std::endl
|
||||
<< "end_header" << std::endl;
|
||||
|
||||
for(Simplify::Vertex &v : Simplify::vertices){
|
||||
p.x = static_cast<float>(v.p.x);
|
||||
p.y = static_cast<float>(v.p.y);
|
||||
p.z = static_cast<float>(v.p.z);
|
||||
f.write(reinterpret_cast<char *>(&p), psize);
|
||||
}
|
||||
|
||||
uint8_t three = 3;
|
||||
for(Simplify::Triangle &t : Simplify::triangles){
|
||||
face.p1 = static_cast<uint32_t>(t.v[0]);
|
||||
face.p2 = static_cast<uint32_t>(t.v[1]);
|
||||
face.p3 = static_cast<uint32_t>(t.v[2]);
|
||||
|
||||
f.write(reinterpret_cast<char *>(&three), sizeof(three));
|
||||
f.write(reinterpret_cast<char *>(&face), fsize);
|
||||
}
|
||||
|
||||
logWriter(" done!\n");
|
||||
|
||||
f.close();
|
||||
}else{
|
||||
std::cerr << "Cannot open " << InputFile.value << std::endl;
|
||||
}
|
||||
|
|
Ładowanie…
Reference in New Issue