# Development and operation of the CGEM Inner Tracker for the BESIII
experiment

Alberto Bortone  on behalf of CGEM-IT Working Group

###### Abstract

The extension of data acquisition for the BEijing Spectrometer (BESIII)
experiment until at least 2030 has resulted in upgrades to both the
accelerator and the detector. An innovative Cylindrical Gas Electron
Multiplier (CGEM) is under construction to upgrade the inner tracker, which is
suffering from aging. The CGEM Inner Tracker (CGEM -IT) was designed to
restore efficiency and enhance the reconstruction of the secondary vertexes
position. Reconstruction in the magnetic field of 1 T requires an analog
readout and an electronic contribution to the time resolution better than 5
ns. The entire system consists of about 10,000 electronic channels and must
maintain a peak rate of 14 kHz/strip of signal hits for the innermost layer of
the CGEM-IT. The CGEM readout system is based on the innovative TIGER ASIC,
which is manufactured using 110 nm CMOS technology. A special readout chain
consisting of GEM Read Out Cards (GEMROC) was developed for data acquisition.
Two of three layers, equipped with the final electronics, have been operating
in Beijing since January 2020 collecting cosmic ray data in a test facility,
remotely controlled by the Italian group due to the pandemic situation. Since
it was not possible to perform further tests on the cylindrical detector, in
July 2021, a test beam was conducted at CERN with the final electronics
configuration on a small prototype consisting of four triple-GEM planar
detectors. In this presentation, the general status of the project CGEM-IT
will be presented, with particular emphasis on the results of the test beam
data collection.

###### keywords:

keyword one , keyword two

###### PACS:

0000 , 1111

###### MSC:

0000 , 1111

††journal: Nuclear Instruments and Methods in Physics Research A

[a]organization=Turin University, department of Physics,addressline=Via Giuria
1, city=Turin, state=Italy [b]organization=Istituto Nazionale di Fisica
Nucleare, Turin section,addressline=Via Giuria 1, city=Turin, state=Italy

##  1 Introduction

The BEijing Spectrometer (BESIII) is an experimental facility for high-energy
physics at the Beijing Electron Positron Collider (BEPCII), with an extensive
physics program that includes charmonium and charmonium-like states, charmed
mesons and baryons, light hadron spectroscopy, $\tau$ physics, QCD and CKM
parameters, baryon form factors, and new physics by studying rare and
forbidden decays [1].
BEPCII is a two-ring $e^{+}$ $e^{-}$ collider operating in the $\tau$-charm
region with a center of mass energy between 2 and $4.95\text{\,}\mathrm{GeV}$
and a luminosity of up to
${10}^{33}\text{\,}\mathrm{c}\mathrm{m}^{-2}\mathrm{s}^{-1}$[2].
The cylindrical core of the BESIII detector consists of multilayer drift
chamber (MDC) [3], as a charged particle tracker. This detector has been
taking data since 2009, but due to the high luminosity of the experiment, its
performance is degrading, with a gain loss per year of about
$4\text{\,}\mathrm{\char 37\relax}$ for the innermost layers, despite the
countermeasures take [4].
Since BESIII is expected to be in operation for the next ten years, an upgrade
is needed. An innovative solution to replace the aging inner MDC is to use a
CGEM Inner Tracker (CGEM-IT).

##  2 The CGEM-IT detector

The gas electron multiplier (GEM) is a well developed micro pattern gas
detector technology, with a wide range of applications in high-energy physics
[5].
In a triple-GEM detector, the multiplication structure takes place in a
cascade of three GEM foils, allowing a higher gain while keeping the discharge
risk low [6].

Figure 1: The CGEM-IT design

The CGEM-IT consists of three coaxial layers of triple GEM (figure 1). Each
cylindrical detector layer is independently assembled, has an autonomous gas
enclosure, and can be operated in stand-alone mode.
The position reconstruction benefits from the combination of two different
algorithms: the Charge Centroid extracts the average position by weighting the
signal amplitude of the firing strips. The µ-TPC (micro Time-Projection
Chamber) uses the drift gap like a time-projection chamber. The positions of
the primary ionizations in the drift gap are reconstructed by knowing the
drift velocity and the arrival time of the signal at the anode [7]. The
detector aims to achieve an r/$\Upphi$ resolution of about
$130\text{\,}\mathrm{\SIUnitSymbolMicro m}$, a resolution along the beam
direction of about $350\text{\,}\mathrm{\SIUnitSymbolMicro m}$, with a
material budget around 1.5% $X_{0}$ [8, 9].

##  3 The CGEM-IT readout

The generated signal is induced on two arrays of strips to obtain a 2D
readout: $570\text{\,}\mathrm{\SIUnitSymbolMicro m}$ wide $\Upphi$-strips,
parallel to the detector axis, and $130\text{\,}\mathrm{\SIUnitSymbolMicro m}$
wide V strips that have a stereo angle with respect to the $\Upphi$ strips.
The stereo angle is maximized according to the active area of each layer
(between 31.0° to 46.7°), while the pitch,
$650\text{\,}\mathrm{\SIUnitSymbolMicro m}$, is the same for both arrays and
all the layers.
The reconstruction algorithms, the foreseen performance and the BESIII
environment place special demands on the readout functions in terms of charge
measurement, time resolution, and sustainable rate, requiring the development
of a dedicated readout chain [10] .
The signal information is acquired by the front-end ASIC TIGER (Torino
Integrated GEM Electronics for Readout). The TIGER chips are mounted in pairs
on Front-End Boards (FEB) and installed on the detector. Data and ASIC Low
Voltage (LV) are fed through Data Low Voltage Patch Cards (DLVPC) by the GEM
Read Out Cards (GEMROC).
The GEMROC boards receive signals from the BESIII timing and trigger
interface, communicate with the BESIII slow control via the Ethernet interface
and via optical fibers with the GEM-Data Concentrators (GEM-DC). These cards
build the events and communicate with the VME-based BESIII DAQ. The GEMROC
boards also manage the power supply and configuration of the FEBs, as well as
the TIGER output data acquisition. The GUFI (Graphical Frontend User
Interface) software has been developed to control the electronics and manage
the data acquisition, while the CIVETTA (Complete Interactive VErsatile Test
Tool Analysis) software is used for on-run and off-line data validation and
analysis [11].

##  4 Beam test

A beam test was performed with planar triple-GEM detectors to investigate the
system integration and evaluate the performance of the triple GEM detector
equipped with the final TIGER + GEMROC readout chain.
The experiment was performed at the SPS North Area test beam facility at CERN,
on the H4 beam line, with a 150 GeV $\pi$ beam and a $80\text{\,}\mathrm{GeV}$
$\mu$ beam.
About 250M triggers have been acquired with beam incidence angles ranging from
0° to 45°, testing various detector and electronics settings.
The tested hodoscope consisted of four planar GEM detectors and two plastic
scintillators coupled to PMTs in order to provide the trigger. The system was
thoroughly tested by scanning different beam particle incidence angles, HV and
readout electronics settings. Figure 2 shows the preliminary evaluation of the
efficiency as a function of the GEM cumulative voltage. The efficiency
increases with GEM voltage until it reaches a plateau. In analyzing these
data, we discovered a GEMROC firmware issue that can reduce the efficiency. A
fix has been implemented and is currently being tested.

Figure 2: Detector efficiency as a function of the cumulative GEM voltage.

## Acknowledgment

We acknowledge the support of the BESIIICGEM Project (645664) - H2020-MSCA-
RISE-2014, and of the FEST Project (872901) - H2020-MSCA-RISE-2019, which
included INFN and IHEP as well as institutes from Mainz and Uppsala
Universities.

## References

  * Ablikim et al. [2020] M. Ablikim, et al. (BESIII),  Future Physics Programme of BESIII,  Chinese Physics C 44 (2020) 040001. doi:10.1088/1674-1137/44/4/040001.
  * BEPC II collaboration [2009] BEPC II collaboration,  BEPC II: construction and commissioning,  Chinese Physics C 33 (2009) 60–64. doi:10.1088/1674-1137/33/s2/016.
  * Ablikim et al. [2010] M. Ablikim, et al. (BESIII),  Design and Construction of the BESIII Detector,  Nucl. Instrum. Meth. A 614 (2010) 345–399. doi:10.1016/j.nima.2009.12.050. arXiv:0911.4960.
  * Dong et al. [2016] M. Dong, et al.,  Aging effect in the BESIII drift chamber,  Chinese Physics C 40 (2016) 016001. doi:10.1088/1674-1137/40/1/016001.
  * Sauli [2016] F. Sauli,  The gas electron multiplier ( GEM ): Operating principles and applications,  Nucl. Instrum. Meth. A 805 (2016) 2–24. doi:10.1016/j.nima.2015.07.060.
  * Bachmann et al. [2002] S. Bachmann, et al.,  Discharge studies and prevention in the gas electron multiplier (GEM),  Nucl. Instrum. Meth. A 479 (2002) 294–308. doi:10.1016/S0168-9002(01)00931-7.
  * Alexeev et al. [2016] M. Alexeev, et al.,  Development and test of a $\mu$TPC cluster reconstruction for a triple GEM detector in strong magnetic field,  in: 2016 IEEE (NSS/MIC/RTSD), 2016, pp. 1–4. doi:10.1109/NSSMIC.2016.8069914.
  * Marcello et al. [2018] S. Marcello, et al.,  A new inner tracker based on GEM detectors for the BESIII experiment,  International Journal of Modern Physics: Conference Series 48 (2018) 1860119. doi:10.1142/S2010194518601199.
  * Balossino et al. [2022] I. Balossino, et al.,  The CGEM-IT: An upgrade for the BESIII experiment,  Symmetry 14 (2022). doi:10.3390/sym14050905.
  * Amoroso et al. [2021] A. Amoroso, et al.,  The CGEM-IT readout chain,  JINST 16 (2021) P08065. doi:10.1088/1748-0221/16/08/p08065.
  * Bortone [2021] A. Bortone, Deployment of the readout electronics for the BESIII Cylindrical GEM Inner Tracker, Ph.D. thesis, UniTO, 2021.


∎

11institutetext: E.M. Mirkes 22institutetext: University of Leicester,
Department of Mathematics, Leicester, UK
22email<EMAIL_ADDRESS>
OrcID http://orcid.org/0000-0003-1474-1734 33institutetext: A.N. Gorban
44institutetext: University of Leicester, Department of Mathematics,
Leicester, UK
44email<EMAIL_ADDRESS>
OrcID http://orcid.org/0000-0001-6224-1430 55institutetext: J. Levesley
66institutetext: University of Leicester, Department of Mathematics,
Leicester, UK
66email<EMAIL_ADDRESS>
OrcID http://orcid.org/0000-0002-3509-0152 77institutetext: P.A.S. Elkington
88institutetext: Weatherford, Gotham Road, East Leake, Loughborough, LE12 6JX,
UK
Tel.: +44-115-9457880, Fax: +44-115-9405373
88email<EMAIL_ADDRESS>
OrcID http://orcid.org/0000-0001-8415-2260 99institutetext: J. A. Whetton
1010institutetext: Weatherford, Gotham Road, East Leake, Loughborough, LE12
6JX, UK
OrcID http://orcid.org/0000-0003-1527-9094

# Pseudo-Outcrop Visualization of Borehole Images and Core Scans

Evgeny M. Mirkes  Alexander N. Gorban  Jeremy Levesley  Peter A. S. Elkington
James A. Whetton

(Received: date / Accepted: date)

###### Abstract

A pseudo-outcrop visualization is demonstrated for borehole and full-diameter
rock core images to augment the ubiquitous unwrapped cylinder view and thereby
to assist non-specialist interpreters. The pseudo-outcrop visualization is
equivalent to a nonlinear projection of the image from borehole to earth frame
of reference that creates a solid volume sliced longitudinally to reveal two
or more faces in which the orientations of geological features indicate what
is observed in the subsurface. A proxy for grain size is used to modulate the
external dimensions of the plot to mimic profiles seen in real outcrops. The
volume is created from a mixture of geological boundary elements and texture,
the latter being the residue after the sum of boundary elements is subtracted
from the original data. In the case of measurements from wireline
microresistivity tools, whose circumferential coverage is substantially less
than 100%, the missing circumferential data is first inpainted using
multiscale directional transforms, which decompose the image into its
elemental building structures, before reconstructing the full image. The
pseudo-outcrop view enables direct observation of the angular relationships
between features and aids visual comparison between borehole and core images,
especially for the interested non-specialist.

###### Keywords:

Wellbore Microresistivity Image Inpainting

##  1 Introduction

A common starting point for the geological analysis of circumferential images
from boreholes and full-diameter cores is the identification of the types and
orientations of bed boundaries and other discontinuities (Rider and Kennedy,
2011). Analysts determine the orientation of individual features and examine
the ways in which the orientation of successive boundaries evolves over depth
to delineate local structures and to infer depositional environments. In
combination with other data at a variety of scales, image data helps refine
geological models and helps define flow units and baffles in reservoir models.
Image data also correlate with variations in mineralogy, porosity, and fluid
content, making them useful in petrophysical analysis, (Xu and Torres-Verdín,
2013; Fu et al., 2016). Additional textural information may be available from
optical surface scans of full diameter rock cores; correlating these with
borehole images in the same well allows the core and its features to be
oriented in the earth’s frame of reference via orientation data acquired by
the logging tool.

High-resolution borehole images are routinely recorded, yet they are largely
underutilized because of a lack of understanding of the data (Joubert et al.,
2016). A contributing factor is the way in which the data are visualized; note
that similar visualizations are applied to circumferential scans of full-
diameter core. Both types of images are commonly displayed as unwrapped
cylinders in which planar features manifest themselves as sinusoids, and
appropriate software is needed to compute their individual orientations. The
cylinder view has the advantage of being straightforward to generate, but it
requires a mental transformation to visualize the angular relationships that
would be observed in outcrops.

To observe boundary orientations in the earth’s frame of reference, a solid
volume corresponding to the material removed from the borehole must be
constructed and sliced vertically in two planes. Modulating the external
dimensions with a proxy for grain size gives a two-sided pseudo-outcrop view
in which planar features are seen to be planar (rather than sinusoidal), and
the angular relationships between them are observed directly without reference
to other plots. All the information available to the reconstruction is located
on the borehole wall or core surface; the interior has not been sampled. This
differs from the typical three-dimensional (3D) inpainting problem in which
the interior of an object is partially sampled. Aspect ratio is another
differentiating characteristic because borehole images are long relative to
their diameter. Image logs from a typical 0.2-m-diameter well commonly exceed
1,000 m in length; for the data sets examined here, the depth sample increment
is 2 mm and comprises 176 circumferential samples. It is desirable to
visualize data quickly.

In the case of wireline tool microresistivity measurements, current flows from
individual 4-mm-diameter electrodes arranged in rows and inlaid in pads
pressed against the borehole wall (Fig. 1(a)). In small-diameter wells, the
rows overlap circumferentially to give full coverage; but when the well
diameter exceeds approximately 150 mm (Fig. 1(c)), the space between pads
creates longitudinal gaps in the image (Fig. 1(b)). The gaps are bound by
linear margins in the tool’s frame of reference. Any 3D inpainting solution
must handle the gaps directly; alternatively, a two-dimensional (2D)
inpainting solution must be applied as a preprocessing step.

In general there are two classes of algorithms for filling holes in digital
images: texture synthesis and inpainting (Criminisi et al., 2004). The former
has been demonstrated for repeating patterns with some stochasticity and the
latter for linear and curvilinear structures such as object contours. Both
classes of algorithm have been applied to borehole images. Hurley and Zhang
(2011) use the statistical properties of the measured parts to generate
realizations for filling the gaps in microresistivity images. This method can
be appropriate for images dominated by textural features, such as vugs typical
of carbonate rock formations; but achieving continuity across gap boundaries
can be problematic, which makes it less well suited to images of rocks with
complex curvilinear structures. This limitation is addressed by constraining
the selection of realizations with structural information from locally adapted
kernel regression, (Takeda et al., 2007; Zhang et al., 2016). An inpainting
method based on compressed sensing has been shown to successfully reconstruct
full circumferential coverage in borehole images with up to 30% coverage loss
and up to 50% in favourable cases (Assous et al., 2014). The inpainting method
builds on the general idea of minimizing the data required to represent
information, which it seeks to do with the fewest number of multiscale
directional transform coefficients such that reversing the process recovers
the full data, including parts missing from the original samples (Elad et al.,
2005). The images are decomposed into elemental building structures (textures
and piecewise smooth parts) that are combined to represent observed features
(Starck et al., 2005). This method is well suited to images containing bedding
surfaces, fractures, slumps, clasts, and other curvilinear features typical of
clastic rock formations and (depending on the choice of directional transform)
is also able to handle textural elements.

Two-dimensional inpainting techniques have been extended to three dimensions
for a variety of applications. In particular, multipoint geostatistics has
been applied to filling gaps in micro-CT scans of porous rocks by extracting
characteristics from a training image to generate a database of
characteristics from which structures are selected to complete the missing
parts in the target image (Zhang and Du, 2014). Statistical methods have also
been applied to the creation of volumetric textures from 2D samples (Chen and
Wang, 2010; Urs et al., 2014). The compressed sensing approach has also been
extended to three dimensions; new multiscale directional transforms have been
designed that enable capturing the true geometry of 3Dl objects (avoiding the
risk of misrepresentation from slice-by-slice processing) (Woiselle et al.,
2011).

An alternative method is designed to be useful for images in which the primary
structural elements are planar or subplanar at the borehole/core scale. This
is the case in many hydrocarbon-bearing siliciclastic rocks logged by
microresistivity imaging tools. The method will propagate solid textures but
not necessarily in an optimal way. It can be easily combined with any advanced
method of texture analysis if textural elements are the primary interest (such
as in vuggy carbonates). It is tested on several real-life examples.

This paper comprises three sections: Description of the method, Results, and
Conclusion.

##  2 Description of the Method

The principal steps in constructing the 3Dl volume are as follows

  1. 1.

Divide the log into overlapping depth windows. In each window, identify the
dominant set of planar/subplanar geological features and remove them from the
image. This is Level 1 decomposition.

  2. 2.

Identify the dominant planar feature set within the residual image and
subtract it. This is Level 2 decomposition. Repeat the process several times
until all of the subsidiary planar geological features are identified and the
final residual image contains only non-planar features. The image may include
strongly non-planar features and a host of other features that are localized
azimuthally and will be referred to as residual texture.

  3. 3.

Identify functions to represent the shape of the primary geological feature
and each subsidiary feature for all the levels of decomposition, and use these
to construct a synthesised volume.

  4. 4.

Propagate the residual textures through the volume.

Once created, the volume is visualized by, for example, making two
longitudinal cuts through the oriented volume to expose a wedge of material,
which allows examination of the orientations of geological features and their
relationships in the earth’s frame of reference. Additional information from
the image measurement can be encoded as variations in the external dimensions
of the wedge.

Image of a borehole wall in several coordinate systems are presented in Fig.
2. Let us consider one of the overlapping depth windows as a cylinder (Fig.
2(c)). Before the earth frame transformation begins, the raw image-log data is
processed in the usual way and stored in an array $f_{ij}$, where
$i=1,\ldots,H$ are rows corresponding to depth and $j=1,\ldots,N$ are columns
corresponding to circumferential location (Fig. 2(a)). These 2D arrays
correspond to the surface of a cylinder (Fig. 2(c)). The borehole has a
nominal radius $R$ and the step between rows is $h$. Note that the data are
discrete with step $h$ in the vertical direction, and $f_{ij}$ is the
restriction of a scalar field $f$ to the surface of the borehole. Each
position $(i,j)$ in the array $f_{ij}$ corresponds to the space radius-vector
$\vec{r}_{ij}$ from the origin $O$ to the point on the cylinder surface.

Some of the $f_{ij}$ may correspond to missing or bad measurements, so a mask
$m_{ij}$ is used to identify the locations of valid values in which a value of
1 is used if a given $f_{ij}$ is known and 0 is used if $f_{ij}$ is unknown.
At this point it may be advantageous to construct a full circumferential
coverage image by application of the morphological components analysis method
(Assous et al., 2014).

The planar structure shown in Fig. 3 is an important part of the method. Two
angles, $\phi$ and $\psi$, describe the planes relative to a reference
direction (True North in this case), and a distance coordinate identifies the
stack of planes by sequence number. The fragment of the borehole between two
neighboring planes is called a slice. The polar angle $\phi\in[0,\pi/2]$ is
the angle that the normal vector $\vec{n}$ to the planes makes with the
vertical axis. The azimuthal angle $\psi\in[0,2\pi]$ is the angle that the
projection of the normal onto the base of the cylinder makes with True North.
For the degenerate case of zero $\phi$, the value of $\psi$ is set to zero.
Three examples illustrate the idea in Fig. 3 with Fig. 3(a) being a reference
configuration. In Fig. 3(b) an increase in $\phi$ coresponds to a higher angle
of dip, and in (Fig. 3(c)) the angle $\psi$ has decreased such that the planes
have the same dip but the high point of the planes is rotated to a different
position. All slices that have substantially the same angles $\phi$ and $\psi$
are designated as a slice family $S_{\phi\psi}$.

In the Cartesian coordinate system (Fig. 4), the origin $O$ is located at the
centre of the upper base of a cylinder, the $y$-axis is directed towards True
North, and the $x$ and $y$ coordinates of the point $f_{ij}$ are defined as

$x_{ij}=R\sin\alpha_{j},y_{ij}=R\cos\alpha_{j},$ (1)

where $\alpha_{j}$ is the angle between the $y$-axis and the point $f_{ij}$
for the $j$ column and is defined as

$\alpha_{j}=\frac{2\pi j}{N}.$ (2)

The $z$-axis is directed downward, and the $z$-coordinate of the point
$f_{ij}$ is given by

$z_{ij}=h(i-1).$ (3)

Assuming that all slices in the family have the same constant thickness equal
to the spacing $h$ between adjacent members of the slice family, it is
possible to enumerate the slices in the slice family by distance from the
origin $O$. The Cartesian form of the unit normal vector →n is

$x_{n}=\sin\phi\sin\psi,y_{n}=\sin\phi\cos\psi,z_{n}=\cos\phi.$ (4)

The directed length of the projection of the vector $\vec{r}_{ij}$ on the unit
normal vector is the dot product of the unit normal vector $\vec{n}$ and the
vector $\vec{r}_{ij}$

$p_{ij}=\vec{n}\cdot\vec{r}_{ij}=R\sin\phi\cos(\psi-\alpha_{j})+h(i-1)\cos\phi.$
(5)

The enumeration of the slice is then this length divided by the depth of one
slice $h$

$\nu_{ij}=\left[\frac{p_{ij}}{h}\right]=\left[\frac{R}{h}\sin\phi\cos(\psi-\alpha_{j})+(i-1)\cos\phi\right],$
(6)

where $[\cdot]$ is rounding to the nearest integer. This describes the slice
family and the number of the slice within the family.

All points in slice $p$ of an ideal slice family have the same value
$\overline{f(p)}$. The best set of planes is computed by minimising the
variance of residual of the actual data after subtraction of the approximation

$V(\phi,\psi)=var(f,S_{\phi\psi})=\frac{1}{\sum_{ij}m_{ij}}\sum_{m_{ij}=1}{\left(f_{ij}-\overline{f(\nu_{ij})}\right)^{2}}.$
(7)

It can be easily shown that the optimal slice constant $\overline{f(p)}$ is
the mean of all points $f_{ij}$ of this slice

$\overline{f(p)}=\frac{\sum_{\nu_{ij}=p,\,m_{ij}=1}{f_{ij}}}{\sum_{v_{ij}=p}{m_{ij}}}.$
(8)

Three images of real borehole walls were selected for the case study (Fig.
5(a), 5(b), and 5(c)). The initial images, the dominant dip sets, the dominant
dip sets plus 2nd, 3rd and 4th level dips, and sums of dip sets plus residual
texture (the ‘total’ column) are presented in Fig. 5 for all cases. The model
slice family is found by minimising the expression (7) over all polar and
azimuthal angles. Surfaces of this function, which depends on a particular
choice of angles, are presented in Fig. 6 for borehole walls from Fig. 5. The
second (local) minimum on the bottom surface corresponds approximately to the
second dip of the case Fig. 5(c). There appear no problems discretising the
slices when they are parallel to the borehole because of convenient
parameterisation. The Nelder-Mead algorithm (Algorithm 1) (Nelder and Mead,
1965) is used with several random starts to seek the minimum of the residual
(7). The random starts were repeated until the best solution appeared twice
within specified accuracy. Practically, the global minimum of the residual
variance has a large domain of attraction and minimisation is not
significantly difficult or computationally expensive if the depths of the
window exceeds its circumference (the width). If the aspect ratio approaches 1
(a rare situation) then the optimisation landscape becomes more noisy and a
special smoothing is needed for optimisation. If the optimisation landscape is
smooth (Fig. 6) and the domain of attraction for the global minimum is large
then the solution of the optimisation problem is expected to be robust with
respect to small perturbations.

Algorithm 1 Minimization of function $V(\phi,\psi)$ (7)

1:$\alpha\leftarrow 1,\beta\leftarrow 2,\gamma\leftarrow 0.5,\sigma\leftarrow
0.5$ $\triangleright$ Initiate parameters of the Nelder-Mead algorithm

2:$\psi_{1}\leftarrow 0,\phi_{1}\leftarrow 0,V^{*}\leftarrow
2V(0,0),\varepsilon\leftarrow 0.0001$ $\triangleright$ Initiate search
parameters

3:repeat$\triangleright$ Search of the second case of the best solution within
specified accuracy

4: if $V^{*}>V(\psi_{1},\phi_{1})$ then

5: $\psi^{*}\leftarrow\psi_{1},\phi^{*}\leftarrow\phi_{1},V^{*}\leftarrow
V(\psi_{1},\phi_{1})$ $\triangleright$ Save the best solution

6: $\psi_{1}\leftarrow 2\pi r_{1},\phi_{1}\leftarrow\pi r_{2}/2$
$\triangleright$ Generate centre of simplex; $r_{1},r_{2}$ are random numbers

7:
$\psi_{2}\leftarrow\psi_{1}+\pi,\phi_{2}\leftarrow\phi_{1},\psi_{3}\leftarrow\psi_{1},\phi_{3}\leftarrow\phi_{1}+\pi/4$
$\triangleright$ Complete simplex

8: repeat$\triangleright$ the Nelder-Mead search for specified initial point

9: Sort vertices to hold $V(\psi_{1},\phi_{1})\leq V(\psi_{2},\phi_{2})\leq
V(\psi_{3},\phi_{3})$ $\triangleright$ Ordering

10:
$\psi_{c}\leftarrow(\psi_{1}+\psi_{2})/2,\phi_{c}\leftarrow(\phi_{1}+\phi_{2})/2$
$\triangleright$ Centroid calculation

11:
$\psi_{r}\leftarrow\psi_{c}+\alpha(\psi_{c}-\psi_{3}),\phi_{r}\leftarrow\phi_{c}+\alpha(\phi_{c}-\phi_{3})$
$\triangleright$ Calculate reflection point

12: if $V(\psi_{1},\phi_{1})\leq V(\psi_{r},\phi_{r})<V(\psi_{2},\phi_{2})$
then

13: $\psi_{3}\leftarrow\psi_{r},\phi_{3}\leftarrow\phi_{r}$ $\triangleright$
Reflection

14: else if $(V(\psi_{r},\phi_{r})\leq V(\psi_{1},\phi_{1}))$ then

15:
$\psi_{e}\leftarrow\psi_{c}+\beta(\psi_{r}-\psi_{c}),\phi_{e}\leftarrow\phi_{c}+\beta(\phi_{r}-\phi_{c})$
$\triangleright$ Calculate expansion point

16: if $(V(\psi_{e},\phi_{e})\leq V(\psi_{r},\phi_{r}))$ then

17: $\psi_{3}\leftarrow\psi_{e},\phi_{3}\leftarrow\phi_{e}$ $\triangleright$
Expanding

18: else

19: $\psi_{3}\leftarrow\psi_{r},\phi_{3}\leftarrow\phi_{r}$ $\triangleright$
Reflection

20: else$\triangleright$ Now $V(\psi_{r},\phi_{r})\geq V(\psi_{2},\phi_{2})$

21:
$\psi_{co}\leftarrow\psi_{c}+\gamma(\psi_{3}-\psi_{c}),\phi_{co}\leftarrow\phi_{c}+\gamma(\phi_{3}-\phi_{c})$
$\triangleright$ Calculate contraction point

22: if $V(\psi_{co},\phi_{co})\leq V(\psi_{3},\phi_{3})$ then

23: $\psi_{3}\leftarrow\psi_{co},\phi_{3}\leftarrow\phi_{co}$ $\triangleright$
Contracting

24: else

25:
$\psi_{j}\leftarrow\psi_{1}+\sigma(\psi_{j}-\psi_{1}),\phi_{j}\leftarrow\phi_{1}+\sigma(\phi_{j}-\phi_{1}),j=2,3$
$\triangleright$ Shrinking

26: until $\min_{i,j}|\psi_{i}-\psi_{j}|\leq\varepsilon$ and
$\min_{i,j}|\phi_{i}-\phi_{j}|\leq\varepsilon$

27:until $|\psi^{*}-\psi_{1}|\leq\varepsilon$ and
$|\phi^{*}+\phi_{1}|\leq\varepsilon$

After finding the best slice family the residual field is computed by
subtracting the values on the slice family found from the known field $f$

$res_{ij}^{1}=f_{ij}-s_{ij}^{1},\forall m_{ij}=1,$

where $s_{ij}^{1}$ is the value of approximation of the field $f$ at the point
$ij$ by the best slice family. This value is equal to
$\overline{f(\nu_{ij})}$. The secondary slice structure $s^{2}$ is found by
repeating the same procedure for the residual random field $r^{1}$ instead of
$f$. Thus, recursively a sequence of layers $s^{1},s^{2},\ldots,s^{L}$ can be
removed until the statistics of the residual are more or less random (or at a
user-selected level), leaving behind a field that is considered to be residual
texture denoted by $t_{ij}$.

One value of the field is assigned per slice. This is the mean of all values
$f_{ij}$ of this slice $\overline{f(p)}$. This value is considered as a
constant when extended in the modeled volume. For each $k=1,\ldots,L$ the
slice $s^{k}$ is inpainted by this value inside the cylinder. These functions
$s^{k}(x,y,z)$ are called below the boundary elements. The boundary element
$s^{k+1}(x,y,z)$ provides the best approximation of the residual

$res_{ij}^{k}=f_{ij}-\sum_{q=1}^{k}s_{ij}^{q}$ (9)

on the cylinder boundary by the function, which is constant in every slice
from a system of parallel slices. The approximation of the image is a sum of
these boundary elements (9). This sum is not constant in any slice for $k>1$
for a general situation when the original field $f_{ij}$ is not constant in a
slice.

The first level $s^{1}$ in Fig. 5 (the second column) dominates in the
decomposition. Fig. 7 concerns case Fig. 5(b). It represents the information
brought from each level separately for five levels. In this case, the variance
of $s^{k}$ for $k>1$ is much less than the variance of $s^{1}$. The residuals
after subtraction of the approximations of various levels for the same real-
life case are presented in Fig. 8. The variance of the residual decreases at
each step. The final residual is used for generation of the residual texture.
It is clear from Fig. 8 that in this example the residual does not change
significantly after the third step.

To compute a value inside the volume, it is sufficient to identify the
appropriate slices to which the point belongs at each level of this
decomposition and then to compute the value of the residual texture inside the
volume. For this purpose, modelling of residual texture as a 3Dl moving
average field (Francos and Friedlander, 1998; Ojeda et al., 2010) is employed
with calibrating the results to the observed statistics on the boundary of the
borehole. Autocorrelations of this field are evaluated by its boundary values
and are used for continuation of the field inside the volume.

The problems of selection the best number of boundary elements is very similar
to the problem of selection of the number of principal components in the
signal approximation and there are several useful heuristics for its solution
Cangelosi and Goriely (2007). Usually, the first 3-5 boundary elements
approximate the field quite well and the rest is left for the random field of
texture. In Fig. 9 the Fraction of Variance Unexplained (FUV) plots are
presented as the functions of number of layers for the images used in the case
study (Fig. 5). The variance of residual monotonically decreases but this
decrease could be slow (Fig. 8 where the difference between level 4 and level
5 residuals is not obvious).

This version of our algorithm seeks planar dominant features; but in
principle, other geological features might be identified. Subtracting these
from the original data leaves the subsidiary (i.e., less geologically
significant) features which are isolated by an iteration of the method. Each
set of subsidiary features is subtracted from the remaining data, and the
process is repeated to isolate all of the significant subsidiary sets. A
quality of fit test is performed on each slice family as it is identified.
Examples of the sum of two-, and three- and four-slice families are presented
in the corresponding columns of Fig. 5.

To create the synthesized volume representing the rock removed during creation
of the borehole, the main function, the $2,\ldots,L$ subsidiary functions, and
the residual texture are summed in accordance with the expression

$\hat{f}_{ij}=\sum_{l=1}^{L}{s_{ij}^{l}+t_{ij}}.$ (10)

This equation is termed the reconstruction algorithm and is exact on the
boundary. Moreover, for an arbitrary internal point of the cylinder
$p=(x,y,z)$, the angle $\alpha$ is calculated as

$\alpha=\arctan(y/x)$ (11)

(that being an analogue in $\alpha_{j}$ in Fig. 4(a)) and then the number of
the slice in each slice family $k$ is calculated as

$\nu^{k}=\left[\frac{R\sin\phi^{k}\cos(\psi^{k}-\alpha)+z\cos\phi^{k}}{h}\right],$
(12)

where $\phi^{k}$ and $\psi^{k}$ are the angles that specify slice family at
level $k$. Then Eq. (10) is used to calculate the value of the field at this
point.

The model matches the empirical boundary values exactly and preserves the
statistical properties of the residual texture inside the volume. This
internal residual texture is referred to as $t$. If necessary, this model can
be improved by a multidimensional autoregressive moving average approach or
other texture analysis algorithm. There are many geological events which could
be classified as texture. Most of them have specific geometries and may
require specific tools for modelling. All these algorithms may be applied
after extraction of the planar layer structure.

The simple and universal moving average (MA) field is used for the modeling of
residuals. Let the residual field $res_{ij}=res_{ij}^{k}$ (9) on the cylindric
surface be given. A MA field $T(x,y,z)$ is constructed by averaging of 3D
white noise in a sliding window. The sliding window is selected in the form of
a rectangular parallelepiped in such a way that the horizontal and vertical
correlation radii of $T(x,y,z)$ coincide with the correlation radii of
residual $res_{ij}$ evaluated on the surface. After finding the dominant dip
for the positive squared residual field $res_{ij}^{2}$, the field $res_{ij}$
is considered in the corresponding slices $\nu$. The mean $\mu(\nu)$ and
variance $\sigma^{2}(\nu)$ of $res_{ij}$ are evaluated in each slice $\nu$.
The values of $T(x,y,z)$ are shifted and scaled in each slice $\nu$ to have
the same mean $\mu(\nu)$ and variance $\sigma^{2}(\nu)$. The result is used
for the inpainting of the texture.

This approach works satisfactorily if the residuals $res_{ij}^{k}$ are small
enough with respect to the layer structure (Fig. 8). If it is necessary, the
model could be refined by more advanced MA and autoregressive models (Francos
and Friedlander, 1998; Ojeda et al., 2010) or by other universal approaches,
which include wavelet analysis and its various generalizations (Portilla and
Simoncelli, 2000), adapted kernel regression methods (Zhang et al., 2016) and
methods of statistical physics (Wang et al., 2017) among others. If the
residual is not small then the (multi)layer structure
$\sum_{l=1}^{L}s_{ij}^{l}$ is not dominant in Eq. (10) and special tools for
the recognition and modelling of the geological structures may be needed.

The orientation of the structural or sedimentary structures may be not
constant over the window; in addition the well path can also vary, thus one
cannot expect that the sinusoids for a sufficiently long borehole are strictly
parallel. If the change of direction in the sliding window is significant then
a special segmentation procedure is needed. Fix a sliding window of the total
depth $d$, with the internal coordinate $z\in[0,d]$. Find the optimal dominant
dip in a shorter sliding window $[0,z]$ and calculate the variance of residual
(7) for this optimal approximation as a function of $z$: $v(z)$. If this
function does not vary from a constant beyond a preselected tolerance level
$\varepsilon$ then no segmentation is needed. If $v(z)$ deviates from the
constant average value by more than $\varepsilon$ and a significant slope
appears then segmentation is necessary and the point of segmentation is the
break point of the optimal piece-wise linear approximation of $v(y)$. For the
detailed description of segmentation algorithms refer to the classical review
by Keogh et al. (2004). Examples of the function $v(y)$ for the borehole
images are presented in Fig. 10 for two non-uniform images with segment
infrastructure: artificial combination of two segments (Fig. 10(a)) and a real
example with variable structure (Fig. 10(b)) taken from Fig. 5(b) . The
segmentation algorithm utilises same calculation of residuals and search of
optimal dips as the basic approximation algorithm does.

Selection of parallel planar slices for the boundary elements formalizes the
idea about planar geological structures. The choice of slices can be modified.
For example, instead of parallel planes the families of planes can be
considered, which include a given straight line outside the cylinder. In this
case, instead of the 2D optimization problem for function (7) on a hemisphere
of normal vectors $\vec{n}$ one has to solve the four-dimensional optimization
problem in the space of straight lines. Algebraic surfaces provide many other
possible choices of the families of slices but introduction of each
sophisticated construction needs reasonable geological justification.

##  3 Results

Figure 5 shows three examples of surface reconstruction from initial images
with different circumferential coverage states. The starting point for the
upper set of images is image 5a with 65% coverage typical of a well drilled
with a 203-mm (8-inch) bit. Successive slice families are shown in images 5b
to 5e, and image 5f is the reconstruction formed from the sum of slice
families plus residual texture. The reconstructed image has all of the main
geological features of the original and has additionally filled the gaps
present in the initial image. Filling the gaps in this way has introduced some
noise. The middle set of images is from another part of the same well, but
this time the initial image 5a has been inpainted with the morphological
components analysis method of Assous et al. (2014), which reconstructs data
missing in the gaps without introducing noticeable noise. The lower set of
images is from a well drilled with a smaller bit and the circumferential
coverage is 100%. This example includes many fine fractures in addition to the
crossing of planar bed boundaries.

Figure 11 shows horizontal and vertical slices through the solid images
reconstructed from the images in Fig. 5. For each case, inpainted images are
without texture and images with texture are presented.

Figure 12 shows an initial image with surface inpainting displayed in the
conventional cylinder view and also a solid-volume view created by taking two
orthogonal slices. Note that the planar geological features appear as curves
on the conventional view, and as planes in the solid volume view. Relative to
the conventional view, therefore, the angular relationships between geological
features in the solid volume view are more straightforward for nonexpert
interpreters to understand.

Finally, the method has been applied to circumferential scans of fullbore core
(Fig. 13). In this case, the external surface of solid volume has been
modulated with a function of the grey level as a proxy for grain size (an
indicator of depositional environment in sedimentary rocks). In the case of
microresistivity images, a modulation based on resistivity value can be used.

##  4 Conclusions

Image volumes created from circumferential microresistivity borehole images
and optical core scans have been used to create pseudo-rock outcrop
visualizations to help non-expert interpreters understand the earth-frame
angular relationships between geological features. Modulating the external
dimensions of the volume with a proxy for grain size adds information
indicating depositional environment. The method focusses on planar geological
features and textures, but it can be developed to include more complex
structures. The method has been demonstrated on images with 65%
circumferential coverage typical of images acquired with pad-based sensors.

###### Acknowledgements.

The authors acknowledge the help of John Winship, Technical Advisor at
Weatherford Laboratories, UK, for providing core scans.

## References

  * Assous et al. (2014) Assous S, Elkington P, Whetton J (2014) Microresistivity borehole image inpainting. Geophysics 79:D31–D39. doi 10.1190/geo2013-0188.1
  * Cangelosi and Goriely (2007) Cangelosi R, Goriely A (2007) Component retention in principal component analysis with application to cDNA microarray data. Biology Direct 2:2. doi 10.1186/1745-6150-2-2
  * Chen and Wang (2010) Chen J, Wang B (2010) High quality solid texture synthesis using position and index histogram matching. Vis Comput 26:253–262. doi 10.1007/s00371-009-0408-3
  * Criminisi et al. (2004) Criminisi A, Pérez P, Toyama K (2004) Region filling and object removal by exemplar-based image inpainting. IEEE Trans Image Process 13:1200–1212. doi 10.1109/TIP.2004.833105
  * Elad et al. (2005) Elad M, Starck JL, Querre P, Donoho DL (2005) Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA). Appl Comput Harmon Anal 19:340–358. doi 10.1016/j.acha.2005.03.005
  * Francos and Friedlander (1998) Francos JM, Friedlander B (1998) Parameter estimation of two-dimensional moving average random fields. IEEE Trans Signal Process 46:2157–2165. doi 10.1109/78.705427
  * Fu et al. (2016) Fu HC, Zou CC, Li N, Xiao CW, Zhang CS, Wu XN, Liu RL (2016) A quantitative approach to characterize porosity structure from borehole electrical images and its application in a carbonate reservoir in the tazhong area, tarim basin. SPE Reserv Eval Eng 19. doi 10.2118/179719-PA
  * Hurley and Zhang (2011) Hurley NF, Zhang T (2011) Method to generate full-bore images using borehole images and multipoint statistics. SPE Reserv Eval Eng 14:204–214. doi 10.2118/120671-PA
  * Joubert et al. (2016) Joubert JB, Millot P, Montaggioni P, Dymmock S, Andonof L, Kadri N, Torres D (2016) Understanding wireline borehole image workflows from the wellsite to the end user. First Break 34:65–78
  * Keogh et al. (2004) Keogh E, Chu S, Hart D, Pazzani M (2004) Segmenting time series: A survey and novel approach. In: Last M, Kandel A, Bunke H (eds) Data Mining in Time Series Databases. World Scientific, Singapore, pp 1–22
  * Nelder and Mead (1965) Nelder JA, Mead R (1965) A simplex method for function minimization. Comput J 7:308–313. doi 10.1093/comjnl/7.4.308
  * Ojeda et al. (2010) Ojeda S, Vallejos R, Bustos O (2010) A new image segmentation algorithm with applications to image inpainting. Comput Statist Data Anal 54:2082–2093. doi 10.1016/j.csda.2010.03.021
  * Portilla and Simoncelli (2000) Portilla J, Simoncelli EP (2000) A parametric texture model based on joint statistics of complex wavelet coefficients. Int J Comput Vision 40:49–70. doi 10.1023/A:1026553619983
  * Rider and Kennedy (2011) Rider M, Kennedy M (2011) The geological interpretation of well logs, 3rd edn. Rider-French Consulting Limited, Sutherland
  * Starck et al. (2005) Starck JL, Elad M, Donoho DL (2005) Image decomposition via the combination of sparse representations and a variational approach. IEEE Trans Image Process 14:1570–1582. doi 10.1109/TIP.2005.852206
  * Takeda et al. (2007) Takeda H, Farsiu S, Milanfar P (2007) Kernel regression for image processing and reconstruction. IEEE Trans Image Process 16:349–366. doi 10.1109/TIP.2006.888330
  * Urs et al. (2014) Urs RD, Da Costa JP, Germain C (2014) Maximum-likelihood based synthesis of volumetric textures from a 2D sample. IEEE Trans Image Process 23:1820–1830. doi 10.1109/TIP.2014.2307477
  * Wang et al. (2017) Wang H, Wellmann J, Li Z, Wang X, Liang R (2017) A segmentation approach for stochastic geological modeling using hidden Markov random fields. Math Geosci 49:145–177. doi 10.1007/s11004-016-9663-9
  * Woiselle et al. (2011) Woiselle A, Starck JL, Fadili J (2011) 3-d data denoising and inpainting with the low-redundancy fast curvelet transform. J Math Imaging Vision 39:121–139. doi 10.1007/s10851-010-0231-5
  * Xu and Torres-Verdín (2013) Xu C, Torres-Verdín C (2013) Pore system characterization and petrophysical rock classification using a bimodal Gaussian density function. Math Geosci 45:753–771. doi 10.1007/s11004-013-9473-2
  * Zhang and Du (2014) Zhang T, Du Y (2014) A multiple-point geostatistical reconstruction method of porous media using soft data and hard data. Journal of Computational Information Systems 10:7213–7224. doi 10.12733/jcis11741
  * Zhang et al. (2016) Zhang T, Gelman A, Laronga R (2016) Structure and texture-based fullbore image reconstruction. Math Geosci 48:1–21. doi 10.1007/s11004-016-9649-7

Figure 1:  Microresitivity borehole imaging tool, pad section (a), unwrapped
image (b), and circumferential coverage as a function of borehole diameter (c)
Figure 2: The standard 2D borehole image (a), 2D image after subtraction of
the dominant dip sinusoid from the ordinate (b), and cylindrical borehole
surface in 3D coordinate systems (c) Figure 3: Reference slope and orientation
of planes (a), steeper planes with increased $\phi$ (b), and rotated planes
with decreased $\psi$ (c) Figure 4: Cartesian co-ordinate system (a), True
North reference (b), and planar slices (c) Figure 5: Three images for the case
study: Image with 65% circumferential coverage (a), 65% circumferential
coverage, two-dimensionally inpainted (b), and 95% coverage (c) Figure 6: Two-
dimensional plots of the variance of residual (7) for the dominant dip
calculated for the borehole wall images Fig. 5(a), 5(b), and 5(c) (in the same
order) Figure 7: Original image and different levels $s^{k}$ ($k=1,\ldots,5$)
of the decomposition for case Fig. 5(b). The values of $s^{k}$ are rescaled to
the standard maximum and minimum Figure 8: Original image and residuals
$res^{k}$ ($k=1,\ldots,5$) after extraction of the approximations. The values
of $res^{k}$ are rescaled to the standard maximum and minimum Figure 9: FUV
plots for the images presented in Fig. 5 as functions of number of slices used
Figure 10: Dependences of the variance of residual (7) (ordinate) on the depth
of a window $z$ (abscissa) for the dominant dip: artificial combination (a)
and a real case Fig. 5(b) (b)

Figure 11: Vertical cross sections and horizontal cross sections at three
different depths reconstructed from the images in Fig. 5(a) – 5(c) in the same
order, with and without texture Figure 12: Initial 2D inpainted
circumferential image (a) and 3D volume with two orthogonal slices (b) Figure
13: Fullbore core surface scan (a) and pseudo-outcrop visualization (b)