From Point Cloud Data to Building Information Modelling: An Automatic Parametric Workflow for Heritage

Building Information Modelling (BIM) is a globally adapted methodology by government organisations and builders who conceive the integration of the organisation, planning, development and the digital construction model into a single project. In the case of a heritage building, the Historic Building Information Modelling (HBIM) approach is able to cover the comprehensive restoration of the building. In contrast to BIM applied to new buildings, HBIM can address different models which represent either periods of historical interpretation, restoration phases or records of heritage assets over time. Great efforts are currently being made to automatically reconstitute the geometry of cultural heritage elements from data acquisition techniques such as Terrestrial Laser Scanning (TLS) or Structure From Motion (SfM) into BIM (Scan-to-BIM). Hence, this work advances on the parametric modelling from remote sensing point cloud data, which is carried out under the Rhino+Grasshopper-ArchiCAD combination. This workflow enables the automatic conversion of TLS and SFM point cloud data into textured 3D meshes and thus BIM objects to be included in the HBIM project. The accuracy assessment of this workflow yields a standard deviation value of 68.28 pixels, which is lower than other author’s precision but suffices for the automatic HBIM of the case study in this research.


General Framework
In the engineering sector, Building Information Modelling (BIM) is considered both an opportunity and a challenge. BIM entails significant time, economic resources and training of the human team, but it also changes the current hiring processes in the development and achievement of engineering and architectural projects. This effort is rewarded with the amount of skills and results that come from applying the BIM technology. Nevertheless, in the field of historic buildings, the use of BIM is limited because of the particularities and characteristics of the heritage, but its presence in the scientific literature is improving. The Heritage or Historic Building Information Modelling (HBIM) approach, which addresses a comprehensive and interdisciplinary management model in the field of heritage, combines building components edges; the validation of the methodology via accuracy assessment, the interface features and the texture generation. The workflow implementation within the BIM environment is then described. Finally, the results are discussed, and the conclusions and future work are presented.

Related Work
Optimising BIM using parametric tools is an issue addressed in the scientific literature. Khaja et al. [42] used parametric tools such as Dynamo [43], an integrated Visual Programming tool from Revit, to demonstrate its potential for automatic population in BIM models. In this line, Rahmani et al. [44] assessed the environmental aspects of the models. ("Data Acquisition and Visualization for IEQ Assessment: A case study of daylight field measurement -White Rose Research Online") worked on statistical simulation on indoor environmental quality, and Rahmani et al. [44] addressed a performance optimisation scenario based on BIM called 'BPOpt' (BIM-based performance optimization).
Most scientific papers on HBIM from point cloud data build these models through manual processes, as in the case study of the Church of Santa Maria at Portonovo (Ancona, Italy) [27]. Here, the objects do not represent the morphological properties from the point clouds; they are created with existing parametric models, but this approach is already described by Murphy et al. [28]. This is also the case of the work by Bassier et al. [29], who reconstructed the portico of the Saint-Jacobs Church (Leuven, Belgium) from an ideal model-without morphological variations-using 3D meshes for later finite element analysis (Finite Element Method, FEM [45]). Barazzetti et al. [46] worked in this line in Castel Masegra in Sondrio, Italy. Sztwiertnia et al. [40] created the HBIM project of the Wang Stave Church in Karpacz, Poland, and identified the difficulty of creating geometric shapes from point cloud data. Progress on adapting complex shapes to parametric elements is respectively made by Antón et al. [8], and Nieto, Antón and Moyano [5]. The former authors converted the point cloud data into as-built 3D solid objects as IFC entities to become part of the HBIM project, which enables the accurate analysis of the heritage city [9] by considering the geometrical alterations of the assets. The latter [5] used Boolean operations within the ArchiCAD®BIM platform to implement and manage structural deformations on walls. Therefore, this paper advances on the problem of representing diverse aspects such as structural deformations in the parametric objects, materials through textures, process changes in construction techniques and deformations over the course of time.
Procedures to create 3D parametric objects from point clouds should be mentioned. Barazzetti et al. [47] considered a semi-automatic approach, as it required successive steps to reach editable parametric objects starting from non-uniform rational B-splines (NURBS) and converted into solids. Recent work [48] developed a new workflow for mesh modelling towards HBIM, for which Autodesk Revit (BIM) and Dynamo are used to explore new HBIM features. However, several studies [49,50] implement these tools for specific objects, but not for building construction units. Finally, the work by Tommasi et al. [51] brings together diverse software in the market for reverse engineering in BIM applied to historical elements, whereas Pocobelli et al. [52] integrated organised moisture data into spreadsheets and associated it with parametric objects development through Dynamo. In addition, there are different reconstruction methods in the scientific literature [53][54][55][56][57][58][59][60][61], most of which use the Delaunay triangulation [54]. Surface reconstruction is the process of rebuilding geometry from the scanned point cloud. Ideally, the reconstructed surface mesh should have all the mesh faces matching the initial geometry from which the scan was taken.
In this paper, in order to create a new workflow to transfer point cloud data into the BIM environment, a pseudo point cloud structure with controllable density is created from the combination of Rhinoceros [62], Grasshopper [63]-a visual programming language for Rhinoceros-and ArchiCAD. The TLS and SFM data are then converted into textured 3D meshes. The important step is to avoid workflows with complementary software such as Cloudcompare [64], MeshLab [65] or Geomagic [66]. In this workflow, different algorithms are tested for meshing and automatic texture mapping onto the parametric objects. The process consists of converting the numerical model (point cloud) into 'Morph' elements (for walls and vertical surfaces) and meshes (for floors or horizontal surfaces). Next, these 3D objects are automatically transformed into BIM parametric objects (GSM file format), thus inheriting their dimensional and texture properties. In addition, the accuracy and quality of the model's parameters are evaluated to verify their consistency with an appropriate scale.
Regarding BIM operability with TLS or SFM data, the latest versions of Graphisoft ArchiCAD ® [23] and Autodesk Revit [34] allow to import these point clouds in TXT or XYZ file formats. This enables the BIM creation by fitting the morphological alterations under diverse modelling approaches. However, this paper advances on the automatic conversion of point cloud data into a parametric model by defining shapes through mathematical functions [51] under the Rhino+Grasshopper-ArchiCAD software combination. This workflow interconnects algorithms capable of modelling 3D objects so that the latter can be interoperable, which contributes to scientific knowledge in this field.

Case Study: The Main Patio Of The Casa De Pilatos
This research is based on a case study selected for its characteristics: a complex morphology and existing deformations in walls, columns and ledges, visible to the naked eye. This, together with its two-level porticated structure, makes the Casa de Pilatos suitable as a case study for the aim of this research. The Casa de Pilatos is the former Palace of the Adelantados of Andalusia, in the historic centre of Seville, Spain ( Figure 1). The building consists of spaces arranged and transformed over time, among which the main patio stands out. The combination of Mudejar, Gothic and Renaissance styles reveals the architectural richness of this space. Additionally, the patio has overcome diverse modifications due to the expansion needs as a Hispanic-Muslim playground [67]. As mentioned above, the existing pathologies are tilt in walls and columns and deflection in ledges, mainly due to the natural subsidence over time.
Remote Sens. 2020, 11, x FOR PEER REVIEW  4 of 23 format), thus inheriting their dimensional and texture properties. In addition, the accuracy and quality of the model's parameters are evaluated to verify their consistency with an appropriate scale. Regarding BIM operability with TLS or SFM data, the latest versions of Graphisoft ArchiCAD ® [23] and Autodesk Revit [34] allow to import these point clouds in TXT or XYZ file formats. This enables the BIM creation by fitting the morphological alterations under diverse modelling approaches. However, this paper advances on the automatic conversion of point cloud data into a parametric model by defining shapes through mathematical functions [51] under the Rhino+Grasshopper-ArchiCAD software combination. This workflow interconnects algorithms capable of modelling 3D objects so that the latter can be interoperable, which contributes to scientific knowledge in this field.

Case Study: The Main Patio Of The Casa De Pilatos
This research is based on a case study selected for its characteristics: a complex morphology and existing deformations in walls, columns and ledges, visible to the naked eye. This, together with its two-level porticated structure, makes the Casa de Pilatos suitable as a case study for the aim of this research. The Casa de Pilatos is the former Palace of the Adelantados of Andalusia, in the historic centre of Seville, Spain ( Figure 1). The building consists of spaces arranged and transformed over time, among which the main patio stands out. The combination of Mudejar, Gothic and Renaissance styles reveals the architectural richness of this space. Additionally, the patio has overcome diverse modifications due to the expansion needs as a Hispanic-Muslim playground [67]. As mentioned above, the existing pathologies are tilt in walls and columns and deflection in ledges, mainly due to the natural subsidence over time.

The Issue of the Scan-to-HBIM Approach
As mentioned above, BIM platforms allow to import point cloud data directly, so there is no need to use Rhinoceros® [62] at an earlier stage. The import is conducted using the interoperability tool for files with topographic coordinates, and then converting the point clouds into mesh objects. These text files must consist of rows of data containing three numerical entries (X, Y and Z coordinates). Next, the user has to create and align the historic building components referring to the point cloud. Only by using the floor/slab tool is the geometry automatically adapted into mesh-like parametric objects, leaving the option to modify the Z coordinate. In contrast, this procedure is only

The Issue of the Scan-to-HBIM Approach
As mentioned above, BIM platforms allow to import point cloud data directly, so there is no need to use Rhinoceros® [62] at an earlier stage. The import is conducted using the interoperability tool for files with topographic coordinates, and then converting the point clouds into mesh objects. These text files must consist of rows of data containing three numerical entries (X, Y and Z coordinates). Next, the user has to create and align the historic building components referring to the point cloud. Only Remote Sens. 2020, 12, 1094 5 of 22 by using the floor/slab tool is the geometry automatically adapted into mesh-like parametric objects, leaving the option to modify the Z coordinate. In contrast, this procedure is only possible using the wall tool if the point cloud is inserted in the model without any transformation. In this way, the user is able to identify existing deformations, but these are not implicit in the parametric 3D model. In addition, the mentioned parametric objects are included in BIM libraries as standard objects. In order to solve this issue, this paper addresses the automatic creation of parametric objects (as mesh, morph and library objects) from point cloud data under the combination of Rhinoceros, Grasshopper and ArchiCAD. The result is editable mesh parametric objects able to contain metadata and an identification label (ID). Concerning the morph objects, these elements have no geometric limits and are created so that the import of special shapes from other programmes is avoided.

Data Collection
The emergence of new digital technologies has revolutionised the documentation of architectural and cultural heritage. Their origin is the field of reverse engineering given the need to monitor serial manufacturing processes, thus controlling the dimensional quality of the product [68]. From the outset, the digital technology entailed a high economic cost, but the emergence of low-cost software made accurate techniques both affordable and suitable for creating geometrical heritage models.
In this research, diverse surveying methods were used to collect the data, such as classic surveying techniques-laser meter and flexometer to measure interior walls thicknesses-and TLS for the patio. The range and the accuracy of TLS make it ideal for measuring historic buildings [69]. Consequently, this technology was used in this research to carry out the two-stage 3D survey: i) the patio, where three stations (scanner positions) were set to create the XYZ global coordinate system, and ii) the galleries, with four stations in the corners by the galleries. Photographs were also taken using the laser scanner to produce a RGB point cloud ( Figure 2). The TLS device was a Leica ScanStation C10, with 120 m range and a 4-megapixel camera to colourise the cloud. The scan registration (scan alignment) was conducted using Leica Geosystems Cyclone®REGISTER 360 software [70]. The error alignment accuracy reached 3 mm in the patio and 5 mm in the galleries. The data used for this research was the TLS point cloud in .e57 format.
Remote Sens. 2020, 11, x FOR PEER REVIEW 5 of 23 possible using the wall tool if the point cloud is inserted in the model without any transformation. In this way, the user is able to identify existing deformations, but these are not implicit in the parametric 3D model. In addition, the mentioned parametric objects are included in BIM libraries as standard objects. In order to solve this issue, this paper addresses the automatic creation of parametric objects (as mesh, morph and library objects) from point cloud data under the combination of Rhinoceros, Grasshopper and ArchiCAD. The result is editable mesh parametric objects able to contain metadata and an identification label (ID). Concerning the morph objects, these elements have no geometric limits and are created so that the import of special shapes from other programmes is avoided.

Data Collection
The emergence of new digital technologies has revolutionised the documentation of architectural and cultural heritage. Their origin is the field of reverse engineering given the need to monitor serial manufacturing processes, thus controlling the dimensional quality of the product [68]. From the outset, the digital technology entailed a high economic cost, but the emergence of low-cost software made accurate techniques both affordable and suitable for creating geometrical heritage models.
In this research, diverse surveying methods were used to collect the data, such as classic surveying techniques-laser meter and flexometer to measure interior walls thicknesses-and TLS for the patio. The range and the accuracy of TLS make it ideal for measuring historic buildings [69]. Consequently, this technology was used in this research to carry out the two-stage 3D survey: i) the patio, where three stations (scanner positions) were set to create the XYZ global coordinate system, and ii) the galleries, with four stations in the corners by the galleries. Photographs were also taken using the laser scanner to produce a RGB point cloud ( Figure 2). The TLS device was a Leica ScanStation C10, with 120 m range and a 4-megapixel camera to colourise the cloud. The scan registration (scan alignment) was conducted using Leica Geosystems Cyclone® REGISTER 360 software [70]. The error alignment accuracy reached 3 mm in the patio and 5 mm in the galleries. The data used for this research was the TLS point cloud in .e57 format.

Point Cloud Meshing Workflow
The reconstruction algorithm is designed in such a way that an end user who does not need to know or understand the entire workflow should be able to obtain results from it with a low amount of effort. The workflow is between the traditional Scan-to-BIM method, which requires significant manual work, and completely automated AI-oriented methods that take the manual design control role from the architect or any other BIM operator.

Point Cloud Meshing Workflow
The reconstruction algorithm is designed in such a way that an end user who does not need to know or understand the entire workflow should be able to obtain results from it with a low amount of effort. The workflow is between the traditional Scan-to-BIM method, which requires significant manual work, and completely automated AI-oriented methods that take the manual design control role from the architect or any other BIM operator.
This workflow, using Grasshopper and Volvox [71] to process the point cloud and reconstruct surfaces, consists of 3 scripts ( Figure 3). The objects generated are later exported to ArchiCAD through the recent ArchiCAD+Grasshopper live connection.
Remote Sens. 2020, 11, x FOR PEER REVIEW 6 of 23 This workflow, using Grasshopper and Volvox [71] to process the point cloud and reconstruct surfaces, consists of 3 scripts ( Figure 3). The objects generated are later exported to ArchiCAD through the recent ArchiCAD+Grasshopper live connection. The interface requires user input to use cutting objects to remove segments from the cloud. After the subset is selected, further subselection based on proximity takes place. This operation decreases the complexity of the model and increases homogeneity, which is important for subsequent surface reconstruction.

Study Model with an Artificial Point Cloud
During the development of the algorithm, an artificial point cloud with adjustable parameters was created. Point clouds of a floor sector, two walls and a column were created to replicate scenarios such as the horizontal and vertical planes, and volumetric objects scanned, respectively. A splitting or clipping system is then established; the user creates cutting objects (plane/box/sphere) into the Rhino 3D space and the dynamic pipeline reads them in real-time ( Figure 4 and Figure 5).  The interface requires user input to use cutting objects to remove segments from the cloud. After the subset is selected, further subselection based on proximity takes place. This operation decreases the complexity of the model and increases homogeneity, which is important for subsequent surface reconstruction.

Study Model with an Artificial Point Cloud
During the development of the algorithm, an artificial point cloud with adjustable parameters was created. Point clouds of a floor sector, two walls and a column were created to replicate scenarios such as the horizontal and vertical planes, and volumetric objects scanned, respectively. A splitting or clipping system is then established; the user creates cutting objects (plane/box/sphere) into the Rhino 3D space and the dynamic pipeline reads them in real-time (Figures 4 and 5).
Remote Sens. 2020, 11, x FOR PEER REVIEW 6 of 23 This workflow, using Grasshopper and Volvox [71] to process the point cloud and reconstruct surfaces, consists of 3 scripts (Figure 3). The objects generated are later exported to ArchiCAD through the recent ArchiCAD+Grasshopper live connection. The interface requires user input to use cutting objects to remove segments from the cloud. After the subset is selected, further subselection based on proximity takes place. This operation decreases the complexity of the model and increases homogeneity, which is important for subsequent surface reconstruction.

Study Model with an Artificial Point Cloud
During the development of the algorithm, an artificial point cloud with adjustable parameters was created. Point clouds of a floor sector, two walls and a column were created to replicate scenarios such as the horizontal and vertical planes, and volumetric objects scanned, respectively. A splitting or clipping system is then established; the user creates cutting objects (plane/box/sphere) into the Rhino 3D space and the dynamic pipeline reads them in real-time ( Figure 4 and Figure 5).  Several cropping objects can be combined in such a way that if the objects intersect each other, the intersection point set is used ( Figure 6). Otherwise, the union of separately cropped point subsets is considered. After the selection is performed, further subsampling based on proximity in RTree [72] is performed to reduce the point quantity.
The algorithm then analyses the shape and planar deviation of the selected cloud subset. If the subset is within the planarity threshold specified by the operator/user/architect, it is then treated as a quasi-planar surface. The fitting plane is calculated and a Delaunay surface reconstruction method is used to produce a mesh object. The fitting planes vertical tilt angle is used to determine whether it is a vertical or horizontal structure ( Figure 7). Several cropping objects can be combined in such a way that if the objects intersect each other, the intersection point set is used ( Figure 6). Otherwise, the union of separately cropped point subsets is considered.
Remote Sens. 2020, 11, x FOR PEER REVIEW 7 of 23 Several cropping objects can be combined in such a way that if the objects intersect each other, the intersection point set is used ( Figure 6). Otherwise, the union of separately cropped point subsets is considered. After the selection is performed, further subsampling based on proximity in RTree [72] is performed to reduce the point quantity.
The algorithm then analyses the shape and planar deviation of the selected cloud subset. If the subset is within the planarity threshold specified by the operator/user/architect, it is then treated as a quasi-planar surface. The fitting plane is calculated and a Delaunay surface reconstruction method is used to produce a mesh object. The fitting planes vertical tilt angle is used to determine whether it is a vertical or horizontal structure (Figure 7). After the selection is performed, further subsampling based on proximity in RTree [72] is performed to reduce the point quantity.
The algorithm then analyses the shape and planar deviation of the selected cloud subset. If the subset is within the planarity threshold specified by the operator/user/architect, it is then treated as a quasi-planar surface. The fitting plane is calculated and a Delaunay surface reconstruction method is used to produce a mesh object. The fitting planes vertical tilt angle is used to determine whether it is a vertical or horizontal structure (Figure 7). Remote Sens. 2020, 11, x FOR PEER REVIEW 8 of 23 On the other hand, in case the selected subset is a volumetric collection of points, a ball-pivot algorithm and several post-processing mesh cleaning procedures are performed to produce 3D shapes ( Figure 8).
Post-processing of the mesh includes: removing faces adjacent to non-manifold edges, removing duplicate faces, correcting naked edges, filling the holes (as much as possible) and rebuilding normal vectors.

Processing Delaunay Meshes
As Delaunay meshing always outputs convex geometry, post-processing is performed to cull the extra parts which most of the time are distinguished by having longer edges. The threshold input option is given to the user to control the culling. This is particularly helpful in openings ( Figure 9).  On the other hand, in case the selected subset is a volumetric collection of points, a ball-pivot algorithm and several post-processing mesh cleaning procedures are performed to produce 3D shapes ( Figure 8).
Remote Sens. 2020, 11, x FOR PEER REVIEW 8 of 23 On the other hand, in case the selected subset is a volumetric collection of points, a ball-pivot algorithm and several post-processing mesh cleaning procedures are performed to produce 3D shapes ( Figure 8).
Post-processing of the mesh includes: removing faces adjacent to non-manifold edges, removing duplicate faces, correcting naked edges, filling the holes (as much as possible) and rebuilding normal vectors.

Processing Delaunay Meshes
As Delaunay meshing always outputs convex geometry, post-processing is performed to cull the extra parts which most of the time are distinguished by having longer edges. The threshold input option is given to the user to control the culling. This is particularly helpful in openings ( Figure 9).  Post-processing of the mesh includes: removing faces adjacent to non-manifold edges, removing duplicate faces, correcting naked edges, filling the holes (as much as possible) and rebuilding normal vectors.

Processing Delaunay Meshes
As Delaunay meshing always outputs convex geometry, post-processing is performed to cull the extra parts which most of the time are distinguished by having longer edges. The threshold input option is given to the user to control the culling. This is particularly helpful in openings ( Figure 9). On the other hand, in case the selected subset is a volumetric collection of points, a ball-pivot algorithm and several post-processing mesh cleaning procedures are performed to produce 3D shapes ( Figure 8).
Post-processing of the mesh includes: removing faces adjacent to non-manifold edges, removing duplicate faces, correcting naked edges, filling the holes (as much as possible) and rebuilding normal vectors.

Processing Delaunay Meshes
As Delaunay meshing always outputs convex geometry, post-processing is performed to cull the extra parts which most of the time are distinguished by having longer edges. The threshold input option is given to the user to control the culling. This is particularly helpful in openings ( Figure 9).  Additionally, the model in the galleries will consist of point cloud partitions at the same level representing ceilings in each area. Here, their continuity is only maintained through the walls. If sub selected at a certain level (Figure 10), the meshes can have all the intermediary Delaunay edges deleted in one meshing process, i.e., the user can directly and individually generate the parts.
Remote Sens. 2020, 11, x FOR PEER REVIEW 9 of 23 Additionally, the model in the galleries will consist of point cloud partitions at the same level representing ceilings in each area. Here, their continuity is only maintained through the walls. If sub selected at a certain level (Figure 10), the meshes can have all the intermediary Delaunay edges deleted in one meshing process, i.e., the user can directly and individually generate the parts.

Crease Processing
As the user segments the point cloud using objects (e.g., boxes), it is extremely complicated to obtain the points that are specifically in the intersections of surfaces such as walls or ceilings. This is especially common in heritage buildings, since their structure can be deformed, and the resulting point cloud data samples may be dispersed in the crease regions. In order to overcome this issue, the architect or BIM operator is advised to use overlapping selections. This will create meshes as presented in the upper half of ( Figure 11). A separate part of the algorithm is developed to process the crease region between two adjacent meshes so that there are no shared or overlapping faces. An oriented bounding box of the crease is created to establish the working space. The faces that fall into this bounding box from both meshes are culled. Based on the points inside the working space a new mesh is generated again using the Delaunay triangulation with the fitted plane. The crease, based on points from the cloud, is then created. The face normals are checked against the average

Crease Processing
As the user segments the point cloud using objects (e.g., boxes), it is extremely complicated to obtain the points that are specifically in the intersections of surfaces such as walls or ceilings. This is especially common in heritage buildings, since their structure can be deformed, and the resulting point cloud data samples may be dispersed in the crease regions. In order to overcome this issue, the architect or BIM operator is advised to use overlapping selections. This will create meshes as presented in the upper half of ( Figure 11). A separate part of the algorithm is developed to process the crease region between two adjacent meshes so that there are no shared or overlapping faces.
Remote Sens. 2020, 11, x FOR PEER REVIEW 9 of 23 Additionally, the model in the galleries will consist of point cloud partitions at the same level representing ceilings in each area. Here, their continuity is only maintained through the walls. If sub selected at a certain level (Figure 10), the meshes can have all the intermediary Delaunay edges deleted in one meshing process, i.e., the user can directly and individually generate the parts.

Crease Processing
As the user segments the point cloud using objects (e.g., boxes), it is extremely complicated to obtain the points that are specifically in the intersections of surfaces such as walls or ceilings. This is especially common in heritage buildings, since their structure can be deformed, and the resulting point cloud data samples may be dispersed in the crease regions. In order to overcome this issue, the architect or BIM operator is advised to use overlapping selections. This will create meshes as presented in the upper half of ( Figure 11). A separate part of the algorithm is developed to process the crease region between two adjacent meshes so that there are no shared or overlapping faces. An oriented bounding box of the crease is created to establish the working space. The faces that fall into this bounding box from both meshes are culled. Based on the points inside the working space a new mesh is generated again using the Delaunay triangulation with the fitted plane. The crease, based on points from the cloud, is then created. The face normals are checked against the average An oriented bounding box of the crease is created to establish the working space. The faces that fall into this bounding box from both meshes are culled. Based on the points inside the working space a new mesh is generated again using the Delaunay triangulation with the fitted plane. The crease, based on points from the cloud, is then created. The face normals are checked against the average normals of the larger meshes to determine which face belongs to which mesh. In the end, as shown in the lower half of Figure 11, a clean crease between the two meshes is produced.
In order to aid this process, collision events are tested for all potential pairs and a connectivity graph is drawn onto the 3D model ( Figure 12). This shows the operator-architect which pairs of meshes should be checked for crease processing.
Remote Sens. 2020, 11, x FOR PEER REVIEW 10 of 23 normals of the larger meshes to determine which face belongs to which mesh. In the end, as shown in the lower half of Figure 11, a clean crease between the two meshes is produced.
In order to aid this process, collision events are tested for all potential pairs and a connectivity graph is drawn onto the 3D model ( Figure 12). This shows the operator-architect which pairs of meshes should be checked for crease processing.

Interface
Using the Human UI for Grasshopper [73] a custom interface was developed for the architectoperator to make use of the algorithm. The interface ( Figure 13) is used for both input values that will determine how the algorithm works and output values that are useful for the user. This includes controls for point cloud subsampling, proximity control for point culling, deviation to treat the selection as planar or volumetric, among others. Some parts of the interface are updated based on the subselection of points. If the subset is determined as planar, an input for Delaunay mesh edge average length multiplier is given to control how much of the long edge faces will be culled. If the subset is determined as volumetric, the input is given for the ball-pivot radius multiplier. Optional smoothing of the mesh may be also performed, which will not increase the mesh face or vertex count.

Interface
Using the Human UI for Grasshopper [73] a custom interface was developed for the architect-operator to make use of the algorithm. The interface ( Figure 13) is used for both input values that will determine how the algorithm works and output values that are useful for the user. This includes controls for point cloud subsampling, proximity control for point culling, deviation to treat the selection as planar or volumetric, among others. Some parts of the interface are updated based on the subselection of points. If the subset is determined as planar, an input for Delaunay mesh edge average length multiplier is given to control how much of the long edge faces will be culled. If the subset is determined as volumetric, the input is given for the ball-pivot radius multiplier. Optional smoothing of the mesh may be also performed, which will not increase the mesh face or vertex count.

Automatic Creation Textures
Point cloud data can also contain information about the colour of the objects scanned. This data, included in the points, is also transferred through the algorithm. However, the ArchiCAD-Grasshopper live connection does not allow the transfer of colour or texture information with morphs so far. This is because texturing in ArchiCAD is not carried out by manual colour coding of individual faces of the morph object but by mapping a texture file onto the object. The latter requires having an actual image file of the texture.
As building components often contain flat-like surfaces (walls, slabs, etc.), a texture image can be used to map them with colour information. The study of the texture creation was conducted on scan data of a wall sector from the Casa de Pilatos Palace. In order to translate the RGB colour information from the point cloud to a texture image file, the following procedures have been implemented: (1) a fitted plane is calculated from the point cloud; (2) the points are projected onto the plane; (3) given that the points are located on the same plane, a bounding rectangle is found to act as the canvas of the future texture image; and (4) the Voronoi [56] division method is used to divide the space into smaller cells, since the point distribution is not regular-diverse point densities exist. Each cell represents an area to be filled with the colour driven by the point which created that cell ( Figure 14).

Automatic Creation Textures
Point cloud data can also contain information about the colour of the objects scanned. This data, included in the points, is also transferred through the algorithm. However, the ArchiCAD-Grasshopper live connection does not allow the transfer of colour or texture information with morphs so far. This is because texturing in ArchiCAD is not carried out by manual colour coding of individual faces of the morph object but by mapping a texture file onto the object. The latter requires having an actual image file of the texture.
As building components often contain flat-like surfaces (walls, slabs, etc.), a texture image can be used to map them with colour information. The study of the texture creation was conducted on scan data of a wall sector from the Casa de Pilatos Palace. In order to translate the RGB colour information from the point cloud to a texture image file, the following procedures have been implemented: (1) a fitted plane is calculated from the point cloud; (2) the points are projected onto the plane; (3) given that the points are located on the same plane, a bounding rectangle is found to act as the canvas of the future texture image; and (4) the Voronoi [56] division method is used to divide the space into smaller cells, since the point distribution is not regular-diverse point densities exist. Each cell represents an area to be filled with the colour driven by the point which created that cell ( Figure 14).

Transfer to BIM (ArchiCAD)
The Grasshopper+ArchiCAD live connection allows the bilateral transfer of geometry between them. Since ArchiCAD works in a BIM environment, it is more rigid when representing nonrectilinear complex geometries. As such the meshes generated with the algorithm contain many faces. There are 3 main object types (tools) in ArchiCAD that can deal with complex forms: 1) the Mesh tool operates the Delaunay meshing algorithm, but it only works in XY plane. This implies that only the flat horizontal parts can be represented using this tool, which may be useful for floors and ceilings; 2) the Morph tool is more flexible and corresponds to the mesh representation in any other usual 3D modelling software, however it is slower than the ArchiCAD Mesh tool; and 3) the Library object tool is also flexible in geometry creation. The library objects are GDL entities containing geometry and texture. Notwithstanding, the colour of the vertexes is not preserved through the transfer of mesh geometry to ArchiCAD. One way of overcoming this issue could be treating each individual mesh face as a separate library part entity, as in the work by Pavlo Menshykh [74], but this approach is not suitable for massive object meshes. The result of the algorithm implementation in BIM is shown in Figure 15.

Transfer to BIM (ArchiCAD)
The Grasshopper+ArchiCAD live connection allows the bilateral transfer of geometry between them. Since ArchiCAD works in a BIM environment, it is more rigid when representing non-rectilinear complex geometries. As such the meshes generated with the algorithm contain many faces. There are 3 main object types (tools) in ArchiCAD that can deal with complex forms: 1) the Mesh tool operates the Delaunay meshing algorithm, but it only works in XY plane. This implies that only the flat horizontal parts can be represented using this tool, which may be useful for floors and ceilings; 2) the Morph tool is more flexible and corresponds to the mesh representation in any other usual 3D modelling software, however it is slower than the ArchiCAD Mesh tool; and 3) the Library object tool is also flexible in geometry creation. The library objects are GDL entities containing geometry and texture. Notwithstanding, the colour of the vertexes is not preserved through the transfer of mesh geometry to ArchiCAD. One way of overcoming this issue could be treating each individual mesh face as a separate library part entity, as in the work by Pavlo Menshykh [74], but this approach is not suitable for massive object meshes. The result of the algorithm implementation in BIM is shown in Figure 15.

Transfer to BIM (ArchiCAD)
The Grasshopper+ArchiCAD live connection allows the bilateral transfer of geometry between them. Since ArchiCAD works in a BIM environment, it is more rigid when representing nonrectilinear complex geometries. As such the meshes generated with the algorithm contain many faces. There are 3 main object types (tools) in ArchiCAD that can deal with complex forms: 1) the Mesh tool operates the Delaunay meshing algorithm, but it only works in XY plane. This implies that only the flat horizontal parts can be represented using this tool, which may be useful for floors and ceilings; 2) the Morph tool is more flexible and corresponds to the mesh representation in any other usual 3D modelling software, however it is slower than the ArchiCAD Mesh tool; and 3) the Library object tool is also flexible in geometry creation. The library objects are GDL entities containing geometry and texture. Notwithstanding, the colour of the vertexes is not preserved through the transfer of mesh geometry to ArchiCAD. One way of overcoming this issue could be treating each individual mesh face as a separate library part entity, as in the work by Pavlo Menshykh [74], but this approach is not suitable for massive object meshes. The result of the algorithm implementation in BIM is shown in Figure 15.

Assigning Semantic Data to Objects
The BIM environment is distinguished among 3D modelling approaches by the ability to assign rich semantic data for each individual object. This is useful in the AEC (Architecture-Engineering-Construction) industry practices. Each object contains information about its construction element classification, object identifier, renovation cycle type, structural loads, etc. A similar function is available in Rhinoceros environment, where each object is able to include associated user-defined Key+Value attributes. These attributes can be populated and have assigned values to them. The information contained in the Rhinoceros object as a user attribute can be later translated to an ArchiCAD property. In the ArchiCAD environment, however, the semantic data are not simple; some are inherent to the specific object type, others are general and some (after the ArchiCAD 20 update) can be customised by the user in order to represent a certain kind of object associated phenomena. To establish a workflow of the translation of semantic data from Rhinoceros attributes to ArchiCAD data, a script was created to read the custom made ID field of the Rhinoceros object and use it as an input for the object attributes when creating the object (Morph/Library) in ArchiCAD. This approach succeeds to transfer the semantic data; therefore, the other parameters can be set to be transferred in the same manner.

Texture Creation
The aforementioned automatic method enables texture image creation regardless of the point cloud composition. As can be seen in the given example (previous Figures 14, 16 and 17), a portion of the point cloud was subtracted due to the bench in front of the wall. This led to a discontinuity of the cloud particles in that place. Nevertheless, the final image is complete.
Remote Sens. 2020, 11, x FOR PEER REVIEW 13 of 23 The BIM environment is distinguished among 3D modelling approaches by the ability to assign rich semantic data for each individual object. This is useful in the AEC (Architecture-Engineering-Construction) industry practices. Each object contains information about its construction element classification, object identifier, renovation cycle type, structural loads, etc. A similar function is available in Rhinoceros environment, where each object is able to include associated user-defined Key+Value attributes. These attributes can be populated and have assigned values to them. The information contained in the Rhinoceros object as a user attribute can be later translated to an ArchiCAD property. In the ArchiCAD environment, however, the semantic data are not simple; some are inherent to the specific object type, others are general and some (after the ArchiCAD 20 update) can be customised by the user in order to represent a certain kind of object associated phenomena. To establish a workflow of the translation of semantic data from Rhinoceros attributes to ArchiCAD data, a script was created to read the custom made ID field of the Rhinoceros object and use it as an input for the object attributes when creating the object (Morph/Library) in ArchiCAD. This approach succeeds to transfer the semantic data; therefore, the other parameters can be set to be transferred in the same manner.

Texture Creation
The aforementioned automatic method enables texture image creation regardless of the point cloud composition. As can be seen in the given example (previous Figure 14, Figure 16 and Figure  17), a portion of the point cloud was subtracted due to the bench in front of the wall. This led to a discontinuity of the cloud particles in that place. Nevertheless, the final image is complete.
The population size sampled for texture creation can be adjusted to obtain low to high resolution images. However, the resolution here should not be misinterpreted with the pixel resolution of the image. The pixel resolution is controlled separately using an input parameter that stands for the height pixel amount. The length derives from the ratio of the bounding rectangle.  The BIM environment is distinguished among 3D modelling approaches by the ability to assign rich semantic data for each individual object. This is useful in the AEC (Architecture-Engineering-Construction) industry practices. Each object contains information about its construction element classification, object identifier, renovation cycle type, structural loads, etc. A similar function is available in Rhinoceros environment, where each object is able to include associated user-defined Key+Value attributes. These attributes can be populated and have assigned values to them. The information contained in the Rhinoceros object as a user attribute can be later translated to an ArchiCAD property. In the ArchiCAD environment, however, the semantic data are not simple; some are inherent to the specific object type, others are general and some (after the ArchiCAD 20 update) can be customised by the user in order to represent a certain kind of object associated phenomena. To establish a workflow of the translation of semantic data from Rhinoceros attributes to ArchiCAD data, a script was created to read the custom made ID field of the Rhinoceros object and use it as an input for the object attributes when creating the object (Morph/Library) in ArchiCAD. This approach succeeds to transfer the semantic data; therefore, the other parameters can be set to be transferred in the same manner.

Texture Creation
The aforementioned automatic method enables texture image creation regardless of the point cloud composition. As can be seen in the given example (previous Figure 14, Figure 16 and Figure  17), a portion of the point cloud was subtracted due to the bench in front of the wall. This led to a discontinuity of the cloud particles in that place. Nevertheless, the final image is complete.
The population size sampled for texture creation can be adjusted to obtain low to high resolution images. However, the resolution here should not be misinterpreted with the pixel resolution of the image. The pixel resolution is controlled separately using an input parameter that stands for the height pixel amount. The length derives from the ratio of the bounding rectangle.  The population size sampled for texture creation can be adjusted to obtain low to high resolution images. However, the resolution here should not be misinterpreted with the pixel resolution of the image. The pixel resolution is controlled separately using an input parameter that stands for the height pixel amount. The length derives from the ratio of the bounding rectangle.

Accuracy Evaluation
Once the meshing algorithm produced the geometry, it was worth assessing the meshing accuracy against the original subsampled point cloud. This was carried out in two sectors of the patio. First, the 3D meshing of a wall sector, whose point cloud had a total of 24,399 points, and its mesh had 3943 polygons and 2033 vertexes. Second, using a sample small mullion column in the patio, with a point cloud consisting of 11,493 points, and a mesh with 1704 faces and 854 points. The geometry comparison was performed as per [8,9], using CloudCompare software to compute the cloud-to-mesh distances (Figures 18 and 19).

Accuracy Evaluation
Once the meshing algorithm produced the geometry, it was worth assessing the meshing accuracy against the original subsampled point cloud. This was carried out in two sectors of the patio. First, the 3D meshing of a wall sector, whose point cloud had a total of 24,399 points, and its mesh had 3943 polygons and 2033 vertexes. Second, using a sample small mullion column in the patio, with a point cloud consisting of 11,493 points, and a mesh with 1704 faces and 854 points. The geometry comparison was performed as per [8,9], using CloudCompare software to compute the cloud-to-mesh distances ( Figure 18 and Figure 19).  The scale in both histograms is constant to ease comprehension. The abscissa axis represents the distance intervals between the two 3D objects, and the ordinate axis shows the number of points in each interval. The 'zero' deviation value on the X-axis represents the correspondence between the objects compared-3D mesh and point cloud [8,9]. In other words, it is the interval in which the distance between the mesh and the point cloud is null. A total of 2000 distance intervals were created for the calculation; the higher the number of intervals is, the more accurate the number of points in the 'zero' value. The standard deviation (σ) was the statistical tool considered to assess the accuracy of the meshing process. In this paper, σ was 68.28 points in the wall sector, and 28.13 vertexes in the sample mullion column.

Discussion of Results
Firstly, the choice of programmes and algorithm structure should be discussed. Given that 3D modelling applications such as Rhinoceros and ArchiCAD are not intended for point cloud processing, the operations performed in this research may be slower than if performed in specific point cloud processing software. However, the live connection explored in this research (Rhinoceros+Grasshopper-ArchiCAD) allows for linking visual programming to BIM in order to create BIM entities in a complex case study, including flat-like and volumetric building components. Attempts to make large amount of calculations at once may lead to crashes; therefore, in order to achieve a more effective procedure, the algorithm is designed to work on a step-by-step basis through point cloud subselection and processing. As a result, considering the aforementioned software limitations, processing time is reduced, thus revealing the established workflow as an agreement between speed and efficiency.
Secondly, it is worth discussing the different levels of parametrisation. BIM software normally identify their tools and objects as being parametric, however that level of parametrisation is normally bound to those specific objects. For instance, a table can have the drawer sizes parametrically adjusted based on the tabletop dimensions. This is a parametric model but limited to the inherent parameters of the given object. On the other hand, a more global parametrisation suggests that there is a parameter-based connection not only inside one entity but also between different objects. In this regard, the methodology presented here achieves various levels of parametrisation. Examples of object level, multi-object level and global level parametrisations are given below: • Surface reconstruction, which performs several calculations in order to obtain different parameters for a selected subset of the point cloud.

•
The detection of colliding pairs for crease processing described in section 2.6 is a multi-object parametrisation.

•
Global parameters: to determine the resolution of point clouds or the BIM objects type choice that the reconstructed mesh will be translated to.
With a view to avoid confusion with the terms parametric (model) and parametrisation, the BIM entities in this paper (mesh/morph/GDL objects) should not be considered as purely parametric, but as elements which depend on parameters at a certain level.
Thirdly, the performance of the workflow and the user interaction are discussed. The performance depends on the meshing method. Delaunay algorithm in Grasshopper is generally faster and more reliable than the ball-pivot method for point cloud 3D meshing, since the latter was not completely stable in this research. This algorithm caused the programme to crash occasionally, which rarely happened when cropped subsampled clouds contain less than 5000 points. The workflow was tested on an HP Pavilion gaming laptop for performance.
The performance obviously depends on the resolution the user is working with. The higher it is, the smaller portions are advised to be taken for meshing at a time. In this case study, the average collection of points selected for meshing after tolerance culling ranged from 2000 to 10,000 points. In this case, the workflow experienced negligible latency and the calculations took less than 0.5 seconds to compute. The only process that would make the user wait for a short lapse could be the initial point cloud down sampling. On the contrary, transferring the generated meshes into ArchiCAD takes longer time (several seconds for a mesh with 2000 to 5000 points). The more objects that are present in ArchiCAD through the live connection, the longer it will take at each subsequent step. The above led to create a different script for the transfer stage only.
Next, this section is divided into three parts: user recommendations, methods for surface reconstruction and accuracy assessment of the research outcomes.

User Recommendations
The algorithm has been tested in this research; several recommendations can therefore be provided to ease the workflow:

•
Every time the geometry is created, it is not necessary to send it to ArchiCAD in order to save the data. The meshes are already saved in Rhinoceros, which means that the transfer to ArchiCAD can be conducted every 5-6 steps or more.

•
If working with architectural models, it is highly advisable to position the point cloud in such a way that the least amount of subsampling operations is required. Moreover, it is advisable to orient the cloud to the Cartesian system in Rhinoceros.

•
It is easier to control the cropping objects in planar views (top, front, side) rather than in perspective.

•
In certain cases, as in this case study, there may be long continuous surfaces such as walls. In the case of heritage objects, the deviation from the fitted plane of those long continuous elements may be extremely high if taken as one piece. In that case, the object may be reconstructed from several pieces that can be determined by the user. This will ensure a better and faster workflow, especially concerning the crease processing algorithm, resulting in seamless and smoother component joints.

Different Methods of Source Reconstruction
It is worth discussing the effectiveness of the aforementioned methods of reconstructing surfaces [59][60][61][62][63][64][65][66][67]. Certain recently developed algorithms are capable of reconstructing surface efficiently and with high quality of the output mesh. However, there is no implementation of such techniques into the Grasshopper visual scripting environment. Creating Grasshopper components of those techniques is out of the scope of this research; thus, only the methods available at hand were used. Some mesh post-processing operations were developed to test the following set of methods. The algorithm is working only with point positions, so the surface reconstruction can be performed from unoriented point clouds. A point cloud subset was selected to test different methods available for meshing: the fountain in the centre of the patio.

Delaunay Meshing
The Delaunay meshing method is available within Grasshopper and is probably the fastest. It is capable of processing large quantities of points within reasonable time frames. The flat-like parts of the point clouds are reconstructed via Delaunay triangulation with some post-processing.

Ball-Pivot
The ball-pivot algorithm [75] is one of such available components used in this paper. It has several advantages and drawbacks. As the radius of the pivot ball remains constant, the accuracy of the yielded model highly depends on the uniformity of the point cloud. If the point cloud has varying densities, which is generally true, the algorithm will fail to connect more distant points and will overrun the closer ones. For this reason, a proximity culling of the points is performed prior to running the ball-pivot algorithm on the point cloud. Although this will not increase the density in the undersampled regions, it will decrease the oversampled region density, thus bringing the overall densities closer to the mean. In order to determine the ball-pivot radius, random subsampling is performed to test the average distance to the nearest 3 neighbouring points. The average of the values is then taken as the base radius. The user is given a numeric control of a multiplier of this value to refine the radius in the end.

Marching Cubes
Marching Cubes [76] is another surface reconstruction algorithm available within Grasshopper via the Cocoon developed by David Stasiuk. The resulting reconstructed surfaces were not just thin meshes but more like volumetric envelopes of the objects. In addition, the Marching Cubes algorithm required more parameter inputs from the user, which is far from being automatic.

Voxelisation
Volvox includes the voxelisation meshing algorithm. This method was also tested during this research. Although the resulting mesh highly resembled the initial object from which the points were scanned, this method is inapplicable to heritage reconstruction since the final mesh does not represent the scanned points accurately, and may lack crucial data.
Finally, Delaunay meshing was chosen to reconstruct the surfaces in flat regions and ball-pivot was used for the volumetric ones.

Meshing Accuracy
In view of the results of the meshing accuracy evaluation, it is worth noting that, according to Antón et al. [8], high values of the standard deviation account for high accuracy in the 3D meshing process. The values obtained in this research (68.28 and 28.13 vertexes) are lower than the value considered optimal by those authors; i.e., the accuracy is also lower. The extreme distance values between the 3D objects are 0.04 m and 0.057 m, respectively. However, these values may be outliers, since the starting populated intervals are those below approximately 0.02 m and 0.04 m, respectively. It is also worth highlighting the symmetry in the data deriving from the accuracy assessment of the wall sector, whereas the case of the mullion column shows a more irregular distribution (Figures 18  and 19). This could indicate a lower 3D meshing accuracy of the algorithm in volumetric components. Consequently, all the above are limitations of this work, but also an opportunity to improve the capabilities of the meshing algorithm in the future. Nevertheless, it should be noted that the aim of the work by the aforementioned authors was to obtain the highest possible accuracy in their 3D heritage models for HBIM by weighing the geometry generation and the mesh triangle amount, which has a direct impact on processing time, computational resources and file size. Thus, the standard deviation calculated in this paper can be accepted for the purpose of modelling the patio of Casa de Pilatos in HBIM. Certainly, the model produced does not include complex decorative details yet, but the automatic algorithm is able to represent the overall shape, volume and deformations of the historical building components under the HBIM environment.

Conclusions
The scientific advance of this work consists of automating and validating the use of BIM in heritage buildings with complex structures. Thus, researchers, architects and BIM operators can identify geometric properties within the platform and represent, evaluate and inventory both structural and pathological deformations of the assets.
In this paper, a new workflow is created to transfer point cloud data into HBIM using the Rhinoceros+Grasshopper-ArchiCAD combination, thus avoiding complementary software such as Cloudcompare, MeshLab or Geomagic. Here, starting from TLS and SFM data, a pseudo point cloud structure with controllable density is created so that the point clouds are converted into textured meshes.
Although according to Facundo et al. [77] Autodesk Revit currently leads the ranking of use in BIM publications in heritage, this software does not solve the problems of surface modelling [52] by itself. In this sense, the Rhino+Grasshopper-ArchiCAD combination demonstrates the possibilities of a new working scenario based on interoperability.
There is a huge potential of parametric modelling, which was described by Lee et al. [78] in 2006. HBIM platforms need algorithmic designs to represent complex structures, and workflows like the one developed in this research can become the next step in this knowledge area. It is in the hands of BIM software developers to introduce these workflows to increase operability as the capabilities of computers to manipulate massive point cloud data improve.
In this work, diverse algorithms such as Marching Cubes, Volvox, Delaunay and ball-pivot have been tested to prove their effectiveness in modelling objects containing geometrical alterations. This is an experimental work to achieve the conversion of TLS or SFM point cloud data into HBIM, thus constituting semantically enriched parametric objects. The 3D meshing accuracy assessment of the algorithm reveals its suitability to represent the overall features and geometrical alterations of the elements modelled.
As mentioned above, not many surface reconstruction algorithms are available at hand in graphical algorithm editors by the time this work was conducted. In the case of the advent of better surface reconstruction tools, the workflow may adopt those new methods for better performance. The ArchiCAD-Grasshopper live connection itself is a relatively new development; further functionalities are therefore expected to be developed in the near future, which will improve the proposed workflow. This may include the possibility of transferring textures directly into the objects, saving generated objects in ArchiCAD in a dynamic manner, among others. This workflow will be used in future case studies and more refinement will be carried out to improve the algorithmic performance.