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