## Variants on ROAM Priority Computations

The method I suggest anyone try first is to do the full view projection
of the wedgie thickness segment for the center of a triangle.  The
reason for this is not to speed up the computation.  I suggest it
because it is fast to write the code, with little room for mistakes, and
it is good for doing sanity checks on faster or more accurate methods
you might try later, such as the method you are using, or the ones I've
outlined.

[Seumas]
> TriZDist = the camera space Z distance of the center of the triangle's
> hypotenuse (I pick this rather than the tri center since the
> hypotenuse center must be computed to find the coordinates of the
> children, if passing triangle vertices on the stack in a recursive
> function).

Any point will work, and your choice is fast.  I choose the center-point
because it gets distorted the least in the worst case compared to the
rest of the representatives you might choose.  Again, I was going for
reliability in that suggestion for the first thing to try, not speed.

The method I suggest for fast direction-sensitive priority computations
is something like the first-order approximation that I posted a few days
ago.  This gives a tighter measurement of the actual errors than using
thickness>Z*chicken_sacrifice, but takes more number crunching (the
expensive stuff like square roots is replaced with fast look-up tables,
and you only end up doing part of a projection, but is is obviously more
expensive than a multiply and compare).  The chicken_sacrifice method
trades away quality to gain speed.  In other words, you might be using
20-50% more triangles to get the same accuracy.  This is fine if you
have the triangle budget to burn, of course.  But you also take more
time because you are dealing with more triangles for the same quality.
If your bottleneck is the priority computation, then this is a win, but
with frame-coherence acceleration this is not the bottleneck, so the
higher-quality priority is worth it.

Before I begin answering your questions requesting ROAM paper clarification,
let me suggest a simpler metric that is fairly robust and can get you going
quickly.  Take the mid-point of your triangle, and create two points, one at
the top of the wedgie and the other at the bottom.  If your centerpoint is
(x,y,z), and the wedgie thickness is e_T, then the two points are
(x,y,z+e_T) and (x,y,z-e_T).  Send these two points through your view
transform, perspective and all, into pixel space.  Let's call these two
points (sx0,sy0) and (sx1,sy1).  Next compute the distance between them,
se=sqrt((sx1-sx0)^2+(sy1-sy0)^2).  Call se your priority and you will
*almost* have a working system.  In a recent post I described how to deal
with the occasional child that has a higher priority than its parent--you
will need that fix for this simple bound estimate.  The main reason for
using this is that it is very close to what you would get with the more
robust bound, but is really fast to code and can be used to debug (do a
sanity check) on the more complex things you might try later.

I actually switch between two different priority computations,
one slow and robust and the other fast but sometimes wrong.

A priority ideally is a tight upper bound on how far any point
on the triangle is from where it should be (with respect to the
finest-resolution surface), measured in pixels.  This priority
can be tweeked of course to add detail near an object sitting on
a terrain, for example, so it sits right.  Note that the pixel
distance errors are zero if you are looking straight down at
a terrain triangle in the middle of the screen, regardless of how
bad the height error is.  The error goes up as the triangle goes
towards the edges of the screen, or if it is rotated so you aren't
looking straight down.  The error becomes largest when you
look at the triangle edge-on, i.e. when the triangle is right
on a silhouette edge.  This gives higher priority exactly where
it counts visually.

The slow priority computation takes your precomputed maximum
height error for a triangle and produces guaranteed upper bounds
on the screen error in pixels.  This is described in fair detail in the
ROAM paper (http://www.llnl.gov/graphics/ROAM/).  In short,
you take the maximum numerator divided by the minumum
denominator for the three triangle vertice's perspective-transform
error.  This works in all cases except when the denominator goes
to zero (this is avoided by setting priority to max when the triangle
ventures inside the near clipping plane).

David Smith wrote:

> [snip]  I want
> to understand the principle before I understand the
> speed up.
>
> "Let (p,q,r) be the camera-space coordinates of a point
> w(v) without perspective."
>
> I think I understand this, basically take a point
> an pass it through my camera position and orientation
> transform.

Yes, translate the camera position to the origin, then rotate so that p is
the screen-right direction, q is the screen-down direction, and r is the
screen-into direction (for a right-hand coordinate system).

>
>
> You then talk about the perspective transform letting
> s = (p/r, q/r)  but you never use s. Curious, why?

I define s(v) two paragraphs before the sentence you quote.  I mention it
in this sentence so you know where the formula hat(dist)(v)=...etc... is
coming from.  You don't actually need s in the code you write.

>
>
> Now,
>
> "Let (a,b,c) be the camera-space vector corresponding
> to world-space thickness vector (0,0,et)....."
>
> Does this last paragraph of section 6.2 mean substitute
> the three vertices of the triangle with the additional
> error, et, tacked onto their height component converted
> into screen space and substituted into equation (3)?
>
> I am getting a couple of negative distances and that can't
> be right.  So I still must be missing something.
>
> -DaveS

The sentence before the one you mention is the key to understanding the
formula.  I indicate the screen-space distortion at a corner vertex is
bounded by projecting the thickness segment--this is the simple method I
suggested at the beginning of this message, applied now at the triangle
corners.  What we need to do to get the robust bound is split up this
projection into two steps, (1) into camera space (without perspective), then
(2) bound in screen space.  Let's say your world-space (x,y,z) goes to
camera-space (x',y',z').  Take the two endponts of a thickness segment,
(x0,y0,z0) and (x1,y1,z1), and these transformed to (x0',y0',z0') and
(x1',y1',z1').  What we need in the formula is the vector between these,
(a,b,c)=(x1'-x0',y1'-y0',z1'-z0').  Note that (a,b,c) will be the same
regardless of where you move the thickness segment, as long as its direction
and length remain constant.  In effect you only need to perform the rotation
part of your camera transform (and some scaling), ignoring the translation,
on the vector (0,0,e_T).  So the bottom line is that (a,b,c) depends only on
your rotation matrix (and perhaps scaling) and the wedgie thickness, and is
the same for all three triangle vertices.

So now for three vertices you have three (p,q,r)'s and one (a,b,c).  In your
code you compute 2/(r^2-c^2) three times (once per vertex), and take the max
of the three.  You also compute ((ar-cp)^2+(br-cq)^2) three times and take
the max.  You then take the first max times the square root of the second
max as your distortion bound and priority.

The one major caveat here is that you need to make sure r^2>c^2, or you
could get divide-by-zero or negative "bounds".  If
r^<=c^2+tiny_positive_number, then set the priority to an artificially large
value.  Also make sure your near clipping plane is as far forward as your
application allows (this makes good sense for other reasons as well, like
good use of finite depth-buffer precision).  This way the "artificially
prioritized" triangles quickly get culled outside the view frustum.  You
never need to compute the priority of a triangle labeled OUT, and you and up
without any artificially-prioritized triangles after a few splits.  NOTE:
what I describe here is better than what I put in the paper--the method
there causes terrain intersecting the near clipping plane to be subdivided
to the finest resolution at the expense of all the other visible terrain
(which I consider to be a bug :).

The fast method, not described in the paper, is a good approximation
to the slow method when the triangle is small and far from the
eyepoint.  It is a first-order approximation taken at the center of
the triangle.  You take the derivative of the projected pixel error as
the length of your height error increases from zero.  You then
scale this derivative by the actual height error to get the approximate
projected error.  Here is a totally cryptic code fragment implementing
this:

x0=q->cp; y0=q->cp; z0=q->cp;
x1=m*x0+m*y0+m*z0+m;
y1=m*x0+m*y0+m*z0+m;
z1=m*x0+m*y0+m*z0+m;
w1=m*x0+m*y0+m*z0+m;
m32=m/(w1*w1);
dx=m/w1-m32*x1;
dy=m/w1-m32*y1;   // FIXED 06-13-01: Thanks Dimitris P.!
dz=m/w1-m32*z1;   // ditto
e2=q->e;
r=e2*e2*(dx*dx+dy*dy+dz*dz)*100.0;
ri=(int)floor(r*(double)QCELL_RITAB);
if (ri<0) ri=0; if (ri>=QCELL_RITAB) ri=QCELL_RITAB-1;
iq=qcell_ritab[ri];

This is maybe ten times as fast as the slow method, so it is worth
decyphering.  The variables are:

q -- the bintree triangle data structure
cp[0..2] aka (x0,y0,z0) -- the triangle center point
m[0..3][0..3] -- the 4x4 view transform matrix
(dx,dy,dz) -- the derivative of the projected errors
q->e -- the "wedgie thickness" precomputed height error
r -- the projected error squared, scaled up a bit
ri -- r converted to fixed point int for table lookup
iq -- finally the priority, obtained from a table lookup
that does the square root and perhaps a log() and scaling.
This is an int in a range from 0..QIMAX (e.g. QIMAX=2047).

Now you've got an iq priority index, and are ready to queue
your triangle.

Hi Kallen,

Glad to hear you've made such good progress on your ROAM
implementation.

You wrote:

> I have my ROAM implementation almost complete but I'm stuck on the priority
> computation part. My implementation is straight from the paper except I
> support paging of bin_trees.
>
> Your paper has a formula that goes dist(v)=| |s(v)-sT(v)| |. Now to my
> understanding this a typical LOD distortion calculation.  Do I use this
> formula for distortion? I assume not as it does not use the wedgie
> thickness' I spent so long computing.  If its not distortion then what is it
> used for or is it used at all?

The formula you mention is the exact error measure per surface point v
(the distance in pixels between the position of the point on the finest-
resolution terrain and the position of the point on the approximate
terrain triangulation T).  This is not used directly in an implementation
(since you would have to use it an infinite number of times to get
a guaranteed bound...).  Instead you typically compute a conservative
upper bound on this error measure.  The wedgie bounds you precomputed
help with this.  Ideally you would find the thickest part of the projected
wedgie and use that.  Unfortunately the thickest part of the wedgie can
be anywhere on the wedgie boundary.  So we end up with the formula
with the p,q,r's and a,b,c's in it to get a conservative upper bound on the
distortion.  Why a conservative upper bound you ask?  This is the only
kind of guarantee you can get that you don't get significant and
unacceptable errors if you are unlucky.  If you are willing to relax this
requirement a bit, you can get much faster priority computations where
you sometimes get unlucky or just don't know for sure how bad the
distortions are.  In fact in my implementation I typically turn on
a kind of first-order approximation for all triangles that are fairly
small and far away from the camera point (which is the vast majority
of triangles you are processing).  My code is rather cryptic and takes
some decoding and deduction, but here is what my computation
looks like:

int qcell_priority(qcell q)
{
int iq,k,ri;
qter qt;
qcell qp;
qlos ql;
qlosref qlr,qlrx;
vmatmat m;
double r,x0,y0,z0,x1,y1,z1,w1,dx,dy,dz,*cp,m32,e2;

qt=q->qt;
if (q->l>=qt->tm_l) return 0;

// update line-of-site list
while (q->losref_list) qlosref_destroy(q->losref_list);
if (qp=q->qp) {
for (qlr=qp->losref_list;qlr;qlr=qlr->qlr1)
{ ql=qlr->ql; if (qlos_overlap(ql,q)) qlrx=qlosref_create(q,ql); }
}else{
for (ql=qt->los_list;ql;ql=ql->ql1)
{ if (qlos_overlap(ql,q)) qlrx=qlosref_create(q,ql); }
}
if (q->losref_list) return QTER_QUEUEN-1;

if (q->l<3) return qcell_priority_good(q);
if (q->flags&QCELL_IN0) {
m=qt->vm->m; cp=q->cp;
x0=cp; y0=cp; z0=cp;
x1=m*x0+m*y0+m*z0+m;
y1=m*x0+m*y0+m*z0+m;
z1=m*x0+m*y0+m*z0+m;
w1=m*x0+m*y0+m*z0+m;
m32=m/(w1*w1);
dx=m/w1-m32*x1;
dy=m/w1-m32*y1;    // FIXED 06-13-01: Thanks Dimitris P.!
dz=m/w1-m32*z1;    // ditto
e2=q->e;    // fixed 06-01-13
r=e2*e2*(dx*dx+dy*dy+dz*dz)*100.0;
ri=(int)floor(r*(double)QCELL_RITAB);
if (ri<0) ri=0; if (ri>=QCELL_RITAB) ri=QCELL_RITAB-1;
iq=qcell_ritab[ri];
}else iq=QTER_QUEUEN-1;
return iq;
}

I'll give a bit of guidance as to what some of this is.  The data structure for

a bintree triangle is a qcell.  This is attached to a terrain data structure
qter.  The priority computation proceeds by checking for line-of-sight
corrections (ignore this), falling back on the fully conservative, expensive
and slow priority computation (qcell_priority_good()) if the triangle is large,

setting the priority to max if the wedgie is not completely behind the near
clipping plane, and otherwise computing the first-order approximation
I mentioned.

The 4x4 matrix m is the view transform.  I perform a differential transform
(compute
the derivative of the projection for the center point of the triangle), giving
vector
(dx,dy,dz).  I then scale this vector by the wedgie thickness and compute the
magnitude squared.  I use a table look-up to convert this to a priority index.
QCELL_RITAB is 1meg, and QTER_QUEUEN is 2048.

>
>
> Then there is the second large formula that uses (p,q,r) and (a,b,c). But I
> thought this formula was for computing an upper bound.  Is this the formula
> for computing distortion?  If not then why must I have an upper bound.
>
> I'm confused can you help me clarify which of these formulas is the
> distortion.
>
> Lastly, in implementing the deferred priority computation lists your first
> sentence states
> "Given a velocity bound on the viewpoint".  What does this mean?

Really what you are shooting for here is a prediction of when a triangle will
be split or merged at the earliest.  You need to estimate or bound two
quantities
to do this.  One is the change of the crossover priority over time (the
crossover
priority is the priority that, if given as a tolerance, would produce the
output
mesh you obtain--think of running the queue optimization with only a tolerance
as the termination condition and not a triangle count).  The second quantity is

the priority of the triangle over time.  If you think of a 1-D plot of these
two
quantities over time, you will get two cones widening to the right that
eventually
intersect.  The point of intersection is the earliest time that a merge or
split
could happen.  So everything comes down to bounding crossover and triangle
priorities over time.

The crossover priority we just measured empirically and observed that less
than 1% change occurs per frame in our application.  For the triangle priority,

just consider all the possible camera motions that could possibly happen over
time.  If your camera can jump around arbitrarily, then no bound is possible.
Typical motion is constrained, however, to some finite maximum speed of
the camera position.  This is enough to allow you to compute a cone-shaped
bound on priority over time for a given wedgie (by "priority" here I really
mean
the conservative upper bound on priority or the differential estimate, whatever

you choose to use).  Think for example of fixing your reference frame to be
that of the camera, and think of the motion as occurring entirely by the
wedgie.
Then the wedgie will move by at most the maximum velocity per second, so
and point in the wedgie will lie in a sphere whose radius grows linearly over
time.  You can work out the math on what this means for the wedgie bound
or differential projection (it is a pretty simple relationship).

A slight modification you might consider to the reasoning above is to allow
very high rates of change to camera orientation, but still have a finite bound
on the velocity of the camera position.  In this case you consider a different
projected error: project the errors onto a sphere.  Now the errors are
rotationally
invariant.  The reasoning above continues to work if you now always imagine
your camera as pointing straight towards the wedgie.

## Contact Information

Mark A. Duchaineau -- duchaineau1@llnl.gov

Updated June 13, 2001