1 /****************************************************************************
2 * MeshLab o o *
3 * A versatile mesh processing toolbox o o *
4 * _ O _ *
5 * Copyright(C) 2005 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23
24 #ifdef WIN32
25 #include <windows.h>
26 #include <Psapi.h>
27 #endif
28 #include <cstdio>
29 #include "Src/MyTime.h"
30 #include "Src/MemoryUsage.h"
31 #include "Src/MarchingCubes.h"
32 #include "Src/Octree.h"
33 #include "Src/SparseMatrix.h"
34 #include "Src/CmdLineParser.h"
35 #include "Src/PPolynomial.h"
36 void DumpOutput( const char* format , ... );
37 void DumpOutput2( char* str , const char* format , ... );
38
39 #include "Src/PointStream.h"
40
41 #include "Src/MultiGridOctreeData.h"
42 #include "filter_screened_poisson.h"
43 #include <QtScript>
44
45
46
DumpOutput(const char * format,...)47 void DumpOutput( const char* format , ... )
48 {
49 char buf[4096];
50 va_list marker;
51 va_start( marker, format );
52
53 vsprintf(buf,format,marker);
54 va_end( marker );
55
56 qDebug(buf);
57 }
DumpOutput2(std::vector<char * > & comments,const char * format,...)58 void DumpOutput2(std::vector< char* >& comments , const char* format , ... )
59 {
60 char buf[4096];
61 va_list marker;
62 va_start( marker, format );
63
64 vsprintf(buf,format,marker);
65 va_end( marker );
66 qDebug(buf);
67 }
68
69
70 #if defined( _WIN32 ) || defined( _WIN64 )
PeakMemoryUsageMB(void)71 double PeakMemoryUsageMB( void )
72 {
73 HANDLE h = GetCurrentProcess();
74 PROCESS_MEMORY_COUNTERS pmc;
75 return GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ? ( (double)pmc.PeakWorkingSetSize )/(1<<20) : 0;
76 }
77 #endif // _WIN32 || _WIN64
78
79 #if defined( _WIN32 ) || defined( _WIN64 )
to_seconds(const FILETIME & ft)80 inline double to_seconds( const FILETIME& ft )
81 {
82 const double low_to_sec=100e-9; // 100 nanoseconds
83 const double high_to_sec=low_to_sec*4294967296.0;
84 return ft.dwLowDateTime*low_to_sec+ft.dwHighDateTime*high_to_sec;
85 }
86 #endif // _WIN32 || _WIN64
87
88
89 template< class Real >
90 struct OctreeProfiler
91 {
92 Octree< Real >& tree;
93 double t;
94
OctreeProfilerOctreeProfiler95 OctreeProfiler( Octree< Real >& t ) : tree(t) { ; }
startOctreeProfiler96 void start( void ){ t = Time() , tree.resetLocalMemoryUsage(); }
printOctreeProfiler97 void print( const char* header ) const
98 {
99 tree.memoryUsage();
100 #if defined( _WIN32 ) || defined( _WIN64 )
101 if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
102 else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
103 #else // !_WIN32 && !_WIN64
104 if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
105 else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
106 #endif // _WIN32 || _WIN64
107 }
dumpOutputOctreeProfiler108 void dumpOutput( const char* header ) const
109 {
110 tree.memoryUsage();
111 #if defined( _WIN32 ) || defined( _WIN64 )
112 if( header ) DumpOutput( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
113 else DumpOutput( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
114 #else // !_WIN32 && !_WIN64
115 if( header ) DumpOutput( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
116 else DumpOutput( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
117 #endif // _WIN32 || _WIN64
118 }
dumpOutput2OctreeProfiler119 void dumpOutput2( std::vector< char* >& comments , const char* header ) const
120 {
121 tree.memoryUsage();
122 #if defined( _WIN32 ) || defined( _WIN64 )
123 if( header ) DumpOutput2( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
124 else DumpOutput2( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
125 #else // !_WIN32 && !_WIN64
126 if( header ) DumpOutput2( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
127 else DumpOutput2( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
128 #endif // _WIN32 || _WIN64
129 }
130 };
131
132
133
134 // Constructor usually performs only two simple tasks of filling the two lists
135 // - typeList: with all the possible id of the filtering actions
136 // - actionList with the corresponding actions. If you want to add icons to your filtering actions you can do here by construction the QActions accordingly
137
FilterScreenedPoissonPlugin()138 FilterScreenedPoissonPlugin::FilterScreenedPoissonPlugin()
139 {
140 }
141 template <class Real>
142 class PoissonParam
143 {
144 public:
145 int MaxDepthVal;
146 int MaxSolveDepthVal;
147 int KernelDepthVal;
148 int MinDepthVal;
149 int FullDepthVal;
150 Real SamplesPerNodeVal;
151 Real ScaleVal;
152 bool ConfidenceFlag;
153 bool CleanFlag;
154 bool DensityFlag;
155 Real PointWeightVal;
156 int AdaptiveExponentVal;
157 int BoundaryTypeVal;
158 bool CompleteFlag;
159 bool NonManifoldFlag;
160 bool ShowResidualFlag;
161 int CGDepthVal;
162 int ItersVal;
163 Real CSSolverAccuracyVal;
164
165 bool VerboseFlag;
166 int ThreadsVal;
167 bool LinearFitFlag;
168 float LowResIterMultiplierVal;
169 float ColorVal;
170
PoissonParam()171 PoissonParam()
172 {
173 MaxDepthVal=8;
174 MaxSolveDepthVal=-1;
175 KernelDepthVal=-1;
176 MinDepthVal=0;
177 FullDepthVal=5;
178 SamplesPerNodeVal =1.5f;
179 ScaleVal=1.1f;
180 ConfidenceFlag=false;
181 CleanFlag=false;
182 DensityFlag=false;
183 PointWeightVal = 4.0f;
184 AdaptiveExponentVal=1;
185 BoundaryTypeVal=1;
186 CompleteFlag=false;
187 NonManifoldFlag=false;
188 ShowResidualFlag=false;
189 CGDepthVal=0;
190 ItersVal=8;
191 CSSolverAccuracyVal=1e-3f;
192
193 VerboseFlag=true;
194 ThreadsVal=omp_get_num_procs();
195 LinearFitFlag = false;
196 LowResIterMultiplierVal=1.f;
197 ColorVal=16.0f;
198 }
199 };
200
201
202 template< class Real>
GetPointStreamScale(vcg::Box3<Real> & bb,float expFact)203 XForm4x4<Real> GetPointStreamScale(vcg::Box3<Real> &bb, float expFact)
204 {
205 qDebug("bbox %f %f %f - %f %f %f ",bb.min[0],bb.min[1],bb.min[2],bb.max[0],bb.max[1],bb.max[2]);
206 Real scale = bb.Dim()[bb.MaxDim()] * expFact;
207 Point3m center = bb.Center();
208 for( int i=0 ; i<3 ; i++ ) center[i] -= scale/2;
209 XForm4x4< Real > tXForm = XForm4x4< Real >::Identity() , sXForm = XForm4x4< Real >::Identity();
210 for( int i=0 ; i<3 ; i++ ) sXForm(i,i) = (Real)(1./scale ) , tXForm(3,i) = -center[i];
211 return sXForm * tXForm;
212 }
213
214
215 template< class Real >
216 class MeshModelPointStream : public OrientedPointStreamWithData< Real, Point3m >
217 {
218 CMeshO &_m;
219 size_t _curPos;
220 public:
MeshModelPointStream(CMeshO & m)221 MeshModelPointStream( CMeshO &m):_m(m),_curPos(0)
222 {
223 vcg::tri::RequireCompactness(m);
224 }
~MeshModelPointStream(void)225 ~MeshModelPointStream( void ){}
reset(void)226 void reset( void ) { _curPos =0;}
nextPoint(OrientedPoint3D<Real> & pt,Point3m & d)227 bool nextPoint( OrientedPoint3D< Real >& pt, Point3m &d)
228 {
229 if(_curPos>=_m.vn)
230 return false;
231 Point3m &nn = _m.vert[_curPos].N();
232 Point3m tp = _m.Tr * _m.vert[_curPos].P();
233 Point4m np = _m.Tr * Point4m(nn[0],nn[1],nn[2],0);
234
235 pt.p[0] = tp[0];
236 pt.p[1] = tp[1];
237 pt.p[2] = tp[2];
238 pt.n[0] = np[0];
239 pt.n[1] = np[1];
240 pt.n[2] = np[2];
241
242 d[0]=Real(_m.vert[_curPos].C()[0]);
243 d[1]=Real(_m.vert[_curPos].C()[1]);
244 d[2]=Real(_m.vert[_curPos].C()[2]);
245
246 ++_curPos;
247 return true;
248 }
249
250 };
251
252 template< class Real >
253 class MeshDocumentPointStream : public OrientedPointStreamWithData< Real, Point3m >
254 {
255 MeshDocument &_md;
256 MeshModel *_curMesh;
257 size_t _curPos;
258 size_t _totalSize;
259 public:
MeshDocumentPointStream(MeshDocument & md)260 MeshDocumentPointStream( MeshDocument &md):_md(md),_curMesh(0),_curPos(0)
261 {
262 _totalSize=0;
263 MeshModel *m=0;
264 do
265 {
266 m=_md.nextVisibleMesh(m);
267 if(m!=0)
268 {
269 vcg::tri::RequireCompactness(m->cm);
270 _totalSize+=m->cm.vn;
271 }
272 } while(m);
273 qDebug("TotalSize %i",_totalSize);
274 }
~MeshDocumentPointStream(void)275 ~MeshDocumentPointStream( void ){}
reset(void)276 void reset( void ) { _curPos =0; _curMesh=0;}
nextPoint(OrientedPoint3D<Real> & pt,Point3m & d)277 bool nextPoint( OrientedPoint3D< Real >& pt, Point3m &d )
278 {
279 Point3m nn(0,0,0);
280 // do
281 {
282 if((_curMesh==0) || (_curPos >= _curMesh->cm.vn) )
283 {
284 _curMesh = _md.nextVisibleMesh(_curMesh);
285 _curPos = 0;
286 }
287
288 if(_curMesh==0)
289 return false;
290 if(_curPos < _curMesh->cm.vn)
291 {
292 nn = _curMesh->cm.vert[_curPos].N();
293 Point3m tp = _curMesh->cm.Tr * _curMesh->cm.vert[_curPos].P();
294 Point4m np = _curMesh->cm.Tr * Point4m(nn[0],nn[1],nn[2],0);
295 // Point3m tp = _curMesh->cm.vert[_curPos].P();
296 // Point3m np = nn;
297 pt.p[0] = tp[0];
298 pt.p[1] = tp[1];
299 pt.p[2] = tp[2];
300 pt.n[0] = np[0];
301 pt.n[1] = np[1];
302 pt.n[2] = np[2];
303 d[0]=Real(_curMesh->cm.vert[_curPos].C()[0]);
304 d[1]=Real(_curMesh->cm.vert[_curPos].C()[1]);
305 d[2]=Real(_curMesh->cm.vert[_curPos].C()[2]);
306
307 ++_curPos;
308 }
309 }
310 assert(nn!=Point3m(0,0,0));
311 return true;
312 }
313
314 };
315
316
317
318
319 template< class Real , int Degree , BoundaryType BType , class Vertex >
_Execute(OrientedPointStream<Real> * pointStream,Box3m bb,CMeshO & pm,PoissonParam<Real> & pp,vcg::CallBackPos * cb)320 int _Execute(OrientedPointStream< Real > *pointStream, Box3m bb, CMeshO &pm, PoissonParam<Real> &pp, vcg::CallBackPos* cb)
321 {
322 typedef typename Octree< Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator;
323 typedef typename Octree< Real >::template InterpolationInfo< false > InterpolationInfo;
324 typedef OrientedPointStreamWithData< Real , Point3D< Real > > PointStreamWithData;
325 typedef TransformedOrientedPointStreamWithData< Real , Point3D< Real > > XPointStreamWithData;
326 Reset< Real >();
327 std::vector< char* > comments;
328
329 XForm4x4< Real > xForm = GetPointStreamScale(bb,pp.ScaleVal);
330 XForm4x4< Real > iXForm = xForm.inverse();
331 DumpOutput2( comments , "Running Screened Poisson Reconstruction (Version 9.0)\n" );
332 double startTime = Time();
333
334 OctNode< TreeNodeData >::SetAllocator(MEMORY_ALLOCATOR_BLOCK_SIZE);
335 Octree< Real > tree;
336 OctreeProfiler< Real > profiler( tree );
337 tree.threads = pp.ThreadsVal;
338 if( pp.MaxSolveDepthVal<0 ) pp.MaxSolveDepthVal = pp.MaxDepthVal;
339
340
341
342 qDebug("Using %i threads\n",pp.ThreadsVal);
343 // int kernelDepth = KernelDepth.set ? KernelDepth.value : Depth.value-2;
344 if(pp.KernelDepthVal<0) pp.KernelDepthVal =pp.MaxDepthVal-2;
345 if( pp.KernelDepthVal>pp.MaxDepthVal )
346 {
347 printf("kernelDepth cannot be greateer Depth.value\n");
348 return false;
349 }
350
351 int pointCount;
352
353 Real pointWeightSum;
354 std::vector< typename Octree< Real >::PointSample >* samples = new std::vector< typename Octree< Real >::PointSample >();
355 std::vector< ProjectiveData< Point3D< Real > , Real > >* sampleData = NULL;
356 DensityEstimator* density = NULL;
357 SparseNodeData< Point3D< Real > , NORMAL_DEGREE >* normalInfo = NULL;
358 Real targetValue = (Real)0.5;
359
360 // Read in the samples (and color data)
361 {
362 profiler.start();
363 // PointStream* pointStream;
364
365 // char* ext = GetFileExtension( In.value );
366 // if( Color.set && Color.value>0 )
367 // {
368 // sampleData = new std::vector< ProjectiveData< Point3D< Real > , Real > >();
369 // if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< Real , Point3D< Real > , float , Point3D< unsigned char > >( In.value );
370 // else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::PlyProperties , 6 , ColorInfo< Real >::ValidPlyProperties );
371 // else pointStream = new ASCIIOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::ReadASCII );
372 // }
373 // else
374 // {
375 // if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStream< Real , float >( In.value );
376 // else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStream< Real >( In.value );
377 // else pointStream = new ASCIIOrientedPointStream< Real >( In.value );
378 // }
379 // delete[] ext;
380 sampleData = new std::vector< ProjectiveData< Point3D< Real > , Real > >();
381 XPointStreamWithData _pointStream( xForm , ( PointStreamWithData& )*pointStream );
382 pointCount = tree.template init< Point3D< Real > >( _pointStream , pp.MaxDepthVal , pp.ConfidenceFlag , *samples , sampleData );
383
384 #pragma omp parallel for num_threads( pp.ThreadsVal )
385 for( int i=0 ; i<(int)samples->size() ; i++ ) (*samples)[i].sample.data.n *= (Real)-1;
386
387 DumpOutput( "Input Points / Samples: %d / %d\n" , pointCount , samples->size() );
388 profiler.dumpOutput2( comments , "# Read input into tree:" );
389 }
390
391 DenseNodeData< Real , Degree > solution;
392
393 {
394 DenseNodeData< Real , Degree > constraints;
395 InterpolationInfo* iInfo = NULL;
396 int solveDepth = pp.MaxSolveDepthVal;
397
398 tree.resetNodeIndices();
399
400 // Get the kernel density estimator [If discarding, compute anew. Otherwise, compute once.]
401 {
402 profiler.start();
403 density = tree.template setDensityEstimator< WEIGHT_DEGREE >( *samples , pp.KernelDepthVal , pp.SamplesPerNodeVal );
404 profiler.dumpOutput2( comments , "# Got kernel density:" );
405 }
406
407 // Transform the Hermite samples into a vector field [If discarding, compute anew. Otherwise, compute once.]
408 {
409 profiler.start();
410 normalInfo = new SparseNodeData< Point3D< Real > , NORMAL_DEGREE >();
411 *normalInfo = tree.template setNormalField< NORMAL_DEGREE >( *samples , *density , pointWeightSum , BType==BOUNDARY_NEUMANN );
412 profiler.dumpOutput2( comments , "# Got normal field:" );
413 }
414
415 if( !pp.DensityFlag ) delete density , density = NULL;
416
417 // Trim the tree and prepare for multigrid
418 {
419 profiler.start();
420 std::vector< int > indexMap;
421
422 constexpr int MAX_DEGREE = NORMAL_DEGREE > Degree ? NORMAL_DEGREE : Degree;
423 tree.template inalizeForBroodedMultigrid< MAX_DEGREE , Degree , BType >( pp.FullDepthVal , typename Octree< Real >::template HasNormalDataFunctor< NORMAL_DEGREE >( *normalInfo ) , &indexMap );
424
425 if( normalInfo ) normalInfo->remapIndices( indexMap );
426 if( density ) density->remapIndices( indexMap );
427 profiler.dumpOutput2( comments , "# Finalized tree:" );
428 }
429
430 // Add the FEM constraints
431 {
432 profiler.start();
433 constraints = tree.template initDenseNodeData< Degree >( );
434 tree.template addFEMConstraints< Degree , BType , NORMAL_DEGREE , BType >( FEMVFConstraintFunctor< NORMAL_DEGREE , BType , Degree , BType >( 1. , 0. ) , *normalInfo , constraints , solveDepth );
435 profiler.dumpOutput2( comments , "# Set FEM constraints:" );
436 }
437
438 // Free up the normal info [If we don't need it for subseequent iterations.]
439 delete normalInfo , normalInfo = NULL;
440
441 // Add the interpolation constraints
442 if( pp.PointWeightVal>0 )
443 {
444 profiler.start();
445 iInfo = new InterpolationInfo( tree , *samples , targetValue , pp.AdaptiveExponentVal , (Real)pp.PointWeightVal * pointWeightSum , (Real)0 );
446 tree.template addInterpolationConstraints< Degree , BType >( *iInfo , constraints , solveDepth );
447 profiler.dumpOutput2( comments , "#Set point constraints:" );
448 }
449
450 DumpOutput( "Leaf Nodes / Active Nodes / Ghost Nodes: %d / %d / %d\n" , (int)tree.leaves() , (int)tree.nodes() , (int)tree.ghostNodes() );
451 DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) );
452
453 // Solve the linear system
454 {
455 profiler.start();
456 typename Octree< Real >::SolverInfo solverInfo;
457 solverInfo.cgDepth = pp.CGDepthVal , solverInfo.iters = pp.ItersVal , solverInfo.cgAccuracy = pp.CSSolverAccuracyVal , solverInfo.verbose = pp.VerboseFlag , solverInfo.showResidual = pp.ShowResidualFlag , solverInfo.lowResIterMultiplier = std::max< double >( 1. , pp.LowResIterMultiplierVal );
458 solution = tree.template solveSystem< Degree , BType >( FEMSystemFunctor< Degree , BType >( 0 , 1. , 0 ) , iInfo , constraints , solveDepth , solverInfo );
459 profiler.dumpOutput2( comments , "# Linear system solved:" );
460 if( iInfo ) delete iInfo , iInfo = NULL;
461 }
462 }
463
464 CoredFileMeshData< Vertex > mesh;
465
466 {
467 profiler.start();
468 double valueSum = 0 , weightSum = 0;
469 typename Octree< Real >::template MultiThreadedEvaluator< Degree , BType > evaluator( &tree , solution , pp.ThreadsVal );
470 #pragma omp parallel for num_threads( pp.ThreadsVal ) reduction( + : valueSum , weightSum )
471 for( int j=0 ; j<samples->size() ; j++ )
472 {
473 ProjectiveData< OrientedPoint3D< Real > , Real >& sample = (*samples)[j].sample;
474 Real w = sample.weight;
475 if( w>0 ) weightSum += w , valueSum += evaluator.value( sample.data.p / sample.weight , omp_get_thread_num() , (*samples)[j].node ) * w;
476 }
477 Real isoValue = (Real)( valueSum / weightSum );
478 // if( samples ) delete samples , samples = NULL;
479 profiler.dumpOutput( "Got average:" );
480 DumpOutput( "Iso-Value: %e\n" , isoValue );
481
482 profiler.start();
483 SparseNodeData< ProjectiveData< Point3D< Real > , Real > , DATA_DEGREE >* colorData = NULL;
484 if( sampleData )
485 {
486 colorData = new SparseNodeData< ProjectiveData< Point3D< Real > , Real > , DATA_DEGREE >();
487 *colorData = tree.template setDataField< DATA_DEGREE , false >( *samples , *sampleData , (DensityEstimator*)NULL );
488 delete sampleData , sampleData = NULL;
489 for( const OctNode< TreeNodeData >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) )
490 {
491 ProjectiveData< Point3D< Real > , Real >* clr = (*colorData)( n );
492 if( clr )
493 (*clr) *= (Real)pow( pp.ColorVal , tree.depth( n ) );
494 }
495 }
496 tree.template getMCIsoSurface< Degree , BType , WEIGHT_DEGREE , DATA_DEGREE >( density , colorData , solution , isoValue , mesh , !pp.LinearFitFlag , !pp.NonManifoldFlag , false /*PolygonMesh.set*/ );
497 DumpOutput( "Vertices / Polygons: %d / %d\n" , mesh.outOfCorePointCount()+mesh.inCorePoints.size() , mesh.polygonCount() );
498 profiler.dumpOutput2( comments , "# Got triangles:" );
499 }
500
501 // FreePointer( solution );
502
503 cb(90,"Creating Mesh");
504 mesh.resetIterator();
505 int vm = mesh.outOfCorePointCount()+mesh.inCorePoints.size();
506 for(auto pt=mesh.inCorePoints.begin();pt!=mesh.inCorePoints.end();++pt)
507 {
508 Point3D<Real> pp = iXForm*pt->point;
509 vcg::tri::Allocator<CMeshO>::AddVertex(pm,Point3m(pp[0],pp[1],pp[2]));
510 pm.vert.back().Q() = pt->value;
511 pm.vert.back().C()[0] = pt->color[0];
512 pm.vert.back().C()[1] = pt->color[1];
513 pm.vert.back().C()[2] = pt->color[2];
514 }
515 for (int ii=0; ii < mesh.outOfCorePointCount(); ii++){
516 Vertex pt;
517 mesh.nextOutOfCorePoint(pt);
518 Point3D<Real> pp = iXForm*pt.point;
519 vcg::tri::Allocator<CMeshO>::AddVertex(pm,Point3m(pp[0],pp[1],pp[2]));
520 pm.vert.back().Q() = pt.value;
521 pm.vert.back().C()[0] = pt.color[0];
522 pm.vert.back().C()[1] = pt.color[1];
523 pm.vert.back().C()[2] = pt.color[2];
524 }
525
526 std::vector< CoredVertexIndex > polygon;
527 while(mesh.nextPolygon( polygon ))
528 {
529 assert(polygon.size()==3);
530 int indV[3];
531 for( int i=0 ; i<int(polygon.size()) ; i++ )
532 {
533 if( polygon[i].inCore ) indV[i] = polygon[i].idx;
534 else indV[i]= polygon[i].idx + int( mesh.inCorePoints.size() );
535 }
536 vcg::tri::Allocator<CMeshO>::AddFace(pm, &pm.vert[indV[0]], &pm.vert[indV[1]], &pm.vert[indV[2]]);
537 }
538 cb(100,"Done");
539
540 // if( colorData ) delete colorData , colorData = NULL;
541
542
543 if( density ) delete density , density = NULL;
544 DumpOutput2( comments , "# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-startTime , tree.maxMemoryUsage() );
545
546 return 1;
547 }
548
549
550
551 template <class MeshType>
PoissonClean(MeshType & m,bool scaleNormal,bool cleanFlag)552 void PoissonClean(MeshType &m, bool scaleNormal, bool cleanFlag)
553 {
554 vcg::tri::UpdateNormal<MeshType>::NormalizePerVertex(m);
555
556 if(cleanFlag) {
557 for (auto vi = m.vert.begin(); vi != m.vert.end(); ++vi)
558 if (vcg::SquaredNorm(vi->N()) < std::numeric_limits<float>::min()*10.0)
559 vcg::tri::Allocator<MeshType>::DeleteVertex(m,*vi);
560
561 for (auto fi = m.face.begin(); fi != m.face.end(); ++fi)
562 if( fi->V(0)->IsD() || fi->V(1)->IsD() || fi->V(2)->IsD() )
563 vcg::tri::Allocator<MeshType>::DeleteFace(m,*fi);
564 }
565
566 vcg::tri::Allocator<MeshType>::CompactEveryVector(m);
567 if(scaleNormal)
568 {
569 for(auto vi=m.vert.begin();vi!=m.vert.end();++vi)
570 vi->N() *= vi->Q();
571 }
572 }
573
HasGoodNormal(CMeshO & m)574 bool HasGoodNormal(CMeshO &m)
575 {
576 for(auto vi=m.vert.begin();vi!=m.vert.end();++vi)
577 if(vcg::SquaredNorm(vi->N()) < std::numeric_limits<float>::min()*10.0)
578 return false;
579
580 return true;
581 }
582
applyFilter(const QString & filterName,MeshDocument & md,EnvWrap & env,vcg::CallBackPos * cb)583 bool FilterScreenedPoissonPlugin::applyFilter( const QString& filterName,MeshDocument& md,EnvWrap& env, vcg::CallBackPos* cb)
584 {
585 bool currDirChanged=false;
586 QDir currDir = QDir::current();
587
588 if (filterName == "Surface Reconstruction: Screened Poisson")
589 {
590 //check if folder is writable
591 QTemporaryFile file("./_tmp_XXXXXX.tmp");
592 if (!file.open())
593 {
594 currDirChanged=true;
595 QTemporaryDir tmpdir;
596 QDir::setCurrent(tmpdir.path());
597 Log("Warning - current folder is not writable. Screened Poisson Merging needs to save intermediate files in the current working folder. Project and meshes must be in a write-enabled folder. Please save your data in a suitable folder before applying.");
598 //errorMessage = "current folder is not writable.<br> Screened Poisson Merging needs to save intermediate files in the current working folder.<br> Project and meshes must be in a write-enabled folder.<br> Please save your data in a suitable folder before applying.";
599 //return false;
600 }
601
602 PoissonParam<Scalarm> pp;
603 pp.MaxDepthVal = env.evalInt("depth");
604 pp.FullDepthVal = env.evalInt("fullDepth");
605 pp.CGDepthVal= env.evalInt("cgDepth");
606 pp.ScaleVal = env.evalFloat("scale");
607 pp.SamplesPerNodeVal = env.evalFloat("samplesPerNode");
608 pp.PointWeightVal = env.evalFloat("pointWeight");
609 pp.ItersVal = env.evalInt("iters");
610 pp.ConfidenceFlag = env.evalBool("confidence");
611 pp.DensityFlag = true;
612 pp.CleanFlag = env.evalBool("preClean");
613
614 bool goodNormal=true, goodColor=true;
615 if(env.evalBool("visibleLayer") == false)
616 {
617 PoissonClean(md.mm()->cm, pp.ConfidenceFlag, pp.CleanFlag);
618 goodNormal=HasGoodNormal(md.mm()->cm);
619 goodColor = md.mm()->hasDataMask(MeshModel::MM_VERTCOLOR);
620 }
621 else
622 {
623 MeshModel *_mm=0;
624 while(_mm=md.nextVisibleMesh(_mm)) {
625 PoissonClean(_mm->cm, pp.ConfidenceFlag, pp.CleanFlag);
626 goodNormal &= HasGoodNormal(_mm->cm);
627 goodColor &= _mm->hasDataMask(MeshModel::MM_VERTCOLOR);
628 }
629 }
630
631 if(!goodNormal)
632 {
633 this->errorMessage = "Filter requires correct per vertex normals.<br>"
634 "E.g. it is necessary that your <b>ALL</b> the input vertices have a proper, not-null normal.<br> "
635 "Try enabling the <i>pre-clean<i> option and retry.<br><br>"
636 "To permanently remove this problem:<br>"
637 "If you encounter this error on a triangulated mesh try to use the <i>Remove Unreferenced Vertices</i> filter"
638 "If you encounter this error on a pointcloud try to use the <i>Conditional Vertex Selection</i> filter"
639 "with function '(nx==0.0) && (ny==0.0) && (nz==0.0)', and then <i>delete selected vertices</i>.<br>";
640 return false;
641 }
642
643 MeshModel *pm =md.addNewMesh("","Poisson mesh",false);
644 md.setVisible(pm->id(),false);
645 pm->updateDataMask(MeshModel::MM_VERTQUALITY);
646 if(goodColor) pm->updateDataMask(MeshModel::MM_VERTCOLOR);
647
648 if(env.evalBool("visibleLayer"))
649 {
650 Box3m bb;
651 MeshModel *_mm=0;
652 while(_mm=md.nextVisibleMesh(_mm))
653 bb.Add(_mm->cm.Tr,_mm->cm.bbox);
654
655 MeshDocumentPointStream<Scalarm> documentStream(md);
656 _Execute<Scalarm,2,BOUNDARY_NEUMANN,PlyColorAndValueVertex<Scalarm> >(&documentStream,bb,pm->cm,pp,cb);
657 }
658 else
659 {
660 MeshModelPointStream<Scalarm> meshStream(md.mm()->cm);
661 _Execute<Scalarm,2,BOUNDARY_NEUMANN,PlyColorAndValueVertex<Scalarm> >(&meshStream,md.mm()->cm.bbox,pm->cm,pp,cb);
662 }
663 pm->UpdateBoxAndNormals();
664 md.setVisible(pm->id(),true);
665 md.setCurrentMesh(pm->id());
666 if(currDirChanged) QDir::setCurrent(currDir.path());
667 return true;
668 }
669 return false;
670 }
671
672
673
674
675
676 MESHLAB_PLUGIN_NAME_EXPORTER(FilterScreenedPoissonPlugin)
677