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@lists.osgeo.org

http://lists.osgeo.org/mailman/listinfo/grass-dev

--

Dr. Benjamin Ducke, M.A.

{*} Geospatial Consultant

{*} GIS Developer

benducke@fastmail.fm