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