[GRASS-dev] GSOC wk4 checkin Horizon Based Voxel Interpolation

[Benjamin Ducke]

Hi Tim,

On 07/13/2013 01:37 AM, Tim Bailey wrote:
> Tim Bailey July 12, 2013
>
> GSOC week 4 check in
>
> Horizon based voxel interpolation
>
> This was an interesting week. After thorough reconsideration of the
> foundational documents for this project as well as reviewing about a
> dozen other three d interpolation projects, I have settled on a
> methodology that I believe resolves many of the outstanding issues. The
> vertical to horizontal anisotropy problems are a rabbit hole that I do
> not need to resolve during this process. Rather than get into the weeds
> of benefits and liabilities of different interpolation operators, I
> think that the best way to move forward is by deterministic flood
> population of voxels. Clearly different disciplines have different
> methods and needs for spatial approximations. In my estimation, the end
> users decisions about what the appropriate geostatistical or logical
> operator to apply should take place between generation of conforming
> horizons, and the operation that generates the r3 map.
>
> The most important point of this is the population of the r3 voxels will
> proceed through a deterministic operator, which can be accomplished with
> existing tools. The various approximations such as geostatistical
> operations can be accomplished before the r3 map is initiated. The r3
> map can essentially be flood filled from points by selecting all the
> points in a tier and applying v.voronoi, followed by v.to.rast. A data
> artifact from this process is an arbitrary step pattern caused by sparse
> data. If you want to smooth transitions you can introduce intermediate
> horizon approximations based on the logical model of the operators
> choosing.

The Voronoi diagram seems like a good choice, as it is a simple
method that makes the results of the interpolation easy to
predict and interpret. However, instead of using v.voronoi
and v.to.rast, you should consider creating your own interpolator
(although you can use v.voronoi/v.to.rast to build a simple
prototype at this stage).

A basic implementation of Voronoi in voxel space is very simple:
You just measure the 3D distances between the empty voxel and all
3D sample points. You then assign the voxel the value of the closest
point. This avoids the overhead involved in using v.voronoi/v.rast
and should be a lot more CPU and disk space efficient.

 From there, you can go on to solve the other issues by using
a _weighted_ Voronoi approach:

1. Simply exaggerate the Z data to account for horizontal anisotropy:
if distances along the Z axis get larger, than your Voronoi cells will
"flatten". Experiment to find a good default scaling for "Z" and add
an option for the user to set the Z scaling coefficient manually.

2. Similarly: You can add increased weight to horizon categories,
depending on the number of cores in which they appear ("evidence
weight").This will suppress the influence of horizons that appear
in only a few samples.

3. You could also define a maximum distance from the closest
sample point, beyond which a voxel will simple remain NULL.
This will reflect reality better, as it will take care of the
edge effect of the Voronoi diagram (cells towards the edges
of the study region tend to become really large).

All of these are really simple to implement, but it might take
some experimenting and good design of parameters to come up
with good default weights and a weighting scheme of distance
vs. evidence weight.

Your basic Voronoi formula could then be:

I=D^w*E

Where

I = "influence": this is calculated for each empty voxel
and each sample point and the empty voxel is assigned the
value of the point with the highest score for I
D = distance from sample point
w = distance weight (by default 2.0)
E = evidence weight

The distance is calculated as:

D=sqrt((x1-x2)^2+(y1-y2)^2+((z1-z2)*m)^2)
With
m=Z exaggeration factor

In the simple case, the evidence weight would be:

E=sum(sample points with same horizon ID)


You will certainly want to play with different
settings for w,E and m -- but that should be fun!

Another question is what to do with voxels that
are not empty, i.e. that have sample points falling
into them. These could either be preserved, i.e.
you convert them to voxels first and then perform the
interpolation only for non-empty voxels; or your could
ignore them during the interpolation (and then use them
for validating the quality of the interpolation result).
The latter approach has the advantage that you do not
need to think about what to do when several points fall
into the same voxel.

Best,

Ben



>
>
> I have a proof of concept dataset which I manually shepherded through
> what is at this point a ten step transformation between a conforming
> horizon descriptions, and an r3 map. Next week I plan on continuing
> working on the generic wrapper that manages the transformation of
> horizon intervals to voxels. Also I will get my project wiki page up on
> very soon.
>
> Have a good weekend
>
> Tim
>
>
>
>
> _______________________________________________
> grass-dev mailing list
> grass-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-dev
>



-- 
Dr. Benjamin Ducke, M.A.
{*} Geospatial Consultant
{*} GIS Developer

   benducke at fastmail.fm