Out of Plumb Assessment for Cylindrical-Like Minaret Structures Using Geometric Primitives Fitting

Abstract

Cultural heritage documentation and monitoring represents one of the major tasks for experts in the field of surveying, photogrammetry and geospatial engineering. Cultural heritage objects in countries like Iraq and Syria have suffered from intentional destruction or demolition during the last few years. Furthermore, many heritage sites in the mentioned places have an added religious value, and were either destroyed or are still in danger. Mosques, churches and shrines typically include one or multiple tower structures, and these towers or minarets are in many cases cylindrical-like objects. Because of their tall and relatively thin body, and adding in their age of construction, observing their inclination or out of plumb is of high importance. Accordingly, it is highly necessary for the continuous monitoring and assessment of their preservation and restoration. In this paper, we suggest an out of plumb assessment procedure using a geometric primitives least squares fitting technique, namely, cylinders, cones, and 3D circles. The approach is based on reconstructing a dense point cloud of the minaret tower which is scaled to reality by control points. Accordingly, the out of plumb is computed by fitting one of the mentioned 3D primitives to the minaret point cloud where its major axis orientation is computed. Two experimental tests of heritage objects in Iraq are presented: the lost heritage of the minaret al Hadbaa in the city of Mosul (1173 AD) and an existing inclined minaret of the religious shrine of Imam Musa AlKadhim in Baghdad (1058 AD). The results show the efficiency of the suggested methodology where the out of plumb is computed as 0.45m±1cm for the shrine minaret and 1.90m±10cm for the model of the minaret al Hadbaa.

1. Introduction

The preservation of cultural heritage objects requires regular monitoring and rehabilitation work. Normally, experts in the field of surveying, photogrammetry and geospatial engineering take the lead in such monitoring because of the high precision and reliability needed to conduct the measurements. Valuable cultural heritage sites in countries like Iraq and Syria have suffered from intentional destruction or demolition during the last few years. Mosques, churches and shrines are at the top of the list of heritage objects that need attention because of their added religious and spiritual value to the people living there [1]. These religious cultural heritage objects typically include one or more tower structures, and these towers or minarets are in many cases cylindrical-like objects. Minarets are vertical elongated structures with a relatively small base diameter. This, in addition to their construction material types, non-robust foundations, and their age of construction, makes observing their inclination or out of plumb of high importance for their preservation and maintenance.

Technically, out of plumb measuring and monitoring techniques can be categorized as:

-

Geodetic techniques: these techniques are the most accurate, where several control targets are placed on specific preselected locations of the study object. Total station instruments are used to fix the coordinates of the control points and baselines with millimetric accuracy in a predesigned network [2]. Frequently, GPS is also used for absolute positioning of the control network. This technique is based on measuring sparse target control points on the structure surface body (Figure 1a). The cylindrical axis is then defined at multiple heights from the best fit circles of the points. Other old techniques included using theodolite instruments, whereby angles are halved to define the cylindrical axis (Figure 1b) at multiple heights.

In recent years, geodetic techniques for tall structures like towers and skyscrapers have started using a combined sensor-assisted approach [3]. As an example, Leica Geosystems has developed the Core Wall Control Survey System (CWCS), which uses GPS sensors combined with high-precision inclination sensors and total stations to deliver precise and reliable coordinates [4].

-

Laser scanning techniques TLS: these techniques are applied normally from terrestrial-based time of flight (TOF) ranging scanners or phase-based laser scanners (Figure 1d). The TLS approach [5,6,7,8,9,10] has the advantage of high positioning accuracy, while it depends on the intrinsic noise of the instrument and how good the coregistration of the multiple scanning point clouds is. Compared to image-based approaches, the TLS is higher in cost, heavier, and requires a level of proficiency.

-

Photogrammetric techniques: these techniques are applied using image-based methods (Figure 1c). These techniques [11,12] range from centimetric to sub-centimetric accuracy depending on the type of camera, image network preplanning, accessibility around the object, illumination, etc.

In this paper, we use photogrammetric image-based methods to create the minaret point clouds, because they have the following advantages:

• They have lower cost compared to the geodetic and laser scanning approaches.
• They are portable and can be applied in situ without disturbing site visitors.
• The available state-of-the-art software tools for automated image-based 3D modeling [13,14,15].
• The possibility to reconstruct lost heritage models using archived and crowdsourced images and videos.

Accordingly, this ease of use, highly detailed reconstruction of the objects, and the permitted accuracy mentioned above comprise the advantages of using this photogrammetric technique for the out of plumb assessment of heritage minarets. In Section 2, below, we introduce three possible techniques of the geometric primitive fitting for the out of plumb calculations. Then, two experimental tests will be discussed in Section 3 for the out of plumb assessment of cultural heritage minaret structures.

2. Out of Plumb Computation Method

To calculate the out of plumb inclination of a cylindrical-like minaret body, we suggest a geometric primitive fitting, namely a cylinder, cone or a 3D circle fitting.

2.1. Cylinder Fitting

There are different methods of least squares fitting of a cylindrical shape [16,17,18]; either linear- or nonlinear-based solutions. Normally, seven parameters are used to define a cylinder (Figure 2):

• Point of origin lying on the cylinder axis, defined as xo, yo, zo.
• Cosine direction normal a, b, c of the cylinder axis.
• Radius of the cylinder R.

In this subsection, we will introduce a modified approach in computing the unknown parameters xo, yo, zo, R and normal cosine direction $\left[\mathrm{a},\mathrm{b},\mathrm{c}\right]$ of the cylinder axis using nonlinear least squares adjustment as follows:

-

The first derivation step is to apply the cross product between the vector derived from the coordinate differences and the cosine direction as:

$$\mathrm{E}=\left[\begin{array}{c}\mathrm{a}\\ \mathrm{b}\\ \mathrm{c}\end{array}\right]\times \left[\begin{array}{c}\mathrm{x}-\mathrm{xo}\\ \mathrm{y}-\mathrm{yo}\\ \mathrm{z}-\mathrm{zo}\end{array}\right]\to \mathrm{d}=\parallel \mathrm{E}\parallel$$

(1)

where × refers to cross product and ∥ ∥ refers to norm.

Therefore, the objective function is to minimize the difference between the distances from the fitting points to the cylinder axis and its radius R as illustrated in Figure 2.

Or

$$\mathrm{EF}=\min .\sum {\left(\mathrm{di}-\mathrm{R}\right)}^2$$

(2)

The observation equation $\mathrm{V}=\mathrm{B}\Delta -\mathrm{F}$ can be formed for $\mathrm{n}$ points as follows:

$$\underset{{\mathrm{V}}_{\mathrm{n}\ast 1}}{\underbrace{\left[\begin{array}{c}{\mathrm{v}}_1\\ {\mathrm{v}}_2\\ ⋮\\ {\mathrm{v}}_{\mathrm{n}}\end{array}\right]=}}\underset{{\mathrm{B}}_{\mathrm{n}\ast 7}}{\underbrace{\left[\begin{array}{ccccccc}\frac{\partial {\mathrm{F}}_1}{\partial \mathrm{a}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{b}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{c}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{R}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{xo}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{yo}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{zo}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{a}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{b}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{c}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{R}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{xo}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{yo}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{zo}}\end{array}\right]}}\underset{{\Delta }_{7\ast 1}}{\underbrace{\left[\begin{array}{c}\Delta \mathrm{a}\\ \Delta \mathrm{b}\\ \Delta \mathrm{c}\\ \Delta \mathrm{R}\\ \Delta \mathrm{xo}\\ \Delta \mathrm{yo}\\ \Delta \mathrm{zo}\end{array}\right]}}-\underset{{\mathrm{F}}_{\mathrm{n}\ast 1}}{\underbrace{\left[\begin{array}{c}{\mathrm{F}}_1\\ {\mathrm{F}}_2\\ ⋮\\ {\mathrm{F}}_{\mathrm{n}}\end{array}\right]}}$$

(3)

where

• Δ: vector of corrections to fitting parameters.
• V: vector of residual errors.

Furthermore, to avoid the rank deficiency in the partial derivatives B matrix, two constraint equations C₁ and C₂ are added to the system as follows:

1st constraint: the cosine direction vector a, b, and c should be normalized as $({\mathrm{a}}^2+{\mathrm{b}}^2+{\mathrm{c}}^2=1)$, Therefore, the constraint equation matrix C₁ will be:

$${\mathrm{C}}_1=\left[\begin{array}{ccccccc}2\mathrm{a} & 2\mathrm{b} & 2\mathrm{c} & 0 & 0 & 0 & 0\end{array}\left], {\mathrm{g}}_1=\right[1-{\mathrm{a}}^2-{\mathrm{b}}^2-{\mathrm{c}}^2\right]$$

(4)

2nd constraint: the second constraint equation that should be added is based on the dot product between the cosine direction a, b, and c of the cylinder axis and the origin center point 0, or:

$$\mathrm{a}.\mathrm{xo}+\mathrm{b}.\mathrm{yo}+\mathrm{c}.\mathrm{zo}=0$$

(5)

Therefore:

$${\mathrm{C}}_2=\left[\begin{array}{ccccccc}\mathrm{xo}& \mathrm{yo}& \mathrm{zo}& \mathrm{a}& \mathrm{b}& \mathrm{c}& 0\end{array}\right], {\mathrm{g}}_2=-\left[\mathrm{a}.\mathrm{xo}+\mathrm{b}.\mathrm{yo}+\mathrm{c}.\mathrm{zo}\right]$$

(6)

Accordingly, Helemert’s method is used to solve the constrained problem as follows [19]:

$$\left[\begin{array}{cc}\mathrm{N}& {\mathrm{C}}^{\mathrm{t}}\\ \mathrm{C}& 0\end{array}\right]\left[\begin{array}{c}\Delta \\ {\mathrm{k}}_{\mathrm{c}}\end{array}\right]=\left[\begin{array}{c}\mathrm{t}\\ \mathrm{g}\end{array}\right] \mathrm{where} \mathrm{C}=\left[\begin{array}{c}{\mathrm{C}}_1\\ {\mathrm{C}}_2\end{array}\right], \mathrm{g}=\left[\begin{array}{c}{\mathrm{g}}_1\\ {\mathrm{g}}_2\end{array}\right], \mathrm{N}={\mathrm{B}}^{\mathrm{t}}\mathrm{B}, \mathrm{and} \mathrm{t}={\mathrm{B}}^{\mathrm{t}}\mathrm{F}$$

(7)

Because of the large number of fitting points, it is expected to have some amount of noisy points, and as a result, a possibly error-prone estimation of the parameters. Therefore, the detection of blundered points using a voting method is applied in the form of the well-known “RANdom SAmple Consensus” (RANSAC) [20,21]. In the RANSAC algorithm, we assume the observational data to include both: inliers and outliers. Inliers can be described as observations without blunders with respect to a specific model, while outliers do not fit that model in any condition. It should be noted that the RANSAC algorithm can find adjustment models even if 50% of the input data set is affected by noise points. On the other hand, there is no guarantee of reaching the optimal solution if we limit the number of iterations.

2.2. Cone Fitting

A similar approach of least squares fitting can be applied to the cone shape [16]. Normally, eight parameters are used to define a cone (Figure 3):

-

Point of origin lying on the cone axis and defined as xo, yo, zo.

-

Cosine direction numbers a, b, c of the cone axis (pointing to the apex).

-

Orthogonal distance from the point of origin to cone s.

-

The cone’s apex semi-angle θ.

The least squares adjustment technique will be applied to minimize the computed distance d between the cone axis and the fitting points as:

$$\mathrm{d}={\mathrm{f}}_{\mathrm{i}}\cos \mathsf{\theta }+{\mathrm{g}}_{\mathrm{i}}\sin \mathsf{\theta }-\mathrm{s}$$

(8)

where

$${\mathrm{f}}_{\mathrm{i}}=\sqrt{{\mathrm{u}}^2+{\mathrm{v}}^2+{\mathrm{w}}^2}$$

(9)

$$\begin{array}{l}\mathrm{u}=\mathrm{c}\left({\mathrm{y}}_{\mathrm{i}}-\mathrm{y}\right)-\mathrm{b}\left({\mathrm{z}}_{\mathrm{i}}-\mathrm{z}\right)\\ \mathrm{v}=\mathrm{a}\left({\mathrm{z}}_{\mathrm{i}}-\mathrm{z}\right)-\mathrm{c}\left({\mathrm{x}}_{\mathrm{i}}-\mathrm{x}\right)\\ \mathrm{w}=\mathrm{b}\left({\mathrm{x}}_{\mathrm{i}}-\mathrm{x}\right)-\mathrm{a}\left({\mathrm{y}}_{\mathrm{i}}-\mathrm{y}\right)\end{array}$$

(10)

$${\mathrm{g}}_{\mathrm{i}}=\mathrm{a}\left({\mathrm{x}}_{\mathrm{i}}-\mathrm{x}\right)+\mathrm{b}\left({\mathrm{y}}_{\mathrm{i}}-\mathrm{y}\right)+\mathrm{c}\left({\mathrm{z}}_{\mathrm{i}}-\mathrm{z}\right)$$

(11)

Therefore, the objective function is to minimize the distances from the fitting points to the cone surface, as illustrated in Figure 3.

Or

$$\mathrm{F}=\min .\sum {\left(\mathrm{di}\right)}^2$$

(12)

The observation equation $\mathrm{V}=\mathrm{B}\Delta -\mathrm{F}$ can be formed for n points as follows:

$$\underset{{\mathrm{V}}_{\mathrm{n}\ast 1}}{\underbrace{\left[\begin{array}{c}{\mathrm{v}}_1\\ {\mathrm{v}}_2\\ ⋮\\ {\mathrm{v}}_{\mathrm{n}}\end{array}\right]=}}\underset{{\mathrm{B}}_{\mathrm{n}\ast 8}}{\underbrace{\left[\begin{array}{cccccccc}\frac{\partial {\mathrm{F}}_1}{\partial \mathrm{a}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{b}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{c}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{s}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{xo}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{yo}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{zo}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathsf{\theta }}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{a}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{b}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{c}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{s}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{xo}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{yo}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{zo}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathsf{\theta }}\end{array}\right]}}\underset{{\Delta }_{8\ast 1}}{\underbrace{\left[\begin{array}{c}\Delta \mathrm{a}\\ \Delta \mathrm{b}\\ \Delta \mathrm{c}\\ \Delta \mathrm{s}\\ \Delta \mathrm{xo}\\ \Delta \mathrm{yo}\\ \Delta \mathrm{zo}\\ \Delta \mathsf{\theta }\end{array}\right]}}-\underset{{\mathrm{F}}_{\mathrm{n}\ast 1}}{\underbrace{\left[\begin{array}{c}{\mathrm{F}}_1\\ {\mathrm{F}}_2\\ ⋮\\ {\mathrm{F}}_{\mathrm{n}}\end{array}\right]}}$$

(13)

Similar geometric constraints mentioned in the cylinder fitting problem are used to run a full rank least squares adjustment (Equation (7)) of the best cone fitting.

After the calculation of the fitting parameters either for the cylinder primitive or the cone, the angular orientation of the cylinder axis is computed as follows:

$$\begin{array}{l}\mathsf{\alpha }={\cos }^{-1}(\mathrm{a}/\sqrt{{\mathrm{a}}^2+{\mathrm{b}}^2+{\mathrm{c}}^2}\\ \mathsf{\beta }={\cos }^{-1}(\mathrm{b}/\sqrt{{\mathrm{a}}^2+{\mathrm{b}}^2+{\mathrm{c}}^2}\\ \mathsf{\gamma }={\cos }^{-1}(\mathrm{c}/\sqrt{{\mathrm{a}}^2+{\mathrm{b}}^2+{\mathrm{c}}^2}\end{array}$$

(14)

However, for such vertical primitive objects, we only need the γ inclination angle to compute the out of plumb as follows (Figure 2b):

$$\mathrm{Out} \mathrm{of} \mathrm{plumb} = \mathsf{\gamma }\times \mathrm{height}$$

(15)

Whereas $\mathsf{\alpha }$ and $\mathsf{\beta }$ are used to determine the out of plumb direction if needed. The standard deviation of the out of plumb in Equation (15) can be computed using the propagation of errors technique. However, tests showed that the computed standard deviations of the cosine directions are of small magnitude, and the out of plumb accuracy is mainly related to the estimated height accuracy. Accordingly, the estimated accuracy of the out of plumb in the following experimental tests are assumed equivalent to the accuracy of the primitives fitting of the point clouds.

2.3. Circle Fitting in 3D space

Fitting a circle in 3D space (Figure 4) needs a multistep solution approach, as follows [16,22]:

• Compute the best fit least squares plane of the data as explained previously.
• Rotate the points such that the least squares plane is the XY plane.
• Rotate data points onto the XY plane.
• Compute the 2D circle fit in the XY plane.
• Rotate back to the original orientation in 3D space.

2.3.1. Fitting a 2D Circle

The least squares adjustment of a circle best fit in a 2D space is shown in Figure 5 as follows:

• Given: XY points, i = 1:n
• Required: the adjusted radius r, and circle center ${\mathrm{x}}_{\mathrm{c}},{\mathrm{y}}_{\mathrm{c}}$

Circle equation can be formulated as:

$${\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{c}}\right)}^2+{\left(\mathrm{y}-{\mathrm{y}}_{\mathrm{c}}\right)}^2-{\mathrm{r}}^2=0$$

(16)

The observation equation will be as follows in the form $\mathrm{V}=\mathrm{B}\Delta -\mathrm{F}$

$$\left[\begin{array}{c}{\mathrm{v}}_{\mathrm{x}1}\\ {\mathrm{v}}_{\mathrm{y}1}\\ {\mathrm{v}}_{\mathrm{x}2}\\ ⋮\\ {\mathrm{v}}_{\mathrm{yn}}\\ {\mathrm{v}}_{\mathrm{zn}}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial {\mathrm{F}}_1}{\partial {\mathrm{x}}_{\mathrm{c}}}& \frac{\partial {\mathrm{F}}_1}{\partial {\mathrm{y}}_{\mathrm{c}}}& \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{r}}\\ ⋮& ⋮& ⋮\\ \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial {\mathrm{x}}_{\mathrm{c}}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial {\mathrm{y}}_{\mathrm{c}}}& \frac{\partial {\mathrm{F}}_{\mathrm{n}}}{\partial \mathrm{r}}\end{array}\right]\left[\begin{array}{c}{\mathsf{\delta }\mathrm{x}}_{\mathrm{c}}\\ {\mathsf{\delta }\mathrm{y}}_{\mathrm{c}}\\ \mathsf{\delta }\mathrm{r}\end{array}\right]-\left[\begin{array}{l}{\mathrm{r}}^2-{\left({\mathrm{x}}_1-{\mathrm{x}}_{\mathrm{c}}^{\mathrm{o}}\right)}^2-{\left({\mathrm{x}}_1-{\mathrm{y}}_{\mathrm{c}}^{\mathrm{o}}\right)}^2\\ ⋮\hspace{1em}\hspace{1em}\hspace{1em}⋮\hspace{1em}\hspace{1em}\hspace{1em}⋮\hspace{1em}\hspace{1em}\hspace{1em}⋮\\ {\mathrm{r}}^2-{\left({\mathrm{x}}_{\mathrm{n}}-{\mathrm{x}}_{\mathrm{c}}^{\mathrm{o}}\right)}^2-{\left({\mathrm{x}}_{\mathrm{n}}-{\mathrm{y}}_{\mathrm{c}}^{\mathrm{o}}\right)}^2\end{array}\right]$$

(17)

The partial derivatives of the function to the unknowns will be as follows:

$$\begin{array}{l}\frac{\partial {\mathrm{F}}_1}{\partial {\mathrm{x}}_{\mathrm{c}}}=-2\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{c}}^{\mathrm{o}}\right)\\ \frac{\partial {\mathrm{F}}_1}{\partial {\mathrm{y}}_{\mathrm{c}}}=-2\left(\mathrm{y}-{\mathrm{y}}_{\mathrm{c}}^{\mathrm{o}}\right)\\ \frac{\partial {\mathrm{F}}_1}{\partial \mathrm{R}}=-2{\mathrm{r}}^{\mathrm{o}}\end{array}$$

(18)

Then a least squares adjustment is performed as $\Delta ={({\mathrm{B}}^{\mathrm{t}}\mathrm{B})}^{-1} \left({\mathrm{B}}^{\mathrm{t}}\mathrm{F}\right)$.

2.3.2. Rodrigues Rotation Formula

To apply the 2nd step of rotating the points to the best fit plane of points, we can utilize the Rodrigues rotation formula to rotate 3D points and get their XY coordinates in the coordinate system of the plane. This is simply a problem of rotating from normal vector $\mathrm{n}1$ to normal $\mathrm{n}2$, as shown in Figure 6.

To compute the rotation matrix, we will apply the following steps:

• Find the axis and angle of rotation using cross product and dot product respectively. The axis of rotation k is a cross product (×) between plane normal $\mathrm{n}1$ and the normal of the new XY coordinates. Thus, $\mathrm{n}2={\left(0,0,1\right)}^{\mathrm{t}}$ and $\mathrm{k}=\mathrm{n}1\times \mathrm{n}2$. Furthermore, $\mathsf{\theta } = {\cos }^{-1}(\frac{\overrightarrow{\mathrm{n}}1.\overrightarrow{\mathrm{n}}2}{\parallel \overrightarrow{\mathrm{n}}1\parallel .\parallel \overrightarrow{\mathrm{n}}2\parallel })$
• Find the rotation matrix M using exponential map:

  $$\mathrm{M}={\mathrm{I}}_{3\ast 3}+\hat{\mathrm{k}}\sin \left(\mathsf{\theta }\right)+{\hat{\mathrm{k}}}^2\left(1-\cos \mathsf{\theta }\right)$$

  (19)

  where

  -

  Skew—symmetric matrix $\hat{\mathrm{k}}= \left[\begin{array}{ccc}0& \mathrm{k}\left(3\right)& \mathrm{k}\left(2\right)\\ \mathrm{k}\left(3\right)& 0& \mathrm{k}\left(1\right)\\ \mathrm{k}\left(2\right)& \mathrm{k}\left(1\right)& 0\end{array}\right]$

  -

  $\mathrm{k}=\left(\mathrm{n}1\times \mathrm{n}2\right)$, then normalized.

The 3rd step after rotation of points is to find the best fit circle Cpoints in 2D space ($\mathrm{xc},\mathrm{yc},\mathrm{r}$). Finally, rotate back to 3D space by taking n2 as the plane normal and $\mathrm{n}1={\left(0,0,1\right)}^{\mathrm{t}}$

$$\mathrm{C}=\mathrm{Cpoints}-\left[\mathrm{xc},\mathrm{yc},0\right]$$

(20)

$$\mathrm{C} = \mathrm{C}\ast {\mathrm{M}}^{\prime } + \mathrm{Points}\_\mathrm{mean}$$

(21)

3. Results and discussion

Two experiments were conducted at Iraqi cultural heritage sites using the suggested geometric primitive fitting. The first experiment was carried out on the inclined minaret of the Imam AlKadhim shrine in Baghdad, and the second experiment was carried out on the lost al Hadbaa minaret in Mosul.

3.1. First Experiment

The Imam Musa AlKadhim shrine was first constructed in 1058 AD (450 Hijri), and included a mosque and one minaret; in 1098 AD, two further minaret structures were added. Then, during the rule of caliph Al Mustansir Billah in 1226 AD, the whole shrine structure was completed [23,24]. In 1508 AD, Shah Ismail Safavi added two more minarets to become four, as is still the case in our current time, and as shown in Figure 7 with two domes beside the four minaret structures.

Each minaret has a foundation 11 m in the ground and an elevation 28 m above the ground, which means a total of 39 m from the base to top. In 2013, an out of plumb increase was observed in the minaret structures, in addition to some cracks in the walls, both inside and outside. This indicated a problem with the site foundations and the soil due to the groundwater and the renovations previously applied at the shrine site.

Therefore, there was a need to ensure the visitors’ safety and to find suitable construction engineering solutions to resolve the problems with the foundations and the groundwater. Furthermore, monitoring the minarets’ stability, calculating the out of plumb, and fixing the wall cracks in the shrine structures was necessary. This was achieved by an extensive project which included TLS scanning of the whole shrine site in addition to geodetic monitoring of the minaret inclinations [24].

Conventionally, geodetic monitoring techniques are applied through multiple epoch measurements of the questionable minarets of the heritage site. Total station measurements are applied to selected critical points on the cylindrical body of each minaret, which comprise circular planes at multiple selected elevations. Then for each circular plane/section, the center coordinates are determined by applying the least squares adjustment for fitting a 2D circle.

The monitoring of the shrine site was carried out from May 2014 to August 2015 using the described geodetic methods, and with photogrammetric techniques, as well. The applied photogrammetric techniques were based on manually measuring the interest points of the minaret using an SLR camera at individual strips around the four minarets [24].

For both measuring techniques, 2500 monitoring object points were measured in the designed multi planes to compute their center points, defining the vertical axis of the minaret structures. It was observed that two minaret structures were stable, while the inclination was critical for the two other minarets [24].

In this paper, the inclination of one of the suspected renovated minaret structures was determined after three years of treatment to check its stability and whether there is still a critical out of plumb in its structure.

3.1.1. Close-Range Photogrammetry

First, a pre-analysis of the camera network around the minaret structure was carried out by reconnaissance and checking accessibility spaces. Unfortunately, it was only possible to capture the image from the ground, and no drone photogrammetry was permitted at the shrine site during the daytime. Therefore, it was decided to capture the images in a highly overlapped strip of images. Images were captured using a DSLR Canon EOS Mark III camera with a fixed focal length of 24 mm and an image frame size of 5760 × 3840 pixels. An average capture distance of 6 m was designated in order to guarantee a depth accuracy of 5 mm with an average ground sample distance (GSD) of 2 mm, as shown in Equation (22) [25].

$$\mathrm{Z}=\sqrt{3\mathrm{fB}\frac{{\mathsf{\sigma }}_{\mathrm{Z}}}{\mathrm{p}}}$$

(22)

where

• σz = precision at object space [mm]
• f = focal length [mm]
• B = base line (step size) [mm]
• p = pixel pitch (width) [mm]

However, higher parts of the minaret were expected to have larger depth accuracy and GSD.

Practically, a total of 450 images were captured in two rings around the minaret structure in different portrait and landscape orientations. To avoid the long processing time required to create the dense point cloud in our laptop machine, and to perform an efficient imaging network configuration for future image capture, we applied the minimal camera network approach presented in [26].

The minimal camera network filtering technique is based on the concept of having at least three images (cameras) viewing the object point simultaneously. The approach starts by labeling the derived sparse point cloud of the object as either over-covered or fair-covered. Over-covered points are the points that appear in more than three images, while fair-covered points, refer to the points that typically appear in three cameras. Therefore, the images are considered redundant and are filtered out if they exclusively contribute points that are covered by more than three cameras, as shown in Figure 8d. The camera filtering is allowed if it maintains the desired small Base/Distance ratio. The filtering is repeated until no more redundant cameras are involved in imaging over-covered points. Accordingly, the number of images around the shrine minarets is reduced from 450 to 169 images (Figure 8a,b) and oriented using the Metashape software tool [13].

The quality of the reconstructed minaret object is evaluated by computing the variance-covariance matrix of the tie points using a separate code. The quality computations are applied through the bundle adjustment of the images and the calculation of the mean of the error ellipsoid axes for each tie point. Then every point is colored with respect to the computed mean error, as shown in Figure 8c.

Table 1. Camera specifications and image orientation quality of the first experiment.

Image Orientation Quality		Adjusted Camera Parameters
GSD [mm]	1.5	Focal length [mm]	23.88
Pixel size [µm]	6.0	p.p in x [mm]	0.113261
Maximum reprojection error [pixel]	1.36	p.p in y [mm]	0.131599
Median reprojection error [pixel]	0.23	Radial lens distortion k1	−0.133579
Mean reprojection error [pixel]	0.37	Radial lens distortion k2	0.092059

summarizes the camera specifications and the imaging network orientation quality.

Mathematically, a minimum of three control points is sufficient to geo-reference the imaging network and, subsequently, the derived point cloud. However, adding more control points is necessary to reliably check the quality of the orientation and the produced point cloud. In this experiment, it was difficult to fix more than four control points because of the limitations regulated by the shrine authorities during the field work.

The root mean squared error (RMSE) of the four control points we used is computed as: σx = 1.2 mm, σy = 2.4 mm and σz = 1.6 mm, respectively.

Afterward, the highly dense 3D point cloud of ≈3 million points was reconstructed, as shown in Figure 9a. A density analysis showed that the average number of points was 15,300 ± 4150 points/m², which is a high level of density that will certainly give a reliable 3D representation of the minaret (Figure 9b).

3.1.2. Geometric Primitive Fitting of the Shrine Minaret

As mentioned, a geometric cylinder and cone fitting is applied to the minaret point cloud to estimate the out of plumb of the minaret axis. Figure 9c shows the best fitting cylinder, and in Figure 9d, the best fitting cylinder is fused with the minaret point cloud. Using the formulas in Section 2.1, the least squares solution is converged to optimal values (Figure 10).

Geometrically, it is clear that the minaret body is not a perfect cylinder, because whenever the points are higher, the diameter is slightly decreased. Therefore, a more efficient fitting is applied either by assuming a cone minaret body for the whole point cloud or by fitting multiple cylinders in elevated adjacent sections.

Originally, the out of plumb computations using Equation (15) are applied using a cylinder and a cone fitting for the whole point cloud, and the results are listed as follows in

Table 2. Fitting the whole point cloud with a cylinder and a cone primitive.

	xo [m]	yo [m]	zo [m]	a	b	c	Out of Plumb [m]
Cylinder fitting	136.94	114.89	31.37	0.019216	−0.019097	0.99963	0.43
Cone fitting	136.94	114.89	31.37	0.020305	−0.019289	0.99961	0.44

, where the computed out of plumb is shown:

To evaluate the goodness of the primitive fitting, we compute the distance between every point and both the cylinder and the cone surfaces, respectively. The results are shown where the mean distance is ±2.3 cm from the fitting cylinder surface (Figure 11a), while it is ±1.4 cm from the fitting cone surface (Figure 11b).

The third fitting option is to use a cylinder or a cone fitting at multiple elevations. Accordingly, multiple height sections with a width of 1 m are selected from the point cloud and followed by a cylindrical or conical fitting applied for each level.

Table 3. The best fitting cylinder parameters of the multi sections along the minaret structure.

Section	xo [m]	yo [m]	zo [m]	a	b	c	R [m]	Out of Plumb [m]
1st	136.84	114.99	25.72	0.0177	−0.0207	0.9996	1.608	-
2nd	136.86	114.97	26.69	0.0148	−0.0265	0.9995	1.592	0.03
3rd	136.87	114.96	27.69	0.0095	−0.0049	0.9999	1.594	0.04
4th	136.88	114.95	28.67	−0.0007	−0.0132	0.9999	1.585	0.06
5th	136.89	114.93	29.69	0.0193	−0.0176	0.9997	1.577	0.07
6th	136.91	114.91	30.68	0.0231	−0.0249	0.9994	1.574	0.10
7th	136.93	114.89	31.67	0.0251	−0.0241	0.9994	1.571	0.14
8th	136.96	114.87	32.67	0.0262	−0.0223	0.9994	1.568	0.17
9th	136.98	114.85	33.68	0.0210	−0.0169	0.9996	1.567	0.20
10th	137.00	114.83	34.68	0.0201	−0.0223	0.9995	1.563	0.23
11th	137.02	114.81	35.67	0.0205	−0.0154	0.9997	1.557	0.25
12th	137.05	114.79	36.69	0.0288	−0.0294	0.9992	1.552	0.29
13th	137.08	114.77	37.68	0.0305	−0.0241	0.9992	1.544	0.33
14th	137.11	114.74	38.67	0.0182	−0.0308	0.9994	1.534	0.37
15th	137.13	114.71	39.67	0.0163	−0.0229	0.9996	1.525	0.40
16th	137.15	114.69	40.71	0.0228	−0.0093	0.9997	1.515	0.43
17th	137.17	114.68	41.41	0.0074	−0.0139	0.9999	1.530	0.45

illustrates the seven parameters of the fitting cylinders in the multi elevation sections shown in Figure 12a. The last column in Table 2 shows the computed out of plumb inclination between every section and the minaret base using the planimetric coordinates of the involved cylinder centers as in Equation (23):

$$\mathrm{out} \mathrm{of} \mathrm{plumb}= \sqrt{{\left({\mathrm{X}}_{\mathrm{top}}-{\mathrm{X}}_{\mathrm{bottom}}\right)}^2+{\left({\mathrm{Y}}_{\mathrm{top}}-{\mathrm{Y}}_{\mathrm{bottom}}\right)}^2}$$

(23)

It is shown that the computed horizontal shift between the minaret base and its top is 45 cm, which represents the out of plumb inclination of the minaret. The result is similar to the field survey measurements implemented by professional surveyors.

As previously stated, Table 3 shows that the radius of the cylindrical minaret tower decreases from bottom to top, and fitting the whole minaret body with one cylinder or cone would not be so accurate. Similar to Figure 11, the multi cylinder fitting is evaluated by computing the distances between the point cloud and the fitting cylinders. Then the standard deviation of 1σ is computed as ±0.8 cm, indicating the high precision of the multi cylinder fitting to the point cloud, as shown in Figure 13.

A thorough verification of the computed out of plumb of Table 3 shows that the inclination is not truly linear, but somehow curvilinear. The best fitting polynomial is found to be cubic with a norm of residuals of 2 cm as shown in Figure 14:

$$\mathrm{y} = {\mathrm{p}}_1{\mathrm{x}}^3+{\mathrm{p}}_2{\mathrm{x}}^2+{\mathrm{p}}_3\mathrm{x}+{\mathrm{p}}_4$$

(24)

where

• y: the out of plumb in [m].
• x: difference of elevations [m].
• p₁ = −7.84 × 10⁻⁵, p₂ = 0.003, p₃ = 0.006, and p₄ = 0.037

Obviously, the minaret was built from bricks line by line, using the techniques available at the time of construction, which are not accurate according to our current standards; this resulted in the slightly curvilinear plumb axis of the minaret.

Finally, a 3D circle fitting is applied to the selected multiple sections of the minaret point cloud and the computed out of plumb is 45 cm with a fitting residuals standard deviation of ±6 cm. The results of the four fitting scenarios show that the best approach is to apply multiple cylindrical sections fitting.

3.2. Second Experiment

The Al-Hadbaa leaning Minaret (Figure 15a) is one of the remarkable monuments in Iraq which was destroyed in 2017 (Figure 15b) during the liberation of the old city of al-Mosul from ISIS [1,27]. The minaret was established in the year 1173 AD (568 Hijri) during the Seljuk rule of Nur ad-Din Zenki.

The minaret is a part of the construction of the Great Mosque or al-Nuri mosque and is known as Al-Hadbaa which means “the hunchback” in the Arabic language. Since 2003, the image of the minaret has decorated one of the banknotes in the Iraqi currency.

According to [28], the cylindrical body of the minaret, which has an elevation of 47.25 m, starts at a level of 17.31 m with a diameter of 2.62 m. The minaret consisted of three parts: the foundations; a prismoid base with dimensions of 12 m, 2 m below the ground; and the cylindrical-like minaret. The minaret is built out of brick with plaster materials, while the internal stairs inside the minaret are built of stone. The effect of the weather and the weak construction materials caused the leaning of the minaret body over the course of hundreds of years after it was established, with a lean of 2.19 m at a height of 47.25 m, as measured in 1981 [28].

Similar to the lost minaret of Aleppo in Syria [29], and based on the archived multimedia records of the minaret, engineering descriptions in the literature, and published 3D models in SketchUp 3D warehouse [30] and Sketchfab [31], a 3D image-based model was created (Figure 16). The model was scaled to reality based on the published dimensions of the minaret base of 9 × 9 meters [28]. Subsequently, a 3D point cloud was created, and the out of plumb computations were performed similarly to the procedure in the first experiment. It should be noted that the point cloud (≈700,000 points) was of moderate accuracy, since it was mainly based on low-resolution crowdsourced images available to the public. Furthermore, since this minaret object is lost, no field measurements could be applied, and the comparisons are based on the out of plumb reported previously.

From the geometric description of the minaret and the available images, it is obvious that the minaret body was not a completely regular cylinder or cone, since at higher points, the diameter is slightly decreased and leans with a curvilinear axis.

Therefore, we apply two kinds of fitting, similar to in the previous experiment: first, by assuming a conical body for the whole minaret; and second, by assuming multiple cone sections with width of 1 m for fitting.

The out of plumb computations of 1.9 m, as shown in

Table 4. Cone fitting parameters of the whole point cloud of al Hadbaa minaret.

	xo [m]	yo [m]	zo [m]	a	b	c	Out of Plumb [m]
Cone fitting	104.76	104.58	31.477	0.016572	−0.063281	0.99786	1.9

, result from using Equation (15) for the whole point cloud cone fitting; the least squares results are shown in Figure 17:

To evaluate the goodness of the primitive fitting, we compute the distance between every point and the cone surface. The results are shown in Figure 18a, where the median of the residuals between the cone surface and the points is −27 cm and the standard deviation is ±9.7 cm of 1σ.

Meanwhile, the out of plumb is computed as 2.3 m using the multi conical section fitting (Figure 18b), with a higher precision for a zero median and ±5 cm standard deviation of 1σ. This computed out of plumb is in agreement with the measured deviation reported in [28] using field surveying methods in 2007.

Lastly, a 3D circle fitting is applied to multiple selected sections of the minaret point cloud, and the computed out of plumb is 2.27 m with a standard deviation of the fitting residuals of ±4.8 cm.

4. Conclusions

In this paper, we presented a geometric primitive fitting approach for calculating the out of plumb of cylindrical-like cultural heritage minaret towers. Three types of primitive fittings for the minaret point clouds are suggested, namely: cylinder, cone and 3D circle.

In the first experiment, the approach is applied with respect to the minaret at the holy shrine of Imam AlKadhem by following the photogrammetric engineering techniques of different tasks. The tasks included:

• Preplanning, image capture with a professional camera.
• Automated image orientation using state-of-the-art software tool.
• Using field survey control points for reality scaling.
• Creating a dense point cloud for 3D modeling of the minaret.
• and finally, to compute the out of plumb amount by primitive fitting of the point clouds.

Error analysis showed (Figure 8c) that the oriented images were able to meet centimetric accuracy of the modeled points of the minaret. Furthermore, the RMSE of the control points was less than half a centimeter. In addition to cleaning the noise from the point cloud, which was of minor magnitude, the RANSAC technique was also applied during the fitting calculations to avoid the effect of the noisy points if they were present on the reliable fitting parameters.

The geometric primitive fitting calculations showed the existence of an out of plumb inclination of 45 cm between the top and the bottom of the shrine minaret body, as illustrated in Table 1 and Table 2. Statistical analysis showed that multi cylinder fitting to the point cloud of the minaret body was accurate to within 1 cm at 68% probability (1σ), as shown in Figure 13. Furthermore, because the minaret was built from bricks line by line, it was found that the out of plumb inclination of the minaret body is better fitted to a cubic polynomial rather than a straight line (Figure 14). These results have been reported to the shrine administration authorities in order to handle the proper structural treatment and monitoring of the minaret.

The second experiment was applied with respect to a digital model of the lost heritage site of the al Hadbaa minaret (Figure 15). A cone fitting is applied, rather than a cylinder fitting, because of the conical-like shape of the minaret. The computed out of plumb using one cone fitting for the whole point cloud was 1.9 m, with an accuracy of 10 cm (Figure 18a), whereas the computed multiple sections of cone fitting and the 3D circle fitting both resulted in 2.3 m out of plumb (Figure 18b). Using the cone primitive has a higher precision of ±5 cm, which is closer to the reported inclination results in 2007.

Accordingly, the proposed primitive fitting showed reliable estimation of the out of plumb of the cylindrical-like minarets, which can be automated and integrated with state-of-the-art software tools.

For future work, we will consider the out of plumb computation of other geometrically shaped objects, like squared minarets, using the proposed technique of primitive fitting. Similarly to our second experimental test of the al Hadbaa minaret, the lost heritage of the great mosque minaret of Aleppo in Syria is a recognizable future case study of square-like minarets.