first commit

alphaLas/WriteOBJ.m

@@ -0,0 +1,26 @@
+% Write simple object files from coordinates and faces indices
+% Copyright ?2017 Ecole polytechnique f?d?rale de Lausanne - EPFL
+% Laboratory of Interdisciplinary Performance in Design - LIPID
+% @license GNU GPLv3 https://www.gnu.org/licenses/gpl-3.0.en.html
+
+
+% Author: Giuseppe Peronato, <[email protected]>
+% Latest modified on 18.07.2017
+
+
+function WriteOBJ(pts,tri,fname,material)
+fid = fopen(fname, 'wt');
+if size(material) > 0    
+    fprintf(fid, 'mtllib ./material.mtl\n');
+end
+for i = 1:size(pts,1)
+   fprintf(fid, 'v %12.4f %12.4f %4.4f\n', pts(i,:));
+end
+fprintf(fid, 'o Object\n');
+if size(material) > 0    
+    fprintf(fid, 'usemtl %s\n',material);
+end
+for i = 1:size(tri,1)
+   fprintf(fid, 'f %g %g %g\n', tri(i,:));
+end
+fclose(fid);

alphaLas/alphaLAS.m

@@ -0,0 +1,82 @@
+% Create a mesh from a LAS dataset using ?-shapes
+% Copyright ?2017 Ecole polytechnique f?d?rale de Lausanne - EPFL
+% Laboratory of Interdisciplinary Performance in Design - LIPID
+% @license GNU GPLv3 https://www.gnu.org/licenses/gpl-3.0.en.html
+
+% Author: Giuseppe Peronato, <[email protected]>
+% Latest modified on 19.07.2017
+
+% Based on a script by Matthew Parkan - EPFL LASIG
+% Original concept implemented in solar energy assessments:
+% "3D-modeling of vegetation from LiDAR point clouds and assessment of its impact on fa?ade solar irradiation"
+% http://dx.doi.org/10.5194/isprs-archives-XLII-2-W2-67-2016
+
+% Dependencies:
+% - Digital Forestry Toolbox (DFT) by Matthew Parkan
+%   https://github.com/mparkan/Digital-Forestry-Toolbox
+% - WriteOBJ function by Giuseppe Peronato
+
+clc
+clear 
+close all
+opengl hardware
+
+% Set parameters
+
+% DFT directory
+addpath(genpath('C:\Users\giupe250\Documents\MATLAB\mparkan-Digital-Forestry-Toolbox-01e26d8\scripts\')) %path to DFT scripts
+
+% IO
+idir = 'C:\Users\giupe250\Data\LiDAR_High\Vegetation_tiles\cleaned_100cm\'; %directory with input *.las files
+odir = 'C:\Users\giupe250\Data\LiDAR_High\Vegetation_tiles_OBJ\cleaned_100cm\'; %directory with output *.obj files
+prefix = 'tile'; %optional prefix to output files
+
+% Alpha-shape https://en.wikipedia.org/wiki/Alpha_shape
+radius = 0.75; %alpha-radius
+holethreshold = 0; % area of the largest hole to fill
+regionthreshold = 1; % minimum volume
+
+% LiDAR
+class = 5; % LiDAR class: class 5  usually corresponds to high vegetation
+
+
+%Script
+files = dir(strcat(idir,'*.las')); % list the *.las files in the input directory
+
+%iterate over each las file
+for i = 1:length(files) %loop over all files
+    filename = strrep(files(i).name,'.las',''); %filename without extension
+    display(sprintf('Converting file %s', filename));
+
+    opt.input.file.points = strcat(idir,filename,'.las'); %input *.las fle
+
+
+    pc = LASread(opt.input.file.points); %pointcloud
+
+    xyz = [pc.record.x, pc.record.y, pc.record.z];
+    classification = pc.record.classification;
+
+    clear pc;
+
+    idxl_veg = ismember(classification, class); %class 5  usually corresponds to high vegetation
+
+    shp = alphaShape(xyz(idxl_veg,:), ...
+        radius, ...
+        'HoleThreshold', holethreshold, ...
+        'RegionThreshold', regionthreshold);
+
+    [tri, pts] = boundaryFacets(shp); %Extract the boundaries of the hulls
+
+    WriteOBJ(pts,tri,strcat(odir,prefix,filename,'.obj'),"Vegetation");
+
+    clear opt;
+    %clear xyz;
+    clear classification;
+    %clear idxl_veg;
+    clear shp;
+    clear pts;
+    clear tri
+
+
+    display(sprintf('Conversion of file %s completed.\n', filename));
+end

alphaXYZ/alphashape.cpp

@@ -0,0 +1,166 @@
+// This script reads an XYZ file and gives as output an OFF file containing the visible alpha shape facets
+// It is based on the visible_alpha_shape_facets_to_OFF example file
+// The compiled executable can be run in command line using as arguments the parameter alpha and the path to the XYZ file
+// Contact: Giuseppe Peronato <[email protected]>
+// Latest modification 12/03/2019
+
+
+#include <chrono>  // for high_resolution_clock
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Alpha_shape_3.h>
+#include <CGAL/Alpha_shape_cell_base_3.h>
+#include <CGAL/Alpha_shape_vertex_base_3.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+
+#include <cassert>
+#include <fstream>
+#include <list>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel               Gt;
+typedef CGAL::Tag_false                                                    Alpha_cmp_tag;
+//We use CGAL::Default to skip the complete declaration of base classes
+typedef CGAL::Alpha_shape_vertex_base_3<Gt, CGAL::Default, Alpha_cmp_tag>   Vb;
+typedef CGAL::Alpha_shape_cell_base_3<Gt, CGAL::Default, Alpha_cmp_tag>     Fb;
+typedef CGAL::Triangulation_data_structure_3<Vb, Fb>                       Tds;
+typedef CGAL::Delaunay_triangulation_3<Gt, Tds, CGAL::Fast_location>       Triangulation_3;
+//Alpha shape with ExactAlphaComparisonTag set to False (note that the tag is also
+//set to false for Vb and Fb)
+typedef CGAL::Alpha_shape_3<Triangulation_3, Alpha_cmp_tag>                Alpha_shape_3;
+typedef Gt::Point_3                                                       Point;
+
+// Record start time
+auto start = std::chrono::high_resolution_clock::now();
+
+int main(int argc, char* argv[])
+{
+	// Initialize parameters
+	double alpha = 0.5;
+	char* infile = "input.xyz";
+
+	// Check the number of parameters
+	if (argc < 2) {
+		// Tell the user how to run the program
+		std::cerr << "Run in command line using as arguments alpha and the path to the XYZ file.\nE.g. alphashape.exe 0.75 infile.xyz\n";
+		system("pause");
+		return 0;
+	}
+	else {
+		alpha = atof(argv[1]);  // alternative strtod
+		infile = argv[2];
+	}
+  std::vector<Point> points;
+
+//read input
+  std::ifstream myfile(infile);
+
+  // new lines will be skipped unless we stop it from happening:    
+  myfile.unsetf(std::ios_base::skipws);
+
+  // count the newlines with an algorithm specialized for counting:
+  int n = std::count(
+	  std::istream_iterator<char>(myfile),
+	  std::istream_iterator<char>(),
+	  '\n');
+
+  std::ifstream is(infile);
+  double x, y, z;
+  for (int i=0;i<n;++i)
+  {
+	  is >> x >> y >> z;
+    points.push_back( Point(x, y, z) );
+  }
+
+  std::cerr << points.size() << " points read!\n";
+
+// set alpha shape
+  std::cerr << "alpha = " << alpha << "\n";
+  Alpha_shape_3 as(points.begin(), points.end(), alpha, Alpha_shape_3::REGULARIZED);
+  std::cerr << as.number_of_solid_components() << " number of solid components\n";
+
+// collect alpha-shape facets accessible from the infinity
+  // marks the cells that are in the same component as the infinite vertex by flooding
+  boost::unordered_set< Alpha_shape_3::Cell_handle > marked_cells;
+  std::vector< Alpha_shape_3::Cell_handle > queue;
+  queue.push_back( as.infinite_cell() );
+
+  while(!queue.empty())
+  {
+    Alpha_shape_3::Cell_handle back = queue.back();
+    queue.pop_back();
+
+    if ( !marked_cells.insert(back).second ) continue; //already visited
+
+    for (int i=0; i<4; ++i)
+    {
+      if (as.classify(Alpha_shape_3::Facet(back, i))==Alpha_shape_3::EXTERIOR &&
+          marked_cells.count(back->neighbor(i))==0)
+        queue.push_back( back->neighbor(i) );
+    }
+  }
+
+  // filter regular facets to restrict them to those adjacent to a marked cell
+  std::vector< Alpha_shape_3::Facet > regular_facets;
+  as.get_alpha_shape_facets(std::back_inserter( regular_facets ), Alpha_shape_3::REGULAR );
+
+  std::vector<Alpha_shape_3::Facet> filtered_regular_facets;
+  BOOST_FOREACH(Alpha_shape_3::Facet f, regular_facets)
+  {
+    if ( marked_cells.count(f.first)==1 )
+      filtered_regular_facets.push_back(f);
+    else
+    {
+      f = as.mirror_facet(f);
+      if ( marked_cells.count(f.first)==1 )
+        filtered_regular_facets.push_back(f);
+    }
+  }
+
+// dump into OFF format
+  // assign an id per vertex
+  boost::unordered_map< Alpha_shape_3::Vertex_handle, std::size_t> vids;
+  points.clear();
+
+  BOOST_FOREACH(Alpha_shape_3::Facet f, filtered_regular_facets)
+  {
+    for (int i=1;i<4; ++i)
+    {
+      Alpha_shape_3::Vertex_handle vh = f.first->vertex((f.second+i)%4);
+      if (vids.insert( std::make_pair(vh, points.size()) ).second)
+        points.push_back( vh->point() );
+    }
+  }
+
+  // writing
+  char* outfile = infile;
+  sscanf(infile, "%[^.]", outfile);  // foo.bar => foo
+  sprintf(outfile, "%s.off", outfile); // foo.txt <= foo
+
+  std::ofstream output(outfile);
+  output.precision(10);
+  output << "OFF\n " << points.size() << " " << filtered_regular_facets.size() << " 0\n";
+  std::copy(points.begin(), points.end(), std::ostream_iterator<Point>(output, "\n"));
+  BOOST_FOREACH(const Alpha_shape_3::Facet& f, filtered_regular_facets)
+  {
+    output << 3;
+
+    for (int i=0;i<3; ++i)
+    {
+      Alpha_shape_3::Vertex_handle vh = f.first->vertex( as.vertex_triple_index(f.second, i) );
+      output << " " << vids[vh];
+    }
+    output << "\n";
+  }
+
+  // Record end time
+  auto finish = std::chrono::high_resolution_clock::now();
+
+
+  std::chrono::duration<double> elapsed = finish - start;
+
+  std::cout << "Elapsed time: " << elapsed.count() << " s\n";
+
+
+  return 0;
+}