Skip to content

Commit 02a1be1

Browse files
jmmullerjmmuller
committed
staticlidar: add e57 format
1 parent b7afa13 commit 02a1be1

File tree

3 files changed

+216
-20
lines changed

3 files changed

+216
-20
lines changed

‎MMVII/CMakeLists.txt‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,14 @@ if((NOT OpenMP_CXX_FOUND) AND (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))
7575
message(STATUS "OpenMP not found: did you install libomp-dev ?")
7676
endif()
7777

78+
# libE57Format
79+
option(USE_LIBE57FORMAT "If ON, e57 point cloud file format will be available" OFF)
80+
if(USE_LIBE57FORMAT)
81+
set(LIBE57_INCLUDE_DIR "" CACHE PATH "Path to libe57format includes")
82+
set(LIBE57_LIBRARY "" CACHE FILEPATH "Path to libe57format .a library")
83+
find_package( XercesC 3.2 REQUIRED )
84+
endif()
85+
7886
#######################################################
7987
## Compilation toolchain
8088
#######################################################
@@ -143,6 +151,19 @@ add_library(P2007 "${mmv2_libsrcs}" "${mmv2_headers}")
143151
target_include_directories(P2007 PRIVATE "${mmv2_include_dir};${mmv1_include_dir};${EIGEN3_INCLUDE_PATH};${PROJ_INCLUDE_DIRS};${GDAL_INCLUDE_DIRS}")
144152
target_link_libraries(P2007 PRIVATE MMVII_compiler_flags Threads::Threads ${PROJ_LIBRARIES} ${GDAL_LIBRARIES})
145153

154+
155+
if(USE_LIBE57FORMAT)
156+
target_compile_definitions(P2007 PUBLIC USE_LIBE57FORMAT=true)
157+
if(NOT EXISTS "${LIBE57_INCLUDE_DIR}")
158+
message(FATAL_ERROR "Path to libe57format includes is not correct: ${LIBE57_INCLUDE_DIR}")
159+
endif()
160+
if(NOT EXISTS "${LIBE57_LIBRARY}")
161+
message(FATAL_ERROR "Path to libe57format .a library is not correct: ${LIBE57_LIBRARY}")
162+
endif()
163+
target_include_directories(P2007 PRIVATE "${LIBE57_INCLUDE_DIR}")
164+
target_link_libraries(P2007 PRIVATE "${LIBE57_LIBRARY}" XercesC::XercesC)
165+
endif()
166+
146167
# Force parallel compilation with MSVC
147168
if(MSVC)
148169
target_compile_options(P2007 PRIVATE "/MP")

‎MMVII/README.md‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ Some external tools need to be present on your system for **MMVII** to run prope
3131
- **[PROJ](http://trac.osgeo.org/proj/)** for coordinate system conversion and coordinate reference system transformation
3232
- **[PROJ additional data](https://download.osgeo.org/proj/)** grids for coordinates tranformations (optional, see doc)
3333
- **[GDAL](https://gdal.org/)** for image files handling
34+
- **[libE57Format](https://github.com/asmaloney/libE57Format)** for e57 point cloud file format (optional)
3435
- **[ccache](https://ccache.dev/)** for recompilation optimization (optional)
3536
- **[OpenMP](https://www.openmp.org/)** multi-platform parallel programming (optionnal)
3637
- **[Doxygen](https://www.doxygen.nl/)** documentation generator (optional)
@@ -91,6 +92,11 @@ Under Linux (Ubuntu) distribution the installation procedure is as follows:
9192
echo 'export PATH=/home/src/micmac/MMVII/bin:$PATH' >> ~/.bashrc
9293
```
9394

95+
- Optional: to use libe57format
96+
- clone, compile and install as explained in https://github.com/asmaloney/libE57Format
97+
- in MMVII cmake, set USE_LIBE57FORMAT then fill LIBE57_INCLUDE_DIR and LIBE57_LIBRARY with the E57-install/ data.
98+
99+
94100
## Windows
95101
Under Windows the installation procedure is as follows:
96102

‎MMVII/src/ImportFormat/ImportStaticScan.cpp‎

Lines changed: 189 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "MMVII_Geom3D.h"
44
#include "../Mesh/happly.h"
55
#include <functional>
6+
#include <E57Format/E57SimpleReader.h>
67

78
/**
89
\file importStaticScan.cpp
@@ -31,7 +32,11 @@ public :
3132

3233

3334
void readPlyPoints(std::string aPlyFileName);
35+
#if (USE_LIBE57FORMAT)
36+
void readE57Points(std::string aE57FileName);
37+
#endif
3438
void convertToThetaPhiDist();
39+
void convertToXYZ();
3540
void getAnglesMinMax();
3641
void estimatePhiStep();
3742
void computeLineCol();
@@ -48,6 +53,11 @@ private :
4853
bool mNoMiss; // are every point present in cloud, even when no response?
4954

5055
// data
56+
bool mHasCartesian; // in original read data
57+
bool mHasIntensity; // in original read data
58+
bool mHasSpherical; // in original read data
59+
bool mHasRowCol; // in original read data
60+
5161
tREAL8 mDistMinToExist; // boundary to check if no response points are present
5262
std::vector<cPt3dr> mVectPtsXYZ;
5363
std::vector<tREAL8> mVectPtsIntens;
@@ -99,12 +109,26 @@ cPt3dr cart2spher(const cPt3dr & aPtCart)
99109
return {theta, phi, dist};
100110
}
101111

112+
cPt3dr spher2cart(const cPt3dr & aPtspher)
113+
{
114+
tREAL8 dhz = aPtspher.z()*cos(aPtspher.y());
115+
tREAL8 x = dhz* cos(aPtspher.x());
116+
tREAL8 y = dhz* sin(aPtspher.x());
117+
tREAL8 z = aPtspher.z()*sin(aPtspher.y());
118+
return {x, y, z};
119+
}
120+
121+
102122
using namespace happly;
103123

104124
void cAppli_ImportStaticScan::readPlyPoints(std::string aPlyFileName)
105125
{
126+
StdOut() << "Read ply file " << aPlyFileName << "..." << std::endl;
106127
mVectPtsXYZ.clear();
128+
mVectPtsTPD.clear();
107129
mVectPtsIntens.clear();
130+
mVectPtsCol.clear();
131+
mVectPtsLine.clear();
108132
try
109133
{
110134
PLYData aPlyF(aPlyFileName,false);
@@ -119,13 +143,19 @@ void cAppli_ImportStaticScan::readPlyPoints(std::string aPlyFileName)
119143
}
120144
}
121145

146+
mHasCartesian = true;
147+
mHasIntensity = false;
148+
mHasSpherical = false;
149+
mHasRowCol = false;
150+
122151
// try to fill points attribute with any "i", "I" or "*intens*" property
123152
auto aVertProps = aPlyF.getElement("vertex").getPropertyNames();
124153
auto aPropIntensityName = std::find_if(aVertProps.begin(), aVertProps.end(), [](const std::string &s){
125154
return (ToLower(s)=="i") || (ToLower(s).find("intens") != std::string::npos);
126155
});
127156
if (aPropIntensityName!= aVertProps.end())
128157
{
158+
mHasIntensity = true;
129159
mVectPtsIntens = aPlyF.getElement("vertex").getProperty<tREAL8>(*aPropIntensityName);
130160
}
131161

@@ -136,6 +166,86 @@ void cAppli_ImportStaticScan::readPlyPoints(std::string aPlyFileName)
136166
}
137167
}
138168

169+
#if (USE_LIBE57FORMAT)
170+
void cAppli_ImportStaticScan::readE57Points(std::string aE57FileName)
171+
{
172+
StdOut() << "Read e57 file " << aE57FileName << "..." << std::endl;
173+
mVectPtsXYZ.clear();
174+
mVectPtsTPD.clear();
175+
mVectPtsIntens.clear();
176+
mVectPtsCol.clear();
177+
mVectPtsLine.clear();
178+
try
179+
{
180+
e57::Reader reader( aE57FileName, {});
181+
MMVII_INTERNAL_ASSERT_tiny(reader.IsOpen(), "Error: unable to open file " + aE57FileName)
182+
StdOut() << "Image2DCount: " << reader.GetImage2DCount() << "\n";
183+
StdOut() << "Data3DCount: " << reader.GetData3DCount() << "\n";
184+
MMVII_INTERNAL_ASSERT_tiny(reader.GetData3DCount()==1, "Error: File should have exactly 1 scan for now")
185+
e57::E57Root fileHeader;
186+
reader.GetE57Root( fileHeader );
187+
/*StdOut() << fileHeader.formatName << " =? " << "ASTM E57 3D Imaging Data File" << std::endl;
188+
StdOut() << fileHeader.versionMajor << " =? " << 1 << std::endl;
189+
StdOut() << fileHeader.versionMinor << " =? " << 0 << std::endl;
190+
StdOut() << fileHeader.guid << " =? " << "Zero Points GUID" << std::endl;*/
191+
e57::Data3D data3DHeader;
192+
reader.ReadData3D( 0, data3DHeader );
193+
// data3DHeader.indexBounds is not correct
194+
const uint64_t cNumPoints = data3DHeader.pointCount;
195+
e57::Data3DPointsFloat pointsData( data3DHeader );
196+
auto vectorReader = reader.SetUpData3DPointsData( 0, cNumPoints, pointsData );
197+
const uint64_t cNumRead = vectorReader.read();
198+
MMVII_INTERNAL_ASSERT_tiny(cNumPoints==cNumRead, "Error: cNumPoints!=cNumRead")
199+
200+
mHasCartesian = pointsData.cartesianX && pointsData.cartesianY && pointsData.cartesianZ;
201+
mHasIntensity = pointsData.intensity;
202+
mHasSpherical = pointsData.sphericalAzimuth && pointsData.sphericalElevation && pointsData.sphericalRange;
203+
mHasRowCol = pointsData.columnIndex && pointsData.rowIndex;
204+
205+
if (mHasCartesian){
206+
mVectPtsXYZ.resize(cNumRead);
207+
for (uint64_t i=0;i<cNumRead;++i)
208+
mVectPtsXYZ[i] = {pointsData.cartesianX[i], pointsData.cartesianY[i], pointsData.cartesianZ[i]};
209+
}
210+
if (mHasSpherical){
211+
mVectPtsTPD.resize(cNumRead);
212+
for (uint64_t i=0;i<cNumRead;++i)
213+
mVectPtsTPD[i] = {pointsData.sphericalAzimuth[i], pointsData.sphericalElevation[i], pointsData.sphericalRange[i]};
214+
}
215+
if (mHasIntensity){
216+
mVectPtsIntens.resize(cNumRead);
217+
for (uint64_t i=0;i<cNumRead;++i)
218+
mVectPtsIntens[i] = pointsData.intensity[i];
219+
}
220+
if (mHasRowCol){
221+
mVectPtsLine.resize(cNumRead);
222+
mVectPtsCol.resize(cNumRead);
223+
mSL_data.mMaxLine = 0;
224+
for (uint64_t i=0;i<cNumRead;++i)
225+
{
226+
mVectPtsLine[i] = pointsData.rowIndex[i];
227+
if (pointsData.rowIndex[i]>mSL_data.mMaxLine)
228+
mSL_data.mMaxLine = pointsData.rowIndex[i];
229+
}
230+
mSL_data.mMaxCol = 0;
231+
for (uint64_t i=0;i<cNumRead;++i)
232+
{
233+
mVectPtsCol[i] = pointsData.columnIndex[i];
234+
if (pointsData.columnIndex[i]>mSL_data.mMaxCol)
235+
mSL_data.mMaxCol = pointsData.columnIndex[i];
236+
}
237+
}
238+
239+
vectorReader.close();
240+
241+
}
242+
catch (const std::runtime_error &e)
243+
{
244+
MMVII_UserError(eTyUEr::eReadFile, std::string("Error reading E57 file \"") + aE57FileName + "\": " + e.what());
245+
}
246+
}
247+
#endif
248+
139249
void cAppli_ImportStaticScan::convertToThetaPhiDist()
140250
{
141251
mNoMiss = false;
@@ -166,6 +276,37 @@ void cAppli_ImportStaticScan::convertToThetaPhiDist()
166276
<< mSL_data.mPhiMin << ", " << mSL_data.mPhiMax << "\n";
167277
}
168278

279+
280+
void cAppli_ImportStaticScan::convertToXYZ()
281+
{
282+
mNoMiss = false;
283+
mVectPtsXYZ.resize(mVectPtsTPD.size());
284+
size_t aNbPtsNul = 0;
285+
for (size_t i=0; i<mVectPtsTPD.size(); ++i)
286+
{
287+
mVectPtsXYZ[i] = spher2cart(mVectPtsTPD[i]);
288+
if (mVectPtsTPD[i].z()<mDistMinToExist)
289+
aNbPtsNul++;
290+
}
291+
StdOut() << aNbPtsNul << " null points\n";
292+
mNoMiss = aNbPtsNul>0;
293+
294+
// get min max theta phi
295+
cWhichMinMax<int, tREAL8> aMinMaxTheta;
296+
cWhichMinMax<int, tREAL8> aMinMaxPhi;
297+
for (const auto & aPtAng: mVectPtsTPD)
298+
{
299+
aMinMaxTheta.Add(0,aPtAng.x());
300+
aMinMaxPhi.Add(0,aPtAng.y());
301+
}
302+
mSL_data.mThetaMin = aMinMaxTheta.Min().ValExtre();
303+
mSL_data.mThetaMax = aMinMaxTheta.Max().ValExtre();
304+
mSL_data.mPhiMin = aMinMaxPhi.Min().ValExtre();
305+
mSL_data.mPhiMax = aMinMaxPhi.Max().ValExtre();
306+
StdOut() << "Box: theta " << mSL_data.mThetaMin << ", " << mSL_data.mThetaMax << " phi "
307+
<< mSL_data.mPhiMin << ", " << mSL_data.mPhiMax << "\n";
308+
}
309+
169310
void cAppli_ImportStaticScan::estimatePhiStep()
170311
{
171312
// find phi min diff
@@ -197,6 +338,8 @@ void cAppli_ImportStaticScan::estimatePhiStep()
197338

198339
void cAppli_ImportStaticScan::computeLineCol()
199340
{
341+
if (mHasRowCol)
342+
return; // nothing to do
200343
tREAL8 aColChangeDetectorInPhistep = 100;
201344

202345
// compute line and col for each point
@@ -387,51 +530,76 @@ tREAL8 cAppli_ImportStaticScan::doVerticalize()
387530
int cAppli_ImportStaticScan::Exe()
388531
{
389532
mPhProj.FinishInit();
390-
readPlyPoints(mNameFile);
533+
std::string aPostFix = ToLower(LastPostfix(mNameFile));
534+
if (aPostFix=="ply")
535+
{
536+
readPlyPoints(mNameFile);
537+
} else if (aPostFix=="e57") {
538+
#if (USE_LIBE57FORMAT)
539+
readE57Points(mNameFile);
540+
#else
541+
MMVII_INTERNAL_ASSERT_tiny(false, "Error: missing libe57format");
542+
#endif
543+
} else
544+
MMVII_INTERNAL_ASSERT_tiny(false, "Error: unknown file format for "+mNameFile);
545+
546+
StdOut() << "Read data: ";
547+
if (mHasCartesian)
548+
StdOut() << mVectPtsXYZ.size() << " cartesian points";
549+
if (mHasSpherical)
550+
StdOut() << mVectPtsTPD.size() << " spherical points";
551+
if (mHasIntensity)
552+
StdOut() << " with intensity";
553+
if (mHasRowCol)
554+
StdOut() << " with row-col";
555+
StdOut() << "\n";
556+
557+
if (mHasCartesian && !mHasSpherical)
558+
{
559+
cRotation3D<tREAL8> aRotFrame = cRotation3D<tREAL8>::RotFromCanonicalAxes(mTransfoIJK);
560+
// apply rotframe to original points
561+
for (auto & aPtXYZ : mVectPtsXYZ)
562+
{
563+
aPtXYZ = aRotFrame.Value(aPtXYZ);
564+
}
565+
convertToThetaPhiDist();
566+
} else if (!mHasCartesian && mHasSpherical) // mTransfoIJK not used if spherical
567+
{
568+
convertToXYZ();
569+
}
391570

392-
StdOut() << "Got " << mVectPtsXYZ.size() << " points.\n";
393571
MMVII_INTERNAL_ASSERT_tiny(!mVectPtsXYZ.empty(),"Error reading "+mNameFile);
394-
if (!mVectPtsIntens.empty())
572+
if (mHasIntensity)
395573
{
396-
StdOut() << "Intensity found.\n";
397574
MMVII_INTERNAL_ASSERT_tiny(mVectPtsXYZ.size()==mVectPtsIntens.size(),"Error reading "+mNameFile);
398575
}
399576

400-
StdOut() << "Sample:\n";
577+
StdOut() << "Cartesian sample:\n";
401578
for (size_t i=0; (i<10)&&(i<mVectPtsXYZ.size()); ++i)
402579
{
403580
StdOut() << mVectPtsXYZ.at(i);
404-
if (!mVectPtsIntens.empty())
581+
if (mHasIntensity)
405582
StdOut() << " " << mVectPtsIntens.at(i);
406583
StdOut() << "\n";
407584
}
408585
StdOut() << "...\n";
409586

410-
cRotation3D<tREAL8> aRotFrame = cRotation3D<tREAL8>::RotFromCanonicalAxes(mTransfoIJK);
411-
// apply rotframe to original points
412-
for (auto & aPtXYZ : mVectPtsXYZ)
413-
{
414-
aPtXYZ = aRotFrame.Value(aPtXYZ);
415-
}
416-
417-
convertToThetaPhiDist();
418-
419587
// check theta-phi :
420588
StdOut() << "Spherical sample:\n";
421589
for (size_t i=0; (i<10)&&(i<mVectPtsTPD.size()); ++i)
422590
{
423591
StdOut() << mVectPtsTPD[i];
424592
StdOut() << "\n";
425593
}
426-
StdOut() << "...\n";
594+
StdOut() << "..." << std::endl;
427595

428596
estimatePhiStep();
429-
430597
computeLineCol();
431598

432599
//fill rasters
433600
/*fillRaster<tU_INT1>("totoMask.png", [this](int i){auto aPtAng = mVectPtsTPD[i];return (aPtAng.z()<mDistMinToExist)?0:255;} );
434-
fillRaster<tU_INT1>("totoIntens.png", [this](int i){return mVectPtsIntens[i]*255;} );
601+
if (mHasIntensity)
602+
fillRaster<tU_INT1>("totoIntens.png", [this](int i){return mVectPtsIntens[i]*255;} );
435603
fillRaster<tU_INT4>("totoIndex.png", [](int i){return i;} );
436604
fillRaster<float>("totoDist.tif", [this](int i){auto aPtAng = mVectPtsTPD[i];return aPtAng.z();} );
437605
fillRaster<float>("totoTheta.tif", [this](int i){auto aPtAng = mVectPtsTPD[i];return aPtAng.x();} );
@@ -444,7 +612,7 @@ int cAppli_ImportStaticScan::Exe()
444612
for (size_t i=0; (i<10)&&(i<mVectPtsXYZ.size()); ++i)
445613
{
446614
StdOut() << mVectPtsXYZ.at(i);
447-
if (!mVectPtsIntens.empty())
615+
if (mHasIntensity)
448616
StdOut() << " " << mVectPtsIntens.at(i);
449617
StdOut() << "\n";
450618
}
@@ -490,7 +658,8 @@ int cAppli_ImportStaticScan::Exe()
490658
mSL_data.mRasterMask = mPhProj.DPStaticLidar().FullDirOut() + mSL_data.mStationName + "_" + mSL_data.mScanName + "_mask.tif";
491659

492660
fillRaster<tU_INT1>(mSL_data.mRasterMask, [this](int i){auto aPtAng = mVectPtsTPD[i];return (aPtAng.z()<mDistMinToExist)?0:255;} );
493-
fillRaster<tU_INT1>( mSL_data.mRasterIntensity, [this](int i){return mVectPtsIntens[i]*255;} );
661+
if (mHasIntensity)
662+
fillRaster<tU_INT1>( mSL_data.mRasterIntensity, [this](int i){return mVectPtsIntens[i]*255;} );
494663
fillRaster<float>(mSL_data.mRasterDistance, [this](int i){auto aPtAng = mVectPtsTPD[i];return aPtAng.z();} );
495664

496665
fillRaster<tU_INT4>("titiIndex.png", [](int i){return i;} );

0 commit comments

Comments
 (0)