C or C++ bicubic surface fit code ?
>mag<
Does anyone know of C or C++ source for
a bicubic spline surface fitter ? Or something
even better ? I need a surface fitter able
to fit a sparse, unordered 2D data set and the
surface must fit through the data points.
Some may say the nature of the fitter depends on a number of
factors, but I've seen a Fortran implementation
of a bicubic spline surface fit which looks
exactly like what I need. As long as I wind up with
a reasonable (use your imagination) surface, I can
probably use it. I need vertices, quads or triangles though,
not just sets of spline coefficients.
Any email, pointers, references etc would
be greatly appreciated. :)
many thanks
Mike
ras...@bayarea.net
Pierre Alliez
>a bicubic spline surface fitter ? Or something
>even better ? I need a surface fitter able
>to fit a sparse, unordered 2D data set and the
>surface must fit through the data points.
Do you have a coarse seed triangular mesh ?
You may try the Loop algorithm if you have, see example :
>mag<
I'm looking for C/C++ bivariate surface fit code. I've posted a similar
request previously but the lack of responses lead me to believe I had
not well-stated my request, or perhaps it is indeed something difficult
to find, so I apologize if this is redundant for some, and many thanks
to those who are able to reply, esp. via email as my ISP removes
messages
every day or so.
Puzzle:
Given a set of data points randomly positioned over a 2D area and
of random value (or height, as you prefer)....
** How do I generate a regular mesh (say, 50x50 etc) of polygons
or triangles such that the surface is
A) a bicubic interpolation *through* the points with the
"tightest"
curve, and/or "least" warped by the choice of coefficients
* again, I'm looking for that regular 50x50 mesh... I did
get
a couple of emails from people prescribing tesselators,
surface
subdivision etc, but that isn't what I had requested or
what I'm
looking for since, unfortunately, it doesn't fit the data
as prescribed :)
and/or
B) well approximated as near the points as possible with "good"
curvature, and perhaps by selection of some limit of error,
eg. possibly using least squares and I could only guess as
what kind of rational spline to use, whether bicubic or
whatever..
I am guessing that "local" interpolation/approximation is
computationally
less expensive, and if I had to choose between local/global, I guess I
would
choose local in this case. Note that the data set could vary from a
very few
points (may 10-20) to a fairly large number (1000-2000), and I would
guess
that a global solution would freeze up even a high end NT with a global
matrix to solve of that size, but that's just a guess :)
I've just barely started figuring out the basics of NURBS, and have
bought the book "The Nurbs Book" which seems great, but I've only
finished the first two chapters. My math is rusty, but I'll probably
eventually
figure it out, but I'm hoping to actually not become a Nurbs expert
unless
I really have to :P I guess that sounds lazy, but with all the other
tasks on
my plate, it's one more apple I guess I would rather save for later to
savor
rather than to sweat over :)
I hope this is a better informed request than my previous one :)
many thanks for email, pointers etc
Mike
ras...@bayarea.net
Sorry for the cross-posting, but the range of people who would
have encountered this problem isn't strictly limited to any of the
groups posted to.
bob_wil...@my-deja.com
ras...@bayarea.net wrote:
>
> Does anyone know of C or C++ source for
> a bicubic spline surface fitter ? Or something
> even better ? I need a surface fitter able
> to fit a sparse, unordered 2D data set and the
> surface must fit through the data points.
You might try looking at Dave Eberly's site:
<http://www.magic-software.com/>
Check under "interpolation".
Bob
Sent via Deja.com http://www.deja.com/
Before you buy.
bob_wil...@my-deja.com
Is there a reason you have to form a regular rectangular mesh?
Forming the Delaunay triagulation and then interpolating that
triangulated network would be less costly in both time and memory.
I have some triangulation and linear interpolation code in C that
you might be interested in.
If you need to do a spline interpolation, again I would encourage
you to check out Dave Eberly's site. Paul Bourke has some helpful
source code too, I think:
http://www.swin.edu.au/astronomy/pbourke/
If you're interested in my code, please reply by email.
>mag<
bob_wil...@my-deja.com wrote:
> Is there a reason you have to form a regular rectangular mesh?
Hmmmm... good question actually.. The main reason is that a
gridded mesh "looks" nicer if viewed only as lines, although I'm
supporting both textured/shaded polys and lines... The only other
reason which may not be a valid one is that I want to be able to
take the resultant mesh and easily smooth it to zero around a
specified border. For example, if the ground plane has an image
of a disk, I would smooth (only near the matching border) the
mesh so that it is zero at the boundaries of the disk... not sure if
that is the right way to do it... there's a couple of other reasons,
but could be equally invalid... So, to be honest, I don't really
know if it is essential or not :P
> If you need to do a spline interpolation, again I would encourage
> you to check out Dave Eberly's site. Paul Bourke has some helpful
> source code too, I think:
> http://www.swin.edu.au/astronomy/pbourke/
thanks for the code references. I've looked at Dave's stuff which is
really cool, but I'm not sure I see the answer in there. It's on my
list of "to do" items though to try out all of his interpolation code :)
thanks again
Mike
P.G.Hamer
> Hi :)
>
> bob_wil...@my-deja.com wrote:
>
> > Is there a reason you have to form a regular rectangular mesh?
>
> Hmmmm... good question actually.. The main reason is that a
> gridded mesh "looks" nicer if viewed only as lines, although I'm
> supporting both textured/shaded polys and lines...
In which case calculating the curved surface based on Delaunay
triagulation, then calculating the points on the grid before display
(of flat plates?) sounds reasonable. I think some of the netlib
TOMS routines already mentioned will calculate the values on
such a grid.
Of course a bicubic interpolation from the grided points would not
produce /precisely/ the same curved surface as that produced
from the Delaunay triagulation; or /exactly/ pass through the
original set of points. However the differences would probably
be unobservable in the graphical presentation anyway.
Peter
>mag<
Thanks for the reply !
"P.G.Hamer" wrote:
> ...... However the differences would probably
> be unobservable in the graphical presentation anyway.
Will that be true even if I display the original data points ? :)
One of the criteria I have is to display the original data points with
the surface.
Thanks to the groups for all of the interesting replies.
One mentioned that since I will convert the spline surface
back to triangles or quads that I'm going to have the flat facets
and not the nice smooth that one would usually want from
a spline... true, but in this case, I'm not going to do the final
rendering with splines, I just want the surface smoothly fit
to the data points, although not linearly , which some keep suggesting
for some reason. The idea is that I can always increase the mesh
resolution to get a smoother and smoother surface, and as the resolution
increases, it approaches the curvature of the spline from which it is
derived. I'm not sure if I'm approaching this the wrong way :)
But, hopefully I can derive a solution from all of the great suggestions
from all you guys.
thanks
Mike
bob_wil...@my-deja.com
> > If you need to do a spline interpolation, again I would encourage
> > you to check out Dave Eberly's site. Paul Bourke has some helpful
> > source code too, I think:
> > http://www.swin.edu.au/astronomy/pbourke/
>
> thanks for the code references. I've looked at Dave's stuff which is
> really cool, but I'm not sure I see the answer in there. It's on my
> list of "to do" items though to try out all of his interpolation
code :)
You ought to follow up on your "to do" list. Dave does have a nice
spline interpolation routine that you can apply to the triangulation.
Determining the heights along the regular grid becomes just a
matter of plugging in an X & Y value, determining the triangle in
which a given coordinate resides and using the parameters of that
triangle to solve for Z.
>mag<
bob_wil...@my-deja.com wrote:
> You ought to follow up on your "to do" list. Dave does have a nice
> spline interpolation routine that you can apply to the triangulation....
Of course, you are correct :)
The only reason it was on the back-burner is because I'm still trying
to understand splines well enough to know what I would be looking
at... :P I finally figured out how to compute the geometry vector
and basis matrix and so I'm now trying to understand B-Splines. At
least now I better understand what's involved to some degree (or
is that to degree 3 ?) :)
Your suggestion of plugging in x and y makes intuitive sense, and
I guess it will make perfect sense as soon as I understand how splines
are attached to one another, since I need to know that in order to
get to the point of being able to plug in x and y over the entire surface
:)
But in the mean time, perhaps you can point me at a specific function in
Dave's magic box ? :) The only one I see at the moment is meshintp,
which seems to use quadratics and also seems to require partial derivatives
in the setup, and it's only C1 continuous. Am I wrong to assume that
a bicubic fit is C2 and I can get by with only specifying endpoints and
tangents for any given patch ? This is all relatively new to me :)
many thanks
Mike
P.G.Hamer
> Hi, :)
> Thanks for the reply !
>
> "P.G.Hamer" wrote:
>
> > ...... However the differences would probably
> > be unobservable in the graphical presentation anyway.
>
> Will that be true even if I display the original data points ? :)
> One of the criteria I have is to display the original data points with
> the surface.
>
That depends on the size of your grid. You could always do a calculation
to evaluate the maximum error (for any of your original points) fo a given
grid.
OTOH I suspect that we are very poor at perceiving the `true' position
of a point in a 3D picture. So you may well be OK for all realistically
sized grids.
Peter
BTW if you ever intend to display contour lines Wahba had a very
good paper some years ago. You might start at
http://www.stat.wisc.edu:80/~wahba/
P.G.Hamer
> BTW if you ever intend to display contour lines Wahba had a very
> good paper some years ago. You might start at
> http://www.stat.wisc.edu:80/~wahba/
I should have said:
This paper discusses how to create contour lines free from `features'
not present in the data. It is not simply about the mechanics of drawing
contours.
Peter
>mag<
Thanks for the reference, but unfortunately I am unable to access that URL,
the page generates a server error (at least just now it did).
What I'm doing now is something I found in the Siggraph 99 papers
which just uses 1D textures which map pixels to the texture according
to their respective "height". I guess this is about as good as I can
get if I use textures.
I thought about another method which would generate separate contour
"lines" which hover just about the surface, but that seems to be a sticky
problem. i.e. I would have to determine which direction to offset the
contour line, and by how much (just enough so that it doesn't get obscured
by the surface mesh itself, and not so much that it becomes nearly
irrelevant).
I would definitely be interested in any references which might have
solved that problem.
In the mean time, I would be interested to know if anyone thinks the
1D texture method is inaccurate... but my intuition tells me that it cannot
be any more inaccurate than the surface itself.
BTW, my initial mesh is only 50x50, but I'm allowing users to extend
the mesh to any size so that users with higher end systems can have
as fine a resolution as they wish to try...
Thanks again for all the replies :)
Mike
PS
I'm slightly stuck on Berstein functions... I could go on faith as to
why they are calculated the way they are and just go forward with
my reading about splines, but if it isn't too complex a matter I would
like to understand how/why they are chosen the way they are :)
Thanks again :)
bob_wil...@my-deja.com
> Of course, you are correct :)
> The only reason it was on the back-burner is because I'm still trying
> to understand splines well enough to know what I would be looking
> at... :P I finally figured out how to compute the geometry vector
> and basis matrix and so I'm now trying to understand B-Splines. At
> least now I better understand what's involved to some degree (or
> is that to degree 3 ?) :)
>
> Your suggestion of plugging in x and y makes intuitive sense, and
> I guess it will make perfect sense as soon as I understand how splines
> are attached to one another, since I need to know that in order to
> get to the point of being able to plug in x and y over the entire
surface
> :)
>
> But in the mean time, perhaps you can point me at a specific function
in
> Dave's magic box ? :) The only one I see at the moment is meshintp,
> which seems to use quadratics and also seems to require partial
Hello Mike,
I think I understand what you're trying to do, but let me
repeat it back to you just to make sure. You want to:
1) take an arbitrary set of 2D points (each with an "offset"
value which will be graphically represented as a height),
2) triangulate these points,
3) derive splines for each triangle such that:
a) the the original points are on the spline surfaces,
b) the surfaces are continuous across triangle edges and
c) the first derivative is continuous across all edges.
4) slice an dice this smooth surface into a regular square
grid, and
5) plot the mesh.
If that's it, then I think Dave Eberly has written most of
the necessary code in his qdintp2.cpp and qdintp2.h files.
That'll get you up to (but not including) the "slicing and
dicing". Check out Dave's explanation in meshintp.pdf.
He cites two sources by David F. Watson. I found the
"Contouring..." book at my local library. Very helpful.
http://www.magic-software.com/gr_intp.htm
// Watson, D.F., 1981, Computing the n-dimensional
// Delaunay Tessellation with Application to Voronoi
// Polytopes:
// The Computer Journal, 24(2), p. 167-172.
//
// Additional information can be found in
//
// CONTOURING: A Guide to the Analysis and Display
// of Spatial Data, by David F. Watson, Pergamon Press,
// 1992, ISBN 0 08 040286 0
> and it's only C1 continuous. Am I wrong to assume that
> a bicubic fit is C2 and I can get by with only specifying endpoints
and
> tangents for any given patch ? This is all relatively new to me :)
To be perfectly honest, I don't know what C1 or C2 continuity
means. First and second derivative continuity perhaps? Sorry,
but you've stumbled into one of my vast areas of ignorance. Can
anyone else elucidate? If you're dealing with a triangulation,
then I don't believe the term bicubic applies. (Feel free to
correct me, anyone.)
P.G.Hamer
> Hi,
> Thanks for the reference, but unfortunately I am unable to access that URL,
> the page generates a server error (at least just now it did).
Still seems to work for me.
> I thought about another method which would generate separate contour
> "lines" which hover just about the surface, but that seems to be a sticky
> problem. i.e. I would have to determine which direction to offset the
> contour line, and by how much (just enough so that it doesn't get obscured
> by the surface mesh itself, and not so much that it becomes nearly
> irrelevant).
> I would definitely be interested in any references which might have
> solved that problem.
You might try the technique used in several graphics packages (eg
the free gnuplot) where a 2D contour plot is shown below the 3D surface;
eg on the z=0 plane.
Peter
P.G.Hamer
surface modelling in
http://iris.usc.edu/Outlines/Papers.html
especially as it considers discontinous functions.
Peter
>mag<
Thank you for the great references. I've been reading through
them and have found some good hints.
I've just complete the chapter on Beziers, so I think I have
a handle on what to do if I decided to try to implement
a bi-cubic Bezier patch. I'm still not sure if Delaunay
triangulation is going to help or hinder my getting a regular
grid of points which all lie in a bicubic surface which well
represents the original data set. I need to read up on Delaunay
triangulation and see what it does besides just finding a set
of triangles to map into the data points (which is my current
presumption, although this is merely conjecture without having
investigated it).
Several replies keep referring to Dave Eberly's libraries, but
just so I don't have to repeat myself again, I am under the
impression that it does not take a random data set and convert
it into a regular gridded mesh (eg. 50x50 or whatever resultant
size somebody wants) but I may be misunderstanding Dave's
documentation. I am a big fan of his stuff however, and I have
no doubt he knows a great deal more than I could pretend to
know, but all I can gather is from my current knowledge and
with that, from what I can glean from his documentation. Also,
his one (are there more?) surface fit algorithm is based on
quadratics, not cubics, but I would assume that with a little
more insight I could convert that to a bi-cubic fit.
So maybe I should break the question into 3 parts... even though
all are related.
1) given a random set of data points which is desired to fit
through a bi-cubic surface, where do I begin ?
A) triangulation ?
B) some kind of regression and intermediate interpolation
to determine approximate regular grid points ?
C) other ?
2) If the answer to (1) is (A) triangulation, then I presume
I go directly to surface fit, which is farther back in
the books I'm reading so I can't even guess at this yet but
I would guess that triangulation will fall out at the end
just because the points start off in some regular grid.
If the answer to (1) is (B) then I would make a similar
assumption, and still assume there is a surface fit to do
3) If the answer is "they're all related" in order to do a bi-cubic
fit, then I can see I've got still more reading to do :)
This has kind of been my presumption so far, and that
in this case I'm guessing that triangulation falls out from
the surface fit in the end, and that the surface fit generates
coefficients which I then apply to somehow determine
the regular grid points.
Does anyone know if any of these are on the right track ? :P
Sorry if I'm being a bit dense here... again, splines, nurbs etc
are kind of new to me.
Thanks :)
Mike
>mag<
tangents in the geometry matrix. I hate to use a page reference
since there are probably a dozen versions of the book, but
I have Foley & van Dam, Edition 2, the C Version.
On page 484, reference (11.12) you see a geometry
matrix Gh = [P1 P4 R1 R4]^T (sorry for the crappy notation)
and the implementation is carried out by simply computing
(for R1 and R4 which are tangent vectors) by substituting
derivative of the basis function and plugging in t. When I
read everything through the first time, the geometry matrix
made intuitive sense, but when I read through the second
time, it didn't make sense that a geometry vector
could randomly go between presuming use of a function
and it's derivative. I don't know why this bothers me, perhaps
because I'm not reading it correctly.
So, I'm just now going on faith that "it just works"..
Questions:
1) can I presume I could plug in all four as points in this
example and not use tangents ?
2) it would seem that position in the matrix doesn't matter
since content doesn't seem to matter
3) exactly what are the limitations of what can be used
in the geometry matrix in this case ?
Since any polynomial for these examples is defined by C * T
(where * is "dot") then C is just coefficients, and of course
it must contain geometry information along with any other
information relative to the transformation... . so F & vD
break C into M * G where M is the basis matrix and G is
the geometry matrix. I know the answer to my above
questions lie in what comes next, but I'm missing it....
These matrices are re-combined as T * M * G so that
the blending functions are given by B = T * M.
So, we want to find the basis matrix given a predetermined
T matrix, so we substitute absolutes for the
T matrix and take it's inverse and use that as the new
blending matrix Mb...
Ok... that gets us up to (11.19) in F&vD which is just
a figure depicting a sample of this matrix.
But, the book goes on with T * Mb * G and it seems
to me that Mb is partially based on derivatives of part
of the T matrix, so how can the two (T & Mb) function
correctly as shown ?
What am I not understanding ?
I've logged this as "it just works" and gone on through
Bezier's, but would sure like to understand this
mystery :) I think then I might be able to solve my
original problems, or at least to understand the coming
text.
A final note:
What finally brought a grin to my face was discovering
that hand computing the general equation for a cubic
Bezier agreed with F & vD after their jumping through
hoops with matrices... all the time I thought I was computing
the polynomials which would *then* be broken into C * T to
be used to determine the final blending functions. My head
is spinning :P
Thanks to you all for your patience and great input.
Peter Spellucci
given a set a random points (x_i,y_i,z_i) with z_i representing some f(x_i,y_i)
fit them by a C2-bicubic surface, then represent this surface by a regular
grid .
(solution: in principle as follows: 1) compute triangulation. 2)compute
C2 -continuous fit (i am not a specialist, but darkly remember this will
have order much higher than bicubic). 3) compute surface values on a regular
covering grid. 4) compute
bicubic C2 spline from these data (will of course result in a different
final surface).
just a question:
why not taking the following, much simpler approach:
1) choose a regular grid covering your data points. (there must be at
least as much data points as there are grid points) 2)set up a bicubic
(tensor model) C2 spline on the regular grid. with natural end conditions,
it is uniquely determined by the values on the regular grid. 3) now make a
linear least squares fit for fitting your data points by this so defined surface.
this allows you to use standard software. directly does what you want (with
the exception that it may not interpolate your data, but this
property is also destroyed by step 3 above).
peter
P.G.Hamer
> Hi Peter
> Thank you for the great references. I've been reading through
> them and have found some good hints.
>
> I've just complete the chapter on Beziers, so I think I have
> a handle on what to do if I decided to try to implement
> a bi-cubic Bezier patch. I'm still not sure if Delaunay
> triangulation is going to help or hinder my getting a regular
> grid of points which all lie in a bicubic surface which well
> represents the original data set. I need to read up on Delaunay
> triangulation and see what it does besides just finding a set
> of triangles to map into the data points (which is my current
> presumption, although this is merely conjecture without having
> investigated it).
AFAIK Delaunay triangulation not only produces a set of non-
overlapping triangles, it also picks nicely shaped triangles. That
is, it avoids choosing long thin triangles wherever possible.
> So maybe I should break the question into 3 parts... even though
> all are related.
>
> 1) given a random set of data points which is desired to fit
> through a bi-cubic surface, where do I begin ?
Several posters, including myself, have suggested that you
use bi-cubic >patches< with continuity conditions across the
joins. This is not a bicubic surface over the whole area, just
over each patch [ie triangle in a Delaunay triangulation]. Code
is in netlib/toms somewhere around the high 600's. I think
its O(n) in computational complexity.
I think others are suggesting splines over the whole surface.
Which also seems to be closer to the area you are reading in
at the moment.
> A) triangulation ?
Essential if you want to fit bicubic patches to the triangles! After
calculating the patches you can use them to obtain the values
at your chosen regular grid points.
> B) some kind of regression and intermediate interpolation
> to determine approximate regular grid points ?
I assume that you meant `to determine the values at your chosen
regular grid points'
Yes, this is another way. You might try the Generalised Additive Model.
The drawback is that it corresponds to [something like] a least-squares
fit to your points.
If your points have measurement errors this could well be better than
trying to fit a surface throught them -- but I understood that your
points were error-free.
If the points are error free, and your function is nice and smooth, then
this approach is quite reasonable. I would expect the errors of surface
fit at the data points to be larger than using method (A) though. Maybe
not significantly bigger.
> C) other ?
Set up your regular grid and do a penalised regression fit to the
data points. Here you have more unknowns than knowns- hence an
under-etermined problem. But techniques such as GCV (generalised
cross validation) can handle this sort of situation. The 1D case is
explored in detail in a paper by Wahba called something like
`the automated french curve' and also in her SIAM monograph
on splines (see http://www.siam.org).
But I think this is getting /far/ to complicated for your problem, and
still smooths rather than interpolates the surface through your points.
In another posting Peter Spelluci suggests that you choose the
grid size so that a least-squares has no more unknowns than
knowns. IMHO this isn't what you want; as the grid would be
to coarse your visualisations (as far as I understand your
problem).
I'm very nervous about querying Peter's expertise, but I suspect
if you aim for a just-determined system (by having the grid as fine
as possible while avoiding under-determinacy) you might end up with
a poor fit away from the points [cf why you don't usually interpolate
N points with a (N-1)th order polynomial].
> 2) If the answer to (1) is (A) triangulation, then I presume
> I go directly to surface fit, which is farther back in
> the books I'm reading so I can't even guess at this yet but
> I would guess that triangulation will fall out at the end
> just because the points start off in some regular grid.
Do your data points start off on a regular grid? So far I've assumed
that they are positioned at random. Please elaborate.
> If the answer to (1) is (B) then I would make a similar
> assumption, and still assume there is a surface fit to do
> 3) If the answer is "they're all related" in order to do a bi-cubic
> fit, then I can see I've got still more reading to do :)
> This has kind of been my presumption so far, and that
> in this case I'm guessing that triangulation falls out from
> the surface fit in the end, and that the surface fit generates
> coefficients which I then apply to somehow determine
> the regular grid points.
You again say that `the surface fit generates coefficients which
I then apply to somehow determine the regular grid points'. This
doesn't fit with my picture of what you're doing. Which is to
determine the values of the interpolating surface >on a regular
grid chosen by you to satisfy visualisation criteria<.
How does a surface fit lead you to determine the grid structure?
Rather than the values on the already chosen grid.
> Does anyone know if any of these are on the right track ? :P
> Sorry if I'm being a bit dense here... again, splines, nurbs etc
> are kind of new to me.
I think that you are getting confused between your options. There
are several approaches, all of which you would be well advised to
use public domain code for, rather than writing your own. [Or if you
insist, all of which you should first get working with PD code before
you try to implement all the detailed algorithms yourself.]
The bicubic patch approach several of us have suggested has little
to do with the general literature on Beziers splines, nurbs, etc. It
will
be described in the ACM Transactions on Mathematical Software
article referenced in the netlib/toms code.
Peter
>mag<
A very interesting and thoughtful response :)
Peter Spellucci wrote: (some clipped for brevity by >mag<)
> the problem posed:
> given a set a random points (x_i,y_i,z_i) with z_i representing some f(x_i,y_i)
> fit them by a C2-bicubic surface, then represent this surface by a regular
> grid .
> (solution: in principle as follows: 1) compute triangulation.
So, by this do you mean to triangulate the existing points as they
are and not yet to worry about whether they are grid-aligned ?
> 2)compute
> C2 -continuous fit (i am not a specialist, but darkly remember this will
> have order much higher than bicubic).
I've just realized in the last few hours that it may not be obvious to
someone viewing the mesh whether it is C1 or C2, but I'll stick
with possible C2 fit if possible just in case, and in case there is
a solution. If I understand Foley & van Dam, some bicubic patches
are C2 continuous, but not all certainly.
> 3) compute surface values on a regular
> covering grid.
This is the part that's confusing me the most... I had presumed
I would have to somehow combine (2) and (3) because from
what I know (so far, and that isn't very long...) I need a regular
grid to compute the fit :P and yet I know intuitively that somehow
making a fit of data "as is" will provide the polynomials to better
compute the points of the regular grid...
> 4) compute
> bicubic C2 spline from these data (will of course result in a different
> final surface).
But wasn't that what we did above in (2) ? :P
Ok... so, the fit isn't necessarily related to computing the splines
from which the final mesh is constructed... I guess it's time
to walk through Dave Eberly's quadratic code just to see how
that works... if it does a fit from randomly placed data, that is...
> just a question:
> why not taking the following, much simpler approach:
> 1) choose a regular grid covering your data points. (there must be at
> least as much data points as there are grid points)
The grid may be 50x50 or larger which is only to determine
final viewable mesh resolution, and the data points could be
anywhere from 1 to several thousand and have no relationship
to the grid other than the fact that they are contained within
the boundaries of the grid base... not sure if that helps answer
that question...
But, assuming I don't have enough data points to cover the grid,
exactly how do I make a grid cover my data points ? What do
I put in those "cells" not containing data points ? What do I
do with the cells containing multiple data points ? (I've presumed
so far just to average multiple data points in a cell)
> 2)set up a bicubic
> (tensor model) C2 spline on the regular grid. with natural end conditions,
> it is uniquely determined by the values on the regular grid.
I haven't studied tensor models yet, but now that you've pointed it out
I'll add that to my pile of math books :)
> 3) now make a
> linear least squares fit for fitting your data points by this so defined surface.
> this allows you to use standard software. directly
by "standard software" what do you mean exactly ? I'm guessing that
what you mean is to convert the least squares fitted (and gridded) data
directly to a viewable mesh ?
> does what you want (with
> the exception that it may not interpolate your data, but this
> property is also destroyed by step 3 above).
Is that true ? That is interesting. I will have to definitely examine
this approach.
Many thanks for your reply :)
Mike
P.G.Hamer
netlib/toms/684
COMPUTATION OF A RECTANGULAR GRID FROM SCATTERED DATA
USING A C2 NONIC POLYNOMIAL REPRESENTATION OVER TRIANGLES
<sorry for the caps, but this is FORTRAN>
You may also be interested to look at
netlib/dierckx
DIERCKX is a package of Fortran subroutines for calculating smoothing
splines for various kinds of data and geometries, with automatic knot
selection.
This library is also called FITPACK, but is independent of the FITPACK
library by Alan Cline.
Note in particular surfit.f
c given the set of data points (x(i),y(i),z(i)) and the set of positive
c numbers w(i),i=1,...,m, subroutine surfit determines a smooth bivar-
c iate spline approximation s(x,y) of degrees kx and ky on the rect-
c angle xb <= x <= xe, yb <= y <= ye.
c if iopt = -1 surfit calculates the weighted least-squares spline
c according to a given set of knots.
c if iopt >= 0 the total numbers nx and ny of these knots and their
c position tx(j),j=1,...,nx and ty(j),j=1,...,ny are chosen automatic-
c ally by the routine. the smoothness of s(x,y) is then achieved by
c minimalizing the discontinuity jumps in the derivatives of s(x,y)
c across the boundaries of the subpanels (tx(i),tx(i+1))*(ty(j),ty(j+1).
c the amounth of smoothness is determined by the condition that f(p) =
c sum ((w(i)*(z(i)-s(x(i),y(i))))**2) be <= s, with s a given non-neg-
c ative constant, called the smoothing factor.
c the fit is given in the b-spline representation (b-spline coefficients
c c((ny-ky-1)*(i-1)+j),i=1,...,nx-kx-1;j=1,...,ny-ky-1) and can be eval-
c uated by means of subroutine bispev.
Peter
>mag<
didn't get attached to the thread. It will make no sense probably to most
of those
outside of this thread...
---------------
>mag<
Well, I guess Fortran wins. I've tried to solve the problem by looking
for C/C++ (per my subject line) and have already a lot of netlib stuff and
others downloaded, but I'm not a Fortran guy, and it somehow didn't
make sense to go back and re-learn Fortran, at least at the time. But if
that's where the answers lie, then that's where I'll be :)
Initially, I started this thread exactly because my company has some old
Fortran code and I didn't want to have to dig into that stuff, esp. since
the app we're building is all VC++. It made no sense to try to maintain
old code that nobody knew anything about :)
BTW, I did try running f2c on some of those Fortran files, but where
the conversion worked, the output was nearly as cryptic as the input :)
and more often, f2c wouldn't do the conversion... just a side note :)
thanks again to all
Mike
P.G.Hamer
> Thanks Peter,
> Well, I guess Fortran wins. I've tried to solve the problem by looking
> for C/C++ (per my subject line) and have already a lot of netlib stuff and
> others downloaded, but I'm not a Fortran guy, and it somehow didn't
> make sense to go back and re-learn Fortran, at least at the time. But if
> that's where the answers lie, then that's where I'll be :)
Sad but sometimes true!
Of course the ACM TOMS article will explain the method used. So you
could find a C/C++ Delaunay triangulation routine and recode the patch
fitting (using a linear maths library if/where appropriate) and mesh
generation.
AFAIR the core of the patch fitting algorithm was some tens of equations
generated by a symbolic maths package. I would expect these to be
pretty obvious in the Fortran code, and pretty easy to convert to C/C++.
There is also the possibility that a search based on keywords from
the Fortran documentation might find something in C/C++.
> BTW, I did try running f2c on some of those Fortran files, but where
> the conversion worked, the output was nearly as cryptic as the input :)
> and more often, f2c wouldn't do the conversion... just a side note :)
Old Fortran code is notorious for being non-standard. Different compilers
had different `extensions' and often non-standard intrinsic functions.
OTOH when f2c'ing core algorithms (ie not i/o based user interfaces)
I've usually managed to get things to work in a hour or so. Only needing
to change a handfull of lines (not counting duplicates of `the same'
change.
Peter
FYI
Abstract of the ACM TOMS paper
http://www.acm.org/pubs/citations/journals/toms/1990-16-3/p246-preusser/
Algorithm writers home page
http://www.fhi-berlin.mpg.de/grz/pub/preusser.html
C Delaunay triangulation code?
http://www.netlib.org/voronoi/
Dave Eberly
news:38C765A7...@bayarea.net...
>
> Several replies keep referring to Dave Eberly's libraries, but
> just so I don't have to repeat myself again, I am under the
> impression that it does not take a random data set and convert
> it into a regular gridded mesh (eg. 50x50 or whatever resultant
> size somebody wants) but I may be misunderstanding Dave's
> documentation.
A couple of my files can take a random data set and fit
with a surface that contains the original points.
* qdintp2.{h,cpp}: Takes random data of form (x,y,f(x,y)),
computes a Delaunay triangulation of the (x,y) portions,
then computes a globally C1 quadratic fit with local
control. You have the freedom to choose the first-order
partial derivatives at the input points. These can be
estimated by measuring variation with neighboring
points in the triangulation. Or you can just go with zero
for all the partial derivatives.
* tpspline.{h,cpp}: Uses thin plate splines for fitting points
of the form (x,y,f(x,y)). The fit has C-infinity smoothness.
> Also,
> his one (are there more?) surface fit algorithm is based on
> quadratics, not cubics, but I would assume that with a little
> more insight I could convert that to a bi-cubic fit.
You say you want a "bi-cubic fit". The term bicubic is used
to denote a tensor product patch. That would imply that you
already have your data points in a rectangular lattice. This
does not appear to be the case, otherwise you would have
just fit with the standard bicubic patch. Are you looking for
a mesh of triangle patches, each patch being cubic?
My impression from the thread is that you are trying to
1. build a surface that contains the random points
2. resample to get a rectangular grid of points
3. fit the rectangular grid to get a bicubic patch surface
If this is the case, then I do not understand why you should
be concerned about the order of continuity in step 1. In
fact, it is not clear why you should be concerned about
the original surface containing the random points. Somehow
I get the impression you believe that by having the original
surface contain the data points that you have a "better fit"
than with a surface that does not contain the points. If so,
how are you measuring the goodness of fit?
--
Dave Eberly
ebe...@magic-software.com
http://www.magic-software.com
>mag<
Dave Eberly wrote:
> A couple of my files can take a random data set and fit
> with a surface that contains the original points.
Hey :) The master speaks !
> * qdintp2.{h,cpp}: Takes random data of form (x,y,f(x,y)),
> computes a Delaunay triangulation of the (x,y) portions,
> then computes a globally C1 quadratic fit with local
> control. You have the freedom to choose the first-order
> partial derivatives at the input points. These can be
> estimated by measuring variation with neighboring
> points in the triangulation. Or you can just go with zero
> for all the partial derivatives.
So how do I know the first derivative of a function I haven't
seen yet ? I'm sure you know from reading the thread that
I'm just recently trying to understand this stuff :)
> * tpspline.{h,cpp}: Uses thin plate splines for fitting points
> of the form (x,y,f(x,y)). The fit has C-infinity smoothness.
Again, I'm a little stuck on the notion of supplying an input
function :P
> You say you want a "bi-cubic fit". The term bicubic is used
> to denote a tensor product patch.
But, I understood bi-cubic to denote a cubic polynomial or rational
equation over 2 orthogonal curve sets... I see "cubic curves" vs
"bicubic patches" and have not yet seen a definition of "bicubic",
not in Foley & Van Dam, not in Watt, not even in "The Nurbs Book"...
so that probably means it's somewhere back in my old math
books :P Not disputing you, but trying to understand what goes
where and why....
> That would imply that you
> already have your data points in a rectangular lattice.
Does it still imply that w.r.t. the home-made definition I have
above ? Or am I even close ?
> This
> does not appear to be the case,
Exactly :)
> otherwise you would have
> just fit with the standard bicubic patch.
yes, I'm trying to figure out how to get a random set of
data points somehow converted eventually to an n x m
mesh such that the mesh "fits" the data, hopefully
with C2 continuity. I've made a lot of assumptions and
some modifications to my assumptions as well along the
thread to try to seek the answers as you've also probably seen.
I think I said several times I'm trying to figure this stuff
out and some folks have given me a little more credit
than I deserve I think presuming I know what I'm doing...
For example, someone suggested looking at your code,
and I found the interpolation folder easily enough, and
read the PDF's, but from what I could tell, the PDF's
referenced 3 unique solutions for that folder, and it wasn't
clear to me that it was fitting "n" random points to a grid, etc....
I guess the point is, yes I'm looking for code, but I'm also
trying to understand as I go along. I've just finished reading
the Bezier chapters and on to B-Splines, and while the chapters
on surfaces regarding Bezier's seem "intuitive", there's still
a lot of holes in the process for me.
> Are you looking for
> a mesh of triangle patches, each patch being cubic?
That sounds right... what I've had my radar locked onto
was a "local" cubic fit, with the final result being the
triangulated mesh, where the points of the triangles are
coincident with a regularly spaced 2D grid as seen
"from above"... :) Now whether I've got my terminology
right I don't know, but the "fuzzy" spot on the radar
seems like the right direction and sounds like what you're saying :)
What I think I understand of the material I've read so far on
splines etc, a cubic patch would (by definition ?) provide better
continuity than a quadratic patch, or am I misunderstaning the
concept ? i.e. a polynomial or rational equation of degree 3 vs
a polynomial of degree 2.. is that a fair comparison ?
> My impression from the thread is that you are trying to
> 1. build a surface that contains the random points
> 2. resample to get a rectangular grid of points
> 3. fit the rectangular grid to get a bicubic patch surface
Or any order thereof so that I wind up with the "regular"
triangulated mesh I mentioned above :) Whatever works
and gives the best fit :)
> If this is the case, then I do not understand why you should
> be concerned about the order of continuity in step 1.
I don't think I am, but I may not have said it correctly...
My primary interest is that the surface fit as well as possible
through the actual data points so when I display the surface
with little mini-spheres to represent the data that the spheres
actually look like they belong to the surface rather than
hovering somewhere near about... does that make sense ?
> In
> fact, it is not clear why you should be concerned about
> the original surface containing the random points. Somehow
> I get the impression you believe that by having the original
> surface contain the data points that you have a "better fit"
> than with a surface that does not contain the points. If so,
> how are you measuring the goodness of fit?
Well, I'm building the surface from sampled data.... my input
data is random... there isn't an original surface :P But I can bet
I'm totally not understanding what you're saying here :)
Well, thanks Dave for your feedback.
From what you've said in this post, I'm guessing I should put my
books aside for a little while and play with your qdintp2.* code,
... if I can figure out what to do with the "function" inputs :)
Any hints as to what I should do there ?
Many thanks :)
Mike
ras...@bayarea.net
>mag<
Thanks to all for your fantastic insight and patience :)
They say you can't teach an old dog new tricks, but
I'm determined to learn this stuff... I've been pushing
pixels since 1982, but all in the low level pipe optimizations....
a "span head" ... But a couple of years ago, one of
the startups I was in (among quite a number) went
belly up and I'd had my fill of spans... it became more
and more like "spam".. :) had that nasty ol' taste to it...
The only "kewl" things were regarding steroscopy or
game physics and GUI's, and there's only so much one
could do at home with stereoscopy when I decided to
go it alone, and it takes a lot to go head-to-head with
a game company these days... I had always been a moderate
net "hacker" (old school version, i.e. would try any kind
of net code experiments, wrote a tcp stack, a kernel here
and there etc) but I needed a project doable for one man...
It was all I could afford (actually I couldn't but tried anyway)
So, I diverted and took a trip around the block to see what
I was missing... the great findings in AI, fuzzy logic, neural nets,
genetic algorithms... and then this contract came along... they
had this ancient drawing package attached to a really kewl
machine and they initially just wanted some cleanup work and
a few new GUI gadgets, but when I looked at their surface code,
it was terrible.. so, on the side I'm trying to figure out the
surface fit scene well enough to at least go back and give them
something that is "definable" with some tangible characteristics...
That's what pulled me to start this thread... They aren't asking
for this work, and may not even want it if I do it, but to me
"it's there" and a great little challenge that isn't pushing pixels
by the span or writing yet another curvey GUI controlly thingy
that will get lost in the dust...
:)
Many thanks to all... sorry if that history was boring, but
I figure some of you deserved to know where I was coming
from with your hard earned attention to trying to help me.
Dave Eberly
> So how do I know the first derivative of a function I haven't
> seen yet ?
It is even "worse" than that. You have a set of points
(x_i,y_i,z_i) that are sampled from an unknown function f
with z = f(x,y). You know that z_i = f(x_i,y_i) but you do not
know the function let alone its derivatives. As I understand
your problem statement, you want
(1) a function z = g(x,y) so that g(x_i,y_i) = z_i,
(2) g is a C^2 function (second derivatives are continuous)
(3) g is as "close" to f as possible.
Additionally you might want
(4) g has "local control"
meaning that if you perturb a data point, the function g
only changes locally about that point. Far away from the
perturbed point, g remains the same.
The qdintp2 code gives you (1) and (4) but not (2). My
documentation does not discuss (3).
The tpspline code gives you (1) and (2) but not (4). My
documentation does discuss (3). Thin plate splines are
set up to minimize bending energy ("closeness" measured
by an integral norm on the second derivatives).
> > * tpspline.{h,cpp}: Uses thin plate splines for fitting points
> > of the form (x,y,f(x,y)). The fit has C-infinity smoothness.
>
> Again, I'm a little stuck on the notion of supplying an input
> function :P
The algorithm does not require you to supply an input function.
You feed it the input data points. I mentioned that the points
are of the form (x,y,f(x,y)) solely to stress that you are thinking
of the z-component as a function of x and y, not that you have
to provide some function f to the code.
> But, I understood bi-cubic to denote a cubic polynomial or rational
> equation over 2 orthogonal curve sets... I see "cubic curves" vs
> "bicubic patches" and have not yet seen a definition of "bicubic",
bicubic = bi + cubic
'bi' refers to 2 (as in 2 independent variables)
'cubic' refers to the degree 3
A bicubic polynomial has 2 independent variables. The largest
degree term for either variable is 3. Note that if you were to
restrict the bicubic polynomial P(s,t) to a line s = a*t+b, you
will get a degree 6 polynomial in t.
Both triangle and rectangle patches can be represented as
bicubic polynomials. Bezier triangles can be used to build
a triangle patch in terms of barycentric coordinates (u,v,w)
where u+v+w=1. The resulting bicubic polynomial is in
terms of u and v (w = 1-u-v, so w is not an independent
variable).
Most folks tend to use 'bicubic patch' to mean 'bicubic
rectangle patch', but as you can see that leads to confusion
since earlier I indicated you can get 'bicubic' with triangle
or rectangle patches.
Check out "Curves and Surfaces for Computer Aided Geometric
Design" by Gerald Farin for a more detailed discussion.
A couple of your statements:
% statement 1
> yes, I'm trying to figure out how to get a random set of
> data points somehow converted eventually to an n x m
> mesh such that the mesh "fits" the data, hopefully
> with C2 continuity.
% statement 2 (reply to my "Are you looking for a mesh
% of triangle patches, each patch being cubic?)
>
> That sounds right... what I've had my radar locked onto
> was a "local" cubic fit, with the final result being the
> triangulated mesh, where the points of the triangles are
> coincident with a regularly spaced 2D grid as seen
> "from above"
Note that both of these indicate a slightly different formulation
than what I mentioned for z=g(x,y). Your initial points are not
necessarily on a rectangular grid, but you appear to want
a rectangular grid to build an interpolating function. I believe
that you have to select between the following two alternatives.
Alternative 1. Construct g according to (1) and (2) and call
it a day. You have a C^2 function that exactly contains the data.
You can use z=g(x,y) to evaluate at new locations whether or
not they are on a rectangular grid.
Alternative 2. Construct some g according to (1) without
regard to order of continuity. Resample that function on
a rectangular grid. Now fit the resampled points with a
mesh of bicubic rectangle patches. The corresponding
function z = h(x,y) is C^2.
> What I think I understand of the material I've read so far on
> splines etc, a cubic patch would (by definition ?) provide better
> continuity than a quadratic patch, or am I misunderstaning the
> concept ?
A cubic provides "higher order continuity" than a quadratic.
By using the term "better continuity" implies that you require
a certain order of continuity. Perhaps you should revisit
this assumption. Does your final goal really require C^2?
> > My impression from the thread is that you are trying to
> > 1. build a surface that contains the random points
> > 2. resample to get a rectangular grid of points
> > 3. fit the rectangular grid to get a bicubic patch surface
>
> Or any order thereof so that I wind up with the "regular"
> triangulated mesh I mentioned above :) Whatever works
> and gives the best fit :)
>
> > If this is the case, then I do not understand why you should
> > be concerned about the order of continuity in step 1.
>
> I don't think I am, but I may not have said it correctly...
> My primary interest is that the surface fit as well as possible
> through the actual data points so when I display the surface
> with little mini-spheres to represent the data that the spheres
> actually look like they belong to the surface rather than
> hovering somewhere near about... does that make sense ?
It sounds to me like either qdintp2 or tpspline will satisfy your
goals. Your last paragraph indicates that the resampling on
a rectangular grid is to produce a set of points that appear
to be on the same "surface" that the original points are on.
Either of the two suggested algorithms should produce this
visual effect.
>mag<
Dave Eberly wrote:
> ...<lots of very kewl things>....
Well, I have learned a lot from this thread and this last
post from Dave seals the immediate tasks I have to try :)
I appreciate all of your replies very much.
Thanks,
Mike
ras...@bayarea.net
Muhammad Ali
> In article <38BDFCEF...@bayarea.net>,
> ras...@bayarea.net wrote:
> > Hi,
> > I'm looking for C/C++ bivariate surface fit code. I've posted a
> similar
> >
> > request previously but the lack of responses lead me to believe I had
> > not well-stated my request, or perhaps it is indeed something
> difficult
> > to find, so I apologize if this is redundant for some, and many thanks
> > to those who are able to reply, esp. via email as my ISP removes
> > messages
> > every day or so.
> >
> > Puzzle:
> > Given a set of data points randomly positioned over a 2D area and
> > of random value (or height, as you prefer)....
> > ** How do I generate a regular mesh (say, 50x50 etc) of polygons
> > or triangles such that the surface is
> > A) a bicubic interpolation *through* the points with the
> > "tightest"
> > curve, and/or "least" warped by the choice of coefficients
> > * again, I'm looking for that regular 50x50 mesh... I did
> > get
> > a couple of emails from people prescribing tesselators,
> > surface
> > subdivision etc, but that isn't what I had requested or
> > what I'm
> > looking for since, unfortunately, it doesn't fit the
> data
> > as prescribed :)
> > and/or
> > B) well approximated as near the points as possible with "good"
> > curvature, and perhaps by selection of some limit of error,
> > eg. possibly using least squares and I could only guess as
> > what kind of rational spline to use, whether bicubic or
> > whatever..
> triangulated network would be less costly in both time and memory.
> I have some triangulation and linear interpolation code in C that
> you might be interested in.
> you to check out Dave Eberly's site. Paul Bourke has some helpful
> source code too, I think:
> http://www.swin.edu.au/astronomy/pbourke/
hi i am really interested in knowing more about your code as i am trying to do a local bicubic surface interpolation of (x,y,z) data set that would give me the resulting control points as well as the knot vecotors in the u and v direction. Hope you could help