Visualize terrain models

[joa-quim]

The final idea would be to able to see terrain models like in

http://w3.ualg.pt/~jluis/mirone/images/marte_planar.jpg

but never mind about the details with the texture and the points. Other more basic things to be solver before it. I'm attaching a simple code with questions in comments. It needs GMT.jl installed (we know that it works fine both in Win and OSX, no linux testing so far)

visu_GMT.jl.txt

[SimonDanisch]

It looks like your x,y vectors are basically unit ranges, which means you can display them as a heightfield. Unlike e.g. Matlab, I don't do anything to focus on the object, or normalize Z automatically when it's too large. So you need to do this yourself currently. One way to do this is via the so called model matrix, which encodes arbitrary transformations. E.g.:

using GLAbstraction
visualize(any_vizz_type, model=translationmatrix(Vec3f0(0.5, 10, 0))*scalematrix(Vec3f0(0.01))*rotatiomatrix(axis, deg2rad(180f0))

With the heightfield visualization, you can also pass the xy extent:

using GMT, GLVisualize, GLAbstraction, Colors
fname = "bathy_1m.nc";
G = gmt("gmtread -Tg " * fname)
Z = G.z
maxi = maximum(Z)
mini = minimum(Z)
xmin,xmax = G.x[1], G.x[end]
ymin,ymax = G.y[1], G.y[end]
#center them
xwidth_half = (xmax-xmin)/2f0
xmin, xmax = -xwidth_half, xwidth_half

ywidth_half = (ymax-ymin)/2f0
ymin, ymax = -ywidth_half, ywidth_half
zwidth_half = maxi-mini
Z = (Z - mini) ./ (maxi-mini) #normalize
w,r = glscreen();
colormap = RGBA{U8}[
    RGBA{U8}(0.0,0.08203125,0.96875,1)
    RGBA{U8}(0.03125,0.9375,0.984375,1)
    RGBA{U8}(0.25,0.80078125,0.36328125,1)
    RGBA{U8}(0.87890625,0.91015625,0.1015625,1)
    RGBA{U8}(0.75,0.125,0.0,1)
 ]
robj = visualize(Z, :surface, grid_min=Vec2f0(xmin, ymin), grid_max=Vec2f0(xmax, ymax), color=colormap)
view(robj)
r()

This yields this result:

Better camera control is on its way and I plan to automatically adjust the camera on an object in GLPlot.

Accepting ranges and Float64 arrays is on its way in https://github.com/JuliaGL/GLVisualize.jl/pull/43 . Although, it will always end up as Float32. If the areas are really that large with that much detail in the visualization, I'm pretty sure that GLVisualize couldn't handle it right now. For that we would need tiles+culling. In it's simple form quite easy to implement.

[joa-quim]

Simon, that's fantastic thanks. A couple things more. Float64 for the z arrays is overkill Float32 should be more than enough. It's for the vector coordinates that one might need an higher resolution. Any chances to have an axis system plotted as well? I see also that the surface was shaded illuminated but can see any instruction to do so. What is the shading direction, can we change it?

If you have a Windows machine and want to try my other (big) project (that I would love to port to Julia) http://w3.ualg.pt/~jluis/mirone it has included the 3D Fledermaus free viewer that has many features that would be lovely to see in GLVisualize. I can give you a quick recipe on how see it working. The data format used by that viewer even uses Int16 from the z array.

[SimonDanisch]

Mirone

I'm on windows, yes ;)

Axis are slowly on their way. Right now you can already do:

robj = visualize(...)
view(robj)
view(visualize(boundingbox(robj).value, :grid))

If you like some other camera, the best idea is to just implement it yourselve, since I have other more pressing to do items on my list. It shouldn't be too hard with this script:

function camtranslation(arrow, speed)
    arrow == :up    && return Vec3f0(-speed,0,0)
    arrow == :down  && return Vec3f0(speed,0,0)
    arrow == :left  && return Vec3f0(0,speed,0)
    arrow == :right && return Vec3f0(0,-speed,0)
    Vec3f0(0)
end
function camrotation(buttons, speed)
    GLFW.KEY_A in buttons && return Vec3f0(0,0,-speed)
    GLFW.KEY_D in buttons && return Vec3f0(0,0,speed)
    GLFW.KEY_W in buttons && return Vec3f0(0,speed,0)
    GLFW.KEY_S in buttons && return Vec3f0(0,-speed,0)
    Vec3f0(0)
end

function create_cam(window, lookatvec=Vec3f0(0), eyeposition=Vec3f0(2))
    @materialize buttonspressed, arrow_navigation, window_size = window.inputs
    translation = map(camtranslation, arrow_navigation, Signal(0.1f0)) # map --> applies callback to signal
    theta = map(camrotation, buttonspressed, Signal(0.1f0))

    PerspectiveCamera(
        window_size,
        eyeposition,
        lookatvec,
        theta,
        translation,
        Signal(0f0),  # Zoom
        Signal(41f0), # Field of View
        Signal(1f0),  # Min distance (clip distance)
        Signal(100f0) # Max distance (clip distance)
    )
end
w.cameras[:my_cam] = create_cam(w)
view(robj, method=:my_cam)

I must admit, that I found a bug in GLWindow when writing this up and haven't found the source yet. Seems like I'm missing some events, so that the keys are not removed correctly, which means sometimes you move in the wrong direction. I'm on it!

Best, Simon

[joa-quim]

OK, to use the fledermaus viewer from within Mirone just open one grid by drag-n-drop on to the Mirone icon (almost everything in Mirone opens by drag-n-drop). Use for example the bathy_1m.nc. Next you can play with the illumination by clicking the icon with a hill and a Sun on top. To call the viewer hit the icon eye (next to the GoogleEarth icon).

Two nice things in it (besides the fact that it is extremely efficient even for big grids). By dragging the vertical arrow one can change the vertical exaggeration. Also, moving the mouse displays the x,y,z coordinates at lower left corner.

One further question. Do you think that your graphical infrastructure could be used to migrate Mirone to Julia? I don't mean right now but 'with time'. I've searched for a Matlab alternative but I didn't find anything convincing (given my GUI programming skills). GTK.jl could be one but simple can't convince myself that GTK is that a good bat.

Thanks for your patience and trials of my codes.

Joaquim

[SimonDanisch]

Miron kinda works, but crashes at random tasks.

What would be a really large grid? Would be fun to see how far we can push GLVisualize. Actually, as long as the grid fits into GPU memory, we could use some tricks to keep a smooth interaction. If it's bigger, we'd need to start streaming, which is a planned but bigger project.

You can actually build 2D/3D GUIs in GLVisualize. I hope they will gradually get more and more appealing. Have you seen: http://julialang.org/blog/2015/10/glvisualize/ That's basically the interaction that is possible right now. The cubecamera is closest to what you're looking for (the little cube steering the camera). As you can see, drag and drop is already present ;) If you properly integrate GMT into FileIO and GLVisualize, you could already drag and drop it easily (in GLPlot). You could use sliders to change all kind of values in a visualisation, e.g. the light properties. Light can be changed via this keyword argument:

light  = Vec3f0[Vec3f0(1.0,1.0,1.0), Vec3f0(0.1,0.1,0.1), Vec3f0(0.9,0.9,0.9), Vec3f0(20,20,20)]
# [diffuse_color, specular_color, ambient_color, position]
# I actually duplicate the light, using -position, to make the back not pitch black.
# For a visualization, signals can be used, too:
using Reactive, GLAbstraction
pos = Signal(Vec3f0(20))
light = map(pos) do p
    Vec3f0[Vec3f0(1.0,1.0,1.0), Vec3f0(0.1,0.1,0.1), Vec3f0(0.9,0.9,0.9), p]
end
view(visualize(..., light=light))
@async renderloop()
push!(pos, newpos)

I see if I can get you a slider example. A simple one is here: https://github.com/JuliaGL/GLVisualize.jl/blob/master/examples/image_processing.jl

Best, Simon

[joa-quim]

Hmm, it's not supposed that Mirone crashes (though I know it may happen). If it is easily reproducible please tell me or file an issue in https://github.com/joa-quim/mirone but before that please update Mirone by doing "Help -> Check for updates". You may have to repeat the process until "Help -> About" says 25 November.

I learned about http://julialang.org/blog/2015/10/glvisualize/ in your today's post in Julia users and will give a look in FileIO to learn what needs to be done to integrate GMT in it, but we have also work to do in GMT itself so things will take time.

And, beside sliders, what about other type of buttons like pushbuttons, checkboxes, radiobuttons, etc. will it be possible to have them too in GLVisualize?

[joa-quim]

Ah, I forgot to reply what is big grid. Mirone being 32 bits plus the RAM fragmentation means that something bigger than ~400 MB starts to be in within its limits. But that is big enough to slow down even fledermaus.

This one is about 100 MB http://w3.ualg.pt/~jluis/misc/tejo10_geo.grd.zip and should give already a good idea of a relatively big file (but real big files can go over GBs or more). The extension is .grd but the file is netCDF too and GMT opens it fine.

[SimonDanisch]

Okay I made a little test and I can easily display a 10.000x10.000 grid with Float32 (400mb) on my crappy onboard gpu, as long as I adapt the resolution of the mesh I produce (it's produced on the GPU anyways, so it's very fast to change it's resolution). It will need some work to adapt the resolution view dependant, so that you can zoom in and actually see more. I will not start with that work before I haven't merged next2, just letting you know the possibilities. Well, you can display arbitrary objects and build signals that depend on clicking them, which means you can build arbitrary complex buttons. So it's just a matter of implementing them.

[SimonDanisch]

I can't download that file, by the way!

will give a look in FileIO to learn what needs to be done to integrate GMT in it

Should be relatively straight forward, but I guess there is no rush ;) Just implement load/save and register the file formats that GMT can load! I think the docs in FileIO should explain that. If not, don't hesitate to ask.

By the way, do you have some other interesting data, like magnetic fields, tsunamie data or what not? I'd like to have some better examples for this kind of data. Would be nice to have some terrain surface + some additional data overlayed, via lines, arrows or meshes.

[joa-quim]

Thanks for letting me know these news. It's interesting the idea that we will be able to create our own buttons (with all its pros and cons). Please continue your work on next2, I'm counting also on the documenting part of it.

[joa-quim]

Oops, sorry for the bad link (something happened and a sub-dir was moved to another location). The correct one is: http://w3.ualg.pt/%7Ejluis/mirone/misc/tejo10_geo.grd.zip I have some others at http://w3.ualg.pt/~jluis/mirone/data-links.html namely, the http://w3.ualg.pt/%7Ejluis/mirone/misc/EMAG2.grd.bz2 has the magnetic anomalies of the world. Regarding the tsunami data, that is more tricky. What I do is either output specific time frames of the simulation in individual files, or in a single 3D netCDF file. But to get nice pictures one must later combine the water height (the waves) with land models and because the two sets have very different amplitudes (~1 m for wave and hundreds to thousands meters for land) the image construction must be done separately. In Mirone I have dedicated tool to do just this. A quick (not so nice) example from data that I happen to have here now, We can work out an example with aerial image draped on a terrain model (Mirone does that and sends the result to IView3D) About FileIO, I didn't have much luck. Besides not understanding how I would implement the data loader (the registering looks simple) I had only errors. For example the bathy_1m.nc file is in netCDF4 and in this version there was a merge with HDF so FileIO should be able to read it with the HDF package, but not

julia> file = query("bathy_1m.nc")
FileIO.File{FileIO.DataFormat{:HDF5}}("bathy_1m.nc")
julia> file = load("bathy_1m.nc")
WARNING: UndefVarError(:load)
 in load at C:\j\.julia\v0.4\FileIO\src\loadsave.jl:77
 in load at C:\j\.julia\v0.4\FileIO\src\loadsave.jl:42
ERROR: UndefVarError: load not defined
 in load at C:\j\.julia\v0.4\FileIO\src\loadsave.jl:77
 in load at C:\j\.julia\v0.4\FileIO\src\loadsave.jl:42

also reading classical formats error for me (note: I'm using my previously installed ImageMagick)

julia> file = load("Capture.png");
WARNING: ErrorException("no decode delegate for this image format `Capture.png' @ error/constitute.c/ReadImage/555")
 in error at C:\j\.julia\v0.4\ImageMagick\src\libmagickwand.jl:149
 in readimage at C:\j\.julia\v0.4\ImageMagick\src\libmagickwand.jl:228
ERROR: no decode delegate for this image format `Capture.png' @ error/constitute.c/ReadImage/555
 in error at C:\j\.julia\v0.4\ImageMagick\src\libmagickwand.jl:149
 in readimage at C:\j\.julia\v0.4\ImageMagick\src\libmagickwand.jl:228

[joa-quim]

and still on tsunamis. I knew I'd seen this before. It would be lovely to be able to do this https://software.intel.com/en-us/articles/3d-fluid-simulation-sample

[SimonDanisch]

Oh yeah, that's exactly what I'm after in the long run! I see if I can use this to show off a first OpenCL/OpenGL interop example for GLVisualize. On the first look it seems to be a pretty big code basis, but it's maybe just noise and I can reduce it to the OpenCL kernel.

[SimonDanisch]

Best would be, if you could give me some example code were you extract the data and how you would like to visualize it. For magnetic fields I assume you would use vectorfields, like this one: This is for dense 3D arrays of velocities. With the next2 branch, position + rotation and all sorts of attribute combinations would be also possible.

I see you visualize seismic data as particles (http://w3.ualg.pt/~jluis/mirone/images/FlederSismic.jpg), which should be really straight forward with GLVisualize.

[SimonDanisch]

julia> file = load("bathy_1m.nc")
WARNING: UndefVarError(:load)

I can reproduce this. This must be a bug then, I'll see what I can do!

[SimonDanisch]

(note: I'm using my previously installed ImageMagick)

ImageMagick seems to have problems with that. I'm still not sure how to resolve that. It should work when you install it with build...

[SimonDanisch]

Uh, I can't find a load function in HDF5.jl. I just assumed @timholy put it there since he registered HDF5 to FileIO. @timholy, is HDF5 not supported yet, or is it working in some mysterious way, which broke at some point?

[joa-quim]

For the magnetic case. I'm the type of Geophysicist that is more interested in the magnetic anomaly intensity (the vector modulus). That is so because the magnetic fields originated in the crust are much smaller than the that of the Earth core and what is really useful is the projection of the crust originated magnetic field on the direction of core field. That said, for people that work on the nucleus things are the other way around (crust fields are noise for them). But it would be a nice example to have an Earth core vector field example. I should be able to generate that data with GMT alone.

The earthquakes, yes they are represented with particles (size proportional to magnitude). The thing to address is, again, that x,y coordinates are angles (long,lat) and z in meters (km) so some vertical scaling must be done and one that does not completely smash the terrain surface if we want to see both.

[joa-quim]

Hi Simon, One more puppy. If you install this LIDAR reader https://github.com/joa-quim/Laszip.jl and, for example, download this file http://www.liblas.org/samples/Mount%20St%20Helens%20Nov%2020%202004.laz than you can load it with

xyz = laz2xyz("Mount St Helens Nov 20 2004.laz");

and I'm sure it will be beautiful when plotted as a point cloud. What I do in Mirone is to color the points in function of their height.

[SimonDanisch]

You mean sort of like this:

I'm positively surprised, that GLVisualize is able to display this.. This is the anti-aliased circle kind of type, which I thought would be way to slow for that many particles...

[SimonDanisch]

This is with the next2 branch...Also I need to find a way, to allow people to define arbitrary functions for the color.

[joa-quim]

Very cool. You may compare the speed with ivew3d in Mirone. Just drag the laz file on Mirone and than bush the "View with Fledermaus" button on the fig that pops out. You may need to update Mirone ("Help -> Check for updates") because I recently made I change to use the internally shipped iview3d.

[SimonDanisch]

The performance of Fledermaus seems worse, although I think I'm doing quite a bit less... So it's not that easy to compare ;)

[SimonDanisch]

I installed iView4 since the update run into an infinite loop it seems. Is iView4 supposed to be slower?

[joa-quim]

Updating should not go into a loop. What it does is to compare your current files with a list of updated files. When the check sum differs the new one is downloaded. There is currently one permanent failure because the IT guys of my University are not able to allow downloading a file without extensions. So, repeating the "Check for updates" will keep downloading one particular file, but you can ignore this.

Regarding the iView4D, it's not that it's slower (it's not) but the problem is that it expects another format (v7 instead of v6) of the ancient Fledermaus (now QPS) format. Add to this that the format is not open and what works is because I managed to hack the binary format. Summary, new v7 still works but complains a lot of missing projection info. That is why I'm forcing to use the old viewer that has the additional bonus of weighting ~5 MB instead of ~200 for iView4D.

[joa-quim]

Hi Simon,

How do we build an image like that point cloud of St Helen's Mount? I tried to use issue #180 as a first guess but fro me it doesn't work. All I get is an empty window.

[SimonDanisch]

You might need a GLAbstraction.center!(w)! From what I remember, your values where pretty out of range for the default camera settings.

[joa-quim]

Well, I believe I need that ... and the rest of the script. Sorry, Simon, but without a manual I'm not able to start using GlViz