Photogrammetric, Geometrical, and Numerical Strategies to Evaluate Initial and Current Conditions in Historical Constructions: A Test Case in the Church of San Lorenzo (Zamora, Spain)

Abstract

Identifying and quantifying the potential causes of damages to a construction and evaluating its current stability have become an imperative task in today’s world. However, the existence of variables, unknown conditions and a complex geometry hinder such work, by hampering the numerical results that simulate its behavior. Of the mentioned variables, the following can be highlighted: (i) the lack of historical information; (ii) the mechanical properties of the material; (iii) the initial geometry and (iv) the interaction with other structures. Within the field of remote sensors, the laser scanner and photogrammetric systems have become especially valuable for construction analysis. Such sensors are capable of providing highly accurate and dense geometrical data with which to assess a building’s condition. It is also remarkable, that the latter provide valuable radiometric data with which to identify the properties of the materials, and also evaluate and monitor crack patterns. Motivated by this, the present article investigates the potential offered by the combined use of photogrammetric techniques (DIC and SfM), as well as geometrical (NURBs and Hausdorff distance) and numerical strategies (FEM) to assess the origin of the damage (through an estimation of the initial conditions) and give an evaluation of the current stability (considering the deformation and the damage).

1. Introduction

The conservation of built heritage is today considered a fundamental aspect of modern society. Their artistic, cultural, and intrinsic value make these constructions extremely important. Complementary to this, the lack of the building‘s own mechanical values and the characteristic behavior of its masonry, the complex interaction between components, and the lack of documentation, make the analysis of such constructions remarkably difficult. Currently, and derived from these considerations, numerous regulations propose the integration of different approaches among which are [1]: (i) the study of the construction’s history; (ii) inspection; (iii) monitoring; and (iv) structural analysis.

Regarding the numerical calculations, the static graphic [2] and limit analyses [3] traditionally provided the necessary tools to study the stability and bearing capacity of historical structures [4]. However, such numerical strategies have among their drawbacks the difficulty to evaluate damages [1].

In contrast with these models, the Finite Element Method (FEM) has been widely used for the evaluation of historical buildings at different levels; from complex and large constructions through macromodelling techniques [5], to the use of micromodelling strategies [6], where the units are independently discretized, or homogenized [7]. However, the large number of involved variables, as well as interaction with other structures, conditions the results.

It is in the field of built constructions where remote sensors and especially photogrammetric and laser scanner systems have proven great worth for their analysis [3,8,9,10]. These sensors are able to provide accurate and dense geometric and radiometric values with which to assess these buildings, as well as obtaining the data through non-intrusive means. Despite this, the data they provide (in form of dense and accurate point clouds) is largely untapped, since it is only used for the construction of simplified CAD models [10].

On one hand, the present article introduces two novel robustness parameters (based on geometrical components) in order to increase their applicability, obtained from the symmetrical Hausdorff distance [11]. These parameters, called Global Hausdorff metric (GHm_s) and Local Hausdorff metric (LHm_s), help ascertain whether the variables or simulated conditions improve or worsen the numerical results, in comparison with the real deformation provided by the photogrammetric and laser scanner systems.

On the other hand, the article introduces a methodology based on a Non-Uniform Rational B-Splines (NURBs) modelling strategy, with the purpose of providing an accurate geometrical model (with the current deformation and damage) for the evaluation of the current stability of the construction. This strategy is able to take advantage of the properties provided by the Structure from Motion products: (i) density; (ii) accuracy; and (iii) photorealistic texture, within a numerical environment.

In order to confirm the feasibility of the proposed geometrical strategies (GHm_s, LHm_s and NURBs modelling), they are applied to a case study: the dome of the church of San Lorenzo in Sejas de Aliste (Zamora, Spain). This construction, built in brick masonry, has suffered severe structural damages, shown through significant deformation, cracking and plastic hinges that reduce its bearing capacity. It seems necessary to perform a structural evaluation in order to design efficient restoration actions.

The article is organized as follows: Section 1 consists of an introduction and brief state of the art, Section 2 describes the different image-based techniques that were employed; Section 3 is made up of the description of the construction, the current deformation, damage, and the numerical aspect through the FEM; Section 4 describes two robustness indices based on geometrical discrepancies, a manual calibration of the model and a complementary strategy to evaluate the current stability of the construction (considering the complex geometry and the presented cracking); and finally, Section 5 shows the conclusions.

2. Image Based Approaches: Digital Image Correlation and Structure from Motion

The great diversity of approaches today, along with their flexibility, place image-based procedures as a suitable solution for the analysis of constructions [3,9,12], materials [13,14], and pathologies [8].

The different methodologies that comprise this approach, particularly in the field of numerical evaluation of constructions, highlight: (i) Digital Image Correlation (DIC); and (ii) image-based modelling procedures. While the former provides mechanical data of materials and constructive solutions (in the form of displacement and strains), the latter allows the definition of a dense, accurate, and photorealistic geometrical model of the construction. Their combination provides relevant information for the numerical analysis of the structure.

2.1. 2D Digital Image Correlation

A wide variety of methodologies has been developed and used to study material and union behavior. Some of these are [14,15]: (i) Moiré interferometry; (ii) Holography interferometry; (iii) Shearography; and (iv) Digital Image Correlation.

These methodologies prove to have important advantages, compared to traditional methods based on strain gauges or LVDT’s (Linear Variable Differential Transformer) such as their non-invasiveness and their full-field data information. In comparison, traditional methods provide only local information and require direct contact with the tested material. Within this wide range of techniques, the use of Digital Image Correlation (DIC) stands out.

To characterize the materials used in the dome, various compression tests were performed separately on each material (three in each material) during the experimental campaign that was carried out. Considering the procedure defined by [16], an extra specification, such as the mortar joint material (made by gympsum mortar), was considered.

In order to verify the flexibility and accuracy of the shown method, one standard sensor was used for both the DIC and the SfM: a digital reflex camera Canon 500D. However, in contrast to the image-based modelling strategy, DIC requires the preparation of the analyzed specimen, following the approach defined by [16]: (i) MIG (Mean Intensity Gradient) evaluation [17] of the speckle pattern; (ii) camera pose estimation [18]; and (iv) camera calibration [19].

Once the specimen has been correctly pre-processed, different images were captured during the test (Figure 1). Also, concerning the test setup, a large focal distance and working distance were used in order to minimize the geometrical distortion, out of plane displacements (approximation to a telecentric lens system) [20], depth of field, and light conditions.

Figure 1. (a) Detail view of the brick and speckle pattern applied during the Digital Image Correlation (DIC) test; (b) Histogram of the speckle pattern.

Figure 1. (a) Detail view of the brick and speckle pattern applied during the Digital Image Correlation (DIC) test; (b) Histogram of the speckle pattern.

The basic principle of DIC is the tracking (or matching) of the different areas of the images which were captured during the test (before and after deformation occurs), called subsets. As an initial approximation of this tracking, a correlation coefficient (generally the Zero Normalized Cross Correlation) [20] is used. Later, this initial approximation is optimized by the use of a non-linear strategy (such as the Inverse Compositional Gauss Newton method) [20] which allows the evaluation of the displacement suffered by the subset along the different captured images (Figure 2) [16].

Complemented to this optimization process an interpolation process (based on splines) is used with the aim of obtain sub-pixel accuracy [20]. Considering multiple subsets in the image, their analysis can provide a full-field displacement. Later, the strains suffered by the specimen during the test, which allow the evaluation of its mechanical properties, can be obtained by a direct relationship between the obtained displacement on the measurement point and the initial length of the virtual extensometers [16]. A total of three virtual extensometers were placed on the ROI: (i) A-A’ and B-B’ in the longitudinal direction; and (ii) C-C’ in the transversal direction.

Figure 2. Digital Image Correlation general outline. In red the reference subset, in blue the initial seed, and in yellow, the final location of the subset.

Figure 2. Digital Image Correlation general outline. In red the reference subset, in blue the initial seed, and in yellow, the final location of the subset.

Concerning accuracy, there are different studies [15,21,22] that endorse the DIC’s precision for the assessment of the material’s mechanical properties. For a standard test configuration, the accuracy may range from values of 0.01–0.04 pixels. Considering a conservative threshold at 0.1 pixels, and an acceptable accuracy for the test of 0.01 mm (from which the critical pixel size is set at 0.1 mm), the test’s configuration is shown in (Table 1).

Table 1. Summary of the different properties set during the Digital Image Correlation (DIC) test carried out with a Canon EOS 500D and a macrolens system 70–300 mm.

Values Adopted during the Digital Image Correlation (DIC) Test
Aperture	7.1
Focal length (mm)	200
Working distance (mm)	2700
Pixel size (mm)	0.063
Acquisition frequency (Hz)	0.33

Table 1. Summary of the different properties set during the Digital Image Correlation (DIC) test carried out with a Canon EOS 500D and a macrolens system 70–300 mm.

Values Adopted during the Digital Image Correlation (DIC) Test
Aperture	7.1
Focal length (mm)	200
Working distance (mm)	2700
Pixel size (mm)	0.063
Acquisition frequency (Hz)	0.33

Once the stress-strain curve has been obtained (by a relationship between the stress applied by the compression press and the strains obtained by DIC) (Figure 3b), it is possible to extract the mechanical properties of the materials as follows (Table 2): (i) the Young Modulus was considered as the ratio between one third of the maximum force achieved and the mean strain provided by the longitudinal extensometers (A-A’ and B-B’); (ii) for the Poisson ratio, the relationship between the strains provided by the longitudinal extensometers and those obtained by the transversal extensometer (C-C’) was taken into account; and (iii) the compression strength was considered as the maximum pressure supported by the specimen.

Figure 3. Results after the experimental campaign (2D DIC). (a) Deformation measurement, expressed in pixels, between two captures and positioning of the virtual extensometers; (b) Stress-strain curve obtained with the virtual extensometer A-A’.

Figure 3. Results after the experimental campaign (2D DIC). (a) Deformation measurement, expressed in pixels, between two captures and positioning of the virtual extensometers; (b) Stress-strain curve obtained with the virtual extensometer A-A’.

Table 2. Mechanical properties obtained by the performed DIC test.

Mechanical Properties Obtained by the DIC Test
	Clay brick	Gypsum mortar
E (GPa)	3.10 ± 0.30	1.15 ± 0.06
ν (-)	0.22 ± 0.05	0.23 ± 0.02
fc (MPa)	7.80 ± 0.90	2.12 ± 0.10

Table 2. Mechanical properties obtained by the performed DIC test.

Mechanical Properties Obtained by the DIC Test
	Clay brick	Gypsum mortar
E (GPa)	3.10 ± 0.30	1.15 ± 0.06
ν (-)	0.22 ± 0.05	0.23 ± 0.02
fc (MPa)	7.80 ± 0.90	2.12 ± 0.10

To assess the accuracy of the previously mentioned procedure, a comparative study was carried out, between the strain rate applied by the compression press (with an average value of −1.77 × 10⁻⁶ s⁻¹), and the one obtained by the different performed DIC tests and the different virtual extensometers used (with an average value of −1.93 × 10⁻⁶ s⁻¹). The results obtained demonstrate the accuracy and suitability of the applied configuration and algorithms, with an estimated precision of 0.056 pixels (which correspond to an approximate value of 3.53 μm). This value proves to be lower than the previously shown critical value.

2.2. Image-Based Modelling: Structure from Motion

In recent years the image-based modelling strategy, called Structure from Motion (SfM), has positioned itself as an attractive alternative to laser scanning systems. Its flexibility—as it can be integrated into different types of platforms (e.g., UAV [9])—low-cost, and qualities of the point cloud (high density, photorealistic texture, and accuracy) place the solution at a vantage position in the evaluation of historical buildings [23].

This technique integrates within its operating structure the advantages of computer vision (automation and flexibility) and photogrammetry (accuracy and reliability) [23] to obtain high density three dimensional models whose accuracy can compete with those of the laser scanner system [24,25].

For this case study, a standard SfM strategy is applied, comprising the following stages: (i) automatic extraction and keypoint matching by applying the Affine-Scale Invariant Feature Transform (ASIFT) algorithm [26]; (ii) automatic hierarchical orientation of images; and (iii) dense model generation through the MicMac algorithm. For further details on this methodology see [12]. Concerning the photogrammetric network a convergent protocol was used, combining a total of 32 cameras with high overlap (around 90%) and throwing a mean GSD (Ground Sample Distance) of 1.61 mm. Complementary to these, different circular targets (along the lower part of the pendentive) were used to scale the model (this measurement were taken by a total station using a radiation approach).

As a result of the implementation of the above-mentioned methodology it is possible to obtain a dense and photorealistic texture point cloud (Figure 4a). Afterwards, applying CAD conversion techniques (meshing, surface parameterization, extrusions, revolutions, etc.), or even generating true-ortho-images, increase further the applicability of the obtained product. More precisely, they help to accurately build CAD models suitable for subsequent numerical simulations, as well as complementary products, which analyze patterns of deformation and cracking for the pathological characterization of the structure [9] (Figure 4b).

Concerning the total error, associated with this point cloud, a quadratic error propagation was used Equation (1). Into this approach, two sources were considered: (i) the error coming from the bundle adjustment of the photogrammetric network; and (ii) the error corresponding to the scaling process Equation (2).

$${\epsilon }_t=\sqrt{{{\epsilon }_p}^2+{{\epsilon }_s}^2}$$

(1)

$${\epsilon }_s=\sqrt{2{{\epsilon }_i}^2+{{\epsilon }_m}^2}$$

(2)

where ε_t represents the total error; ε_s the scale error; ε_p the error associated with the photogrammetric network; ε_i the origin error stablished as $\sqrt{2}*(\frac{pixelsize}{\gamma }),$ where γ is the subpixel accuracy of the target detection algorithm (estimated in 0.5); and ε_m the error associated with the total station.

As a result a budget error of ε_t = 4.38 mm was obtained (with values of ε_P = 3.22 mm, ε_i = 1.14 mm, ε_m = 2.50 mm, and ε_s = 2.97 mm).

Figure 4. (a) 3D model obtained by the proposed methodology; (b) Detail view of the most damaged section through the texture model.

Figure 4. (a) 3D model obtained by the proposed methodology; (b) Detail view of the most damaged section through the texture model.

3. Structural Evaluation of San Lorenzo’s Dome

3.1. San Lorenzo Church

The church is built with irregular masonry walls (slabs of slate) fixed with lime mortar and at the corners finished with granite masonry. The parish church of San Lorenzo in Sejas de Aliste is located in the region of Aliste, Zamora province (Spain), 32 m long and 17 m wide, it belongs to the family of temples with transept crossing, Latin cross-shaped floor plan, and transept and nave at different heights (Figure 5).

Figure 5. San Lorenzo church: (a) Orthophoto of the main façade through the methodology proposed; (b) Orthophoto of the west façade (chancel) of the construction; (c) Floor plan-view of the church, red color indicates the damaged area of the dome.

Figure 5. San Lorenzo church: (a) Orthophoto of the main façade through the methodology proposed; (b) Orthophoto of the west façade (chancel) of the construction; (c) Floor plan-view of the church, red color indicates the damaged area of the dome.

The transept crossing is the most representative element of this temple. Its importance in the building is highlighted in the interior through the semi-elliptical dome that shelters the whole crossing. Its eight ribs marked with bands stand out. The transept is highlighted in the outside as well, covering the dome with a hipped roof that rises above the nave and transept height. This roof is built with a pavilion-shaped chestnut-timber framing, with regularly placed rafters that lean on the main beams and bear the load of the roof, made up of curved tiles and wooden roof boards (Figure 6a). Overall stability is obtained by use of tie beams at the top of the bearing walls, which collect the loads of the rafters and the hip rafters or main beams. Angle-ties, placed at 45° in each corner under the hip rafters prevent the transversal deformation of the tie beams (Figure 6a).

Figure 6. (a) Constructive section of the church’s transept; (b) Transversal section of the dome geometry (initial state estimated by the Structure from Motion (SfM) point cloud) with dimensions in meters.

Figure 6. (a) Constructive section of the church’s transept; (b) Transversal section of the dome geometry (initial state estimated by the Structure from Motion (SfM) point cloud) with dimensions in meters.

Concerning the dome, the construction has an estimated diameter of 6.72 m and a total height (measured from the pendentive) of 2.63 m. This structure was built with traditional tile brick and gypsum mortar, reaching a total thickness of 5.00 cm. It resembles, from a construction point of view, a Catalan vault (Figure 6b). It is also worth mentioning the presence of an infill (basically composed of a mixture of sand, clay and fragments of bricks with a medium compaction) at the support of the dome. This infill reaches a total height of 0.75 m and an average thickness of 0.65 m, and its presence contributes to the stability of the construction.

3.2. Present Damage and Deformation

The characterization of both, deformations and cracking patterns, is key to understanding the structure in terms of stability and safety. The high density, accuracy, and photorealistic texture of the point cloud obtained by the proposed methodology (Section 2.2) can address this task foregoing any need of physical contact with the structure. Through evaluation of the obtained product, it is possible to obtain a hypothesis for the origin of the damage.

It is worth noting that there is widespread damage in the area enclosed by three ribs (corresponding to the southern part of the dome). This area has two main cracks, in the parallel direction, which are interconnected through the presence of two plastic hinges. At its maximum, there is a deflection of 19.70 cm (compared to the initial estimated model) (Figure 7).

These structural pathologies seem to be related to the presence of assymetric loads acting on part of the dome’s shell (Figure 7a). More specifically the current damages, which are located under the south wing, can be attributed to a failure of the timber structure.

On one hand, the evaluation of photogrammetric products (which are the result of the previously defined SfM strategy) allows an estimation of the possible causes of the dome’s damage. However, it is required to have numerical strategies to verify these assumptions and assess the current state of the construction. For the present case study, and considering the hypothesis of failure of the timber structure several numerical analyses were performed: (i) numerical evaluation of the timber structure for the worst load case: snow; (ii) evaluation of the dome’s stability under self-weight; and (iii) numerical evaluation of the interaction timber-dome as a result of a timber failure.

Figure 7. Results of the visual inspection over the different photogrammetric products: (a) Surface comparison between the initial proposed model and the most deformed one estimated by the SfM point cloud; (b) Damage inspection in the orthophoto, in green the main observed cracks, in blue the secondary cracks, in yellow the material removal.

Figure 7. Results of the visual inspection over the different photogrammetric products: (a) Surface comparison between the initial proposed model and the most deformed one estimated by the SfM point cloud; (b) Damage inspection in the orthophoto, in green the main observed cracks, in blue the secondary cracks, in yellow the material removal.

Previous investigations carried out by [16], verified the stability of the timber structure for the most adverse load case: the presence of snow. Yielding a maximum deflection of 2.35 cm, it proves to be insufficient in order to interact with the dome. Considering these results, it is possible to conclude that the interaction between the cover and the dome seems to be linked with the presence of pathological agents (mainly moisture and biological organisms) which reduce the bearing capacity of the timber structure until it fails and rests on the dome.

3.3. Numerical Simulation of the Initial State of the Dome: Self-Weight and South Wing Support

Understanding the degradation mechanisms present in the construction requires a geometrical model of its initial state, a material characterization, as well as its boundary conditions, and load assessment. According to these, several numerical simulations (through non-linear static analysis) were performed in order to understand the causes and the construction’s initial conditions. Several improvements, regarding the geometric and mechanical aspects, are introduced in comparison to the previous investigations performed by [16]: (i) consideration and modelization of the infill-dome interaction; and (ii) account of gypsum as union material.

Regarding the mechanical aspect, a macromodelling strategy of the masonry was followed. This technique blends the bricks, mortar joints as well as the brick-mortar interface into one continuum assuming homogeneous material properties (Table 3). Also, the recommendations exposed by reference [6,27] were considered. An initial estimation of the masonry’s Young modulus was estimated using the formulas displayed in [28], setting the initial Young modulus at 2.54 GPa. However, further visual inspections showed the presence of an erratic masonry with low overlap between units. In accordance with this, a reduction of the initial Young modulus was considered for subsequent simulations (half of the initial estimated), yielding a final value of 1.22 GPa, analogous to those used in similar studies [7,29].

Table 3. Mechanical properties adopted for the macromodelling of the masonry.

Mechanical Properties for the Masonry Structure
Em (GPa)	Young Modulus	1.22
δm (kg/m3)	Density	1800.00
νm (-)	Poisson coefficient	0.25
ft,m (MPa)	Tensile strength	0.16
fc,m (MPa)	Compressive strength	1.60
dt,m (mm)	Ductility index in tensile	0.093
dc,m (mm)	Ductility index in compression	1.6
βc,m (-)	Shear retention	0.2

Table 3. Mechanical properties adopted for the macromodelling of the masonry.

Mechanical Properties for the Masonry Structure
Em (GPa)	Young Modulus	1.22
δm (kg/m3)	Density	1800.00
νm (-)	Poisson coefficient	0.25
ft,m (MPa)	Tensile strength	0.16
fc,m (MPa)	Compressive strength	1.60
dt,m (mm)	Ductility index in tensile	0.093
dc,m (mm)	Ductility index in compression	1.6
βc,m (-)	Shear retention	0.2

For the numerical simulation of the infill, a Morh-Coulomb failure criterion was considered, with its mechanical properties set according to the visual inspection (medium compaction) and the recommendations shown by [7,29,30] (Table 4).

Table 4. Mechanical properties adopted for the infill simulation.

Infill Mechanical Properties
Ei (GPa)	Infill Young Modulus	0.80
δi (kg/m3)	Infill density	1800.00
G (GPa)	Infill shear modulus	Ei/2
fi (MPa)	Cut-off tension	0.02
Φi (deg)	Infill friction angle	39
ci (MPa)	Infill cohesion	1 × fi

Table 4. Mechanical properties adopted for the infill simulation.

Infill Mechanical Properties
Ei (GPa)	Infill Young Modulus	0.80
δi (kg/m3)	Infill density	1800.00
G (GPa)	Infill shear modulus	Ei/2
fi (MPa)	Cut-off tension	0.02
Φi (deg)	Infill friction angle	39
ci (MPa)	Infill cohesion	1 × fi

Concerning the load (for the numerical evaluation of the interaction between timber structure and the dome), a value of 8000 N was considered, resulting from the combination of different loads: (i) 650 N/m² for the arabic tiles and wooden board; and (ii) 400 N/m² for the snow load. Finally, the numerical model (for both simulations), had a total of 46,181 high order solid elements (CTE30) [31] (Figure 8a).

It is possible to observe that in the absence of external loads acting on the dome, the structure seems to be stable under its own weight (Figure 8b). In spite of this, considering the support of the south wing (roof tile, boards, and rafters), the dome begins to present damage (cracking) and its deformation tendency (Figure 8c) seems to be similar to the one shown in the photogrammetric model (Figure 7a).

However, in terms of deformation, considering for this purpose six control nodes along the damage area (Figure 8a), the model exhibits high rigidity. This suggests that the initial mechanical conditions are inadequate to reproduce the damage and deformation presented in the dome (Table 5).

Table 5. Comparison between the obtained and expected displacement of the control nodes in the numerical model (initial considerations).

Control Node	Displacement Obtained (mm)	Displacement Expected (mm)
54	1.95	148.00
20256	0.48	46.00
56	1.57	198.00
21125	0.54	52.00
443	2.12	196.00
64123	0.57	25.00

Table 5. Comparison between the obtained and expected displacement of the control nodes in the numerical model (initial considerations).

Control Node	Displacement Obtained (mm)	Displacement Expected (mm)
54	1.95	148.00
20256	0.48	46.00
56	1.57	198.00
21125	0.54	52.00
443	2.12	196.00
64123	0.57	25.00

Figure 8. (a) Isometric view of the mesh and the control points (nodes) used for the numerical simulations; (b) First principal stress distribution, expressed in N/mm² for the self-weight case; (c) First principal stress distribution, expressed in N/mm², for the numerical model which considers the asymmetric load.

Figure 8. (a) Isometric view of the mesh and the control points (nodes) used for the numerical simulations; (b) First principal stress distribution, expressed in N/mm² for the self-weight case; (c) First principal stress distribution, expressed in N/mm², for the numerical model which considers the asymmetric load.

The high discrepancies shown in the previous numerical simulation suggest the need for an optimization of the mechanical properties. However, performing such an optimization requires inevitably having robustness indices to quantify the level of improvement/worsening introduced by the different variable’s variations.

Exploiting the advantages offered by the SfM or laser scanner systems, two novel robustness parameters (based on geometric discrepancies) are proposed: (i) a global parameter, based on the similarity between the numerical and real model; and (ii) a local index which provides data about the geometrical variations introduced by the new variables considered in different areas of the construction.

4. SfM, NURBS Modelling, Global and Local Hausdorff Metrics: Geometrical Strategies to Improve the Knowledge about the Initial and Current State of the Constructions

4.1. Global and Local Hausdorff Metrics as Geometric Accuracy Indices

The Hausdorff distance or Hausdorff metric is used in a wide range of fields, such as point cloud [32] and meshes [33] comparison, object recognition [34], and image comparison and matching [35]. This metric proves to be a robust strategy for the similarity evaluation of two compact and non-empty sub-sets within a metric space. It is formulated as follows Equations (3) and (4):

$$d(y,X)={\min }_{x\in X}{||y-x||}_2$$

(3)

$$d_H(Y,X)={\max }_{y\in Y}d(y,X)$$

(4)

where ||.||₂ stands for the Euclidean norm; min the minimum value (distance); max the maximum distance; X and Y are the two compact sub-sets defined by the numerical and photogrammetric nodes; and x and y the considered points inside these sub-sets.

It is worth mentioning that, considering the previously defined concept of Hausdorff distance, the value of the norm does not have a symmetrical nature; it is therefore different in each direction (d_H(X,Y) ≠ d_H(Y,X)). For that reason, the symmetrical Hausdorff distance d_SH Equation (5) is used as metric comparison to avoid potential errors of geometrical similarity. This way a more robust solution is provided for geometry comparison.

$$d_{SHi}=\max \{d_H(y,X),d_H(x,Y)\}$$

(5)

where d_SHi is the symmetrical Hausdorff distance; of sub-set i, between models (numerical and photogrammetric); and x and y are two points that respectively belong to sub-sets X and Y.

On the other hand, understanding the global structural behavior of the analyzed construction inevitably requires several numerical analyses in order to adapt the simulated behavior to the real one. It is necessary to take into account the consideration that new conditions or new values of variables may worsen or improve the global and/or local result of the structure. It is therefore possible to define, out of the previously shown comparison metric Equation (3), two novel geometrical indices of robustness that represent improvements or worsening in the new numerical simulations in comparison to a reference model, considering the different variations of the variables or conditions: Global Hausdorff metric Equation (6); and Local Hausdorff metric Equation (7).

$${GHm}_s=\left(\frac{{\displaystyle \sum _{i=1}^n}{d}_{SH}(\mathrm{i})-{\displaystyle \sum _{i=1}^n}{d}_{S{H}_b}(\mathrm{i})}{{\displaystyle \sum _{i=1}^n}{d}_{S{H}_b}(\mathrm{i})}\right)\times 100$$

(6)

$${LHm}_s(\mathrm{i})=\frac{d_{SH}(\mathrm{i})}{d_{S{H}_{ref}}(\mathrm{i})}$$

(7)

where GHm_s represent the Global Hausdorff metric index and LHm_s the Local Hausdorff metric index, d_SH(i) the symmetrical Hausdorff distance to cluster i considered for the model; d_SHb(i) the symmetrical Hausdorff distance for cluster i of the base model (the model that results from the geometrical discrepancies between the initial model and the photogrammetric one); and d_SHref(i) the symmetrical Hausdorff distance from cluster i to the reference one (which may be the base model).

On one hand, GHm_s is able to provide a global value, expressed in percentage, for the improvement/worsening of the numerical simulation model in comparison to the model that was considered as base model. On the other hand, LHm_s provides a comparison of the variations between the numerical model and the reference model at a local level (values lower than one indicates a local improvement and values higher than one, a worsening).

For this case study, the reference model was considered to be the base model, obtained by the application of Equation (3) between the photogrammetric model and the non-deformed numerical model (Figure 9).

Figure 9. Graphical distribution of the different considered symmetrical Hausdorff distance (d_SH) (expressed in m) for the base model.

Figure 9. Graphical distribution of the different considered symmetrical Hausdorff distance (d_SH) (expressed in m) for the base model.

Finally, and considering GHm_s and LHm_s as the robustness indexes, a manual calibration was carried out, according to the established lower and upper bounds (Table 6) (Figure 10). According to [36], which provides a range of mechanical properties for historical masonry constructions, the upper and lower bounds were established with a safety factor of 1.35, since nowadays only visual inspection and geometrical survey are available (without an extensive experimental campaign).

Table 6. Parameters and variables considered during the manual calibration stage.

	Variable	Initial Value	Upper Bound	Lower Bound	Update Value
ft,j (MPa)	Masonry tensile strength	0.16	0.20	0.05	0.13
Ei (GPa)	Infill Young Modulus	0.80	1.00	0.05	0.50
Em (GPa)	Masonry Young Modulus	1.22	0.89	1.33	0.90
β (-)	Shear retention factor	0.20	0.01	0.20	0.15

Table 6. Parameters and variables considered during the manual calibration stage.

	Variable	Initial Value	Upper Bound	Lower Bound	Update Value
ft,j (MPa)	Masonry tensile strength	0.16	0.20	0.05	0.13
Ei (GPa)	Infill Young Modulus	0.80	1.00	0.05	0.50
Em (GPa)	Masonry Young Modulus	1.22	0.89	1.33	0.90
β (-)	Shear retention factor	0.20	0.01	0.20	0.15

Figure 10. (a) First principal stress distribution, expressed in N/mm2 of the updated model; (b) Geometrical accuracy, in terms of Local Hausdorff metric (LHm_s) of the updated model; in green, values where the geometrical model improves the results, in orange values where no improvements are carried out and in red, areas where the updated numerical model displays worse behavior.

Figure 10. (a) First principal stress distribution, expressed in N/mm2 of the updated model; (b) Geometrical accuracy, in terms of Local Hausdorff metric (LHm_s) of the updated model; in green, values where the geometrical model improves the results, in orange values where no improvements are carried out and in red, areas where the updated numerical model displays worse behavior.

Noteworthy is the presence of a red area (Figure 10). The said phenomenon is associated with the presence of an offset in the spatial distribution of the plastic hinge in comparison to the photogrammetric one. Considering the results provided by the GHm_s and LHm_s indices (Figure 10b) (Table 7), a mild improvement in the geometrical similarity between the photogrammetric and numerical model (Figure 9b) is observable compared with the initial conditions (Figure 10a) and previous studies carried out on the dome [16] (presence of an infill, independent oculus, and manual calibration of the mechanical properties).

On one hand, the obtained numerical results, with a value of GHm_s of 7.40%, are insufficient to study the current stability of the dome based on an initial state model. The discrepancies, derived from the large number of currently unknown variables, call for the use of additional sensors as well as additional experimental campaigns (in laboratory and in field tests).

On the other hand, the causes of the current damage and deformation correspond to the initial one: a local failure of the timber structure (south wing) could be the cause of pathological agents acting on the wood (moisture and biological agents).

In order to understand the current stability of the construction it is required to evaluate it with the actual deformation and damage (cracks). Motivated by this, and given the geometrical and radiometric properties provided by the SfM systems, a geometric strategy is defined below.

Table 7. Comparison between expected and predicted displacement of the considered control nodes.

Control Node	Displacement Obtained (mm)	Displacement Expected (mm)
54	16.56	148.00
20256	22.80	46.00
56	26.10	198.00
21125	31.52	52.00
443	46.19	196.00
64123	38.92	25.00

Table 7. Comparison between expected and predicted displacement of the considered control nodes.

Control Node	Displacement Obtained (mm)	Displacement Expected (mm)
54	16.56	148.00
20256	22.80	46.00
56	26.10	198.00
21125	31.52	52.00
443	46.19	196.00
64123	38.92	25.00

4.2. Analysis of the Current Stability of the Construction Based on a SfM and NURBs Approach

It should be stressed that the structural evaluation of historical constructions not only implies the assessment of the damage’s causes, but also requires a thorough understanding of the current stability (considering the actual deformation and damage), in order to take efficient restoration actions on the construction and to predict its integrity in case of different events (e.g., earthquakes). With the aim of improving the knowledge of the current stability, with respect to previous studies (Section 3.3 and [16]), a new approach is needed.

Although the point clouds obtained by the previously defined SfM approach, rich in geometric (density and accuracy) and radiometric (photorealistic texture) features, accurately represent the actual state of the construction, it is required to have additional strategies capable of exporting these properties into a numerical environment. The resulting mesh (triangulation of the SfM point cloud) has significant shortcomings to be considered as a suitable CAD/CAM model. Among its deficiencies, the following stand out [37]: (i) High density/resolution, which implies a large number of triangles and (ii) inadequate shapes.

Under the said framework, a methodology able to exploit these features based on the Non-Uniform Rational B-Splines (NURBs) and enhanced by the integration of structural pathologies (such as cracks and lack of material) is proposed. It follows the workflow shown below (Figure 11).

Figure 11. Proposed workflow for the study of the current stability of the construction.

Figure 11. Proposed workflow for the study of the current stability of the construction.

Considering the point cloud as the starting point, this product is firstly meshed by a standard Delaunay triangulation. Usually, these meshes present a non-manifold structure, which implies a low quality product with non-natural triangles which hinder the NURBs’ generation. In order to minimize this drawback we use a topological reconstruction, which generates a maniflod mesh, based on the approach defined by [38].

Once the mesh has been correctly defined, a region clustering (boundary contours) was carried out, comprising two critical stages. In a first stage the boundary vertex (limits of the construction, lack of material and the absence of an oculus piece) of the mesh are extracted. Later, in a second stage, the cracks are integrated into these boundary contours, through a projection of the observed crack path on the SfM point cloud, as established [9].

After that, a correct representation of these regions by means of NURBS patches is required. For this purpose and in order to build a regular base on which to correctly estimate the parameters of the different regions, a quadrilaterization of the mesh is carried out. This procedure is based on the combination of Morse theory and Spectral mesh analysis according to [37]. This methodology guarantees a complete quadrilateral description of the mesh, with a C¹ (tangential continuity) between neighboring patches, ensuring a continuity along the edges.

Since the construction’s surface needs to be fitted using different NURBs regions (quadrilaterial patches), a regularization process of these regions is necessary. This procedure comprises several steps [37]: (i) selection of one random border in the considered path and its opposite; (ii) border’s regularization using B-Splines with a lambda density; and (iii) matching between points by means of the Fast Marching Method.

Finally, the points obtained by the regularization procedure are used as control points to fit each quadrilateral path to a NURBs’ surface. It is worth mentioning, that, in construction elements such as arches, vaults or domes, the acquisition of its extrados (as a point cloud) is not possible in most of the cases, due to the presence of several setbacks (e.g., presence of infill, poor lighting conditions or lack of accessibility). Therefore, this previously shown strategy restricted the analysis of such constructions with the membrane theory (Mindlin-Reissner or Kirchhoff-Love theories). These theories limit the geometry of the numerical model’s different elements to a minimum recommended size of ten times the construction’s thickness. It implies as well the waste of the geometrical potentialities offered by the previously shown methodology.

Based on what is remarked above, a complementary strategy is proposed; able to estimate the construction’s extrados based on its intrados geometry. This methodology is made up of the following stages: (i) decorative elements removal; (ii) normal estimation of the points by means of eigenvalue analysis of the covariance matrix [39]; (iii) translation of each point along the normal direction (with a value equivalent to the construction’s thickness); (iv) point cloud meshing based on the Poisson approach and (v) projection of the cracks along its orthogonal direction. As a result, an accurate geometrical model of the construction is obtained with which to evaluate its actual stability (Figure 12).

Regarding the numerical aspect, and for the present case study, an incremental static non-linear FEM was carried out [40]. The material properties and the modelling strategy remain the same as those estimated for the initial model (considering the most appropriate ones) in Section 3.3, including a discrete model of the cracks.

For the present case study, this cracking is modelled considering the residual transversal stiffness (shear strength) through Equation (8). Concerning the normal stiffness, only a contribution in compression was considered, dismissing any contribution to the tensile regime according to Equation (9) (Table 8).

$$K_{t,c}=\frac{G_b{G}_m}{h_m(G_b-G_m)}{\beta }_{crack}$$

(8)

$$K_{n,c}=\frac{E_b{E}_m}{h_m(E_b-E_m)}$$

(9)

where K_t,c and K_n,c represent the tangential and normal stiffness respectively; G_b and G_m the shear modulus of brick and mortar, respectively; E_b and E_m the Young modulus of brick and mortar, respectively; h_m the mortar thickness; and β_crack the shear retention factor.

Figure 12. (a) Isometric view of the considered mesh model; (b) Discrepancies, expressed in mm, between the Non-Uniform Rational B-Splines (NURBs) and the photogrammetric models.

Figure 12. (a) Isometric view of the considered mesh model; (b) Discrepancies, expressed in mm, between the Non-Uniform Rational B-Splines (NURBs) and the photogrammetric models.

Table 8. Mechanical properties considered for the interaction between macroblocks (cracks).

Mechanical Properties of the Cracks
hm (mm)	Mortar thickness	15.00
Gb (N/mm2)	Brick’s shear modulus	1.27
Gm (N/mm2)	Mortar’s shear modulus	0.47
βcrack	Shear retention factor	βc,m
Kt,c (N/mm3)	Tangential stiffness	121.88
Kn,c (N/mm3)	Normal stiffness (compression)	49.74

Table 8. Mechanical properties considered for the interaction between macroblocks (cracks).

Mechanical Properties of the Cracks
hm (mm)	Mortar thickness	15.00
Gb (N/mm2)	Brick’s shear modulus	1.27
Gm (N/mm2)	Mortar’s shear modulus	0.47
βcrack	Shear retention factor	βc,m
Kt,c (N/mm3)	Tangential stiffness	121.88
Kn,c (N/mm3)	Normal stiffness (compression)	49.74

Finally a mesh for the numerical simulation is provided, with a total of 45,350 elements, clustered in: 45,196 high order solid elements and 154 high order interface elements.

For the stability analysis, all the loads acting on the dome (self-weight, infill pressure and asymmetric load) were considered. Afterwards, the estimated safety factor was established as the ratio between the current load and the collapse load obtained in the numerical simulation (Figure 13).

According to the study carried out in the Section 4.1 and the inspection of the SfM point cloud, the collapse mechanisms are mainly due to the formation of plastic hinges in the tensile regime.

The complexity of the model and the uncertainties associated with the variables (e.g., soil properties) require the study of the influence of different mechanical variables in the global stability of the construction, through parametrics analysis.

For these analyses, only the most important mechanical properties (to the tensile regime) were considered, namely: (i) Young modulus; (ii) tensile strength; and (iii) shear retention factor.

It can be observed, that the stability of the dome is mainly conditioned by the mechanical properties of the masonry, rather than the mechanical properties of the infill. Therefore, a safety factor (considering the initial mechanical properties) of 1.23 was established (Figure 13c).

However, it is worth mentioning that only the most important cracks were taken into account. By following a discrete strategy, minor and diffuse cracks were not considered. For this reason further investigation, integrating complementary approaches is necessary in order to obtain a better estimation of the actual stability.

Figure 13. (a) Parametric analysis of different tensile strengths and shear retention factors; (b) Parametric analysis of different masonry and infill’s Young modulus; (c) Maximum principal stress (σ₁), expressed in N/mm², at collapse of the initial considered model.

Figure 13. (a) Parametric analysis of different tensile strengths and shear retention factors; (b) Parametric analysis of different masonry and infill’s Young modulus; (c) Maximum principal stress (σ₁), expressed in N/mm², at collapse of the initial considered model.

5. Conclusions

Based in the already established photogrammetric techniques of Digital Image Correlation (DIC) and image-based modelling (SfM), and complemented with geometrical (NURBs modelling and Hausdorff distance) and numerical methodologies (FEM), the strategies defined and used in the article allow the needs of structural evaluation of historical constructions to be met.

On one hand, two novel geometric quality indices are introduced and defined, called Global Hausdorff metric or GHm_s and Local Hausdorff metric or LHm_s. They allow to assess globally (GHm_s) and locally (LHm_s) the robustness, in geometric terms, of the obtained numerical model in comparison to the point cloud (deformed shape) of the construction. These indices can calibrate the different variables, based on the geometrical similarity between models acting on the numerical simulation.

On the other hand, with the aim of evaluating the actual stability of the construction and exploiting the geometrical and radiometric components of the obtained products (SfM point clouds), a modelling strategy based on NURBs is proposed. This strategy is able to profit from these properties to obtain an accurate geometrical model (with the actual deformation and damage), that serves as a basis for subsequent numerical analysis.

In order to validate these parameters and modelling methodology, it was applied to a real case study: the dome of the church of San Lorenzo in Sejas de Aliste (Zamora, Spain). Several simulations were carried out to understand the degradation process between the initial and the current state, and to corroborate the viability of the defined robustness parameters with a value of 7.40%, for the GHm_s. When studying the current construction’s stability, through the modelling strategy defined in the article, the results reveal a damaged construction with an estimate safety factor of 1.23.

However, the complexity of the model, the initial state, the absence of comprehensive knowledge of the different construction stages and the need of more experimental campaigns hinder the numerical results and the correct estimation of the safety factor. Taking this into account, further research will focus on the following aspects: (i) dynamical tests; and (ii) a robust calibration procedure (e.g., Non-Linear Square Minimization) based on the geometrical indices defined to enhance the numerical simulation of the dome.

Complementary to this, concerning the used image-based procedures, their potential includes: (i) flexibility (these may be used in the evaluation of mechanical properties of materials as well as geometrical models); (ii) wide range of applications, its use may be extended to other types of constructions such as tunnels or bridges; (iii) non-contact and non-destructive techniques; (iv) low associated cost; and (v) abundance of geometric and radiometric data. However, the methodology has some limitations: (i) the lack of geometrical information in non-visible areas, requiring complementary sensors such as electric tomography or ground penetration radar; and (ii) the model’s accuracy, with several millimeters of error, restricting this strategy to constructions with large deformations.