> faces/triangles in GTS are made from edges, so, an extra "layer".  agree re bug
> hunt, running your points through both implementations should provide useful
> info.
>

I found and fixed the kd-tree bug. So now it is correctly returning the nearest
vertex to the test point.

So my algorithm tests all triangles that use the nearest vertex, and then
calculates the minimum distance to each triangle, keeping the smallest of those.

It seems that this part is now not returning quite the correct minimum distance,
so I have to focus on that part next.

Post a reply to this message

> So my algorithm tests all triangles that use the nearest vertex, and then
> calculates the minimum distance to each triangle, keeping the smallest of those.
>

Ah, I think I figured out the problem.

Instead of testing the distance to all triangles that share the nearest vertex,
what I *should* be doing is finding a bounding box around the test point based
on the distance to the nearest vertex (probably plus some tolerance), then test
all triangles that overlap that bounding box (which I should be able to use the
bounding box tree to find quickly).

Post a reply to this message

On 7/13/22 3:28 PM, jceddy wrote:
> 
>> Either there is a bug in my kd-tree implementation, or there is someth
ing
>> fundamentally incorrect with the approach. Will have to do some more t
esting to
>> narrow down which it is.
>>
> 
> Looks like the kd-tree is not always returning the nearest mesh vertex 
to the
> test point. I create a test that compares the minimum distance found us
ing
> simulated annealing vs. the one found using the nearest-vertex test, an
d then
> added some code to find the nearest vertex by using a brute-force check
 of *all*
> the vertices in the mesh, which I run in the case the simulated-distanc
e found a
> better solution than the kd-tree nearest-vertex approach.
> 
> I fairly quickly ran into a case where kd-tree is not returning the cor
rect
> vertex, so now need to figure out why that is the case.
> 
     There is something fundamentally wrong with your approach: the 
closest approach is not necessarily on a face that includes the closest 
vertex. If you look at the image I posted to p.b.i, when searching for 
the closest point to the cross in the bottom, the closest vertex is the 
red one, but the closest face is the orange one, which doesn't include 
the red vertex…

         Jerome

PS: Sorry I sent that to your email instead of the group.
-- 
mailto:jeb### [at] freefr
http://jeberger.free.fr
Diaspora*: https://framasphere.org/u/jeberger

Post a reply to this message

>      There is something fundamentally wrong with your approach: the
> closest approach is not necessarily on a face that includes the closest
> vertex. If you look at the image I posted to p.b.i, when searching for
> the closest point to the cross in the bottom, the closest vertex is the
> red one, but the closest face is the orange one, which doesn't include

>

Yeah, I figured that out. :)

I am now finding the nearest vertex, then creating a bounding box around the
sample point that is large enough to contain the nearest vertex (plus a small
tolerance factor). Then I traverse the bounding box tree for the mesh and find
all leaves that overlap the sample bounding box, and calculate the minimum
distance for each leaf's triangle. The minimum of all of those is the minimum
distance to the mesh.

It performs fairly well (actually best performance of all my methods so far,
although still slows down quite a bit around detailed portions of the mesh where
it ends up sampling more triangles), and the output is great.

The output actually looks really good even if I turn max_gradient way down on
the isosurface (the actual max gradient is like 451, but I was turning it down
to 4 to get quick test renders).

Going to post some example outputs soon, then turn my attention to cleaning up
the code.

My current plan is to add a virtual method to object (probably called
"proximity") with a default behavior of ray sampling, or even just spitting out
a "not implemented for this object type" error, and then start overriding it for
the various object types.

Then will wire up an internal proximity function (I like proximity better than
"minimum distance"), and a proximity pattern type, as well, I think.

Oh, I also want to tinker around with doing an inverse translation on the sample
point and then translating the result vector back to measure it...I think it
might work the way I want now that I am not sampling the wrong set of triangles,
and if I can get it to work properly will allow me to remove a bunch of code I
added to Mesh to generate transformed versions of the vertex data and bounding
box tree.

Post a reply to this message

> Oh, I also want to tinker around with doing an inverse translation on the sample
> point and then translating the result vector back to measure it...I think it
> might work the way I want now that I am not sampling the wrong set of triangles,
> and if I can get it to work properly will allow me to remove a bunch of code I
> added to Mesh to generate transformed versions of the vertex data and bounding
> box tree.

I got the inverse-transform stuff to work correctly.

Cleaned up the code, removed all the simulated annealing stuff I had added,
changed "minimum_distance" to "proximity", implemented the Proximity() virtual
function in ObjectBase (which throws an error if it's called for objects that
don't have it implemented), implemented the Proximity() override for Mesh,
cleaned up pattern parser code for "proximity" pattern, and changed it to call
object->Proximity().

Proximity returns negated values for points inside the object.

The proximity pattern does not attempt to scale or clamp values, so it will look
like the values are just 0.0-1.0 repeating to the user...it is up to them to do
any required scaling that makes sense. This is something that we might want to
give some options for eventually, but it's outside the scope of what I am doing
ATM.

Updated my fork on github with all of the changes.

Next step (after some more testing) will be to implement Proximity() for other
object types.

Post a reply to this message

> Next step (after some more testing) will be to implement Proximity() for other
> object types.

Implemented Proximity() now for:

Mesh
Sphere
Box
Triangle
Plane
Disc
Cone (includes cylinder)

For non-solid objects like Triangle, Disc, Cone (if open), the proximity will
always be positive, since there is no well-defined "inside", so for testing
isosurfaces, I set threshold to 0.001 instead of 0.0, so it actually renders a
little thickness around the object.

I also implemented fn_nearest_point, that returns the nearest point calculated
by the Proximity() code.

So for proximity, you define the function like this:

#declare fn_X = function { proximity { shape } }

For nearest_point, you define the function like this:

#declare fn_N = function { nearest_point { shape } }

proximity functions return a double
nearest_point functions return a 3d vector position <x, y, z>

Attached a test render for an open cone, with some test points arrayed around
it.

Post a reply to this message

Attachments:
Download 'test_new.png' (91 KB)

Preview of image 'test_new.png'

> Implemented Proximity() now for:
>
> Mesh
> Sphere
> Box
> Triangle
> Plane
> Disc
> Cone (includes cylinder)

Found a bug in the bounding box tree traversal code I wrote that made it so
Mesh::Proximity() was checking any triangle whose BBox overlapped the sample
BBox in the x direction, fixing it to narrow down by y and z directions as well
sped it up considerably. XD

Anyone want to hazard a methodology for implementing Proximity() for Blob,
Bicubic Patch, or perhaps for any of the CSG shapes? (or any shape I haven't
listed above for that matter?)

Post a reply to this message

> Anyone want to hazard a methodology for implementing Proximity() for Blob,
> Bicubic Patch, or perhaps for any of the CSG shapes? (or any shape I haven't
> listed above for that matter?)

I have done a few more of the easier ones (like Torus, etc.).
I ended up implementing the Bicubic Patch one using similar estimation logic as
is used to create the triangle mesh approximation. It is doing proximity to all
the triangles in the mesh and returning the nearest point overall...could be
slow for very detailed patches.

Been thinking about Blob, but haven't worked up the guts to attempt it yet. It
seems like you should be able to come up with something you can feed to the root
solver, but I'm still not clear on how to do it.

Post a reply to this message

> Been thinking about Blob, but haven't worked up the guts to attempt it yet. It
> seems like you should be able to come up with something you can feed to the root
> solver, but I'm still not clear on how to do it.

For Height Field, I implemented something similar to what I did for Mesh, where
I put all of the points in a kd-tree, then I look at all points with x- or z-
distance within the distance of the sample point to the nearest height field
vertex, plus one, triangulate the field mesh defined by those points, and check
for nearest point on any of those triangles. (Unless, that is, the y-coordinate
of the sample point is above the max-y or below the min-y of the height field,
in which case I just use the nearest vertex as the nearest point on the height
field). I haven't handled water-level yet.

Here is the code:

DBL HField::Proximity(Vector3d &pointOnObject, const Vector3d &samplePoint,
TraceThreadData *threaddata) {
 if (Data->kdTree == nullptr) {
  return MAX_PROXIMITY_DISTANCE;
 }

 Vector3d transformedPoint = samplePoint;
 Vector3d diff;
 if (Trans != nullptr) {
  MInvTransPoint(transformedPoint, transformedPoint, Trans);
 }
 DBL px = transformedPoint[X];
 DBL py = transformedPoint[Y];
 DBL pz = transformedPoint[Z];

 Vector3d nearest = Data->kdTree->nearest_point(transformedPoint);

 if (py < Data->min_y || py > Data->max_y) {
  diff = nearest - transformedPoint;
 }
 else {
  DBL dist = (nearest - transformedPoint).length();
  DBL xmin = px - dist;
  DBL zmin = pz - dist;
  DBL xmax = px + dist;
  DBL zmax = pz + dist;
  int minX = max(0, int(xmin) - 1);
  int maxX = min(Data->max_x, int(xmax) + 2);
  int minZ = max(0, int(zmin) - 1);
  int maxZ = min(Data->max_z, int(zmax) + 2);

  dist = MAX_PROXIMITY_DISTANCE;
  for (int z = minZ; z < maxZ - 1; z++)
  {
   for (int x = minX; x < maxX - 1; x++)
   {
    Vector3d p1(x, Data->Map[z][x], z);
    Vector3d p2(x + 1, Data->Map[z][x + 1], z);
    Vector3d p3(x + 1, Data->Map[z + 1][x + 1], z + 1);
    Vector3d p4(x, Data->Map[z + 1][x], z + 1);

    Vector3d trianglePt;
    DBL test = NearestOnTriangle(trianglePt, transformedPoint, p1, p2, p3);
    if (test < dist) {
     dist = test;
     nearest[X] = trianglePt[X];
     nearest[Y] = trianglePt[Y];
     nearest[Z] = trianglePt[Z];
    }

    test = NearestOnTriangle(trianglePt, transformedPoint, p1, p3, p4);
    if (test < dist) {
     dist = test;
     nearest[X] = trianglePt[X];
     nearest[Y] = trianglePt[Y];
     nearest[Z] = trianglePt[Z];
    }
   }
  }

  diff = nearest - transformedPoint;
 }

 if (Trans != nullptr) {
  MTransDirection(diff, diff, Trans);
 }

 pointOnObject = samplePoint + diff;
 return diff.length();
}

I am getting really high max gradient values (resulting in isosurface holes)
around the highest and lowest points in the field, and I'm not clear on whether
that is due to a problem in the code, or maybe like due to the orientation of
the triangles relative to the camera ray there or something, so thought I would
ask if anyone could take a look at this and let me know if they see any glaring
issue.

Cheers,
Joe

Post a reply to this message

>  if (py < Data->min_y || py > Data->max_y) {
>   diff = nearest - transformedPoint;
>  }

After I posted this, I realized that this is not correct and might be causing
the issue, since it will be incorrect for sample points "above" or "below" flat
triangles at min_y or max_y.


I think I can still optimize for points above or below the y bounding planes. I
should only have to check triangles within one "step" of the nearest vertex.

So I modified the code thusly:

 if (py < Data->min_y || py > Data->max_y) {
  minX = int(nearest[X]) - 1;
  maxX = int(nearest[X]) + 2;
  minZ = int(nearest[Z]) - 1;
  maxZ = int(nearest[Z]) + 2;
 }

Here is the new version of the method I am currently testing:

DBL HField::Proximity(Vector3d &pointOnObject, const Vector3d &samplePoint,
TraceThreadData *threaddata) {
 if (Data->kdTree == nullptr) {
  return MAX_PROXIMITY_DISTANCE;
 }

 Vector3d transformedPoint = samplePoint;
 Vector3d diff;
 if (Trans != nullptr) {
  MInvTransPoint(transformedPoint, transformedPoint, Trans);
 }
 DBL px = transformedPoint[X];
 DBL py = transformedPoint[Y];
 DBL pz = transformedPoint[Z];

 Vector3d nearest = Data->kdTree->nearest_point(transformedPoint);

 DBL dist = (nearest - transformedPoint).length();
 DBL xmin = px - dist;
 DBL zmin = pz - dist;
 DBL xmax = px + dist;
 DBL zmax = pz + dist;
 int minX = max(0, int(xmin) - 1);
 int maxX = min(Data->max_x, int(xmax) + 2);
 int minZ = max(0, int(zmin) - 1);
 int maxZ = min(Data->max_z, int(zmax) + 2);

 if (py < Data->min_y || py > Data->max_y) {
  minX = int(nearest[X]) - 1;
  maxX = int(nearest[X]) + 2;
  minZ = int(nearest[Z]) - 1;
  maxZ = int(nearest[Z]) + 2;
 }

 dist = MAX_PROXIMITY_DISTANCE;
 for (int z = minZ; z < maxZ - 1; z++)
 {
  for (int x = minX; x < maxX - 1; x++)
  {
   Vector3d p1(x, Data->Map[z][x], z);
   Vector3d p2(x + 1, Data->Map[z][x + 1], z);
   Vector3d p3(x + 1, Data->Map[z + 1][x + 1], z + 1);
   Vector3d p4(x, Data->Map[z + 1][x], z + 1);

   Vector3d trianglePt;
   DBL test = NearestOnTriangle(trianglePt, transformedPoint, p1, p2, p3);
   if (test < dist) {
    dist = test;
    nearest[X] = trianglePt[X];
    nearest[Y] = trianglePt[Y];
    nearest[Z] = trianglePt[Z];
   }

   test = NearestOnTriangle(trianglePt, transformedPoint, p1, p3, p4);
   if (test < dist) {
    dist = test;
    nearest[X] = trianglePt[X];
    nearest[Y] = trianglePt[Y];
    nearest[Z] = trianglePt[Z];
   }
  }
 }

 diff = nearest - transformedPoint;

 if (Trans != nullptr) {
  MTransDirection(diff, diff, Trans);
 }

 pointOnObject = samplePoint + diff;
 return diff.length();
}

Sorry for the SPAM...sometimes just putting my code in front of other people
forces me to thing more deeply about it and I can figure out issues that I
hadn't previously.

Cheers,
Joe

Post a reply to this message