Contour Plot algorithm for curved contours

I have been working on plotting contours for a set of 2d points with z-
values. After some research here and on the web, I started with a
Delaunay triangulation where the known points made up the vertices of
each triangle. Using basic linear interpolation I divided each
triangle into polygons for each of the contour lines that crossed the
triangle.

So far so good, I can see the contours and everything makes sense, but
it is too pointy since I'm using linear interpolation in the
triangles.

Next I want to try another interpolation method that would give me
curved contour lines, like inverse distance perhaps.

my main concerns are:
- performance
- understanding the interpolation concept: again, linear was very
basic and very easy to understand.

If any other interpolation method can be based on TINs and give curved
contours that would also work since I already the triangulation part
working.

My background is more to the development side, not so much on graphics
and would appreciate any ideas, suggestions or input you might have on
this project.

Thanks,

dar

<darukan@gmail.com> wrote in message
news:14dcb1a8-e882-4681-bfa4-428e24a3dd35@d77g2000hsb.googlegroups.com...
> If any other interpolation method can be based on TINs and give curved
> contours that would also work since I already the triangulation part
> working.

http://www.geometrictools.com/LibFou...rpolation.html

Look at the code/documentation for the section: "Piecewise quadratic
interpolation of unordered data of the form (x,y,f(x,y)). The code uses
Delaunay triangulation to order the data. The interpolation algorithm is
by Zoltan J. Cendes and Steven H. Wong (a reference is given in the
header file). The resulting interpolation is globally C1 and has local
control."

--
Dave Eberly
http://www.geometrictools.com

Dave Eberly:
> <darukan@gmail.com> wrote in message
> news:14dcb1a8-e882-4681-bfa4-428e24a3dd35@d77g2000hsb.googlegroups.com...
> > If any other interpolation method can be based on TINs and give curved
> > contours that would also work since I already the triangulation part
> > working.
>
> http://www.geometrictools.com/LibFou...rpolation.html
>
> Look at the code/documentation for the section: "Piecewise quadratic
> interpolation of unordered data of the form (x,y,f(x,y)). The code uses
> Delaunay triangulation to order the data. The interpolation algorithm is
> by Zoltan J. Cendes and Steven H. Wong (a reference is given in the
> header file). The resulting interpolation is globally C1 and has local
> control."
>
> --
> Dave Eberly
> http://www.geometrictools.com

Get rid of yours triangles at all. There are better methods, not well-
known, as I could notice. Let's try to imagine a flat table and set of
nails, hammered into it. Each nail's head has its own elevation above
table's surface. take a thin flaexible shhet of steel (tinfoil) and
press it to these nails. After this op each point of this tinfoil will
have some elevation above table's surface. This value one may threat
as interpolation of z-values to the whole 2D surface. To perform this
op one need to solve biharmonic equotation- not a great problem, if
you have not too many points up to 500 approximatelly.

"Andrew_M" <mats@smartfills.com> wrote in message
news:858f6bed-cd9d-4a9e-8f84-e6ff85954b2a@s50g2000hsb.googlegroups.com...
> Get rid of yours triangles at all. There are better methods, not well-
> known, as I could notice.

What is your measure of "better"?

> Let's try to imagine a flat table and set of
> nails, hammered into it. Each nail's head has its own elevation above
> table's surface. take a thin flaexible shhet of steel (tinfoil) and
> press it to these nails. After this op each point of this tinfoil will
> have some elevation above table's surface. This value one may threat
> as interpolation of z-values to the whole 2D surface. To perform this
> op one need to solve biharmonic equotation- not a great problem, if
> you have not too many points up to 500 approximatelly.

How does one go about computing the z-value for this surface given
an (x,y) that is not one of the original data points? Is this an
expensive operation, one that could be noticeable when you have
a very large number of (x,y) to evaluate? I am interested in your
description of the details of solving the biharmonic equation for this
problem.

--
Dave Eberly
http://www.geometrictools.com

Dave Eberly:
> "Andrew_M" <mats@smartfills.com> wrote in message
> news:858f6bed-cd9d-4a9e-8f84-e6ff85954b2a@s50g2000hsb.googlegroups.com...
> > Get rid of yours triangles at all. There are better methods, not well-
> > known, as I could notice.
>
> What is your measure of "better"?

A better metod uses the only info what actually exists. The input is
set of points with z-values, therefore triagularization is odd. One
shoud make interpolation without it

>
> > Let's try to imagine a flat table and set of
> > nails, hammered into it. Each nail's head has its own elevation above
> > table's surface. take a thin flaexible shhet of steel (tinfoil) and
> > press it to these nails. After this op each point of this tinfoil will
> > have some elevation above table's surface. This value one may threat
> > as interpolation of z-values to the whole 2D surface. To perform this
> > op one need to solve biharmonic equotation- not a great problem, if
> > you have not too many points up to 500 approximatelly.
>
> How does one go about computing the z-value for this surface given
> an (x,y) that is not one of the original data points? Is this an
> expensive operation, one that could be noticeable when you have
> a very large number of (x,y) to evaluate? I am interested in your
> description of the details of solving the biharmonic equation for this
> problem.
>
> --
> Dave Eberly
> http://www.geometrictools.com

Course, it is not a trick to calculate z-value if (x,y) is one of
original points. Detailed description take much more than one page, so
if you really interesting in this method you'd better contact me
directly. One needs compute what happens with sheet of tinfoil under
pressure of set of forces, pressing at original points. N original
points- N forces. This is one of typical problems of the theory of
elasticity. To compute these forces one needs solve system of linear
equations. This stage is crucial- I performed this op for 400 hundred
points, but I'm awaiting problems, in N is too large 1000 and more.
But, it means only that this algorithm must be modified for large
numbers of original points. When one knows these forces, the next is
very simple. No triangles, one may compute result everywhere- not only
inside of a convex polygon, drawn around set of original points. You
may download http://www.smartfills.com/Html/GradientFillsDemo.zip .
Unzip and run. Press button number 2, 3 or 4 and you'll see how this
method works. Randomly moving set of vertices, so described above op
is performed to R, G and B color channels independly and you'll see
isohypses for each color channel (when button number 1 is pressed, you
may see result of original interpolation and no isohypses- true
color).

"Andrew_M" <mats@smartfills.com> wrote in message
news:75ec64a3-b837-44ca-8427-dde15b3605fc@z66g2000hsc.googlegroups.com...
> A better metod uses the only info what actually exists. The input is
> set of points with z-values, therefore triagularization is odd. One
> shoud make interpolation without it

The very nature of this problem is that you are trying to generate
information where there is none (i.e. interpolation). Triangulation
is one way to add structure to allow you to interpolate. Your idea
is essentially "thin plate splines". Nothing guarantees that the
sampled data comes from a surface described by a thin plate
spline, so in fact your proposed method also adds structure where
there is none.

I do not think there is any "better" here. The OP must choose an
interpolation algorithm. If the results are acceptable, then all is
well. What you do with the interpolating surface is relevant for
for talking about "better". This was the reason for my asking
how you evaluate your surface. A value of the thin plate
spline at a point depends on *all* the input data points, and
it is relatively expensive to compute (compared to the C1
quadratic surface I proposed). Moreover, you might want the
value to depend only on sample points "near" the point at which
you want the value. Thin plate spline does not allow you to
"cull" in this manner. The C1 quadratic surface does--it has
local control.

> To compute these forces one needs solve system of linear
> equations. This stage is crucial- I performed this op for 400 hundred
> points, but I'm awaiting problems, in N is too large 1000 and more.

This is the classical problem with thin plate splines. The system
of equations is so large that you cannot store all required data
in the computer memory. Most likely you have to page in/out the
data from disk as you use Gaussian elimination on the rows of the
matrix, keeping only a small number of rows in memory at one time.

> No triangles, one may compute result everywhere- not only
> inside of a convex polygon, drawn around set of original points.

This, effectively, becomes an extrapolation process. In most practical
applications, once your input point is "far" from the sampled data, the
interpolation results are not meaningful.

--
Dave Eberly
http://www.geometrictools.com

Dave Eberly:
> "Andrew_M" <mats@smartfills.com> wrote in message
> news:75ec64a3-b837-44ca-8427-dde15b3605fc@z66g2000hsc.googlegroups.com...
> > A better metod uses the only info what actually exists. The input is
> > set of points with z-values, therefore triagularization is odd. One
> > shoud make interpolation without it
>
> The very nature of this problem is that you are trying to generate
> information where there is none (i.e. interpolation). Triangulation
> is one way to add structure to allow you to interpolate. Your idea
> is essentially "thin plate splines". Nothing guarantees that the
> sampled data comes from a surface described by a thin plate
> spline, so in fact your proposed method also adds structure where
> there is none.

I do not see your point of view. The only way if "the sampled data
doesn't come from a surface described by a thin plate spline" is if z-
value is real altitude of sone 3D object (nice sample- bottom of a
ship). Otherwice if z is something like temperature, air pressure,
concentration of some element an so on- not actual vertical
coordinate, this method can be used. As a matter of fact, this method
is limited by the only one factor- difference between two each z must
be much less tnan distance between corresponding points.

>
> I do not think there is any "better" here. The OP must choose an
> interpolation algorithm. If the results are acceptable, then all is
> well. What you do with the interpolating surface is relevant for
> for talking about "better". This was the reason for my asking
> how you evaluate your surface. A value of the thin plate
> spline at a point depends on *all* the input data points, and
> it is relatively expensive to compute (compared to the C1
> quadratic surface I proposed). Moreover, you might want the
> value to depend only on sample points "near" the point at which
> you want the value. Thin plate spline does not allow you to
> "cull" in this manner. The C1 quadratic surface does--it has
> local control.

Remember, how 1D splines arise. They are result of some extremal
problem- to find function with minimal curvature if F(x[i])=Y[i]. It
means that at least in some aspect splines are THE BEST (optimal,
whichever term do your prefer) instrument for 1D interpolation. Same
thing with 2D interpolation. Method, what I'm talking about, is direct
analogue - one needs find function with minimal curvature and if
F(x[k], y[k])=Z[k]. This is why I can tell about it better method.
What do you threat as disadvantage (resulting function is not local-
it depends on the whole set of input data, I can confirm it) is the
only property. Same thing with cubic splines, and noone claims that
zigzag is better due to it offers "local control". This is the nature
of 2D event, nothing more.

>
> > To compute these forces one needs solve system of linear
> > equations. This stage is crucial- I performed this op for 400 hundred
> > points, but I'm awaiting problems, in N is too large 1000 and more.
>
> This is the classical problem with thin plate splines. The system
> of equations is so large that you cannot store all required data
> in the computer memory. Most likely you have to page in/out the
> data from disk as you use Gaussian elimination on the rows of the
> matrix, keeping only a small number of rows in memory at one time.

I do not use Gaussian method, but, yep, this is the problem. I had no
need to solve system with thousands of rows, but there some iteration
methods are possible even without storing matrix of this system.

>
> > No triangles, one may compute result everywhere- not only
> > inside of a convex polygon, drawn around set of original points.
>
> This, effectively, becomes an extrapolation process. In most practical
> applications, once your input point is "far" from the sampled data, the
> interpolation results are not meaningful.

For 2D case extrapolation and interpolation are same thing. There is
no "in" nor "out". Set of point doesn't splits 2D surface like set of
points splits a line.

Andrew_M wrote:
>
> I do not see your point of view. The only way if "the sampled data
> doesn't come from a surface described by a thin plate spline" is if z-
> value is real altitude of sone 3D object (nice sample- bottom of a
> ship). Otherwice if z is something like temperature, air pressure,
> concentration of some element an so on- not actual vertical
> coordinate, this method can be used. As a matter of fact, this method
> is limited by the only one factor- difference between two each z must
> be much less tnan distance between corresponding points.
>

It's clear that your first sentence is absolutely correct.

--
Kenneth Sloan KennethRSloan@gmail.com
Computer and Information Sciences +1-205-932-2213
University of Alabama at Birmingham FAX +1-205-934-5473
Birmingham, AL 35294-1170 http://KennethRSloan.com/

Kenneth Sloan:
> Andrew_M wrote:
> >
> > I do not see your point of view. The only way if "the sampled data
> > doesn't come from a surface described by a thin plate spline" is if z-
> > value is real altitude of sone 3D object (nice sample- bottom of a
> > ship). Otherwice if z is something like temperature, air pressure,
> > concentration of some element an so on- not actual vertical
> > coordinate, this method can be used. As a matter of fact, this method
> > is limited by the only one factor- difference between two each z must
> > be much less tnan distance between corresponding points.
> >
>
> It's clear that your first sentence is absolutely correct.
>
Thanks for comfirming my ideas. Method what I was talking about can be
used not only for interpolation but for approximation as well. The
only problem is to make some practical approach to full procedure-
lage dataset produces huge matrix, what causes lots of problem. This
is interesting task, and I'm currently working with. One small
snapshot may be downloaded from http://www.smartfills.com/Html/images/isohypses.png
.. Random set of point with random z-values