okey dokey.

added two methods for Unknown. This is probably unsafe (information destroying, not just unable to solve) for nonlinear or multivariable equations, but should fit the purpose at hand okay.


Code:


then


Code:

Sadly, I'm hoping to avoid scipy v0.11.x if possible. I hope that doesn't make things too challenging, as I really appreciate your algorithmic help on this. I tried doing some tricks with flipping around axis choices when I detect a bad choice, but that didn't help.

Edit: I will try this new method and report back.

Edit #2: Here are my first results with epsilon == 1e-4:


Code:

look at the new code. it's scipy free.

[edit]
patch to make collect_additive_terms free "safe", even if it still can't solve everything.


Code:


New behavior (good):

Code:


Old behavior (bad):

Code:


But doesn't matter if you only have 1 variable


[edit2]
to change the broadcasting problem:

change

Code:


to

Code:

Sadness; the problem remains.


Code:

What version of numpy are you running? I'm not using anything fancier than http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#advanced-indexing

If all else fails, you can convert that to a series of 4 sequential assigns.

[edit]

Also, I've uploaded Unknown.py to
https://github.com/elfprince13/ModelDecomposition

I'll probably just do the assigns. By the way, don't I need to collapse 4 times, once each for each of the lambdas and once for the unknown Cartesian coordinate? Also, I am thinking of doing this process only if projection yields a degenerate triangle, for the sake of a shorter code path. Finally, I'll be opening a repo for this code soon.

Nope! Just collapse once to be able to solve the m*z+b=0 for z=-b/m, the __call__ automatically collapses all variables passed to it.


And it might be worthwhile to use timeit on a batch of these for profiling purposes to compare our two codepaths - I still suspect that this will be performant against your python version because of the improved used of vectorization.

Incidentally, double check the end of the big chunk of linear algebra code. I made a copy + paste error and only the second of these two lines is necessary:

Code:



As far as repos go, it might be worthwhile to maintain a shared codebase (with separate forks, but same project history) on this project in case we want to merge our projects down the road, given the closely related endgoals.

Hmm, I'm confused about that. It doesn't seem to me that that would help the lambda1, lambda2, and lambda3 = 1. - lambda1 - lambda2 being represented with Unknowns. And outputting the final rval for "successful" solutions seems to confirm this:

Code:

I also notice that every solution is either -524.6335275 or 0.0, which worries me, considering that neither of those are valid coordinates for this model.

Are you sure you're __call__ing properly in the code where you replaced the vectorized assign that wasn't working?
This:

Code:

Replaces the unknowns just fine.

Also, can you maybe pastebin your current code/some input matrices/coordinate sets so I can debug?


[edit]

Looked at pastebin.

Code:


Should be:

Code:


The < epsilon bit was a typo on my part, but the rval thing is the proper replacement for the

Code:

or

Code:

that weren't working

[edit 2]
You dropped your line defining Pplanar somewhere along the way. It should be erroring at you unless you're running it interactively and it's now a constant global.

Nope, that line is there. And I have the proper solved and real_z stuff there now. :/ I'm still getting huge numbers of unsolveable things that seem to actually match the degeneracies I was seeing. In fact, that's pretty logical to me, considering I'm still passing your version of this function unknown variables that could be anywhere along a line in-plane with the triangle, essentially. I figured out the condition that causes this (below) and will try to do some axis-flipping before I even get to your or my version of this function:
Code:
(tv is tri)

I don't think axis flipping will help a real degenerate triangle, because linear algebra (a singular matrix is a singular matrix), and I'm quite surprised that mine is having trouble with any correct geometry, since all of the transformations should be isometries.

Can you do a full run overnight and see if the output from mine matches or is more or less complete from your version of the function? It should at least resolve issues where the triangle projection was producing a line.

Analytically/geometrically speaking, can you describe the condition that is causing problems? My spatial reasoning/mental code interpretation module has switched off for the night

The problem isn't the tri; all of the tris that we're giving the geom_cart2bary() in any of its forms are fine. The problem is the way I iterate over triangles. I measure the triangle's projection onto each axis, and iterate over the two longest axes, to make sure the generated voxel patch will have no holes. For each point along the plane generated by iterating along those two axes, it asks geom_cart2bary() to find the third coordinate on the triangle. It then generates a voxel there.

However, there are rare triangles that have lines parallel to one of the three axes, AND when I pick the long axes, the projection of the triangle onto the plane I pick is a line. When I give your function or mine points in that plane, they do not intersect the triangle's plane at all! Which means everything goes squiffy.



In fact, doing axis flipping DOES help, as I demonstrate in this screenshot. Previously-missing regions turn into proper solids or slatted solids. The slats are because it's not iterating with a high enough resolution in the new axis chosen. I plan to find a way to resolve that.



I can.

See above.

Maybe I'm just misinterpreting the process, but it seems like it would make more sense to iterate in triangle-coordinate land, and project afterwards? That way you know your target surface will always be a square with a particular resolution, rather than a line or arbitrarily small triangle?

Yeah, that would work very nicely indeed, and would probably remove a "thin triangle" issue I've been struggling with, wherein certain very thin triangles disappear entirely, leading features made up of thin triangles to not show up. How would you recommend that I go about that?

What about something like this:

For a triangle with edges of length l1, l2, l3. subdivide l1, l2, and l3 each into something like m= ceil(ln / voxel_side_length) line segments (i.e m+1 vertices), and generate a "grid" (or rather several overlapping grids) in barycentric coordinates, and for each point on the grid, calculate the interior normal vector.

Any voxel containing 40% or more of the interior normal vector gets filled.

Select output color/texture using a weighted function that counts based on the total length of interior nomal vector contained by that voxel from a given source texture.

What does the internal normal vector mean in the context of a planar triangle (or just a triangle, since all triangles are planar)?

Counterproposal: Pick the two longest sides as the "axes" for iteration. Iterate on voxel_side/2 steps or so, picking lambda l1 and l2 based on that "grid" (of course, it's not necessarily a square grid). Generate l3, x, y, and z from that. If the point is inside the triangle, rasterize to the voxel grid and fill the corresponding voxel. Use current method of scanning over a quadrilateral section of the texture to compute the average color inside that quad section of the triangle.

In virtually all computer graphics software, vertex order specifies which side of the triangle is the front and which is the back, and then back face can then be safely not drawn. The normal vector pointing out from the "in/back" side of the triangle is the "interior normal vector" for a closed mesh.

I worry that iterating the barycentric "grid" over only 2 axes will lead to empty spaces on triangles where all 3 sides are relatively close to the same length.

Filling based on normals rather than faces also should give a better approximation of the surface, because you're filling voxels that would be significantly occupied by a skin of a certain thickness, rather than simply anywhere a triangle passes through a teensy corner of a cube.

Example: Blue circle is filling by normal/volume, red is filling by surface.


You can of course tune the ratios needed to fill a voxel, and even have different ratios depending on whether you're looking at the base or tip of the normal vector (if you haven't filled yet, accept 30% coverage to prevent a hole, but prefer 60%+ near the outside).

That all makes perfect sense to me, thanks for the great (and diagrammatic) explanation. So you're suggesting iterating over three axes, correct? Projecting lines perpendicular to each side across the triangle, and then generating a normal at every point of intersection? Or more specifically:

1) Iterate L1 over the length of one side
2) Iterate L2 over the length of another side
3) Iterate L3 over the length of the last side
4) For every L1L2, L2L3, or L1L3 intersection:
4a) Compute the normal (cross the vectors from the edges to the points?)
4b) Figure out what voxel(s) that normal lies inside
4c) If >x% of that normal lies inside any voxel, fill the voxel based on a patch of the texture on the triangle surface

Edit: Just for fun, the Metropolitan Museum. Note the missing triangles.

Yes So in terms of complexity it's really more like iterating over 2 axes 3 times than it is like iterating over 3 axes



Models typically store vertices in counter-clockwise order, so I *think*
given Cartesian coordinate vectors for v1,v2,v3, the direction of the exterior normal vector is given by (v2-v1) x (v3 - v2), which means the interior normal vector is given by (v3 - v2) x (v2 - v1) (this is essentially just the right-hand rule as applied to solenoids and magnetic fields). Obviously, you can double check this if the results seem fishy.


Yep! Though since voxels may have many different normals extending into them, you may want to do some kind of weighted decision making to choose the textures.

Very cool!

That all makes sense. I got stuck when I got to the point where I'll be taking points on the edges of edges, casting them into (or out of!) the triangle perpendicular to those edges, and finding the points of intersection. Any insights?


Code: