Statistics in C++ (Mean, Median and Standard Deviation)

Just a few simple statistics. I found this function extremely useful for my software DASOS. I used it for calculating features that may characterise trees. :)

#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include <cmath>

typedef struct Statistics {
   double mean;
   double median;
   double stdev;
} Statistics;

Statistics getMeanMedianStd(
        const std::vector<  double > &i_vector
        )
{
   Statistics stats;
   std::vector <  double  > vector(i_vector);
   std::sort(vector.begin(),vector.end());
   double sumD = std::accumulate(vector.begin(),vector.end(),0.0);
   stats.mean = sumD/double(vector.size());
   std::vector < double  > Diff(vector.size());
   std::transform(vector.begin(),vector.end(),Diff.begin(),
                  std::bind2nd(std::minus< double  > (),stats.mean));
   stats.stdev  = std::inner_product(Diff.begin(),Diff.end(),
                                      Diff.begin(),0.0);
   stats.median = vector[std::floor(double(vector.size())/2.0)];
   stats.stdev = std::sqrt(stats.stdev/vector.size());
   return stats;
}

int main(int argc, char *argv[])
{
    std::vector < double  > values({10.0,12.0,11.0,4.5,10,13,22,2,1});
       for(unsigned int i=0; i < values.size(); ++i)
       {
          std::cout << values[i] << " " ;
       }
       std::cout << "\n";
       Statistics stats = getMeanMedianStd(values);
       std::cout << "Mean               = " << stats.mean   << "\n"
                 << "Median             = " << stats.median << "\n"
                 << "Standard Deviation = " << stats.stdev  << "\n";
    std::cout << "   ***   EXIT   ***\n";
    return EXIT_SUCCESS;
}

PS: for the median if the array contains an even number of values, then the bigger one is chosen between the two of them.

How to add metrics to DASOS

DASOS is our open source software for managing full-waveform (FW) LiDAR data [1] (http://miltomiltiadou.blogspot.co.uk/2015/03/las13vis.html).  Nevertheless it has a limited number of metrics, but it's design make it easy for people to add their own. So this blogpost aims to explain how to add your own metrics derived from the voxelised FW LiDAR data. It is advised to forward me (mmiltoo(at)gmail(dot)com) the new metrics classes in order to add the to future releases of DASOS and your name will also be added to the contributors.

There are two steps for generating your own metrics:
Step 1. Create a new Class for your new metric. To avoid confusion let's name our new class Metric. Please save the Metric.h file into the director ./include/Maps and the Metric.cpp file into the ./src/Maps directory to keep the files tidy. The Metrics class will inherit from the the base class Map. The header file (Metric.h) is standard and it should always be as the following (please replace all METRIC/Metric with the name of your actual new metric):

#ifndef METRIC_H
#define METRIC_H
#include "Map.h"
//-------------------------------------------------------------------------
/// @file Metric.h
/// @author <your name>
/// @version 1.0
/// @date <date generated>
/// @class Metric
/// @brief
//-------------------------------------------------------------------------


class Metric: public Map
{
public:
    //-------------------------------------------------------------------------
    /// @brief default constructor
    //-------------------------------------------------------------------------
    Metric(
            const std::string i_name,
            Volume *i_obj
            );
    //-------------------------------------------------------------------------
    /// @brief default destructor
    //-------------------------------------------------------------------------
    ~Metric();

private:
    //-------------------------------------------------------------------------
    /// @brief method that creates the Map
    //-------------------------------------------------------------------------
    void createMap();
};

#endif // METRIC_H

The only method that needs to be implemented is the createMap() which is a compulsory virtual function. The .cpp file should be as follow:

#include "Metric.h"

//-----------------------------------------------------------------------------
Metric::Metric(
        const std::string i_name,
        Volume *i_obj
        ):
    Map(i_name,i_obj)
{
}


//-----------------------------------------------------------------------------
void Metric::createMap()
{
   // Loop through all the voxels and generate the metric of interest
   // the variable m_noOfPixelsX, m_noOfPixelsY and m_noOfPixelsZ gives you 
   // the the number of voxels in x,y,z axis respectively
   for(unsigned int x=0; x<m_noOfPixelsX; ++x)
   {
      for(unsigned int y=0; y<m_noOfPixelsY; ++y)
      {
         for(unsigned int z=0; z<m_noOfPixelsZ; ++z)
         {
            // in m_mapValues all the values of the metrics are saved
            // the following command assigns the value 0 at the (x,y) position 
            // of the map
            m_mapValues[getIndex(x,y)]=-0.0f;

            // the voxelised 3D volume is the m_object variable and some 
            // useful examples of using it are shown below:

            // get the length of the voxel
            float voxelLength = m_object->getVoxelLen(); 
            // get the intesity at voxel (x,y,z)
            float intensity = m_object->getIntensity(x,y,z);
            // check whether the value of the voxel at (x,y,z) is considered to 
            // be inside or outside the scanned object. This checks whether the 
            // intensity value is above the boundary threshold. 
            bool isInside = m_object->isInside(x,y,z);            
         }
      }
   }
}

//-----------------------------------------------------------------------------
Metric::~Metric()
{}

Step 2: Link the new Metric Class with the rest of the program. This is done by modifying the .cpp file of MapsManager class.  Here it is shown the 4 additions that needs to be done in order to link your new metric with the rest of the program.

#include "MapsManager.h"
#include "ThicknessMap.h"
#include "FirstPatch.h"
#include "LastPatch.h"
// 1st ADDITION: include the header file of the new class
#include "Metric.h
// end of 1st ADDITION
#include <map>
#include <algorithm>

//-----------------------------------------------------------------------------
MapsManager::MapsManager():m_map(0),
    m_FWMetrics({"THICKNESS",
                 "LOWEST_RETURN",
                 "LAST_PATCH",
                 "FIRST_PATCH"
   // 2nd ADDITION: add the a name for your new metric
   // please do not forget the comma at the beginning
                 , "METRIC"
   // end of 2nd ADDITION
                })
{
   // The types should aggree with the fw metrics list
   m_types =
   {
      {"THICKNESS",3},
      {"FIRST_PATCH",9},
      {"LAST_PATCH",11}
    // 3rd ADDITION: give a number to your metric. This number must be unique
    // please do not forget the comma at the beginning
      , {"METRIC", 12}
   // end of 3rd ADDITION
   };
}

//-----------------------------------------------------------------------------
const std::vector<std::string> MapsManager::getNamesOfFWMetrics()const
{
   return m_FWMetrics;
}

//-----------------------------------------------------------------------------
void MapsManager::createMap(
        mapInfo *m_infoOfMap
        )
{
   if (m_map!=0)
   {
      delete m_map;
      m_map=0;
   }

   std::string s(m_infoOfMap->type);
   std::transform(s.begin(), s.end(), s.begin(), toupper);
   switch (m_types[s])
   {
  
   case 3:
      std::cout << "Density map\n";
      m_map = new DensityMap(m_infoOfMap->name,m_infoOfMap->obj);
      break;
  
  
   case 9:
      std::cout << "Length of first continues patch of non empty voxels\n";
      m_map = new FirstPatch(m_infoOfMap->name,m_infoOfMap->obj);
      break;
 
   case 11:
      std::cout << "Length of last continues patch of non empty voxels\n";
      m_map = new LastPatch(m_infoOfMap->name,m_infoOfMap->obj);
      break;

   // 4th ADDITION: link your metric class with the rest of the program
   case 12:
       std::cout << "Brief Discreption of the new Metric\n";
       m_map = new Metric(m_infoOfMap->name,m_infoOfMap->obj);
       break;
   // end of 4th ADDITION
   default:
      std::cout << std::string (s) << " is not a valid type of map";
      break;
   }
   // create and save the map
   if(m_map!=0)
   {
      m_map->createAndSave(m_infoOfMap->thres,m_infoOfMap->samp);
      delete m_map;
      m_map=0;
   }
}


//-----------------------------------------------------------------------------
MapsManager::~MapsManager()
{
   if(m_map!=0)
   {
      delete m_map;
   }
}

Once those steps are done you should be able to call your new metrics from the main program:

 ./DASOS -las myLasFile.LAS -map METRIC metric.asc

If you add the above class you should get a black asc file for all LAS files because all the values of the map are set to zero.

I hope you find that useful and if you have any questions please contact us at the following Google group:
https://groups.google.com/forum/#!forum/dasos---the-native-full-waveform-fw-lidar-software

Work Cited:

[1] 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

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

Bingo Game in C++

Once I was responsible of organising a bingo night and I wrote the following C++ code to make the game possible and allow myself to participate as well!

Since we were all programmers there, we enjoyed the idea of having our C++ code for playing the Bingo. By the end, I got slightly accused of customising the game in the favour of the event's organisers. Well, random is random and it wasn't my fault that we won the first two prizes!

The idea of the code is simple:
1. Firstly, an array of 100 elements is created and initialised to contained the values from 0 to 99 inclusive
2. Secondly, the array is shuttle using a random seed.
3. Then, the game starts! Every time the user presses 'enter' the next number inside the array is printed on the screen
4. and if there is a winner, the user may press 'b' to break the loop and end the game

The code of the Bingo is here:

#include < vector >
#include < iostream >
#include < time.h >
#include < stdlib.h >

int main(void)
{
    // create and initialise an array containing the values 0-99
    unsigned int legth = 100;
    std::vector < unsigned short int > m_bingoNumbers(legth);
    for(unsigned int i=0; i < m_bingoNumbers.size(); ++i)
    {
       m_bingoNumbers[i] = i;
    }

    // shuffle the array
    srand(time(NULL));
    for(unsigned int i=0; i< legth*2; ++i)
    {
       unsigned int rand1 = rand()%100;
       unsigned int rand2 = rand()%100;
       unsigned short int temp = m_bingoNumbers[rand1];
       m_bingoNumbers[rand1] = m_bingoNumbers[rand2];
       m_bingoNumbers[rand2] = temp;
    }

    // print numbers one by one 
    char b;
    for(unsigned int i=0; i< m_bingoNumbers.size(); ++i)
    {
       std::cin >> std::noskipws >> b;
       if(b=='b')
       {
           break;
       }
       if(b=='\n')
       {
          std::cout << " " << m_bingoNumbers[i];
       }
    }
    std::cout << "\n";
    return 0;
}

Here, there is also an example on how to compile and run the game:

$: g++ main.cpp -o bingo
$: ./bingo

 99
 4
 2
 50
 80
 41b

Finding all files of a given extension from a directory in C++

I have just found this useful recently and my blog it's a nice place to back it up for future reference. :)
The following small script takes as input an extension and a directory and prints all the files inside that directory that have that extension.

#include < iostream >
#include< stdio.h >
#include< cstdlib >
#include< iostream >
#include< string .h >
#include< fstream >
#include< sstream >
#include< dirent .h >
#include < vector >
// < extension > < directory >
int main(int argc, char *argv[])
{
    if(argc!=3)
    {
       std::cerr << "Too few arguments. Please include the following\n"
                 << "<.extension> \n";
       return EXIT_FAILURE;
    }
    std::string  dirStr(argv[2]);
    std::string extension(argv[1]);
    DIR *dir;
    struct dirent *ent;
    unsigned int count(0);
    if ((dir = opendir (dirStr.c_str())) != NULL) {
      /* print all the files and directories within directory */
      while ((ent = readdir (dir)) != NULL)
      {
        std::string current(ent->d_name);
        if(current.size()>extension.size())
        {
           std::string ext = current.substr(current.length()-extension.length());
           if(ext==extension)
           {
              std::cout << current << "\n";
              count++;
           }
        }
      }
    }
    std::cout << count << " files with extension " << extension << " found\n";
    std::cout << "   ***   EXIT   ***\n";
    return EXIT_SUCCESS;
} 

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.