use mesh for custom walls

[stephanschulz]

Hello and thanks for sharing your code with the public.

I am usually working with openframeworks.cc which is a c++ toolkit for artists. It has a library that makes interfacing with voro++ easy. Here a link: https://github.com/patriciogonzalezvivo/ofxVoro

I am trying to extend it, to allow people to use a 3D mesh as a custom wall object. I am following the original torus example: http://math.lbl.gov/voro++/examples/torus/ but had no luck too far. I guess I am not understanding what happens in this function: https://github.com/chr1shr/voro/blob/master/examples/walls/torus.cc#L57

I adapted 'bool point_inside(double x,double y,double z) ' to successfully checking if a point is inside or outside my 3D mesh structure. But inline bool cut_cell_base(vc_class &c,double x,double y,double z) is a mystery to me.

Here an image showing that I can use ray tracing to:

• check if a point is inside/outisde
• find the intersection point of a ray going from double x,double y,double z to the origin of the 3D mesh. It looked to me that 'bool point_inside(double x,double y,double z) ' generates a normal vector. But I am not sure from where to were and if the magnitude is important.

I was hoping you have some advice. Thank you.

[chr1shr]

When the Voronoi cell for a particular point is being constructed, the code works in a coordinate system where the origin is centered on that point.

In the basic Voronoi cell construction routine, the code computes displacement vectors (x,y,z) from a particle to its neighbor. The call c.nplane(x,y,z,w_id) cuts the Voronoi cell by the plane

v . (x,y,z) = 0.5*rsq

where v is a 3D position vector, the "." is a scalar product, and rsq=x*x+y*y+z*z. This is the perpendicular bisector between this particle and its neighbor.

In addition, the other version of the function, c.nplane(x,y,z,rsq,w_id) allows for rsq to be specified independently.

For the curved walls like in the torus example, the code computes the closest point on the wall from the particle under consideration. It applies a plane cut to the Voronoi cell based on the tangent plane of the wall at this closest point.

I think it's feasible to adapt the code to handle the case that you mention, but there are a few issues:

• You'll probably need to compute tangent planes at each mesh point. You could look at neighboring mesh points and fit a plane through them.
• Currently the voronoicell class only handles convex polyhedra. This is usually the case for Voronoi cells, but this representation has limitations on accuracy for non-convex domains like the hyperboloid in your example. It can work okay in some situations, and the torus is one example of a non-convex domain. But, if your domain has sharp non-convex internal angles (e.g. a V-shaped domain) then you may get some strange results near those positions.

A full solution for non-convex domains is actually rather complicated, both mathematically and computationally. To illustrate this consider the following 2D example:

Suppose you were trying to compute Voronoi cells inside the red boundary. Then the Voronoi cell for point p would actually be split into two, unless you altered the definition of the Voronoi cell in some way. Nevertheless, you can still come up with good approximations in many cases.

[stephanschulz]

Thank you for taking the time to answer my question, it is very much appreciated.

I am using a ray tracing libraryto find the intersection on the mesh, when drawing a ray from the cell center to the 3d object's center. This also allows me to know if a point is inside or outside of the 3d Mesh object. The library also returns the normal vector at this intersection. I guess this should then allow me to get the perpendicular at that point and pass those x,yz in to c.nplane(x,y,z,w_id) https://github.com/chr1shr/voro/blob/master/examples/walls/torus.cc#L70

I will give it a try and report back.

[stephanschulz]

Thanks again for your input. I am pretty certain i am now correctly generating a tangent vector to the intersection point and it's normal vector. As mentioned I get the intersection's normal vector as a return value from the ray tracing library.

To insure that it is really a perpendicular vector or the normal vector and that all other values are correct i draw a couple of the results. In the video you see 2 planes with their normal and in red tangent vectors. But simply feeding the tangent vector in to return c.nplane(tangentNorm.x,tangentNorm.y,tangentNorm.z,w_id); does not seem to effect the cells.

I am also noticing that the x,y,z point that is feed in to inline bool cut_cell_base(vc_class &c,double x,double y,double z) is not necessary a cell center or cell vertices. I compared these x,y,z values with all cell center or cell vertices and do not get a match. Not sure what "particles" it uses to trigger those functions.

https://user-images.githubusercontent.com/1054816/102909462-fbedd780-4446-11eb-9ee8-a5e952383030.mov For now i switch to a non concave 3d mesh to avoid the problems you discussed earlier.

[chr1shr]

It's not sufficient to call nplane with the components of the tangent vector. The displacement of the plane also needs to be calculated.

Suppose n=(nx,ny,nz) is your normal vector and it's normalized so that ||n||=1. It should also point outward. Let p be the position vector of the particle, and c be the position of the closest point on the surface. Then, if v is a vector in the coordinate system relative to the particle, the plane is specified by

v . n = n . (c-p)

If you compare this to the form listed in my previous post, you'll see that you need to call nplane with (nx,ny,nz,rsq,w_id) where rsq=2* (n.(c-p)).

There is some flexibility with scaling. For example, you'll get the same result if you pass in (2*nx,2*ny,2*nz,2*rsq,w_id). Thus it may not be necessary to normalize n.

[stephanschulz]

Thanks again. I translated your suggestions in to openframeworks' way of using c++.

  ofVec3f p = ofVec3f(x,y,z);
        ofVec3f cc = ofVec3f(pointOnSurface);
        ofVec3f n = ofVec3f(normalOnSurface).getNormalized() * -1; //tangentNorm
        float rsq = 2 * ( n.dot(cc - p) );

        cout<<"----"<<endl;
        cout<<"p "<<p<<endl;
        cout<<"cc "<<cc<<endl;
        cout<<"n "<<n<<endl;
        cout<<"rsq "<<rsq<<endl;

        return c.nplane(n.x,n.y,n.z,rsq,w_id);

It seems to produce something almost correctly. I reduced the number of cells to see more clearly what's what. I can see the cells being cut at the surface plane (yeah) but the cells still overshoot at the other side of the 3d mesh. ???

https://user-images.githubusercontent.com/1054816/102932484-5fd7c680-446e-11eb-9278-5446eca33568.mov

[stephanschulz]

i guess i am doing some right since the cells closer to the center of this hourglass object work pretty well.

https://user-images.githubusercontent.com/1054816/102933539-69622e00-4470-11eb-9acc-19abdb287a4f.mov

[chr1shr]

This looks like it's doing the right thing, although since the surface is only approximated with single plane cut for each point, there's a limit to its accuracy.

You could try cutting the Voronoi cells by multiple tangent planes on the surface (and not just the tangent plane at the closest point). Within the cut_cell_base function you can make multiple nplane calls.

For a convex domain, you can get arbitrarily accurate by applying more tangent planes. But for a non-convex domain, you might find that applying too many tangent planes will remove parts of the Voronoi cell that you want to keep.

[stephanschulz]

Within the cut_cell_base function you can make multiple nplane calls.

Interesting. I guess those extra plane cuts should be done with cell points on the other side of 3D mesh. Maybe somehow get the vertices that make up a the cell's faces. Not sure how to loop through different faces. But there is a function that seems to then let me extract all x,y,z points per face. https://github.com/chr1shr/voro/blob/master/src/cell.cc#L2221

Then i could get a point on that face; maybe the centroid and check it with the 3d mesh surface as i have been doing so far.

Do you know which function i need to call to get a point to cell's faces, or how to get the differences vertices per face?

BTW happy holidays :)

[stephanschulz]

I will give this a try. https://stackoverflow.com/a/43664448 http://math.lbl.gov/voro++/examples/polygons/ where line 50 shows how to "Loop over all faces of the Voronoi cell"

Sorry i am noticing i am using this thread a bit like a note book. I hope you are ok with this.

[xarthurx]

@chr1shr Is the neighbor tracking really enabled? https://github.com/chr1shr/voro/blob/46db70fec291439178401a3c06fe32177578bbd4/src/cell.hh#L340

https://github.com/chr1shr/voro/blob/46db70fec291439178401a3c06fe32177578bbd4/src/cell.hh#L351

From the source code, it seems even neighbour id is taken, it wasn't really used, is it?

[chr1shr]

The code is correct, and this is not a bug. There are two variant classes for computing Voronoi cells. The voronoicell class does not do neighbor tracking, and it uses simplified data structures. The voronoicell_neighbor class does neighbor tracking and has additional data structures for holding this information.

Since the voronoicell class does not do neighbor tracking, the p_id argument to nplane is ignored. Thus, in these inline functions one might as well set it to zero rather than pass in a value.

If you look at the corresponding functions in voronoicell_neighbor you will see that p_id is passed to nplane.

[xarthurx]

@chr1shr Thanks for the quick reply. I'm trying to implement the same functionality as this issue addressed, and some of the source code confused me a bit.

 52:                 // This template takes a reference to a voronoicell or
 53:                 // voronoicell_neighbor object for a particle at a vector
 54:                 // (x,y,z), and makes a plane cut to to the object to account
 55:                 // for the toroidal wall
 56:                 template<class vc_class>
 57:                 inline bool cut_cell_base(vc_class &c,double x,double y,double z) {
 58:                         double orad=sqrt(x*x+y*y);
 59:                         double odis=orad-mjr;
 60:                         double ot=odis*odis+z*z;
 61: 
 62:                         // Unless the particle is within 1% of the major
 63:                         // radius, then a plane cut is made
 64:                         if(ot>0.01*mnr) {
 65:                                 ot=2*mnr/sqrt(ot)-2;
 66:                                 z*=ot;
 67:                                 odis*=ot/orad;
 68:                                 x*=odis;
 69:                                 y*=odis;
 70:                                 return c.nplane(x,y,z,w_id);
 71:                         }
 72:                         return true;
 73:                 }
 74: 

I copied the code from the torus example above for reference.

1. The term particle refers to the center of voronoi cells, correct? I'm a bit confused at the beginning about the use of particle, vertex, particle at vector (x,y,z), point, etc.
2. Any reason why the voronoicell class doesn't contain the its center? but we need to pass both the voronoicell class as well as its particle location (x,y,z) to the cut_cell_base class?
3. In the above code, it seems you actually MODIFY the value of x, y, z. Is this the transformation from the global coordinate to the local coordinate? If not, what does it actually do? Modifying variables inside a function in this way is a bit strange and hard to understand... Could you please explain a bit more?

P.S. The current method for container initialization, wall initialization, etc, is based [0,0,0] as the coordinate origin. This is fine for the examples shown on the website, even for the torus example, as the walls are all parametrically defined. However, for the customized wall using a mesh (vertices defined in arbitrary global positions), the user needs to translate it to the original position, or compute a bounding box of that mesh and initialize a container there. It might be helpful to consider in the future upgrade of this circumstance.

[xarthurx]

OK, some updates here. Hope these demos proofs of using general mesh as the bounding wall to generate Voronoi data structures. The demos below show how different types of meshes (convex, non-convex) generate different results.

Experiments

For all the results, the point sets include the offset barycentric centres of the triangles from the mesh (as the outer layer), and a regular grid (as the inside filling) with different densities. Each cell is cut by three planes, derived from the three vertex normals of the triangle which the closest point of the cell centre locates (We can increase the cutting plane of course, but for these examples, this is enough).

Sphere

Bunny

Original

simplified

Human Body (very dense)

Thoughts

As you can see from the images, for sharp non-convex areas (bunny ears, human hands, feet, etc.) the result is not following the geometry, this is highly related to the number of cutting planes, which can be fixed by adding more planes. But it is good enough for my own use. Adding more planes will also increase the computation speed (below).

Advice for the future development team

Only the cells that are in the outer layer needs to be cut by the custom wall. If the library allows the user to define the particles cut/uncut by the wall, it will speed up a lot when both the point set and the # of walls are large (O(MxN) at the moment, where M ~ # of tri in the mesh bound, and N ~ # of particles in the point set)

Though I can create a custom hash map inside the wall class for that purpose, I would recommend the implementation on the library side.