Hanyi Wang

Shiting Yi

Mukul M Sharma

A computationally efficient approach to modeling contact problems and fracture closure using superposition method *

Keywords: Fracture Closure, Contact Problem, Asperities, Contact Stress, Fracture Width, Superposition

The shape of fractures in an elastic medium under different stress distributions have been well studied in the literature, however, the fracture closure process during unloading has not been investigated thoroughly. The fracture surface is normally assumed to be perfectly smooth so that only the width and stress intensity at the fracture tip are critical in the analysis. In reality, the creation of a fracture in a rock seldom produces smooth surfaces and the resulting asperities on the fracture surfaces can impact fracture closure. Correctly modeling this fracture closure behavior has numerous applications in structural engineering and earth sciences. In this article, we present an approach to model fracture closure behavior in a 2D and 3D elastic media. The fracture surface displacements under arbitrary normal load are derived using superposition method. The contact stress, deformation, and fracture volume evolution can be estimated in a computationally efficient manner for various fracture surface properties. When compared with the traditional integral transform methods used to model fracture closure, the superposition method presented in this study produces comparable results with significant less computation time. In addition, with the aid of parallel computation, large scale fracture closure and contact problems can be successfully simulated using our proposed dynamic fracture closure model (DFCM) with very modest computation times.

Introduction

Fracture widths in a solid elastic medium are governed by the stress field around the crack tip and by parameters that describe the resistance of the material to the imposed stress. Thus, the analysis of stresses near the crack tip constitutes an essential part of fracture mechanics. For brittle materials exhibiting linear elastic behavior, methods of elasticity are used to obtain stresses and displacements in cracked bodies. [START_REF] England | Some two-dimensional punch and crack problems in classical elasticity[END_REF] studied 2D crack problems using complex variable techniques and the crack displacement was expressed as a double integral over the fracture length. [START_REF] Sneddon | The Distribution of Stress in the Neighborhood of a Crack in an Elastic Solid[END_REF] investigated the distribution of stresses produced in the interior of an elastic solid by the opening of an internal crack under the action of pressure applied to its surface. The solutions are derived from Westergaard's stress function (1939) and the critical crack propagation pressure under uniform load is obtained based on the rate of energy dissipation that proposed by [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF]. [START_REF] Shah | Stress intensity factor for an elliptical crack under arbitrary normal loading[END_REF] presented solutions of stresses and the stress intensity factor for an elliptical crack embedded in an elastic solid and subjected to complex pressure distributions expressible in terms of a 2D polynomial function. The foundation of this work is based on the potential function proposed by [START_REF] Segedin | Some Three-Dimensional Mixed Boundary Value Problems in Elasticity[END_REF]. Bell (1979) derived the solution for stresses and displacements due to arbitrarily distributed normal and tangential loads acting on a circular crack in an infinite body, in the form of series of Bessel function integrals. Understanding the effect of stress and load conditions on crack geometry is fundamental to solving crack propagation problems in solid materials. These stresses can be a result of pressured liquid injection, such as magma driven dykes [START_REF] Rubin | Propagation of magma filled cracks[END_REF], sub-glacial drainage of water [START_REF] Tsai | A model for turbulent hydraulic fracture and application to crack propagation at glacier beds[END_REF], underground CO 2 or waste repositories [START_REF] Abou-Sayed | Evaluation of Oily Waste Injection Below the Permafrost in Prudhoe Bay Field[END_REF] or intentionally induced hydraulic fractures for stimulating hydrocarbon reservoirs [START_REF] Economides | Reservoir Stimulation[END_REF] or for geothermal energy exploitation [START_REF] Nemat-Nasser | Mechanics of elastic and inelastic solids Vol(5): Hydraulic fracturing and geothermal energy[END_REF]. Although hydraulic fracturing can be dated back to the 1930s [START_REF] Grebe | Increasing crude production 20,000,000 bbl. from established fields[END_REF], it has been in the last twenty years that this practice has become commonplace. Combining the analytic solutions of stresses and displacements with a material balance on the injected fluid, leads to the development of extensively used pressurized fracture propagation models, such as the KGD model [START_REF] Geertsma | A rapid method of predicting width and extent of hydraulic induced fractures[END_REF][START_REF] Khristianovich | Formation of vertical fractures by means of highly viscous fluids[END_REF], the PKN model [START_REF] Nordgren | Propagation of vertical hydraulic fractures[END_REF][START_REF] Perkins | Widths of hydraulic fractures[END_REF], and the radial model [START_REF] Abe | Growth rate of a penny-shaped crack in hydraulic fracturing of rocks[END_REF]. Besides using analytic methods (e.g., complex potential function method and the integral transform method) to relate load, stress and displacement, numerical methods such as discrete element methods (McClure et al. 2016;[START_REF] Sesetty | Simulation of Simultaneous and Zipper Fractures in Shale Formations[END_REF][START_REF] Zhao | Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method[END_REF], cohesive zone methods [START_REF] Bryant | Arbitrary Fracture Propagation in Heterogeneous Poroelastic Formations Using a Finite Volume-Based Cohesive Zone Model[END_REF][START_REF] Wang | SPE-168600-PA. Poroelastic and Poroplastic Modeling of Hydraulic Fracturing in Brittle and Ductile Formations[END_REF], and extended finite element methods [START_REF] Dahi-Taleghani | Numerical Modeling of Multi-stranded Hydraulic Fracture Propagation : Accounting for the Interaction between Induced and Natural Fractures[END_REF]Gordeliy and Peirce, 2013;[START_REF] Wang | Numerical Modeling of Non-Planar Hydraulic Fracture Propagation in Brittle and Ductile Rocks using XFEM with Cohesive Zone Method[END_REF][START_REF] Wang | Numerical Investigation of Fracture Spacing and Sequencing Effects on Multiple Hydraulic Fracture Interference and Coalescence in Brittle and Ductile Reservoir Rocks[END_REF]) can be used to solve stress and displacement fields that are coupled with material balance and lubrication equations to model fluid driven fracture propagation. Even though these numerical methods enable us to model fracture propagation with complex interactions under fully coupled physics, just as analytic models, they assume smooth crack surfaces. This assumption is legitimate for modeling fracture initiation and propagation when internal load exerts sufficient pressure on the fracture surface to prop it open and maintain a finite fracture aperture. However, during the unloading process, the fracture gradually closes and the stress acting normal to the fracture plane results in fracture apertures approaching the scale of the surface roughness, and fracture surfaces can no longer be treated as perfectly smooth.

With the assumption of a pressurized fracture with perfectly smooth surfaces imbedded in an elastic medium, the opposite surfaces will come into contact all at once when the internal pressure acting on fracture surfaces drops to zero, and the fracture is then completely mechanically closed. In reality, the creation of cracks seldom produces perfectly smooth surfaces. It is well established that the fracture properties of brittle heterogeneous materials such as rock, concrete, ceramics, wood, and various composites strongly depend on the microstructure and on the damage process of the material. Both these parameters are also at the source of the roughness of crack surfaces [START_REF] Morel | Scaling of crack surfaces and implications for fracture mechanics[END_REF]. This is especially the case for large scale subsurface fractures, such as those generated by hydraulic fracturing. [START_REF] Van Dam | Analysis of Hydraulic Fracture Closure in Laboratory Experiments[END_REF] presented scaled laboratory experiments on hydraulic fracture closure behavior. Their work shows that the roughness of the fracture surfaces appears as a characteristic pattern of radial grooves and it is influenced by the externally applied stress. They also observed up to a15% residual aperture long after shut-in. [START_REF] Fredd | Experimental Study of Hydraulic Fracture Conductivity Demonstrates the Benefits of Using Proppants[END_REF] used sandstone cores from the East Texas Cotton Valley formation to show sheared fracture surface asperities that had an average height of about 0.09 inches. [START_REF] Warpinski | Examination of a Cored Hydraulic Fracture in a Deep Gas Well (includes associated papers 26302 and 26946)[END_REF] reported fracture surface asperities of about 0.04 and 0.16 inches for nearly homogeneous sandstones and sandstones with coal and clay rich bedding planes, respectively. Field measurements [START_REF] Warpinski | Evaluation of a downhole tiltmeter array for monitoring hydraulic fractures[END_REF]) using a down-hole tiltmeter array indicated that the fracture closure process is a smooth, continuous one which often leaves 20%-30% residual fracture width.

Fractures play a decisive role in determining production and its decline trend. In natural fractures, the most important properties that affect the fluid flow are the fracture width, normal stress, contact area and the roughness of the fracture surfaces. These properties are all interdependent and directly affect each other. In propped fractures, the stress distribution after closure can substantially influence the crushing or embedment of proppants and the resulting conductivity [START_REF] Warpinski | Stress Amplification and Arch Dimensions in Proppant Beds Deposited by Waterfracs[END_REF]. In addition, the fracture closure behavior significantly impact our estimation and correct interpretation of in-situ stress through diagnostic fracture injection tests (McClure et al. 2016;Wang and Sharma 2017b) or pump-in and flow back tests [START_REF] Raaen | Stress determination from hydraulic fracturing tests: the system stiffness approach[END_REF][START_REF] Savitski | Revisiting Microfrac In-situ Stress Measurement via Flow Back -A New Protocol[END_REF]. The changes of fracture compliance during closure is the key input parameter for flow back analysis [START_REF] Fu | Estimating Effective Fracture Pore Volume From Flowback Data and Evaluating Its Relationship to Design Parameters of Multistage-Fracture Completion[END_REF]).

An extensive literature exists on the relationship between fracture aperture, conductivity, internal pressure, and effective stress. These studies include analytical models [START_REF] Greenwood | Contact of nominally flat surfaces[END_REF][START_REF] Gangi | Variation of whole and fractured porous rock permeability with confining pressure[END_REF]Brown and Scholz, 1985;[START_REF] Cook | Natural joints in rock: mechanical, hydraulic, and seismic behavior and properties under normal stress[END_REF][START_REF] Adams | Contact modeling-forces[END_REF][START_REF] Myer | Fractures as collections of cracks[END_REF] and numerical approaches [START_REF] Pyrak-Nolte | Single fractures under normal stress: the relation between fracture specific stiffness and fluid flow[END_REF][START_REF] Lanaro | A random field model for surface roughness and aperture of rock fractures[END_REF][START_REF] Kamali | Fracture closure and conductivity decline modeling-Application in unpropped and acid etched fractures[END_REF] as well as experimental studies [START_REF] Barton | Strength, deformation and conductivity coupling of rock joints[END_REF]Brown et al. 1998;[START_REF] Marache | Experimental and modelled mechanical behavior of a rock facture under normal stress[END_REF][START_REF] Matsuki | Estimation of closure of a fracture under normal stress based on aperture data[END_REF]. However, all these studies focus on detailed investigations of the closure of fractures modeled as two parallel plates with rough surfaces, and ignore the impact of the fracture geometry and the evolution of the stress redistribution during fracture closure.

Even though analytic methods can give closed from solutions for fracture displacement and stress distribution under arbitrary load, it requires a priori knowledge of load distribution (i.e., the arbitrary load can be written as a function of local coordinates). However, in contact problems, the total pressure/stress acting on the opposite surfaces of fracture is a dynamic process and coupled with fracture displacement, so it is not known in advance and does not have a general form that can be represented by smooth analytic functions. To resolve this issue, Wang and Sharma (2017a) presented an integral transform method to model contact problem with analytic solutions that relate load and displacement, the solution integrals have to be divided into piecewise segments to approximate the distribution of possible discrete internal loads, however, this method is not very efficient because the number of integrals that need to be evaluated grows exponentially with increasing number of segments. Using numerical methods (such as finite element method, finite volume method, boundary element method, etc.) to model fracture closure and contact problems also face challenges, because in many cases there is a sharp contrast of scales between fracture geometry (e.g., tens or hundreds of meters for subsurface fractures) and surface asperities (normally in the range of µm or mm). For a numerical model to capture both large scale fracture deformation and small scale contact behavior during fracture closure, the mesh or element size has to be extremely small (equivalent to the scale of the fracture aperture), which makes this effort very computationally expensive. In addition, contact problems are often non-linear and require a large number of iterations to converge, this makes numerical modeling of fracture closure processes more challenging, especially when the fracture size is large or a large number of fractures need to be modeled simultaneously.

In this article, we propose a new method and general algorithms to model the dynamic behavior of fracture closure on rough fracture surfaces and asperities. Analytic solutions that relate fracture width with arbitrary normal net pressure/stress using superposition are derived. Coupling analytic solutions with a general contact law that up-scales surfaces roughness and asperities, large scale 2D and 3D fracture closure behavior can be modeled efficiently. Different methods of modeling fracture closure behavior and their computation efficiency are also compared and investigated.

Fracture Surface Displacement with Arbitrary Normal Load

When the internal pressure of a fracture declines, the fracture aperture diminishes accordingly until the fracture faces come into contact when the asperities of opposing fracture walls touch each other. This induces an additional normal load (i.e., contact stress) acting on the area of contact and changes the subsequent fracture closure behavior as the internal pressure declines further. Here we are not seeking a general solution that relates displacement to the load, because the load itself is dynamically related to fracture geometry, material properties and surface roughness during the fracture closure process. Instead, we determine the displacement based on superposition, by adding the influence of discretized fracture segment/surface areas that subject to constant computed loads. In this section, we present a method to determine fracture aperture under an arbitrary normal load in a two and three dimensional elastic medium using superposition.

1 Two Dimensional Fracture

The problem of a pressurized line crack in an infinite, elastic domain is defined by the following conditions:

The fracture in Cartesian coordinate x-y is shown in Fig. 1.

is the pressure acting on the crack surface. All displacements, and stresses, are zero at infinity. Linear elasticity and plane strain are assumed.

Fig.1 2D Crack ( ) in an infinite domain

The stress equilibrium equations in a 2D domain can be expressed as

These equations are satisfied if we assume where is termed the Airy Stress Function. Hooke's Law that relate stress and strain:

where is Young's modulus and is Poisson's ratio, and , are normal and shear strain respectively. For plane stress condition, and for plane strain condition, .

A convenient equation that represent three strains are defined in terms of derivatives of only two displacements is given by Substituting Eq.( 5) into Eq.( 6), and using Eq.( 7), the stress equilibrium can be expressed as, [START_REF] Westergaard | Bearing Pressures and Cracks[END_REF] found an Airy Stress Function of complex numbers can be the solution for the stress field in an infinite plate containing a crack. He discussed several Mode I crack problems that could be solved using:

where and Im denote the real and imaginary part, and z is a complex number. Then, the stresses can be calculated by taking derivatives of the Airy Stress Function:

Where,

For plane strain and symmetric loading on y=0, , and the displacement field at x-direction, and y-direction, can be determined:

The use of Airy Stress Function is a powerful technique for solving 2D equilibrium elasticity problems. After finding the Airy Stress Function for a specific problem, all the stress field and displacement can be determined. Considering a case of normal traction p from x=b to c is applied to a line crack |x|≤ a in an infinite elastic body (shown in Fig. 2). The complex number Z for Westergaard solution for Mode I fracture is given by (Tara et al. Finally, the width of a fracture, , along x-direction can be calculated by set y=0 in Eq.( 14):

where is plane strain Young's modulus.

For a special case, when b=-a and c=a, Eq.( 17) is reduced to the analytical solution for a constant internal pressure:

If the fracture from -a to a is divided into a total number of n line segments, and in the i th ( ) segment, there exists a uniform pressure that acts on the opposite fracture surface. If n is large enough, the pressure distribution can be approximated to any distribution of . The solution of fracture width for each individual segment with constant pressure is given in Eq.( 17). Use superposition by adding the influence of each individual segment over the entire fracture, the final fracture width profile for any given pressure distribution can be obtained:

2 Three Dimensional Fracture

For the three-dimensional case assume that a crack is created in the interior of an infinite elastic medium and that it is "penny-shaped" with a radius of R, as shown in Fig. 3. The fracture is axially symmetry about the z-axis, and the internal pressure is a function a radius r. In the absence of body force, the equations of elastic equilibrium in a cylindrical coordinate system reduce to [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF]:

At the fracture surface (

), the boundary conditions are:

Fig. 3 Penny shaped fracture in cylindrical coordinates.

Through integral transforms, the fracture aperture for a given internal pressure,

, can be obtained as a double integral [START_REF] Sneddon | The Distribution of Stress in the Neighborhood of a Crack in an Elastic Solid[END_REF]:

where is a dummy variable, is fracture width along r-direction. Eq.( 24) is not suitable for contact problems because is itself coupled with , which is not known in the first place. To determine based on superposition, we first need the solution of when is constant over a circular area. Let's assume p(r) is imposed within a radius : Let be the inner integral of Eq.( 24):

Eq.( 26) has to be evaluated piecewise for different ranges of r. When

When

The outer integral of Eq.( 24) also needs to be evaluated in a piecewise manner. When When Eq.( 29) and ( 30) give the fracture width distribution under the load given by Eq.( 25). The integrals in the solutions can't be evaluated analytically without resorting to elliptic integrals of the second kind [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. With the above derivation, the solution of fracture width as a function of fracture radius, for a constant pressure applied in the circular region of , can be obtained using superposition, as shown in Fig. 4. Using the solution of the constant pressure acting on minus the solution of the constant pressure acting on , we get the fracture width distribution for a constant pressure acting on :

For a special case, when a=0 and b=R, Eq.( 31) reduces to the well-known analytical solution of a penny-shaped crack for a constant internal pressure:

Fig. 5 shows a penny-shaped fracture discretized into a total number of n circular segments along the fracture radius. In each segment, there exists a uniform pressure that acts on the opposite fracture surface. If n is large enough, the pressure distribution can be arbitrary, . The solution for fracture width for an individual segment is given in Eq.( 31). Again, we use the superposition principle by adding the influence of each individual segment over the entire fracture to get the final solution of fracture width for any given pressure distribution:

There are two main advantages of using superposition. One is that the mathematical form of the pressure distribution, , is not required a priori, the other is the pressure distribution does not necessarily need to be continuous, which suits the nature of contact problems.

Solution Dependent Contact Load

If the fracture surface is perfectly smooth, according to Eq.( 18) and Eq.( 32), the entire fracture surface area will come into contact all at once when the internal pressure drops to zero. However, as discussed previously, depending on material properties, heterogeneity and loading conditions, the generated crack or fracture surface is always rough with asperities distributed across the surface area. So fractures can retain a finite aperture after mechanical closure due to a mismatch of asperities on the fracture walls, this is especially the case for subsurface fractures that are stimulated in hydrocarbon and geothermal reservoirs. Scaled laboratory experiments show that the created fracture surfaces are rough (shown in Fig. 6) Detailed measurements and modeling of surface roughness and mechanical properties of asperities for every fracture encountered is not practical and, many times not even possible, so it is desirable to study and investigate specimen samples, and upscale the influence of surface microscopic structure to some macroscopic empirical relationships (i.e., a contact law) that relates fracture width and the associated contact stress. [START_REF] Zangerl | Normal stiffness of fractures in granitic rock: a compilation of laboratory and in-situ experiments[END_REF] compiled laboratory and in-situ experiments of fracture closure behavior for different rock samples and showed that the effective normal stress and fracture normal closure (the displacement from zero effective normal stress condition) is highly non-linear, as shown in Fig. 7. is the contact stress acting normal on the fracture, and is the effective normal stress at which the aperture is reduced by 90%. is the fracture aperture and, is the fracture aperture when the contact stress is equal to zero. When , the contact stress is zero because the fracture walls are mechanically detached. It should be emphasized that the contact width is determined by the tallest asperities, and the strength, spatial and height distribution of asperities are reflected by the contact reference stress (e.g., if the tallest asperities on two fracture samples are the same, then they should have the same , but the one with a higher median asperity height or higher Young's modulus will have a higher value of , provided other properties are the same). The influence of all the rock properties, modulus, asperities patterns, density, distribution, etc., can be up-scaled into two contact parameters, and , using Eq.( 34) to match the (e.g., curves in Fig. 7) displacement vs loading stress curves for fractured rock samples

Coupling Superposition Solutions with Contact Law to Model Fracture Closure

In this section, we present the algorithm to determine fracture aperture and contact stress distribution with a given fracture geometry, rock properties, fracture internal pressure and parameters for the contact law. When the fracture closes on rough surfaces and asperities as the internal pressure drops, the total pressure acting on the fracture surfaces becomes the internal pressure plus contact stress. The total pressure distribution impacts surface displacement and in turn, determines the contact stress, so the contact stress and fracture width distribution need to be solved simultaneously. Fig. 8 shows a schematic of the algorithm used for solving the fracture contact problem based on a Picard iteration [START_REF] Berinde | Picard iteration converges faster than Mann iteration for a class of quasi-contractive operators[END_REF]). During the first step, if the contact stress is not known, it is assumed to be zero, and the resulting fracture width is smaller than the actual width if some portion of the fracture tip is already in contact. Based on this initial guess of the fracture width, we can calculate the corresponding contact stress and get the updated total pressure distribution. Convergence is achieved when the absolute change of total pressure between two consecutive iterations is smaller than a specified tolerance. The total net pressure/stress, acting on the fracture surface is determined by: where is fluid pressure inside the fracture, is the far field stress that is perpendicular to the fracture face. From

Eq.( 34), we can observe that the relationship between contact stress and fracture aperture is highly non-linear, and small changes of fracture aperture can lead to very large increases in contact stress. This poses great challenges for a coupled iterative model. The Picard relaxation factor, , avoids over correction in each iteration, and has an value between 0 and 1. The proposed algorithm of modeling fracture closure behavior is a general algorithm, where different contact laws and different methods that relate displacement and loading conditions can be applied.

Validation of Superposition Implementation

In this section, we examine and validate the superposition method implemented in this study by comparing the fracture width distribution with analytic solutions under uniform internal pressure. The fracture is divided into a number of segments along the fracture length or fracture radius, and a pressure is imposed on each segment to represent a uniform internal pressure condition. First, a 2D plan strain case is examined as shown in Fig. 9. The fracture length is 50 m, Young's modulus is 20 GPa, Poisson's ratio is 0.25 and the internal pressure is 5 MPa. The results from the superposition method, using Eq.( 19), with a different number of segments are compared with the analytic solution, using Eq.( 18). It can be observed that when the total number of segment is only 6 or 8, the superposition method overestimates fracture width, but when the total number of segment increases to 10, the results from the superposition method agree very well with the analytical solution.

Next, we examine a 3-D penny-shaped fracture with the same Young's modulus, Poisson's ratio and internal pressure as the 2-D case, and set the fracture radius to be 50 m. The results from superposition method, using Eq.( 33), with different numbers of segments are compared with the analytical solution, using Eq.( 32),are shown in Fig. 10. Results from the superposition method match very well with the analytical solution. It seems that the superposition method for a 3-D pennyshaped fracture is less sensitive to the total number of segments or segment size than the 2-D fracture.

Fracture Closure and Contact Stress

In this section, the dynamic fracture closure process with declining internal pressure and the impact of surface roughness properties (represented by up-scaled contact parameters and ) will be examined for both 2-D and 3-D fracture geometries.

1 2-D Fracture Closure

1.1 Uniform Far-Field Stress

A Base Case was set up to demonstrate the method. The fracture has a tip-to-tip length of 20 m (a =10 m in Fig. 2) with a far-field stress of 20 MPa. The surrounding rock has a Young's modulus of 20 GPa, Poisson's ratio of 0.25. On the fracture surface, the contact width is 2 mm and the contact reference stress is 5 MPa. The initial fluid pressure inside the fracture is 25 MPa. The length of the discretized segment for all the following simulated cases is set as 0.2 m. Using the proposed algorithm (Fig. 8), the evolution of the fracture width profile and contact stress distribution can be determined, as the fluid pressure inside the fracture gradually declines. The results are shown in Fig. 11 and Fig. 12. To demonstrate the impact of fracture roughness and surface asperities on fracture closure behavior, the case without surface asperities (fracture surface is completely planar and smooth) is also included. The result shows that at a relatively high fracturing fluid pressure, the fracture asperities have negligible impact on fracture width distribution, and the contact stress is always concentrated at the tip of the fracture, where the contact stress is much higher than in the middle of the fracture. We can also see that the fracture surfaces do not contact each other like parallel plates. In fact, the fracture closes on rough surfaces starting from the tip, and closes progressively all the way from the edges to the center of the fracture. As the fluid pressure continues to decline, more and more of the fracture surfaces come into contact and this changes the subsequent fracture closure behavior. At lower fluid pressures, contact stresses start to counter-balance the far-field stress and the fracture becomes stiffer and less compliant. If the fracture faces were perfectly parallel and smooth, the fracture width would have collapsed to zero when the fluid pressure dropped to 20 MPa (i.e., the far-field stress). However, because of the surface asperities and mismatches, the fracture still retains a residual fracture width even when the fluid pressure inside the fracture drops below the far-field stress. Fig. 13 shows the fracture width and the corresponding contact stress for different values of contact width with a contact reference stress of 5 MPa using Eq.( 34). As expected, when the fracture width is larger than or equal to the contact width, the contact stress is zero. However, when the fracture width near the fracture tip is smaller than the contact width, the contact stress and fracture width follow a non-linear relationship, resembles that of Fig. 7. As the fracture width gets smaller and smaller, more force is needed to further close the fracture (less compliant fracture). Fig. 14 shows the evolution of fracture average width for different contact width while other input parameters remain the same as the Base Case. As can be observed, when the pressure is relatively high (above 23 MPa), the fracture average width and fluid pressure follows a linear relationship, since the impact of asperities is negligible. However, as the pressure continues to decline, a larger contact width leads to an earlier departure from the linear relationship. This is because the rough fracture walls will come into contact at a higher pressure if the contact width is larger, and the larger contact width results in larger residual fracture width when the fluid pressure drops to 15 MPa. 34). It can be observed that for the same contact width, the higher the contact reference stress, the more rapid the increase of contact stress as the fracture width decreases. Physically, the contact reference stress represents how hard and compliant the fracture surface asperities are. Fig. 16 shows the evolution of fracture average width for different contact reference stress while other input parameters remain the same as the Base Case. The results reveal that the contact reference stress does not have much impact on the pressure at which the fracture asperity and pressure depart from a linear relationship, but it does impact the fracture width evolution after most of fracture surface has come into contact. The lower the contact reference stress, the smoother the change in fracture average width as pressure declines. A lower contact reference stress also leads to smaller residual fracture width. There is compelling field evidence to indicate that hydraulic fracture closure occurs on rough fracture surfaces in the field. Fig. 17 shows the normalized tilt-meter response plotted against wellbore pressure during the shut-in period of well-2B at three different stations from the GRI/DOE M-site. Each tilt-meter response is normalized by dividing the maximum tilt (displacement, also an indication of fracture width) measured at that instrument during the test. The tilt-meter demonstrates that soon after shut-in, the measured displacement declines linearly with pressure (roughly constant fracture stiffness). After the wellbore pressure declines to a certain level, the measured displacement vs pressure departs from a linear relationship.

Even though these data were measured after weeks of shut-in, we can still observe the existence of a residual fracture width supported by surface asperities and mismatches even after the fracturing fluid pressure drops below the minimum in-situ stress. This field measurement is consistent with the general trend shown in Fig. 14 and Fig. 16. If the fracture geometry and rock properties can be obtained with certainty, the tilt-meter measurement, along with diagnostic fracture injection tests (DFITs), can be used to infer the contact parameters (i.e., and ) using our proposed model. The 2D fracture closure model presented here can be applied to both PKN and KGD fracture geometry, the difference is that for a PKN type fracture geometry (fracture length fracture height), a plane strain condition is assumed in the direction of fracture length, while for KGD type fracture (fracture height fracture length), the plane strain condition is assumed in the direction of fracture height.

1.2 Non-Uniform Far-Field Stress

2D fracture models have been used for decades with reasonable success. However, in certain formations, neither of these two models (KGD and PKN) can be used successfully to model hydraulic fracture propagation. Today, pseudo-threedimensional (P3D) models are the most commonly used hydraulic fracture model. P3D models, are in essence an extension of the classical PKN, but unlike the PKN model, the height of the P3D fracture is not limited to the reservoir thickness, as shown in Fig. 18, the fracture is allowed to grow vertically into the adjacent layers, and the fracture width distribution is still based on the 2D plane strain conditions (no strain in the x-direction) of the PKN model. Using the superposition method of Eq.( 19), the fracture width distribution under an arbitrary far-field stress can be obtained by first using Eq.( 35) to get the total net pressure/stress distribution, and then the dynamic fracture closure process follows the algorithm that presented in Fig. 6. Fig. 20 shows the fracture width distribution across three layers and the corresponding contact stress at different internal fluid pressure. Again, we can notice that at relative high pressure, contact stresses are only concentrated at the fracture edges (top and bottom) and the shape of fracture width is roughly elliptic. As pressure continues to decline, more of the fracture surface comes into contact and the fracture width distribution becomes more influenced by far-field and contact stresses. As expected, the fracture has a larger width in the middle layer where the far-field stress is the smallest, and a smaller fracture width in the bottom layer, where the far-field stress is the highest. Because the contact stress is closely coupled with fracture width, the bottom layer has the highest contact stress. Fig. 21 shows the resulting evolution of fracture average width. There is clear departure from the linear relationship between fracture average width and fluid pressure when the fluid pressure drops close to 24 MPa, which is far-field stress in the bottom layer.

2 3-D Fracture Closure

Radial or penny-shaped fractures are often created if the fracture is only propagating in a single, relative thick rock layer or the stress difference between the penetrated rock layers is small. The algorithm for modeling dynamic fracture closure behavior presented in Fig. 8 is a general approach. To model radial fracture closure, the fracture width distribution under an arbitrary normal load is calculated using Eq.( 33). Assuming a radial fracture with a radius of 20 m, and all the other input parameters remaining the same as the previous Base Case, Fig. 22 and Fig. 23 show the evolution of fracture average width for different contact width and contact reference stress , respectively. Similar to the results obtained for a 2D fracture geometry, larger contact width leads to earlier departure from the linear relationship between fracture average width and fluid pressure. Also the contact reference stress strongly affects the fracture width evolution after a certain portion of the fracture surface already in contact, but has less influence on the pressure at which the contacting asperities begin to influence the overall fracture closure behavior.

Comparisons of Computation Efficiency

To demonstrate the computational advantage of using the superposition method to determine fracture surface displacement and dynamic fracture closure, the computation speed for both 2D and 3D fracture problems are examined using both the superposition method and the piecewise integral transform (Appendix). Unlike the superposition method, the integral transform method calculates fracture surface displacement through nested integration operators. Fig. 24 shows the fracture width profile calculated from these two different methods for 2D and 3D fractures using the parameters from the previous validation cases. The results show that even with a relatively small number of discretized segments (i.e., the fracture length or radius of 50 m is divided into 10 segments), both superposition and integral transform methods provide accurate fracture width profile predictions.

Fig.24 Comparison of fracture width calculation for different methods with 10 discretized segments

To test the computational efficiency of superposition and integral transform methods, dynamic fracture closure cases under uniform far-field stress are simulated. For a 2D fracture, the Base Case from Section 6.1.1 is used. For a 3D fracture, the radial case from section 6.2 is used. In all scenarios, the internal pressure drops from 25 MPa to 15 MPa with pressure intervals of 0.1 MPa. Fig. 25 and Fig. 26 show the total computation time needed to complete each simulation case for 2D and 3D fracture closure, respectively, with different numbers of discretized segments. The results indicate that the total computation time and the total number of discretized segments follows a linear relationship on a log-log plot, regardless of which method is used. In addition, we can observe that the superposition method dramatically improved the computational speed for both 2D and 3D cases (by 3 to 4 orders of magnitude). From our previous analysis, we know that the superposition method yields accurate results even with relatively coarse discretization. However, to better capture dynamic fracture closure behavior and the sharp contact stress profiles, finer discretization may be needed. So to make our dynamic fracture closure model more suitable for handling large scale fractures, parallel computation has been implemented. To test the benefits of parallelization, both 2D and 3D cases are simulated on a workstation with 14 CPUs and the results are presented in Fig. 27. For the 3D cases, we notice that parallelization reduces computation time for all cases, but the difference becomes smaller when the total number of discretized segments drops below 50. When the total number of discretized segments is above 500, parallelization increases 3D computation speed roughly 10 times. However, for 2D cases, parallelization has a completely different impact on computation efficiency. As can be observed, parallelization increases computation speed only when the total number of discretized segments is above 500, and the more the total number of discretized segments, the more the speedup that can be achieved with parallelization. When the total number of discretized segments is below 500, it actually takes more time to complete a 2D fracture closure simulation with parallelization than just using a single CPU. This is because calculating 2D fracture geometry is relative fast for a small number of discretized segments, and under such circumstances, the amount of time required to coordinate parallel tasks between different CPUs (message passing overhead) surpasses that of doing useful calculations, so parallelization becomes detrimental and increases the total computation time. One way to mitigate this issue is to optimize the number of parallelized CPUs based on the total number of discretized segments through trial and error, but in general, using a single CPU is adequate for a 2D fracture closure problem. Overall, the dynamic fracture closure model (DFCM) presented in this study provides us an approach to model fracture closure process accurately and efficiently. By combining this superposition method with parallel computation, even large scale fracture closure and contact problems can be successfully simulated within minutes.

Conclusions

Understanding the role of stress and loading conditions on fracture geometry and the impact of fracture surface asperities and roughness is fundamental to solving fracture closure and contact problems. A wide range of applications, such as underground CO 2 or waste repositories, stimulating hydrocarbon reservoirs, geothermal energy exploitation, measurement of minimum in-situ stress from fracture injection test are impacted by such fracture closure problem solutions. The shape of fractures under certain loads in an elastic medium and the corresponding distribution of stress have been well studied in the literature, however, the fracture closure process during unloading has not been investigated thoroughly. Fracture closure problems are challenging because one needs to capture both large scale fracture deformation and small scale contact deformation. In addition, the non-linear nature of contact deformation problems poses serious convergence issues in iterative solutions. These issues make the use of fully-coupled numerical models to simulate fracture closure computationally prohibitive.

In this study, the fracture surface displacements under arbitrary normal load are derived using superposition method. A dynamic fracture closure model (DFCM) for both 2D and 3D fracture geometry is proposed. This model is capable of simulating the fracture closure process due to fluid depressurization under the influence of rough fracture walls and asperities. When compared with previous work using integral transforms, the new model dramatically improves the computational speed (1000 to 10,000 fold increase in speed) with no compromise on accuracy. With the aid of parallel computation, large scale fracture closure and contact problems can now be successfully simulated using our proposed model within a very reasonable time frame.

1,E-01 1,E+00

Fig. 2

 2 Fig.2 Partial load on a 2D line crack with a constant pressure p applied on

Fig. 4

 4 Fig.4 Solution of a constant pressure applied on

Fig. 5

 5 Fig.5 Illustration of discretizing fracture radius into a number of segments with uniform pressure within each segment.

Fig. 6

 6 Fig.6 Example of fracture surface roughness pattern generated by a hydraulic fracture (van Dam et al. 2000)

Fig. 7

 7 Fig.7 Compiled laboratory and in-situ experiments showing highly non-linear normal closure behavior (Modified from Zangerl et al. 2008) To quantify the non-linear fracture closure behavior, Willis-Richards et al. (1996) proposed a relationship to estimate contact stress for a certain fracture aperture, based on the work of Barton et al. (1985):

Fig. 8

 8 Fig.8 Solution flow diagram for solving coupled contact stress and fracture width distribution.

Fig. 9

 9 Fig.9 Validation of superposition implementation for 2D fracture with different numbers of segments.

Fig. 11

 11 Fig.11 Fracture width evolution with and without surface asperities at different fluid pressure under uniform far-field stress.

Fig. 13

 13 Fig.13The relationship between contact stress and fracture width for different .

Fig. 15

 15 Fig.15 The relationship between contact stress and fracture width for different

Fig. 17

 17 Fig.17 Normalized tilt-meter data from shut-in of Well-2B from the GRI/DOE M-site. Data Courtesy of Norm Warpinski.

Fig. 18 Fig. 19

 1819 Fig.18 Illustration of pseudo-3D fracture model (Adachi et al. 2010)Despite the existence of abundant literature on P3D fracture propagation, few studies have investigated fracture closure behavior in multi-layered formations. To demonstrate the capability of our model, a case with a fracture penetrating three rock layers with different far-field stresses will be presented. Fig.19shows our example case where the fracture is contained in three layers, and both of the top and bottom layer has a higher far-field stress than the middle layer. The total fracture height is 20 m and all the other input parameters remain the same as the previous Base Case.

Fig. 20

 20 Fig.20 Fracture width and the corresponding contact stress distribution at different fluid pressure under a non-uniform far-field stress.

Fig

  Fig.22 Fracture average width evolution for different during closure with 3D fracture geometry.

  Fig.25 Computation speed comparison for 2D fracture closure process

Fig. 27

 27 Fig.27 Computation speed comparison with and without parallelization.

Appendix: Fracture Surface Displacement with Arbitrary Normal Load using Nested Integral Transform [START_REF] England | Some two-dimensional punch and crack problems in classical elasticity[END_REF] showed that the width of a 2D fracture with any arbitrary pressure/stress profile is given by: where is a dummy variable, y is the distance to the middle of the fracture, is the half fracture height (PKN model) or half fracture length (KGD model). Let be the net pressure/stress distribution along the y-direction, which is the total pressure/stress acting on the fracture surfaces minus the far field stress in the fracture opening direction. To obtain and , has to be divided into even and odd functions, such that , with being even and being odd. Then and are determined by:

where is a dummy variable. If the is symmetric to the middle of the fracture , then . For example, assuming a uniform pressure inside the fracture, , , From Eq.(A2):

Substituting Eq.(A4) into Eq.(A1), we have, Eq. ( A5) is the well-known equation for a 2D static pressurized crack. If is the fracture half-height, then Eq.(A5) describes the fracture width distribution for the PKN fracture geometry, and if is the fracture half-length, Eq.(A5) describes the fracture width distribution for the KGD fracture geometry. When an analytic form of is not available a priori (e.g., coupled with other physics) or is not a continuous function, numerical integration to calculate the fracture width has to be implemented in a discretized manner along the fracture height or length. Assume a fracture half-height or half-length is discretized into a number of n segments. In each segment, there exists a uniform net pressure/stress that acts on the fracture surfaces. Then the function of in the m th segment with net pressure is calculated as:

Similarly, the integration to obtain also needs to be integrated piecewise. The function of in the m th segment is calculated as: Substitute Eq.(A6) into Eq.(A7):

We define the integration operator I: Re-write Eq.(A8) using the definition from Eq.(A9), we have Notice that is the summation of the integration operator I that operated on different bounds and integrands, we put all the terms in the R.H.S of Eq.( A10) in to a matrix, called MatrixI(k,j), with k rows and j columns. Then Eq.( A10) becomes

Re-arrange , we have Based on Eq.(A12), an algorithm can be developed to calculated fracture width within the m th segment of a 2D fracture, and this process has to be repeated from the 1 st segment to the n th segment to obtain the entire fracture width profile.

The radial crack opening due to an excess pressure distribution in a linearly elastic material of infinite extent can be expressed as [START_REF] Sneddon | Crack Problems in the Classical Theory of Elasticity[END_REF]:

where is the length of fracture radius, and are dummy variables and is the normalized radius, which is defined as where is the local fracture radius. Assume a fracture is discretized into a number of n segments along radial direction, if belongs to the m th segment where , then

We define integration operators R1 and R2:

We can rewrite Eq.(A15) using Eq.( A16) and (A17),to get:

From Eq.(A18), we can observe that can be expressed by the summation of two 1D matrices:

where Based on Eq.(A19) to Eq.(A21), an algorithm can be developed to calculated fracture width within the m th segment of a radial fracture, and this process has to be repeated from the 1 st segment to the n th segment to obtain the entire fracture width profile.


Sumera Sattar

Haris Ahmad Khan

Khurram Khurshid
email: khurram.khurshid@ist.edu.pk

OPTIMIZED CLASS-SEPARABILITY IN HYPERSPECTRAL IMAGES

Keywords: Dimensionality reduction, hyperspectral image visualization, class separability

Image visualization techniques are mostly based on three bands as RGB color composite channels for human eye to characterize the scene. This, however, is not effective in case of hyper-spectral images (HSI) because they contain dozens of informative spectral bands. To eliminate redundancy of spectral information among these bands, dimensionality reduction (DR) is applied while at the same trying to retain maximum information. In this paper, we propose a new method of information-preserved hyperspectral satellite image visualization that is based on fusion of unsupervised band selection techniques and color matching function (CMF) stretching. The results show consistent, edge-preserved and pre-attentive feature less images with high class separability. Different visualization techniques are compared to demonstrate the effectiveness of our scheme that can prompt an important advancement in the field.

INTRODUCTION

Hyper-spectral images contain wide range of spectral information in adjacent narrow bands. However, as three bands are used in display systems, rendering hyper-spectral images directly on conventional display systems is not feasible. HSI are used in various applications [START_REF] Stein | Hyperspectral imaging for intelligence, surveillance, and reconnaissance[END_REF][START_REF] Ardouin | A demonstration of hyperspectral image exploitation for military applications[END_REF] thus having capability of a good meaningful visualization is a prime requirement. Fundamentally there can be two ways to display HS data. First, based on classification where we display each class with different color thus resulting in a good but unnatural looking image. The second is to map N components to three for display but it results in loss of information. Band selection techniques for HS image visualization are categorized into supervised and unsupervised methods. Supervised methods like minimum estimated abundance covariance (MEAC) depend on spectral signature of data [START_REF] Huqqani | Comparative Study of Supervised Classification of Urban Area Hyperspectral Satellite Imagery[END_REF]. On the other hand, unsupervised methods do not rely on any prior knowledge but they are computationally expensive. For example, linear combination of the data has been utilized to acquire three new bands in method called Principal component analysis (PCA) [START_REF] Jia | Segmented principal components transformation for efficient hyperspectral remote-sensing image display and classification[END_REF]. But in certain cases, its performance drops down due to badly correlated data. Therefore, a well-organized and suitable means of visualization of the HSI is needed [START_REF] Schowengerdt | ‖Remote Sensing-Models and Methods for Image Processing‖[END_REF].

A lot of work has been done in fusion of data. Some bands are highly informative during fusion and weak features are suppressed. Recently saliency based visualization (SBV) of HSI method has been proposed [START_REF] Ahmad Khan | Saliency Based Visualization of Hyper-Spectral Images‖ IEEE Trans[END_REF] in which spatial and spectral information about all bands was taken into account, but this approach doesn't possess class separability property. So an efficient, low complexity, computationally efficient, highly informative and more class distinctive method are required. To address the previously stated issues and to make utilization of this rich source of data, we propose a new HSI visualization technique by trying to overcome the problems using different useful band selection techniques.

The rest of this paper is arranged as follows: different band selection systems are presented in section 2. In section 3 proposed visualization approaches is presented with two schemes/variations. The experimental results are presented in section 4.

LITERATURE REVIEW

The visualization using bilateral filtering is based on calculating the weight of each pixel in band, the weight of each pixel in whole HSI is individually calculated, depending on the relative contribution of the pixel in the image [START_REF] Kotwal | Visualization of Hyperspectral Images Using Bilateral Filtering[END_REF]. One bit transforms (1BT) method is investigated based on using all bands of HSI data for selecting the appropriate bands for RGB show, determination of bands is performed completely over 1BT representations by means of typically Boolean operation [START_REF] Demir | A low-complexity approach for the color display of hyperspectral remote-sensing images using one-bit transform-based band selection[END_REF]. The segmented principle components transformation (PCT) is the fusion method, three methods for dividing the bundle of bands into subgroups is presented. First principle component of every subgroup is utilized for display [START_REF] Anastassopoulos | Fusion of hyperspectral data using segmented pct for color representation and classification[END_REF]. A Pixel-level fusion [START_REF] Tsagaris | Multispectral image fusion for improved RGB representation based on perceptual attributes[END_REF] is proposed based on perceptual attributes that display an image with perceptually high quality. Typical DR methods to display RGB image principle components analysis is investigated [START_REF] Tyo | Principal components-based display strategy for spectral imagery[END_REF], [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF]. Noise-adjusted principal component analysis (NAPCA) is performed in two steps. The first step is noise whitening and other step is to perform ordinary PCA to the noise-whitened data because of the unity of the noise whitened data in noise variance [START_REF] Roger | A faster way to compute the noise-adjusted principal components transform matrix[END_REF]. Other transformation-based color display approaches are proposed by Jacobson et al [START_REF] Jacobson Andm | Design goals and solutions for display of hyperspectral images[END_REF], [START_REF] Jacobson | Linear fusion of image sets for display[END_REF], where to produce naturallooking images fixed weighting linear spectral envelopes are used that displays the same materials in different data sets consistently but they will be incapable to show some different objects distinctively because there is no adaptation to specific image information. Hyvarinen et al. [START_REF] Hyvärinen | Independent Component Analysis[END_REF] proposed the idea to decompose a group of multivariate signals into statistically independent component sources with minimum loss of object information. Wang et al [START_REF] Wang | ‖Independent Component Analysis-Based Dimensionality Reduction with Applications in Hyperspectral Image Analysis‖[END_REF] investigated DR-ICA that can capture information using virtual dimensionality to approximate the quantity of dimensions needed. Saliency based visualization [START_REF] Ahmad Khan | Saliency Based Visualization of Hyper-Spectral Images‖ IEEE Trans[END_REF] Investigated weighted fusion method by calculating the saliency in each band of HSI. Both spatial and spectral information are taking into account. They applied the saliency map technique on PCA image by dividing data into three equal parts. For final visualization fused PCA, CMF and saliency based tri-color images. But these techniques don't ensure the different class information in the final image. Supervised band selections methods are used for visualization e.g. Minimum estimated abundance covariance (MEAC) [START_REF] Yang | An efficient method for supervised hyperspectral band selection[END_REF] depend on spectral signature and many other mentions in [START_REF] Su | Hyperspectral Image Visualization Using Band Selection[END_REF].

METHODOLOGY

In this paper, we investigate a new methodology for displaying HS image considering two major goals: class separability and consistent rendering [START_REF] Jacobson | Linear fusion of image sets for display[END_REF]. To accomplish these objectives, it is essential to discover effective information handling steps. We selected three unsupervised band selection methods on the basis of their effectiveness in remote sensing that include Principal Component Analysis (PCA) [START_REF] Jia | Segmented principal components transformation for efficient hyperspectral remote-sensing image display and classification[END_REF], Maximum Noise Fraction (MNF) [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF] and Independent Component Analysis (ICA) [START_REF] Wang | ‖Independent Component Analysis-Based Dimensionality Reduction with Applications in Hyperspectral Image Analysis‖[END_REF]. The goal here is to show significant information as distinctive as possible for high-class separability. PCA uses highest Eigen values to select the three principle components from hundreds of band. In MNF [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF] goal is to pick bands in a way that maximize their signal to-noise (SNR) ratio by minimizing the noise fraction using a linear transformation. MNF considers the spatial information of each band for covariance matrix that makes MNF much more effective and reliable to image quality. The third method used is ICA [START_REF] Hyvärinen | Independent Component Analysis[END_REF] that measures data, statistical independency using mutual information as a criterion. It captures data that can't be retained or preserved by second-order statistics based dimension reduction methods like PCA and MNF. These three methods are applied in 2 ways:

Scheme-I: a) True Color b) CIR

This scheme is designed in two ways namely true color and near infrared color. Here the original HSI is segmented into a) visible range (R, G, B) and b) near infrared (R, G, NIR). Then we apply three DR methods to each spectrum. From 9 selected bands further 3 uncorrelated bands are selected for tri-stimulus composite display. The resultant image for true color gives the natural appearance while CIR version is helpful to display vegetation in red color. Flowchart of scheme-I is shown in Fig. 1.

Scheme-II

In scheme-I we selected only specific spectrum for visualization which results is loss of important information present in other bands. In this method we use all bands for image visualization. Keeping focus on our goals, first we spectrally segment the original data cube into three groups: visible range, NIR and middle infrared. In order to get highly informative featured bands, we further sub-group each spectral group into three equal parts and applied PCA on each sub-group to get principle components (PC). Through correlation metric Eq. (1) [START_REF] Su | Hyperspectral Image Visualization Using Band Selection[END_REF] from all PC's we selected the three most distinct, highly structured maximum energy bands.

( ) ( ) ( ) ( ) ( ) ( )

( ) ( )

By adding MNF, PCA, ICA images, the fused image ( ) gives inconsistent visualization. We also apply CMF stretching on original HSI to get (

On the basis of Perceptual based class separability, we prefer to assign different weights (found empirically) for tri-stimulus display. Resultant image ( ) has the ability to discriminate different data. The proposed approach gives image with maximum achievable distinction in the RGB components and therefore, enhanced illustration and visual discrimination of end-members.

RESULTS AND DISCUSSION

For performance evaluation, we select 3 statistical assessment parameters given in [START_REF] Kotwal | Visualization of Hyperspectral Images Using Bilateral Filtering[END_REF] where Entropy describes amount of information in image calculated as Eq.( 4). Here is the probability density of the intensity level x in the image.

= ∑ ( ) (x) (4)
Mean describe average gray value and Variance ( ) measures the deviation of gray values of the pixels from the image mean which describes image contrast. After removing water absorption atmospheric effected bands, we performed the experiments on AVIRIS Cuprite (C) and PAVIA CENTER (P) and compared the result with different previous methods. First image data was collected over the Cuprite mining site where the five minerals of interest are exposed alunite (A), buddingtonite (B), calcite(C), kaolinite (K), and muscovite (M) labeled in Fig. 3(a). For better subjective assessment of minerals we zoom-in our result in Fig. 3 (b) that clearly shows different class data with different contrast which is not seen in any previous methodology. Qualitative results are given in Fig. 2 for 6 different methods (3 proposed and 3 state of the art). We can see that none of the images is able to discriminate the different class information present in ROSIS Pavia center image except of our proposed method (scheme-II) that displays image with special contrast, while other methods are not able to give such details. As original data contain vegetation our CIR scheme displays it in red color because infrared range is sensitive for vegetation. Quantitative assessment in Table-1 well proves that our results are highly informative as images with high variance have better contrast and appealing for visualization.

Fig 1 :

 1 Fig 1: Scheme-I (a) True Color (b) CIR Display

Fig. 2 :

 2 Fig.2: Row-wise: AVIRIS Cuprite, ROSIS Pavia Centre Column from left to right: Segmented PCA [7], CMF Stretched, SBV, Scheme1a, Scheme 1b, Scheme 2 Optimized Class Separability

Fig. 3 :

 3 Fig.3: (a) Spatial position of five minerals [15] (b) Result by proposed scheme II

Table 1 :

 1 Quantitative Results for Cuprite (C) and Pavia (P)

	Method	Entropy	Variance	Mean
		C	P	C	P	C	P
	SPCA	6.73 6.46 633.8	434.2 96.33 98.55
	CMF	6.31 7.28 670.5 820.0 78.53 75.76
	SBV	7.15 7.32 594.4	794.4 95.00 90.87
	Scheme-I(a) 6.42 7.01 534.3 598.34 93.23 95.32
	Scheme-I(b) 6.64 7.11 523.8 543.69 81.25 87.04
	Scheme-II 7.45 7.67 920.2	849.5 98.37 99.55


Stefan Kinauer

Riza Alp Güler

Siddhartha Chandra

Iasonas Kokkinos

Structured Output Prediction and Learning for Deep Monocular 3D Human Pose Estimation

Keywords:

In this work we address the problem of estimating 3D human pose from a single RGB image by blending a feed-forward CNN with a graphical model that couples the 3D positions of parts. The CNN populates a volumetric output space that represents the possible positions of 3D human joints, and also regresses the estimated displacements between pairs of parts. These constitute the 'unary' and 'pairwise' terms of the energy of a graphical model that resides in a 3D label space and delivers an optimal 3D pose configuration at its output. The CNN is trained on the 3D human pose dataset 3.6M, the graphical model is trained jointly with the CNN in an end-to-end manner, allowing us to exploit both the discriminative power of CNNs and the top-down information pertaining to human pose. We introduce (a) memory efficient methods for getting accurate voxel estimates for parts by blending quantization with regression (b) employ efficient structured prediction algorithms for 3D pose estimation using branch-and-bound and (c) develop a framework for qualitative and quantitative comparison of competing graphical models. We evaluate our work on the Human 3.6M dataset, demonstrating that exploiting the structure of the human pose in 3D yields systematic gains.

Introduction

Human pose estimation has made rapid progress thanks to deep learning, as witnessed by the improvements reported on large-scale benchmarks [START_REF] Toshev | Deeppose: Human pose estimation via deep neural networks[END_REF][START_REF] Pishchulin | Deepcut: Joint subset partition and labeling for multi person pose estimation[END_REF][START_REF] Insafutdinov | Deepercut: A deeper, stronger, and faster multi-person pose estimation model[END_REF][START_REF] Cao | Realtime multi-person 2d pose estimation using part affinity fields[END_REF][START_REF] Newell | Stacked hourglass networks for human pose estimation[END_REF][START_REF] Wei | Convolutional pose machines[END_REF][START_REF] He | Mask r-cnn[END_REF]. In this work we focus on the more challenging task of 3D human pose estimation from a single monocular image, which can have many applications in human-computer interaction, augmented reality, and can eventually lead to addressing generic 3D object pose estimation.

Recovering 3D information from a single 2D image is clearly ill-posed, given that different 3D scenes can project to the same 2D image. However, exploiting task-specific prior knowledge can increase the probability of the more plausible scenes. All leading approaches to 3D human pose estimation, such as [START_REF] Tome | Lifting from the deep: Convolutional 3d pose estimation from a single image[END_REF][START_REF] Chen | 3d human pose estimation= 2d pose estimation+ matching[END_REF][START_REF] Bogo | Keep it smpl: Automatic estimation of 3d human pose and shape from a single image[END_REF][START_REF] Mehta | Monocular 3d human pose estimation in the wild using improved cnn supervision[END_REF], rely on incorporating prior knowledge about the structure of the 3D human body. Two-stage approaches, e.g. [START_REF] Chen | 3d human pose estimation= 2d pose estimation+ matching[END_REF][START_REF] Bogo | Keep it smpl: Automatic estimation of 3d human pose and shape from a single image[END_REF][START_REF] Guler | Human joint angle estimation and gesture recognition for assistive robotic vision[END_REF], firstly detect joint positions in 2D and subsequently lift joints into 3D by relying on prior knowledge about the 3D human pose. The advantage of such approaches is that they can exploit large datasets constructed for the prediction of 2D landmarksthe disadvantage is that errors in the 2D stage can propagate to the 3D predictions and can often not be recovered from. Inherently 3D approaches [START_REF] Pavlakos | Coarse-to-fine volumetric prediction for single-image 3d human pose[END_REF][START_REF] Burenius | 3d pictorial structures for multiple view articulated pose estimation[END_REF] discretize the depth variable and train a CNN to score every possible combination of position and depth with respect to the presence of a joint -one can understand that the CNN learns to use the scale of the joint to guess its depth. This was recently shown in [START_REF] Pavlakos | Coarse-to-fine volumetric prediction for single-image 3d human pose[END_REF] to deliver results that are largely superior over previous 2-stage approaches. More recent works have delivered further improvements by fusing the 2D and 3D streams [START_REF] Tekin | Fusing 2d uncertainty and 3d cues for monocular body pose estimation[END_REF].

In directly regressing the pose from the input image, the aforementioned approaches do not explicitly impose constraints that exploit the dependencies between the human joints. In our understanding, what is missing from existing works is a method to exploit the structure of the human pose in a natively 3D setup. Authors in [START_REF] Tekin | Structured prediction of 3d human pose with deep neural networks[END_REF] acknowledge this deficiency of contemporary methods, and propose to use a stacked denoising autoencoder to learn these dependencies implicitly. Other approaches to combining structured prediction with deep learning have recently been successfully pursued in 2D human pose estimation e.g. [START_REF] Tompson | Joint training of a convolutional network and a graphical model for human pose estimation[END_REF][START_REF] Yang | End-to-end learning of deformable mixture of parts and deep convolutional neural networks for human pose estimation[END_REF], while current approaches to incorporating structure in feedforward CNNs for pose estimation rely on cascading, or stacking the outputs of CNNs in 2D [START_REF] Newell | Stacked hourglass networks for human pose estimation[END_REF][START_REF] Wei | Convolutional pose machines[END_REF], which can become prohibitive when done in 3D, due to the increased memory and computation load. In this work we develop novel techniques that allow us to 'explicitly' capture the dependencies between human joints via an energy function that consists of unary and pairwise terms, and thereby pursue this direction in the arguably harder 3D setting.

Our contribution consists in showing that one can combine a volumetric representation with a structured model that imposes constraints between the relative positions of parts. Rather than relying exclusively on a feed-forward architecture, we show that one can append a structured prediction algorithm that propagates information on a graph that represents the part positions, and still remain computationally efficient. In particular, we train a CNN to not only populate the 'unary' 3D score maps for each part, but also to provide estimates of the relative positions of parts in 3D. This coupling of the part positions results in an optimization problem, that we treat as yet-another layer of a deep network, and thereby add functionality to the network.

We can describe our contributions as (i) allowing for a high resolution in the estimation of the pose without increasing the computation/memory budget (Sec. 2.1) (ii) constructing the coupling term so as to make optimization tractable while working with a high resolution in the depth coordinates (Sec. 2.2) (iii) exploiting rich connectivity, i.e. allowing the edges in the underlying graph of the model to form loops, by using fast approximate inference that combines Branch-and-Bound with ADMM (Sec. 2.3) and (iv) using Deeply Supervised Network (DSN) [START_REF] Lee | Deeply-supervised nets[END_REF] training to accelerate training, by forcing the 2D and 3D inputs to structured prediction to reflect part of the solution, while also leveraging on larger 2D datasets during training (Sec. 2.4).

We evaluate our approach on the Human3.6M dataset. Our results indicate that blending local information about the 3D pose of parts with information about the 3D displacements of parts yields systematic improvements over approaches that rely on either of these two cues alone, while incurring a marginal computational overhead, in the order of a fraction of a second. Fig. 1: We consider the task of 3D human pose estimation from a single RGB image. Our approach involves a fully convolutional neural network that provides a 'bottom-up' estimate of the 3D positions of parts and their relative displacements, and a structured prediction layer that combines them into a coherent estimate of the pose. The whole architecture is trained end-to-end, allowing us to optimize the CNN outputs with the respect to the subsequent pose estimation algorithm.

Method

We start by formulating our approach in terms of a structured prediction problem, and then provide the details about the individual components of our proposed approach. We represent the pose Φ in terms of the concatenation of the 3D coordinates of N individual parts φ i Φ = {φ 1 , . . . , φ N }.

(

) 1
Given an image I, we score a candidate pose in terms of a graphical model that considers individual properties of parts, as well as properties of some of their pairwise combinations:

S I (Φ) = N i=1 U i (φ i ) + i,j∈E P i,j (φ i , φ j ), (2)
where U stands for unary and P for pairwise potentials, and E is the set of edges used in our graphical model. The unary and pairwise terms are delivered by the CNN, while the structured prediction layer couples the parts through the optimization of Eq. 2.

If we consider a generic cost function, this can be challenging even for simple cases, let alone for the 3D pose space we are working with. Our main technical contributions aim at making the construction and optimization of Eq. 2 tractable while still exploiting the structure of the output space.

Quantized Regression for Depth Estimation

One of the main challenges in constructing a volumetric CNN is that the amount of memory and computation scales linearly in the granularity of the depth quantization, requiring to tradeoff accuracy for speed/memory. The root of the problem is that the underlying quantity is continuous, but plain regression-based models may be neither sufficiently accurate, nor expressive enough to capture the uncertainty and multimodality of the depth value caused by depth ambiguity, or occlusion.

Instead, we follow recent successful developments in object detection [START_REF] Girshick | Fast r-cnn[END_REF][START_REF] Ren | Faster r-cnn: Towards real-time object detection with region proposal networks[END_REF], dense correspondence estimation [START_REF] Guler | Densereg: Fully convolutional dense shape regression in-the-wild[END_REF], and pose estimation [START_REF] Papandreou | Towards accurate multi-person pose estimation in the wild[END_REF] where a combination of classification and regression is used to attack the image-based regression problem. We use a first classification stage to associate a confidence value with a set of non-overlapping depth intervals, corresponding to a coarse quantization of the depth value. If we have N classes and a depth range of, say D units, the k-th class is associated with a quantized depth of q k = k D N . This however may be at a very coarse depth resolution. We refine this coarse estimate by combining it with the results af a regression layer that aims at recovering the residual between the ground-truth depth values and their quantized depth estimates.

As shown in Fig. 2 this strategy allows us to 'retarget' the voxels to 3D positions that lie closer to the actual part positions, without requiring the exhaustive sampling of the 3D space. In particular a voxel v lying at the k-th depth interval will become associated with a novel 3D position of part i,

p v i = k D N + r i (v)
, where r i (v) is the residual regressed by our network for the i-th part type at voxel v.

The value of the associated unary terms, U i (p v i ), is obtained in terms of the inner product between a joint-specific weight vector, w i and a feature vector extracted from the CNN's output at the 2D position associated with voxel v. Fig. 2: Unary 3D coordinates via quantized regression. To efficiently regress the unary 3D coordinates, we use a divide and conquer strategy. We begin by quantizing the 3D space into voxels. We estimate the score of each joint belonging to each of these voxels using a classifier. Finally we regress a residual vector per voxel which indicates the offset between the center of the voxel and the continuous 3D position of each joint. Left: Sigmoid function on classified voxels and regressed residual vectors (in black) for two joints. Right: Regressed residual vectors for all joints.

Efficient Optimization with Quadratic Pairwise Terms

Having described how the unary terms are constructed in our model, we now turn to the pairwise terms and the resulting optimization problems. The expression for the pairwise term in Eq. 2, P i,j (φ i , φ j , I) would suggest constructing a six-dimensional function. Instead, as in the Deformable Part Model paradigm [START_REF] Felzenszwalb | Pictorial structures for object recognition[END_REF], we use a pairwise term that penalizes devations from a nominal displacement µ i,j :

P i,j (φ i , φ j , I) = - 3 d=1 c d (φ d,i -φ d,j -µ d,i,j ) 2 , (3)
where the c d parameters allow us to calibrate the importance of the different dimensions. These parameters are forced to be positive, while the expression in 3 corresponds to the log-probability under an axis-aligned Gaussian model, centered at the predicted part position. We note that as in [START_REF] Chen | Articulated pose estimation by a graphical model with image dependent pairwise relations[END_REF][START_REF] Sapp | Cascaded models for articulated pose estimation[END_REF], µ i,j is image-dependent, and in our case is the output of a sub-network which is trained end-to-end. This enables us to capture dependencies between parts, where an estimate of their ideal displacement is combined with the local evidence provided by their unary terms. One important advantage of the pairwise terms is that since they encode the relative position of parts they are often easier to model, since e.g. the distance between human joints is much more predictable than the actual positions of the joints. As such they can simplify the overall problem.

Another crucial advantage of the particular form of the pairwise term is that by virtue of being in the form of a quadratic cost function, it can easily be bounded from above and below using interval arithmetic -in particular, we rely on the 3D Branch-and-Bound algorithm introduced recently in [START_REF] Kinauer | Monocular surface reconstruction using 3d deformable part models[END_REF] to efficiently search over optimal combinations of parts in 3D. A brute-force, dynamic programming-type algorithm for solving this task would require a quadratic number of operations, since it would need to compare pairs of points. Our implementation has a low-constant linear complexity for the construction of per-part KD-tree data structures, and logarithmic best-case complexity for the subsequent optimization. In practice optimization requires less than a tenth of a second on a CPU, while further accelerations could be obtained through GPU-based implementations.

Network Connectivity: from star-shaped to loopy graphs

The Branch-and-Bound (BB) algorithm we use for efficient inference only accommodates a star-shaped graph topology. This can be problematic if one wants to model human pose in terms of a tree-structured graph, or introduce loops to capture more constraints. For this we employ master-slave type approximate inference techniques that allow us to use BB for slave problems and coordinate them through a master. In particular we rely on the Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF][START_REF] Martins | Dual decomposition with many overlapping components[END_REF][START_REF] Boussaid | Fast and exact: Admm-based discriminative shape segmentation with loopy part models[END_REF] which matches the continuous nature of the pose estimation problem [START_REF] Boussaid | Fast and exact: Admm-based discriminative shape segmentation with loopy part models[END_REF]. The approach to subdivide difficult problems into smaller and easier ones has before been seen in [START_REF] Komodakis | Mrf optimization via dual decomposition: Messagepassing revisited[END_REF]. The authors introduced Dual Decomposition to optimize MRF-type energies, outperforming former state of the art of "tree-reweighted message passing" algorithms. Later works on ADMM like [START_REF] Boussaid | Fast and exact: Admm-based discriminative shape segmentation with loopy part models[END_REF] borrow from developments outlined in [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] to reach convergence in a lower number of iterations. Loopy graphs are subdivided into easier to handle trees and coordinated via a master problem, which turns out to be updating the dual variables.

The method we outline below uses approximate inference to obtain solutions in ω(T log N ) operations, where T is a low constant in the order of tens, N is the number of voxels, and log N is the cost of re-solving the slave subproblems. The ω (best-case) notation relates to the (exact) Branch-and-Bound algorithm, which also empirically has typically this performance. Even though the ADMM-based results are now only approximately optimal, the cost function being optimized reflects more accurately the problem structure, which can positively affect accuracy.

We consider the case where the set of graph edges in Eq. 2 corresponds to a graph with loops. Denoting by R ⊂ 1...K the subset of point indices belonging to more than one star graph, our optimization problem can equivalently be rewritten as follows:

max S(Φ) = N i=1 S i (Φ i ) s.t. Φ i (r) = u(r) ∀r ∈ R, (4)
where S i is a set of loop-free subproblems, defined so that S(Φ) = N i=1 S i (Φ i ) for a common solution Φ. The consistency is enforced by the 'master', to whom the 'slave' subproblems S i deliver their solutions Φ i -obtained through Branch-and-Bound. In particular a relaxation to the constraints is updated and used to reset the problem solved by the slaves -at each step the relaxation becomes tighter and at convergence consistency is guaranteed. Dual Decomposition relaxes the constraints in Eq. 4 by introducing a Lagrange Multiplier λ i (r) for each agreement constraint. ADMM augments this with a quadratic constraint violation penalty resulting in an augmented Lagrangian function:

A(Φ, u, λ) = N i=1 (S i (Φ i )+ r∈R λ i (r), Φ i (r) )- r∈R (( N i=1 λ i (r))u(r)- ρ 2 N i=1 (Φ i (r)-u(r)) 2 )
(5) where ρ is a positive parameter that controls the intensity of the augmenting penalty. The quadratic term ensures rapid convergence by acting like a regularizer of the solutions found across different iterations. To maximize the augmented Lagrangian, ADMM iteratively performs the following steps:

Φ t+1 i = argmax Φi A(Φ i , u t , λ t ) (6)
u t+1 = argmax u A({Φ t+1 i }, u, λ t ) (7) λ t+1 i (r) = λ t i (r) -ρ(Φ t+1 i (r) -u t+1 (r)) (8)
In words, the slaves efficiently solve their sub-problems and update the master about Φ i , then the master sets u t+1 (r), and the current multipliers λ t+1 i (r), and communicates them back to the slaves for the next iteration. Unlike [START_REF] Boussaid | Fast and exact: Admm-based discriminative shape segmentation with loopy part models[END_REF] who used dynamic programming to efficiently solve the slave problems, here we combine ADMM with the Branch-and-Bound algorithm. Interestingly, both of the additional terms contributed by the master problem to the slave problems, λ i (r)u(r), (Φ i (r) -u(r)) 2 can be easily bounded using interval arithmetic, allowing for a straightforward incorporation into the original Branch-and-Bound method. With these changes we have observed similar convergence behavior as the one reported in [START_REF] Boussaid | Fast and exact: Admm-based discriminative shape segmentation with loopy part models[END_REF]; In typically 15-20 (sometimes even less) ADMM iterations the slaves converge to a consistent pose estimate.

Deeply Supervised 2D-and 3D-Learning

We have observed substantial simplifications in the learning procedure by employing Deeply Supervised Network (DSN) [START_REF] Lee | Deeply-supervised nets[END_REF] training. In particular we use loss functions that directly operate on the unary and pairwise terms, before these are integrated through structured prediction. We empirically observed that this substantially accelerates and robustifies learning, by helping the network come up with good 'proposals' to the subsequent combination stage. As discussed in Sec.2.1, the unary coordinates are obtained by adding the quantization and regression signals. Rather than expect this result to be correctly obtained only by back-propagation from the last layer, we also associate a classification and regression problem with each 2D image position. We associate every pixel with a set of discrete labels corresponding to quantized depth values. For each joint we learn a different classification function; we consider a voxel as being positive if the respective joint is within certain proximity to the center of the 3D volume. We train this classifier using the cross-entropy loss. We also regress residual vectors between voxel centers and groundtruth joints using an L1 loss which is only active when a voxel is close enough to 3D landmarks. For the pairwise terms, we regress vectors that point from each 3d joint to others. Similar to the unary coordinates, we regress these quantities in a fully-convolutional manner. The smooth L1 loss for the pairwise offsets between a specific joint and the rest of the joints is only active on pixels within certain proximity to the specific joint.

Training with a Structured Loss Function

Having outlined our cost function and our optimization algorithm, we now turn to parameter estimation. Our graphical model is defined in Eq. 2, and the pairwise terms are described in Eq. 3. As outlined in the preceding sub-sections, our network generates the unary terms U i (p v i ), the nominal displacements µ i,j and the 3D coordinates φ i . In this section we describe training of all these parameters, as well as the calibration parameters c in Eq. 3, using a structured loss function [START_REF] Joachims | Cutting-plane training of structural svms[END_REF][START_REF] Pepik | Multi-view and 3d deformable part models[END_REF][START_REF] Boussaid | Fast and exact: Admm-based discriminative shape segmentation with loopy part models[END_REF] that reflects the geometric nature of the problem we want to address. Once our loss function is defined, back-propagation can be used to update all of the underlying network parameters. While authors in [START_REF] Zhang | Improving object detection with deep convolutional networks via bayesian optimization and structured prediction[END_REF] use an Intersection-over-Union (IoU) based structured loss for the task of detection, given that in this setup we have access to continuous ground truth values that naturally capture the underlying geometry of the problem, we opt for simplicity and use a more straightforward structured loss function. Given that Φ denotes the 3D coordinates for a candidate configuration of parts (Eq. 1), and Φ denotes the groundtruth 3D coordinates, we use the Mean Euclidean Distance, ∆( Φ, Φ) = 1 P P p=1 φ p -φp 2 as a loss for our learning task, penalizing the 3D displacement of our estimated landmarks from their ground truth positions. As in standard structured output prediction, we use this loss to induce a set of constraints in pose space:

S( Φ) > S(Φ) + ∆( Φ, Φ) ∀Φ, (9)
requiring that the score of the ground truth configuration should be greater than the score of any other configuration by a margin depending on how far the particular configuration is from the ground truth. Since this cannot hold in general, we introduce slack variables ξ: ξ(Φ) = max(S(Φ) + ∆( Φ, Φ) -S( Φ), 0). Thus, the slack variables represent the violations of the constraints in Eq. 9, and our goal here is to learn the model parameters that minimize the slack variables.

Standard training of structural SVMs [START_REF] Joachims | Cutting-plane training of structural svms[END_REF][START_REF] Pepik | Multi-view and 3d deformable part models[END_REF][START_REF] Boussaid | Fast and exact: Admm-based discriminative shape segmentation with loopy part models[END_REF] typically finds the most violated configuration given by Φ * = argmax Φ (S(Φ) + ∆( Φ, Φ) -S( Φ)) and tries to reduce the violation of this configuration by updating the model parameters appropriately via the cutting-planes or Franke-Wolfe algorithm. In this work we use the standard stochastic gradient algorithm to minimize these slack variables. We do so by first finding K most violated configurations for each input sample (K is a hyper-parameter which affects the convergence speed; we set K = 20 based on experiments on a validation set). We then compute the sub-gradients of the model parameters with respect to each of these violated constraints and back-propagate them through the network.

3 Experimental Evaluation

Network Architecture

In our experiments we use a fully-convolutional 151 layer ResNet, with weights initialised from a model pre-trained on MPII for 2D body pose estimation [START_REF] Insafutdinov | Deepercut: A deeper, stronger, and faster multi-person pose estimation model[END_REF]. Both 3D and 2D branches of our network are implemented as single-level convolution layers branching from the last layer of the ResNet. The input images to the system are cropped and rescaled to a fixed size of 320x320; the downsampling factor of our network is 16, leading to a cube of 20x20x20 dimensions for 3D unary detection and residual regression branches and 20x20 spatial dimensions for the 2D branches.

Dataset

We use the largest available 3D human pose dataset Human3.6M (29) to train and evaluate our approach. The dataset consists of 3.6 million video frames of daily life activities performed by actors whose 3D joint locations are recorded by motion capture systems. Following the recent works in the literature, we have used frames from subjects S1, S5, S6, S7 and S8 for training and S9 and S11 for testing. We have used frames from all 4 cameras and all 15 actions in our training and testing in an action-agnostic manner.

We have sub-sampled the videos at 10 frames per second. Several videos that suffer from drift of the groundtruth joints are removed from the dataset. Due to the projective geometry, it is not possible to obtain "groundtruth data-cubes" from 3D poses. In particular, we cannot assume 3D points project to 2D points according to an orthogonal projection model. To cope with this issue, Pavlakos et al [START_REF] Toshev | Deeppose: Human pose estimation via deep neural networks[END_REF] create a data-cube using image coordinates for x and y dimensions and real-world coordinates for the z dimension (distances relative to the root node). At test time the depth of the root node and the intrinsic camera parameters are used to obtain 3D pose estimates.

Unlike their approach, which requires knowledge of the root node's z-coordinate at test time, we estimate z coordinates such that the ratio of standard deviations of real-world and projected coordinates in x, y dimensions is preserved in the z dimension. This approximation naturally introduces some reconstruction error, but leads to a system that estimates pose up-to a similarity transform agnostic to the distance of the person to the camera and the intrinsic camera parameters.

Joint Training with 2D Pose

Our network is initialized with ResNet parameters obtained by training for 2D joint localization on the MPII dataset, but we observe that including samples from MPII as training samples increases performance -apparently not doing so results in the network forgetting about 2D joint localization. As in [START_REF] Sun | Compositional human pose regression[END_REF] we modify the labelled joints of the Human3.6m dataset in order to be able to utilize the 2D data. In particular we include a joint of "thorax" between shoulders that is connected to the "neck" and discarding "chin" and "abdomen" joints. The resulting skeleton structure is identical to the one of MPII. We have verified that two identical networks trained with baseline and MPII-type label structures lead to equivalent evaluation scores, thus it is fair to compare to existing methods. The active losses for an MPII sample are 2D detection and X and Y pairwise offset values, while the 3D position estimates are ignored.

Results

Since our groundtruth comes in the form of projected coordinates, we can obtain the 3D pose only up-to a similarity transform. We report "reconstruction error", which is measured as the mean euclidean distance to the ground truth, after applying Procrustes analysis.

We experimented with a number of graph topologies and notice that performance depends on the graph structure: center star describes the graph topology where all joints are connected to one central root node at the human's torso. It performs better than "unary only", indicating that the body center "knows" something about the other body parts. stick figure is a graph that directly corresponds to the human skeleton, i.e. the wrist is connected to the elbow, the elbow is connected to the shoulder, and so on. Clearly the shoulder knows better where the elbow has to be than the root node in the torso. This structure clearly performs better than "center star". extended stick figure is an extension to "stick figure", containing all its edges plus additional connections between the elbows of left and right arm, left and right knee, head to shoulders and torso to knees. This shows that additional loops boost performance, stabilizing against outliers or false evidence. 2-hop follows the human skeleton like "stick figure" and adds connections from every joint to its indirect (2-)neighbours in the skeleton. This connects, for example, hand with shoulder and ankle to hips and left to right knee. "2-hop" performs best, helping to resolve occlusions and improving accuracy. Our experiments, reported in Table 1 clearly indicate that the 2-hop graph topology outperforms all of the other structures that we experimented with. This indeed justifies using approximate inference (ADMM), since these results require employing a loopy graph.

In Table .2 we compare the performance of our method to existing methods. Our results indicate that (a) our quantizatoin + regression-based unary network already delivers excellent results, at the level of the current state of the art. (b) Structured prediction yields an additional, quite substantial boost.

We note that there are some methods that only use a single camera or only S-11 frames as test samples and the rest of the videos for training. In order to compare our approach to such works, we present our results per camera and per subject in Tab.3. Our results show that we are also outperforming the very recent work of [START_REF] Sun | Compositional human pose regression[END_REF], who uses only S-11 as test set and obtains 48.3, which is inferior with respect to our S-11 result(45.91), even though we have not used S-9 for training.

We provide qualitative results in Fig. 3, demonstrating cases where the ADMM inference clearly increases the pose estimation performance. Figure 4 shows some example images from the LSP dataset in the left column, augmented with the inferred body skeleton. The other three columns illustrate the plausible 3D structure as inferred by our approach.

Conclusion

In this work we have introduced an efficient method for 3D human pose estimation from 2D images. To this end, we augment the functionality of existing deep learning networks by adding a final layer that optimizes an energy function with variables in three dimensions. We have shown our method to deliver state-of-the-art 3D human pose estimation results, and intend to explore ways of extending our method to other tasks, such as general 3D object pose estimation. Another avenue of investigation is to further expand our optimization approach to handle even more general pairwise potentials. Furthermore we intend to make our optimization module software publicly available where it can be easily adopted for other applications in the computer vision community.

Fig. 3 :

 3 Fig. 3: Examplar pose estimates by ADMM inference: Blue indicates the ground truth pose, whereas red and green is the solution obtained from "unaries alone" and ADMM respectively.

Table 1 :

 1 Comparison of average reconstruction errors for different graph topologies.

		Directions Discussion Eating Greeting Phoning	Photo	Posing Purchases
	UNARY alone	49.69	49.45	47.77	50.69	54.80	57.35	43.76	44.11
	center star	49.41	49.26	47.35	49.93	50.97	56.12	43.62	43.43
	stick figure	49.13	49.19	47.15	49.70	50.50	55.57	43.53	43.59
	extended stick figure	49.16	49.07	47.35	49.82	50.67	55.45	43.60	43.57
	2-hop	48.89	48.75	47.07	49.40	49.82	55.31	43.30	43.47
		Sitting	Sit. Down Smoking Waiting Walk Dog Walking Walk Tog. Average
	UNARY alone	65.39	95.76	53.53	46.27	51.53	41.59	49.52	53.48
	center star	61.50	78.09	52.51	45.88	50.63	41.08	49.41	51.42
	stick figure	60.14	79.46	51.52	45.74	50.59	40.73	49.33	51.12
	extended stick figure	59.94	78.51	51.42	46.01	50.39	40.89	49.32	51.08
	2-hop	60.48	78.20	51.69	45.63	50.16	40.74	49.17	50.87
					Average error
			Yasin et al. [36]		108.3
			Rogez et al. [37]		88.1
			Tome et al. [8]		70.7
			Pavlakos et al. [13] 3		53.2
			(Ours)Unary		53.48
			(Ours)ADMM		50.87

Table 2 :

 2 A comparison of our approach to methods that report reconstruction error in literature.

Table 3 :

 3 55.62 51.24 56.10 55.22 54.54 S-11 51.14 42.86 47.83 41.90 45.91 Average 53.72 47.64 52.59 49.57 Reconstruction errors for videos for specific cameras and test subjects in the Human 3.6M dataset.

	cam1 cam2 cam3 cam4 Average
	S-9

Acknowledgements: This work has been funded by the European Horizon 2020 programme under grant agreement no. 643666 (I-Support).