Interactive Seismic Interpretation
with Piecewise Global Energy Minimization
Thomas Höllt∗

Johanna Beyer∗

Fritz Gschwantner†

Philipp Muigg‡

King Abdullah Universityof Science and Technology

King Abdullah Universityof Science and Technology

VRVis Research Center

SimVis GmbH

Helmut Doleisch‡

Gabor Heinemann§Markus Hadwiger∗

SimVis GmbH

Heinemann OiL GmbHKing Abdullah Universityof Science and Technology

A BSTRACT
Increasing demands in world-wide energy consumption and oil depletion of large reservoirs have resulted in the need for exploringsmaller and more complex oil reservoirs. Planning of the reservoir
valorization usually starts with creating a model of the subsurface
structures, including seismic faults and horizons. However, seismic
interpretation and horizon tracing is a difficult and error-prone task,
often resulting in hours of work needing to be manually repeated.
In this paper, we propose a novel, interactive workflow for horizon
interpretation based on well positions, which include additional geological and geophysical data captured by actual drillings. Insteadof interpreting the volume slice-by-slice in 2D, we propose 3D seismic interpretation based on well positions. We introduce a combination of 2D and 3D minimal cost path and minimal cost surface
tracing for extracting horizons with very little user input. By processing the volume based on well positions rather than slice-based,
we are able to create a piecewise optimal horizon surface at interactive rates. We have integrated our system into a visual analysis
platform which supports multiple linked views for fast verification,
exploration and analysis of the extracted horizons. The system iscurrently being evaluated by our collaborating domain experts.
Index Terms:
I.3.8 [Computing Methodologies]: Computer
Graphics—Applications; I.4.6 [Image Processing and Computer
Vision]: Segmentation—Edge and feature detection
1

I NTRODUCTION

Fossil fuels are today’s most important source of energy. According
to the International Energy Outlook 2010 [36] by the U.S. Energy
Information Administration, roughly 58% of the marketed energy
in 2007 were oil and gas products. It is expected that in 2035 oiland gas will still account for more than 50% of the marketed energy. The combination of dwindling oil in place and increasingenergy demand in general results in rising prices and justifies the
exploration of smaller or more complex reservoirs.
Planning of reservoir valorization usually starts with creating amodel of the subsurface structures, the seismic interpretation. A
seismic reflection volume, or seismic cube, is acquired by sending
sound waves, for example triggered by explosions, into the earth
and recording the emerging echoes. At the boundary of two subsurface layers, a seismic horizon, a part of the induced wave proceeds,
while another part is reflected. The reflected waves are recorded
over time by a grid of geophones arranged on the ground, result∗ e-mail:{thomas.hollt|johanna.beyer|markus.hadwiger}@kaust.edu.sa† e-mail:fritz.gschwantner@vrvis.at‡ e-mail:{muigg|doleisch}@simvis.at§ e-mail:gheinemann@heinemannoil.com

ing in a set of 1D traces containing a reflection value for each time
step. The values for a single time step of all geophones arrangedin the grid form a time-slice. The grid itself defines the spatial resolution. In the volume, a set of traces with equal position in thegrid for either the x- or y- coordinates form the crossline and inline
slices, respectively. The traces and time-slices are inscribed along
the z-axis (see Figure 3).
The extraction of seismic horizons is one of the main tasks wheninterpreting the seismic cube to build a model of the subsurface
structures. Horizons are represented in the volume by bands of locally extremal values, whereas most other structures like faults areprimarily defined by their interaction with horizons.
The seismic interpretation of the volume can be used for thetime-depth conversion, a process where the time axis is converted
to actual spatial depth values. Estimating the sonic speed for eachlayer, the layers can then be stretched or compressed using the associated sonic speed resulting in a depth-converted volume.
In most cases, the seismic interpretation is not based solely on
the seismic cube. Usually, physical drillings provide additionaldata, resulting in well logs. A well log is a detailed record of acertain property of the subsurface structures at the drill hole. Logs
can be either geological or geophysical. The former are acquired byvisually or chemically inspecting samples brought to the surface,but also contain live drilling parameters such as the speed at which
the drill bit deepens the borehole (i.e., the rate of penetration), thelatter by lowering a measurement device into the borehole. A comprehensive survey of a drilling usually contains several well logs,
geological as well as geophysical. Typical properties besides therate of penetration are for example weight, porosity or resistivity.
In the following, we use the term well log synonymously for a single log as well as for a complete survey of one drill hole, containing
several logs. The well log data is especially useful when performing time-depth conversion, as it contains actual information about
the depth location of subsurface structures.
Interpreting the seismic cube is a cumbersome, time-consuming
process. The data is dense, hard to visualize in 3D, noisy, and ambiguous. It may happen that after hours of interpretation work it
becomes apparent during time-depth conversion that large parts ofthe volume were misinterpreted due to a single wrong decision inwhat seems to be a branching horizon. In current practice, which ismainly based on manual inline-, and crossline-slice inspection, thisis very tedious, as it requires manual correction of multiple slices.
Additionally, these slices need to be adjusted according to well log
information which usually is only available for very few slices, andmight be located between interpreted slices where no interpretation
has been performed.
In this paper, we propose a novel workflow for seismic horizon
interpretation based on well positions. We triangulate the spatial
positions of all available well logs, which divides the volume into a
set of triangular prisms. Instead of interpreting the volume slice-byslice in 2D, we propose performing full 3D seismic interpretation,
going from prism to prism. The part of a horizon intersecting agiven prism can be extracted by specifying as little as a single seedpoint, using global energy minimization techniques to compute a
horizon surface of minimal “cost”. We combine very fast 2D minimal cost path tracing for determining horizon intersections along
the faces of the prisms resulting from well log triangulation withsubsequent 3D minimal cost surface computation, which is constrained by the 2D contours on the prism faces. By combining the
optimal surfaces from the individual prisms, we receive a piecewise
globally optimized horizon. Extracting the horizon in a piecewisemanner allows for an interactive algorithm which gives the user fullcontrol over adjusting the surface by setting an arbitrary number of
constraints, forcing the surface to pass through specified locations.
We provide a framework with multiple linked 2D and 3D views, toenable easy interaction with the data and verification of the resulting surfaces. We have integrated our system into the SimVis visual
analysis framework [11].
2

R ELATED W ORK

The importance of oil and gas for today’s societies has resulted in aconsiderable amount of previous work in the area of seismic interpretation and visualization, as well as commercial software packages, such as HydroVR [20], Petrel [34] or Avizo Earth [37].
One important area is seismic interpretation and seismic horizon
extraction. Pepper and Bejarano [28] give a good overview overthese seismic interpretation techniques. A lot of work has focusedon fully automatic approaches that require the definition of someparameters, but afterwards work without user interaction. Keskes
et al. [17] and Lavest and Chipot [19] show an abstract outline ofsuch algorithms for fully automated 3D seismic horizon extraction
and surface mesh generation. Tu et al. [35] present an automatic approach for extracting 3D horizons, but based on grouping 2D traces.
Faraklioti and Petrou [14] and later Blinov and Petrou [3] present
real 3D surface reconstruction using connected component analysis
parameterized by the local waveform and layer direction. However,processing the complete volume at once is a lengthy process, thus
parameter tuning is inconvenient. Additionally, optimal parameters
are not necessarily equal for all features in the volume.
Castanie et al. [9] propose a semi-automatic approach. Horizons
are traced one by one from a user-defined seed point. Interactivity in these kind of methods is limited as the horizon extraction iscostly, resulting in long waiting times between seeding.
Patel et al. [26] present a technique for quick illustrative rendering of seismic slices for interpretation. They use transfer functions
based on precomputed horizon properties. Their method, however,
only works on 2D slices. In a subsequent publication [27] they extended their illustration technique for rendering of 3D volumes tobe used for knowledge-assisted visualization of seismic data. Thismethod still relies on slice-based precomputation.
Borgos et al. [4], as well as Patel et al. [25], propose interactive workflows for the surface generation based on an automatic
preprocessing step. Both extract surface patches from the volume
in a preprocessing step by extrema classification and growing, respectively. This step lasts for several hours. In a second, interactive
step the user then builds the horizon surface using the pre-computed
patches. Where in the work of Borgos et al. the surface patches
function as simple building blocks, Patel et al.’s approach is moresophisticated. The extracted patches are subdivided and stored hierarchically. The user starts the surface generation with a single seed
patch and by climbing up the hierarchy adds connected patches.
The authors themselves target their approach to quick extraction ofhorizon sketches which then need to be refined in a second step.
Lampe et al. [18] present a technique to deform and render volumes along curves. They show two applications, one of which isthe visualization of seismic data along well logs.
We propose an interactive seismic interpretation workflow exploiting global energy minimization inspired by segmentation tech-niques employed in computer vision. Our method does not require
any pre-computation and allows modification of computed surface
patches by adding constraints. By partitioning the volume into
prisms, the time waiting for global optimizations is minimized. By
applying the surface extraction on smaller chunks (i.e, prisms) instead of the complete volume at once, we maintain an interactive
workflow for extracting piecewise global optimal horizon surfaces.
Typical image segmentation systems based on energy minimization can be divided into two classes based on the types of constraints they provide. Approaches based on graph cuts [5, 7, 8]
work by dividing the image in fore- and background elements. The
boundary is then defined implicitly between neighboring fore- andbackground image elements. Consequentially, constraints can be
set by defining specific elements as fore- or background. The second type of energy minimization algorithms for segmentation arealgorithms like Live Wire/Intelligent Scissors [22, 24], which trace
a line through a set of given points (the constraints) which is supposed to bound the feature which is to be extracted. Here constraints are explicit, giving the user full control over the boundary.
Depending on the application, both types of constraints have theiradvantages. For seismic interpretation it is crucial to be able to directly interact with the boundary (representing the seismic horizon),making explicit constraints a much better choice than implicit ones.
While graph cuts adapt naturally to 3D [2, 6], or higher dimensions in general (the cut in an n-dimensional graph is always of dimension n-1), the application of Live Wire/Intelligent Scissors will
always produce a line. These lines can be used as a boundary for 2D
objects, but not 3D objects, which are bounded by a surface. Poonet al. [29] exploit this property to track vessels in 3D. There are
a number of publications dealing with the extension of Live Wire
to 3D to create minimal cost surfaces. Falcão and Udupa [13, 12]
propose automatic tracing in 2D using a user-defined trace on orthogonal slices as constraints. Their approach, however, requiresconsiderable user supervision, especially when trying to segment
objects with a topology different from a sphere. Schenk et al. [33]
introduce a combination of Live-Wire and shape-based interpolation. Poon et al. [30, 31] present a technique similar to Falcão et
al. [12, 13] which can handle features with arbitrary topology, by
defining inner and outer contours. However, all these methods are
based on creating some kind of network of lines computed using theclassical 2D Live Wire algorithm, which has several drawbacks; A
set of minimal cost paths in a 3D image does not necessarily assemble a minimum cost surface. Another major problem is that theplacement of the key paths is crucial for the quality of the resulting surface, especially if the topology changes between slices asmight be the case for concave objects. Additionally, from a user
perspective it can be much harder to edit a network of minimal cost
paths. For example, a constraint set in one slice will not necessarily
propagate into the next slices, etc.
The only method we know of to find minimum cost surfaces in
3D with explicit constraints like in the Live Wire approach was proposed recently by Grady [15]. This work formulates the minimal
cost surface problem as an extension of the original 2D minimal
cost path problem and shows that it can be solved with linear programming. In later work, Grady [16] presents an optimized solution for regular 6-connected lattices using minimum-cost circulation network flows. This approach results in the real minimal costsurface and is also independent of topology. The algorithm requires
a single arbitrarily shaped closed contour as input and allows discretionary closed contours as constraints.
3

W ORKFLOWAs input to our system we assume a seismic volume, where thex- and y-dimensions correspond to actual spatial coordinates, andthe z dimension corresponds to the travel time of reflected seismic
waves acquired by geophysics during reflection seismology.

cb

3.1 Current Practice
The workflow of current seismic interpretation systems is typically
based on the manual or semi-automatic interpretation of every n-th
inline- and crossline-slice of a seismic volume, where n can range
from skipping several slices each step to processing each slice individually, one by one. Generally, in these systems, multiple horizons
are extracted for each slice. After all desired slices have been interpreted (i.e., horizons have been extracted), the gaps in between
are filled using automatic growing and surface interpolation. In
a subsequent step the geological model is built, based on the extracted seismic horizons. For building the geological model, thez-axis of the seismic volume (i.e., the time axis) has to be transformed to the spatial dimension z (i.e., the depth of the volume),
the so called time-depth conversion. Therefore, the seismic volume data is transformed such that the previously extracted horizons
fit and correspond to the additional spatial data provided by welllogs. This is an iterative process. Often errors done during the interpretation phase will show in this step, requiring refinement of theextracted horizon. Because this current workflow is heavily basedon 2D slice processing, this often results in multiple slices that need
to be updated manually. Additionally, in this workflow the selectedaxis-aligned slices used for extracting horizons do not necessarily
correspond to the well positions, making it much harder to fit horizon surfaces to actual well log data.
3.2 Proposed Workflow
Figure 1 displays our integrated and interactive workflow for seismic interpretation with focus on the linked 2D and 3D views forsemi-automatic horizon tracing. Our system has been integrated
into SimVis [11], a visual analysis framework with support forlinked views, advanced data selection techniques, and data importand export functions.
In contrast to other seismic interpretation systems, we propose
an interactive workflow that is based on well positions rather thanaxis-aligned slices. Conceptually, we always work on the subvolume that is defined by three well positions that make up a triangle
on the geological surface and form a triangular prism through theVolume View
shows:
volume renderinghorizons
well log positions
prismsinline / crosslines

horizonupdates

actions:
prism selection

prism
selection

Main Interpretation View
shows:
unfolded prism sides
horizons
well log data
actions:
prism selection
horizon seeding
constraint placement
trigger 3D tracing

constraints

Inline / Crossline Slice Views
shows:
inlines / crosslines
horizons

horizonupdates

actions:
constraint placementnt

Figure 1: Linked-view interaction diagram. The different views of thesystem are coupled to each other, triggering actions upon interaction,and updating the other views accordingly.

a

b

c

a

a

Figure 2: The original prism (left) unfolded into a single slice (right).

seismic volume along the z-direction. Therefore, the well logs that
we use in our system are created by drillings that are roughly orthogonal to the surface, resulting in a unique (x, y)-coordinate foreach well log. This approach was proposed by our collaborating
domain experts. The main advantage of this approach is that thehorizon can easily be adjusted directly at the well position, whichis the point where errors in the interpretation become visible when
adjusting the horizon to the ground truth data from the well log
during time/depth-conversion of the seismic volume. Additionally,
if horizons need to be adjusted during the depth-conversion of thevolume, our semi-automatic method does not require the user tomanually adjust several different slices, but can automatically update only those triangles that are adjacent to the corresponding welllog. In the best case, a single constraint placed at the well position
can be enough to correct all adjacent prisms. After extracting a horizon surface patch for one well log prism, the horizon can then be
extended by processing the neighboring well log prisms, creating
a set of horizon surface patches that constitute the single horizon
surface. The result from one well log will therefore be transferredto the neighbors functioning as constraints for the two- and threedimensional tracing, eliminating the need for seed points in all but
the first prism.
In our approach, the volume is divided into a set of prisms bytriangulating the well positions. The triangles are computed using
Delaunay triangulation, guaranteeing non-intersecting triangles andresulting in the maximum area covered by the well logs, i.e. theirconvex hull. Additionally, we take advantage of the fact that Delaunay maximizes minimum angles. Small angles would result in very
close edges and limit the freedom of the 3D tracing.
The interpretation itself is done triangle (prism) by triangle (prism),
instead of slice by slice. To extend the horizon trace to areas notcovered by wells, additional points of interest (e.g. the corners ofthe volume) can be added to the triangulation.
To compute the seismic horizons we have developed an algorithm that combines the advantages of semi-automatic 2D and 3D
global minimal cost path algorithms. We separate the horizon extraction algorithm into a 2D minimal cost path tracing on prism
faces orthogonal to the (x, y) plane, and a subsequent 3D minimal cost surface computation, which is initialized by the previously
computed 2D contour. For more details on the tracing process refer
to Section 4.1. Implementation details on the 2D part can be foundin Section 4.1.1, 3D tracing is described in Sections 4.1.2 and 4.1.3.
Dividing the tracing into an easy to handle and fast 2D tracing,followed by the slower 3D tracing constricted by the result from the2D trace has the advantage that the user can very quickly produce
a large amount of constraints for the 3D tracing (i.e. the boundary)in the 2D step. This improves the quality of the 3D tracing, minimizing the need to set additional 3D constraints and thus avoids
lengthy recomputation. In our system, the horizon tracing is started
by the user selecting a triangle (i.e., three well logs) in the volume
or timeslice view. Next, the sides of the corresponding prism areshown as a single unfolded slice in the main seismic interpretation
view (see Figure 2). Starting with a single user-defined seed point,a horizon is automatically traced on the sides of each prism. We
employ a minimum cost path algorithm similar to the Live Wirealgorithm presented by Mortensen et al. [22]. Due to the cyclic nature of the prism sides, only a single user-provided input point isFigure 3: Our framework for interactive seismic interpretation using piecewise global energy minimization. The 3D view in (a) shows the volume
alongside well positions, and a complete horizon surface. In the main interpretation view (b), the unfolded sides of a prism are shown overlaid
with well logs. The slice views (c) and (d) show inline and crossline slices, respectively, (e) the timeslice view with the well positions. A horizon
covering the complete volume was traced by adding prism vertices on the corners of the volume. The amplitude at the horizon vertices is plottedon the horizon using a cool to warm colormap.

required, which acts as start as well as end point for the traced 2D
contour. The minimum cost path can be computed in real time fortypical prisms, allowing the initial trace to be visualized instantly inthe main interpretation view as well as in the volume view. In thelatter, a surface inside the edited prism is interpolated right away togive a rough idea of the 3D surface.
Seismic data exhibits a high degree of uncertainty as well as alarge amount of noise. Therefore, we have added the possibility
for the user to add constraints to the automatic tracing algorithm
by specifying additional points that have to be part of the extracted
seismic horizon. New constraints instantly update the horizon andcan be continuously modified by the user.
When a satisfying result for the horizon is reached on the prism
faces, the user can initialize the 3D tracing. The points of the traceon the prism sides are used as constraints for the minimum costsurface algorithm to compute the horizon patch on the inside of theprism. The 3D minimum cost surface algorithm is much more computationally demanding than its 2D counterpart. Thus computation
can take several seconds, but usually does not exceed more than 30
seconds for typical prism sizes. After the computation is finished,
the interpolated surface in the volume view will be replaced by thecorrect minimum cost surface.
Adding constraints requires re-computation of the horizon patch.
However, the comparably small size of the prisms in combination
with the large number of constraints gained from the trace on theprism faces minimize the need for adjusting the 3D traces considerably. If, however, adjustments are needed, the fact that we use a real
3D tracing method over a network of 2D traces is a big advantage
over slice-based approaches, as constraints can be set on the exact
desired position and not only on key slices and every constraint has
impact on the whole surface patch instead of just the current inline
or crossline trace. Additionally, interpreted horizons can be imported and exported in a format readable by the industry standard
software Petrel [34].
We believe that a boundary-constricted minimum cost approach,as described above is the optimal solution for horizon tracing interms of horizon quality, as constraints define horizon elements directly rather than subsurface elements bounded by the horizon. If
needed, the horizon can virtually be traced manually by placing
constraints at every desired grid point, yet most of the time only
very sparse constraints are needed. On the other hand, a global
approach is very costly, making interactively tracing on complete

volumes impossible. By dividing the volume in prisms, though, weare able to retain an interactive workflow with short waiting times
between the 3D tracing of single prisms.
Figure 3 illustrates the integrated nature of our system. All ourviews are linked. If one view is updated, all corresponding views
are updated as well. The 3D volume view (Figure 3a) shows all interpreted horizon parts in addition to the well tops in context of theseismic cube. Embedding the horizon surfaces in the volume rendering allows quick verification of the tracing results. This view canalso be used to select a prism for editing. Figure 3b shows the maininterpretation view. Here the unfolded sides of the selected prism
as well as the corresponding well logs are displayed. For every log,
multiple attributes gathered from the drillings can be visualized as
1D curves. This view is used to set the initial seed point for eachhorizon and allows to modify the boundary of the horizon patch inthe current prism.
Additionally, we have linked the 2D inline and crossline slice views(Figure 3c/d). When the currently displayed inline or crossline slice
intersects the horizon surface, the intersection contour is displayed
immediately. These views can also be used to add constraints for3D tracing. Figure 3e shows a slice view from the top, the so called
timeslice view. Well positions and triangles are displayed and canbe selected here.
4 S EISMIC H ORIZON S URFACE T RACING
Our system is based on automatic global energy minimization ofseismic horizon patches, according to a cost function and userspecified constraints, leading to a piecewise optimized surface. Due
to noise and the ambiguous nature of seismic data in general, theminimal cost path or minimal cost surface might not immediately
match the desired horizon after initialization. Therefore, we allow
the user to add an arbitrary number of constraint points to guide thetracing algorithm, forcing contours and surfaces through selected
points on the horizon.
When adapting techniques like intelligent scissors for seismic interpretation, horizons are defined by the boundaries produced by thealgorithms. As indicated before, precise control over the boundaryis necessary in cases of false classification. This can only be guaranteed by explicit constraints. However, in contrast to graph cuts,
which extend naturally to 3D, a 3D extension of intelligent scissors is not straightforward. Whereas implicit boundary approaches
always operate on image elements, regardless of dimensionality,explicit boundary approaches operate on the boundary elements
which differ between dimensions. Thus, lines in 1D are boundedby points, areas in 2D are bounded by closed curves, volumes in 3D
are bounded by closed surfaces, and so on. Until recently, there wasno real extension of 2D minimal cost path algorithms with explicit
boundary constraints to 3D minimal cost surfaces. Most proposed
methods were based on simple networks of 2D minimal cost paths.
A set of minimal cost paths in a 3D image, yet does not necessarily
define a minimum cost surface. A simple example that illustrates
this is the minimal surface connecting two circles in 3D [16]. Connecting opposing points on both circles with minimal energy results
in a cylinder. However, the minimal cost surface bounded by thesame two circles is a curved shape called a catenoid.

4.1.1 2D Contour Creation

3. The initial surface z0 fitting the boundary acquired in the firststep needs to be inscribed in the graph.

The initial closed contour for each prism is computed as a minimalcost path on the prism sides. For any prism, the three sides are
re-sampled from the volume using texture mapping with tri-linear
interpolation, and stitched together into a single image. This image, superimposed with selected data of the cornering well logs, isshown in the main interpretation view and also used to construct
the graph used for the minimum cost path computation. Figure 4
shows the mapping of an image to primal and dual graph elements.
In contrast to standard approaches in image segmentation [16, 23],
where the minimal-cost path/surface computation is done on the
dual graph, resulting in a boundary in between image elements, theseismic horizons we want to extract are represented in the data bybands of extremal values. Thus, for seismic interpretation we do
not compute the minimal-cost path in between the voxel grid but
directly on the grid. The underlying graph then has to be the primal graph, instead of the dual graph. The resulting minimal cost
path will be a set of edges connecting pixels in the image instead ofdividing pixels.
We have implemented the minimal cost path computation using
Dijkstra’s algorithm [10]. The algorithm is applied on the primal
graph constructed as a 4-connected lattice according to Figure 4,
based on the stitched image of the prism sides. Edge weights aredefined using the cost function defined in Section 4.1.4. The minimal path has to be bounded by at least two nodes, otherwise it will
not contain any edges. To be able to compute the minimal cost pathusing just a single seed point, the cyclic nature of the prism sides
has to be taken into account when creating the graph. In addition to
the regular lattice arcs connecting each node on the right-most sideof the graph with its counterpart on the left-most side and vice versa
are added. To coerce the path around the prism, all arcs directly tothe right of the seed, but pointing to the left are removed. That way,
the single seed point functions as start- and endpoint of the path andit is ensured that the path completely encircles the prism.
Dijkstra’s algorithm cannot handle constraints directly. However, constraints can be added simply by dividing the problem. Instead of trying to compute the complete path at once, the path canbe coerced through any node (the constraint) in the graph by computing the minimal cost path in segments. Every constraint willthen bound two connected segments, one as a start- and one as atarget-node.

4. The MCNF through z0 is computed.

4.1.2 3D Graph Initialization and MCNF

5. The minimal-cost surface is constructed from the set of saturated edges of the circulation, bounded by the initial contour.

MCNF describes the problem of putting as much flow as possible
through the edges of a network at minimal cost, without adding or
draining flow with dedicated source or sink nodes. Therefore, acapacity and a cost is assigned to every edge. As the graph does notcontain a source or a sink, the incoming and outgoing flow at every
node must be equal, thus flow can only circulate in the graph. The
cost is always prioritized over the capacity of edges. If a cycle adds
cost to the total flow network, it will not be used, regardless of its
capacity. The capacity of a cycle is defined by the lowest capacityof the edges it consists of and the total flow in the graph is definedby the accumulated flow of all cycles.
To compute the maximum flow in a prism, the dual graph is created by using all voxels contained in the prism, according to thevolume-graph scheme shown in Figure 4. The graph is constructedin two steps: first the basic graph structure of the dual is created(Figure 5a) and costs and capacities are assigned to each edge (Figure 5b). In this step, a cost of 1 is used for all edges (not indicatedin the figure). The capacities of the dual edges correspond to thecosts of the primal facets and are computed using Equation 3 inSection 4.1.4. The edges contained in the outer facets in the graph,have to be treated separately. These edges do not have all four adjacent voxels and thus the cost cannot be computed using the costfunction. We initialize these edges with infinite capacity. Edge ca-

4.1 Minimal Horizon Surface Computation
Our approach for computing horizon surfaces using global energy
minimization is based on previous work of Grady [16]. In that
paper, Grady formulates the minimization problem for boundary
constricted 3D minimum cost surface computation and shows thatminimal-cost circulation network flow (MCNF) is an optimized solution for the dual of this particular problem. By circulating the
maximum possible amount of flow in the dual graph through any
initial surface z0 fulfilling the boundary constraints, eventually a
set of edges will be saturated. This forms a cut, which corresponds
to the minimal cost surface in the primal. The process is similar
to the maximum flow, minimum cut method used for graph cuts,with the type of constraints being the main difference. Graph cuts
define fore- and background elements in the image as constraints,
which then form the source and sink nodes in the graph, whereas
for MCNF, with no source or sink nodes present, the constraints areadded as the boundary of the initial surface z0 .
Our algorithm for computing minimal horizon surfaces consists ofthe following main steps:
1. A closed contour bounding the desired surface within a prism
has to be computed. This is done in a flattened representationof the prism’s sides (see Figure 2).
2. Next, the volume contained in the prism has to be represented
as a weighted, directed graph.

6. If the surface is not satisfactory it can be forced through desired points in the volume by placing additional constraints.
These constraints are then added as additional (inner) boundary patches in step 3. Steps 4 and 5 need to be reapplied.
Below, the steps of our algorithm are explained in more detail.
Pixel

PrimalDual 2D

Dual 3D

facet

cube

2D Lattice

node

edge

facet
edge

facet
node

3D only:

cube

edge

node

Figure 4: Correspondence between picture elements and items in
the primal and dual graph.

aDual Graph
Primal Graphb
Voxels - grey value
indicates amplitudeccapacity denoted byedge thickness∞ on sides-ncost & capacity
initialization

z surfaceinitialization

Dual Edgeon projected
Dual Facet

Prism

-n

--n
n

-n

extract min
cost surface

MCNF
MC
n
-n

2D sliceused forillustration

cost 1 for all edgeseminimal path/surface candidates
in primal corresponding to saturated dual edges

0Primal Nodes
belonging toboundary

Dual Noded2 alternate
valid paths

cut backpointing arc for dual edges
corresponding to initial surface

saturated
dual
edges

cost -n for dual edges
corresponding to initial
surface to induce flow
initial surface
in primal graphdead end path
separate patch
minimalpath/surface

2D horizon functions
as surface boundary

g

f

h

Figure 5: Tracing pipeline.

pacities are indicated in Figure 5b by the thickness of the edge. At
this point there will be no flow in the graph, because all edges havepositive cost. In the second step, the initial surface z0 is added to thegraph. The initial surface is a set of facets in the primal, boundedby the closed contour obtained from the prism sides.In Figure 5c,the boundary is indicated by the two purple primal nodes, the initial
surface by the path between these nodes. To induce flow throughthe initial surface, negative costs are assigned to the dual edges corresponding to the primal facets of z0 . To allow all possible paths,
the cost is the negative of the sum of all positive costs in the graph,such that every path through z0 is negative and thus reducing total cost for the circulation. Additionally, we allow flow throughthe surface only in one direction by cutting all edges in the opposite direction. This guarantees that no cycles going back and forththrough the surface are generated. Here it becomes obvious why weassign infinite capacity to the outer edges in the first step. All allowed cycles will go through the initial surface and the inner edges
in one direction and come back on the (much fewer) edges on thefaces of the prism. If the capacity of these edges were delimited, thealgorithm would most likely terminate, delivering a set of saturated
edges on the outside of the prism which would not be helpful forfinding the minimal-cost surface.
For most image segmentation approaches, the initial 2D trace
will be performed on one slice in the volume, resulting in a planarboundary. In that case, the part of the slice bounded by the obtained
2D trace is the simplest match for the initial surface z0 . The 2D
trace in our approach, does not result in a single planar boundary,however, a very simple initial surface can be constructed dividing
the surface into four planar parts. First we cut the prism above thehighest point in the trace with an orthogonal slice, resulting in atriangle shaped boundary in that slice. This triangle functions as
the cover of our surface. The other three parts are on each sideof the prism, bounded by the initial trace on the bottom, the coverslice on the top and the prism edges on the sides. Figure 5g showsa sketch of an initial surface fitting a 2D horizon path on the sides
of a prism, the cover is depicted in light, the sides in dark cyan, the2D trace on the prisms sides by the purple curve.
For the actual MCNF computation we use the LEMON C++
graph library [1], which can compute the flow for all edges in thegraph. This is used to identify the saturated edges in the graph.
4.1.3 Horizon Extraction
The set of saturated edges does not translate directly to the desired
minimal cost surface. In real world datasets, usually a large amountof edges not contributing to the minimum cost surface are saturated.
E.g., in our application, tracing a horizon in a prism of a thousand
points per layer can easily result in over 10,000 saturated edgesof which 2,000 contribute to the actual surface. Figure 5d shows
the cases of saturated dual edges which have to be removed. Dualedges/primal facets not connected to the initial boundary could simply be avoided by growing the surface from the initial boundary.
However, dead end paths/surface patches, connected to the surface
would not be removed by growing.
Instead we prune facets by checking that all four of the facet’s
edges are either shared with another facet candidate, or are partof the initial boundary. By iterating over all facets and discarding
facets which do not fit these criteria, the set of surface candidates
is gradually reduced. Once a pass over all facet candidates finishes without finding any more facets that need to be discarded, theremaining facets form the minimal-cost surface and the algorithm
returns (Figure 5e/h). If two alternate surfaces bounded by the initial contour exist (top left part of the path in Figure 5d), both must
be of the same cost and thus are equivalent solutions. As the user
can always force one of these we just use the first complete closed
surface found by the algorithm.
4.1.4 The Cost Function
We define edge and facet weights using a cost function containing
two components.
The first component defines the snappyness of the surface toridge and valley lines/surfaces in the image/volume. It calculates
the deviation of the grey values of all current voxels (i.e., the voxels corresponding to the current edge or facet whose cost is to be
defined) from a specified target value. We allow tracing of minima, maxima and zero crossings. Before tracing the kind has to bechosen, resulting in a suitable target value (e.g., for a [−1.0..1.0]scaled volume −1.0 for minima, 0.0 for zero crossings and 1.0 for
maxima).
The second component of the cost function is important for thesmoothness of the traced line/surface. It is defined by the similarityof the voxels belonging to the edge/facet in question. The higher
the difference of voxel grey values, the higher the cost. If needed,
the relative weight of the two cost function components can be adjusted.
The dimensionality of the boundary items in the primal graphvaries from dimension to dimension. In the dual, however thebounding item is always an edge. Thus we define the cost using
the voxels adjacent to an edge in the dual graph.

The snappyness component is defined by:m

w1 = ∑ |t − f (vk )| ,

(1)

k=1

where f (v) is the grey value of the image element corresponding tonode v, t is the target amplitude, and m is the number of adjacentimage elements (i.e. two pixels in 2D, four voxels in 3D). Accordingly,m
1 m
w2 = ∑ f (vk ) − ∑ f (vl )
(2)
m l=1
k=1
describes the smoothness component of the cost function andw = cw1 + (1 − c)w2

(3)

the complete cost function using Equations 1 and 2, including a
constant c ∈ [0..1] for linear blending between w1 and w2 .
5 S EISMIC H ORIZON V ISUALIZATION
Visualization of the extracted horizon surface is done simultaneously in all 2D and 3D views. The 2D inline and crossline views
display a slice of the volume with overlaid horizon data. All extracted horizons are displayed as 2D contours, with the currently
active horizon and prism being highlighted. The current horizon
can be modified at any time by clicking on and dragging the contour. After the first step of the minimal cost path tracing algorithm,
the inner part of the horizon is linearly interpolated, whereas after
3D surface tracing, the actual minimal surface is displayed in the2D view.
The 3D view combines direct volume rendering of the seismic
cube with the visualization of the horizon surface geometry and thewell logs. Volume rendering is performed in a single-pass GLSL
raycaster. For correct visibility of the horizon surfaces and welllogs we simply adjust the depth values of the ray stop-positions
used in raycasting, and draw the geometry prior to the volume [32].
Additionally, the 3D view enables picking of well log prisms andof already extracted horizon surfaces and is linked to the seismic
interpretation-, inline-, and crossline view. The 3D view is usedfor visual inspection of the ongoing seismic interpretation and supports plotting of several horizon properties (i.e., amplitude, deviation from the target amplitude, cost, and the distance to a second
horizon) directly on the horizon’s surface geometry. This additional
information is computed on-the-fly during rendering using a combination of vertex and fragment shaders. Deviations and cost arecomputed as scalar values, scaled to [0..1] and mapped onto a blue
to red colormap presented in [21] for easy interpretation. The amplitude can either be visualized as simple grey values or with thesame blue to red scheme used for the other scalar values. However,here blue is mapped to minus one and red to one. Surfaces withcost and distance to another surface, respectively, can be seen in the
results section (Figure 6).
6 R ESULTS
To support an interactive workflow for seismic interpretation, we
have integrated our system into the SimVis visual analysis framework. The system was recently installed at our domain expert collaborators for evaluation. Figure 6 shows 3 horizons extracted froma 240x240x1509 voxels seismic dataset using the proposed technique. The horizon was extracted in the area covered by 8 wells,
covering roughly one fourth of the dataset. All the horizons weretraced using a single seed point plus between five and eight constraints. All constraints were set on the sides of the prism, noneon the inside. For details see Table 1. For comparison, we plotted
the distance to the same horizons manually traced by our domain
experts onto the surfaces using the cold to warm color map presented in [21]. For the first two examples, the distance was at most
2.39 pixels. The second horizon showed a very strong reflection
in large parts of the volume, resulting in the shorter editing timeabcConstraints
1+8
1+5
1+5

Edit t/min
<5
<2
<5

Tracing t/s
61
53
74

Max Dist.
2.39px2.27px4.30px

Table 1: Details for the horizons shown in Figure 6. Number of constraints (seed + additional constraints), editing time, tracing time andmaximum deviation from manually traced ground truth data.

as well as faster tracing. For the third example, the maximum distance was over four pixels. However, it should be noted that thearea of largest distance in the last example seems to actually be ofhigher accuracy with our approach. Figure 6f shows a cutout of
the slice views showing the horizon trace (supposed to be a maximum trace) with our method in magenta and the manually traced
horizon in blue. The manually traced horizon is further away fromthe area of maximum amplitude to allow for an overall smoothersurface. Outside of the area of interest the manual interpretation
was done very sparsely resulting in a rough approximation of thehorizon only. With a (semi-)automatic approach like ours the horizon can be traced exactly over the complete volume at low cost. It
should be noted that with placement of additional constraints, thedistance to the manually traced surface can be further reduced. In
theory, it would even be possible to completely mimic the manualinterpretation process by placing constraints on every grid point.
Rendering of the interpreted horizons alongside the volume ina 3D view happens in real time. Rendering a dataset of size
240x240x1509 with three interpreted surfaces and additional wellgeometry results in more than 30 fps on a Geforce 280 GTX at fullscreen resolution (1920x1080 pixels).
7 C ONCLUSIONS
We have presented an interactive workflow for rapid interpretationof seismic volumes using a combination of 2D and 3D cost minimization techniques. The workflow guides the user through the interpretation process using a combination of multiple linked 2D and
3D views. The prism-based workflow, instead of working along
inline and crossline directions only, was received very well by ourcollaborators. It bounds the size of the working volume, and allows incrementing the piecewise optimal horizon surface bit by bit,by subsequently appending additional optimal surface patches thatwere created by processing neighboring well log prisms. Decoupling the interaction in 2D and the eventual computation in 3D was
also perceived very well by the experts, as it simplifies user interaction and time expenditure considerably. Visualization of the extracted horizons, superimposed with additional information and incontext of the original 3D volume, gives instantaneous feedback on
the quality of the seismic interpretation. Additionally, our approachis scalable to large seismic volumes, because the time needed forcomputing minimal-cost paths and surfaces depends on the size ofthe selected prism instead of the volume size.
ACKNOWLEDGEMENTS
This research project was funded in parts by the Austrian Research Funding
Agency (FFG), in scope of the SEIVIS project (No. 818063). The datasets
are courtesy of Heinemann Oil GmbH.

R EFERENCES[1] LEMON library for efficient modeling and optimization in networks:
http://lemon.cs.elte.hu.
[2] C. Armstrong, B. Price, and W. A. Barrett. Interactive segmentation ofimage volumes with live surface. Computers & Graphics, 31(2):212–
22, 2007.
[3] A. Blinov and M. Petrou. Reconstruction of 3-d horizons from 3D seismic datasets. IEEE Transactions on Geoscience and Remote
Sensing, 43(6):1421–1431, 2005.
[4] H. G. Borgos, O. Gramstad, G. V. Dahl, P. L. Guern, L. Sonneland,
and J. F. Rosalba. Extracting horizon patches and geo-bodies from3d seismic waveform sequences. SEG Technical Program Expanded
Abstracts, 25(1):1595–1599, 2006.

Figure 6: Three horizon patches computed using our algorithm. (a),(b) and (c) show the distance in pixels from a manual trace of thesame horizon. (d), (e) and (f) show detailed slice view cutouts ofthe surface indicated by the green boxes in (a), (b) and (c). Manual
traces are in cyan, traces computed with our method in magenta. (d)shows a region of the first horizon, with > 2 pixel difference, where
the manual trace is clearly better than the automatic trace, manualand automatic traces in (e) are very similar. (f) shows a part of thehorizon with > 4 pixel difference. Here, the automatic trace is closer
to the actual ridgeline.
[5] Y. Boykov and M.-P. Jolly. Interactive organ segmentation using graphcuts. In MICCAI ’00: Proceedings of the Third International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 147–175, 2000.
[6] Y. Boykov and M.-P. Jolly. Interactive graph cuts for optimal boundary& region segmentation of objects in n-d images. In ICCV 2001. Proceedings of the Eighth IEEE International Conference on Computer
Vision, pages 105–112 vol.1, 2001.
[7] Y. Boykov and V. Kolmogorov. An experimental comparison of mincut/max-flow algorithms for energy minimization in vision. In EMMCVPR ’01: International Workshop on Energy Minimization Methods
in Computer Vision and Pattern Recognition, pages 359–374, 2001.
[8] Y. Boykov, O. Veksler, and R. Zabih. Fast approximate energy minimization via graph cuts. Pattern Analysis and Machine Intelligence,
IEEE Transactions on, 23(11):1222–1239, 2001.
[9] L. Castanie, B. Levy, and F. Bosquet. Volumeexplorer: Roaming large
volumes to couple visualization and data processing for oil and gas
exploration. In Proceedings of IEEE Visualization Conference ’05,
pages 247–254, 2005.
[10] E. Dijkstra. A note on two problems in connexion with graphs. Numerische Mathematik, 1(1):269–271, 1959.
[11] H. Doleisch, M. Gasser, and H. Hauser. Interactive feature specification for focus+context visualization of complex simulation data. In
Proceedings of the 5th Joint IEEE TCVG - EUROGRAPHICS Symposium on Visualization (VisSym 2003), pages 239–248, 2003.
[12] A. Falcão and J. Udupa. A 3d generalization of user-steered live-wiresegmentation. Medical Image Analysis, 4(4):389–402, 2000.
[13] A. X. Falcão and J. K. Udupa. Segmentation of 3d objects using live
wire. In Medical Imaging 1997: Image Processing, pages 228–235,
1997.
[14] M. Faraklioti and M. Petrou. Horizon picking in 3d seismic data volumes. Machine Vision and Applications, 15:216–219, 2004.
[15] L. Grady. Computing exact discrete minimal surfaces: Extending andsolving the shortest path problem in 3d with application to segmentation. In Computer Vision and Pattern Recognition, 2006 IEEE Com-

puter Society Conference on, volume 1, pages 69–78, 2006.
[16] L. Grady. Minimal surfaces extend shortest path segmentation methods to 3d. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 32(2):321–334, 2010.
[17] N. Keskes, P. Zaccagnino, D. Rether, and P. Mermey. Automatic extraction of 3-d seismic horizons. SEG Technical Program Expanded
Abstracts, 2(1):557–559, 1983.
[18] O. D. Lampe, C. Correa, K.-L. Ma, and H. Hauser. Curve-centric
volume reformation for comparative visualization. IEEE Transactions
on Visualization and Computer Graphics, 15(6):1235–1242, 2009.
[19] P. Lavest and Y. Chipot. Building complex horizons for 3-d seismic.
SEG Technical Program Expanded Abstracts, 12(1):159–161, 1993.
[20] E. M. Lidal, T. Langeland, C. Giertsen, J. Grimsgaard, and R. Helland.
A decade of increased oil recovery in virtual reality. IEEE Computer
Graphics and Applications, 27:94–97, 2007.
[21] K. Moreland. Diverging color maps for scientific visualization. In
Proceedings of the 5th International Symposium on Visual Computing,pages 92–103, 2009.
[22] E. N. Mortensen and W. A. Barrett. Intelligent scissors for image
composition. In SIGGRAPH ’95: Proceedings of the 22nd annual
conference on Computer graphics and interactive techniques, pages
191–198, 1995.
[23] E. N. Mortensen and W. A. Barrett. Interactive segmentation
with intelligent scissors. Graphical Models and Image Processing,60(5):349–384, 1998.
[24] E. N. Mortensen, B. Morse, W. A. Barrett, and J. K. Udupa. Adaptive
boundary detection using ‘live-wire’ two-dimensional dynamic programming. In Computers in Cardiology 1992. Proceedings., pages
635–638, 1992.
[25] D. Patel, S. Bruckner, I. Viola, and M. E. Gröller. Seismic volume visualization for horizon extraction. In Proceedings of the IEEE Pacific
Visualization Symposium 2010, pages 73–80, 2010.
[26] D. Patel, C. Giertsen, J. Thurmond, J. Gjelberg, and M. E. Gröller. The
seismic analyzer: Interpreting and illustrating 2d seismic data. Visualization and Computer Graphics, IEEE Transactions on, 14(6):1571–
1578, 2008.
[27] D. Patel, Ø. Sture, H. Hauser, C. Giertsen, and M. E. Gröller.
Knowledge-assisted visualization of seismic data. Computers &Graphics, 33(5):585–596, 2009.
[28] R. Pepper and G. Bejarano. Advances in seismic fault interpretation
automation. Search and Discovery Article 40170, Poster presentation
at AAPG Annual Convention, pages 19–22, 2005.
[29] M. Poon, G. Hamarneh, and R. Abugharbieh. Live-vessel: Extending
livewire for simultaneous extraction of optimal medial and boundary
paths in vascular images. In Lecture Notes in Computer Science, Medical Image Computing and Computer-Assisted Intervention (MICCAI),
pages 444–451, 2007.
[30] M. Poon, G. Hamarneh, and R. Abugharbieh. Segmentation of complex objects with non-spherical topologies from volumetric medical
images using 3d livewire. In SPIE Medical Imaging, volume 6512-31,
pages 1–10, 2007.
[31] M. Poon, G. Hamarneh, and R. Abugharbieh. Efficient interactive 3d
livewire segmentation of objects with arbitrarily topologies. Computerized Medical Imaging and Graphics, 32(8):639–650, 2008.
[32] H. Scharsach, M. Hadwiger, A. Neubauer, and K. Bühler. Perspective
isosurface and direct volume rendering for virtual endoscopy applications. In Eurovis 2006, pages 315–322, 2006.
[33] A. Schenk, G. Prause, and H. O. Peitgen. Efficient semiautomatic segmentation of 3d objects in medical images. In MICCAI ’00: Proceedings of the Third International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 186–195, 2000.
[34] Schlumberger Information Solutions. Petrel seismic to simulation
software: www.slb.com/services/software/geo.
[35] P. Tu, A. Zisserman, I. Mason, and I. Cox. Identification of eventsfrom 3d volumes of seismic data. In IEEE International Conference
on Image Processing, ICIP ’94, volume 3, pages 309–313, 1994.
[36] U.S. Energy Information Administration. International energy outlook, 2010.
[37] Visualizaton Sciences Group. Avizo earth: www.vsg3d.com/
avizo/earth.