Implicit Objects and the Marching Cubes Algorithm (implementation in C++)

Hello,

today I am going to talk about algebraic representation of Objects and the Marching Cubes Algorithm. This tutorial's code is available at:
https://github.com/Art-n-MathS/ImplicitObjects

The code is part of the open source software DASOS, which is able to reconstruct forest from full-waveform LiDAR data. For more information about DASOS please look at: http://miltomiltiadou.blogspot.co.uk/2015/03/las13vis.html

Please note that some of the explanation were copied from my paper:
Miltiadou, M, Warren M.A, Grant M., Brown M., 2015, Alignment of Hyperspectral Imagery and full-waveform LiDAR data for Visualisation and Classification purposes, 36th International Symposium of Remote Sensing of the Enviroment, ISPRS Archives
Available at: http://www.int-arch-photogramm-remote-sens-spatial-inf-sci.net/XL-7-W3/1257/2015/isprsarchives-XL-7-W3-1257-2015.pdf

Implicit Representation of objects

Numberical implicitization was introduced by Blinn in 1982. Numerical implicitization allows definition of complex objects, which can easily be modified, without involving large number of triangles.

A function f(X) defines an object, when every n-dimensional point X the following conditions apply:

f(X) = a, when X lies on the surface of the object
f(X) > a, when X lies inside the object
f(X) < a, when X lies outside the object

where,
X = an n-dimensional point. It usually it is in Euclidean space defining the x, y, z poisitions of the point
f(X) = the function that defines the object of our interest. In the coding example is the function of the sphere
a = the isosurface of the object, which defines the boundaries of the object. f(X) is equal to a iff X lies on the surface of the object.

Marching Cubes Algorithm:

Even though numerical implicitization is beneficial for reducing storage memory and for various resolution renderings of the same object, visualising numerical/algebraic objects is not straight forward, since they contain no discrete values. This problem can be addressed either by ray-tracing or by polygonisation. In 1983, Hanraham suggested a ray-tracing approach, which used the equation derived from the ray and the surface intersection to depict the algebraic object into an image. The Marching Cubes algorithm was later introduced for polygonising implicit objects.

The Marching cubes algorithm constructs surfaces from implicit object using a search table. Assume that f(X) defiens an algebraic object. At first the space is divided into cubes, named voxels.Each voxel is defined by eight corner points, which lie either inside of outside the object. By enumerating all the possible cases and linearly interpolating the intersections along the edges, the surface of the implicit object is constructed (Lorensen and Cline, 1987).

In 2021 a research was published that investigates the performance of 6 different data structures for extracting an iso-surface (creating a 3D polygonal model) from voxelised data (Miltiadou et, al. 2021). It was shown that Integral Volumes was the fastest from the six, but it consumes the most memory. Each data structure has its owns prons and cons. Additionally, in that paper a new category of data strucutre "Integral Trees" were introduced. 


For more information about the Marching Cubes algorithm please look at the following blogpost:
http://paulbourke.net/geometry/polygonise/

Example Output:

The output .obj file, visualised in Meshlab

Code

This tutorial's code is available at: https://github.com/Art-n-MathS/ImplicitObjects

In this example, a sphere is defined using an function, polygonised using the Marching Cubes Algorithm and then exported into an .obj file.

Requirements:
- gmtl library
- c++11
- qmake

Compile:
$: qmake
$: make clean
$: make

Run:
$: ./ImplicitObjects

References

Miltiadou, M.; Campbell, N.D.F.; Cosker, D.; Grant, M.G. A Comparative Study about Data Structures Used for Efficient Management of Voxelised Full-Waveform Airborne LiDAR Data during 3D Polygonal Model Creation. Remote Sens. 2021, 13, 559. https://doi.org/10.3390/rs13040559

Blinn, J.F, 1982. A Generalization of Algebraic Surface Drawing. ACM Trans.Graph, pp. 235-256

Hanrahan, P., 1983. Ray tracing algebraic surfaces. ACM SIGGRAPH Computer Graphics, Vol 17, No. 3.

Lorense, W. E., & Cline, H.E., 1987. Marching Cubes: A high resolution 3D surface construction algorithm ACM Siggraph Computer Graphics, Vol 21, No 4 pp. 163-169


Miltiadou, M, Warren M.A, Grant M., Brown M., 2015, Alignment of Hyperspectral Imagery and full-waveform LiDAR data for Visualisation and Classification purposes, 36th International Symposium of Remote Sensing of the Enviroment, ISPRS Archives

DOI: 10.5194/isprsarchives-XL-7-W3-1257-2015

Parsing a file into a words (C++ Programming)

It usually seems difficult to parse files, so I found this way which you can parse a file into words within a few lines and it works perfect. So keeping it here makes it more accessible when I will next need it. ;p

Let's assume that you have the following file that need to be parsed:

My name is Milto
and I like Capoeira! :)

Then the following code will read the file and print each word at separated line.

#include < iostream >
#include < fstream >
#include < iterator >
#include < string >
#include < vector >

int main(void)
{
   std::string filename = "inputt.xt";
   std::ifstream mystream(filename.c_str());
   if(!mystream)
   {
      std::cerr << "File \"" << filename << "\" not found.\n";
   }
   std::istream_iterator< std::string > it(mystream);
   std::istream_iterator< std::string > sentinel;
   
   std::vector < std::string > words(it,sentinel);

   for(unsigned int i=0; i< words.size();++i)
   {
      std::cout << words[i] << "\n";
   }
   
   return 0;
}

The output result will be the following:

$: g++ parsing.cpp -o out
$: ./out
My
name
is
Milto
and
I
like
Capoeira!
:)

You may download this example code from here:
https://www.dropbox.com/s/2chkdvjf2m0uhcp/parseFileIntoWords.zip?dl=0

DASOS - Open Source Software for processing full-waveform LiDAR and hyperspectral Images

What is DASOS?

DASOS is an open source software that aims to ease the way of handling the full-waveform LiDAR data. Its name was derived from from the Greek word δάσος (=forest). For any publication using the software please cite the following paper: Open source software DASOS: efficient accumulation, analysis, and visualisation of full-waveform lidar


The aim of this software is to enhance the visualisations and classifications of forested areas using coincident full-waveform (fw) LiDAR data and hyperspectral images. It uses either full-waveform LiDAR only or both datasets to generate metrics understandable for foresters and 3D virtual models.

Influenced by Persson et al (2005), voxelisation is an integral part of DASOS. The intensity profile of each full-waveform pulse is accumulated into a voxel array, building up a 3D density volume. The correlation between multiple pulses into a voxel representation produces a more accurate representation, which confers greater noise resistance and it further opens up possibilities of vertical interpretation of the data. The 3D density volume is then aligned with the hyperspectral images using a 2D grid similar to Warren et al (2014) and both datasets are used in visualisations and classifications.

There are three main functionalities of DASOS:

• the generation of 2D metrics aligned with hyperspectral Images 
• construction of 3D polygonal meshes and
• the characterisation of objects using feature vectors. 

Aligned Metrics

Here is a list of the available metrics:

From FW LiDAR (LAS1.3 or Pulswaves file formats):
- Non-Empty Voxels
- Density
- Thickness
- First Patch
- Last Patch
- Lowest Return
- Average Height Difference (works as an edge detection algorithm)
- AGC intensity


A visual explanation of the available full-waveform LiDAR metrics is given below:

There are also the following Hyperspectral metrics (derived from .bil files & .igm files)
- Hyperspectral Mean
- NDVI
- A single hyperspectral band

Polygonal Meshes

The following example is from New Forest and generated at one meter resolution.The first one was generated using FW LiDAR data and three user defined bands of the hyperspectral images. The second one was generated with only the FW LiDAR data.

The following video was also rendered in Maya using a polygon exported from DASOS:

https://www.youtube.com/watch?v=ZsnB9822xmU&feature=youtu.be 

List of Feature Vectors

This is useful for characterising object inside the 3D space (e.g. trees). For each column of the voxelised FW LiDAR, information around its local area are exported. The exported format is .csv to easy usage in common statistical packages like R and matlab.

In the following images there are two output examples. The top on contains processed information about the data and the second file contains the raw intensity values. Additionally each line is a feature vector. 

Related Links

The full user guide is available at:
https://github.com/Art-n-MathS/DASOS/blob/master/DASOS_userGuide_v2.pdf

The windows executable and the code are available at:
https://github.com/Art-n-MathS/DASOS

To add your own customised 2D metrics, you may read the following tutorial:
http://miltomiltiadou.blogspot.co.uk/2016/07/how-to-add-metrics-to-dasos.html

For any questions regarding the usage of the software please use the following Google Group:
https://groups.google.com/forum/#!forum/dasos---the-native-full-waveform-fw-lidar-software

Updates about DASOS could be found @_DASOS_ on twitter

For more information about the algorithms please refer to the related paper:
https://www.researchgate.net/publication/334069759_Open_source_software_DASOS_efficient_accumulation_analysis_and_visualisation_of_full-waveform_lidar

Related Papers:

- Open source software DASOS: efficient accumulation, analysis, and visualisation of full-waveform lidar

-  Alignment of Hyperspectral Imagery and full-waveform LiDAR data for visualisation and classification purposes

- Detection of dead standing Eucalyptus camaldulensis without tree delineation for managing biodiversity in native Australian forest.

- Detecting Dead Standing Eucalypt Trees from Voxelised Full-Waveform Lidar Using Multi-Scale 3D-Windows for Tackling Height and Size Variations

- Improving and optimising visualisations of full-waveform LiDAR data

- Reconstruction of a 3D Polygon Representation from full-waveform LiDAR data

- A Comparative Study about Data Structures Used for EfficientManagement of Voxelised Full-Waveform Airborne LiDARData during 3D Polygonal Model Creation

Acknowledgements

Special thanks are given to my industrial supervisor Dr. Michael Grant, who has supported me during my entire studies.

I would also like to thanks my initial and current supervisors Dr. Matthew Brown and Dr. Neill D.F Cambpell and all the people who occasionally got involved: Dr. Mark Warren, Susana Gonzalez Aracil, Dr.  Daniel Clewley and Dr. Darren Cosker.

This project is funded by the Centre for Digital Entertainement and Plymouth Marine Laboratory.

The data used were collected by the NERC Airborne Research and Survey Facility (ARSF). Copyright is held by the UK Natural Environment Research Council (NERC).


Please note that this text was taken and modified from the EngD thesis of Milto Miltiadou, which was submitted to University of Bath in 2017.

LAS1.3 file Reader in C++

Hello,

today I am going to talk about a part of my code I wrote for my research. Well, my research is about visualising and classifying Airborne Remote Sensing data. About a year ago, I was provided a python script (https://github.com/pmlrsg/arsf_tools) that reads the LAS file but working in python was very very slow, so eventually I had to write my own LAS1.3 reader in C++.

The code for reading a LAS1.3 file has now been released as an open source code under the GNU General Public Licence, version 3, and it is available at:
https://github.com/Art-n-MathS/LAS1.3Reader.

In this blogpost I am planning to explain how you can use the code and also how the code works.

Please note that this work is supported by the Centre for DIgital Entertainment at the University of Bath, and Plymouth Marine Laboratory. 

1. Data Specifications

The data that I used to test the program were provided by ARSF NERC and they were collected using a small footprint Leica ALS50-II LiDAR system.

The file specifications of LAS1.3 files were published in 2010 and they are available here: http://www.asprs.org/a/society/committees/standards/LAS_1_3_r11.pdf

The point data record format is assumed to always be  of format 4 and it is also assume that waveform packets exists. Otherwise the reader prints an error and terminates.

2. Compile and run the Program

Before compiling the program, make sure that c++11 and gmtl libraries are installed on your computer.

Further in the main.cpp file replace the <DIRECTORY OF LAS FILE> with the name of the LAS.13 file of your interest. Unfortunately due to license restrictions I couldn't provide a sample file. But I may get to do it later.

Compile the program using the makefile:

make

and run the progrma:

./LASReader

The output is the following one. It first print the related pulse information, then all the waveform samples and at the end the associated discrete values:

1199653 waveforms found
1511062 discrete points found
There are 2477840 Discrete Without Waveforms
----------------------------------------------------------
the pulse manager has : 1199653 pulses
Point                            4.3787e+05 1.0486e+05 16.316
Return Number                    

Number of returns for this pulse 

Time                             3.8836e+05
Scan Angle                       
Classification                   

Temporal Sample Spacing          2
AGC gain                         
Digitiser Gain                   0.017291
Digitiser Offset                 0
No. of Samples                   256
Sample Length                    0.29979
Return Point Location            21.816
Point in Waveform                3.2702
Offset                           0.021667 0.0059378 0.29887
Origin                           4.3787e+05 1.0486e+05 19.576
Waveform Samples: ( x , y , z , I )
( 4.3787e+05 , 1.0486e+05 , 19.576 , 15 )
( 4.3787e+05 , 1.0486e+05 , 19.277 , 16 )
( 4.3787e+05 , 1.0486e+05 , 18.978 , 14 )
( 4.3787e+05 , 1.0486e+05 , 18.679 , 14 )
( 4.3787e+05 , 1.0486e+05 , 18.381 , 14 )
( 4.3787e+05 , 1.0486e+05 , 18.082 , 14 )
...
...
( 4.3787e+05 , 1.0486e+05 , -55.739 , 14 )
( 4.3787e+05 , 1.0486e+05 , -56.038 , 14 )
( 4.3787e+05 , 1.0486e+05 , -56.337 , 13 )
( 4.3787e+05 , 1.0486e+05 , -56.636 , 14 )

Associated discrete points (x , y  , z , I)
( 4.3787e+05 , 1.0486e+05 , 16.316 , 141

3. Classes Information

There are three classes in the LAS reader: Las1_3_handler, PulseManager and Pulse.The Las1_3_handler takes as input a LAS1.3 file, read the binary file according to the specifications, saves all the information read inside a PulseManager and returns the PulseManager. (For more information on how to read binary files please look at: http://miltomiltiadou.blogspot.co.uk/2013/09/read-binary-files-into-structures-in-c.html)

The following code shows how to use the Las1_3_handler class. Please note that the Pulse Manager is dynamic allocated and the user is responsible for realising the memory afterwards.

Las1_3_handler lasHandler("/local1/data/scratch/mmi/2010_098_FW/classified_manually/LDR-FW10_01-201009822.LAS");
PulseManager *p = lasHandler.readFileAndGetPulseManager();
// do stuff with p
delete p;

The PulseManager contains and manages multiple Pulses. Each Pulse contains the point data record information, the associated waveform packet and all the discrete returns, which are the peak reflectances of that Pulse. Nevertheless, due to the speed of the flight most of the time waveforms are only saved for about half of the pulses. The position and intensity of the discrete returns with no waveform associated are saved into m_discretePoints and m_discreteIntensities inside the PulseManager. So far the only thing you can do is to get the number of Pulses and print all the information of a Pulse of your interest.

std::cout << "the pulse manager has : " << p->getNumOfPulses() << " pulses\n";
p->printPulseInfo(104510);

An example of output is shown in Section 2. Further if the input pulse number doesn't exist a error message is printed instead.

4. Future Work

To extend the code and add more stuff is considerable easy. The only thing that need to be pointed out is the way of calculating the position of the waveform samples. To reduce the amount of memory used only the origin and offset of the waveform samples are saved. So, you have to calculate the position of each sample before used and this can be done inside the Pulse class as follow:

gmtl::Vec3f tempPosition = m_origin;
for(unsigned short int i=0; i< m_noOfSamples; ++i)
{
   std::cout << "( " << tempPosition[0] << " , " << tempPosition[1] << " , "
             << tempPosition[2] << " , " <<  (int) m_returns[i] <<  " )\n";
   tempPosition-=m_offset;
}

More information about coding are given in the Doxygen Documentation. Also feel free to contact me if you have further questions.

I would also like to give you an insignt of my future posts and publications with the following demo of visualising fw LiDAR data with hyperspectral iamges:

The approach of generating polygon meshes from full-waveform LiDAR data has been presented at RSPSoc Conference 2014 and if you are interested you can either download the extended abstract from here: http://rspsoc2014.co.uk/, or wait for one of my following blogpost where I will explain that in more details. Please don't forget to reference if you use any of this work.

Quick Point Cloud Search in C++ using unordered_multimap

 The main idea of the algorithm as explained to one of my previous post is:
- At first, the points are imported into a grid.
- Then, for every point of our interest
    - find the square the point lies in
    - loop through all the points that lie inside that square only
    - return the closest point

Here I implemented a 2D points cloud search for aligning the hyperspectral images with the LiDAR data[1] in the open source software DASOS [2], but it can easily be extended to work for 3D.

At first we need to find the limits of the points cloud. So we need to loop through all the points m_points, and find the minimum and maximum values.

 // Loop through all the points and find the minimum and max
if(i_pointsCloud.size() < 1)
{
   std::cout << "Warning: Points cloud too small\n";
   m_min.m_x = 0.0f; m_min.m_y = 0.0f;
   m_max.m_y = 0.1f; m_max.m_y = 0.1f;
   return;
}
m_min.m_x = m_points[0].m_x;
m_min.m_y = m_points[0].m_y;
m_max.m_x = m_points[0].m_x;
m_max.m_y = m_points[0].m_y;
for(unsigned int i=1; i < m_points.size(); ++i)
{
   if (m_min.m_x > m_points[i].m_x)
   {
      m_min.m_x=m_points[i].m_x;
   } else if (m_max.m_x < m_points[i].m_x)
   {
      m_max.m_x=m_points[i].m_x;
   }
   if (m_min.m_y > m_points[i].m_y)
   {
      m_min.m_y=m_points[i].m_y;
   } else if (m_max.m_y < m_points[i].m_y)
   {
      m_max.m_x=m_points[i].m_y;
   }
}

We also need to calculate how many squares our grid will have in each dimension. In the example below i_samplingRate is the average number of expected points per square and m_nX, m_nY is the number of squares in the x,y axes respectively.

float rate = (m_max.m_x-m_min.m_x)/(m_max.m_y-m_min.m_y);
float numOfSquares = m_points.size() / i_samplingRate;
m_nY = sqrt(numOfSquares/rate);
m_nX = ceil(numOfSquares/m_nY);

So, after we calculate the number of square per axis, how do we facilitate unordered_multimap for quick points cloud search?
Each element of an unordered_multimap has a key associated with its value. The key value of each square (x,y) on the grid will be equal to x + y*m_nX. All the points of the points cloud are saved into a 1D array (m_points), so the value of the element is the index of the point of our interest.

for (unsigned int i=0; i < m_points.size(); ++i)
{
   unsigned int sqX = double((m_points[i].m_x-m_min.m_x)/
                             (m_max.m_x-m_min.m_x))*m_nX;
   unsigned int sqY = double((m_points[i].m_y-m_min.m_y)/
                             (m_max.m_y-m_min.m_y))*m_nY;
   m_map.emplace(getKeyOfSquare(sqX,sqY),i);
}

So far, we have a grid with all the points associated to its squares. So, the next step is to be able to find the square (x,y) that a given point(i_point) lies inside.

const unsigned int x = float((i_point.m_x-m_min.m_x)/
                             (m_max.m_x-m_min.m_x)*(float)m_nX);
const unsigned int y = float((i_point.m_y-m_min.m_y)/
                             (m_max.m_y-m_min.m_y)*(float)m_nY);

Then we can get the key of that square as explained above (key = x+y*m_nX) and using that key we can query the unordered_multimap and get all the indices of the points that lie inside that square. Then by looping through all the those point, the closest point to i_point can be found.

unsigned int closestPointIndex = 0;
float minDis = sqrt(pow(m_points[0].m_x-i_point.m_x,2.0)+
                    pow(m_points[0].m_y-i_point.m_y,2.0));
auto itsElements = m_map.equal_range(i_point.m_x + i_point.m_y * m_nX);
for (auto it = itsElements.first; it != itsElements.second; ++it)
{
   const float distance = sqrt(pow(m_points[it->second].m_x-i_point.m_x,2.0)+
                                pow(m_points[it->second].m_y-i_point.m_y,2.0));
   if(distancesecond;
      minDis = distance;
   }
}
closestpoint = m_points[closestPointIndex];

You may download the full version of the program here:
https://github.com/Art-n-MathS/QuickPointsCloudSearch

Please note that Unordered_multimap is a part of C++2011. So, to compile and run the program you need to add the tag -std=c++11 as follow:

g++ -std=c++11 main.cpp GridSearch.cpp -o myProgram
./myProgram

The output of the program is the following:

All the points of the point cloud are:

(0,0) (0,1) (0,2) (0,3) (0,4) (0,5) (0,6) (0,7) (0,8) (0,9) (1,0) (1,1) (1,2) 
(1,3) (1,4) (1,5) (1,6) (1,7) (1,8) (1,9) (2,0) (2,1) (2,2) (2,3) (2,4) (2,5) 
(2,6) (2,7) (2,8) (2,9) (3,0) (3,1) (3,2) (3,3) (3,4) (3,5) (3,6) (3,7) (3,8) 
(3,9) (4,0) (4,1) (4,2) (4,3) (4,4) (4,5) (4,6) (4,7) (4,8) (4,9) (5,0) (5,1) 
(5,2) (5,3) (5,4) (5,5) (5,6) (5,7) (5,8) (5,9) (6,0) (6,1) (6,2) (6,3) (6,4) 
(6,5) (6,6) (6,7) (6,8) (6,9) (7,0) (7,1) (7,2) (7,3) (7,4) (7,5) (7,6) (7,7) 
(7,8) (7,9) (8,0) (8,1) (8,2) (8,3) (8,4) (8,5) (8,6) (8,7) (8,8) (8,9) (9,0) 
(9,1) (9,2) (9,3) (9,4) (9,5) (9,6) (9,7) (9,8) (9,9)

The closest point to (0,0.2) is (0,0)

The closest point to (5.2,5.9) is (5,6)

*** Note ***
There is a chance to miss the closest point and that may happen when the point of interest is very close to the edges of the cell that its contained. This can be resolved by checking whether the distance between the point of interest and its closest point that within the same cell is smaller than the distance between the point of interest and the edges of its cell.

PS: if you are interested on how unordered_multimap works in general, please have a look at one of my previous posts, which focuses on that: http://miltomiltiadou.blogspot.co.uk/2014/03/unorderedmultimap-for-quick-search-in.html

Reference:

[1] Miltiadou, M., Warren, M. A., Grant, M., & Brown, M. (2015). Alignment of hyperspectral imagery and full-waveform LiDAR data for visualisation and classification purposes. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences-ISPRS Archives.

DOI: 10.5194/isprsarchives-XL-7-W3-1257-2015

[2] Miltiadou, M., Grant, M. G., Campbell, N. D., Warren, M., Clewley, D., & Hadjimitsis, D. G. (2019, June). Open source software DASOS: Efficient accumulation, analysis, and visualisation of full-waveform lidar. In Seventh International Conference on Remote Sensing and Geoinformation of the Environment (RSCy2019) (Vol. 11174, p. 111741M). International Society for Optics and Photonics.

DOI: 10.1117/12.2537915

Pass an array pointer as an input to a function C++.

That's a really short post, but that's the second time I find it useful so I decided to keep it on my blog for future reference.

So let's assume that you have pointer to an array and you would like to use this array into a function. Here is an example of how to correctly send the pointers to avoid memory leaks:

#include < iostream >

void modifyData(unsigned int **i_array,unsigned int short i_len)
{
    (*i_array)[0] = 1;
    (*i_array)[2] = 2;
}

int main(int /*argc*/, char **/*argv*/)
{
  unsigned int short len = 5;
  unsigned int *mydata = new unsigned int[len];
  // initialise the values of the array
  for(unsigned int i=0; i < len; ++i)
  {
      mydata[i] = 0;
  }
  // print the original values of the array
  for(unsigned int i = 0; i < len; ++i)
  {
      std::cout << mydata[i] << " ";
  }
  std::cout << "\n";
  // call function that will modify the data
  modifyData(&mydata,len);
  // print the new values of the array
  for(unsigned int i = 0; i < len; ++i)
  {
      std::cout << mydata[i] << " ";
  }
  std::cout << "\n";
  delete []mydata;

  return 0;
}

Unordered_multimap

Hash tables are structures that allow you to search elements using a key. In my recent publication, I am using the data structure that was named "Voxel Hashing" to test the efficiency of hash table while voxelizing/interpreting full-waveform LiDAR data for 3D polygonal mesh creation [1]. This approach is the one selected to be used on the open source software DASOS [2] since it does not store empty voxels and even though it was not the fastest in the processes of 3D polygonal mesh creation, it is able to find the value of a voxel at constant time O(n) using the Hash function. The C++ implementation of the unordered_multimap is used and this blog post explains how to use it. 

Unordered_multimap:

As mentioned before, it is used for quick search and it allows you to save more than one values with  the same key, while unordered_map doesn't. How does it work? Simply, for every key there is a bucket where all the values associated with that key are saved. 

So, unordered_multimap is in the same header file with unordered_map, so the header file that need to be included is:

#include < unordered_map >

To create a map you need to define the type of the two components of the map. For example:

std::unordered_multimap < unsigned int , std::string > myMap;

if you would like to give initial values in the map you can do it as follow:

// this gives the initial values to the constructor
std::unordered_multimap < unsigned int , std::string > myMap ({{100,"a"},{120,"b"}});

// while the following on initialises the map and then adds the values
std::unordered_multimap < unsigned int , std::string > myMap ={{100,"a"},{120,"b"}};

If you would like to add items you can do it using the command "emplace":

// add more elements
myMap.emplace(100,"c");
myMap.emplace(100,"d");
myMap.emplace(200,"e");

You may like to loop through all its elements as follow:

//loop through all its elements
std::cout << "These are all the elements saved into the map:\n";
for (auto& x: myMap)
{
   std::cout << "(" << x.first << " , " << x.second << ")\n";
}

And in my opinion the most important feature is to be able to loop through all the elements with the same key. Here is an example of how you may do it:

// print all the elements with Key = 100
 unsigned int key = 100;
 std::cout << "These are all the elements with key value: " << key << "\n";
 auto itsElements = myMap.equal_range(key);
 for (auto it = itsElements.first; it != itsElements.second; ++it)
 {
    std::cout << "(" << it->first << " , " << it->second << ")\n";
 }

By the end, another small useful feature is the ability to query the unordered_multimap and get the number of elements that exists with the same key. This is achieved using the command count() as follow:

// counts how many entries exist with the same key
std::cout << "There are " << myMap.count(key) << " entries with key equal to " << key << "\n";

To sum up here is the entire program, which defines an unordered_multimap with some initial values, adds 3 more elements, prints all the elements inside the map and finally prints all the elements with key 100.

#include < iostream >
#include < unordered_map >

int main(int /*argc*/, char **/*argv*/)
{
   // define and initialise map
   std::unordered_multimap < unsigned int , std::string > myMap ({{100,"a"},{120,"b"}});

   // add more elements
   myMap.emplace(100,"c");
   myMap.emplace(100,"d");
   myMap.emplace(200,"e");

   //loop through all its elements
   std::cout << "These are all the elements saved into the map:\n";
   for (auto& x: myMap)
   {
      std::cout << "(" << x.first << " , " << x.second << ")\n";
   }

   // print all the elements with Key = 100
   unsigned int key = 100;
   std::cout << "These are all the elements with key value: " << key << "\n";
   auto itsElements = myMap.equal_range(key);
   for (auto it = itsElements.first; it != itsElements.second; ++it)
   {
       std::cout << "(" << it->first << " , " << it->second << ")\n";
   }

   // counts how many entries exist with the same key
   std::cout << "There are " << myMap.count(key) << " entries with key equal to " << key << "\n";


   return 0;
}

And the output of the program is the following:

These are all the elements saved into the map:
(200 , e)
(100 , d)
(100 , c)
(100 , a)
(120 , b)
These are all the elements with key value: 100
(100 , d)
(100 , c)
(100 , a)
There are 3 entries with key equal to 100

References

[1] Miltiadou, M.; Campbell, N.D.F.; Cosker, D.; Grant, M.G. A Comparative Study about Data Structures Used for Efficient Management of Voxelised Full-Waveform Airborne LiDAR Data during 3D Polygonal Model Creation. Remote Sens. 2021, 13, 559. https://doi.org/10.3390/rs13040559

[2] Miltiadou, M., Grant, M. G., Campbell, N. D., Warren, M., Clewley, D., & Hadjimitsis, D. G. (2019, June). Open source software DASOS: Efficient accumulation, analysis, and visualisation of full-waveform lidar. In Seventh International Conference on Remote Sensing and Geoinformation of the Environment (RSCy2019) (Vol. 11174, p. 111741M). International Society for Optics and Photonics. 

Link: https://www.researchgate.net/publication/334069759_Open_source_software_DASOS_efficient_accumulation_analysis_and_visualisation_of_full-waveform_lidar