CobWeb ― a toolbox for automatic tomographic

image analysis based on machine learning

techniques: application and examples
Swarup Chauhan1,2*, Kathleen Sell2,5, Frieder Enzmann2, Wolfram
Rühaak3, Thorsten Wille4, Ingo Sass1, Michael Kersten2
1Institute of Applied Geosciences, University of Technology, Darmstadt 64287, Germany
2 Institute for Geosciences, Johannes Gutenberg-University, Mainz 55099, Germany
3 Federal Institute for Geosciences and Natural Resources (BGR), Hannover 30655, Germany
4
APS Antriebs-, Prüf- und Steuertechnik GmbH, Götzenbreite 12, Göttingen-Rosdorf 37124, Germany
5igem – Institute for Geothermal Ressource Management, Berlinstr. 107, Bingen 55411, Germany

Abstract. In this study, we introduce CobWeb 1.0 which is a graphical user interface tailored
explicitly for accurate image segmentation and representative elementary volume analysis of
digital rock images derived from high resolution tomography. The CobWeb code is a work
package deployed as a series of windows executable binaries which use image processing and
machine learning libraries of MATLAB®. The user-friendly interface enables image segmentation
and cross-validation employing K-means, Fuzzy C-means, least square support vector machine,
and ensemble classification (bragging and boosting) segmentation techniques. A quick region of
interest analysis including relative porosity trends, pore size distribution, and volume fraction of
different phases can be performed on different geomaterials. Data can be exported to ParaView,
DSI Studio (.fib), Microsoft® Excel and MATLAB® for further visualisation and statistical
analysis. The efficiency of the new tool was verified using gas hydrate-bearing sediment samples
and Berea sandstone, both from synchrotron tomography datasets, as well as Grosmont carbonate
rock X-ray micro-tomographic dataset. Despite its high sub-micrometer resolution, the gas hydrate
dataset was suffering from edge enhancement artefacts. These artefacts were primarily normalized
by the dual filtering approach using both non-local means and anisotropic diffusion filtering. The
desired automatic segmentation of the phases (brine, sand, and gas hydrate) was thus successfully
achieved using the dual clustering approach
.
Keywords:
X-ray computed tomography, machine learning, graphical user interface, edge enhancement, LSSVM, REV, gas hydrates, Berea
sandstone, Grosmont carbonate

* Corresponding author
E-mail: schauhan@uni-mainz.de
Swarup Chauhan designed and developed the graphical user interface (GUI) and performed Analysis
Kathleen Sell performed XCT experiments of gas hydrate (GH) and performed image filtration on GH XCT datasets
Frieder Enzmann conducted the gas hydrate XCT experiments
Wolfram Rühaak helped in optimizing the GUI
Throsten Wille shared his expertise on XCT experiments and supported the research study
Ingo Sass supervised the research studies, proof read and edited the manuscript
Michael Kersten supervised the research studies, proof read and edited the manuscript
 2

1 Introduction
Despite the availability of both commercial and described in Section 2. Section 3 highlights the
open source software for digital rock physics (DRP) image segmentation techniques used in Section 4
analysis as compiled in Figure 1, an ideal tool for for REV, relative porosities, and pore size
accurate automatic image analysis at ambient distribution analysis.
computational performance is difficult to pinpoint.
The best practice so far among researchers is to 2 CobWeb 1.0
alternate between different available software tools
and to synthesize the different datasets using home- 2.1 Salient Features
brew workflows. Porosity and particularly The word Cobweb means “a tangled three-
permeability can vary dramatically with small dimensional spider web”, i.e. something resembling
changes in segmentation, as significant features can a cobweb in delicacy or intricacy (Oxford
be lost when thresholding greyscale tomography Dictionaries). According to Marriam-Webster, it
images to binary images, even if using the most may also mean something that entangles obscures
advanced data acquiring techniques like or confuses, as is the philosophy of machine
synchrotron tomography (1). Our new CobWeb 1.0 learning - elegant, sophisticated yet stochastic and
visualisation and image analysis toolkit addresses confusing. This inspired us to name our software
some of the challenges of selecting representative tool CobWeb. The first version enables reading and
elementary volume (REV) for X-ray computed processing of (reconstructed) XCT files in both tiff
tomography (XCT) datasets reported earlier (1-7). It and raw formats. The graphical user interface (GUI)
is customized to perform image analysis and is embedded with visual inspection tools for
accurate greyscale phase segmentation of zooming in/out, cropping, color, and scale, and
reconstructed high resolution XCT and synchrotron assisting in the visualisation and interpretation of
tomographic datasets. As the only one currently 2D and 3D stack data. Noise filters such as non-
available, it is based on machine learning techniques local means, anisotropic diffusion, median and
of excellent performance for segmentation analysis contrast adjustments are implemented to increase
as detailed previously (8, 9). This software tool signal-to-noise ratio. The user can choose from a
package was developed on a MATLAB® series of five different segmentation algorithms,
workbench and can be used as Microsoft Windows namely K-means, Fuzzy C-means (unsupervised),
standalone executable (.exe) files or as a least square support vector machine (supervised),
MATLAB® plugin. In this paper, we demonstrate bragging and boosting (enable classifiers) for
exemplarily 3D tomographic REV analysis of Berea accurate automatic segmentation and cross-
sandstone, Grosmont carbonate rock, and gas validation. Relevant material properties like relative
hydrate-bearing sediment datasets. For the latter porosities, pore size distribution, volume fraction
geomaterial, Sell et al. (10) highlighted problems (pore, matrix, mineral phases) can be quantified and
with the edge enhancement (ED) effect and visualized as graphics output. The data can be
recommended image morphological strategies to exported into different file formats such as
compensate for this artefact. In this paper, we Microsoft® Excel (.xlsx), MATLAB® (.mat),
suggest a strategy to eliminate ED artefacts using ParaView (.vkt) and DSI studio (.fib). The current
the same dataset with the machine learning version is supported for Micosoft® Windows PC
approach. The respective MATLAB code is operating systems (Windows 7 and 10).
provided in the appendix. The salient features of
CobWeb 1.0 and its overall framework are
 3

Figure 1. Market survey of the currently available commercial software (a) and open source software (b) assisting in digital rock
physics analysis with features as indicated in the legend.
 4

Figure 2. Screenshots of the CobWeb GUI. XCT stack of Grosmont carbonate rock is shown as an example of representative elementary
volume analysis. The top panel displays the XCT raw sample, the K-means segmented ROI, and the porosity of single slice No. 10. The
bottom plot shows pore size distribution of the complete REV stack, the relative porosity and volume fraction, respectively.
 5

Figure 3. The general workflow of the CobWeb software tool, where the arrow denotes the series in which different
modules are compiled and executed. A separate file script is used to generate .dll binaries and executables.

history subpanel is a uilistbox that displays errors,

2.2 Window Panels
processing time/status, processing instruction, files
The main GUI window panel divides into three generated/exported and executed callbacks. The
main parts (Figure 2), the tool menu strip, the control subpanel is an assemblage of uibuttons to
inspector panel, and the visualisation panel. The tool initialize the XCT data analysis process and the
strip contains menus for zoom in and out, pan, rotate, progress bar. The visualisation panel is where the
point selection, color bar, legend bar, and results are displayed in several resized windows,
measurement scale functionalities. The inspector which can be moved, saved and deleted. The window
panel is divided into subpanels where the user can panels displayed in the visualisation module are
configure the initial process settings such as embedded with uimenu and submenu to export, plot
segmentation schemes (supervised, unsupervised, and calculate different variables like porosity, PSD,
ensemble classifiers), filters (contrast, non-local volume fraction, entropy, or receiver operational
means, anisotropic filter, fspecial), and distance characteristics. To get the desired user functionalities,
functions (link distance, Manhattan distance, box MATLAB® internal uilibaries were used and where
distance) to assist segmentation and geometrical inadequate. Therefore, numerous specific adaptions
parameter selection for image analysis (REV, were adopted from Yair Altman’s undocumented
porosity, PSD, volume fraction). The display Matlab website and the Matlab File Exchange
subpanel records and displays the 2D video of the community. Specifically, the GUI Layout Toolbox
XCT stack and the histogram of the raw image(s). The from David Sampson has been used to configure the
 6

CobWeb GUI layout; the uitable, which uses the input comprises a desired cluster number, the given
MATLAB java-component was designed using image resolution, and the slice number. Afterwards,
uitable customisation reoprt provided by (10) . the chosen slice is displayed in 2D format, and in a
separated window. Following this, an option of ROI
2.3 Overall Framework selection is proposed to be accepted or denied. If the
An overview of the different modules of the user accepts, the REV will be cropped and updated to
CobWeb toolkit are compiled in Figure 3, and the the main structure and the cropped ROI is displayed
arrow displayed indicates the series in which they are in a separate window as a slice.
executed. The advantage of using MATLAB® is the Based on the option selected in the uitable, the
access to the structure and respective variables, which respective unsupervised and supervised loop is
are used for further investigations. As a stand-alone, initialized. If LSSVM or Bragging and Boasting is
the GUI can be executed on different PC and HPC chosen as the segmentation scheme, a right click
clusters without any license issues. The CobWeb 1.0 uimenu is initialized with the options of pixel
framework can be broadly divided into three modules selection, training, and testing. On pressing the Pixel
Section option, the subroutine uPixelSel() initiates a
2.3.1 Control Module uitable window representing the columns Clusters,
Features, X-Coordinate, and Y-Coordinate, and
In the control module, the CobWeb menu creates requiring user inputs in the respective columns. The
the main figure panel, assembles the size/position of user can explore interesting features (pores, minerals,
the panel and subpanel windows, initializes the matrix, noise/specks) in the 2D image windows and
control buttons and generates a main structure. collect the data using zoom in, zoom out, and data
Ideally, any button can be activated after the GUI is cursor tools. Once the respective features’ X,Y
displayed, but an exception will be displayed in the coordinates are fed in the uitable, the data has to be
history window, highlighting the next step. That step arranged and exported for training and testing. This is
is to load data where the Load function checks the file fulfilled by pressing the export button which initiates
properties, loads the data in tiff and raw formats, the subrountine uExportTable(). The export
creates and displays 2D video of the selected stack, subroutine collects a total of 36 (6 x 6) pixel values in
saves the video file in the current folder, and updates the perimeter of the X,Y coordinates of the respective
the respective variables to the main structure. The features given in the uitable. Thereafter, with the
Stop function ends the execution. However, when the training and testing options, respective models
processing is inside a loop, the Stop function breaks (LSSVM, Ensemble Classifiers) are trained using the
the loop after the i-th iteration. pixel values of the representative slice. The
classification is then performed on the 3D stack, and
2.3.2 Analysis Module
the main structure is updated. In the case of
The Start function is in the analysis category and unsupervised techniques based on the option selected
is a densely nested function, where the bullet points for image filtering, segmentation, and distance
and the indented bullet points represent the outer and function, the complete stack is processed. For FCM,
the inner nested loops, respectively (Figure 3). the user is given an option to choose the membership
Initially, data is gathered and a sanity check is criteria (7) between the range of one to two (decimal
performed to evaluate whether all the options are values).
correctly selected. If the conditions are not satisfyied, The progress can be monitored in the progress
an exception alert pops up in the History panel, bar, the color of the control buttons (red to grey) and
highlighting the error and offering an alternative the History window, which gives the related
process. The second loop is the image modification information on processing time, segmentation
loop, where initially the user input is required. This scheme, and filter options executed.
 7

which can save current plots to different image

2.3.3 Visualisation Module
formats (bmp, matlab, tiff, Jpeg, ect.) and closes the
The visualisation module consists of the plot current window.
function and volume rendering function. The plot
function is a densely nested function, where the bullet 3 Tomography Datasets Used for Evaluation of
points and the indented bullet points in represent the CobWeb 1.0
outer and the inner nested loops, respectively. The 3.1 Gas Hydrate-Bearing Sediments
execution of the Plot() callback accesses the main
structure, identifies the representative slice number, The in-situ synchrotron-based tomography
and plots the segmented 2D image in a new window. experiment and post-processing of synchrotron data
The displayed window is embedded with right click conducted to resolve the microstructure of gas
menu and submenus, namely: hydrate-bearing (GH) sediments is given in detail by
(11-12), and (9). In brief, the tomographic scans were
 Porosity → Porosity, Pore Size Distribution acquired with a monochromatic X-ray beam energy
 Validate → Entropy, Receiver Operation of 21.9 KeV at Swiss Light Source synchrotron
Characteristics, 10-fold Cross Validation facility (Paul Scherrer Institute Villigen, Switzerland)
 Export Stack → Paraview Format using the TOMCAT beamline (Tomographic
When choosing the desired options mentioned above, Microscope and Coherent Radiology Experiment)
their respective subroutines (uPoreSzVol, uCalVal, (13). Each tomogram was reconstructed from
uExport) are started and the desired results are plotted sinograms by using the gridded Fourier
as shown in Figure 2. The user has the option to transformation algorithm (14). Adjacent to this, a
choose a single slice or the complete stack. The three-dimensional stack was derived resulting in an
visualized parameter distributions display in a image size of 2560 x 2560 x 2160 voxels with a voxel
separate window and can further be exported using resolution of 0.74 μm and 0.38 μm at 10-fold and 20-
right-click menu into Excel, ASCII, or MATLAB fold optical magnification, respectively.
formats.
The Export Stacks option bundles the complete 3.2 Grosmont Carbonate Rock
stack in ParaView code format (.vkt files) using the
The datasets of the Grosmont carbonate rock
subroutine function Vtkwrite. When the VolRender
were obtained from the GitHub FTP server
button is activated, the subroutine checks for the size
(http://github.com/cageo/Krzikalla-2012) provided
of the stack and offers an option to perform volume
for the benchmark project reported by Andrä et al.
reduction of the greyscale values. If ignored,
(15, 16). They acquired the Grosmont carbonate rock
depending on the computer RAM capacity, it may
for their benchmark test from the Grosmont
take a relatively long time to plot a high quality
Formation in Alberta, Canada. The Grosmont
volume rendered figure of the 3D stack. If accepted,
Formation was deposited during the Late Devonian
the pixel information is reduced 10-fold, which
and is divided into four facies members, LG UG-1,
fastens the plotting process, but the image quality is
UG-2, and UG-3 from the bottom up. The sample was
hampered. Therefore, CobWeb 1.0 offers the option
taken from UG-2 facies and is mostly composed of
to export the stack and visualize it using ParaView or
dolomite and karst breccia (17, 18). Laboratory
DSI studio.
measurements of porosity and permeability reported
2.3.4 Miscellaneous Options in Andrä et al. (16) are 0.21 (21 %) and 150 mD ─
470 mD, respectively. The Grosmont carbonate
CobWeb 1.0 has a clear button to reset the dataset was measured at the high-resolution X-
CobWeb by delegating back to the main panel. This ray computer tomographic facility of the University
corresponding function (Clear) clears the figure axis, of Texas with an Xradia MicroXCT-400 instrument
the history window panel, and resets the progress bar. (ZEISS, Jena, Germany). The measurement was
The plot windows embed, save, and close fbuttons performed using 4x objective lenses, 70 kV
 8

polychromatic X-ray beam energy and a 25 mm CCD obtained using a Helium pycnometer 1330
detector. The tomographic images were reconstructed (Micromeritics Instrument Corp., Germany) and a
from the sinograms using proprietary software and Pascal 140+1440 Mercury porosimeter (Thermo
corrected for the beam hardening effect as typical for Electron Corporation, Germany) as described by
lab-based polychromatic cone-beam X-ray Giesche (20). The permeability reported in the same
instruments (18). The retrieved image volume was benchmark test (16) ranges between 200 mD and
cropped to a dimension of 1024 x 1024 x 1024 with a 500 mD. Machel and Hunter (17) reported for this
voxel resolution of 2.02 μm. sample a mineral composition of Ankerite, Zircon, K-
feldspar, Quartz, and Clay using a polarized optical
3.3 Berea Sandstone Rock microscope and a scanning electron microscope. The
The Berea sandstone dataset was also obtained synchrotron tomographic scans of Berea sandstone
from the GitHub FTP server provided for the were also acquired at the SLS TOMCAT beamline
benchmark project reported by Andrä et al. (15, 16). (13). The beam energy was monochromatized to 26
The Berea sandstone sample plug was acquired from keV for optimal contrast with an exposure time of 500
Berea sandstone TM Petroleum Cores (Ohio USA). ms. This resulted in a 3D tomographic stack of
Porosity values of around φ = 0.20 (20 %) were dimension 1024 x 1024 x 1024 and voxel resolution
of 0.74 μm.

4 Results and Discussion 4.1.1 Dual Filtering of Gas Hydrate-Bearing

Sediments
4.1 Image Pre-Processing
Due to ED artefacts affecting the quality of the
XCT and synchrotron tomographic datasets were
hydrate-bearing sediment tomograms, the data had to
used to validate the functionality of CobWeb 1.0. A
undergo a data post-processing routine to enhance the
total of 12 region of interests (ROIs) were thus
image quality. Details of the image enhancement are
investigated to determine suitable REVs. Figure 4
published in Sell et al. (9). In brief, several image
shows the ROIs selected for the Berea, Grosmont and
enhancement techniques were tested in preliminary
gas hydrate samples. Image pre-processing is one of
studies including filters and filter combinations to
the essential and precautionary steps before image
gain best-fit results for further numerical simulations.
segmentation (21, 22). Image filtering techniques
To our knowledge, the combination of the non-local
such as blur, background intensity variation and
means and the anisotropic diffusion filter, both
contrast help in reducing artefacts. Image denoising
implemented in Avizo (ThermoFisher Scientific),
filter such as median filter, non-local means filter, and
works best for the given GH data.
anisotropic diffusion filter can assist in lowering the
phase misclassification and improving the The concept of the anisotropic diffusion (AD)
convergence rate of automatic segmentation schemes. filter is to smooth out noise in predefined areas of an
In the case of XCT volume stack of Berea sandstone, image, but also to stop at sharp edges representing
the 3D reconstructed raw images (10243) had boundaries between phases. This way, edges and
sufficiently high resolution and contrast, thus they did sharp boundaries between phases are preserved, and
not show any noticeable change on using the above image noise is significantly reduced (23, 24). A
mentioned filters. However, with XCT images comparison of the current voxel with its six neighbors
(10243) of the Grosmont carbonate rock, non-local takes place, and diffusion is fulfilled when the
means filtering yielded better visualisation and threshold stop criterion is not exceeded. If the
performance results compared to anisotropic difference between one voxel and its six adjacent
diffusion filter. neighbors exceeds the given value no diffusion takes
place. Another option to control the diffusion process
of the filter is to reduce or increase the diffusion time.
The parameter number of iterations defines how often
the algorithm will be used on the data. The bigger this
 9

number, the more blurred the resulting image. The supervised ML techniques rely on features
Smoothing is performed by applying a Gaussian also termed as feature vectors (FVs). The FVs are sets
filter. For our investigations, the threshold stop of instances that represent descriptive information on
criterion was set to the value 22,968 as this is the which ML algorithm is used to train the classification
approximated transition of the grain phase to hydrate. model. They further identify these features in an
AD was run on a CPU device with five iterations. unknown dataset and group them into respective
classes. Least square support vector machine
The non-local means filter (NLM) is a windowed
(LSSVM) is one such supervised ML technique,
version of the non-local means algorithm (25, 26).
which in recent years has emerged as a reliable
The main aim is to de-noise data based on comparing
technique to segment digital rocks images (7). Khan
voxels for similarities in a selected window in which
et al. (30) provides concise description and
a new weight for a voxel is assigned. After a Gauss
MATLAB® code snippet for the implementation of
kernel was run on the weighted values, the new value
the LSSVM library on XCT images, whereas
will be assigned replacing the former grey values. The
Chauhan et al. (8) validated its best performance and
filter is most efficient if the image is affected by white
accuracy in comparison to other common ML
noise. In Avizo the parameter window size, the local
techniques. In practice, an FV is a group containing
neighborhood, and the similarity value can be
subsets of different pixel values. For example, the FV
customized. Furthermore, the NLM filter is also an
of class four is a group encapsulating pixel values
appropriate tool for salt-and-pepper de-noising
corresponding to the pore, matrix rock, and noise. The
caused by image sensor defects (27). For this study,
pixel values were selected from a single 2D slice
the NLM filter was run in 3D mode on a CPU device.
representative of the REV. This FV was used for
The search window was set to 21 and the local
training the classification model. The training
neighborhood to 6 at a similarity value of 0.71.
performance was monitored using a 10 K-fold cross-
4.2 Phase Segmentation validation technique (31-33).

4.2.1 Grosmont Carbonate and Berea Sandstone 4.2.2 Gas Hydrate-Bearing Sediments ― Dual
Rock Samples Clustering

The K-means algorithm was used for the The edge enhancement (ED) effect was
segmentation of REV analysis of Berea and Grosmont significant in all reconstructed slices of the GH
rocks. K-means is one of the simplest unsupervised dataset. The ED effect was seen around the quartz
machine learning (ML) algorithms commonly used to grains mostly, with high and low pixel intensities
address clustering problems (28, 29, and 7). The K- adjacent to each other. The high intensity pixel (EDH)
means algorithm iteratively calculates the Euclidean values were very close to GH pixel values, while the
distance between the data points (pixel value) to its low intensity pixel (EDL) values showed a variance
nearest centroid (cluster). The algorithm converges between noise and brine phase pixel values.
when the objective function, i.e. the mean square root Therefore, immediate segmentation performed on the
error of Euclidean distance, reaches the minimum. pre-filtered GH datasets using CobWeb 1.0 resulted
This is when no further pixel is left to be assigned to in misclassification. Further parameterising and
the nearest centroid (cluster). However, the K-means tuning the unsupervised (K-means) and supervised
algorithm has the tendency to terminate without (LSSVM) modules of CobWeb 1.0 specifically,
identifying the global minimum of the objective distance function (i.e., functions euclidean distance
function. Therefore, running the algorithm several sqeuclidean, sum of absolute differences cityblock,
times is recommended to increase the likelihood that and mandist) and different permutation and
the global minimum of the objective function will be combination between kernel type, bandwidth and
identified. The performance of the K-means cross-validation parameters, showed significant
algorithm is strongly governed by the initial choice of improvement, but the segmentation was still not
the cluster centres (7). optimal. The aim was to eliminate the ED features
 10

completely without altering the phase distribution  Step 6: K-means clustering

between GH and the brine. This prompted K-means segmentation with three classes was
development of a GH-specific workflow as explained performed on the rescaled dataset to obtain the
below. The appendix provides the MATLAB® script final result.
for this workflow comprised of 6 steps:
4.3 Representative Elementary Volume Analysis
 Step 1: Filtering and REV selection
Four REVs of size 4 × 7003 were cropped The representative elementary volume (REV)
from the raw (16 bit) data stack. These REVs can be defined as the smallest volume, which should
were dual-filtered using NLM and ASD filters ideally represent the average effective macroscopic
(see section 4.1.1). Figure 6 shows the pre- behavior of the geomaterial. As a result, the transport
filtered raw dataset. of the effective parameters (mass, momentum,
 Step 2: K-means clustering energy) mathematically modelled within the REV
In this step, K-means segmentation is becomes independent of the sample size (34).
performed on the REVs to label the phases Figure 5 schematically explains the relationship
into different classes. The class sizes were between porosity and the volume of the porous media.
varied between three to twenty, and it was thus In a small REV (region I), high fluctuation in porosity
established that class seven captured all the is caused by heterogeneity at the pore scale. As
desired phases (noise, edge enhancement low volume increases (region II), porosity starts to
intensities (EDL), brine, quartz, edge normalize above some Vmin value within a small
enhancement high intensities (EDH), GH). standard deviation around a constant value of
 Step 3: Indexing porosity. The porosity measured in this region is
The pixel indices corresponding to the scale-independent, and an accurate representation of
respective classes (desired phases) were a large-scale system. The increase in REV value
extracted from the segmented slice(s). above a Vmax may result in increase/decrease in
Thereafter, using these pixel indices as porosity related to increases in heterogeneity,
reference, corresponding pixel values were associated with `macroscopic’ volume features
extracted from the 16-bit raw images. The (region III) (34). For heterogeneous porous media,
obtained pixel values represent noise, EDL, porosity theoretically lies in between region I and
brine, quartz, EDH, and GH phases in the raw region III depending on the effective parameter under
images. Then, histogram distribution of the investigation; however, the determination of ideal
pixel values in each phase was plotted. The region II for a real heterogeneous system may be
skewness of the histogram plots was inspected difficult and subjective (2-6).
visually and mean and standard deviation for In particular while performing permeability
each histogram were calculated. If an overlap tensor simulation using XCT data, the size of the
of pixel intensities was found in the different minimum REV should be assessed based not only on
histograms (phases), step four was repeated. porosity but also on geometrical parameters such as
 Step 5: Rescaling raw REV pore size distribution, void ratio, and coordinate
The pixel values corresponding to the phases, number (5, 6). For this study, we investigated
i.e brine, quartz, and GH, were replaced by different ROIs and REV sizes between 3003 to
their mean values, with an exception for EDH 5003 resolution, and established that an REV of size
pixel values. The latter were replaced with the 471 x 478 x 480 suited best. Figure 7 shows the REVs
mean value of quartz. These assignments lead of Berea sandstone and Grosmont carbonate rock and
to optimal segregation of the phase boundaries their respective geometrical parameters.
in the raw dataset and finally to the
elimination of the ED effect.
 11

Figure 4. The most suitable REVs of Berea sandstone and Grosmont carbonate rock shown in panel and gas hydrate-bearing sediment datasets shown
in the panel b.

Figure 5. Schematic representation of the relationship between porosity () and volume (V) of porous media. Bachmat and Bear (1986).
 12

Figure 6. 2D filtered, rescaled, and segmented slices of gas hydrate REV1 dataset
 13

Figure 7. Top panel shows surface plot of REVs Berea sandstone and Grosmont carbonate (size 471x478x480) using visualisation software ParaView.
Middle plot shows the relative porosity (%) trend for Berea sandstone and Grosmont carbonate REVs samples. Bottom plot shows the pore size
distribution of Berea sandstone and Grosmont carbonate. XCT images were segmented using K-means. In the case of Grosmont, a non-local means
filter was used
 14

Figure 8. The top panel shows relative porosity trend analysis of gas hydrates, the middle and bottom panel show the geometrical pore size
distribution of the respective REVs. The analysis was performed using CobWeb 1.0.

calculated based on the morphological scheme

4.4 Estimation of Relative Porosity and Pore Size suggested in Rabbani et al. (36). The mean relative
Distribution porosity value of Berea sandstone is 17.3 ± 2.6 %,
In the case of the Grosmont carbonate and Berea whereas for Grosmont carbonates the mean porosity
sandstone, the respective REVs were segmented value is lower (10.5 ± 2.3 %), as shown in Figure
using K-means and LSSVM, and the PSD was 7. The regression coefficient value of R2 = 0.092 for
the Berea sandstone porosity trend indicates that
 15

porosity remains constant throughout the REV sizes independent heterogeneities. However, there is high
chosen and is therefore consolidated for scale- variance compared with the mean PSD values. The
independent heterogeneities. In the case of the exact reason is unknown, but it may be due to the
Grosmont carbonate rock, the chosen REV size was drastic increase and decrease of the quartz grains, as
the best found out of five others explored, which can be noticed in Figure 6. The first and last 2D slices
consolidate again for scale-independent of ROI 1 in show either non-isotropic or isotropic
heterogeneities. The average pore size distribution distribution of quartz grains, which might have
thus obtained was 6.70 μm ± 0.68 μm and contributed to the respective high and low standard
14.21 μm ± 0.66 μm for Berea and Grosmont plug deviation seen in the porosity distribution. Figure 9
samples, respectively. shows the surface and volume rendered plots of REV
1 and REV 2. Due to the high accuracy of
Similarly, the porosity and PSD of the four GH
segmentation, the quartz grain, brine and GH
REVs were analyzed using CobWeb 1.0 except for
boundaries are clearly segregated and the ED effect
segmentation, which was performed using a different
is completely eliminated.
workflow as discussed above. Figure 8 shows the
comparison of the porosity trends of different GH
REVs. The selected REVs consolidate for the scale

Figure 9. Segmented REVs of gas hydrate sample displayed as surface and volume rendered plots. as analyzed using CobWeb 1.0 and exported to
VTK format using CobWeb 1.0 ParaView plug-in. Quartz grain phase is represented in green color, gas hydrate in red, and in blue the liquid brine
phase.
 16

segment data without compromising data quality at a

5 Conclusions and Outlook
reasonable speed makes it a favourable tool for REV
With CobWeb 1.0, this paper introduces a analysis.
visualisation and image analysis toolkit dedicated to
However, CobWeb 1.0 is somewhat limited
representative elementary volume analysis of digital
regarding its volume rendering capabilities, which
rocks. It was developed on the MATLAB®
will be one of the features to improve with the next
framework and can be used as a MATLAB® plugin or
version. The volume rendering algorithms
as a standalone executable. It relies on robust image
implemented thus far are not as sophisticated as those
segmentation schemes based on machine learning
in ParaView or DSI studio codes, which rely on the
(ML) techniques, which can be tested and cross-
OpenGL marching cube scheme. At present, the
validated in parallel. Dedicated image preprocessing
densely nested loop structure appears to be the best
filters such as the non-local means, anisotropic
choice for systematic processing. However, in future
diffusion, averaging and the contrast enhancement
versions vectorisation and indexing approaches
functions help to reduce artefacts and increase the
(bsxfun, repmat) will have to be tried and a check for
signal to noise ratio. Petrophysical and geometrical
significant changes in processing speed.
properties such as porosity, pore size distribution and
MATLAB®―Java synchronisation will be explored
volume fractions can be computed with ease on a
further to configure issues related to multi-threading
single representative 2D slice of a complete REV 3D
and visualisation (Java OpenGL).
stack. This was tested further using synchrotron
datasets of the Berea Sandstone, a gas hydrate- In the science segments, the file readers and
bearing sediment and a tomography dataset of the subroutines will be improved to analyse and overlay
Grosmont carbonate rock. The gas hydrate dataset, scanning electron microscope data with XCT data to
despite it’s nanoscale resolution, was infested with enhance mineral identification. A module CrackNet
strong edge enhancement (ED) artefacts, which (crack network) is planned which will explicitly
causes discrepancies in different modelling tackle segmentation of cracks, fissures in
approaches. A combination of dual filtering and dual geomaterials using machine learning techniques, and
clustering approaches is proposed to completely a mesh generation plugin (stl format) for 3D printing.
eliminate the ED effect in the gas hydrate sediments Pore network extraction and skeletonisation schemes
and the code is attached as an appendix. The REV such as modified maximum ball algorithm (37) and
studies performed on Berea sandstone, Grosmont medial axis transformation (38) will be considered
carbonate rock and GH sediment using CobWeb1.0 such that the data can be exported to an open-source
shows relative porosity trends with very low linear pore network modelling package (39).
regression values of 0.092, 0.1404, 0.0527
respectively. CobWeb 1.0’s ability to accurately

6 Acknowledgements
We thank Heiko Andrä and his team at Fraunhofer ITWM, Kaiserslautern, Germany, for providing us with the
synchrotron tomography dataset of the Berea sandstone. This study was funded within the framework of the
SUGAR (Submarine Gashydrat Ressourcen) III project by the Germany Federal Ministry of Education and
Research (BMBF grant 03SX38IH). The sole responsibility of the paper lies with the authors.
 17

7 References
1. Leu, L., Berg, S., Enzmann, F., Armstrong, R.T., and Kersten, M., (2014). Fast X-ray micro-tomography
of multi-phase flow in Berea sandstone: A sensitivity study on image processing. Transp. Porous Media
105, 451-469.
2. Zhang, D., Zhang, R., Chen, S., and Soll, W. E., (2000). Pore scale study of flow in porous media: Scale
dependency, REV, and statistical REV. Geophysical Research Letters, 27, 1195–1198.
3. Gitman, I. M., Gitman, M. B., and Askes, H., (2006). Quantification of stochastically stable representative
volumes for random heterogeneous materials. Archive of Applied Mechanics 75(2–3), 79–92
4. Razavi, M. R., Muhunthan, B., and Al Hattamleh, O., (2007). Representative elementary volume analysis
of sands using X-ray computed tomography. Geotechnical Testing Journal, 30, 212–219.
5. Al-Raoush, R., and Papadopoulos, A., (2010). Representative elementary volume analysis of porous
media using X-ray computed tomography. Powder Technology 200, 69–77.
6. Costanza-Robinson, M. S., Estabrook, B. D., and Fouhey, D. F., (2011). Representative elementary
volume estimation for porosity, moisture saturation, and air–water interfacial areas in unsaturated porous
media: Data quality implications. Water Resources Research, 47, W07513.
7. Chauhan, S., Rühaak, W., Khan, F., Enzmann, F., Mielke, P., Kersten, M., Sass, I., (2016a). Processing
of rock core microtomography images: Using seven different machine learning algorithms, Computers &
Geosciences, 86, 120-128.
8. Chauhan, S., Rühaak, W., Anbergen, H., Kabdenov, A., Freise, M., Wille, T., and Sass, I., (2016b). Phase
segmentation of X-ray computer tomography rock images using machine learning techniques: an accuracy
and performance study, Solid Earth, 7, 1125-1139, doi:10.5194/se-7-1125-2016.
9. Sell, K., Saenger, E.H., Falenty, A., Chaouachi, M., Enzmann, F., Kuhs, W.F., and Kersten, M., (2016).
On the path to the digital rock physics of gas hydrate-bearing sediments – processing of in situ
synchrotron-tomography data. Solid Earth, 7, 1243-1258.
10. Altman, Y., (2016). Blog entry: https://undocumentedmatlab.com/blog/uitable-customization-report. Last
check: 03/08/2011.
11. Chaouachi, M., Falenty, A., Sell, K., Enzmann, F., Kersten, M., Haberthür, D., and Kuhs, W. F., (2015).
Microstructural evolution of gas hydrates in sedimentary matrices observed with synchrotron x-ray
computed tomographic microscopy, Geochem. Geophy. Geosy., 16, 1711–1722.
12. Falenty, A., Chaouachi, M., Neher, S. H., Sell, K., Schwarz, J.-O., Wolf, M., Enzmann, F., Kersten, M.,
Haberthur, D., and Kuhs, W.F., (2015). Stop-and-go in situ tomography of dynamic processes – gas
hydrate formation in sedimentary matrices, Acta Crystallogr. A, 71, s154,
doi:10.1107/S2053273315097740.
13. Stampanoni, M., Groso, A., Isenegger, A., Mikuljan, G., Chen, Q., Bertrand, A., Henein, S., Betemps, R.,
Frommherz, U., Böhler, P., Meister, D., Lange, M., and Abela, R., (2006). Trends in synchrotron-based
tomographic imaging: the SLS experience. In: Proceedings of SPIE, The International Society for Optical
Engineering, 6318, pp. M1–M14.
14. Marone, F., and Stampanoni, M., (2012). Regridding reconstruction algorithm for real-time tomographic
imaging, J. Synchrotron Radiat., 19, 1029–1037.
15. Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M.,
Madonna, C., Marsh, M., Mukerji, T., Saenger, E.H., Sain, R., Saxena, N., Ricker, S., Wiegmann, A.,
Zhan, X., (2013): Digital rock physics benchmarks-Part I: Imaging and segmentation. Computers &
Geosciences, 50(0), 25-32.
16. Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M.,
Madonna, C., Marsh, M., Mukerji, T., Saenger, E.H., Sain, R., Saxena, N., Ricker, S., Wiegmann, A.,
Zhan, X., 2013. Digital rock physics benchmarks-part II: Computing effective properties. Computers &
Geosciences, 50(0), 33-43.
 18

17. Machel, H.G. and Hunter X., (1994). Facies models for middle to late Devonian shallow-marine
carbonates, with comparisons to modern reefs: a guide for facies analysis. Facies 30: 155.
doi:10.1007/BF02536895.
18. Buschkuehle, B.E., Hein, F.J., Grobe, M., (2007). An overview of the geology of the Upper Devonian
Grosmont carbonate bitumen deposit, Northern Alberta, Canada. Natural Resources Research, 16, 3–15.
19. Jovanović, Z., Khan, F., Enzmann, F., and Kersten, M., (2013). Simultaneous segmentation and beam-
hardening correction in computed microtomography of rock cores. Computers. Geosciences. 56, 142-15
20. Giesche, H., (2006). Mercury porosimetry: a general (practical) overview. Particle & Particle Systems
Characterization, 23 (1), 9–19.
21. Iassonov, P., Gebrenegus, T., and Tuller, M., (2009). Segmentation of X-ray computed tomography
images of porous materials: a crucial step for characterization and quantitative analysis of pore structures.
Water Resources Research, 45, W09415.
22. Schlüter, S., Sheppard, A., Brown, K., and Wildenschild, D., (2014). Image processing of multiphase
images obtained via X-ray microtomography: A review. Water Resources Research, 50, 3615-3639
23. Kaestner, A., Lehmann, E., and Stampanoni, M., (2008). Imaging and image processing in porous media
research, Adv. Water Resour., 31, 1174–1187.
24. Porter, M. L. and Wildenschild, D., (2010). Image analysis algorithms for estimating porous media
multiphase flow variables from computed microtomography data: A validation study, Computat.Geosci.,
14, 15–30.
25. Buades, A., Coll, B., and Morel, J.-M., (2005a). A non-local algorithm for image denoising, IEEE
Computer Society Conference on Computer Vision and Pattern Recognition, San Diego, CA, USA, 60–
65.
26. Buades, A., Coll, B., and Morel, J. M. (2005b): A review of image denoising algorithms, with a new one,
Multiscale Model. Sim. 4, 490–530.
27. Sarker, S., Chowdhury, S., Laha, S., and Dey, D., (2012). Use of non-local means filter to denoise image
corrupted by salt and pepper noise, SIPI, J Signal & Image Processing: An International Journal, 223–
235.
28. MacQueen, J., (1967). Some Methods for classification and Analysis of Multivariate Observations.
University of California Press, 281-297 pp.
29. Jain, A.K., (2010). Data clustering: 50 years beyond K-means, Pattern Recognition Letters, 31(8), 651-
666.
30. Khan, F., Enzmann, F., and Kersten, M., (2016). Multi-phase classification by a least-squares support
vector machine approach in tomography images of geological samples, Solid Earth, 7, 481-492,
doi:10.5194/se-7-481-2016.
31. Larson, S. C., (1931). The shrinkage of the coefficient of multiple correlation. Journal of Educational
Psychology, 22(1), 45.
32. Kohavi R., (1995). A study of cross-validation and bootstrap for accuracy estimation and model selection.
In Proceedings of International Joint Conference on AI. pp. 1137–1145.
33. Dietterich, T.G., (1998). Approximate statistical tests for comparing supervised classification learning
algorithms. Neural computation, 10(7), 1895-1923.
34. Wiącek, J, and Molenda, M., (2016). Representative elementary volume analysis of polydisperse granular
packings using discrete element method. Particuology, 27:88-94.
35. Bachmat, Y., and Bear, J., (1986). Microscopic Modelling of Transport Phenomena in Porous Media. 1:
The Continuum Approach. Transport in Porous Media 1, 1986, 213-240.
36. Rabbani, A., Jamshidi, S., and Salehi S., (2014). An automated simple algorithm for realistic pore network
extraction from micro-tomography images, Journal of Petroleum Science and Engineering, 2014, Volume
123, Pages 164-171, ISSN 0920-4105, https://doi.org/10.1016/j.petrol.2014.08.020.
 19

37. Arand, F., Hesser, J., (2017). Accurate and efficient maximal ball algorithm for pore network extraction,
In Computers & Geosciences, Volume 101, 2017, Pages 28-37.
38. Katz, R.A & Pizer, S.M, (2003). International Journal of Compter Vision 55: 139.
http://doi.org/10.1023/A:1026183017197
39. Gostick, J., Aghighi, M. Hinebaugh, J., Tranter, T., Hoeh, M.A., Day, H., Spellacy, B., Sharqawy, M.,
Bazylak, A., Burns, A., Lehnert, W., and Putz, A., (2016). OpenPNM: A Pore Network Modeling Package.
Computing in Science & Engineering 18(4), p 60-74.

8 Appendix A: MATLAB snippet for removal for Edge Enhancement Effect in gas hydrate
datasets

8.1 Gas Hydrate Segmentation

8.2 Step 1
The Dual Clustering approach, by which first the 16bit gas hydrate was filtered using Anisotropic
diffusion (ASD) filter, and then with non-local means (NLM) filter, to minimize/normalize the
edge enhancement (ED) artefacts.

8.3 Step 2
 read slice by slice 3D prefiltered raw data
 for this example the reading is restricted
 to only four slices (700x700x4); it can be changed using nZ variable

mfname= 'Xe_17w_8_ROI4_ADS_NLM';
ifname=[mfname,'.raw'];
nX=700;
nY=700;
nZ=4;
ldim = 1;
xDi=[nX nY]';
grenzwert=0;
clusterS =7;
ifid=fopen(ifname,'r');
M=zeros(nX,nY,nZ,'uint16');
SeData = nZ-ldim;
SeData = 1:1:SeData;
dim = size(M);

for k=ldim:nZ
disp(sprintf('Reading slice no. %d....',k));
s=sprintf('Slice: % d', k');
S=fread(ifid, [xDi(1) xDi(2)], 'uint16');
M(:,:,k)=S;
%figure; imagesc(M(:,:,k)); colorbar;
end

Reading slice no. 1....

Reading slice no. 2....
Reading slice no. 3....
Reading slice no. 4....

8.4 Display image

figure; imagesc(M(:,:,1)); colorbar;

title('Raw Prefiltered ASD NLM')
 20

8.5 Concatenate raw data into single array

%* concatenate array will be used in step three

M = M(:,:,ldim:nZ);
rawM = double(M(:));

8.6 Perform k-means clustering

Here, clustering is restricted to class 7 optimal to enable clustering of all the available features:
for ii = 1:2
R=double(M(:,:,ii));
[r,c,v]=find(R>grenzwert);
cyl=R>grenzwert;
R1=cyl.*R;
[m, n, w]=find(R1);
G = kmeansK(w,clusterS);
S=sparse(r,c,G,size(R,1),size(R,2));
M_seg=full(S);
SegImg(:,:,ii)=M_seg;
%figure; imagesc(SegImg(:,:,ii)); colormap(parula(5)); colorbar;
%title('K-means prefiltered');
end
 21

8.7 Display image

figure; h = imagesc(SegImg(:,:,1)); colormap(jet(max(h.CData(:))));

title('K-means prefiltred');

8.8 Step three

The purpose is to |index out| pixel values of different phases:
% noise
% edge enhanced low (EDL)
% liquid
% quartz
% edge enhanced high (EDL)
% gas hydrate
% from the concatenated raw images matrix using segmented class values
% thereafter compare their histogram % _as sanity check_
% to identify if any overlapping boundaries

8.9 Index noise pixels

rangeNl = 0;
indN = find(h.CData(:)==rangeNl);
rawO = rawM(indN);

8.10 Plot histogram noise

[cN, countN] = hist(rawO, 10);

%figure; bar(countN, cN);
%title('noise')
 22

8.11 Index EDL pixels

rangeNu = 2;
indD = find(h.CData(:)>rangeNl & h.CData(:)<=rangeNu);
rawD = rawM(indD);

8.12 Plot histogram noise

[cD, countD] = hist(rawD, 100);

%figure; bar(countD, cD);
%title('Edge Enhanced low noise')

8.13 Index liquid pixels

rangeLl = 1;
rangeLu = 3;
indL = find(h.CData(:)>=rangeLl & h.CData(:)<=rangeLu);
if min(SegImg(indL))==rangeLl & max(SegImg(indL))==rangeLu
rawL = rawM(indL);
min_rawL = min(rawL);
max_rawL = max(rawL);
Avg_rawL = mean(rawL);
else
fprintf('min and max for liquid dont match.....\n')
return
end

8.14 Plot histogram liquid

[cL, countL] = hist(rawL, 100);

figure; bar(countL, cL);
title('Liquid')
 23

8.15 Index EDH pixels

rangeE = 5;
indE = find(h.CData(:)==rangeE);
if min(SegImg(indE))==rangeE
rawE = rawM(indE);
min_rawE = min(rawE);
max_rawE = max(rawE);
Avg_rawE = mean(rawE);
else
fprintf('min and max for EDH dont match.....\n')
return
end

8.16 Plot histogram EDH

[cE, countE] = hist(rawE, 10);

figure; bar(countE, cE);
title('Edge Enhanced high noise')
 24

8.17 Quartz index phases

rangeQu = 4;

8.18 Quartz

indQ = find(h.CData(:)==rangeQu);
if min(SegImg(indQ)) == rangeQu
rawQ = rawM(indQ);
min_rawQ = min(rawQ);
max_rawQ = max(rawQ);
Avg_rawQ = mean(rawQ);
else
fprintf('min and max for quartz dont match.....\n')
return
end
%indQ = find(h.CData(:)>=rangeQl & h.CData(:)<=rangeQu);

8.19 Plot histogram quartz

[cQ, countQ] = hist(rawQ, 100);

figure; bar(countQ, cQ);
title('Quartz')
 25

8.20 Gas Hydrate

rangeMl =6;
rangeMu =7;
%indM =find(h.CData(:)>=rangeMu);
indM = find(h.CData(:)>=rangeMl & h.CData(:)<=rangeMu);
if min(SegImg(indM))==rangeMl & max(SegImg(indM))==rangeMu
rawMu = rawM(indM);
min_rawMu = min(rawMu);
max_rawMu = max(rawMu);
Avg_rawMu = mean(rawMu);
elseif min(SegImg(indM))==rangeMu & max(SegImg(indM))==rangeMu
rawMu = rawM(indM);
min_rawMu = min(rawMu);
max_rawMu = max(rawMu);
Avg_rawMu = mean(rawMu);
else
fprintf('min and max for gas hydrate dont match.....\n')
return
end

8.21 Plot Histogram Gas Hydrate

[cM, countM] = hist(rawMu, 100);

figure; bar(countM, cM);
title('Methane')
 26

8.22 Step 4 - rescaling the raw images

First min-max and mean of respective phases are calculated for the respective (raw) phases
(obtained above), which thereafter are replaced by their respective mean values.
%average values

8.23 With an exception to GH dataset

Where EDH (raw pixels) are replaced with averaged quartz values:
%as they are in close vicinity to quartz pixel values
M_replace = M(:);
min_li = min_rawL;
max_li =max_rawL;
avg_li =Avg_rawL;
min_Qz = min_rawQ;
max_Qz = max_rawQ;
avg_Qz = Avg_rawQ;
min_EDH = min_rawE;
max_EDH = max_rawE;
avg_EDH = Avg_rawE;
min_GH = min_rawMu;
max_GH = max_rawMu;
avg_GH = Avg_rawMu;
%--------------------------------------------------------------------------
% indxes of liquid pixels
%--------------------------------------------------------------------------
Ind_rep_L = find(M_replace>=min_li & M_replace <= max_li);
% replacement by average liquid value
if min(M_replace(Ind_rep_L))==min_li & max(M_replace(Ind_rep_L)==max_li)
M_replace(Ind_rep_L)=avg_li;
else
fprintf('min and max for liquid dont match.....\n')
return
end
%--------------------------------------------------------------------------
% indxes of quartz pixels
%--------------------------------------------------------------------------
 27

Ind_rep_Q = find(M_replace>= min_Qz & M_replace<= max_Qz);

% replacement by average quartz value
if min(M_replace(Ind_rep_Q))==min_Qz & max(M_replace(Ind_rep_Q))==max_Qz
M_replace(Ind_rep_Q)= avg_Qz;
else
fprintf('min and max for quartz dont match.....\n')
return
end
%--------------------------------------------------------------------------
% indxes of EDH pxels
%--------------------------------------------------------------------------
Ind_rep_E = find(M_replace>= min_EDH & M_replace<=max_EDH);
%

Replace by average quartz values

if min(M_replace(Ind_rep_E))==min_EDH & max(M_replace(Ind_rep_E))==max_EDH
M_replace(Ind_rep_E)= avg_Qz;
else
fprintf('min and max for EDH dont match.....\n')
return
end
%--------------------------------------------------------------------------
%indexes of gas hydrate pixels
%--------------------------------------------------------------------------
Ind_rep_M = find(M_replace>=min_GH & M_replace<=max_GH);
% replacement by average gas hydrate value
if min(M_replace(Ind_rep_M))== min_GH & max(M_replace(Ind_rep_M))==max_GH
M_replace(Ind_rep_M)= avg_GH;
else
fprintf('min and max for methane dont match.....\n')
return
end
%--------------------------------------------------------------------------

8.24 Reshape rescaled array

%--------------------------------------------------------------------------
M_replaced = reshape(M_replace,[dim(1), dim(2), dim(3)]);
clear M_replace;
figure; imagesc(M_replaced(:,:,1));
title('Rescaled Raw');
 28

8.25 Step 5
K-means clustering is performed on the rescaled images to obtain segmetation in three classes:
clusterS =3;
initialcenters = [avg_li,avg_Qz,avg_GH];
for ii = 1:dim(3)
R=double(M_replaced(:,:,ii));
[r,c,v]=find(R>grenzwert);
cyl=R>grenzwert;
R1=cyl.*R;
[m, n, w]=find(R1);
G = kmeans(w,clusterS,'Distance','sqeuclidean','start',initialcenters');
S=sparse(r,c,G,size(R,1),size(R,2));
M_seg=full(S);
SegImg(:,:,ii)=M_seg;
%figure; imagesc(SegImg(:,:,ii)); colormap(parula(5)); colorbar;
end
figure; imagesc(SegImg(:,:,1)); colormap(parula(5)); colorbar;
title('Gas Hydrate ROI');
fclose('all');
 29

Published with MATLAB® R2015b