IsoSurface plots for Wavefunctions on specific grid

[louisponet]

Hi,

I'm trying to visualize some wavefunctions inside a crystal, and other related things. Now the issue I'm running into is that these wavefunctions are defined on a grid, which is spanned by the 3 primary vectors of the crystal. So for example the primary vectors I'm working with are: a=(2.41,0.0, 3.54),b=(-1.21, 2.09,3.54),c=(-1.21,-2.09,3.54) in cartesian coordinates, which results in the wavefunctions being defined as a 81x81x81 grid of values, inside the 3D space spanned by these vectors.

Now I tried running through the code for the visualize(vol,:iso,isovalue=...) function, but I get lost somewhere after the definitions of hull::AABB and other data parts (namely, how it progresses afterwards to define the points of the vertices in the mesh). Is there an easy way to get to some kind of functionality where I could also specify the vectors along which my volume data is defined to the visualize() function, or manually specify them in this data definition?

Ultimately the problem I'm working on comes down to analyzing how the dipole of an energy band of the bandstructure of a crystal moves around while the momentum wavevector changes. It would be awesome to be able to make some animation showing the way the electron cloud moves around, and it's effect on the bandstructure and other properties of the material.

There is also some idea brewing in my mind to start a side project writing some kind of crystal visualization package for julia, similar to the program VESTA which is what our group is using now. Your package would be an awesome backbone for this, it basically already has everything that would be required (if I get the above to work)!

grtz

[SimonDanisch]

Sounds like a great project! If you figure out how to represent the transformation from the unit vector to your axis vectors and express it in a 4D matrix, you should be able to just pass that as the model (matrix) argument! Too tired to figure out the right formular for that myself right now ;)
Let me know how it goes!

[louisponet]

Ok, seems like I figured out the correct model transform matrix. I verified it by taking one of your volume examples and it looks like it's the same, but with the box now in the shape of my crystal unit cell. This is the matrix: Mat4f0([2.41004 -1.20502 -1.20502 0.0;0.0 2.08716 -2.08716 0.0;3.53817 3.53817 3.53817 0.0;0.0 0.0 0.0 1.0]).

Edit: it works now, normalizing the wavefunction mesh and taking isorange 0.2f0 made things better (I had issues with holes appearing etc)!
Edit2: weird left/right mouse movement I mentioned is related to the upvector of the camera not changing, which is fixed by the issue below

[louisponet]

Hi Simon,

I've been making a bit of progress on getting some isosurface and crystal structure viewer to work. This is where I'm at so far (it's not a bad start):

I was wondering about a couple of things.

I implemented some marching cubes style isosurface extractor, and the I create a GLNormalMesh from the resulting vertices. I did this because I didn't fully understand your isosurface implementation inside the volume shader, or at least I couldn't get it to work in this case. One of the issues I guess is that the data is both positive and negative (which I display with red and blue in this case), together with the fact that even when using the correct model transformation matrix, it still results in slices of white that run through the isosurface. I'm quite happy with the solution I have, however it's not indexed yet (I was bored of implementing the marching cubes for now), meaning that I'm just using a 1:length(vertices) array as indices. This is probably not very performant. Is there a mesh that works in the pure oldschool GLDrawArray(....,..,GLTriangles) way?

Secondly, do you have any pointers to where I can search for why there are these weird stripes appearing (in the picture I aimed it so that it was just right to see them, a bit more to one side and it's mostly whitish hue, a bit more to the other and it gets to a more primary colour hue)?
Together with this question, how exactly does the rendering pipeline work, and how can I influence it? For example, the two "atoms" are hyperspheres, if I _view(visualize(...)) their meshes after those of the isosurfaces they don't show at all, even though my isosurface is supposed to have a 0.6 alpha channel. The same situation is happening with my :linesegments that I use to mark my boundaries with. However, I can't seem to make these appear behind the other meshes no matter what I try. Any ideas?

If it's handy I can post the code I'm using, I don't want to overclutter things.

[SimonDanisch]

Hi, really sorry that I didn't answer at all...
The isovalue shader should be easy to adapt for what you need. If you can give me some code I can look at the artifacts - both for iso surface rendering and the one you created above!
For the rest of the issues, it's also probably easiest if you supply me with the code and I can see what the problem is!

[louisponet]

No worries about the reply, it's nothing game braking atm!
Ok, so it's a little bit tricky to give you minor code to work with.

The easiest thing to start with is I guess the artifacts above, which I thought might have come from the way normals are generated, but when I disabled them in my shaders nothing changed.
I've made this pastebin that gives a working example of the marching_cubes mesher that I implemented, but it looks horrendous due to me trying to optimize the inner loop for speed. The smallest possible working example of this system is given here: https://pastebin.com/0sKpb643. It should be a "select all and run" experience, which in the end should show you a sphere that is visualized. The marching cubes is just the classic algorithm from paul bourke ported to julia.

For getting the isosurface shader working in my situation might be a bit more tricky, since it involves the actual wavefunction files, a reader function for them. If that is not an issue I can compile another pastebin for you to copy-paste. Let me know!

EDIT I forgot to check again the issues that I talked about. This is code purely for replicating the lines in the picture.

[SimonDanisch]

Okay, I had a chance to take a look at this. Since the artifact appears only with transparent colors, I'm guessing that the mesh has overlapping vertices. When loading it into microsoft 3D builder, it also told me the mesh is broken...
Why did you store a shader into the Mesh? I'm surprised that it doesn't throw an error, maybe I should open an issue, that it should throw for unsupported attributes.

Did you realize, that Meshing.jl contains marching cubes & marching tetrahedra?
https://github.com/JuliaGeometry/Meshing.jl/blob/master/src/marching_cubes.jl

Best,
Simon

[louisponet]

So I didn't understand GLNormalColorMesh it seems, I thought it might've had a thing where all the extra keywords would get pushed through to the final RenderObject with visualize(..).

I didn't know about the meshing.jl! I tried using both it's marching_cubes and marching tetrahedra methods to render an iso mesh using this code:

using Meshing,GeometryTypes,ColorTypes,GLVisualize

sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
    sqrt(sum(dot(v,v))) - 1 # sphere
end

# meshing marching_cubes test 
m = marching_cubes(sdf,0.1)
m_test = GLNormalColorMesh(vertices=m.vertices,faces=m.faces,color=RGBA(1.0,0.0,0.0,0.6))

screen = glscreen()
_view(visualize(m_test),camera=:perspective)
renderloop(screen)

#Meshing marchingTetrahedra test
m = Meshing.GLNormalColorMesh(sdf)
m_test = GLNormalColorMesh(vertices=m.vertices,faces=m.faces,color=RGBA(1.0,0.0,0.0,0.6))
screen=glscreen()
_view(visualize(m_test),camera=:perspective)
renderloop(screen)

It seems that the weird lines also appears when I use these methods. I checked the meshes and for sure I have overlapping vertices, because how can I otherwise form GLTriangles. I don't think this should be an issue though? Could it be something related to how normals are generated when calling GLNormalColorMesh(..)?

Cheers

[SimonDanisch]

oh, i didn't meant vertices, i meant overlapping triangles :-D

[louisponet]

Nope, I checked it using my code (x.wfcs[1] is just some values+points), and running this small for loop

vertices_pos,vertices_neg = marching_cubes(x.wfcs[1],0.4);

triangles = [[vertices_pos[i],vertices_pos[i+1],vertices_pos[i+2]] for i=1:Int(length(vertices_pos)/3)]
for i=1:length(triangles)
  test = find(x->x==triangles[i],triangles)
  if length(test)>1
    println(i)
  end
end

this didn't result in anything printed, and also my vertices_pos/neg are not empty.

[SimonDanisch]

well that tests only for an exact overlap ;) do you have another 3d program in which you can test this? i currently dont have anything installed that lets me change the transparency :(

[louisponet]

Ok, that's true. No this is beyond my graphical skills, what program would I use for something like this? How can I load something from julia into an external program?

What I tried is coloring the Triangles red/blue one by one using:

using Meshing,GeometryTypes,ColorTypes,GLVisualize

sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
    sqrt(sum(dot(v,v))) - 1 # sphere
end

# meshing marching_cubes test 
m = marching_cubes(sdf,0.1)
it_colors = [[RGBA(1.0,0.0,0.0,0.6) for i=1:3];[RGBA(0.0,0.0,1.0,0.6) for i=1:3]]
m_test = GLNormalVertexcolorMesh(vertices=m.vertices,faces=m.faces,color=[it_colors[(i-1)%6+1] for i = 1:length(m.vertices)])
m_test.color
screen = glscreen()
_view(visualize(m_test),camera=:perspective)
renderloop(screen)

I'm not sure if there are really overlapping triangles.
What's interesting is that if you put the camera exactly at the point where the stuff happens, and then translate it using right mouse, each time a next "ring" passes from back to front the color darkens or the other way around.

If I put glEnable(GL_CULL_FACE) in the function startin from
https://github.com/JuliaGL/GLAbstraction.jl/blob/c66428f33f048a56c33c57c8cbaf745487ae9bbd/src/GLRenderObject.jl#L72

the lines disappear, it now only has a milky hue to it, however i'm not so sure if this is useful...

[SimonDanisch]

Good question, maybe Blender!

But GL_CULL_FACE is a good point! I've disabled it for stupid reasons, as I found out later... I was actually pretty sure that I enabled it on newer versions of GLVisualize. I will take a look!
I'm guessing that this will just hide the issue though. Anyways, transparency isn't really working correctly right now, but I wouldn't expect such artifacts...

[louisponet]

Ok I discovered renderdoc. Pretty awesome app for capturing glStuff ... what I did was this:

1. Download renderdoc (linux/windows)
2. create an executable:

#!/opt/julia0.6/julia
using Meshing,GeometryTypes,ColorTypes,GLVisualize

sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
    sqrt(sum(dot(v,v))) - 1 # sphere
end;

# meshing marching_cubes test 
m = marching_cubes(sdf,0.1);
it_colors = [[RGBA(1.0,0.0,0.0,0.6) for i=1:3];[RGBA(0.0,0.0,1.0,0.6) for i=1:3]];
center = sum(m.vertices)/length(m.vertices);
normals = [normalize(center-vert) for vert in m.vertices];
m_test = GLNormalVertexcolorMesh(vertices=m.vertices,faces=m.faces,color=[it_colors[(i-1)%6+1] for i = 1:length(m.vertices)],normals=normals);

screen = glscreen()
_view(visualize(m_test),camera=:perspective)
renderloop(screen)

change to your julia exec ofc and chmod +x it.
3. open renderdoc, path to executable -> what you just created
4. run, wait until it shows a window, rotate till you see the weirdness
5. trigger capture
6. close window, double click captured thing.
7. open colour pass 1 -> glDrawElements with a large nr in the parentheses.
8. open tab "mesh output"
9. check weird vertices.... (Around vidx 700 things go weird, normals are often NaN as well)

In my implementation I don't have any rogue vertices or NaN normals. If you want I can send you a log file.

cheers

[louisponet]

Turns out it has probably nothing to do with the meshes. I guess I'll have to write my own rendering to get this sorted out (if I want to).

Anyway thanks for the help!

cheers

[SimonDanisch]

Would be awesome if we could just fix it :D I do actually still have a branch of glvisualize with an approximated order independent transparency rendering pass. I can see if I can merge that!
That's probably the most appropriate way of fixing this! Please keep me posted, I can help you if you want to take a stab at this!

[louisponet]

Ok, sounds good, if you merge it let me know!

I would love to try and figure it out, I'm slowly getting a better and better understanding of the whole framework of packages, but it's pretty hard. In any case I have a lot of things to do this month so I'm not sure how much I can focus on this. Probably November is when I will continue to work on my package and implement some more lowlevel rendering/api calls.

However if by then we figured it out from your end that would be even better :D! If I can do any testing or something, that I will have time for so let me know.

[louisponet]

Hi Simon it's been a while.

Did you find your branch with the correct transparancy rendering? In the meantime, in light of trying to get the isosurface shader to work for my case, I played around a bit today and I think this is a minimal working example that shows the issues.

using GLVisualize,GLWindow,GeometryTypes
using JLD

points = load("data.jld")["data"]
cell = Mat3f0([2.41004 -1.20502 -1.20502;0.0 2.08716 -2.08716;3.53817 3.53817 3.53817])
cell /= norm(cell)

window = glscreen()
vol = visualize(points,:iso,isovalue=0.12f0,isorange=0.02f0,model=Mat4f0([cell [0.0,0.0,0.0];[0.0 0.0 0.0 1.0]]))
_view(vol,window)
renderloop(window)

The data.jld file is found in the attached .zip archive.

As you can see there are some strange openings/bands in the isosurface. I also have to excessively finetune the isorange to get something that at least has some full surface areas. Any ideas on how to optimize this for my usecase?

For reference, this is how it is supposed to look (minus the atoms):

data.zip

[louisponet]

Hi Simon,

I've finally had some time to play around with this cheap tricks for OpenGL transparency thing. I've come up with this quick and dirty way of doing things (please, if there is a better place/way to do this give me some hints): In GLRender.jl, I changed the list render function to this:

 function render(list::Vector{RenderObject{Pre}}) where Pre
    isempty(list) && return nothing
    # first(list).prerenderfunction()
    first(list).prerenderfunction()
    vertexarray = first(list).vertexarray
    program = vertexarray.program
    glUseProgram(program.id)
    glBindVertexArray(vertexarray.id)
    for renderobject in list
        Bool(Reactive.value(renderobject.uniforms[:visible])) || continue # skip invisible
        # make sure we only bind new programs and vertexarray when it is actually
        # different from the previous one
        if renderobject.vertexarray != vertexarray
            vertexarray = renderobject.vertexarray
            if vertexarray.program != program
                program = renderobject.vertexarray.program
                glUseProgram(program.id)
            end
            glBindVertexArray(vertexarray.id)
        end
        for (key,value) in program.uniformloc
            if haskey(renderobject.uniforms, key)
                if length(value) == 1
                    gluniform(value[1], renderobject.uniforms[key])
                elseif length(value) == 2
                    gluniform(value[1], value[2], renderobject.uniforms[key])
                else
                    error("Uniform tuple too long: $(length(value))")
                end
            end
        end
        if haskey(program.uniformloc,:pass)
            pass_loc  = program.uniformloc[:pass][1]
            if haskey(renderobject.uniforms,:transparent)
                glBindVertexArray(vertexarray.id)
                
                glDisable(GL_CULL_FACE)
                glDepthFunc(GL_LESS)
                gluniform(pass_loc,Cint(1))
                renderobject.postrenderfunction()
                
                gluniform(pass_loc, Cint(2))
                glEnable(GL_CULL_FACE)
                glCullFace(GL_FRONT)
                glDepthFunc(GL_ALWAYS)
                renderobject.postrenderfunction()
                
                gluniform(pass_loc, Cint(3))
                glDepthFunc(GL_LEQUAL)
                renderobject.postrenderfunction()
                
                gluniform(pass_loc, Cint(2))
                glCullFace(GL_BACK)
                glDepthFunc(GL_ALWAYS)
                renderobject.postrenderfunction()
                
                gluniform(pass_loc, Cint(3))
                glCullFace(GL_BACK)
                glDisable(GL_CULL_FACE)
                glDepthFunc(GL_LEQUAL)
                renderobject.postrenderfunction()
            else
                gluniform(pass_loc, Cint(0))
                renderobject.postrenderfunction()
            end
        else
            renderobject.postrenderfunction()
        end
        # we need to assume, that we're done here, which is why
        # we need to bind VertexArray to 0.
        # Otherwise, every glBind(::GLBuffer) operation will be recorded into the state
        # of the currently bound vertexarray
    end
    glBindVertexArray(0)
    return
end

Then I added the following to the standard.frag shader :

uniform int pass;

void main(){
    vec4 color = get_color(color, o_uv);
    if (pass == 1){
        color = vec4(vec3(color.xyz),0*color[3]);
    }
    else if (pass==2){
        color = vec4(vec3(color.xyz),0.75*color[3]);
    }
    else if (pass==3){
        color = vec4(vec3(color.xyz), (color[3]-0.75*color[3])/(1.0-0.75*color[3]));
    }
    {{light_calc}}

    write2framebuffer(
        color,
        o_id
    );
}

To play around you can run the following after adding the lines (the data.jld is attached):

using GLVisualize,GLWindow,GeometryTypes,ColorTypes
using JLD

window = glscreen()
wfc_mesh = load("data.jld")["wfc_mesh"]
cell_mesh   = load("data.jld")["cell_mesh"]
_view(visualize(cell_mesh))
_view(visualize(wfc_mesh,transparent=true))
renderloop(window)

This seems to be working very well with regards to transparency. However there is some issue with the depth now, things that are supposed to be behind the transparent surfaces show up either always behind or always in front, depending on the rendering order. Any ideas why this is?

[louisponet]

Changing the last pass:

gluniform(pass_loc, Cint(3))
glCullFace(GL_BACK)
glDisable(GL_CULL_FACE)
glDepthFunc(GL_LEQUAL)
renderobject.postrenderfunction()

to

gluniform(pass_loc,Cint(0))
renderobject.postrenderfunction()

Seems to work perfect when you make sure to render the transparent things after the opaque ones.

[SimonDanisch]

This looks nice, although I suspect it to be more expensive compared to the OIT transparency pass I was using... Sorry, this is a bit of a low priority for me right now, until I have the higher level API figured out.

[louisponet]

Yes I understand, no problem. It probably is quite expensive because of the amount of render passes, I guess it depends on the usecase which is more suitable. Probably your OIT is better for general uses.

[SimonDanisch]

Btw, the OIT shaders live in: https://github.com/JuliaGL/GLWindow.jl/tree/sd/screencam

[louisponet]

Oh I didn't see that branch before! I'll have a look at it.

[louisponet]

Is the is_transparent_pass uniform used anywhere? I can't seem to locate it in any of the packages.

[SimonDanisch]

Not that I know of !