option to export cloth from cloth simulator filter

[mikoontz]

Hello,

I was directed here from an issue I opened in the {lidR} package (r-lidar/lidR#244 (comment)). I'll copy that text here:

Here's my feature request: Would it be possible to extract the "cloth" that is used to classify the point cloud into ground/non-ground in the cloth simulator filter algorithm as a raster object?

My current workflow is to use the lidR::classify_ground() function with the algorithm = argument as csf() (with a few tuned parameter values) to classify ground/non-ground, then to use lidR::grid_terrain() with the algorithm = argument set to tin() to create a digital terrain model. During the interpolation process, there are sometimes some DTM artifacts associated with the uneven coverage of classified "ground" points within a scene, and the cloth itself might be a better representation of the terrain underneath the vegetation compared to the spatial interpolation between the "ground" points. I was finding the smoother cloth to be a better representation of the DTM when using the cloth simulator filter algorithm as implemented in CloudCompare (and then exporting the cloth itself as a raster layer).

Thanks for your thoughts and for considering it!

It looks like this is possible in the C++ code (https://github.com/Jean-Romain/RCSF/blob/7de8509dbe8ec413a77a8ead1e43cc5789d6e2de/src/CSF.cpp#L176-L177. Would it be straightforward to make this an option for {RCSF} to export the cloth as a {raster} object? (Then perhaps that option could be forwarded on to the use of the lidR::csf() algorithm.)

[Jean-Romain]

See jianboqi/CSF#27

It is probably doable but not trivial. I won't do it in a close future for sure.

[mikoontz]

Right on, and thanks for the link to the issue you opened in the original CSF code. If you'd be open to creative solutions or PRs, I might try to think about how this might be achievable with the CSF code as written. Maybe writing to a temporary file is okay? Or writing to disk in a more permanent way with it being left up to the user to re-import the cloth (in whatever format it's in) as a raster/data.table/point cloud?

[Jean-Romain]

It is ok but one requirement: RCSF is a standalone package with no dependency (expect Rcpp). So you can't add raster as a dependency in RCSF. If you want to get a RasterLayer it must be in lidR. Thus you must find a trick to achieve that. It won't be trivial. Maybe something like

In RCSF

CSF <- function(cloud, 
    sloop_smooth = FALSE, 
    class_threshold = 0.5, 
    cloth_resolution = 0.5, 
    rigidness = 1L, 
    iterations = 500L, 
    time_step = 0.65,
    save_cloth = "path/to/file")
{
}

In lidR

csf() is an algorithm. It works only within classify_ground() (the new name of lasground() in version 3.0.0). So there is no way to use it to return a RasterLayer. You could maybe create a new function csf_cloth() that run RCSF::CSF() and read back the file.

algo = csf()
algo(las)
#> Error : The 'csf' algorithm has not been called in the correct context. Maybe it has been called alone but it should be used within a lidR function. 

[mikoontz]

Hmm, I see what you mean. That makes sense that this would be a two step process: 1) add the ability to write the cloth to disk in the RCSF package, then 2) modify lidR so that calling the csf() function in the classify_ground() context triggers writing the cloth to disk.

What if you don't need to return a raster when calling csf() in the lidR package, but rather just trigger a side effect of RCSF writing to disk? It seems like that might be a simpler change?

So this line (copied below) in lidR::csf() might include an additional argument save_cloth which is the path to the file where the RCSF package will write the cloth to disk (as you've outlined above).

gnd <- RCSF:::R_CSF(cloud, sloop_smooth, class_threshold, cloth_resolution, rigidness, iterations, time_step, save_cloth)

If save_cloth is NA (the default), then the cloth is not written to disk and the lidR::csf() function performs as it currently does within the classify_ground() function (or lasground() for lidR <v3.0.0). If save_cloth is a filepath, then the cloth is written to disk when the RCSF:::R_CSF() function is called, but there is no change to what is currently returned when the lidR::csf() function is called within classify_ground().

If you think this seems reasonable, I can see about implementing it. Or if you think this way would be no good, then I'm curious to hear that too!

[Jean-Romain]

And how do you get back a RasterLayer. It is up to the user to deal with the file?

[mikoontz]

Yes, that'd be my first instinct-- to leave it up to the user to deal with the file as they'd like.

[Jean-Romain]

The problem with your option is that if you only want the cloth, the side effect is the classification of the points. I prefer a new function like csf_cloth() with an API like.

raster = csf_cloth(las, csf())

[mikoontz]

Yes, that's true. For my use case, I would definitely always want the point classifications too. But that of course doesn't mean that everyone would want that.

In your example above, do you imagine that the csf() function would be re-written as I suggested (with RCSF package modified first, then csf() modified to allow for the side effect of writing the cloth to disk) but that the new csf_cloth() function would be used in place of classify_ground()?

So if a user wants to classify ground/non-ground for their point cloud only, they'd use classify_ground(las, csf()).

If a user wants to both classify ground/non-ground for their point cloud AND export the cloth to disk, they'd use classify_ground(las, csf(save_cloth = "path/to/file"))

And if a user wants to just get the cloth without modifying the point cloud at all, they'd use csf_cloth(las, csf(save_cloth = "path/to/file")) which could write the cloth to disk, then read it in as a raster object (or some such)?

[Jean-Romain]

In your example above, do you imagine that the csf() function would be re-written as I suggested (with RCSF package modified first, then csf() modified to allow for the side effect of writing the cloth to disk) but that the new csf_cloth() function would be used in place of classify_ground()?

Yes

So if a user wants to classify ground/non-ground for their point cloud only, they'd use classify_ground(las, csf()). If a user wants to both classify ground/non-ground for their point cloud AND export the cloth to disk, they'd use classify_ground(las, csf(save_cloth = "path/to/file"))

Yes if the exported file is something that can be trivially be used. Not a raw text file. I think, no matter how, it should alway write a RasterLayer

And if a user wants to just get the cloth without modifying the point cloud at all, they'd use csf_cloth(las, csf(save_cloth = "path/to/file")) which could write the cloth to disk, then read it in as a raster object (or some such)?

In that case no need to provide a path. It must be saved in a temporary file and be loaded back as a RasterLayer. And here we need a trick to make everything consistent.

Honestly it is probably not the hardest part to export the cloth. Introducing this feature consistently is more difficult imo.

[mikoontz]

Copy that. I agree that a consistent implementation is important. If you definitely want to use a RasterLayer (or a stars object?), it looks like you'd have to be willing to add more dependencies to lidR, since I'm not seeing them in the DESCRIPTION file https://github.com/Jean-Romain/lidR/blob/master/DESCRIPTION. If you're up for that, I think this is totally possible to implement in a consistent way.

It looks like the nodes of the cloth are written to disk as a tab delimited text file in the Cloth.cpp code:
https://github.com/Jean-Romain/RCSF/blob/8ac13c36b841e0c69e7cbc52ebb31594e9bdd853/src/Cloth.cpp#L381-L383

That should be easy enough to read back into R, then convert the regular grid to a raster layer using raster::rasterFromXYZ(), then assign the proper coordinate reference system to be the same as the input las.

All of this hinges on first modifying RCSF to allow a save_cloth= argument, which will then write that text file to disk. I think it makes more sense to then let lidR handle the process of reading that file and converting to a RasterLayer, possibly then writing it back to disk as a geoTIFF. That will keep RCSF lightweight, and allow lidR to handle all the other important geospatial issues (like adding the right crs to the geoTIFF).

Modifying RCSF

I think that implmenting the change to RCSF wouldn't be too challenging (as you said). My guess is that a third argument needs to be passed to the do_filtering() method here:
https://github.com/Jean-Romain/RCSF/blob/0cddbbe4b6cf158ea1c5f31f6b55dd233e37f8ab/src/R_CSF.cpp#L57.

That third argument corresponds to the exportCloth boolean argument in the C++ code: https://github.com/Jean-Romain/RCSF/blob/7de8509dbe8ec413a77a8ead1e43cc5789d6e2de/src/CSF.cpp#L103-L105.

It looks like we'd also consider modifying the call to the .saveToFile() method in order to specify a known file location that we can use in the R code:
https://github.com/Jean-Romain/RCSF/blob/7de8509dbe8ec413a77a8ead1e43cc5789d6e2de/src/CSF.cpp#L177

The default (calling the .saveToFile() method with no arguments) will save the .txt file in the home directory as "cloth_nodes.txt", which maybe is fine (or better? it's certainly simpler):

https://github.com/Jean-Romain/RCSF/blob/7de8509dbe8ec413a77a8ead1e43cc5789d6e2de/src/Cloth.cpp#L367-L374

Modifying lidR

I could take a first pass at modifying lidR to help deal with the .tsv file written by RCSF.

Something like (borrowing from the current csf() function:

  f = function(las, filter)
  {
    . <- X <- Y <- Z <- NULL
    assert_is_valid_context(LIDRCONTEXTGND, "csf")
    las@data[["idx"]] <- 1:npoints(las)
    cloud <- las@data[filter, .(X,Y,Z, idx)]
    gnd <- RCSF:::R_CSF(cloud, sloop_smooth, class_threshold, cloth_resolution, rigidness, iterations, time_step, save_cloth)

begin psuedocode

If save_cloth is FALSE (the default):

• Proceed to RCSF:::R_CSF() as before, without any calls to save the file (R_CSF() will have to interpret "save_cloth=FALSE" appropriately).

If save_cloth is a file path "path/to/file":

• Let RCSF write the .txt file to disk in the default location: the home directory as a tab delimited text file named "cloth_nodes.txt".
• Use read.tsv() to read the file into R.
• Use raster::rasterFromXYZ() to make it into a RasterLayer.
• Assign the crs of the RasterLayer to be the crs of the las object.
• Use raster::writeRaster() to write the RasterLayer to the provided "path/to/file" as a geoTIFF
• unlink("cloth_nodes.txt")

If save_cloth is TRUE, as might be the case if someone uses the proposed csf_cloth() function and simply wants a RasterLayer returned without writing to disk:

• Let RCSF write the .txt file to disk in the default location: the home directory as a tab delimited text file named "cloth_nodes.txt".
• Use read.tsv() to read the file into R.
• Use raster::rasterFromXYZ() to make it into a RasterLayer.
• Assign the crs of the RasterLayer to be the crs of the las object.
• Use raster::writeRaster() to write the RasterLayer to a geoTIFF file to the home directory with a name "cloth_nodes.tif"
• unlink("cloth_nodes.txt")
• In the csf_cloth() function, use raster::raster() to read the "cloth_nodes.tif" file
• unlink("cloth_nodes.tif")

end pseudocode

    idx <- cloud$idx[gnd]
    return(idx)
  }

  class(f) <- LIDRALGORITHMGND
  return(f)
}

I also remembered that one of your implemented algorithms was able to be run on its own (i.e., not as the value to be passed to the algorithm= argument within the context of a different function). I looked around a bit and then found it-- it's dalponte2016(). Would you rather see the functionality of your proposed csf_cloth() be instead implemented as an option to call csf() outside of the call to classify_ground() for consistency with dalponte2016() (which can also be called within the context of segment_trees())? I am fine either way, but this seemed potentially inconsistent so I wanted to bring it up.

From ?dalponte2016:

Because this algorithm works on a CHM only there is no actual need for a point cloud. Sometimes the user does not even have the point cloud that generated the CHM. \code{lidR} is a point cloud-oriented library, which is why this algorithm must be used in \link{segment_trees} to merge the result with the point cloud. However the user can use this as a stand-alone function like this:

 chm = raster("file/to/a/chm/")
 ttops = find_trees(chm, lmf(3))
 crowns = dalponte2016(chm, ttops)()

[mikoontz]

The distinction between save_cloth=TRUE and save_cloth="path/to/file" doesn't feel quite polished enough (because there seems to be an unnecessary writing of a RasterLayer to disk, and then reading that same RasterLayer back into R, so we should treat this as a rough draft.

[Jean-Romain]

Two comments:

csf()() similary to dalponte()() is the most elegant solution and was my first idea but it cannot work. Not going through details but if you read the code and you understand well function factories you will find why.

Not stars dependency. I don't know this package. It looks interesting but I don't have time to learn it. raster does the job and is already a dependency of lidR.

If save_cloth is TRUE, as might be the case if someone uses the proposed csf_cloth() function and simply wants a RasterLayer returned without writing to disk:

You will face a technical issue that you will have difficulty to solve. But it can be solved. When you will be there we will discuss about that.