added -thin_with_voxel filter to LASlib

[tiagodc]

Hi Jean, I added a custom filter to LASlib. It performs a systematic grid filter based on a fixed voxel size at read time. It works the same way as -thin_with_grid, except that it keeps one random point per (cubic) voxel, instead of pixel, on the loaded cloud.

I integrated it into LASlib so it won't mess with any structures in the rlas package, being easy to use as a filter argument in read.las function and subsequently lidR functions, e.g. -thin_with_voxel 0.025.

I plan on releasing a refactored version of TreeLS later on this year. Dealing with TLS and MLS point clouds requires some robust 3D filters synced with data streaming, and instead of building functions inside a more specific package just for that, it would be nice to keep this functionality in a broader scope , i.e. rlas.

Cheers!

[Jean-Romain]

Interesting! And useful! This is a great PR :smile: However I cannot accept this PR as is.

1. You must explain me your code a bit. I will have to maintain that code so I must understand every single line. The code is not hard and I could understand it easily but I prefer that you also give me few key that will be helpful.
2. Don't modify deprecated.r in this PR. If I made a mistake open an issue or another PR and please explain what I did wrong.
3. This is the most important point. Do not modify the code of LASlib manually or at least the less as possible. I do agree that you must modify the parser but I'm sure you can put your class LAScriterionThinWithVoxel in a separate file such as rlas_extended_lasfilter.cpp with the proper includes. Let me explain why.

To compile LASlib in R I had to modify several hundred of lines of LASlib. 99% of the modifications are automatic. This is why there is this script laslib2R.sh. 1% is done manually and documented in the same script. Sometime I had to update laslib. So I copy all the file from LASlib I run the script and I make the manual modifications. This is prone to error and a very boring and long job actually. Moreover the script is out of date. I should really spend time to review all the comments. This is why you should not modify LASlib. And if you have to do it anyway, you must document you modifications in this file to enable me to redo the modification in case of LAslib update.

[tiagodc]

Alright! Thanks for all the info, I'll simplify it since some of the code is not used, describe it and minimize the changes on lasfilter.cpp.

Cheers

[tiagodc]

by the way... I haven't changed the LAScriterionThinPulsesWithTime class, I guess you meant the LAScriterion class, right?

[Jean-Romain]

My bad, this is just the github diff that hide large portion of code. I meant the argument parser line 3090. And indeed it makes more sense :wink:

[tiagodc]

ok, I have decoupled most of the code from LASlib and declared LAScriterionThinWithVoxel within lasfilter_voxelgrid.hpp. I reduced the voxelGrid class' code to the essentials. I also documented the changes made to lasfilter.cpp in the laslib2R.sh file. Bellow follows the description of the filter algorithm:

• Usage: read.las(cloud.las, filter="-thin_with_voxel 0.05")
• The voxelGrid object is a registry that stores unique voxel indices, according to a regular cubic spacing (from the example, one voxel is 0.05 x 0.05 x 0.05 units³)
• key voxelGrid members:
  • voxelSideLength: voxel length on each dimension (in cloud units)
  • dynamic_registry: set of unique vector<int> values, which represent the x, y and z voxel index
  • getLength(): calculates the voxel index (int) from an one dimensional coordinate (double). The index is defined as the distance in discrete steps (number of voxelSideLength) from the axis origin (x=0, y=0 or z=0)
  • checkRegistry(double x, double y, double z): takes one 3D point, calculates its voxel indices and check if the voxel is empty. If the voxel is already filled, returns false, otherwise returns true, allowing the point to be streamed further. This way, every voxel allows up to 1 point to be stored, and only the first streamed point to fit in a voxel is kept
  • resetDynamicReg(): clears the dynamic_registry
• Logical process: once called, -thin_with_voxel will instantiate a LAScriterionThinWithVoxel as an empty registry. Every point streamed will be checked for indexation. All first points assigned to a specific voxel will be indexed and subsequently read.

Hope it's all set!

[tiagodc]

By the way, the only change made to lasfilter.cpp was to include a parsing criterion inside LASfilter::parse, as you predicted ;)

[Jean-Romain]

Sounds good to me. I will read the code more carefully, try to find were it could crash and merge otherwise in few days.

You did not only add a new feature here. You also helped me to understand better how the -thin_with_grid works. It works as expected actually but I have never read the code before. Based on your code I'm pretty sure I can easily create a new LAScriterion to filter surface points such as thin_with_surface_point 1. Streaming equivalent to lasfiltersurfacepoint in lidR.

[tiagodc]

Great, I'm glad I could help :)

[Jean-Romain]

You must also modify LASfilter::usage() L1607 of lasfilter.cpp

[tiagodc]

OK, modified

[Jean-Romain]

I fixed the bug of scan angle test differently to ensure compatibility with future version. Support of LAS1.4 put me in the "dependency hell"

[Jean-Romain]

I looked very very quicky into laslib code of thin_with_grid but I looks extremely optimized. I think that it stores a single bit per cell to record the state of the cell. You used 3 x 8 bytes i.e. 192 bit per voxel in a container that relie on a complex algorithm to find the key. It means that you need roughly a full copy of the loaded point cloud coordinates to run your streaming filter. It may work for your own usage but I think is not suitable for more generic usage.

[tiagodc]

Well... I fully agree that the algorithm as is is not memory efficient, I programmed it in a way simple enough for it work fast.

... It means that you need roughly a full copy of the loaded point cloud coordinates to run your streaming filter. ...

About that part I don't fully agree, actually what limits the registry's size is the point cloud's bounding box and the voxel size. For very heavy TLS point clouds we can have up to hundreds of points per cm³, but for every cluster of points there will be only 1 vector<int> in the registry taking memory.

For heavy ALS clouds it shouldn't be a big issue regarding memory also, since those clouds are very sparse (when compared to TLS), thus most theoretical voxels will be empty space and never allocated into the index registry, even when using small voxel sizes.

I benchmarked the filter with different clouds here, and at worst case scenario it consumed roughly 1.2Gb of memory when reading a dense TLS point cloud (~ 165 million points), with extents of 250 x 250 x 200 m on the XYZ directions, respectively, using a voxel size of 2 cm. Read time was less than two minutes for the filtered data set on my laptop.

Summarizing, I agree that memory efficiency of the algorithm can be greatly improved, but as is I would say it is already suitable for more generic usage.

I'll dig deeper into the -thin_with_grid filter's code anyways, but can't guarantee I'll manage to follow Martin's steps that one, as it's quite complex (gotta study harder on my spare time :bowtie:)

[tiagodc]

Well... let me know your thoughts on those comments, I'll do my best to adapt it.

[Jean-Romain]

About that part I don't fully agree, actually what limits the registry's size it is the point cloud's bounding box and the voxel size. For very heavy TLS point clouds we can have up to hundreds of points per cm³, but for every cluster of points there will be only 1 vector in the registry taking memory.

Yes that was I said you have a copy of the loaded point cloud. Not the whole point cloud indeed. But it is anyway an issue. Let assume we want to stream 1 billions points spread on 400 km² to generate a new decimated point cloud we will run out of memory without any doubt.

Summarizing, I agree that memory efficiency of the algorithm can be greatly improved, but as is I would say it is already suitable for more generic usage.

No I don't think it is. As you mentioned it is safe for TLS centered on 0,0 but not for ALS because of integer overflow (and georeferenced TLS as well). And as I mentioned we will run out of memory in some (rare) cases. I really enjoy what you did but we must discuss about an improved implementation. I think the first step consists in understanding how ThinWithGrid deal without having a prior in the bounding box. Then we will be able to adapt.

To me it is a fun and interesting challenge :smile: Let me have a guess on how it works in ThinWithGrid. I think that the object is regularly reset to 0. Thus for a given file, you don't need to know the bbox. You take the first point as offset and you deal with that. A file being max 1 km wide the maximum range of the offset'ed values is only [-1000, +1000] whatever the first point used. When the streamer automatically change the read file the filter is reset to 0. Thus even if you are dealing with a file located 400 km away from the first one your are dealing with numbers that range [-1000, +1000] and a limited number of voxel (you don't store each 400 km² existing pixels/voxels). Just a guess.

[tiagodc]

Alright... I'll dig on that on my development branch as well, and I guess you're further ahead than me on interpreting ThinWithGrid. I figured he used the first point as a reference for the next ones, but still I don't follow how he stores the information on the pixels filled...

need some extra coffee to keep going on that train of thought... haha

[Jean-Romain]

No I'm not ahead in my understanding. It is just a guess.

The only point I'm sure is that a single bit is used to store the state of a cell. Line 1433 there is a bit shift operator U32 pos_x_bit = 1 << (pos_x%32); followed L1434 by a bitwise operator if ((*array)[pos_y][pos_x_pos] & pos_x_bit)

Bit-to-bit computation in brainfucking but it is possible to abtract that with a std::vector<bool>

[Jean-Romain]

Let me make another guess. In LASlib there is a U32** variable with is the matrix that stores the "raster" registry. This matrix can be extended as many time as you need. This is the role of malloc. Each cell of the matrix is a U32 and thus can store the state of 32 cells.

On 400 km² dataset processing, let say 200 x 200 km if the memory is never freed, then, with a grid resolution of 0.1, you have a final grid of (20000/0.1)² = 40.000.000.000 cells. We can use 1 bit per cell we need at least 40.000.000.000 bits = 5 Gb of memory to process a 400 km² ALS dataset. Thus the memory is necessarily freed regularly.

[Jean-Romain]

Ok I feel I understood a bit the code. It is very complex and clever actually.

[Jean-Romain]

Integer overflow

For integer overflow the following should make the job (I renamed voxelSideLength into voxel_spacing)

bool voxelGrid::checkRegistry(double x, double y, double z)
{
    if (voxel_spacing < 0)
    {
      x_anker = I32_FLOOR(x-x_anker / voxel_spacing);
      y_anker = I32_FLOOR(y-y_anker / voxel_spacing);
      z_anker = I32_FLOOR(z-z_anker / voxel_spacing);
      voxel_spacing = -voxel_spacing;
    }

    int nx = getLength(x_anker, x);
    int ny = getLength(y_anker, y);
    int nz = getLength(z_anker, z);

LAScriterionThinWithVoxel(F32 voxel_resolution)
{
      // member resolution is useless
      box.setVoxel(-resolution);
 };

Running out of RAM

That should not happen in lidR. lidR read one chunks at a time. It does not actually stream the whole dataset. At least few tiles could be streamed at once but that do not implies 4 billion of voxels.

Memory optimization

Edited because it was stupid

[tiagodc]

Awesome!! :D

I was actually studying bit-wise operations before and got to a conclusion very similar when reading thin_with_grid again!

I think I'll take on the challenge and try to use simple unsigned long int as containers and manage the voxels by bit-wise operations instead... one cup of coffee more and (hopefully) I'll push a new version tonight.

Thanks for framing all the logic above as well :)

[Jean-Romain]

Forget what I said about memory optimization. It was stupid and it will never work. I speak too fast especially at the end of the workday.

Moreover the main difference between raster based thinning and voxel based thinning is that in the first case 99% of the cells are likely to get a point inside while in the second case only a small percentage of the voxel actually exist. Thus using a std::set was pertinent. It could maybe be optimized but not as much as I was expected. I'll think about that tomorrow but give up with the vector of bool. It cannot work without a prior knowledge on the bbox.

But it does not mean that you cannot optimized a little. Check std::unordered_set if it is more appropriated in our case.

[tiagodc]

Well... I couldn't get quite there, I pushed another approach on this algorithm to my devel branch, if you wanna take a look...

I tried to manage it using nested vector registries + bit-wise operations at the innermost vectors to save some memory. This approach actually improved the speed of the algorithm on small point clouds, but memory management was less stable than the current version using std::set, becoming a problem on large point clouds...

I think I made some progress there, if I dig further into it I might be able to optimize further the voxel containers... but the code is getting more and more complicated this way.

Thanks for the tip on std::unordered_set, I'll check it tomorrow.

[Jean-Romain]

Let try with an unordered_set to figure out if it performs better and also you don't need to create a class voxelGrid only a LAScriterionThinWithVoxel. This way it will look more like the LASlib's criterion. I prefer to be consistent with LASlib style. Then we can consider more complex optimization as an interesting challenge. I read your code I think you understood the logic behind LAScriterionThinWithGrid like me and it is a good basis to go try further.

[tiagodc]

OK, in order to use the std::unordered_set I need to define a custom hash to build the container's keys on inserts.

Any idea on how to define this hash in an efficient way? I could only think of a way using strings, but I'm afraid the algorithm will become less memory efficient with string keys...

[Jean-Romain]

I didn't know that. Ok give up. This is going to be too complex. In few weeks we will step back and have more insight about what we did.

[Jean-Romain]

Alright! All that git stuff drive me crazy. I close this PR! Please make another clean PR from the current branch master. I reshaped your code and changed a bit my mind. Your PR should contain a single commit and you need to modify a single file (lasfilter.cpp). Rewritten like the following the code is very simple and it is much easier to put it into lasfilter.cpp. In you previous PR you added several files so I preferred to use your files to decouple your code and LASlib but reshaped like that its better. Don't worry with las2R.sh I must update it anyway. Change DESCRIPTION to add your name in contributors.

Thanks

class LAScriterionThinWithVoxel : public LAScriterion
{
public:
  inline const CHAR* name() const { return "thin_with_voxel"; };
  inline I32 get_command(CHAR* string) const { return sprintf(string, "-%s ", name()); };
  inline U32 get_decompress_selective() const { return LASZIP_DECOMPRESS_SELECTIVE_CHANNEL_RETURNS_XY | LASZIP_DECOMPRESS_SELECTIVE_Z; };
  inline BOOL filter(const LASpoint* point)
  {
    if(voxel_spacing < 0)
    {
      xoffset = point->get_x();
      yoffset = point->get_y();
      zoffset = point->get_z();
      voxel_spacing = -voxel_spacing;
    }

    I32 nx = I32_FLOOR((point->get_x() - xoffset) / voxel_spacing);
    I32 ny = I32_FLOOR((point->get_y() - yoffset) / voxel_spacing);
    I32 nz = I32_FLOOR((point->get_z() - zoffset) / voxel_spacing);
    std::vector<I32> nVals = {nx, ny, nz};

    return !dynamic_registry.insert(nVals).second;
  };
  void reset()
  {
    voxel_spacing = -voxel_spacing;
    xoffset = 0;
    yoffset = 0;
    zoffset = 0;
    dynamic_registry.clear();
  };
  LAScriterionThinWithVoxel(F32 voxel_spacing)
  {
    this->voxel_spacing = -voxel_spacing;
    xoffset = 0;
    yoffset = 0;
    zoffset = 0;
  };
  ~LAScriterionThinWithVoxel(){ reset(); };

private:
  double voxel_spacing;
  double xoffset;
  double yoffset;
  double zoffset;
  std::set< std::vector<I32> > dynamic_registry;
};

 REprintf("  -thin_with_voxel 0.1\n");

      else if (strcmp(argv[i],"-thin_with_voxel") == 0)
      {
        if ((i+1) >= argc)
        {
          REprintf("ERROR: '%s' needs 1 argument: voxel_spacing\n", argv[i]);
          return FALSE;
        }
        add_criterion(new LAScriterionThinWithVoxel((F32)atof(argv[i+1])));
        *argv[i]='\0'; *argv[i+1]='\0'; i+=1;
      }