Details of subdivision algorithm.

[Leonid Raiz]

Hello,

Is there a publicly accessible document describing the recommended algorithm for subdividing curved triangles? I am familiar with the powerpoint presentation at amf.wikispaces.com but it leaves several questions unanswered. 

1. Hermit interpolation allows easy calculation of mid-point coordinates and tangent direction. However the ppt is silent about recommended approach to compute the normal at edge mid-point. Is averaging normals at end points with subsequent division by the length of result considered to be adequate? 
2. In case when edge directions are given and therefore there is no normals  information in the file (or normals should be ignored) subdivision necessitates computing surface normals for each triangle sharing the edge. I assume that it is done using data related to edges of these triangles but what are the specifics of the recommended approach? 

The answers to these two questions become important if there is a need to perform more than one level of subdivision. I could come up with my own answers to these and similar questions but it would be nice to know the official recommendations without having to pay for F2915 spec on ASTM web site.

Thanks,

- LR

[Hod Lipson]

Hi Leonid,

 

Here is the process I recommend. Jon Hiller independently implemented this in his Open Source code on the AMF wiki. Let me know if this works for you. Note that each curved triangle needs to be subdivided to depth 5, i.e. converted into 1024 planar triangles. Also note that the triangle is assumed not to be “too curved”, i.e. not more than 25% out of plane.

 

I am planning to do an update to the ASTM spec with this information, so let me know if this makes sense.

 

Here are the four key functions I used (I can send you code if you want, let me know)

 

InitializeCurvedTriangle() – This function creates a curved triangle object from information available in the mesh file. This information will include the triangle’s vertices, a normal for each vertex (if specified), and two tangents for each edge (if specified). If either a normal or a tangent is unspecified, it is set to zero. The implementation details of this function depend on the specific overall architecture of your specific AMF data structure.

EstimateTangentsNormals() – This function estimates any unspecified normal and tangent information. For example, unspecified normals can be estimated from tangents, and tangents can be estimated from normals with some additional assumptions.

This function estimates the vectors in the following order:

1.     Normalize all specified normals to length 1

2.     Normalize all specified tangents to length 1

3.     If edge tangent not specified and normal is specified, estimate edge tangent from normal as the vector perpendicular to the normal in the direction of the straight edge

4.     If edge not specified and normal is not specified, estimate edge tangent along linear edge

5.     If a normal is not specified, estimate the normal as perpendicular to two edge tangents

6.     Scale each tangent by the length of the corresponding straight-edge chord

Subdivide(depth) – This function subdivides a curved triangle into four sub curved-triangles. By recursively calling itself with an argument of depth-1 for each subdivided triangle, this function will produce 4ⁿ triangles. When depth equals zero, subdivision terminates and the vertices of the resulting triangles can be used as the final mesh.

HermiteInterpolation(p,v₀,v₁,t₀,t₁,n₀,n₁) – An interpolation function that can interpolate the point v, tangent t, and normal n, at parameter p along a Hermite curve given be two endpoints v₀ and v₁, two end tangents t₀ and t₁, and two end normals n₀ and n₁

 

--hod

 

 

Hod Lipson

Associate Prof. of Mechanical & Aerospace Engineering and Computing & Information Science

Cornell University, 242 Upson Hall, Ithaca NY 14853, USA

Office: (607) 255 1686 Lab: (607) 254 8940 Fax: (607) 255 1222

Email: Hod.L...@cornell.edu

Web: http://www.mae.cornell.edu/lipson

Administrative Assistant:  Joe Rogan jfr...@cornell.edu, (607) 255-6981, Upson 260

Calendar: http://www.mae.cornell.edu/lipson/calendar.htm

--
You received this message because you are subscribed to the Google Groups "STL 2.0" group.
To view this discussion on the web visit https://groups.google.com/d/msg/stl2/-/7OT7yRiQX8AJ.
To post to this group, send email to st...@googlegroups.com.
To unsubscribe from this group, send email to stl2+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/stl2?hl=en.

[Leonid Raiz]

Thanks Hod for your reply.

Even though I could not figure out answers to my questions from your post I did find them by reading Jon Hiller's code.  

Thanks for pointing me in that direction.

- LR

To unsubscribe from this group, send email to stl2+un...@googlegroups.com.

[Hod Lipson]

Glad you found a solution.

 

Can you explain what question you had was  not covered by the instructions below? I would like to include these instructions in version 1.2 of the format, so let me know how they can be improved.

 

--hod

 

 

Hod Lipson

Associate Prof. of Mechanical & Aerospace Engineering and Computing & Information Science

Cornell University, 242 Upson Hall, Ithaca NY 14853, USA

Office: (607) 255 1686 Lab: (607) 254 8940 Fax: (607) 255 1222

Email: Hod.L...@cornell.edu

Web: http://www.mae.cornell.edu/lipson

Administrative Assistant: Craig Ryan, cd...@cornell.edu  (607) 255-0992, Upson 258

1.      Hermit interpolation allows easy calculation of mid-point coordinates and tangent direction. However the ppt is silent about recommended approach to compute the normal at edge mid-point. Is averaging normals at end points with subsequent division by the length of result considered to be adequate? 

2.      In case when edge directions are given and therefore there is no normals  information in the file (or normals should be ignored) subdivision necessitates computing surface normals for each triangle sharing the edge. I assume that it is done using data related to edges of these triangles but what are the specifics of the recommended approach? 

The answers to these two questions become important if there is a need to perform more than one level of subdivision. I could come up with my own answers to these and similar questions but it would be nice to know the official recommendations without having to pay for F2915 spec on ASTM web site.

 

Thanks,

- LR

--
You received this message because you are subscribed to the Google Groups "STL 2.0" group.
To view this discussion on the web visit https://groups.google.com/d/msg/stl2/-/7OT7yRiQX8AJ.
To post to this group, send email to st...@googlegroups.com.
To unsubscribe from this group, send email to stl2+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/stl2?hl=en.

--
You received this message because you are subscribed to the Google Groups "STL 2.0" group.

To view this discussion on the web visit https://groups.google.com/d/msg/stl2/-/ozhxGo9MDLMJ.

To post to this group, send email to st...@googlegroups.com.

To unsubscribe from this group, send email to stl2+uns...@googlegroups.com.

[Leonid Raiz]

[Alexander Zimmermann]

Hello.
Maybe Leonid run into the same problem as I did:

The normal vector at the subdivision point is not defined in detail.

Jon Hiller therefore uses the following code for the normal vector  ns  in
http://sourceforge.net/p/amftools/code/32/tree/trunk/src/Mesh.cpp

1270         ts = (dh1*v0 + dh2*v1 + dh3*t0 + dh4*t1).Normalized();
1271

1272         Vec3D nstemp = (n0*(1-s)+n1*s).Normalized();

1273         Vec3D sd = nstemp.Cross(ts);

1274         ns = (ts.Cross(sd)).Normalized();

But there are two problems with this (I can provide examples for both, if wanted)

1) The vector   (n0*(1-s)+n1*s)    maybe zero.
    Then  nstemp  will be zero and so do  sd  and  ns .

2) Even if  nstemp  is not zero it may be orthogonal to  ts  and therfore  sd  will be zero.
    Then of course  ns  will be zero.


I recommend to calculate  ns  as follows:

1272         Vec3D sd = n0.Cross(t0)*(1-s) + n1.Cross(t1)*s;
1273         ns = (ts.Cross(sd)).Normalized();

I think this will work in any practical case.

But no matter how the calculation will be done, there should be a clear and explicit
definition of the normal vector at the subdivision point, written down in the specification.

Just my two cents.

- AZ

On Sunday, 13. January 2013 00:18:00 UTC+1 Hod Lipson wrote:

 

Can you explain what question you had was  not covered by the instructions below? I would like to include these instructions in version 1.2 of the format, so let me know how they can be improved.

From: st...@googlegroups.com [mailto:st...@googlegroups.com] On Behalf Of Leonid Raiz

Sent: Saturday, January 12, 2013 4:43 PM
To: st...@googlegroups.com
Subject: Re: Details of subdivision algorithm.

 

Even though I could not figure out answers to my questions from your post I did find them by reading Jon Hiller's code.  

[leonid]

I interpret the intent of Jon's code as linear interpolation of normals followed by an adjustment to keep the result perpendicular to interpolated tangent.  Such an intent makes sense to me. However the implementation indeed may be improved to avoid problem (2) pointed out by Alexander.

Specifically, I think the way to do it is: 

Vec3d nstemp = (n0*(1-s) + n1*3);

ns = (nstemp - ts * nstemp.Dot(ts)).Normalized;

nstemp should not be parallel to ts in all reasonable cases and therefore proposed method should not run into problem.

Similarly, issue (1) pointed out by Alexander happens only when n0 and n1 are anti-parallel . I think this case may be ignored because triangles would be too curved.

- LR

--
You received this message because you are subscribed to the Google Groups "STL 2.0" group.

To view this discussion on the web visit https://groups.google.com/d/msg/stl2/-/SE5DmExLc7UJ.

[leonid]

Correction -

Vec3d nstemp = (n0*(1-s) + n1*s);

ns = (nstemp - ts * nstemp.Dot(ts)).Normalized;

[leonid]

On second thought Jon's original formula should produce the same results as mine.  Furthermore I don't follow Alexander's argument regarding problem (2). Cross product of perpendicular non-zero vectors is never zero.

[Alexander Zimmermann]

You are rigth, issue (1) will only happen, when n0 and n1 are anti-parallel,
but the triangle will not be too curved, marginal.

Consider the example

v0 = [0,0,0]
t0 = [0,1,0]
n0 = [-1,0,0]

v1 = [1,0,0]
t1 = [0,-1,0]
n1 = [1,0,0]

The subdivision point is [.5, .25, 0]  , this is 25% away from the straight line,
which is allowed in the specification.

[Alexander Zimmermann]

Consider the following example with

a = sqrt(8)/3

v0 = [0,0,0]
t0 = [1/3,-a,0]
n0 = [a,1/3,0]

v1 = [1,0,0]
t1 = [1/3,-a,0]
n1 = [a,1/3,0]

If I'm not wrong the subdivision point is [1/2, 0, 0] and the  tangent is [4/3, sqrt(2)/3, 0].
This is parallel to the normal  [a,1/3,0] therefore the cross product will be zero.

Of course this is a non rational example and floating point numbers may not work exactly like this,
but a little inaccuracy will flip the normal vector by  180°  .

[Reinoud Zandijk]

On Sat, Jan 19, 2013 at 02:18:41AM -0800, Alexander Zimmermann wrote:
> Consider the following example with
>
> a = sqrt(8)/3
>
> v0 = [0,0,0]
> t0 = [1/3,-a,0]
> n0 = [a,1/3,0]
>
> v1 = [1,0,0]
> t1 = [1/3,-a,0]
> n1 = [a,1/3,0]
>
> If I'm not wrong the subdivision point is [1/2, 0, 0] and the tangent is
> [4/3, sqrt(2)/3, 0].
> This is parallel to the normal [a,1/3,0] therefore the cross product will
> be zero.
>
> Of course this is a non rational example and floating point numbers may not
> work exactly like this,
> but a little inaccuracy will flip the normal vector by 180� .

I've encounted this in the wild a couple of times when rendering. I resolved
it by checking the angle between the normal calculated from the tangents and
the normal calculated from the flat vectors. If it is bigger or equal than 90,
i flip the calculated normal.

The simple averaging of the two end normals to get the interpolated normal
might not be mathematically sound. If you look at the resulting triangle mesh
in the AMF demo code with the curved triangle example, you can see that the
surface of the triangle that is spun up by the two tangents has rimples if one
slices it. My own program also shows the same problem and when i use better
lighting, one can see the various triangles its created from. Since the
tangent calculations are mathematical sound and the curves meet correctly i
suspect the normal calculations since those are used exclusively for the inner
triangle when subdividing since the middle triangle has no known tangents.

What would be a better approximation of the normal in the middle of an edge?
Or is it mathematically sound after all and are we looking at FP resolution
issues?

With regards,
Reinoud

[leonid]

Both Alexander's examples have to do with triangle edges that have inflections. Perhaps the spec should be explicit in guarding against such data. Triangle edges with inflections may indeed lead to subdivided surfaces with ruffles.

> --
> You received this message because you are subscribed to the Google Groups "STL 2.0" group.

[Hod Lipson]

In the specification we will require that curved triangles do not have inflections and are not too curved (i.e. not more than 25% out of plane).
It will be the responsibility of the producer to make sure curved triangles meet this requirement, and if not subdivide them.

Can anyone check the behavior of the alternative method, wher ethe normal in the middle is perpendicular to the tangent in the middle (calculated from the Hermite interpolation function) and the direction to the other vertex? I am wondering if that solves the problem.

--hod

-----Original Message-----
From: st...@googlegroups.com [mailto:st...@googlegroups.com] On Behalf Of leonid
Sent: Saturday, January 19, 2013 9:15 AM
To: st...@googlegroups.com
Subject: Re: Details of subdivision algorithm.

[Alexander]

My counter example (1) has no inflection and is exactly 25% out of plane.

But it's correct, that example (2) only work's with inflections.

Regardless of which formula will be used, please remember to write it down in detail in the specification.

- AZ

Hod Lipson <hod.l...@cornell.edu> schrieb:

--
Diese Nachricht wurde von meinem Android-Mobiltelefon mit K-9 Mail gesendet.

[Alexander Zimmermann]

Did you mean the the third vertex of the triangle (not on the current edge)? I think this will help.

If you want to use this vertex directly, you have to introduce an additional parameter to the HermiteInterpolation method.
My suggested formula
  ns = ts  x  ( (1-s) * (n0 x t0) + s * (n1 x t1) )
avoids this additional parameter and tries to guess the direction to the other vertex with the help of the cross product of tangent and normal at the start and end vertices.

BTW: I'm not sure if inflections may be usefull in some cases, so forbiding it generally might not be what we want.

[leonid]

Using the third vertex may indeed avoid the problems pointed out by Alexander. Unfortunately such an approach would either (a) loose continuity of normals between adjacent triangles or (b) further complicate coding because of the necessity to take adjacency into account. The method based on Hermite interpolation is C1 and works well on sufficiently tessellated surfaces but as Alexander demonstrated the limits of its applicability are not well defined.  Are there known examples of misbehavior for triangles were all normals at vertices have positive dot product with flat triangle normal? May be this condition should be explicitly stated. 

I also not sure if demanding a specific subdivision algorithm in a standard is a good idea in the first place.  May be it is better to ask file produces to output vertex normals but leave it up to file consumers to compete in a way how this information is used to improve overall surface quality.  I suspect it is not realistic to expect CAD programs to produce files in the expectation of specific subdivision algorithms.  Maybe the spec just needs to prescribe that normals data be used to increase the accuracy of surface representation but does not demand a specific formula. At the end of the day flat triangles have O(h^2) accuracy of surface representation while any method taking normals into account is bound to have O(h^4) accuracy and thus produce much better results.  May be the spec should mention Hermitte interpolation as a suggestion only with additional wording describing limits of its applicability.

--
You received this message because you are subscribed to the Google Groups "STL 2.0" group.

To view this discussion on the web visit https://groups.google.com/d/msg/stl2/-/aUAF9zb-cu0J.

[Jonathan Hiller]

Hey folks,

There are two separate issues here:
1) Ambiguous normals. We've identified 2 cases (neither of which
"should" happen according to our soft guidelines - but they will in
the wild of course!) that should be specified to avoid ambiguity.
These could be addressed by the following language:
-The normal at the subdivided center point shall be linearly
interpolated from the two endpoint normals then forced to the nearest
perpendicular direction from the calculated center point edge tangent.
(formula here)
-In the case of an ambiguous normal (describe scenarios here), normals
shall be determined at 1/3 and 2/3 of the edge according to (above).
The centerpoint normal shall then be recalculated according to (above)
using the 1/3 and 2/3 normals instead of the endpoint normals.
-Calculated normals shall always be in the direction facing outward
from the triangle (we can explore ways to formalize this statement,
although it's borderline obvious to include in the spec)


2) The second underlying (bigger) issue is the "ruffling". This is a
known limitation of the simple recursive triangle patch scheme. For
the computer graphics people, we're dealing with c1 continuity
(http://en.wikipedia.org/wiki/Smooth_function#Order_of_continuity).
The 0th and 1st derivatives match (position and angle) but curvature
(2nd derivative) continuity is not enforced. Subdividing triangles
with constant curvature (i.e. on a sphere or cylinder) yields
basically perfect results. Triangular patches with slightly varying
curvature are more than adequate. However triangles with highly
variable 2nd derivatives (for example those with inflection points) do
not fare so well -- hence the guideline on limiting on out-of-plane
curvature.

A possible solution here is to generate all 1024 sub-triangles
directly from the original specified triangle (i.e. 2d Hermite
interpolation or similar). While certainly not as elegant as
recursion, I think this is the only way to always avoid ruffled
surfaces without changing the spec to require significantly more
information.

~Jon

ps. Hod, regarding the alternative method you mentioned: If I
understand correctly this would drop us down to only C0 continuity
since adjacent triangles sharing an edge and sharing endpoint normals
would end up with different midpoint normals - hence the surface would
not be smooth. Do I misunderstand your descriptiuon?

On Sat, Jan 19, 2013 at 2:31 PM, Hod Lipson <hod.l...@cornell.edu> wrote:
>
> In the specification we will require that curved triangles do not have inflections and are not too curved (i.e. not more than 25% out of plane).
> It will be the responsibility of the producer to make sure curved triangles meet this requirement, and if not subdivide them.
>

> Can anyone check the behavior of the alternative method, wher ethe normal in the middle is perpendicular to the tangent in the middle (calculated from the Hermite interpolation function) and the direction to the other vertex? I am wondering if that solves the problem.
>
> --hod
>
>
>

> -----Original Message-----
> From: st...@googlegroups.com [mailto:st...@googlegroups.com] On Behalf Of leonid
> Sent: Saturday, January 19, 2013 9:15 AM

> To: st...@googlegroups.com
> Subject: Re: Details of subdivision algorithm.
>

[leonid]

Here is something else to consider. Jon's two step approach in normal interpolation gave me an idea. What if instead of linear interpolation we would consider rotational interpolation of normal vector from start to and end direction where angle of rotation varies linearly depending on parameter value? Have such an approach been considered? It still creates C1 surfaces and avoids the problem of normal ambiguity. I wonder how it would affect ruffles?

[Martin Wicke]

We had the discussion about the subdivision scheme before. We noted then
(https://groups.google.com/d/msg/stl2/-/Mgdq3OskxiAJ) that the bumpiness
is inherent to the subdivision scheme, and cannot be fixed at all,
unless a subdivision scheme with C2 continuity (almost) everywhere is
used. Modified Butterfly
(http://mrl.nyu.edu/~dzorin/papers/zorin1996ism.pdf) seems to be the
obvious choice here, since it is interpolating, C2,
topology-independent, and easy to implement. It does, however, break the
independence of triangles -- you need information about neighboring
triangles to compute the subdivision. Btw., this is true of all C2
subdivision schemes, there is no free lunch.

The consensus back then was that the independence of triangles was of
higher value than smooth surfaces, and the artifacts are acceptable in
the anticipated use cases (smooth triangles are not going to be very
curved anyway, and the producer is responsible of making sure the
surface subdivided according to AMF rules is close enough to the
original). I wonder whether that has changed. I think something like
Jon's solution will fix the immediate problem, but I am not a huge fan
of the amount of hacks that are required to make this work (for purely
aesthetic reasons, also I don't know how many more we will need).

However, I strongly suggest that the standard prescribe the subdivision
scheme, whichever it is, exactly (including all the little tricks to
make it work). If necessary (and honestly, I think it is necessary), I
would include a full implementation in the standard, like you did for
the rand() function.

If the standard doesn't exactly nail down the subdivision, the surface
geometry is essentially not defined. I've attached Fig. 1 from Zorin et
al.'s paper linked above as a demonstration of what even minor changes
to the subdivision scheme can do to the surface.

Martin

[Alexander Zimmermann]

I've thought in that direction, too. But the resulting formulas require trigonometric functions and look really very ugly.
And there is still one case, when the normal rotates exactly 180° , where this is not unique.

- AZ