🚀 go-pugleaf

------- Start of summary -------
This summarizes (to some extent) the responses that I received 
on my posting "Subject: image reconstruction/irregular grid." 

I still intend to followup with some of you, but it's taking a while
to assimulate all of this information.  If anyone has
anything further to add or any additions/corrections to
the list of references, let me know (I could not find
some of the references in either Purdue's or Texas A&M's
online library system).

CONTENTS
 1. Thanks
 2. Original post
 3. Edited replies
 4. BibTeX summary of references recommended

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1. Thanks

First of all, let me thank all of the people who responded (either
by email or a followup):
 dan@research.nj.nec.com (Dan Ruderman)
 yking@cs.ubc.ca (Yossarian Yggy King)
 PACKER@amarna.gsfc.nasa.gov (Charles Packer)
 sjreeves@eng.auburn.edu (Stan Reeves)
 sloan@cis.uab.edu (Kenneth Sloan)
 rak@crosfield.co.uk (Richard Kirk) 
 watson@maths.uwa.oz.au (David Watson)
 nicholls@stats.ox.ac.uk (Geoffrey Nicholls)
 A4422DAJ@vm.univie.ac.at (Christian Cenker)
 butler@cluster.gps.caltech.edu (Bryan Butler)
 cice@cacs.usl.edu (C. H. Chu)
 ph@HOSTESS.GRAPHICS.CS.CMU.EDU (Paul Heckbert)
 a8131dan@helios.edvz.univie.ac.at (Hans G. Feichtinger)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2. Original Post

I'm looking for some feedback on a problem that I am working on.  It
is basically an interpolation problem in which I would like to use
image reconstruction techniques (or anything similar that is relatively
fast and non-interactive, e.g., thin plate splines). The catch is
that I do not have ground truth or an underlying model.

I have data (a soil property) on a non-uniform grid for an area of a
few hectares.  In one direction, the samples are spaced approx. 5-8
meters apart, but in the other, they are spaced approx 20-30
centimeters apart.  [I'm using a sensor and a slow-moving vehicle to
make these measurements. Since the vehicle does not/cannot always travel
in straight lines, I do not have a regular grid.] As far as I know, 
no one has ever sampled soil properties this densely.

Since I have no ground truth, my approach to testing different
interpolation methods is to find some (spatial?) model of the soil
property so that I can
 1. create a pseudo-realistic example surface,
 2. resample (simulating my current sampling process),
 3. use whatever interpolation techniques that I want to compare,
 4. then compare the errors.

If my example surfaces are the best that I can create, then the "best"
technique should be the one that I use for my real data.

My first thought regarding a model was to find published semivariogram
constants/models for this soil property (hence pseudo-realistic) and
use kriging to create my example surface.  How does this sound?

Since the sampling is so dense, I bawk a little at using the
semivariogram/kriging approach as my underlying model, but this is the
only way that I can think that I may be able to use realistic surfaces
in my testing (I assume that there are a few published reports of
semivariogram fitting for the variable that I am interested in).

Why not just use kriging as my interpolation method? Well, I need
something that can be fully automated.  Many folks more knowledgeable
than myself have insisted that semivariogram fitting simply
cannot be automated.

So, does this approach sound good or bad? Could it be improved upon?

FOR THE IMAGE PROCESSING FOLKS:

Does anyone have any references that present an overview of image
reconstruction from irregular sampling? The only ones that I have
found make an assumption of band-limited images. I suppose what I'm
looking for are methods which do not make very strict assumptions like
this.

To be clear, the space that I'm sampling from is not actually an image
(and not rectangular), but I am hoping that if I pad with zeros or
extend my boundaries, an image reconstruction technique may work (I
hope that there will be no significant edge effects).  Pixels/points
are real-valued and samples are from the spatial domain.

Image reconstruction is an area in which I have practically no
experience, so I would be most grateful for any pointers. Email and I
will summarize (to s.i.p if necessary). If anyone wants the refs that
I mentioned above, I have temporarily put them up for anonymous ftp from
pasture.ecn.purdue.edu:/pub/mccauley/misc/recons.irr-samples.band-limited.

My apologies for this being sort of lengthy, but I'm afraid that more
bandwidth would be wasted if I did not make myself clear the first
time. This is for my PhD research and I haven't gotten to know too
many people locally that I could bounce these ideas off of.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

3. Edited replies:
------------------------------
A fairly standard approach to surface reconstruction in the computer
vision literature is to use regularization, where one tries to select
the surface that "best" fits the data while maintaining some degree
of smoothness/coherence/fit to your expected model. Typically this
is laid out mathematically as a problem of minimizing:

	F[z*] = Omega[z*] + Lambda*SUM(Rho(Mi[z*],Ci))

where
	Omega measures the "smoothness" of the surface z* (fit to model)
	Mi[z*] represents a point on the estimated surface z*
	Ci is a constraint from your data
	Rho(Mi[z*],Ci) is a measure of how well point i on the surface
		matches its corresponding constraint. Thus the sum over
		all the constraints is a measure of fit to data.
	Lambda is a parameter that controls the tradeoff between
		smoothness and fit to data.

A good explanation and demonstration of this technique is in W.E.L. Grimson,
>From Images to Surfaces, MIT Press 1981. Grimson's function minimizes the
quadratic variation of the surface while fitting range data, thus:
	Omega[z*] = INT(INT(Zxx^2 + 2Zxy^2 + Zyy^2)dx)dy
	Rho(Mi[z*],Ci) = (Zi - Ci)^2	{Euclidean distance}

This functional is probably not the best for your purposes, but the
explanation may be helpful from a technique standpoint. For a more
up to date paper that uses a similar approach (and which may have some
good references for you), have a look at Stephenson/Delp, "Viewpoint
invariant recovery of visual surfaces," IEEE Trans. on Patt. Anal.
and Mach. Int., 14(9), 1992.

------------------------------

Is the Landsat image resampling methodology of any relevance here?
We resample from irregular grid to regular output grid. I'm not
a mathematician, but I've spent a non-trivial part of my programming
career on this problem.

[beats me --jdm]

------------------------------

[posted as a followup to s.i.p] 
|> If you want to get a unique mapping from your samples to a continuous
|> surface, then you must make some assumptions that are as restrictive as
|> bandlimitedness.  Bandlimitedness may or may not be the best assumption
|> in your case, but it at least guarantees a one-to-one correspondence
|> between the samples and the continuous surface.  There are certainly
|> other options that would use different basis functions and still 
|> guarantee uniqueness, but I don't think the theory would be nearly as
|> well-developed.

------------------------------

My general impression of your article is that you are pushing very, very
hard on the easy part of the problem.

If you are sampling very densely, compared to previous work, then
interpolation is not your problem.  First, show that really
simple-minded interpolation schemes DO NOT work.  If you can't show, for
example, that piecewise linear approximation of the surface by
triangular facets is dead-wrong, then it is the method of choice.
[triangulate the convex hull of your data poins with the Delaunay
 triangulation; assign the measurement at each site as the z-coordinate
of each vertex in the triangulation, and use linear interpolation within
each triangular domain.  Avoid at all costs extrapolating the surface
beyond the convex hull. Done.

Well, perhaps not done.  Depending on what you plane to *do* with this
surface, you may want to now re-sample the facetted surface on a regular
grid.  Feed these grid points to your favorite surface fitter.

Or...look at the CAD literature for techniques involving higher-order
surface patches.  You may need to estimate derivatives of your function
at the vertices.

Now for the hard parts - how much do you believe your measurements?
That is - *must* the final surface *interpolate* the data points?  Why?
Are you sure?  If not - then you want an approximating scheme, not an
interpolating scheme.  Here's the simplest one I know:  do the same
triangulation as above, but then "relax" the mesh by perturbing each
data point towards the average of its neighbors (neighbors are those
data points joined by an edge in the triangulation).

And...what properties do you want the final surface to have?  Does it
have to be smooth?  How smooth?  Why?

I suppose these numbers were in your article, but I've forgotten -
please remind me:  how many data points do you expect to have?  10, 100,
1000, 1000000?  And, what sort of anisotropies are there in the sampling
scheme?  Are the data points densely arranged along lines (curves) which
are themselves sparsely distributed over the domain?  In that case, it
makes sense to solve the 1D problems separately, and then conbine.

Feel free to poke me if any of the above interests you.

[I probably will --jdm]
------------------------------

We made a very quick image pyramid filter that can generate a whole image
from sparse data.  It produces a smooth surface that should pass through the
points (they are assumed to lie on pixel positions, so you may have to
round your data a bit.  It is not as fancy as the method you suggest, and 
there are some artefacts - for example a ring with a dot at the centre is not
filled with a circularly symmetric function - but otherwise the results are
smooth at the level you might expect.  I'm not sure if it is suited to your
data - I suspect we want many more points to our images than you wnt on your
data, but here it is anyway...

(1). Constructing the pyramid....

	Make a low resolution (LR) image half the size of the high resolution 
	(HR) you started with, where the LR pixel value is the average of all
	the defined values in the corresponding 2*2 pixel region on the HR.  If
	none of the HR pixels are defined, then the LR pixel is undefined.  You
	can reserve a value or use an extra bit plane to say whether a pixel
	is defined or not.

	Repeat the above with boundary conditions until all the pixels are 
	defined.  This will always be the case if you go on until the LR is 
	only one pixel provided something was defined, and that avoids doing
	a lot of tests on the biggest HR levels.

(2). Reconstructing the image...

	Starting with the smallest LR image of the previous set of images and
	the HR image it was derived from, where the HR pixel is undefined we
	replace it with the interpolated value from the surrounding LR pixels.
	Though the interpolation is linear, successive interpolations at 1/4 
	and 3/4 intervals give a parabolic spline but the actual values may
	depend on the register between the defined values and the LR pixel 
	positions.

	The result will be an HR image where all the pixels are defined.
	Repeat the process until you are back to the original resolution.

We use it to generate artificial shadows, shading, filling in irregular gaps
in images, and so on.

------------------------------

Two comments 
             - first could you think of it as a regular lattice with
missing data ? There's certainly work done on this.

[no --jdm]

             - second, from image processing, if you model with a
Markov Random Field on an anisotropic lattice, presumably using a
simple clustering model 

- - represented by a Gibbs potential like

\sum_{i\in S_p}\sum_{j\in n_i} (a_i-a_j)^2 + 
	\sum_{i\in S_d} (a_i-g_i)^2/\sigma^2

where a_i are the estimated values (poss. vector) g_i are data,
and S_p is the set of pixels, which may be a superset of the set
S_d of pixels at which you have data -

and have anisotropic or completely inhomogeneous Gibbs parameters -
then you have alot of scope to do modeling, and are in
well charted statistical teritory. There are some good applied
refernces on this kind of thing. You might click with

[references deleted. see end of this posting. --jdm]

The point is there's nothing inherently regular-lattice about MRF modeling.

------------------------------

[posted as a followup to s.i.p. --jdm] 
|> A group at the univerity of vienna/austria/europe is working on
|> non-uniform sampling and reconstruction of signals, most
|> ly mathematical theory, but they are doing a lot in
|> implementing their algorithms on suns.
|> 
|> if you re interested in their work please contact
|> hans g. feichtinger under a8131dan@vm.univie.ac.at

------------------------------

you might want to check out:

R. J. Renka, Multivariate Interpolation of Large Sets of Scattered
Data, ACM Transactions on Mathematical Software, 14, 1988, 139-148

and you can get the code by ftp'ing to research.att.com.
i've used this to do conversions from planetary images (full disk)
to lat-long projections.  it's pretty slow in 3 dimensions (which
is what i use it for), but i did some tests on the 2-D version,
and it seemed to go pretty fast.

------------------------------

Here are some references on that subject:
Marvasti's book (if you can find a copy at a library or order one)
seems to be a very thorough survey of this topic; though not very
implementation-oriented.

[I couldn't not find this, but it sounds like an interesting
 work. Anyone have a phone number for ordering? --jdm]

You'll also find related work in the fields of geometric modeling and
computer graphics under topic headings like "scattered data interpolation".

I believe that the approach thay many in the geometric modeling
community would take to your problem is to triangulate the set of points
(Delaunay triangulation does a good job) and then use Bezier triangles
with certain continuity constraints to interpolate.
Farin's book should help you get started if you choose that route.

------------------------------

[posted as a followup to s.i.p. --jdm] 

|> Scattered data interpolation can be done by a wide variety of methods
|> For a summary of such methods, ftp
|> cod.nosc.mil  /pub/grid.ps
|>               /pub/grid.tex
|>               /pub/grid.ascii
|> However, because your data is anisotropically scattered it is likely
|> that the only method you will find to be satisfactory is natural
|> neighbor interpolation because this methods uses local subsets of
|> data which are choosen in a manner that compensates for the
|> anisotropic density of the data.  If you, or anybody else, would like
|> a demonstration of this method, send me your most nasty data set and
|> specify the desired regular grid size.  I will grid the data, either
|> on a rectangular or triangular grid as you wish, and then you may
|> display the grid with any conventional contouring package.  If the
|> grid size that you choose is dense enough, the contouring package
|> will not alter the representative surface defined by my grid.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4. BibTeX summary of references recommended

@Book{ marks91,
  author = 	 "Robert J. Marks",
  title = 	 "Introduction to {S}hannon Sampling and Interpolation Theory",
  publisher = 	 "Springer--Verlag",
  year = 	 1991,
  address = 	 "New York, NY",
  purdue = 	 "Engineering Library    511.42 M342i 1991",
  annote = 	 "ISBN 0387973915 (U.S.), 3540973915"
}

@Article{ jerri77,
  author = 	 "Abdul Jerri",
  title = 	 "???",
  journal = 	 "IEEE Proceedings?",
  year = 	 "1977"
}

@Book{ grimson81,
  author = 	 "William Eric Leifur Grimson",
  title = 	 "From Images to Surfaces : A Computational Study of
		  the Human Early Visual System",
  publisher = 	 "MIT Press",
  year = 	 "c1981",
  series = 	 "Artificial intelligence",
  address = 	 "Cambridge, Mass."
}

@InProceedings{,
  author = 	 "Robert L. Stevenson and Edward J. Delp",
  title = 	 "Viewpoint invariant recovery of visual surfaces from
		  sparse data",
  pages = 	 "309--312",
  booktitle = "Proceedings 3rd International Conference on Computer Vision",
  year = 	 "1990",
  organization = "IEEE Computer Society",
  publisher = "IEEE",
  address = 	 "IEEE Service Center, Piscataway, NJ, USA ",
  month = 	 dec,
  annote =       "(IEEE cat n 90CH2934-8)"
}

@Article{besag91,
  author = 	 "J. Besag, and J. York and A. Mollie",
  title = 	 "Bayesian Image Restoration with Two Applications in
		  Spatial Statistics",
  journal = 	 "Ann. Inst. Statist. Math.",
  year = 	 "1991",
  volume = 	 "43",
  pages = 	 "1--59"
}

@Misc{ geman1,
  author = 	 "S. Geman and D. Geman",
  title = 	 "Stochastic Relaxation, Gibbs Distributions, and the
		  Bayesian Restoration of Images",
  howpublished = "??",
  year = 	 "??"
}

@Book{ geman90,
  author = 	 "D. Geman",
  title = 	 "Random Fields and Inverse Problems in Imaging",
  publisher = 	 "Springer-Verlag",
  year = 	 "1990",
  volume = 	 "1427",
  series = 	 "Lecture Notes in Mathematics"
}

@Article{renka88,
  author = 	 "Robert J. Renka",
  title = 	 "Multivariate Interpolation of Large Sets of
		  Scattered Data",
  journal = 	 "ACM Transactions on Mathematical Software",
  year = 	 "1988",
  volume = 	 "14",
  pages = 	 "139--148",
  month = 	 jun
}

@Article{ chen87,
  author = 	 "David Shi Chen and Jan P. Allebach",
  title = 	 "Analysis of Error in Reconstruction of
		  Two-Dimensional Signals from Irregularly Spaced Samples",
  journal = 	 "IEEE Trans. Acoustics, Speech, and Signal Processing",
  year = 	 "1987",
  volume = 	 "ASSP-35",
  number = 	 "2",
  pages = 	 "",
  month = 	 feb,
  note = 	 "(this assume bandlimited images)"
}

@InProceedings{ mitchell87,
  author = 	 "Don P. Mitchell",
  title = 	 "Generating Antialiased Images at Low Sampling Densities",
  volume = 	 "21",
  number = 	 "4",
  pages = 	 "65--72",
  booktitle = "Computer Graphics (SIGGRAPH '87 Proceedings)",
  year = 	 "1987",
  month = 	 jul
}

@InProceedings{ shlomot89,
  author = 	 "E. Shlomot Y. Y. Zeevi",
  title = 	 "A Nonuniform Sampling and Representation Scheme for
		  Images Which are Not Band--limited",
  booktitle = "Proc. IEEE Conf.",
  year = 	 "1989",
  address = 	 "Tel Aviv"
}

@Book{ marvasti87,
  author = 	 "Farokh A. Marvasti",
  title = 	 "A Unified Approach to Zero-Crossings and Nonuniform
		  Sampling of Single and Multidimensional Signals and Systems",
  publisher = 	 "Nonuniform Pub.",
  year = 	 "1987",
  address = 	 "PO Box 1505, Oak Park IL 60304"
}

@InBook{schumaker79,
  author = 	 "L. L. Schumaker",
  title = 	 "Approximation Theory II",
  chapter = 	 "Fitting Surfaces to Scattered Data",
  publisher = "Academic Press",
  year = 	 "1979",
  pages = 	 "203--268"
}

@Article{ franke82,
  author = 	 "R. Franke",
  title = 	 "Scattered Data Interpolation: Tests of Some Methods",
  journal = 	 "Mathematics of Computation",
  year = 	 "1982",
  volume = 	 "38",
  pages = 	 "181--200",
  month = 	 jan
}

@Book{farin88,
  author = 	 "Gerald E. Farin",
  title = 	 "Curves and Surfaces for Computer Aided Geometric
		  Design : A Practical Guide",
  publisher = 	 "Academic Press",
  year = 	 "1990",
  address = 	 "Boston",
  edition = 	 "2nd"
}

@Book{,
  author = 	 "James D. Foley and Andries van Dam and Feiner and Hughes",
  title = 	 "Computer Graphics : Principles And Practice",
  publisher = 	 "Addison-Wesley",
  year = 	 "1990",
  series = 	 "Addison-Wesley Systems Programming Series",
  address = 	 "Reading, Mass.",
  edition = 	 "2nd"
}
------- End of summary -------

-- 
James Darrell McCauley                Dept of Ag Engineering, Purdue Univ
internet: mccauley@ecn.purdue.edu     West Lafayette, Indiana 47907-1146, USA
bitnet: mccauley%ecn@purccvm          UUCP: pur-ee!mccauley