Hi everyone. I’ve been working on an experimental plugin for using Julia as an embedded stage in PDAL. There is a similar feature for using Python already.
I’d like to get a feel for how much people are using Julia for Point Cloud tasks currently and how interested they would be in using it if it were easily accessible in PDAL?
I am no expert on Julia (this project is motivated by our data-science team favouring it for other tasks), so from the guidance of my team members, the most useful packages to base the interface around would be
LasIO.jl
RoamesGeometry.jl
TypedTables.jl
Any thoughts on how to structure this interface would be appreciated. The Python version expects a function with a dictionary of Numpy arrays (one per dimension in the PC) as its only argument and returns the same thing. The metadata/header is also exposed as a global.
The code is here, contributions very welcome! It’s still in a very rough state and is unlikely to work outside of Debian
Hi and welcome! Being able to use Julia in PDAL pipelines would be great. I thought about doing that as well a while ago, but never got around to it.
My understanding is that this is most useful for doing custom steps that are not already available in PDAL. So input and output are probably already covered, and LasIO may not be very helpful, since PDAL’s internal data model is not equal to LAS I believe.
While Python uses a dictionary of Numpy arrays, using a Dict of Arrays in Julia is probably not optimal. I have not looked into PDAL’s data model, but perhaps you can use something similar in Julia? It would be great for performance if PDAL does not have to make any copies.
It is probably also a good idea to ask about the design to the PDAL mailing list, and let them know that you are working on this.
I used to work with a lot of point clouds at FugroRoames with @greghislop and @andyferris and others. We definitely did a lot of point cloud processing there, and a fair bit of it in Julia. TypedTables.jl had its origins in work we did there as well as RoamesGeometry.jl (and StaticArrays.jl, Rotations.jl, CoordinateTransformations.jl, Proj4.jl, PlyIO.jl, Displaz.jl a good bit of Geodesy.jl and various other things!) I almost got involved in PDAL development at the time but it never really caught on.
Certainly PDAL integration could be handy for an organization which already uses PDAL heavily. As I understand it, the PDAL data structure is essentially tabular so it should definitely be reflected into Julia as a Table of some kind. A TypedTables.Table is the obvious choice for this kind of narrow-but-long table without missing data. However there’s many possible Julia table types which share a common interface Tables.jl which allows them to interoperate quite well and share APIs. Depending on your needs you could try other tables, but TypedTables is a good starting point.
By the way, for visualization you may find my old project Displaz.jl interesting. It’s admittedly a bit crusty at this point, but it’s built specifically for lidar point clouds (las/laz files and others) and we used it very heavily and successfully.
I had a quick browse through your embedding code. One thing you should beware of is that you need to carefully GC root a lot of variables. For example the arrays made at PDAL-julia/Invocation.cpp at 9399cc8a63e96ef2d9fbc8fd71077b05bddb18fe · cognitive-earth/PDAL-julia · GitHub look like they may not be protected and could be garbage collected at any moment. Generally you should consider that any jl_ function could trigger a GC and delete any variables you haven’t protected in some way. Protecting variables can be done by
Using JL_GC_PUSH for variables within a function
Adding your variables to some Julia module in which case they’ll be rooted by that module. This can be a bit global and non-threadsafe.
More sophisticated, by using the GC extensions in julia_gcext.h, in particular jl_gc_set_cb_root_scanner can allow the GC to discover roots held in your C++ objects.
Thanks for the feedback @visr, glad you think its a useful idea. My intention with LasIO was to use the LasPoint abstraction, but after investigating further it seems that a TypedTable representation is more appropriate anyway.
I’ve started with a TypedTables interface and will try others if the need arises. So far I have been able to pass the PDAL data structure into a user-supplied function in Julia as a TypedTable which is promising. It may be too early to worry about performance, but how would I go about telling if I am doing any unnecessary copies in this function?
We are indeed already using Displaz when developing, although we are targeting Cesium more generally for visualisation of our datasets.
I am not quite at the point of worrying about freeing memory and GC but thats definitely on the TODO list. I’ve added a comment to remind me to take a look at julia_gcext.h when I get to that.
You can try @allocated (or the allocation tracking in BenchmarkTools which seems somewhat inexplicably more reliable in julia 1.4). Generally that aspect of the code looks reasonable; unless you call copy on your arrays which were passed from the C code, the storage should be shared. I wouldn’t worry about small amounts of allocation to set up the Table itself.
Some other things I noticed:
You can create the Table very slightly more directly (without the Dict temporary) using something like (assuming the C++ code creates names and vals):
julia> names = [:column1, :column2]
2-element Array{Symbol,1}:
:column1
:column2
julia> vals = [[1,2], ["x","y"]]
2-element Array{Array{T,1} where T,1}:
[1, 2]
["x", "y"]
julia> Table(NamedTuple{tuple(names...)}(tuple(vals...)))
Table with 2 columns and 2 rows:
column1 column2
┌─────────────────
1 │ 1 x
2 │ 2 y
Your Symbol conversion looks a bit complex. I’d tentatively suggest making the desired Vector{Symbol} of column names on the C++ side directly. See the C function jl_symbol to make a symbol from a C string (return value from jl_symbol may be cast to a jl_value_t*).
The GC stuff is not about freeing memory, but about protecting it from being freed behind your back.Thus you usually should worry about it immediately during development otherwise you can get weird segfaults or other memory corruption which will really ruin your day!
If you’re used to manual memory management in C or using smart pointers in C++, the required mindset can be rather confusing: in some sense it’s the opposite of these systems: rather than having unexpected leaks, you’ll get things being unexpectedly freed. I assume you’ve seen it already, but do make sure to read Embedding Julia · The Julia Language in detail.
Fun to hear you’re using Displaz for development. It’s definitely a dev tool, not a polished instrument for end users! I don’t have time to develop it a lot these days, but if you or your team are interested in improving it I’m always happy to provide some advice and code review.
This is absolutely correct in the context of PDAL. Generic and efficient point cloud processing is challenging because you’d like to add arbitrary attributes to points as part of the processing pipeline, or remove them to free up memory. This begs for a tabular data structure of some kind rather than the simpler array-of-struct approach which LasIO gets away with.
Unfortunately C++ makes it really hard to make a good tabular data structure which is memory and time efficient, has good syntax, and where the schema isn’t known at compile time. We tried more than once (including an early effort by myself) at Roames and overall we failed to make something that most people wanted to use. To some extent, TypedTables was a reaction to those failures. I still think it’s the correct approach and hits a sweet spot of really convenient syntax and efficiency.
Another interesting thought: given you’re worried about copies, I though I’d point out that you can have your point cloud in “multiple coordinate systems at once” by combining CoordinateTransformations with the lazy MappedArrays package as I just realized over here: API for returning scaled and offset point cloud · Issue #34 · visr/LasIO.jl · GitHub
This works even with fixed-point integer coordinates like are used in LAS. So in principle you could have your LAS files mmaped from disk into memory, represent them as a normal Julia MappedArray of points in a different geographic coordinate system (eg, via Geodesy.jl / CoordinateTransformations.jl). Then do random access on your point cloud, only paying for which part of the point cloud you read or write to.
All this at a very high level of abstraction and convenience.
I’m a little stuck on this. I have quite complex types that I’m passing into Julia (a few arrays and a function), and can’t figure out how to use any of those approaches. I tried the 3rd option but it seems that the julia_gcext.h header is only available in newer versions of Julia (I have 1.0.3 as that is the latest version in Debian stable). When I upgraded to the latest one (1.4.1), it compiled but then segfaults when running.
I also am yet to see any issues around data being unexpectedly freed, is there some way I can trigger this condition reliably?
This could be due to the very issue we’re discussing!
On older Julia, gcext didn’t exist as you can see — in that case, the only real way to proceed is to globally root your temporary objects within some object which is a module-global.
Having said that, I looked a bit closer and I think maybe no julia objects outlive the Invocation::execute() invocation? In that case it’s all a lot easier: just use JL_GC_PUSH to root a few things temporarily and you’ll be fine.
You can always run GC.gc() to trigger GC at an inconvenient point.
In general, what I’d suggest to make rooting easier is to avoid using std::vector to hold Julia objects, but rather use a julia array (jl_alloc_vec_any is simplest and likely efficient enough if it’s holding only a handful of other arrays containing the real data). Then you only need to root the single parent array and the runtime will know about all child objects.
Also, I’m not sure why you’ve got pointer-to-pointer jl_value_t** in the following code. That seems weird and perhaps only “works” because sizeof(jl_value_t*) == sizeof(jl_value_t**)!
Okay thanks, I’ve made this change and everything seems to still work.
I’ve tried using the JL_GC_PUSHARGS macro, and it is segfaulting… See it commented out here. Any ideas on this one?
signal (11): Segmentation fault
in expression starting at no file:0
unknown function (ip: 0x7f611d8c64e5)
jl_apply_generic at /lib/x86_64-linux-gnu/libjulia.so.1 (unknown line)
jl_call at /lib/x86_64-linux-gnu/libjulia.so.1 (unknown line)
_ZN4pdal5jlang10Invocation7executeERSt10shared_ptrINS_9PointViewEENS_12MetadataNodeE at /home/jtfell/PDAL-julia/pdal/./libpdal_plugin_filter_julia.so (unknown line)
_ZN4pdal11JuliaFilter3runESt10shared_ptrINS_9PointViewEE at /home/jtfell/PDAL-julia/pdal/./libpdal_plugin_filter_julia.so (unknown line)
_ZN4pdal11StageRunner3runEv at /usr/local/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal5Stage7executeERNS_14BasePointTableERSt3setISt10shared_ptrINS_9PointViewEENS_13PointViewLessESaIS6_EE at /usr/local/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal5Stage7executeERNS_14BasePointTableE at /usr/local/lib/libpdal_base.so.10 (unknown line)
_ZN42JuliaFilterTest_JuliaFilterTest_test1_Test8TestBodyEv at ./julia_filter_test (unknown line)
_ZN7testing8internal35HandleExceptionsInMethodIfSupportedINS_4TestEvEET0_PT_MS4_FS3_vEPKc at ./julia_filter_test (unknown line)
_ZN7testing4Test3RunEv.part.657 at ./julia_filter_test (unknown line)
_ZN7testing8TestInfo3RunEv.part.658 at ./julia_filter_test (unknown line)
_ZN7testing9TestSuite3RunEv.part.659 at ./julia_filter_test (unknown line)
_ZN7testing8internal12UnitTestImpl11RunAllTestsEv at ./julia_filter_test (unknown line)
_ZN7testing8UnitTest3RunEv at ./julia_filter_test (unknown line)
main at ./julia_filter_test (unknown line)
__libc_start_main at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
_start at ./julia_filter_test (unknown line)
Allocations: 16569645 (Pool: 16566726; Big: 2919); GC: 27
[1] 27023 segmentation fault
Possibly unrelated, but when I try to run the plugin as a stage in PDAL (as a shared library), I get this error for the same build that has the tests passing (the GC stuff above is commented out) and the julia function is identical to the one used in the tests.
signal (11): Segmentation fault
in expression starting at no file:0
unknown function (ip: 0x7fe67c3e0866)
unknown function (ip: 0x7fe67c2d400d)
jl_get_global at /lib/x86_64-linux-gnu/libjulia.so.1 (unknown line)
_ZN4pdal5jlang10Invocation7executeERSt10shared_ptrINS_9PointViewEENS_12MetadataNodeE at /home/jtfell/PDAL/build/./lib/libpdal_plugin_filter_julia.so (unknown line)
_ZN4pdal11JuliaFilter3runESt10shared_ptrINS_9PointViewEE at /home/jtfell/PDAL/build/./lib/libpdal_plugin_filter_julia.so (unknown line)
_ZN4pdal11StageRunner3runEv at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal5Stage7executeERNS_14BasePointTableERSt3setISt10shared_ptrINS_9PointViewEENS_13PointViewLessESaIS6_EE at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal5Stage7executeERNS_14BasePointTableE at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal15PipelineManager7executeENS_8ExecModeE at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal14PipelineKernel7executeEv at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal6Kernel8innerRunERNS_11ProgramArgsE at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal6Kernel11doExecutionERNS_11ProgramArgsE at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal6Kernel3runERKSt6vectorINSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEESaIS7_EERSt10shared_ptrINS_3LogEE at /home/jtfell/PDAL/build/lib/libpdal_base.so.10 (unknown line)
_ZN3App7executeERSt6vectorINSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEESaIS6_EERSt10shared_ptrIN4pdal3LogEE at ./bin/pdal (unknown line)
main at ./bin/pdal (unknown line)
__libc_start_main at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
_start at ./bin/pdal (unknown line)
Allocations: 3248 (Pool: 3240; Big: 8); GC: 0
[1] 24283 segmentation fault
Yeah, this is not the correct usage. JL_GC_PUSHARGS will allocate space on the stack for a known number of arguments, so in your usage, args gets set to an uninitialized pointer to the new space (not julia_args.data()!
Correct use would look more like:
size_t nargs = calculate_nargs();
jl_value_t **args;
JL_GC_PUSHARGS(args, nargs);
for (size_t i = 0; i < nargs; ++i)
args[i] = some_arg(i);
But regardless of that there are still several other reasons why your GC rooting is broken. In particular, the following line will allocate memory in the loop, but you haven’t protected any of the elements of arrayBuffers from being collected:
As I said above, I think you’d find it easier to avoid using std::vector to hold jl_value_t, it’s just going to make correct rooting difficult.
Instead, just use a Julia array. Here’s the pattern you need for prepareData:
jl_value_t *prepareData(PointViewPtr& view)
{
jl_value_t *arrayBuffers = jl_alloc_vec_any(0);
// As soon as you have this `jl_value_t`, you need to protect it
// immediately, before any other `jl_` functions are called.
// Any call to a `jl_*` function may cause arrayBuffers to be freed
// behind your back.
//
// Luckily, anything stored in arrayBuffers will be known to the GC
// via the arrayBuffers root so you don't need to root it separately.
JL_GC_PUSH1(&arrayBuffers);
// ...
for (auto di = dims.begin(); di != dims.end(); ++di)
{
// ...
jl_array_t* array_ptr = jl_ptr_to_array_1d(array_type, data, view->size(), 0);
// Put array_ptr into arrayBuffers immediately to protect from GC:
jl_array_ptr_1d_push(arrayBuffers, array_ptr);
}
// ...
// Critically important: you must pair a POP with every PUSH
JL_GC_POP();
// Also critically important: now that you've POP'd,
// arrayBuffers is no longer rooted, so you should
// avoid calling any other `jl_` function until it is
// rooted again by Invocation::execute()
return arrayBuffers;
}
Thanks so much for that code sample - everything has just come together! I successfully ran a Julia script embedded in a PDAL pipeline (build steps are in the project README if you’re interested).
I am currently relying on a globally installed version of TypedTables.jl locally, how would you suggest I package that (along with other helpful libs eg. RoamesGeometry) up with my shared library? I tried to use the Pkg module a few different ways (calling Pkg methods at the top of my Julia file, using a project.toml file) but couldn’t figure out a simple way to get it running in a venv-style setup.
The way to precisely define dependencies is definitely to have a Project.toml and a Manifest.toml which you distribute and use with Pkg to create the appropriate Julia environment. Precisely how this is done depends on the way you want to set things up.
For a normal Julia library you’d just create a package and add any dependencies to the Project.toml along with acceptable version number ranges for compatibility. PkgTemplates.jl is a great way to create a Julia package with the standard file structure, Project.toml, test harness integration and CI configuration files, templates for documentation, etc etc. So in principle you could just create a normal julia package PDALPipelines.jl or some such, and load this into your embedded Julia.
The problem is that this doesn’t really define a self-contained application and pin down the precise versions of libraries you’d like to deploy with. For that, you want a Manifest.toml which you can use with Pkg.instantiate to create a deployable environment with precise versions of each library installed.
If you had a normal julia application with a “main.jl” file and with project and manifest, you’d set up the environment with Pkg.instantiate and then run Julia in “project mode” with julia --project=path/to/project main.jl where path/to/project contains the project and manifest.
However your case is a bit unusual since you’re embedding. Ideally you should probably compile a custom sysimage from your Manifest.toml for super fast startup. GitHub - JuliaLang/PackageCompiler.jl: Compile your Julia Package is the tool for that. Then your C code can load the custom sysimage as required using the appropriate path with jl_init_with_image.
All this will get much more complicated if you want to allow your users to load arbitrary Julia packages to use from their filtering code. But I think for a first version you can probably just install a few particularly useful libraries and see how that goes.
I have the package compiler working on my local dev machine (debian) and the startup is nice and quick which is great. I think we’ll go with manually loading a useful set of libraries and can rebuild the docker image with extras if and when we need.
Sorry to keep bringing you errors, but I am trying to build an alpine-linux based docker image and am getting this weird error at runtime:
fatal: error thrown and no exception handler available.
#<null>
jl_errorf at /julia-source/src/rtutils.c:77
jl_load_dynamic_library at /julia-source/src/dlload.c:233
jl_get_library_ at /julia-source/src/runtime_ccall.cpp:50
jl_get_library_ at /julia-source/src/runtime_ccall.cpp:39 [inlined]
jl_load_and_lookup at /julia-source/src/runtime_ccall.cpp:61
jlplt_pcre2_compile_8_373 at /usr/local/lib/pdal_jl_sys.so (unknown line)
compile at ./pcre.jl:120
compile at ./regex.jl:72
#occursin#363 at ./regex.jl:172 [inlined]
occursin at ./regex.jl:172 [inlined]
isdirpath at ./path.jl:117 [inlined]
normpath at ./path.jl:337
abspath at ./path.jl:383
abspath at ./path.jl:391
__init__ at ./sysinfo.jl:119
jfptr___init___18703 at /usr/local/lib/pdal_jl_sys.so (unknown line)
jl_apply at /julia-source/src/julia.h:1700 [inlined]
jl_module_run_initializer at /julia-source/src/toplevel.c:74
_julia_init at /julia-source/src/init.c:788
jl_init_with_image__threading at /julia-source/src/jlapi.c:74 [inlined]
jl_init_with_image__threading at /julia-source/src/jlapi.c:63
_ZN4pdal5jlang10Invocation10initialiseEv at /usr/local/lib/libpdal_plugin_filter_julia.so (unknown line)
_ZN4pdal5jlang10InvocationC2ERKNS0_6ScriptENS_12MetadataNodeERKNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEE at /usr/local/lib/libpdal_plugin_filter_julia.so (unknown line)
_ZN4pdal11JuliaFilter5readyERNS_14BasePointTableE at /usr/local/lib/libpdal_plugin_filter_julia.so (unknown line)
_ZN4pdal5Stage7executeERNS_14BasePointTableERSt3setISt10shared_ptrINS_9PointViewEENS_13PointViewLessESaIS6_EE at /usr/local/lib/libpdal_base.so.10 (unknown line)
_ZN4pdal5Stage7executeERNS_14BasePointTableE at /usr/local/lib/libpdal_base.so.10 (unknown line)
_ZN42JuliaFilterTest_JuliaFilterTest_test1_Test8TestBodyEv at ./julia_filter_test (unknown line)
_ZN7testing8internal35HandleExceptionsInMethodIfSupportedINS_4TestEvEET0_PT_MS4_FS3_vEPKc at ./julia_filter_test (unknown line)
_ZN7testing4Test3RunEv.part.0 at ./julia_filter_test (unknown line)
_ZN7testing8TestInfo3RunEv.part.0 at ./julia_filter_test (unknown line)
_ZN7testing9TestSuite3RunEv.part.0 at ./julia_filter_test (unknown line)
_ZN7testing8internal12UnitTestImpl11RunAllTestsEv at ./julia_filter_test (unknown line)
_ZN7testing8UnitTest3RunEv at ./julia_filter_test (unknown line)
main at ./julia_filter_test (unknown line)
unknown function (ip: 0x7ff7ab650230)
unknown function (ip: 0)
unknown function (ip: 0x745f7265746c6965)
It looks like you’re moving your sysimage around, which may have broken the ability to find many of the other .so files which need to go with the Juila distribution (including maybe libpcre.so?)
However I haven’t tried embedding myself so I don’t know precisely what what the problem is here.
You were right, it wasn’t looking in the right place for all the dependency libraries. I have managed to get it to work, but only after doing some questionable file copying and it only works when I’m in the folder where the plugin was built.
and run PDAL with the filter: pdal translate /root/PDAL-julia/pdal/test/data/1.2-with-color.las ./julia-out.las julia \ --filters.julia.script=/root/PDAL-julia/pdal/test/data/test1.jl \ --filters.julia.module="TestModule" \ --filters.julia.function="fff";
This definitely doesn’t seem like the right thing to be doing. I suppose symlinking the deps is better but I don’t understand why it has to be in the ~PDAL-julia/pdal folder to run successfully. Any pointers on this? Or perhaps someone you think would be able to help? Thanks.
I’ve got a working build on alpine-linux (published on Docker Hub), thanks again for the help getting to this point!
I’ve been experimenting with implementing some simple filters and I need to do some spatial lookups. It seems like GridIndex from RoamesGeometry along with AcceleratedArrays.jl is what I’m looking for, but the TypedTables.FlexTable datatype I’m using for IO between PDAL and Julia doesn’t work with these packages.
These are probably obvious, basic questions but I’m still figuring out the ecosystem
According to the TypedTables README, " You can speed up data analysis tasks with acceleration indices, by using the AcceleratedArrays.jl package" but I can’t find an example of doing this
More generally, RoamesGeometry takes the approach of creating Tables with “grouped” columns like “position”, “intensity”, “classification” while I have just attached every dimension (as they are called in PDAL) at the top level eg. “X”, “Y”, “Z”, “Classification”, “Intensity”. Is there some learned wisdom around this choice? I’m leaning towards thinking that this is to allow GridIndex to be written as it is? Its not too late for me to switch to this approach but I just want to understand the tradeoffs I’d be making. There’s definitely an advantage to maintaining the abstractions from PDAL from an ecosystem point of view
I wasn’t really involved in RoamesGeometry. It looks like it currently assumes AbstractVector{<:SVector{3,T}} for some T as the array of points. Arguably this is a bit too specific, but you can probably still get it to work with your use case by constructing a StructArray{SVector{3,Float64}} which should efficiently wrap your separate x,y,z component arrays to present it as a single array without copying.
Well I don’t know about learned wisdom, but the idea is that position is simply much more convenient to manipulate as a vector with three components than as individual components; it’s a higher level of abstraction. For example, suppose you want the distance between two nearby points at indices i and j in the point cloud. It’s just norm(pos[i] - pos[j]) (ok approximately; ignoring for now curved geospatial coordinate systems). No need to fiddle with implementing the Euclidean norm yourself.
It’s the same thing for attributes like color: really the color is one entity which you’d like to manipulate with the standard color space utilities provided Colors.jl. So the appropriate abstraction is to have the color as a single attribute rather than splitting it up.
I suspect PDAL does individual components because it’s just so painful and difficult to implement these kind of powerful abstractions in C++ that people just gave up!
For your plugin, you’ll have to decide whether it’s worth it to your users to have an extra layer of translation between the Julia data structures and the PDAL view of the world. It will definitely make things easier on the Julia side and allow many more standard Julia libraries to “just work” with the translated attributes. The downside is extra complexity in translating between PDAL and Julia.
