DGtalTools 2.1.0
Loading...
Searching...
No Matches
mesh2heightfield.cpp
1
30#include <iostream>
31#include <fstream>
32#include "DGtal/base/Common.h"
33#include "DGtal/helpers/StdDefs.h"
34#include "DGtal/images/ImageContainerBySTLVector.h"
35#include "DGtal/io/writers/GenericWriter.h"
36#include "DGtal/io/readers/MeshReader.h"
37#include "DGtal/io/writers/MeshWriter.h"
38#include "DGtal/images/ConstImageAdapter.h"
39#include "DGtal/kernel/BasicPointFunctors.h"
40#include <DGtal/base/BasicFunctors.h>
41#include "DGtal/math/linalg/SimpleMatrix.h"
42
43#include <DGtal/math/linalg/SimpleMatrix.h>
44#include <DGtal/math/linalg/EigenDecomposition.h>
45
46
47
48#include "CLI11.hpp"
49
50
51
52using namespace std;
53using namespace DGtal;
54
55
56
104typedef PointVector<9, double> Matrix3x3Point;
105typedef SimpleMatrix<double, 3, 3 > CoVarianceMat;
106
110static CoVarianceMat
111getCoVarianceMatFrom(const Matrix3x3Point &aMatrixPt){
112 CoVarianceMat res;
113 for(unsigned int i = 0; i<3; i++){
114 for(unsigned int j = 0; j<3; j++){
115 res.setComponent(i, j, aMatrixPt[j+i*3] );
116 }
117 }
118 return res;
119}
120
121
122std::pair<DGtal::Z3i::RealPoint, DGtal::Z3i::RealPoint>
123getMainDirsCoVar(Mesh<Z3i::RealPoint>::Iterator begin, Mesh<Z3i::RealPoint>::Iterator end )
124{
125 std::pair<DGtal::Z3i::RealPoint, DGtal::Z3i::RealPoint> res;
126 Matrix3x3Point c ;
127 unsigned int nb = 0;
128 //compute centroïd of point cloud
129 Z3i::RealPoint centroid;
130 int cp=0;
131 for(auto it=begin; it != end; it++){
132 Z3i::RealPoint pt=*it;
133 centroid[0] += pt[0];
134 centroid[1] += pt[1];
135 centroid[2] += pt[2];
136 ++cp;
137 }
138 centroid /=cp;
139 //compute matrix like PCL : https://pointclouds.org/documentation/moment__of__inertia__estimation_8hpp_source.html
140 for(auto it=begin; it != end; it++){
141
142 Z3i::RealPoint pt=*it;
143 pt[0] -= centroid[0];
144 pt[1] -= centroid[1];
145 pt[2] -= centroid[2];
146
147 trace.info()<<pt<< std::endl;
148 c[4] += pt[1] * pt[1];
149 c[5] += pt[1] * pt[2];
150 c[8] += pt[2] * pt[2];
151 pt *= pt[0];
152 c[0] += pt[0];
153 c[1] += pt[1];
154 c[2] += pt[2];
155
156 nb++;
157
158 }
159 c[3] = c[1];
160 c[6] = c[2];
161 c[7] = c[5];
162 //normalize matrix
163 c = c/nb;
164 CoVarianceMat covar = getCoVarianceMatFrom(c);
165 SimpleMatrix<double, 3, 3 > eVects;
166 PointVector<3, double> eVals;
167 DGtal::EigenDecomposition<3, double, CoVarianceMat>::getEigenDecomposition (covar, eVects, eVals);
168 unsigned int temp = 0;
169 unsigned int major_index = 0;
170 unsigned int middle_index = 1;
171 unsigned int minor_index = 2;
172 //reorder first and second axis according to their eigen value
173 if (eVals(major_index) < eVals(middle_index))
174 {
175 temp = major_index;
176 major_index = middle_index;
177 middle_index = temp;
178 }
179 if (eVals(major_index) < eVals(minor_index))
180 {
181 temp = major_index;
182 major_index = minor_index;
183 minor_index = temp;
184 }
185 if (eVals(middle_index) < eVals(minor_index))
186 {
187 temp = minor_index;
188 minor_index = middle_index;
189 middle_index = temp;
190 }
191
192 res.first = eVects.column(major_index);
193 res.second = eVects.column(middle_index);
194
195 return res;
196}
197
198template<typename TPoint, typename TEmbeder>
199TPoint
200getNormal(const TPoint &vect, const TEmbeder &emb ){
201 Z3i::RealPoint p0 = emb(Z2i::RealPoint(0.0,0.0), false);
202 Z3i::RealPoint px = emb(Z2i::RealPoint(10.0,0.0), false);
203 Z3i::RealPoint py = emb(Z2i::RealPoint(0.0,10.0), false);
204 Z3i::RealPoint u = px-p0; u /= u.norm();
205 Z3i::RealPoint v = py-p0; v /= v.norm();
206 Z3i::RealPoint w = u.crossProduct(v); w /= w.norm();
207 SimpleMatrix<double, 3, 3> t;
208 t.setComponent(0,0, u[0]); t.setComponent(0,1, v[0]); t.setComponent(0, 2, w[0]);
209 t.setComponent(1,0, u[1]); t.setComponent(1,1, v[1]); t.setComponent(1, 2, w[1]);
210 t.setComponent(2,0, u[2]); t.setComponent(2,1, v[2]); t.setComponent(2, 2, w[2]);
211 return ((t.inverse())*vect);
212}
213
214
215
216
217template<typename TImage, typename TVectorField, typename TPoint>
218void
219fillPointArea(TImage &anImage, TVectorField &imageVectorField,
220 const TPoint aPoint, const unsigned int size, const unsigned int value, const Z3i::RealPoint &n){
221 typename TImage::Domain aDom(aPoint-TPoint::diagonal(size), aPoint+TPoint::diagonal(size));
222 for(typename TImage::Domain::ConstIterator it= aDom.begin(); it != aDom.end(); it++){
223 if (anImage.domain().isInside(*it)){
224 anImage.setValue(*it, value);
225 imageVectorField.setValue(*it, n);
226 }
227 }
228}
229
230
231int main( int argc, char** argv )
232{
233 typedef ImageContainerBySTLVector < Z3i::Domain, Z3i::RealPoint > VectorFieldImage3D;
234 typedef ImageContainerBySTLVector < Z2i::Domain, Z3i::RealPoint > VectorFieldImage2D;
235 typedef ImageContainerBySTLVector < Z3i::Domain, unsigned char > Image3D;
236 typedef ImageContainerBySTLVector < Z2i::Domain, unsigned char> Image2D;
237 typedef DGtal::ConstImageAdapter<Image3D, Z2i::Domain, DGtal::functors::Point2DEmbedderIn3D<DGtal::Z3i::Domain>,
238 Image3D::Value, DGtal::functors::Identity > ImageAdapterExtractor;
239
240
241// parse command line using CLI ----------------------------------------------
242 CLI::App app;
243 std::string inputFileName;
244 std::string outputFileName {"result.pgm"};
245 double meshScale = 1.0;
246 double triangleAreaUnit = 0.01;
247 double widthImageScan = {500};
248 double heightImageScan = {500};
249 bool orientAutoFrontX = false;
250 bool orientAutoFrontY = false;
251 bool orientAutoFrontZ = false;
252 bool orientAuto = true;
253 bool invertNormal = false;
254 bool orientBack = false;
255 bool exportNormals = false;
256 bool setBackgroundLastDepth = false;
257 bool backgroundNormalBack = false;
258 int maxScan {255};
259 int centerX {0};
260 int centerY {0};
261 int centerZ {200};
262
263
264 double nx{0};
265 double ny{0};
266 double nz{1};
267 DGtal::Z3i::RealPoint secDir; //orientation to be used when flag orientAuto is used
268 unsigned int minDiagVolSize {400}; //min size used to automatically resize mesh if included in unit box.
269 app.description("Convert a mesh file into a projected 2D image given from a normal direction N and from a starting point P. The 3D mesh discretized and scanned in the normal direction N, starting from P with a step 1.\n Example:\n mesh2heightfield -i ${DGtal}/examples/samples/tref.off heighfield.pgm \n");
270 app.add_option("-i,--input,1", inputFileName, "mesh file (.off)" )
271 ->required()
272 ->check(CLI::ExistingFile);
273 app.add_option("-o,--output,2", outputFileName, "sequence of discrete point file (.sdp) " );
274 app.add_option("--meshScale,-s", meshScale, "change the default mesh scale (each vertex multiplied by the scale) ");
275 app.add_option("--remeshMinArea,-a", triangleAreaUnit, "ajust the remeshing min triangle are used to avoid empty areas");
276
277 app.add_option("--heightFieldMaxScan", maxScan, "set the maximal scan deep." );
278 app.add_option("-x,--centerX",
279 centerX, "choose x center of the projected image.");
280 app.add_option("-y,--centerY",
281 centerY, "choose y center of the projected image.");
282 app.add_option("-z,--centerZ",
283 centerZ, "choose z center of the projected image.");
284 auto optNx = app.add_option("--nx",
285 nx, "set the x component of the projection direction.");
286 auto optNy = app.add_option("--ny",
287 ny, "set the y component of the projection direction.");
288 auto optNz = app.add_option("--nz",
289 nz, "set the z component of the projection direction.");
290 app.add_flag("--invertNormals,-v", invertNormal, "invert normal vector of the mesh");
291 app.add_option("--width", widthImageScan, "set the width of the area to be extracted as an height field image. (note that the resulting image width also depends of the scale parameter (option --meshScale))" );
292 app.add_option("--height", heightImageScan, "set the height of the area to extracted as an height field image. (note that the resulting image height also depends of the scale parameter (option --meshScale))" );
293 app.add_flag("--orientAutoFrontX", orientAutoFrontX,"automatically orients the camera in front according the x axis." );
294 app.add_flag("--orientAutoFrontY", orientAutoFrontY,"automatically orients the camera in front according the y axis." );
295 app.add_flag("--orientAutoFrontZ", orientAutoFrontZ,"automatically orients the camera in front according the z axis." );
296 app.add_flag("--orientBack", orientBack, "change the camera direction to back instead front as given in orientAutoFront{X,Y,Z} options.");
297 app.add_flag("--exportNormals", exportNormals, "export mesh normal vectors (given in the image height field basis).");
298 app.add_flag("--backgroundNormalBack", backgroundNormalBack, "set the normals of background in camera opposite direction (to obtain a black background in rendering).");
299 app.add_flag("--setBackgroundLastDepth", setBackgroundLastDepth, "change the default background (black with the last filled intensity).");
300
301 app.get_formatter()->column_width(40);
302 CLI11_PARSE(app, argc, argv);
303 // END parse command line using CLI ----------------------------------------------
304
305 orientAuto = optNx->count() == 0 && optNy->count() == 0 && optNy->count() == 0;
306
307 trace.info() << "Reading input file " << inputFileName ;
308 Mesh<Z3i::RealPoint> inputMesh(true);
309 inputMesh << inputFileName;
310 std::pair<Z3i::RealPoint, Z3i::RealPoint> b = inputMesh.getBoundingBox();
311 double diagDist = (b.first-b.second).norm();
312 if(diagDist<minDiagVolSize){
313 meshScale = minDiagVolSize/diagDist;
314 inputMesh.rescale(meshScale);
315 }
316
317 triangleAreaUnit *= meshScale;
318 // get vertex
319 inputMesh.vertexBegin();
320
321 inputMesh.quadToTriangularFaces();
322 trace.info() << " [done] " << std::endl ;
323 int maxArea = triangleAreaUnit+1.0 ;
324 while (maxArea> triangleAreaUnit)
325 {
326 trace.info()<< "Iterating mesh subdivision ... "<< maxArea;
327 maxArea = inputMesh.subDivideTriangularFaces(triangleAreaUnit);
328 trace.info() << " [done]"<< std::endl;
329 }
330
331 std::pair<Z3i::RealPoint, Z3i::RealPoint> bb = inputMesh.getBoundingBox();
332 Image3D::Domain meshDomain( bb.first-Z3i::Point::diagonal(1),
333 bb.second+Z3i::Point::diagonal(1));
334
335 // vol image filled from the mesh vertex.
336 Image3D meshVolImage(meshDomain);
337 VectorFieldImage3D meshNormalImage(meshDomain);
338 Z3i::RealPoint z(0.0,0.0,0.0);
339 for(Image3D::Domain::ConstIterator it = meshVolImage.domain().begin();
340 it != meshVolImage.domain().end(); it++)
341 {
342 meshVolImage.setValue(*it, 0);
343 meshNormalImage.setValue(*it, z);
344 }
345
346 // Filling mesh faces in volume.
347 for(unsigned int i =0; i< inputMesh.nbFaces(); i++)
348 {
349 trace.progressBar(i, inputMesh.nbFaces());
350 typename Mesh<Z3i::RealPoint>::MeshFace aFace = inputMesh.getFace(i);
351 if(aFace.size()==3)
352 {
353 Z3i::RealPoint p1 = inputMesh.getVertex(aFace[0]);
354 Z3i::RealPoint p2 = inputMesh.getVertex(aFace[1]);
355 Z3i::RealPoint p3 = inputMesh.getVertex(aFace[2]);
356 Z3i::RealPoint n = (p2-p1).crossProduct(p3-p1);
357 n /= invertNormal? -n.norm(): n.norm();
358 Z3i::RealPoint c = (p1+p2+p3)/3.0;
359 fillPointArea(meshVolImage, meshNormalImage, p1, 1, 1, n);
360 fillPointArea(meshVolImage, meshNormalImage, p2, 1, 1, n);
361 fillPointArea(meshVolImage, meshNormalImage, p3, 1, 1, n);
362 fillPointArea(meshVolImage, meshNormalImage, c, 1, 1, n);
363 }
364 }
365
366
367 if(orientAutoFrontX || orientAutoFrontY || orientAutoFrontZ)
368 {
369 Z3i::Point ptL = meshVolImage.domain().lowerBound();
370 Z3i::Point ptU = meshVolImage.domain().upperBound();
371 Z3i::Point ptC = (ptL+ptU)/2;
372 centerX=ptC[0]; centerY=ptC[1]; centerZ=ptC[2];
373 nx=0; ny=0; nz=0;
374 }
375
376 if(orientAutoFrontX)
377 {
378 nx=(orientBack?-1.0:1.0);
379 maxScan = meshVolImage.domain().upperBound()[0]- meshVolImage.domain().lowerBound()[0];
380 centerX = centerX + (orientBack? maxScan/2: -maxScan/2) ;
381 }
382 if(orientAutoFrontY)
383 {
384 ny=(orientBack?-1.0:1.0);
385 maxScan = meshVolImage.domain().upperBound()[1]- meshVolImage.domain().lowerBound()[1];
386 centerY = centerY + (orientBack? maxScan/2: -maxScan/2);
387 }
388 if(orientAutoFrontZ)
389 {
390 nz=(orientBack?-1.0:1.0);
391 maxScan = meshVolImage.domain().upperBound()[2]-meshVolImage.domain().lowerBound()[2];
392 centerZ = centerZ + (orientBack? maxScan/2: -maxScan/2);
393 }
394 if (orientAuto)
395 {
396 Z3i::Point ptL = meshVolImage.domain().lowerBound();
397 Z3i::Point ptU = meshVolImage.domain().upperBound();
398 Z3i::Point ptC = (ptL+ptU)/2;
399 int maxScanZ = meshVolImage.domain().upperBound()[2]-meshVolImage.domain().lowerBound()[2];
400 int maxScanY = meshVolImage.domain().upperBound()[1]-meshVolImage.domain().lowerBound()[1];
401 int maxScanX = meshVolImage.domain().upperBound()[0]-meshVolImage.domain().lowerBound()[0];
402 maxScan = std::max(std::max(maxScanX, maxScanY), maxScanZ);
403 trace.info() << "Automatic orientation on volume of bounds: " << ptL << " " << ptU << std::endl;
404 trace.info() << "Computing main dir...";
405 auto dir = getMainDirsCoVar(inputMesh.vertexBegin(), inputMesh.vertexEnd());
406 trace.info() << "[done]"<<std::endl;
407 nx = dir.first[0];
408 ny = dir.first[1];
409 nz = dir.first[2];
410 secDir = dir.second;
411 centerX=ptC[0]-(double)maxScan*nx/2;
412 centerY=ptC[1]-(double)maxScan*ny/2;
413 centerZ=ptC[2]-(double)maxScan*nz/2;
414
415
416 }
417 functors::Rescaling<unsigned int, Image2D::Value> scaleFctDepth(0, maxScan, 0, 255);
418 if(maxScan > std::numeric_limits<Image2D::Value>::max())
419 {
420 trace.info()<< "Max depth value outside image intensity range: " << maxScan
421 << " (use a rescaling functor which implies a loss of precision)" << std::endl;
422 }
423 trace.info() << "Processing image to output file " << outputFileName;
424
425 Image2D::Domain aDomain2D(DGtal::Z2i::Point(0,0),
426 DGtal::Z2i::Point(widthImageScan, heightImageScan));
427 Z3i::Point ptCenter (centerX, centerY, centerZ);
428
429 Z3i::RealPoint normalDir (nx, ny, nz);
430 Image2D resultingImage(aDomain2D);
431 VectorFieldImage2D resultingVectorField(aDomain2D);
432
433
434 for(Image2D::Domain::ConstIterator it = resultingImage.domain().begin();
435 it != resultingImage.domain().end(); it++){
436 resultingImage.setValue(*it, 0);
437 resultingVectorField.setValue(*it, z);
438 }
439 DGtal::functors::Identity idV;
440
441 unsigned int maxDepthFound = 0;
442 Z3i::Point c (ptCenter-normalDir, DGtal::functors::Round<>());
443 DGtal::functors::Point2DEmbedderIn3D<DGtal::Z3i::Domain > embedder(meshVolImage.domain(),
444 c,
445 normalDir,
446 widthImageScan);
447 if (orientAuto){
448 embedder = DGtal::functors::Point2DEmbedderIn3D<DGtal::Z3i::Domain >(meshVolImage.domain(),
449 c,
450 normalDir,
451 secDir,
452 widthImageScan);
453 }
454 int k = 0;
455 bool firstFound = false;
456 while (k < maxScan || !firstFound){
457 embedder.shiftOriginPoint(normalDir);
458 trace.progressBar(k, maxScan);
459 ImageAdapterExtractor extractedImage(meshVolImage, aDomain2D, embedder, idV);
460 for(Image2D::Domain::ConstIterator it = extractedImage.domain().begin();
461 it != extractedImage.domain().end(); it++)
462 {
463 if(resultingImage(*it)== 0 && extractedImage(*it)!=0)
464 {
465
466 if (!firstFound) {
467 firstFound = true;
468 }
469 maxDepthFound = k;
470 resultingImage.setValue(*it, scaleFctDepth(maxScan-k));
471 resultingVectorField.setValue(*it, getNormal(meshNormalImage(embedder(*it)), embedder));
472 }
473 }
474 if (firstFound){
475 k++;
476 }
477 }
478 if (setBackgroundLastDepth)
479 {
480 for(Image2D::Domain::ConstIterator it = resultingImage.domain().begin();
481 it != resultingImage.domain().end(); it++){
482 if( resultingImage(*it)== 0 )
483 {
484 resultingImage.setValue(*it, scaleFctDepth(maxScan-maxDepthFound));
485 }
486 }
487 }
488 bool inverBgNormal = backgroundNormalBack;
489 for(Image2D::Domain::ConstIterator it = resultingImage.domain().begin();
490 it != resultingImage.domain().end(); it++){
491 if(resultingVectorField(*it)==Z3i::RealPoint(0.0, 0.0, 0.0))
492 {
493 resultingVectorField.setValue(*it,Z3i::RealPoint(0, 0, inverBgNormal?-1: 1));
494 }
495 }
496 resultingImage >> outputFileName;
497 if(exportNormals){
498 std::stringstream ss;
499 ss << outputFileName << ".normals";
500 std::ofstream outN;
501 outN.open(ss.str().c_str(), std::ofstream::out);
502 for(Image2D::Domain::ConstIterator it = resultingImage.domain().begin();
503 it != resultingImage.domain().end(); it++){
504 outN << (*it)[0] << " " << (*it)[1] << " " << 0 << std::endl;
505 outN <<(*it)[0]+ resultingVectorField(*it)[0] << " "
506 <<(*it)[1]+resultingVectorField(*it)[1] << " "
507 << resultingVectorField(*it)[2];
508 outN << std::endl;
509 }
510 }
511 return EXIT_SUCCESS;;
512}
DGtal
Definition ATu0v1.h:57