NON-CONTACT PHOTOGRAMMETRIC METHODOLOGY TO EVALUATE THE STRUCTURAL HEALTH OF HISTORICAL CONSTRUCTIONS

Accurate studies of cultural heritage structures usually require the application of combined techniques to understand its structural behaviour. It is in this context, where the article present a set of procedures based on the dual combination of photogrammetric methodologies (image-based modelling and digital image correlation) and the finite element method to understand the structural behaviour of these structures. Through the interpretation of the different obtained products, by the defined approach, it is possible to estimate and evaluate the causes of the structural damage that the analysed infrastructure can suffer and also design subsequent restoration mechanisms, always from a perspective of minimal intrusion to the structure. This case of study aims to validate the methodology presented, using as an example a case of vernacular architecture (a dome) in Sejas de Aliste (Zamora, Spain). * Corresponding author. E-mail address: luisj@usal.es (L. J. Sánchez-Aparicio).


INTRODUCTION
The preservation of architectural heritage is considered a fundamental issue in the cultural life of modern societies.The digitalization and modelling of buildings relevant to the cultural heritage through non-destructive techniques is the first step to form the basis for any kind of decision regarding intervention on cultural heritage.However, their structural analysis is sometimes an aspect that is not considered or only studied separately, being essential for their safeguarding and maintenance.
The conservation of these infrastructures requires understanding their structural behavior and consequently: (i) their boundary conditions,(ii) the mechanical properties of the materials,(iii) the origin of the damage that building suffers and (iv) their vulnerability (Sánchez-Aparicio et al., 2014).
In this process of preserving of cultural heritage structures, procedures based on non-destructive methods are particularly attractive, since they suppose minimum intervention and preservation of the original constructions.It is in this field where geomatic sensors such as terrestrial laser scanner (Armesto-González et al., 2010;Villarino et al., 2014a) and/or digital cameras (Riveiro et al., 2013;Villarino et al., 2014b) have acquired important roles, due to their capacity to acquire accurate geometry without any contact with the structure.However, the terrestrial laser scanner entails some drawbacks such as the difficulty for transport and the lack of flexibility during data acquisition.Wherein this situation photogrammetry can offer several solutions for these drawbacks.Also it is remarkable, in the case of image-based modelling, that this data acquisition system is complemented by the digital image and computer vision approaches, allowing to progress from 2D (images) to 3D (points cloud) in a wide range of applications within the field of cultural heritage structures preservation.
Thanks to the advances in the integration of photogrammetry and computer vision, as well as the numerical methods of finite elements, it is possible to understand the actual behavior of historical constructions with high accuracy.Taking all of this into account, this paper attempts to demonstrate the strong potentialities of image-based modelling, digital image correlation and finite element methods to evaluate the structural behavior and the origin of the damage in these structures.This paper is organized as follows: Section 1 motivates the importance and complexity of assessing structural health of cultural heritage buildings together with a brief state-of art; Section 2 presents two photogrammetric approaches: imagebased modelling and digital image correlation to obtain accurate data for structural analysis; Section 3 is devoted to the preprocessing methodologies to obtain accurate CAD models suitable for finite element simulations, the numerical simulation and also the post-processing strategies to evaluate the robustness of the results; and finally in Section 4, the conclusions are drawn.

Image-based modelling
Nowadays, the tridimensional reconstruction of objects is generally performed using two types of sensors: in one hand, range sensors such as laser scanner or structure light systems; on the other hand, the image-based modelling, namely Structure from Motion (SfM).This approach, the SfM, combines the accuracy and reliability of photogrammetric methods and the automation and flexibility advantages of computer vision (Sánchez-Aparicio et al., 2014).
In the present case of study, an in-house software PW (Photogrammetry Workbench) was developed, which encloses photogrammetric and computer vision algorithms, in a SfM environment.This software applies a multi-stage methodology composed for the following steps (García-Gago et al., 2014): (i) automatic extraction and matching keypoints; (ii) automatic hierarchical orientation of images and (iii) dense surface generation.

Automatic extraction and matching keypoints
The field of architectural documentation, in which often, complex scenarios need to be modeled, classic photogrammetric methods of correspondence as ABM (Area-Based Matching) or LSM (Least Squares Matching), are insufficient, due to the variations in scale, perspective and lightning (García-Gago et al., 2014).It is imperative therefore the use of more sophisticated strategies.
It is in this field where algorithms such as: (i) SUSAN (Smallest Univalue Segment Assimilating Nucleus); (ii) SIFT (Scale-Invariant Feature Transform); (iii) MSER (efficient Maximally Stable Extremal Region) or (iv) SURF (Speeded Up Robust Features), have acquired important roles.However, these new strategies are not unaffected by the perspective variations of different images.In this respect, a variation of SIFT, the ASIFT (Affine Scale-Invariant Transform), has been used and added to the PW tool.
Finally, the matching process is carried out by the employment of the SIFT descriptors, according to the Euclidean distance and secondly filtered by the ORSA (Optimized Random Sampling Algorithm) (Moisa et al., 2012) approach.

Automatic hierarchical orientation
After the extraction and matching of the aforementioned features, the orientation of the images is carried out by a multistage approach.Firstly, a relative orientation of the images is calculated, through the fundamental matrix (Longuet-Higgins, 1987).Secondly, a bundle adjustment of all images is performed by an iterative and least-square process based on the collinearity condition (Tait, 1994).
In cases in which the internal calibration parameters of the camera are unknown, the last step allows the calculation of these parameters, a self-calibration, through the incorporation into the equation as unknowns (Quan, 2010).

Dense surface generation
Given accurate exterior orientation parameters and a sparse point cloud, generated from the keypoints extracted, a robust and reliable automatic dense image matching process is needed.
Nowadays there is a large number of matching algorithm developed by the scientific community, capable of determining the exact 3D coordinate of the object for each pixel.In this field, algorithms like PMVS (Path-based Multi-View Stereo), SGM (Semi-Global Matching), MicMac (Multi-Images Coreespondances, Méthodes Automatiques de Corrélation) or SURE (Surface Reconstruction from imagery) have acquired important roles, and the use of one or other algorithm depends, mainly, on the shooting configuration.
For the present case of study a convergent protocol in the image acquisition has been followed and a MicMac strategy was used.
As a result a dense point cloud model was obtained (Figure 1).

2D Digital image correlation
Although one of the basic requirements to understand the structural behavior of the construction implies the need of an accurate geometric model, this is not the only one.
Complementary to a suitable material model is needed in order to represent the mechanical behavior of the materials.
Nowadays, several approaches based on interferometric and white-light methods, have been used successfully in experimental solid mechanics, to evaluate the material behavior at different stress values.These methodologies compared with the traditional measurement systems, like strain gauges or LVDTs (Linear Variable Differential Transformer), provide some important advantages as a full-field data and contact-free.Among these techniques, non-interferometric methods based on image processing, such as digital image correlation (DIC), have been increasingly used (Ghorbani et al., 2014;Pan et al., 2009b).

Pre-processing stage
In contrast to the image-based modelling strategy, DIC requires the preparation of the analyzed specimen.This preparation involves several steps: (i) a suitable test surface pattern; (ii) geometrical camera calibration; and (iii) pose estimation.
The basic principle of DIC is the tracking (or matching) of the same image subsets located in two digital images of the test specimen surface recorded before and after deformation (Figure 2).As an initial approximation a correlation coefficient is used, later, a non-linear strategy is considered, in order to analysis the deformation suffered by the subset in the reference and deformed image.Considering multiple subsets in the images, the analysis of them, can provided a full-field displacement whose resolution depends of the subset size.It should be noted, that for this purpose, for the case where the natural texture is insufficient, the specimen surface must be covered with a random speckle pattern.This type of patterns provides a surface with random grey levels, whose geometric distribution (size) plays a critical role in the achieved accuracy.
In order to minimize this uncertainty, different paint strategies have been tested and evaluated.For these evaluations a global parameter, called MIG (Mean Intensity Gradient) (Pan et al., 2010a) (Equations 1-2), was used.
Where: f x (x ij ), f y (x ij ) are the gradient along x and y direction respectively W,H are image width and height in pixels It is remarkable that both mean bias error and standard deviation error of the measured displacements are closely related with the MIG coefficient, as show (Pan et al., 2010a).Large MIG values involves large grey variations and therefore a wide number of singular features.Hence, a good speckle patter will have a large MIG value (Table 1) (Figure 3).Since the captured images are the 2D projection of the specimen surface, the accuracy in the estimation of the displacements is conditioned not only by a proper preparation of the specimen, but also by an accurate geometrical distortion correction and perpendicularity between the sensor and the specimen.
Camera calibration is technically not a part of the DIC technique, but in practice and given the image-based DIC nature, some errors can be introduced due to the geometrical distortions of the image captured.For this reason, a geometrical distortion correction approach, proposed by (Jean-Yves, 2010), is employed.
As discussed above, one of the basic requirements for the successful use of the 2D DIC method involves the use of a flat specimen surface and a parallelism between the specimen and the camera sensor, in order to keep the pixel size as constant as possible.Theoretical analysis done by (Meng et al., 2006) indicates that non-parallel angles less than 5º introduces a displacement error smaller than 0.01 pixel.In practice, this uncertainty is unavoidable due to factors such as (Pan et al., 2009a): (i) deviation of the specimen surface from ideal plane, (ii) non-perfect contact between the loading device and the test specimen; and (iii) the Poisson´s ratio of the material tested.In order to minimize this, whenever the use non-telecentric (i.e standard) lens was used, the camera sensor was placed far from the specimen to approximate this to a telecentric imaging system.Also a robust camera pose estimation approach is considered (Schweighofer and Pinz, 2006), in order to identify and correct the pose of the sensor in relation to the specimen.Since this approach need, at least, four coplanar but noncollinear points to determine the camera pose, a chessboard A B C D panel was used.Then, the corners are extracted, using for this the enhanced version of the traditional Harris corner detector proposed by (Jean-Yves, 2010).

Processing stage
The correlation process (processing stage) starts with the ROI (Region Of Interest) discretization into clusters called subsets.
Once the digital image has been split into subsets, a correlation criterion is used, in order to track the different image subsets between the initial state (reference image) and the currently analyzed image (deformed state).Nowadays, several correlation criterions has been designed and used in literature, among which have acquired special attention, given their robustness (Pan et al., 2009a) Where: i,j are the pixel (or subpixel) numbers x,y are the coordinates of the pixel in the reference image x',y' are the coordinates of the pixel in the deformed image f,g are the reference and deformed image functions M is the number of pixels belonging to the subset Δg,Δf are the grey value difference between the point analysed and the subset mean value Once the different subsets are localized, a non-linear optimization algorithm is required.In contrast to the classical non-linear optimization, based on the forward additive Gauss-Newton method, the Inverse Compositional Gauss-Newton (Pan, 2014) was used.This approach gives equivalent results, in terms of accuracy, but with lower computational costs.
Considering the DIC as an image-based approach, the accuracy of this methodology partially depends of the pixel size.Also, if subpixel interpolation schemes are not employed, a displacement accuracy of up to one pixel can be achieved (Gencturk et al., 2014).In order to improve this accuracy, various interpolation schemes have been proposed (Luu et al., 2011): (i) bicubic; (ii)bicubic spline; (iii) bicubic B-spline; (iv) biquintic B-spline.For the present case a biquintic B-spline had been used, since this interpolation scheme mitigates DIC errors form interpolating signals with high frequency content (Luu et al., 2011).
As described before, the implementation of the DIC analysis includes the use of a ROI on the specimen analyzed.The pixels inside this ROI are considered as points to be computed by the abovementioned methodology.Conventional correlation calculation generally starts with the upper-left point of the ROI, after that the calculation analysis is carried out point by point along each row or column.However, this approach presents important drawbacks as (Pan et al., 2009b): (i) presence of geometrical discontinuities (cracks or holes) or irregular shapes; (ii) discontinuous deformations; and (iii) presence of wrong data points that can provide a wrong initial guess of the next point evaluated.
Considering the stated above, a universally applicable reliability-guided (RG-DIC) method, developed by (Pan et al., 2010b), was used.In this approach, the ZNCC coefficient is used to identify the reliability of the point computed.This means, that the correlation stage begins with an initial seed and then, is guided through the neighbors with highest ZNCC coefficient.This implies, that for the full-field displacement (displacement suffered for all subset evaluated), the calculation path was made along the most reliable direction and possible error propagation of the conventional DIC method can be entirely avoided (Pan et al., 2010b).

Post-processing stage
This full field displacement, obtained on the ROI, is the basis for the subsequent mechanical properties evaluation.For this mechanical evaluation an extensometer is needed.This extensometer gives the displacement values between two subset and consequently the strain values in both axis (x and y).The relation with the stress applied affords the stress-strain curves.

Experimental setup
The main goal of the DIC experimental campaign presented is to obtain the mechanical properties (through the stress-strain curves) of the different material that compose the dome.From these curves the mechanical properties of the material can be extracted ( For the case of study evaluated, and due to the lack of materials to be tested, only compression tests were performed on brick, extracted from a part of the dome and mortar specimens, manufactured with similar mechanical properties, as those expected in the dome, considering the guidelines exposed in the regulation UNE-EN 1052-11 (AENOR, 2007) and UNE-EN 772-1 (AENOR, 2011).The image capture system has been composed by a Canon EOS 500D digital camera equipped with a 70-300mm macro lens (fixed in 200mm) and supported by a tripod (Figure 4) (Figure 5) (Table 3).During the processing stage in the DIC workflow a subset size of 40 pixels was considered and a value of 30 pixels for the subset-step (distance between subsets).Once the deformations have been evaluated through DIC, an extensometer is placed between two points in order to obtain the stress-strain curve (Figure 6). Figure 6.Stress-stain curve, for brick, obtained by DIC approach.

Point cloud pre-processing
The image-based modelling strategy shows great potentialities, in terms of accuracy and dense reconstruction, thanks to the robustness of the presented methodologies.However, the generated point cloud may present some noise.In order to improve the quality of this product, a noise reduction based on a approach is presented: (i) normal estimation proposed by (Grimm and Smart, 2011); and (ii) adaptive Moving Least Square proposed by (Dey and Sun, 2005).Finally a mesh reconstruction is used.As a result a suitable mesh for CAD modelling is generated.

CAD modelling
For a suitable structural evaluation and thus a comprehensive knowledge about the damage that can suffer the structure, it is imperative to build an accurate CAD model.Thanks to the presence of superabundant and accurate information in the point cloud provided by the image-based modelling this requirement is possible.
Traditionally this step could be made through three different approaches (Sánchez-Aparicio et al., 2014): (i) orthogonal views; (ii) sections applied along directions and over the mesh; and (iii) Non-Uniform Rational B-Splines (NURBS) generated from the mesh.
For the obtained point cloud, all the strategies defined above, can be applied.Due to the geometrical properties of the structure, which has large deformations, and due to the need of generation a CAD model of the initial state, the second approach out of those mentioned above has been chosen, and complemented with additional modelling strategies defined by (Wang et al., 2012).

Numerical simulations
Nowadays the numerical simulations through finite element method (FEM) have shown high efficiency in the simulation of historical constructions (Villarino et al., 2014a).This numerical strategy requires different type of input data.Some of this input data can be provided by the presented methodology such as: (i) CAD model (image-based modelling); and (ii) the mechanical properties of the material (digital image correlation).
Through the image-based modelling approach defined above, the deformed shape obtained in combination with the use of orthoimages obtained through the orthographic projection of the dense point cloud can be used to evaluate the agents that caused the structural damage.
For the present case of study, the evaluation of both products (deformed shape and pathological inspection through the orthophoto), clearly show a possible cause of the structural damage (Figure 7) (Figure 8).The deformation of Aliste´s dome can be attributed to the presence of an asymmetric load.This eccentric load, attributed to the weight of the rooftop, only affects an area equivalent to one of the roof´s divisions.Further research sheds light on the causes, supported by the foregoing observation: a leak on the roof, which causes the deterioration of the timber structure and its deformation until part of its load rests on the dome.
Complementary to what is exposed above, a material model is needed to adequately represent the masonry behavior as well as a modeling strategy.For the material model of the masonry structure, a non-linear constitutive model was considered, based on Total Strain Crack Model (TSCM) (Diana, 2011).This constitutive model shows a good compromise between stability, in the crack opening, and a moderate computational cost.
Considering the mechanical properties of the different materials obtained through DIC (Section 2.2.4), homogenization formulas (Tensing, 2013) are used, in order to obtain the masonry mechanical properties of the dome (Table 4).
In the case of the timber structure´s mechanical properties, the content of Concerning the modelling strategies, several approaches can be use (Lourenço, 2009): (i) detailed micro-modelling; (ii) simplified micro-modelling; and (iii) macro-modelling.For the present case study, a macro-modelling approach is used, in which the tile units and mortars joints are smeared together as one continuum.As a result of the exposed above, a 3D FEM has been performed, considering both structures, the masonry, that constitute the dome and the timber roof structure.The discretization of the model has been carried out considering a 10-node isoparametric tetrahedral element, contact-target elements and pinned supports in their bases.All these result in a total of 58,682 discrete elements (26,903 for the dome and 31,779 for the timber structure) calculated in the numerical model (Figure 9). Figure 9. FEM adopted mesh for the dome and timber structure.
One the model has been performed; several hypotheses have been taken into account: (i) the structure´s own-weight; (ii) snow load acting on the roof (considering this as the worst-case estimate); and (iii) a damage model due to the presence of an eccentric load (whose origin is attributed to a timber structure degradation process).
As conclusion derived from the numerical simulations, in a normal case (own-weight/snow-case), the structure seems to be stable and safe.In any case, the maximum tensile strength, with a value of 0.2N/mm 2 , is exceeded.Also the maximum displacements predicted were of 3mm which is insufficient to contact the dome (Figure 10).As mentioned previously, the origin of the damage that the dome structure suffers is caused by a part of the timber structure that transfer its load onto.For this reason a third model was evaluated, taking into account an eccentric load acting on the dome (which simulate the timber structure degradation process) (Figure 11).

Model-distance comparison approach: a complementary tool to validate the numerical model
At this point three models have been obtained, namely: (i) numerical model (initial and deformed state); and (ii) imagebased model (considered as "ground truth").The distance between these 3D models, provides an overview of the numerical simulation´s robustness.This robustness evaluation can be used also to assess additional aspects of the initial model considerations and also evaluate different hypothesis of the damage (Figure 12).For this purpose, a model comparison has been evaluated through the strategy builds on the Hausdorff distance (Equation 8).
  max sup inf ( , ),sup inf ( , ) where d H distance between models sup sub-set supremum inf sub-set infimum X components (nodes) of the numerical model Y components (points) of the image-based model With the considerations exposed in the present paper, the damage has been properly located and the load hypothesis seems to be correct and also, in global terms, the movement of the dome.However the dome displacements, are inaccurate (high rigidity) with discrepancies of several centimetres with respect to the image-based model.Such error can be associated with the presence of elastic supports.Further investigations and sensibility analysis are required in order to improve this part and obtain a better simulation of the dome.

CONCLUSIONS
The proposed method shows strong potentialities in the field of built cultural heritage preservation.Through this approach it is possible to solve important drawbacks when the preserving cultural heritage structures.Thanks to the SfM approach it is possible to have a complete tridimensional model of the deformed state, which can be used to estimate the initial state and the global deformation of the structure.In contrast with the traditional damage visual inspections, the use of orthoimages or the dense point cloud can provide the tools needed to evaluate the damages that the structure suffers.Complementary to the exposed above, the use of the DIC technique can provide fullfield information of the specimen tested, in contrast to traditional tools like LVDTs or strain gauges.Said information, not only, can be used to estimate the engineering curves (stressstrain) but can be used to understand the behaviour of different construction solutions which are characteristic in the cultural heritage field.
Among the different advantages it should be remarked the following: (i) non-invasive technology; (ii) flexibility; (iii) low associated cost; and (iv) great efficiency, affording the radiometric, geometric and structural characterization of the analysed object with high accuracy and detail even in complex architectural scenarios.
The dual combination of photogrammetric and finite element approaches is a reliable foundation from which to analyse appropriate restoration actions, following the procedure show above.Thanks to the versatility of the methodologies presented and their complementarity, it is possible to estimate the origin of the damage.On one hand, through the analysis of the different photogrammetric products it is possible to establish different hypothesis (which in many cases are unknown) of the initial state or loads acting on the structure and its verification.On the other hand, a comparison between the numerical model and the image-based one, can provided relevant information about its degradation process.
For the case of study that concerns this article (the Aliste´s dome), the causes of the damage seem to be clear.Nevertheless, further investigations will be required to approximate the obtained numerical data to the real one.Part of these investigations will be focused on a better understanding of the structural behavior of the dome with elastic supports and also several sensitivity analysis, considering as the variables to analysis the mechanical properties of the material and also elastic supports.The possibility to consider the comparison strategy (once the numerical model can predict better the deformations), defined in Section 3.4 into a deterministic framework for FEM model calibration will also be evaluated.

Figure 1 .
Figure 1.Dense point cloud obtained from the image-based modelling strategy.

Figure 2 .
Figure 2. Digital image correlation general outline.In red the reference subset, in blue the initial subset on the deformed image and in yellow the final subset obtained.

Figure 3 .
Figure 3. Speckle patterns tested for the DIC test.

Figure 4 .
Figure 4. Experimental setup during the DIC test.

Figure 7 .
Figure 7. Shape comparison between the initial (in green) and the deformed stage (in blue).In red the load acting on the structure, in black the structural movement and in yellow the hinges observed.

Figure 8 .
Figure 8.An orthoimage (plan view) of the dome: result of the inspection in search for damage (red: deformed zone, blue: cracks and green: material removal).

Figure 10 .
Figure 10.Maximum displacements (snow case), expressed in mm, of the timber structure.The most deformed area is displayed in green-yellow colours.

Figure 11 .
Figure 11.Results of the numerical simulation, in terms of damage (in red colour), considering the degradation process of the timber structure and in yellow the hinges predicted.

Figure 12 .
Figure 12.Comparison between the numerical model (initial state) and the image-based model.In red the maximum discrepancies and in green the minimum.

Table 1 .
Comparison, in term of MIG, between the different speckle patterns tested.

Table 5 .
Table 5 are considered as show (de Oliveira Feio, 2005) into an orthotropic material model.Mechanical properties for the timber structure.