\input header
\input epsf

\title{Subdivision Surfaces analysis and interrogation exercises: Discussion}

\section{Analysis}

\pop Over what shape region is the basis function non-zero, of the ternary scheme whose non-zero coefficients form the pattern
$$\matrix{&&&*&*\cr
&&*&*&*&*\cr
&*&*&*&*&*&*\cr
*&*&*&*&*&*&*&*\cr
*&*&*&*&*&*&*&*\cr
&*&*&*&*&*&*\cr
&&*&*&*&*\cr
&&&*&*}$$

\pop The scheme is over the square grid, and must be `neither', because the rows have 
an even number of entries.

The convex hull of the non-zero entries is octagonal, with sides alternately 1 unit and $3\sqrt2$ units long.

A ternary scheme cannot have any influence beyond a boundary 3/2 times as large as the convex hull of its entries.

Applying a copy of that boundary scaled down (to 1/2 times the size of that convex hull) to each of the entries at the top of the mask indicates that the influence cannot fill that theoretical boundary.  The true boundary has a fractal form, consisting of pieces of line at 45 degrees, touching the scaled convex hull all the way along the sloping sides and in a cantor set along the horizontal and vertical sides.


\filbreak
\pop What is the precision set (constant, linear , quadratic ?) of the binary scheme
whose mask is
$$\matrix{
&&1&&3&&1\cr\ \cr
&3&&4&&4&&3\cr\ \cr
1&&4&&8&&4&&1\cr\ \cr
&3&&4&&4&&3\cr\ \cr
&&1&&3&&1} \qquad/14$$

\pop The scheme is on a triangular grid and is binary.  By symmetry it is primal, and so it has new points at the old ones and also at the midpoints of the edges, which come in three orientations.  The stencils for the vertex vertices and the horizontal edges are
$$\matrix{
&1&&1\cr\ \cr
1&&8&&1\cr\ \cr
&1&&1}\quad/14\qquad \qquad\hbox{ and }
\matrix{
&3\cr\ \cr
4&&4\cr\ \cr
&3}$$

Summing these independently we find that 1+1+1+1+1+1+8 = 14 and 3+4+4+3=14, so that the scheme has constant precision.  We only need to check one of the edge stencils by symmetry.

Multiplying by the linear field of old vertex values
$$\matrix{
-2&&0&&2\cr\ \cr
&-1&&1\cr\ \cr
-2&&0&&2\cr\ \cr
&-1&&1\cr\ \cr
-2&&0&&2\cr\ \cr
&-1&&1}$$
gives the new field
$$\matrix{
-14&&0&&14\cr\ \cr
&-8&&8\cr\ \cr
-14&&0&&14\cr\ \cr
&-8&&8\cr\ \cr
-14&&0&&14}\qquad/14$$
which is not a linear field, and so the scheme has only constant precision.



\filbreak

\pop What is an appropriate univariate scheme to use for the boundaries of a surface defined by Loop's scheme, a binary scheme with mask
$$\matrix{
&&1&&2&&1\cr\ \cr
&2&&6&&6&&2\cr\ \cr
1&&6&&10&&6&&1\cr\ \cr
&2&&6&&6&&2\cr\ \cr
&&1&&2&&1}\qquad/16$$

\pop Multiplying the mask by the extruded cardinal data set
$$\matrix{
&&0&&0&&0\cr\ \cr
&0&&0&&0\cr\ \cr
1&&1&&1&&1\cr\ \cr
&0&&0&&0\cr\ \cr
&&0&&0&&0}$$
we get 
$$\matrix{
&&2&&2&&2\cr\ \cr
&8&&8&&8\cr\ \cr
12&&12&&12&&12\cr\ \cr
&8&&8&&8\cr\ \cr
&&2&&2&&2}\qquad/16$$
and so a suitable boundary scheme is $[2,8,12,8,2]/16$ , or $[1,4,6,4,1]/8$, the cubic spline scheme.  This is interesting, because the Loop scheme is a box-spline with quartic surface pieces, and yet it naturally gets cubic boundaries.

\filbreak
 \pop Analyse the continuity of Loop's scheme at an ordinary vertex

 \pop Spotting that the scheme is a box-spline says that it has
 quartic pieces meeting with $C^2$ continuity.  However, it is worth
 also running through the eigenanalysis.

Let us set up an algebraic data-set for one coordinate, and see how it
is refined.

We observe that every new vertex in or on the first ring of the data
is defined by data in or on the second ring, and, indeed that no
points outside the second ring influence their images after one
refinement.  We only therefore need to consider four types of point:

$$\matrix{
&&d\cr\ \cr
&b&&c\cr\ \cr
a&&b&&d}$$
where $a$ is the centre of the refinement.

\filbreak
Label the original vertex values with upper case letters, $A$ denoting
the central value and $B_i$ denoting the value on the $i^{th}$ spoke.
Label the `C' points with $C_{i+1/2}$, so that the zero subscript is
always on the zeroth spoke.

Label the values after refinement with corresponding lower case letters.

Then the new values are related to the old by
$$\eqalign{
16a &= 10A + \sum_iB_i\cr
16b_i &= 6A + 6B_i + 2B_{i+1} + 2B_{i-1}\cr
16c_{i+1/2} &= 2A +6B_{i+1}+6B_i + 2C_{i+1/2}\cr
16d_i &= A + B_{i+1}+B_{i-1}+10B_i + C_{i+1/2}+D_{i-1/2} + D_i}$$

\filbreak
If we apply the Fourier\footnote{$^3$} {This is not a Fourier
transform but the closely related Discrete Cosine Transform.  The full
Fourier uses $e^{ij\omega}$ in place of the cosines here. DCT is a
useful short cut for manual calculations where we can rely on symmetry
to give the second eigencomponents which are associated with negative
$\omega$ values.}
 transform that 
${B}_i = \sum_{\omega}\hat B_{\omega}\cos(i\omega)$ and 
$\hat b_{\omega} = \sum_i {b}_i\cos(i\omega)$ 
these equations become

$$\eqalign{
16a &= 10A + 6{\hat B}_0\cr
16{\hat b}_{\omega} &= 6A + {\hat B}_{\omega}(6+4\cos(\omega)) \cr
16{\hat c}_{\omega} &= 2A +12{\hat B}_{\omega}\cos(\omega/2) + 2{\hat C}_{\omega}\cr
16{\hat d}_{\omega} &= A + {\hat B}_{\omega}(10+2\cos(\omega)) 
+ 2{\hat C}_{\omega}\cos{\omega/2} +{\hat D}_{\omega}}$$

\filbreak
For $\omega=0$, $\cos{\omega}=\cos{\omega/2}=1$, and these become
$$\eqalign{
16a &= 10A + 6{\hat B}_0\cr
16{\hat b}_{0} &= 6A + 10{\hat B}_{0} \cr
16{\hat c}_{0} &= 2A +12{\hat B}_{0} + 2{\hat C}_{0}\cr
16{\hat d}_{0} &= A + 12{\hat B}_{0} + 2{\hat C}_{0} + {\hat D}_0 }$$
which may be expressed as the matrix equation
$$16\left[\matrix{a\cr{\hat b}\cr{\hat c}\cr{\hat d}}\right]
= \left[\matrix{10&6&0&0\cr
6&10&0&0\cr
2&12&2&0\cr
1&12&2&1}\right]
\left[\matrix{A\cr{\hat B}\cr{\hat C}\cr{\hat D}}\right]$$

The matrix here has eigenvalues 1, 1/4, 1/8 and 1/16

\filbreak 
For $\omega = \pi/3$ , $\cos{\omega}=1/2,\ \cos{\omega/2}=\sqrt3/2$.
We also lose all the $A$ and $a$ terms, so that the equations become

$$\eqalign{
16{\hat b}_{\omega} &= 8{\hat B}_{\omega} \cr
16{\hat c}_{\omega} &= 6\sqrt3{\hat B}_{\omega} + 2{\hat C}_{\omega}\cr
16{\hat d}_{\omega} &= 11{\hat B}_{\omega}  
+ \sqrt3{\hat C}_{\omega}+{\hat D}_{\omega}
}$$

which may be expressed as the matrix equation
$$16\left[\matrix{b\cr c\cr d}\right] =
\left[\matrix{
8\cr
6\sqrt3&2&\cr
11&\sqrt3&1
}\right]
\left[\matrix{B\cr C\cr D}\right]$$
whose eigenvalues are 1/2, 1/8 and 1/16.

\filbreak
For $\omega = 2\pi/3$ , $\cos{\omega}=-1/2,\ \cos{\omega/2}=1/2$.
We again lose all the $A$ and $a$ terms, so that the equations become

$$\eqalign{
16{\hat b}_{\omega} &=  4{\hat B}_{\omega} \cr
16{\hat c}_{\omega} &= 6{\hat B}_{\omega}\cos(\omega/2) + 2{\hat C}_{\omega}\cr
16{\hat d}_{\omega} &= 9{\hat B}_{\omega} 
+ {\hat C}_{\omega}\cos{\omega/2} +{\hat D}_{\omega}}$$
so that the eigenvalues are 1/4, 1/8 and 1/16.

\filbreak
Finally, for $\omega=\pi$, $\cos(\omega) = -1,\ \cos(\omega/2)= 0$

$$\eqalign{
16{\hat b}_{\omega} &= 2{\hat B}_{\omega} \cr
16{\hat c}_{\omega} &= 2{\hat C}_{\omega}\cr
16{\hat d}_{\omega} &= 8{\hat B}_{\omega}  +{\hat D}_{\omega}}$$

has eigenvalues 1/8, 1/8 and 1/16.

\filbreak
Thus the spectrum has: an eigenvalue of 1 from $\omega=0$, two
eigenvalues of 1/2 from $\omega =\pm\pi/3$, two eigenvalues of 1/4
from $\omega=\pm2\pi/3$, and one eigenvalue of 1/4 from $\omega =0$.
All other eigenvalues are smaller, and in fact equal to 1/8 or less.

Thus we have $\mu=\lambda^2$ as needed for bounded curvature.

In fact the eigenvectors associated with $\lambda$ give as the natural
configuration the regular grid, and those associated with $\mu$ give
second order natural configurations which are quadratic functions.

\filbreak
 \pop Analyse the continuity of Loop's scheme at an extraordinary
 vertex of valency 8

 \pop This is more interesting.

The matrix for $\omega=0$ is
$$\left[\matrix{
10&6\cr
6&10\cr
2&12&2\cr
1&12&2&1}\right]/16$$

and for general $\omega>0$ is

$$\left[\matrix{
6+4\cos(\omega)\cr
12\cos(\omega/2)&2\cr
10+2\cos(/omega)&2\cos(\omega/2)&1}\right]/16$$

\filbreak
We need to check the eigencomponents for $\omega=0, \pi/4, \pi/2, 3\pi/4, \pi$.
The first and last of these give one component each, the middle three
pairs of equal eigenvalues.

\filbreak
$$\matrix{
\omega&\cos(\omega)&\cos(\omega/2)&\hbox{eigenvalues}\cr
\ \cr
0&1&1&1,1/4,1/8,1/16\cr
\pi/4&.7071&.9237&.55,1/8,1/16\cr
\pi/2&0&.7071&6/16,1/8,1/16\cr
3\pi/4&-.7071&.3823&3.17/16,1/8,1/16\cr
\pi&-1&0&1/8,1/16}$$

$$\eqalign{\lambda &= 0.55\cr
\mu_0 &= .25 < .55^2 (=.304)\cr
\mu_2 &= .375 > .55^2}$$

\filbreak
All other eigenvalues are less than $\lambda^2$.  Thus for nicely
symmetric data around the extraordinary point there will be a flat
point where the curvature tends to zero, while for saddle-shaped data
there will be a sharp point where the curvature is unbounded.

The next problem is to find an appropriate set of rules so that the
eigenvalues behave, giving $\mu_0 = \mu_2 = \lambda^2$ exactly.  To
achieve this the coefficients need to depend on the valencies of the
vertices to which they refer.




\filbreak
 \pop Does Loop's scheme have lateral artifacts ?

 \pop Look back at the boundary scheme question.  The extruded data gets
refined into a new data set still extruded, and so there are no
lateral artifacts.

It is also possible to divide the scheme by $(1+x)$ to give the quotient
$$\matrix{
&&1&&1\cr\ \cr
&2&&4&&2\cr\ \cr
1&&5&&5&&1\cr\ \cr
&2&&4&&2\cr\ \cr
&&1&&1}$$
without any remainder, and so the scheme can be shown this way also to
be artifact-free.



\section{Interrogations}

 \pop Write an algorithm (in the style of the pseudocode in the notes
 and to the same level of detail) for the intersection of two surfaces
 by subdivision.  Note that it will come out to rather more lines, as
 there is more to be covered.

 \pop Here it is

\cleartabs
\+\cr
\+Pointset:Srfint(Surface:$S_A$, $S_B$)\cr
\+Surface:$S_1$,$S_2$\cr
\+Begin&\cr
\+&If\qquad&Hull($S_A$) overlaps Hull($S_B$)\cr
\+&Then&If\qquad&Simple($S_A$)\cr
\+&    &Then&If\qquad&Simple($S_B$)\cr
\+&    &    &Then&Srfint&(Approx($S_A$),Approx($S_B$)\hskip 150pt &(1)\cr
\+&    &    &Else&$S_1,S_2$ := Halves($S_B$)\cr
\+&    &    &    &Union(Srfint($S_A$,$S_1$),Srfint($S_A$,$S_2$))&&(2)\cr
\+&    &    &Endif\cr
\+&    &Else&If &Simple($S_B$)\cr
\+&    &    &Then&$S_1,S_2$ := Halves($S_A$)\cr
\+&    &    &    &Union(Srfint($S_1$,$S_B$),Srfint($S_2$,$S_B$))&&(3)\cr
\+&    &    &Else&If\qquad&Hull($S_A$) $<$ Hull($S_B$)&(4)\cr
\+&    &    &    &Then&$S_1,S_2$ := Halves($S_B$)\cr
\+&    &    &    &    &Union(Srfint($S_A$,$S_1$),Srfint($S_A$,$S_2$))\cr
\+&    &    &    &Else&$S_1,S_2$ := Halves($S_A$)\cr
\+&    &    &    &    &Union(Srfint($S_1$,$S_B$),Srfint($S_2$,$S_B$))\cr
\+&    &    &    &Endif\cr
\+&    &    &Endif\cr
\+&    &Endif\cr
\+&Else&$\emptyset$\cr
\+&Endif\cr
\+End\cr
\+\cr

I did warn you that this would be quite long.  

Notes:

(1)\qquad This is a special Srfint for intersecting facets.

(2)\qquad The Srfint here could be a sawn-of version which already
knows that its first argument is simple, so it doesn't have to test
again. The Union here is also a version that knows that the gap it has
to fix up is the one between the two second arguments.

(3)\qquad Now the Srfint could be a simpler version which knows that
its second argument is simple, and Union knows that the gap it has to
fix up is the one between the two first arguments.

(4)\qquad The $<$ denotes a test on the maximum size of the two boxes.

\filbreak
 \pop Write an algorithm for finding the first intersection of a given
 straight line (defined by a point $Q$ and a tangent vector $T$) with
 a surface $S$ by subdivision.  `First' means that the value of
 $[P-Q].T$ has the smallest positive value of all points $P$ which are
 both on the line and on the surface.  (It is the one which you see
 when your eye is at $Q$ looking along $T$.)

 \pop This is a nice example for subdivision because it has lots of tests
for culling and also the requirement to find the nearest, so it needs
a priority queue.

\cleartabs
\+\cr
\+Point: NPTL($S$,$Q$,$T$)\cr
\+Point:$R$\cr
\+Point,dist: best so far\cr
\+Begin&\cr
\+&If\qquad&Hull($S$) intersects $[Q,T]$\cr
\+&and&Hull($S$).$T$ $>$ $Q.T$\cr
\+&Then&push $S$ into queue at Hull($S$).$T$\cr
\+&&best so far = no point, very far away\cr
\+&&Until\quad&queue is empty\cr
\+&&Do&Pop $S$ from queue\cr
\+&&&If\qquad&Hull($S$).$T$ $<$ best so far\cr
\+&&&Then&If\qquad&Simple($S$)\cr
\+&&&&Then&$R$ = NPTL(Approx($S$),$Q$,$T$)\cr
\+&&&&&Update best so far with $R$\cr 
\+&&&&Else&$S_1,S_2$ = Halves($S$)\cr
\+&&&&&If\qquad &Hull($S_1$) intersects $[Q,T]$\cr
\+&&&&&and&Hull($S_1$).$T$ $>$ $Q.T$\cr
\+&&&&&Then&push  $S_1$ into queue at Hull($S_1$).$T$\cr
\+&&&&&Endif\cr  
\+&&&&&If &Hull($S_2$) intersects $[Q,T]$\cr
\+&&&&&and&Hull($S_2$).$T$ $>$ $Q.T$\cr
\+&&&&&Then&push  $S_2$ into queue at Hull($S_2$).$T$\cr
\+&&&&&Endif\cr  
\+&&&&Endif\cr  
\+&&&Endif\cr
\+&&Enddo\cr
\+&&return best so far\cr
\+&Else\cr
\+&&return no intersection\cr
\+&Endif\cr
\+End\cr

There is a variant on this which checks whether the thing being pushed
into the queue is going to be worth looking at when it comes off, in
terms of being potentially nearer than the current best so far.  This
is the type of judgement which can only be made on the basis of
experiments with realistic data.  Making the extra test makes things
slower, because it takes time.  It might make them faster by keeping
the queue slightly shorter.


\filbreak
 \pop The silhouette curve is defined as the border curve between
 parts of a surface $S$ which are lit by illumination from a given
 point $Q$ and those parts which are not.

What is the equation linking a point $P$ of the surface, the given
point $Q$ and the surface normal $N$ to $S$ at $P$~?

Write a pseudocode algorithm to compute the silhouette by subdivision.

 \pop Oops.  The question as intended was meant to ignore the
 possibility of one piece of surface shading another.

The equation is just
$$[P-Q].N = 0$$

We are therefore going to need a normal hull to judge whether a piece
of surface might contain a part of the silhouette curve.  That test
needs to compare the inner product of two hulls. This produces an
interval. For two boxes this can be computed by just checking the corners.
(Supplementary questions: prove that this is the case.  Does this
result generalise to hulls which are not just boxes~?).  The hull of
$[P-Q]$ is obtained by just displacing all of the hull of $P$ by the
coordinates of $Q$.

\cleartabs
\+\cr
\+Pointset:silhouette(Surface:$S$,Point:$Q$)\cr
\+Surface:$S_1$,$S_2$\cr
\+Begin&\cr
\+&If\qquad&[Hull($S$)-$Q$].Nhull($S$) spans 0\cr
\+&Then&If\qquad&Simple($S$)\cr
\+&&Then&silhouette(Approx($S$),$Q$)\cr
\+&&Else&$S_1,S_2$ := Halves($S$)\cr
\+&&&Union(silhouette($S_1$,$Q$),silhouette($S_2$,$Q$))\cr
\+&&Endif\cr
\+&Else&$\emptyset$\cr
\+&Endif\cr
\+End\cr
\+\cr


\end









