1 Introduction

With the development of three-dimensional (3D) measurement technology, 3D digital data have become widely used in the fields of cultural heritage preservation and archaeological studies. For example, laser scanners have become more compact and provide faster results with improved accuracy. In addition, software used for point-cloud processing has becoming more sophisticated, which has facilitated the task of 3D measurement. Photogrammetry technology has also matured, and software has become easier to use and more inexpensive. As a result, even non-experts in computer science can now easily obtain 3D data of cultural properties. Accordingly, 3D digital data is now widely used for the preservation of cultural properties and for analysis and visualization (Kamakura et al., 2005; Lu et al., 2011, 2013).

Fig. 1
figure 1

First solar boat of king Khufu. The wooden boat was built about 4600 years ago and has been reassembled using parts that had been buried underground

Full size image

Virtual restoration, which reconstructs the original state of cultural properties in virtual space, is highly useful because it enables us to discover new historical facts. Many cultural properties have been damaged by the effects of aging or disasters. In many cases, physical restoration is difficult owing to the possibility of destroying the cultural properties. To address this risk, various studies have been conducted on virtual restoration, which includes stitching together partial shapes, supplementing lost parts based on the continuity of shapes, and analogizing damaged parts by learning-based approaches (Papaioannou et al., 2002; Huang et al., 2006; Vendrell-Vidal & Sánchez-Belenguer, 2014; Sizikova & Funkhouser, 2017; Hou et al., 2018).

However, virtual restoration of non-rigid cultural properties remains a challenging task. Most of the research on virtual restoration has focused on rigid objects such as stones and pottery. On the contrary, many cultural properties are composed of non-rigid materials such as wood, leather, and cloth. Because 3D data in virtual space can deform arbitrarily, uniquely determining the original state can be ambiguous. In particular, the problem cannot be resolved by 3D data alone without knowledge of the nature of the material and at times the historical background of the objects. Therefore, a method for shape restoration of non-rigid cultural properties is needed.

Fig. 2
figure 2

Deterioration state of Khufu’s boat. The timbers has been distorted by environmental changes and years of display, and many gaps between timbers are observed

Full size image

The present study focuses on two large, wooden Egyptian assets known as the ”solar boat of King Khufu” (Fig. 1) that were constructed approximately 4600 years ago and were discovered underground in 1954. The first boat, which is 42 m in length and weighs an estimated 50 tons, is currently exhibited at the Solar Boat Museum. Although physical restoration was conducted on this boat for many years, it appears to be distorted owing to the effects of long-term display and other factors. Specifically, gaps are present between components, as shown in fig. 2. Comparatively, the second boat has not been well preserved. The conservation team has conducted careful excavation and preservation works (Yoshimura & Kurokochi, 2017) and has determined that physical restoration of the second boat will be substantially more difficult than that of the first boat. Therefore, we plan to restore the first boat before restoring the second boat by analogy. Previous investigations have revealed that the structures and shapes of the two boats are similar. Thus, developing a restoration plan based on 3D data (Fig. 3) will strongly benefit the restoration of the second boat.

Fig. 3
figure 3

Point-cloud of Khufu’s first boat. We measured the boat with the IMAGER5010C laser scanner from 258 positions. The purpose of this study is to restore the outer shape of the main body of the boat. In the experiment, we extracted the outside of the main body and estimated the deformations of the timbers

Full size image
Fig. 4
figure 4

Type of deformation. The deformation of timber is limited to bending in fiber direction and twist along fiber direction. Due to the processing technique and the nature of the wood, bending in the tangential direction is assumed to be negligible

Full size image

Wood is one of the most commonly used materials with the ability to deform. Owing to its favorable adaptation in construction projects, wood has long been used to create many cultural properties. Its anisotropy in strength and deformation differs among varieties (Sadoh, 1996) and processing methods. Because a tree comprises numerous cells that grow in a cyclical manner, wood expands and contracts anisotropically depending on its moisture content in proportion to its orientation to the grain. Due to its structure, trees are typically processed as cut wood for timber which is long along the fiber direction (Dinwoodie, 1981). Moreover, timber is deformed in the fiber direction because wood is resistant to deformation in the fiber direction. Therefore, we assume that the deformation of the timber is limited in bending in the fiber direction and twisting along the fiber direction, as shown in Fig. 4. In the case of ship construction, the curved shape from top to bottom is made by laying out long bent timbers (Tanner, 2018; Nayling & Jones, 2014).

In this study, we propose a method for shape restoration based on non-rigid shape assembly with ruled surface representation for cultural properties composed of multiple timbers. Ruled surface is popular for designing developable surfaces (Pottmann et al., 2008; Bayram et al., 2012; Tang et al., 2016; Kilian et al., 2008; Rose et al., 2007) and architectures (Flöry & Pottmann, 2010; Pottmann et al., 2015). The developable surface means the surface that can be flattened without stretching or compression. The original timber was flat and was deformed to construct the artifacts. The deformation process causes the expansion of the surface. Hence, the deformed timber is not developable but still ruled if bent and twisted only along one direction. The deformation is also limited to bending and expansion to specific degrees. To solve the ambiguity problem in mutually aligning deformable components, we need to represent the deformation with a low-dimensional model based on these limitations. Because the ruled surface effectively represents the deformation of timber, we propose a deformation method referred to as Ruled Surface Free form Deformation (RS-FFD). In addition, we introduce an error metric that considers contact in the non-rigid alignment process to obtain a physically consistent solution.

We make no assumptions about the overall shape after assembly; however, we consider constraints derived from archaeological findings and physical properties of components. It is considered that the components of the boats were fixed and assembled by first drilling holes in the components and then using tenon joints or ropes to join them (Creasman, 2013). In this study, we use the positions of the holes for estimating the constrained deformation. The boats have been previously disassembled, and the parts were stored in a stone chamber while maintain the positional relationship of the components. We assume that all components to be assembled have already been excavated and that the positional correspondences of the components are accurate. It should be noted that because the first solar boat has been well preserved and repaired, the term ”restoration” used here does not imply the repair of damages or defects.

2 Related Work

In recent years, many analyses of cultural properties have been conducted using 3D shape data. Such interdisciplinary research generally focuses on such analyses beyond the measurement and storage of 3D shape data, and the subjects vary from small-scale properties (Malzbender et al., 2001; Gilboa et al., 2013) to large-scale ones (Ikeuchi et al., 2007; Banno et al., 2008) over the world. One of the analysis researches using the shapes of cultural properties is virtual restoration. Since we focused on non-rigid virtual restoration, we introduced relevant studies to the virtual restoration of cultural properties and constrained non-rigid registration.

2.1 Virtual Restoration of Cultural Properties

There is a global demand for the use of 3D data to restore data and estimate the original appearances of various cultural properties. In particular, it is beneficial when it is difficult to physically restore cultural properties. There have been many studies on such virtual restorations in various fields, such as digital heritage and digital archaeology.

Papaioannou et al. (2002) proposed a 3D object restoration method through processing, where a depth map is used after segmenting an object surface. Koller &Levoy (2005) proposed a method where the local features and directions are utilized as metadata and the match between them is calculated to search for adjacent debris members and realize the assembly of partial objects. Also, Huang et al. (2006) proposed a method to calculate the 3D features of surfaces after segmenting object surfaces so as to match the necessary components with correspondences. Matching using the features of the 3D shape of an object surface has often been used for restoration (Brown et al., 2008; Huang et al., 2006; Vendrell-Vidal & Sánchez-Belenguer, 2014). In addition to the feature-based method, there was also a research that proposed the restoration of finer objects by providing guides based on template shapes (Zhang et al., 2015). In addition, methods using spectrum analysis (Altantsetseg et al., 2014), machine learning (Toler et al., 2010; Funkhouser et al., 2011), and the genetic algorithm (Sizikova & Funkhouser, 2017) have been proposed.

Through these methods, it is possible to restore the cultural properties that are composed of fragments, which can be regarded as rigid bodies. However, restoration through these methods is not feasible for the cultural assets that cause non-rigid deformation. It is not the case that the study has not been conducted on cultural properties that have been deformed. A research has previously analyzed a ship’s restoration condition by combining related historical data and shipbuilding knowledge (Jones et al., 2013), where the distortions that occurred in the physical restoration were analyzed based on the drawing information in (Nayling & Jones, 2012).

In many existing studies, virtual restoration was attempted by matching objects based on the feature quantities of 3D shapes. In these methods, virtual restoration is impossible in the situations in which each part is deformed. In addition, it is possible to analyze deformation situations by utilizing them in case detailed records are kept about cultural properties. However, in a situation in which only limited information is left, a new method needs to be devised.

2.2 Constrained Non-rigid Registration

Non-rigid registration is performed by constraining the deformation of a property, where some assumptions that are valid in each problem are considered. There are various characteristics to be preserved, such as the surface shape, volume, elastic energy, and topology. It is also necessary to understand the nature of deformation and to limit the range when expressing the deformation of an object.

Rohlfing et al. (2003) proposed a non-rigid registration method for medical MR image registration. Constraint terms for volume preservation under the assumption that cellular tissue is not compressed, and smoothness constraint with bending energy added. Regarding the CG applications that deform 3D models, there is also a study that proposed natural morphing by imposing constraints on such deformations so that they do not change the volume before and after deformation (Lan et al., 2017). Lin et al. (2013) proposed a method that limits deformation by considering the analogy of a spring in which a triangular mesh is put on an image and an elastic force acts between the vertices of a mesh, making it possible to change the topology and maintain it. Also, Igarashi et al. (2005) proposed an as-rigid-as-possible (ARAP) deformation system that can be operated in real time by designing optimization problems to be linear using the quadric error metric. (Sorkine & Alexa, 2007) applied the ARAP strategy in 3D modeling.

In the alignment task of two-point clouds, a method that falls into the probability density estimation problem of the Gaussian mixture model was proposed (Myronenko & Song, 2010; Horaud et al., 2011). The assumption for the input is that the global shapes are similar among the input. Ma et al. (2016) introduced a process in which the local shape is maintained and realized non-rigid registration for various deformations between various input point clouds. The rigid-as-rigid strategy is also practical in the non-rigid the alignment task (Fujiwara et al., 2011). Numerous non-rigid registration methods have been proposed (Zeng et al., 2010; Chang et al., 2020; Madsen et al., 2020; You et al., 2018; Zampogiannis et al., 2021; Ma et al., 2016; Yang et al., 2018; Hirose, 2021), including learning-based approaches (Wang et al., 2019; Feng et al., 2021; Shimada et al., 2019; Malik et al., 2020), but most of them are aimed at solving the correspondence problem between objects.

Our study is about mutually deforming multiple adjacent components, and there is no reference object to be matched. Since the specific target substance is wood, we introduce a deformation model to express the deformation properties of wood. We aimed at describing the deformation with low-dimensional variables (Salzmann et al., 2007) that reduce the ambiguous solutions in the optimization process as much as possible. In the optimization process, constraints were also imposed on the deformation, taking into account the physical properties of wood. Non-rigid assembly is performed by solving the nonlinear optimization problem. Also, using digital cultural property data, the virtual restoration results were presented in this paper.

3 Ruled Surface FFD

In this study, we propose a constrained deformation method known as the Ruled Surface FFD (RS-FFD), that considers the directionality of the deformation. As previously mentioned, the deformation of timber can be effectively represented by a ruled surface. Our deformation method generates a ruled surface that fits the point-cloud of an object; the point-cloud is deformed indirectly through the ruled surface. Figure 5 shows the processes of RS-FFD. We describe the details of the proposed method in this section. Hereafter, except for in the overlap detection process, all points and objects are assumed to be in 3D space.

3.1 Ruled Surface Representation

We fit a ruled surface to the object by considering its geometric shape, and the anisotropy of the deformation is determined by the grain of the wood. Because the wood is cut in a fixed direction to the grain, the shape of the object is strongly related to the deformation model. We assume that the ruled surface \({{\varvec{r}}}(\beta , \gamma )\) is composed of a base curve \({{\varvec{d}}}(\beta )\) with the parameter \(\beta \in \mathbb {R}\) and the union of lines \(\gamma {{\varvec{g}}}(\beta )\) that move on the curve:

$$\begin{aligned} {{\varvec{r}}}(\beta , \gamma ) = {{\varvec{d}}}(\beta ) + \gamma {{\varvec{g}}}(\beta ), \end{aligned}$$
(1)

where \({{\varvec{g}}}(\beta )\) is the unit vector at a point \({{\varvec{d}}}(\beta )\), \(\gamma \in \mathbb {R}\) represents a position on the line. \({{\varvec{d}}}(\beta )\) is called the directrix, and \({{\varvec{g}}}(\beta )\) is called the generator, as shown in Fig. 6.

The directrix is estimated by fitting a B-spline curve to a medial skeleton of the object, and the generator is estimated from the spatial distribution of the object by applying principal component analysis (PCA) to the point-cloud. The flow of the process is also shown in the figure. Given the ruled surface, every point on the object is associated with the surface. The object deforms through the ruled surface by moving the control points of the spline curve and the twist angles, which is explained later.

Fig. 5
figure 5

Deformation flow. A ruled surface is fitted to the scanned point cloud, and each point is associated with the ruled surface. The point cloud is deformed by changing the parameters of the ruled surface

Full size image
Fig. 6
figure 6

Ruled Surface Estimation. A skeleton is extracted from the input shape data, and the Directrix is obtained by fitting a spline curve to the extracted skeleton. The Generator is estimated from the distribution direction obtained by applying PCA to the input data

Full size image

3.1.1 Directrix Estimation

For estimating the directrix, a medial skeleton of the point-cloud is first found. We adopted the method proposed by Huang et al. (2013), which gives a discrete point set \(X^{\textrm{skl}} = \{ {{\varvec{x}}}_{\phi }^{\textrm{skl}} \in \mathbb {R}^{\textrm{3}} \}\) \((\phi = 1, 2, ..., n_\textrm{s})\) that represents the medial skeleton, where \( n_s \in \mathbb {N}\) is the number of points of the skeleton.

Our method fits a \( \rho \)-order spline curve \({{\varvec{c}}}(\tau ) (\tau \in [\tau _{\rho },\) \( \tau _{n_\textrm{c} - \rho -1}] \subset \mathbb {R})\) to the obtained point set, where \( n_c \in \mathbb {N}\) is the number of control points of the spline curve. For discrete points \(\hat{\tau }\in \mathbb {R}\), \( {{\varvec{c}}} (\hat{\tau }) \) is a discrete \( \rho \)-order spline curve with the B-spline basis function \( b_{i,\rho }(\hat{\tau })\). The open-uniform knot vector \( \tau _i \) and the control points \(P = \{ {{\varvec{p}}}_{i} \in \mathbb {R}^{\textrm{3}} \}\) (\( i = 1, 2, ..., n_\textrm{c}\)) define the curve as shown in the following equation:

$$\begin{aligned} {{\varvec{c}}}(\hat{\tau })= & {} \sum _{i}^{n_\textrm{c}} b_{i,\rho } (\hat{\tau }) {{\varvec{p}}}_i, \end{aligned}$$
(2)
$$\begin{aligned} b_{i,1}(\hat{\tau })= & {} {\left\{ \begin{array}{ll} 1 &{} (\tau _i \le \hat{\tau } \le \tau _{i+1}) \\ 0 &{} (\textrm{otherwise}) \end{array}\right. },\nonumber \\ b_{i,\rho (>1)}(\hat{\tau })= & {} \frac{(\hat{\tau }-\tau _i)}{\tau _{i+\rho -1}-\tau _i}b_{i,\rho -1}(\hat{\tau })\nonumber \\{} & {} + \frac{(\tau _{i+\rho } - \hat{\tau })}{\tau _{i+\rho }- \tau _{i+1}}b_{i+1,\rho -1}(\hat{\tau }). \end{aligned}$$
(3)

Both end points of the principal component direction of the skeleton points are set as the end points of the spline curve, and the discrete points of the spline curve are placed between them at equal intervals. The control points P are estimated by minimizing the distance between \( {{\varvec{c}}}(\hat{\tau }) \) and \(X^{\textrm{skl}}\),

$$\begin{aligned} \bar{P} = \mathop {\mathrm{arg~min}}\limits _{P} \sum _{j}^{n_\textrm{d}} \Vert {{\varvec{c}}}(\hat{\tau }_j) - {{\varvec{x}}}_{\textrm{closest}, j}^{\textrm{skl}}\Vert , \end{aligned}$$
(4)

where \(n_\textrm{d} \in \mathbb {N}\) is the number of discrete points, and \({{\varvec{x}}}_{\textrm{closest}, j}^{\textrm{skl}} \in \mathbb {R}^3\) is the point closest to \({{\varvec{c}}}(\hat{\tau }_j)\). The spline curve \( {{\varvec{c}}}(\hat{\tau }_j)\) is the directrix of the ruled surface.

3.1.2 Generator Estimation

We assume that the generator \({{\varvec{g}}}(\beta )\) is perpendicular to the directrix and parallel to the surface. We applied PCA to the point-cloud of the object to obtain the surface direction. The generator is defined by a principal component vector \( \varvec{v} \in \mathbb {R}^3\) that is nearly perpendicular to the curve \( {{\varvec{d}}}(\beta ) \) and along the surface direction. \( {{\varvec{v}}} \) is not precisely perpendicular to the curve \( {{\varvec{d}}}(\beta ) \) at any points on the curve; \({{\varvec{g}}}(\beta )\) is obtained by Gram-Schmidt orthonormalization, as follows:

$$\begin{aligned} {{\varvec{g}}}(\beta ) = \displaystyle \frac{{{\varvec{v}}} - \left( \displaystyle \frac{d}{d\beta } {{\varvec{d}}}(\beta ) \cdot {{\varvec{v}}}\right) \displaystyle \frac{d}{d\beta } {{\varvec{d}}}(\beta )}{\left\| {{\varvec{v}}} - \left( \displaystyle \frac{d}{d\beta } {{\varvec{d}}}(\beta ) \cdot {{\varvec{v}}}\right) \displaystyle \frac{d}{d\beta } {{\varvec{d}}}(\beta ) \right\| }. \end{aligned}$$
(5)

3.2 Ruled Surface FFD

We propose a deformation method based on the ruled surface representation. We define the orthogonal matrix \({{\varvec{L}}}(\beta )\) that determines the FFD coordinates on a ruled surface. Note that \(\det (L)=1\). \({{\varvec{L}}}(\beta ) \) is represented on a ruled surface as follows:

$$\begin{aligned} {{\varvec{L}}}(\beta ) = \left( \begin{array}{c} {{\varvec{l}}}_1 (\beta )\\ {{\varvec{l}}}_2 (\beta )\\ {{\varvec{l}}}_3 (\beta ) \end{array} \right) = \left( \begin{array}{l} \displaystyle \frac{d}{d\beta } {{\varvec{d}}}(\beta )\\ {{\varvec{g}}}(\beta )\\ {{\varvec{l}}}_1(\beta ) \times {{\varvec{l}}}_2(\beta ) \end{array} \right) . \end{aligned}$$
(6)

\({{\varvec{l}}}_1\), \({{\varvec{l}}}_2\), and \({{\varvec{l}}}_3\) are the orthonormal basis. The local coordinates (s, t, u) of a point \({{\varvec{x}}}_c \in \mathbb {R}^3\) in the point-cloud are defined by the nearest point \({{\varvec{d}}}(\beta _c)\) on the curve as the origin, and the base coordinate \( {{\varvec{L}}} (\beta _c)\), as follows:

$$\begin{aligned} {{\varvec{x}}}_c = {{\varvec{d}}}(\beta _c) + (s,t,u) \cdot {{\varvec{L}}}(\beta _c). \end{aligned}$$
(7)
Fig. 7
figure 7

Update of Ruled Surface. The directrix is updated by moving the control points \(p_i\), and the generator is rotated by the control angle \(\alpha _i\)

Full size image

3.2.1 Twist Deformation

Because the above representation cannot twist the object, we introduce the angle \( \alpha _i \) associated with the control point. The parameter set \(A = \{ \alpha _i \in \mathbb {R} \}\) (\( i = 1, 2, ..., n_\textrm{c}\)) controls the extent of twist required for the control points to affect the surrounding points. To deform the object smoothly, the twisting angle \( \alpha (\beta ) \) associated with the directrix is defined by the weighted average according to the inverse ratio of the distance from the control points, as follows:

$$\begin{aligned} \alpha (\beta ) = \sum _i^{n_c} \alpha _i\frac{{||{{\varvec{d}}}(\beta ) - {{\varvec{p}}}_i||^{-1}}}{\sqrt{\sum _\iota ^{n_c}{||{{\varvec{d}}}(\beta ) - {{\varvec{p}}}_\iota ||}^{-2}}}. \end{aligned}$$
(8)

In this equation, \(\alpha (\beta )\) represents the angle that rotates a local coordinate \( {{\varvec{l}}} _3(\beta ) \) around \( {{\varvec{l}}}_1 (\beta )\), which is the tangential direction of the curve. The generator is redefined with a rotation matrix \({{\varvec{R}}}\) as:

$$\begin{aligned} \hat{{{\varvec{g}}}} (\beta ) = {{\varvec{g}}} (\beta ) {{\varvec{R}}} (\alpha (\beta ), {{\varvec{l}}}_{1} (\beta )). \end{aligned}$$
(9)

The matrix \({{\varvec{R}}} (\theta , {{\varvec{a}}}) \) is the rotation around the axis \( {{\varvec{a}}} = (a_1, a_2, a_3) \) by \(\theta \). The local coordinate is also redefined as follows:

$$\begin{aligned} \hat{{{\varvec{L}}}}(\beta ) = \left( \begin{array}{l} \hat{{{\varvec{l}}}_1} (\beta )\\ \hat{{{\varvec{l}}}_2} (\beta )\\ \hat{{{\varvec{l}}}_3} (\beta ) \end{array} \right) = \left( \begin{array}{l} \displaystyle \frac{d}{d\beta } {{\varvec{d}}}(\beta )\\ {{\varvec{g}}} (\beta ) {{\varvec{R}}}(\alpha (\beta ), \hat{{{\varvec{l}}}_{1}}(\beta )) \\ \hat{{{\varvec{l}}}_1} (\beta ) \times \hat{{{\varvec{l}}}_2} (\beta ) \end{array} \right) . \end{aligned}$$
(10)

3.2.2 Anisotropic Deformation

Timber deformation occurs anisotropically. Here, we introduce the expansion and contraction rate \( m \in \mathbb {R}\), and \( r_\textrm{2} \in \mathbb {R}\) and \( r_\textrm{3} \in \mathbb {R}\) are constant expansion and contraction rates in the radial and tangential directions to the directrix (fiber) direction. The ratio of the fiber direction \( r_\textrm{1} \) is set to 1 and omitted here. m is calculated from the lengths of the initial and deformed curves, and \( r_\textrm{2} \) and \( r_\textrm{3} \) are given as the prior knowledge. Finally, the point \(\hat{{{\varvec{x}}}_c}\) of the object is represented as follows:

$$\begin{aligned} \begin{aligned} \hat{{{\varvec{x}}}_c}&= {{\varvec{d}}}(\beta _c) + s \hat{{{\varvec{l}}}_{1}}(\beta _c)\\&\quad + (1 + r_\textrm{2} m) t \hat{{{\varvec{l}}}_{2}}(\beta _c) + (1 + r_\textrm{3} m) u \hat{{{\varvec{l}}}_{3}}(\beta _c). \end{aligned} \end{aligned}$$
(11)

The spline curve is updated by changing the control points and angles, as shown in Fig. 7. Accordingly, the ruled surface defined above is also updated, and the point-cloud associated with it deforms by recalculating \(\hat{{{\varvec{x}}}_c}\).

4 Non-rigid Assembly with RS-FFD

In this section, we present a method for assembling partial objects using RS-FFD. The framework is a non-rigid alignment of multiple objects by a constrained deformation model. We assume that the order of the adjacent components is known and that approximate correspondences are given. Specifically, it is known that the boat was built by fixing the components with ropes. We can find the positions of the holes, and the closest holes between components can be the corresponding points for alignment.

The solution must be physically consistent because virtual restoration is the recreation of the production process, and accuracy is required over appearance. To prevent spatial overlap of the components, we introduce a distance metric considering the collisions between them such that the distance of the corresponding points increases according to the overlap. We also applied constraints to the amount of deformation with criteria for obtaining practical solutions. The constraints can also narrow the searching space for optimization.

4.1 Overlap-Aware Distance Metric

Here, we define the distance metric for spatially assembling multiple 3D objects. In actual conditions, objects cannot overlap. Therefore, when assembling the components, it is necessary to design a distance metric to avoid spatial overlap. If the shape of the contact surface is correctly obtained, a 3D collision detection method can be applied. In some cases, however, the shape is incomplete, or only the surface shape is obtained, such as the present case. Therefore, based on the assumption that adjacent components are smoothly connected on the surface, the overlap is detected by projecting the surface shape onto two-dimensional (2D) planes.

Fig. 8
figure 8

Overlap detection by winding number. a The overlap between two components near a point is detected by projecting the surrounding points onto a 2D plane. b The point which causes the overlap is penalized in the distance calculation. c The winding number \(n_w \in \mathbb {Z}\) is not equal to zero if a point \({{\varvec{x}}}_m \in \mathbb {R}^2\) is inside the boundary of another component \(B=\{{{\varvec{b}}}_l \in \mathbb {R}^2\}\). \({{\varvec{d}}} \in \mathbb {R}^2\) is an arbitrary vector from a point \({{\varvec{x}}}_m\) that crosses the line between \({{\varvec{b}}}_l\) and \({{\varvec{b}}}_{l+1}\). If \({{\varvec{d}}}\) crosses a boundary line in counter clockwise direction, \(n_w\) increments (\(n_w \leftarrow n_w + 1\)); if it crosses in clockwise direction, \(n_w\) decrements (\(n_w \leftarrow n_w - 1\)). d The winding number is zero if a point \({{\varvec{x}}}_m\) is outside the boundary

Full size image

4.1.1 Overlap Detection on 2D Plane

As previous noted, overlap between two components near a point is detected by projecting the surrounding points onto a 2D plane. Assume that the corresponding point of a point \( {{\varvec{x}}} \in \mathbb {R}^3\) is \(\acute{ {{\varvec{x}}}} \in \mathbb {R}^3\) in another component. A plane is fit to a set of points that are within a distance \(\lambda \) of \(\acute{{{\varvec{x}}}}\). By finding the alpha-shape of the projected points and calculating the concave hull, the boundary points of the component can be obtained. Then, the overlap between components can be detected by calculating the winding number (Sunday, 2021), as shown in Fig. 8. The boundary points of a component are regarded as a closed curve, and the winding numbers for the points in another component are calculated. If a point with a non-zero winding number is detected, point \( {{\varvec{x}}} \) is considered to be overlapped.

4.1.2 Overlap-Aware Distance Metric

When overlap of components is detected near the corresponding points, the distance between the points is penalized, which increases the distance cost between the adjacent components and eliminates physically impossible solutions. Let \( p (> 1) \in \mathbb {R}\) be a penalty constant; the distance \( d ({{\varvec{x}}}, \acute{{{\varvec{x}}}})\) between point \( {{\varvec{x}}} \) and point \(\acute{{{\varvec{x}}}}\) is defined simply as follows:

$$\begin{aligned} d({{\varvec{x}}}, \acute{{{\varvec{x}}}}) = {\left\{ \begin{array}{ll} p \Vert {{\varvec{x}}} - \acute{{{\varvec{x}}}} \Vert &{} (\mathrm {overlap-case}) \\ \Vert {{\varvec{x}}} - \acute{{{\varvec{x}}}} \Vert &{} (\textrm{otherwise}) \end{array}\right. }. \end{aligned}$$
(12)

Using the above distance metric, we define the distance cost for alignment. The target object has \(n_q \in \mathbb {N}\) components \( H = \{\eta _q\} (q = 1, 2, ..., n_q) \). Let \( n_p \in \mathbb {N}\) be the number of pairs of adjacent components and \( n_{<\eta , \acute{\eta }>} \in \mathbb {N}\) be the number of corresponding points between components \( \eta \) and \( \acute{\eta } \). Then, the total distance cost for alignment \(\varepsilon _{\textrm{reg}} \in \mathbb {R}\) is defined as the summation of the distance between all corresponding points as follows:

$$\begin{aligned} \begin{aligned} \varepsilon _{\textrm{reg}}= \frac{1}{n_p} \sum _{<\eta , \acute{\eta }>}^{n_p} \frac{1}{n_{<\eta , \acute{\eta }>}} \sum _{\kappa }^{n_{<\eta , \acute{\eta }>}} d({{\varvec{x}}}_{<\eta , \acute{\eta }>, \kappa }, \acute{{{\varvec{x}}}}_{<\eta , \acute{\eta }>, \kappa }). \end{aligned} \end{aligned}$$
(13)

4.2 Constraints on Deformation

In this section, we describe the constraints imposed to the deformation by considering the physical properties of the wood. The constraints for deformation consist of the following three terms.

  1. (i)

    Constraint on total elasticity

  2. (ii)

    Constraint on uniform elasticity

  3. (iii)

    Constraint on bending

Figure 9 outlines these constraints visually. These constraints are derived from the nature of wood such that it is a biological material composed of aggregation of many cells. The total elasticity is a limitation of scaling and is given by the moisture content and the scaling factor associated with it. Because wood is composed of small cells, it does not expand and contract locally. The uniform elasticity is maintained by constraints of the small differences in local length of the objects. Although wood is a deformable material, its fracture strength is exceeded when deformation occurs at an angle greater than a particular level. Here, we introduce bending constraints to suppress such deformation. Notably, the deformable area for bending changes depending on boundary conditions and processing techniques (Fig. 9).

Fig. 9
figure 9

Constraints on deformation. a The total elasticity is the constraint for overall scaling given by the moisture content. b The uniform elasticity guarantees that the local deformation is uniform over the object. c The smooth bending is introduced to ensure that the fracture strength is not exceeded

Full size image

4.3 (i) Constraint on Total Elasticity

Here, we define the conditions for the total elasticity components. The condition on the length of the component is defined using the length of the directrix \( \xi _\eta \in \mathbb {R}\), which is calculated as the sum of the lengths between the discrete points \(\{ {{\varvec{c}}}_{\eta , j} \}\), as follows:

$$\begin{aligned} \xi _\eta = \sum _j^{n_{d, \eta } - 1} \Vert {{\varvec{c}}}_{\eta ,j} - {{\varvec{c}}}_{\eta ,j+1} \Vert . \end{aligned}$$
(14)

Assuming that the length of the initial state obtained similarly is \( \tilde{\xi _\eta } \), the global length variation is defined as \( \varDelta \xi _\eta = |\frac{\xi _\eta }{\tilde{\xi _\eta }} - 1 |\). We compare \( \varDelta \xi _\eta \) to the threshold of the range condition and set it as a constraint for total elasticity.

4.4 (ii) Constraint on Uniform Elasticity

The uniform elasticity \(\zeta _{\eta ,j}\in \mathbb {R}\) of a discrete point \({{\varvec{c}}}_{\eta ,j}\) is evaluated by the ratio of the distance between adjacent discrete points on the directrix to the total length of the initial directrix.

$$\begin{aligned} \zeta _{\eta ,j} = \frac{\Vert {{\varvec{c}}}_{\eta ,j} - {{\varvec{c}}}_{\eta ,j+1} \Vert }{\xi _\eta } \end{aligned}$$
(15)

Let \( \tilde{\zeta }_{\eta , j} \) be the ratio of the initial states. The local length variation is defined as \( \varDelta \zeta _{\eta , j} = |\zeta _{\eta ,j} - \tilde{\zeta }_{\eta ,j}|\). We compare \( \varDelta \zeta _{\eta , j} \) with a threshold of the range condition and set it as a constraint for uniform elasticity.

4.5 (iii) Constraint on Bending

The constraint on bending is evaluated by the angular difference of the base local coordinate of adjacent discrete points on the directrix from the initial state. The Euler’s angle difference between the base coordinates of adjacent discrete points \( {{\varvec{c}}}_{\eta , j} \) and \( {{\varvec{c}}}_{\eta , j + 1} \) is then evaluated. The basis \( {{\varvec{l}}}_{\eta , j, 3} \), \( {{\varvec{l}}}_{\eta , j, 2} \) and \( {{\varvec{l}}}_{\eta , j, 1} \) are deployed by rotating \( \theta _{\eta , j^*, 3} \), \( \theta _{\eta , j^*, 2} \), \( \theta _{\eta , j^*, 1} \) in this particular order. The rotation matrix from base on \({{\varvec{c}}}_{\eta , j} \) to \({{\varvec{c}}}_{\eta , j + 1}\) can be represented as follows:

$$\begin{aligned} \begin{aligned} {{\varvec{R}}}_{\eta ,j^*}&= \\&{{\varvec{R}}}(\theta _{\eta ,j^*,3}, {{\varvec{l}}}_{\eta ,j,3}) {{\varvec{R}}}(\theta _{\eta ,j^*,2}, {{\varvec{l}}}_{\eta ,j,2}) {{\varvec{R}}}(\theta _{\eta ,j^*,1}, {{\varvec{l}}}_{\eta ,j,1}) . \end{aligned} \end{aligned}$$
(16)

Moreover, the rotation angle can be expressed as Eq. (17) using the local FFD coordinates.

$$\begin{aligned} \begin{aligned} {{\varvec{R}}}_{\eta ,j^*}&= ({{\varvec{L}}}_{\eta ,{j}})^\textrm{T} {{\varvec{L}}}_{\eta ,{j+1}} \end{aligned} \end{aligned}$$
(17)

By solving Eqs. (16) (17) with a limited range of angles, we find the rotation angles \(\theta _{\eta , j^*, l} (l = 1, 2, 3)\). Let \(\tilde{\theta }_{\eta , j^*, l}\) be the angular difference in the initial state. The rotation angle after deformation \(\theta _{\eta , j^*, l} \) is compared to evaluate the angular difference between adjacent discrete points, which is defined as \( \varDelta \theta _{\eta , j^*, l} = |\theta _{\eta ,j^*,l} - \tilde{\theta }_{\eta ,j^*,l}|\).

4.6 Conditioning by Constraints

The constraints described above were implemented as conditions in the optimization. The constraints need to work similar to thresholding such that the cost must significantly increase when the conditions are not satisfied during optimization. We adopted the sigmoid function \(\varsigma _a(x) = \frac{1}{1 + e^{-ax}}\) for the constraints. If the threshold is violated, the cost increases dramatically, whereas the cost is close to zero if the conditions are satisfied.

The cost by the constraint terms for optimization are expressed as,

$$\begin{aligned} \varepsilon _{\textrm{cons}} = w_{\textrm{tel}}f_{\textrm{tel}} + w_{\textrm{uel}}f_{\textrm{uel}} + w_{\textrm{bnd}}f_{\textrm{bnd}}, \end{aligned}$$
(18)

where \(f_{\textrm{tel}}\), \(f_{\textrm{uel}}\), and \(f_{\textrm{bnd}}\) are the constraint terms for total elasticity, uniform elasticity, and bending smoothness, respectively, and \(w_{\textrm{tel}}\), \(w_{\textrm{uel}}\), and \(w_{\textrm{bnd}}\) are the respective weights for each term.

4.7 (i) Constraint Term for Total Elasticity

The evaluation function is defined by the degree of change \( \varDelta \xi _\eta \) from the initial state with a threshold \( \xi _\textrm{th} \in \mathbb {R} \). The sum of all components is the overall evaluation function.

$$\begin{aligned} \begin{aligned} f_\textrm{tel}=&\sum _{\eta }^{n_q} \varsigma _a (\varDelta \xi _\eta - \xi _\textrm{th}) \end{aligned} \end{aligned}$$
(19)

This term controls the ratio of expansion and contraction of the components.

4.8 (ii) Constraint Term for Uniform Elasticity

Uniform elasticity is defined as the sum of the change amount \( \varDelta \zeta _{\eta , j}\) from the initial state of the local length change of the directrix with a threshold \( \zeta _{\textrm{th}} \in \mathbb {R}\).

$$\begin{aligned} \begin{aligned} f_\textrm{uel} =&\sum _{\eta }^{n_q} \sum _{j}^{n_{d,\eta } - 1} \varsigma _a (\varDelta \zeta _{\eta ,j} - \zeta _{\textrm{th}}) \end{aligned} \end{aligned}$$
(20)

This term suppresses the dramatic changes in local length and enables uniform length change.

4.9 (iii) Constraint Term for Bending

The constraint term for bending takes the sum of the angular difference of the local coordinate based on the adjacent discrete points \( \theta _{\eta , j^*, l} \) as the evaluation function. This term limits bending and twisting and achieves smooth deformation.

$$\begin{aligned} \begin{aligned} f_\textrm{bnd} =&\sum _{\eta }^{n_q} \sum _{j}^{n_{d, \eta } - 1} \sum _{l}^3 \varsigma _a (\varDelta \theta _{\eta ,j^*,l} - \theta _{\textrm{th}}) \end{aligned} \end{aligned}$$
(21)

4.10 Optimization

Non-rigid assembly is performed by an optimization framework that minimizes the cost of the distance between the correspondences and the constraints. The variables are the set of control points \( \tilde{P} = \{P_q\}\) (\(q = 1, 2, ..., n_q\)) and the set of twist angles \( \tilde{A} = \{A_q\}\) (\(q = 1, 2, ..., n_q\)).

$$\begin{aligned} \{\tilde{P}, \tilde{A} \} =\mathop {\mathrm{arg~min}}\limits _{\tilde{P}, \tilde{A}} (\varepsilon _{\textrm{reg}} + \varepsilon _{\textrm{cons}}) \end{aligned}$$
(22)

We used the Levenberg–Marquardt method (Levenberg, 1944; Marquardt, 1963) for solving the nonlinear least squares problem.

5 Experiment

5.1 Scanning and Data Processing

Here, we described the pre-processing of the first boat’s scan data. The first boat is currently on display at the Giza Solar Boat Museum, and we measured the boat using a laser scanner (Z+F IMAGER 5010C) at the museum. The overall shape obtained by aligning the partial 3D data (Oishi et al., 2005) is shown in Fig. 10. Only the visible surface of the boat (outer shape and inner measurable parts) was scanned. The point cloud (approx. 1.4 million points) of the entire boat, in which multiple scan data were merged, was extracted for each component, and the component’s data subjected to segmentation and down-sampling was treated as the surface data. The surface data was up to 30,000 points per component.

Fig. 10
figure 10

Partial 3D shape data of the main body. The whole body was decomposed into individual pieces of timber, resulting in 26 pieces of partial data

Full size image

The hull data was manually segmented into each component as pre-processing, as shown in Fig. 10, and the data of each component had an outer shape and an inner partial shape. There were many missing areas, especially at the inner side of the boat, because of the occlusions and limited spaces for scanning. The shapes of the surfaces where the components are adjacent could not be measured, so we only used the outer surface for the experiment.

5.2 Corresponding Points

The components of the first boat’s hull are not fixed by nails but by ropes and tenons. It is thought that a plurality of holes were made so that they can be pulled from the side to the back of the member, where ropes are passed through the holes, and small pieces of wood are inserted into the mortises for fixing. The components that make up the hull have a plurality of holes spaced apart, as shown in Fig. 11. In this research, the positions of these holes were considered to be corresponding points, so we decided to use manually selected points at fixed intervals as corresponding points for alignment.

Fig. 11
figure 11

Fixing method of components. The components are fixed by ropes and tenons. We assume that the positional correspondences of the components are given by the past archaeological inspections

Full size image

5.3 Parameters

  • parameters of the deformation model

The spline curve that describes the object’s deformation was set as dimension 3, the number of control points was 4, and the knot vector was uniform ({0, 0, 0, 1, 1, 1, 1 }). The number of discrete points on the directrix was 100. The scaling ratio of anisotropy was set based on (Sadoh, 1996) so that the ratio of the fiber direction, radial direction, and tangential direction is \( \textrm{r}_1 : \textrm{r}_2 : \textrm{r}_3 = 1: 5: 10 \).

  • parameters of the non-rigid assembly

When calculating the distances between the corresponding points, the range for plane fitting and overlap detection was set to \( \lambda = 0.5 [\textrm{m}] \). This is the same distance as the width of the component, and it was set because it can cover the point clouds of a sufficient number of members for the plane estimation. The constant for the distance penalty was \( p = 5 \). The range conditions for the deformation constraints were 0.5 % for total elasticity, 5 % for uniform elasticity, and \( \theta _\textrm{th} = 0.05 \) for bending, and the sigmoid function parameter was \( a = 100 \). The weight of the constraint terms was \( w _ {\textrm{tel}} = w _ {\textrm{uel}} = w _ {\textrm{bnd}} = 1 ^{(*)}\).

*Weights of constraints (\( w_{\textrm{tel}} \), \( w_{\textrm{uel}} \), \( w_{\textrm{bnd}} \)) : The alignment term unit \( \varepsilon _{\textrm{reg}} \) is a distance [m] that is literally weighted but still represents the Euclidean distance. However, each weight was also introduced to absorb the differences in the dimensions of the constraint terms, where each constraint term is a dimensionless quantity because it is expressed as a simple sum of sigmoid functions. The constraint term plays a thresholding role and linearly increases (\( 0 \rightarrow 1 \)) around the threshold for the combinations of the control point \( {{\varvec{p}}}_i \) and control angle \( \alpha _i \) that cause deformation beyond the threshold.

Since the constraint term is an implementation of the range condition, the term is preferable to be \( \varepsilon _{cons} \approx 0 \). Therefore, w needs to be sufficiently large for the error distance \( \varepsilon _\textrm{reg} \) [m] to make the constraint term function a range condition. In this research, the error is in the order of [cm], so the weight of the constraint term was set as 1 for all the conditions.

5.4 Experiment 1 : Virtual Reconstruction by Non-rigid Assembly

The purpose of the experiment was to examine the possibility that alternative solutions to the overall shape of the boat can be obtained through rerestoring to fill in the local distortions and gaps. Here, the state of the boat at the scanning time was used as the initial state. The components constituting the bottom of the boat were treated as a one-piece object. In addition, optimization was performed by global optimization with the control points and control angles of all the components as variables. Table 1 summarizes the parameters used in the experiment. As for the total elasticity, the upper limit was set with reference to the shrinkage rate at a moisture content of 10%. Also, threshold values were empirically set for local elasticity and smooth bending.

Table 1 parameter setting
Full size table

Figure 12 shows enlarged images of two representative areas where changes in the shape could be observed before and after the restoration. The white area above the deck was not involved in the optimization; the area was just shown for reference. As shown in the upper part of the figure, the gaps between the components decreased. In the lower part of the figure, it can be seen that the gaps between the distorted parts, where the components did not engage with each other, were narrowed by the deformation method.

In Fig. 13, the distances between the corresponding points decreased, and the optimization converged. The proposed method reduced the distance between the components, and the effectiveness of the method in reducing the distortion of the hull was confirmed. However, the distance between the corresponding points did not perfectly match, and the reason for this may be that the corresponding points were manually given. It was also thought that there was a limitation in the representation of the deformation model.

Fig. 12
figure 12

Virtual restoration results of the first boat. a We can observe the significant gaps between components in the 3D data before optimization. b The gaps between the components decreased after optimization while avoiding overlaps

Full size image
Fig. 13
figure 13

Average distance during optimization

Full size image

5.5 Experiment 2 : Reconstruction from Flat Timbers

The experiment was conducted to apply the proposed method by experimentally creating the state before the component assembly. First, we experimentally reproduced the state of the flat timbers from which the components of the first boat were cut out, as shown in Fig. 14. Then, the control points of the directrix were moved to make the skeleton curve of the components a straight line. Afterward, the control points of the original curve were projected onto a 3D line passing through the center of gravity of the original spline curve whose direction vector is the tangent vector at the central discrete point of the spline curve, which made it easy to handle the initial alignment of the components after flattening. The component was deformed by the proposed deformation method according to the projected control points, and the deformation and assembly parameters were the same as those of Ex. 1. It was also confirmed that the flat state was within the deformation condition range (Fig. 14).

The overall optimization results after the flattening of all the components and the distribution of the distances between the corresponding points were shown in Fig. 15. Also, the convergence of the optimization with the average distances of the correspondences was shown in Fig. 16. From the resulting figure, the optimization converges and the possibility of reproducing the boat’s shape can be seen. Moreover, the average distance after the optimization was smaller than that of Ex1. This shows that the optimization of the boat’s flat shapes made it possible to avoid the local minima and reach an optimal solution.

5.5.1 Experiments with and Without Constraints

To see the effect of the constraint terms, we also optimized the case without constraint terms. Figure 17b, e show examples with the constraint terms, and (c) and (f) without. As can be seen from the figure, even if the corresponding points are known, it is clear that the solution is not uniquely determined without the constraint. When we checked the length \(\xi _\eta \) of each component, we found that many of them exceeded the constraint value when no constraints were imposed, as shown in Table 2. In the case of not imposing constraints, the slight difference or number of corresponding points would cause the overall balance to be lost, resulting in incorrect solutions. On the other hand, in the proposed method, it was confirmed that the constraints were correctly designed and functional.

Figure 18 shows results with different parameters fot the constraint terms. Figure 18a shows the result of the no-twist deformation model: B-spline FFD. The deformed model without twist cannot recover the original shape at all. Figure 18b shows the case where the deformation is strongly restricted. In this case, the shape is restored to some extent, but there are many gaps between the components. These results also indicate that the constraint terms for deformation need to be set appropriately.

Fig. 14
figure 14

Flattening timbers. We experimentally reproduced the state of the flat timbers to simulate the actual production process. We flattened each component by moving the control points of the directrix to make the skeleton curve a straight line

Full size image
Fig. 15
figure 15

Result of non-rigid assembly from flat timbers. The figures on the left show the process of shape deformation in iterative optimization. The graphs on the right are histograms of the number of corresponding points against the distance between corresponding points

Full size image
Fig. 16
figure 16

Average distance between corresponding points

Full size image
Fig. 17
figure 17

Comparison with and without constraints. c, f In the case of not imposing constraints, the individual components will be deformed too much, and as a result, the overall shape is incorrect. b, e In the proposed method, the reconstructed shape is close to the current shape but slightly slimmer, indicating that the constraints were correctly designed and functional

Full size image
Fig. 18
figure 18

Comparison with different parameters. a is the result for Experiment 2 without twist parameters in the deformation model. If twist parameters are not used, the shape cannot be restored at all because the deformation of the timber cannot be adequately represented. b is a result with

Full size image
Table 2 Length of components after optimization with and without constraint terms
Full size table

6 Conclusion

In this paper, we proposed a non-rigid shape assembly method for the virtual restoration of the cultural properties composed of wooden materials. The proposed deformation method, Ruled-Surface FFD, is low-dimensional but has sufficient degrees of freedom for representing the deformation of wooden objects. The model allowed anisotropic deformation based on the knowledge of the nature of wood. We also proposed a framework for assembling multiple partial objects using a constrained non–rigid alignment method with an overlap-aware error metric. In the experiment, we applied the proposed method, where we scanned the data of a large wooden structure: the Solar Boat of King Khufu. The results showed that the proposed method works well in regard to fixing the distortions resulting from long-term exhibitions. We also confirmed that the method can restore the boat using the flat timbers from which it was cut. The result showed that the method also works well in this case, and the residual of the alignment was smaller than that of the above-mentioned experiment.

In the future, we aim at evaluating the restoration results from the viewpoint of archaeology. A complementary verification with the second solar boat would be also possible.