extfieldsofstudy:
Computer Science
Physics
Medicine

text:
Efficient partition of integer optimization problems with one-hot encoding

Quantum annealing is a heuristic algorithm for solving combinatorial optimization problems, and hardware for implementing this algorithm has been developed by D-Wave Systems Inc. The current version of the D-Wave quantum annealer can solve unconstrained binary optimization problems with a limited number of binary variables. However, the cost functions of several practical problems are defined by a large number of integer variables. To solve these problems using the quantum annealer, integer variables are generally binarized with one-hot encoding, and the binarized problem is partitioned into small subproblems. However, the entire search space of the binarized problem is considerably larger than that of the original integer problem and is dominated by infeasible solutions. Therefore, to efficiently solve large optimization problems with one-hot encoding, partitioning methods that extract subproblems with as many feasible solutions as possible are required. In this study, we propose two partitioning methods and demonstrate that they result in improved solutions.

The combinatorial optimization problems aim to minimize cost functions defined by discrete variables, and these problems often have significant real-world applications. In general, the cost function of a combinatorial optimization problem can be expressed as the Hamiltonian of a classical Ising model 1 . Therefore, many algorithms for solving combinatorial optimization problems have been inspired by physics. Simulated annealing (SA) 2 is one of the most famous algorithms, employing thermal fluctuations to escape local minima. In contrast to SA, quantum annealing (QA) 3 is a method that exploits quantum fluctuations and the resulting tunneling effect. A popular research topic involves evaluating whether it is more advantageous to employ quantum effects than thermal fluctuations, and numerous studies have been conducted on this topic [4][5][6][7][8][9] . In addition, the recent development of a commercial quantum annealer by D-Wave Systems Inc 10 . has attracted many companies and researchers. The performance of QA has been experimentally studied using the quantum annealer and compared with that of SA [11][12][13] , and several companies have demonstrated the applicability of the annealer to practical problems [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] .
The time-dependent Hamiltonian of QA is given as followŝˆ= where H 0 is the target Hamiltonian representing the cost function, and H q denotes the quantum fluctuation term for which the ground state is trivial. The initial values of the coefficients are set to A(0) = 1 and B(0) = 0, and the system is prepared in the trivial ground state determined by Ĥ q . Then, the strength of the quantum fluctuation is reduced toward zero, and the coefficients are set to A(τ) = 0 and B(τ) = 1 at the end of QA, where τ is the annealing time. The dynamics of the system is described by the Schrödinger equation: where ψ(t) is the state vector of the system, and ℏ is set to 1 for simplicity. Given that the coefficients change sufficiently slowly, the adiabatic theorem 31 ensures that the system remains close to the instantaneous ground state of the time-dependent Hamiltonian. Thus, by setting the annealing time τ to be sufficiently large, the ground states of the target Hamiltonian H 0 can be obtained with high probability. The current version of the D-Wave quantum annealer (D-Wave 2000Q) implements transverse-magnetic-field QA, for which the quantum fluctuation is given as follows:Ĥ , (3) where N q denotes the number of qubits. The quantum annealer can handle a cost function as follows: where the interactions between qubits are restricted to the Chimera graph 32 , which is a 16 × 16 grid of complete bipartite graphs K 4,4 in D-Wave 2000Q. It should be noted that the number of operable qubits is less than N q = 2,048 due to defects in the qubits and connectivities. Due to the limited number of available qubits, large optimization problems cannot be solved directly using the D-Wave quantum annealer. In real settings, large problems are partitioned into subproblems that can be handled by the quantum annealer. The subproblems are iteratively optimized by the quantum annealer, and the optimization result is used to improve the current solution [33][34][35] . A cluster of spins in the subproblem is simultaneously updated in this scheme; this iterative method is a type of large-neighborhood local search algorithm 36 . Although these algorithms can be performed using classical computers, subproblems are fundamentally restricted to tree structures that are solvable in polynomial time by belief propagation or dynamic programming [37][38][39][40] . Therefore, using the quantum annealer is advantageous if it can solve subproblems with many closed loops more efficiently than classical algorithms. Furthermore, solving subproblems that are as large as possible is essential for improving solution accuracy 41 . The size of subproblems that can be embedded into the quantum annealer strongly depends on the quality of the minor embedding, in particular for problems with few interactions. Because subproblems must be iteratively embedded, fast algorithms for embedding larger subproblems are required for exploiting the potential of the quantum annealer. Although complete-graph embedding [42][43][44] can be used for problems with dense interactions, a subproblem-embedding algorithm, that was developed in a previous study 41 , may be effective in improving the solution accuracy of sparse problems.
In addition, the quantum annealer requires the cost function to be represented in the form of a quadratic unconstrained binary optimization (QUBO) problem or Ising model; however, many cost functions in practical problems are defined by integer variables. The binarization of integer variables is generally achieved using one-hot encoding 1 . For example, suppose that we wish to solve the following integer optimization problem with N integer variables {S i } i = 1,2,...,N : is a binary variable that is assigned to component q of S i , and x 1 i q ( ) = indicates that component q is selected for S i . In addition, feasible solutions are constrained to configurations in which exactly one component is selected for each S i . Subsequently, a penalty term is introduced to obtain the following unconstrained form: where the second term formulates the penalty term which is introduced to extract feasible solutions satisfying the , which we call the one-hot constraint, and parameter λ controls the strength of the penalty term. By setting parameter λ to a sufficiently large value, the ground states of the original integer optimization problem (Eq. (5)) are correctly encoded. However, the performance of the D-Wave quantum annealer is significantly affected by noise and intrinsic control errors when λ is larger than necessary. Therefore, to obtain highly accurate solutions, we must explore an appropriate value of λ, which is a tedious task for optimization under the one-hot constraint. In addition, the entire search space of the binarized optimization problem (Eq. (7)) is dominated by infeasible solutions. Figure 1(a) presents the problem graph of Eq. (7). In this figure, vertices and edges represent binary variables and the interactions between them, respectively. Q binary variables 1,2, , are assigned to each S i , and the total number of binary variables is NQ. Although the number of configurations of the binary variables is 2 NQ , the number of feasible solutions is only Q N . Therefore, to efficiently solve large optimization problems under the one-hot constraint using the quantum annealer, partitioning methods are required for extracting subproblems with as many feasible solutions as possible. A simple example of an undesirable partition is depicted in Fig. 1(b). Here, suppose that we wish to improve the current solution presented in Fig. 1(b) and that the three binary variables enclosed by the green rectangle are extracted as the subproblem. In this case, superior www.nature.com/scientificreports www.nature.com/scientificreports/ feasible solutions cannot be explored by optimizing the subproblem because only the current solution in the subproblem satisfies the one-hot constraint. To the best of our knowledge, the partitioning method proposed by Nishimura et al. 30 is the first to focus on the one-hot constraint. This method is applicable to the double-constrained problems , such as the assignment problem and the traveling salesman problem. However, the extracted subproblems still contain infeasible solutions for which parameter λ must be adjusted. In this study, we propose two partitioning methods applicable to problems whose cost function involves a single one-hot constraint, as illustrated in Eq. (7). The first method is similar to the partition proposed by Nishimura et al., while the second method extracts subproblems comprising only feasible solutions and does not require adjusting parameter λ. The performance of the proposed methods is assessed for several Potts models, which are generalized Ising models whose cost function is defined by integer variables 45 . We demonstrate that the proposed methods efficiently obtain superior solutions.

Results
In this section, we propose efficient partitioning methods for solving large integer optimization problems under the one-hot constraint. In addition, we assess the performance of the proposed methods for several Potts models.
Proposed methods. We propose two partitioning methods: a multivalued partition and a binary partition.
These methods are summarized in Fig. 2. Both methods extract a subproblem that involves binary variables assigned to the tentatively selected components for each S i . The resulting subproblems include feasible solutions other than the current solution.
The multivalued partition extracts a subproblem with two or more components for each S i , as illustrated in Fig. 2(a). In addition to the tentatively selected component, the multivalued partition randomly selects one or more components for each S i , and then extracts a subproblem that comprises the binary variables assigned to the selected components. The extracted subproblem involves feasible solutions other than the current solution, and the randomly selected components are explored for each S i by optimizing the subproblem. However, the extracted subproblem still contains infeasible solutions, and the penalty term remains in the cost function of the subproblem. This partitioning method is similar to the partition proposed by Nishimura et al. Although subproblems are embedded using complete-graph embedding in the study by Nishimura et al. 30 , we employed  the subproblem-embedding algorithm developed in our previous paper 41 . Details of achieving a multivalued partition using the subproblem-embedding algorithm are provided in the Methods section.
The binary partition is summarized in Fig. 2(b). In addition to a tentatively selected component, the binary partition randomly selects exactly one component for each S i . Subsequently, new binary variables {y i } i = 1,2,...,N that represent "stay in the tentatively selected component (y i = 0)" or "transit to the randomly selected component (y i = 1)" are introduced for each S i , and a binary subproblem is constructed whose cost function is defined by {y i } i = 1,2,...,N . The cost function of the binary subproblem is derived in the Methods section. Thereafter, a subproblem of the binary subproblem is embedded into the D-Wave quantum annealer by the subproblem-embedding algorithm 41 . Here, the cost function of the binary subproblem does not involve the penalty term because all solutions in the binary subproblem are feasible. Therefore, the binary partition does not require adjusting parameter λ. In addition, a larger number of binary variables can be embedded into the D-Wave quantum annealer because the penalty term, which generates fully connected interactions between x i q ( ) and ′ x i q ( ) , is not involved. Consequently, the number of feasible solutions involved in the embedded subproblem is significantly increased using the binary partition. The binary subproblem can be regarded as one of the simplest cases of optimization under the half-hot constraint 46 . The penalty term of the half-hot constraint is given by and Q/2 components are extracted. The half-hot constraint is proposed to avoid the difficulty caused by the longitudinal magnetic field in the penalty term of the one-hot constraint. This difficulty is avoidable using the binary partition, which may contribute to improving solution accuracy. A disadvantage of the binary partition is that only two components are considered for each integer variable. As demonstrated in the following subsection, this leads to poor performance for the ferromagnetic Potts model.
Performance assessment. The performance of the proposed methods is evaluated for the following four types of Potts models on a cubic lattice with 10 × 10 × 10 integer variables: the ferromagnetic, anti-ferromagnetic, Potts glass 47 and Potts gauge glass 48,49 models. While the ground states of the ferromagnetic and anti-ferromagnetic Potts models are trivial, it is generally difficult to obtain the ground states of the Potts glass and Potts gauge glass models due to competing interactions. The cost function is given by where Q is set to 4, S i ∈ (1, 2, 3, 4), Δ ij ∈ (0, ±1), δ is the Kronecker delta function, and J ij represents the interaction between the nearest neighbors on the cubic lattice with the periodic boundary condition. The cost function is represented in QUBO form using the one-hot constraint as follows: Parameters J ij , Δ ij , and λ in each model are presented in Table 1. Δ ij ≠ 0 generates interactions between different components in the Potts gauge glass model, and Fig. 3 illustrates the local interactions generated by the first term of Eq. (10). Although it is generally difficult to determine an appropriate value of λ a priori, we can derive the lower bound of λ to correctly encode the original optimal solutions for the ferromagnetic and anti-ferromagnetic Potts models. The lower bound strongly depends on whether there exist infeasible solutions that set the first term of Eq. (10) to a smaller value than that of the original optimal solutions. λ > 0 is sufficient if such infeasible solutions do not exist; however, a sufficiently large value of λ is necessary if such infeasible solutions exist. For the ferromagnetic Potts model, the first term of Eq. (10) for an infeasible solution with is lower than that of the original optimal solution (e.g., x x 1, 0 ) by 3N. Because the second term in Eq. (10) increases by Nλ in this infeasible solution, λ > 3 is required. In this study, we set λ = 3.3 because an unnecessarily large value is not preferable, as mentioned in the Introduction section. While for the anti-ferromagnetic Potts model, the original optimal solutions minimize the first term of Eq. (10). Therefore, λ > 0 is sufficient, and we set λ = 1.0, which is the same value as J ij in the first term. However, for the Potts glass and Potts gauge glass models, the lower bound cannot be derived because the original optimal solutions are not trivial. At least, by setting λ > 3, we can restrict energy changes caused by a single-spin flip from the optimal solutions to be a positive value, and λ is set to 3.3 in this study. www.nature.com/scientificreports www.nature.com/scientificreports/ The optimization process demonstrated in this study is illustrated in Fig. 4. The original large problem is partitioned using three partitioning methods: random, multivalued, and binary partitions. The random partition does not address whether an extracted subproblem contains feasible solutions for each S i . A subproblem-embedding algorithm proposed in the literature 41 is used for embedding a subproblem into the D-Wave quantum annealer with defects in qubits and the interactions between them (details on embeddings are provided in the Methods section). After optimizing the embedded subproblem by the D-Wave quantum annealer under the parameter settings provided in Table 2, the variables in the subproblem are replaced by the best solution among the 1,000 solutions obtained using the quantum annealer. Subsequently, a greedy algorithm is executed by a conventional digital computer to recover the one-hot constraint and obtain an exact (local) minimum. In this algorithm, if there exist integer variables violating the one-hot constraint, the constraint is first recovered by extracting the integer variables and selecting exactly one component that minimizes the local energy for each integer variable. Then, an integer variable is randomly selected, and the tentatively selected component is replaced with one that minimizes the local energy. Refining the current solution is completed when all local energies are minimized. Finally, the best solution obtained in the procedure is updated, and the above processes are iterated. We then compare the solution accuracy for the three partitioning methods.   www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 5 represents the energies obtained by the three partitioning methods. The average, maximum, and minimum energies for 16 trials are plotted, and the same 16 initial states are used for each partitioning method. The horizontal axis represents the number of iterations, which is the number of subproblem optimizations performed by the D-Wave quantum annealer. The plot for the multivalued partition is shifted slightly to the left to avoid overlap with other plots. Figure 5(a,b) illustrate the energies obtained for the ferromagnetic and anti-ferromagnetic Potts models, respectively. The ground states of these models are trivial, and the minimum energy is −3 and 0 for the ferromagnetic and anti-ferromagnetic Potts models, respectively. Although the multivalued partition is expected to solve large optimization problems more efficiently than the random partition, the performances of the random and multivalued partitions are almost identical. The performance of the binary partition differs from the other methods, however; it is the lowest for the ferromagnetic Potts model, and the highest for the anti-ferromagnetic Potts model. Figure 5(c,d) present the energies obtained for the Potts glass and Potts gauge glass models, respectively. As expected, superior solutions are obtained with a smaller number of iterations using the multivalued partition rather than the random partition, in particular for the Potts gauge glass model. Of the three partitioning methods, the binary partition shows the highest performance for both the Potts glass and Potts gauge glass models.

Discussion
In this section, we discuss the differences between the three partitioning methods. The following three questions arise from the results presented in the previous section.
• Why is the multivalued partition not superior to the random partition for the ferromagnetic and anti-ferromagnetic Potts models? • Why is the performance of the binary partition the lowest for the ferromagnetic Potts model?   www.nature.com/scientificreports www.nature.com/scientificreports/ • Why does the binary partition exhibit the highest performance for all models except for the ferromagnetic Potts model?
A possible answer to the first question is that, for the ferromagnetic and anti-ferromagnetic Potts models, improved feasible solutions can be obtained through infeasible solutions because there likely exist infeasible solutions whose energy is lower than that of a current feasible solution in a neighbor. Figure 6(a) presents a simple example for the one-dimensional ferromagnetic Potts model. We assume that the binary variable enclosed by the green rectangle is extracted as a one-variable subproblem, which is one of the simplest cases of the random partition. The energy change caused by flipping the extracted binary variable is −2J + λ because two interactions are simultaneously recovered (−2J) and the one-hot constraint is violated (+λ). If λ < 2J, flipping the binary variable decreases the energy despite violating the constraint. It should be noted that λ > J is sufficient to correctly encode the ground states of the one-dimensional ferromagnetic Potts model. This is because the energy of the lowest-energy infeasible states, in which two components are commonly selected for each S i , is −2NJ + Nλ and must be larger than that of the ground states (−NJ). Consequently, if λ is appropriately tuned (J < λ < 2J), the current solution is updated to the infeasible solution by optimizing the subproblem, and superior feasible solutions are obtained via the infeasible solution. Whether energy changes become negative or not in spite of violating the one-hot constraint strongly depends on the number of simultaneously recovered interactions. For the ferromagnetic and anti-ferromagnetic Potts models without competing interactions, many interactions can be simultaneously recovered by violating the one-hot constraint. As a result, the multivalued partition is not effective in improving the solution accuracy for these models. In contrast, for the Potts glass and Potts gauge glass models with competing interactions, the performance of the multivalued partition is superior to that of the random partition.
In answer to the second question, subproblems that can eliminate domain walls are rarely extracted by the binary partition. Figure 6(b) presents one of first excited states, which is commonly observed in the optimization of the ferromagnetic Potts model. The 10 variables in Fig. 6(b) are divided into two domains: five variables S 1 , ..., S 5 are aligned to q = 1 in one domain, while the remaining five variables S 6 , ..., S 10 are aligned to q = 2 in the other domain. The boundary between the domains is referred to as the domain wall. To improve the current solution, an extracted subproblem must contain one of the ground states because the current solution is the first excited state. For example, to align all integer variables {S i } i = 1,2,...,10 to q = 1, component q = 1 must be selected for integer variables S 6 , ..., S 10 . The probability of component q = 1 being selected for S 6 , ..., S 10 is (1/3) 5 = 1/243 because, in addition to the tentatively selected component, the binary partition randomly selects one component for each S i . This probability exponentially decreases with respect to the number of variables, and the extraction of only two components is not suitable for the ferromagnetic Potts model. We further conjecture that the binary partition exhibits poor performance for optimization problems containing ferromagnetically ordered domains, and that the concomitant use of the binary and multivalued partitions may be preferable for such problems.
The answer to the third question is that there exist several binary subproblems that can improve the current solution. Figure 6(c) illustrates the local interactions in the anti-ferromagnetic Potts model. The current solution is one of the first excited states in which the local energy with respect to S 1 and S 4 is not minimized, which we say "the interaction between S 1 and S 4 is broken" in this paper. Assuming that the integer variable S 4 is updated to improve the current solution, there are two binary subproblems that can improve the current solution. Therefore, the disadvantage of the binary partition, which is that only two components are considered for each integer variable, is mitigated for the optimization of the anti-ferromagnetic Potts model. We can thus exploit the advantages of the binary partition, which are that the extracted subproblems contain a larger number of feasible solutions and that the adjustment of parameter λ is not required. This is also the case for the Potts glass and Potts gauge glass models, in which competing interactions generate several binary subproblems that improve the current solution. Figure 6(d) presents a simple example for the Potts gauge glass model. One of the ground states and first excited states are illustrated at the top of Fig. 6(d), where the interaction depicted by a dashed line represents the broken interaction. There exist no configurations that minimize all of the local energies due to the competing interaction between different components, and one interaction is broken even in the ground state. Suppose that we update integer variable S 4 to improve the current solution in the first excited state, then there are two binary subproblems that can improve the current solution, as illustrated at the bottom of Fig. 6(d). One subproblem recovers the interaction between S 1 and S 4 , while the other subproblem recovers the interaction between S 3 and S 4 . The competing interactions generate two binary subproblems that improve the current solution, and each subproblem recovers different interactions. Thus, the disadvantage of the binary partition is mitigated as long as Q is not very large. It should be noted that although the number of binary subproblems that improve the current solution increases as Q is increased for the anti-ferromagnetic Potts model, this number does not increase for the Potts gauge glass model.

Conclusion
In this study, we proposed two partitioning methods to efficiently solve large optimization problems under the one-hot constraint using the D-Wave quantum annealer. The performance of the proposed methods was assessed for the ferromagnetic, anti-ferromagnetic, Potts glass, and Potts gauge glass models. Of the three partitioning methods, the binary partition showed the highest performance for all models except for the ferromagnetic Potts model. The advantages of the binary partition are that it enables embedding a larger number of binary variables and does not require adjusting parameter λ. However, its disadvantage is that only two components are considered for each integer variable. Although this disadvantage leads to poor performance for the ferromagnetic Potts model, the effect is mitigated for optimization problems that have many binary subproblems improving the current solution, such as the anti-ferromagnetic Potts model, and for optimization problems with competing interactions. Although the multivalued partition exhibits a better performance than the random partition for www.nature.com/scientificreports www.nature.com/scientificreports/ Potts glass and Potts gauge glass models, we did not identify problems for which the multivalued partition is most suitable. Future studies should focus on constructing algorithms that can efficiently solve the ferromagnetic Potts model using the binary partition. In addition, the performance of the proposed methods should be assessed for various optimization problems, such as the graph coloring problem whose cost function is represented as the Hamiltonian of the anti-ferromagnetic Potts model.

Methods
Details on partitioning and embedding are provided in this section.
Subproblem-embedding algorithm. In this subsection, we briefly explain the subproblem-embedding algorithm developed in a previous study 41 . This algorithm aims to quickly find minor embeddings of subproblems to efficiently implement large-neighborhood local searches using the D-Wave quantum annealer. Given a large optimization problem whose cost function is represented in QUBO form, this algorithm embeds binary variables one by one into the Chimera graph. After a randomly selected binary variable is embedded into the Chimera graph first, this algorithm embeds a binary variable interacting with the already embedded binary variables into the graph. The latter procedure is iterated until all qubits in the Chimera graph are used. It is generally necessary to assign several qubits to one binary variable and extend a chain to correctly represent interactions between binary variables on the Chimera graph. This algorithm implements Dijkstra's algorithm to greedily determine how to extend chains. Although chains often cannot be extended by only unused qubits, this difficulty can be avoided in the embedding of subproblems because it is not necessary to embed all binary variables. In this case, the algorithm stops attempting to embed the binary variable and attempts to embed other binary variables that can be easily embedded. Thus, extracting and embedding a subproblem is simultaneously implemented in this algorithm, and the resulting subproblem comprises only binary variables that can be easily embedded. Therefore, the computational time is significantly lower than in Cai's algorithm 50 . In addition, this algorithm can easily deal with hardware defects by implementing Dijkstra's algorithm on a Chimera graph with defects.
It is the random partition to directly apply this algorithm to extract and embed a subproblem of integer optimization problems because this algorithm does not address whether an extracted subproblem contains feasible solutions for each S i or not. To combine the multivalued partition with this algorithm, the order of the binary variables embedded into the Chimera graph must be appropriately specified, as described in the following subsection.
Multivalued partition. The multivalued partition requires that the binary variable assigned to the tentatively selected component must be embedded into the Chimera graph. In addition, more than two binary variables should be embedded for each integer variable to distinguish the multivalued and binary partitions. On the other hand, the subproblem-embedding algorithm extracts and embeds a subproblem comprising binary variables that can be easily embedded. Therefore, in order to combine the multivalued partition and the subproblem-embedding algorithm, it is needed to appropriately specify the order of the binary variables embedded into the Chimera graph.
First, an integer variable is randomly selected, and binary variables assigned to the selected integer variable are embedded into the Chimera graph. Then, to determine binary variables to be additionally embedded into the Chimera graph, we select an integer variable as follows: 1. An already embedded binary variable x embedded is selected in the order of being embedded into the Chimera graph. 2. An integer variable S ctr to which x embedded is assigned is selected. 3 .. , criterion 1 is important than criterion 2 to avoid embedding independent integer variables. For the reminder of the binary variables assigned to S i , we prioritize criterion 2 to achieve the multivalued partition. It should be noted that the number of components embedded into the Chimera graph is not uniform for each integer variable because the subproblem-embedding algorithm embeds only binary variables that can be easily embedded. If only one component can be embedded, the integer variable is excluded from the subproblem. Table 3 represents the average number N S (Q embed ) of embedded integer variables with Q embed components per subproblem. The performance is assessed for embedding the Potts gauge glass model on the cubic lattice into D-Wave 2000Q_2 with defects and is averaged over 1,000 trials. 65.7(=14.5 + 8.0 + 43.2) integer variables are embedded into the Chimera graph on average, and the average number of embedded binary variables is as follows: It should be noted that, to distinguish the multivalued and binary partitions, Q embed > 2 is required for most integer variables. All four components are embedded for 65.8% of the integer variables in the subproblem, indicating that we can embed the multivalued subproblem that is distinct from the binary subproblem.
Binary partition. To solve large optimization problems using the binary partition, the cost function of the binary subproblem must be derived from the cost function of the original large problem. The general form of the local energy between S i and S j is given by  . The binary partition extracts a binary subproblem by randomly selecting one component in addition to the tentatively selected component for each integer variable. The local energy of the binary subproblem in QUBO form is given as follows: ij ij i j ii i jj j ij ij ij ij ij where y i ∈ {0, 1}, α i and β i denote the tentatively selected component and randomly selected component for S i , respectively, and y i = 0(y i = 1) indicates "stay in the tentatively selected component α i " "("transit to the other component β i ")" It should be noted that the cost function of the binary subproblem does not contain the penalty term because all solutions in the binary subproblem satisfy the one-hot constraint. The problem graph of the binary subproblem extracted from the three-dimensional Potts model is a cubic lattice with bond dilutions. The density of the interactions in the binary subproblem is lower than that in the multivalued subproblem because the cost function of the binary subproblem does not contain the penalty term, which generates partially fully connected interactions between x i q ( ) and ′ x i q ( ) . The average number of embedded binary variables is 408 when the binary partition is used, and only 225 when the multivalued partition is used. Furthermore, all configurations in the binary subproblem satisfy the one-hot constraint, while the configurations in the multivalued subproblem do not. Therefore, the average number N feasible of feasible solutions involved in the embedded subproblem is considerably increased using the binary partition. Table 4 illustrates log 10 N feasible in a subproblem embedded by the multivalued and binary partitions combined with the complete graph embedding 44 and subproblem-embedding algorithm 41 for the Potts gauge glass model on the cubic lattice.   Table 4. Average number of feasible solutions in an embedded subproblem: log 10 N feasible .


extfieldsofstudy:
Computer Science
Medicine

text:
Sensor Data Acquisition and Multimodal Sensor Fusion for Human Activity Recognition Using Deep Learning.

In this paper, we perform a systematic study about the on-body sensor positioning and data acquisition details for Human Activity Recognition (HAR) systems. We build a testbed that consists of eight body-worn Inertial Measurement Units (IMU) sensors and an Android mobile device for activity data collection. We develop a Long Short-Term Memory (LSTM) network framework to support training of a deep learning model on human activity data, which is acquired in both real-world and controlled environments. From the experiment results, we identify that activity data with sampling rate as low as 10 Hz from four sensors at both sides of wrists, right ankle, and waist is sufficient in recognizing Activities of Daily Living (ADLs) including eating and driving activity. We adopt a two-level ensemble model to combine class-probabilities of multiple sensor modalities, and demonstrate that a classifier-level sensor fusion technique can improve the classification performance. By analyzing the accuracy of each sensor on different types of activity, we elaborate custom weights for multimodal sensor fusion that reflect the characteristic of individual activities.

Introduction
As wearable devices are widely used nowadays, Human Activity Recognition (HAR) has become an emerging research area in mobile and wearable computing. Recognizing human activity leads to a deep understanding of individuals' activity patterns or long-term habits, which contributes to developing a wide range of user-centric applications such as human-computer interaction [1,2], surveillance [3,4], video streaming [5,6], AR/VR [7], and healthcare systems [8,9]. Although activity recognition has been investigated to a great extent in computer vision area [10][11][12], the application is limited to a certain scenario that equips pre-installed cameras with sufficient resolution and guaranteed angle of view. On the other hand, wearable sensor approach allows continuous sensing during daily activities without spatio-temporal limitation [13,14], since devices are ubiquitous and do not require infrastructural support.
Although many commercial fitness trackers or smartwatches adopt wrist-mountable form due to the ease of accessibility, no consensus has been made regarding the positioning of sensors and their data acquisition details. For instance, there are many public datasets available which collect data from Inertial Measurement Units (IMU) attached on different parts of the human body such as ear, chest, arm, wrist, waist, knee, and ankle with various parametric settings [15][16][17][18]. Data acquisition details may vary depending on the type of activities investigated in the research. For example, a sensor on • Build a testbed system and collect activity data in both real-world and lab environments that contain high variations in the signal patterns.

•
Analyze activity data to empirically study the data acquisition details regarding the sensor position and parametric settings, and investigate the characteristic of each sensor modality on detecting different types of ADLs. • Perform sensor fusion on multimodal sensor data by implementing a two-level ensemble model that combines single-modal results from the deep learning activity classifier.
The rest of this paper is organized as follows. Section 2 summarizes previous work closely related to activity recognition and sensor fusion research. Section 3 describes our testbed and experiment protocols used in the study. Section 4 discusses details of our deep learning framework and sensor fusion model developed for data analysis, and Section 5 presents the evaluation results and findings. Section 6 continues the discussion on current limitations and future work. Finally, Section 7 states the final remark.

Related Work
Human activity recognition has been investigated using various sensor modalities [26]. Besides the body-worn sensors, object sensors such as Radio Frequency Identifier (RFID) tags are widely used to infer human activities based on the movements detected in smart home environments [27]. In addition, ambient sensors that capture the change of environments including light, sound, temperature, pressure, and humidity are combined to detect activity events [28]. This section introduces literature related to on-body sensor-based HAR, especially the one investigating the placement of sensors, mobile sensing and public datasets, multimodal HAR and their sensor fusion techniques.

Sensor Positioning
Selecting optimal placement of sensors on different body positions has been widely studied for various purposes of activity recognition that range from gait monitoring or fall detection for elderly subjects to daily living activities of healthy subjects usually with no impairments [29]. By assuming a single IMU on-body sensor, authors of [30] defined the cross-location activity recognition problem and proposed a transfer learning based HAR model, which adapts a new sensor position based on known sensor placements. We limit the scope of study to ADLs and focus on finding feasible sensor location utilizing minimum number of body-worn sensors.
Authors of [19] classify activities into four broad levels of activity groups such as very low level, low level, medium level, high level (e.g., from static to dynamic level) based on the compendium of physical activities that leverages the rate of energy expenditure. Seven accelerometer sensors are placed on wrist, arm, waist, chest, ankle, and knee in addition to the ear-worn sensor. According to the experiment result, knee sensor shows the best performance in recognizing high-level activity (e.g., running), while the next best positions were ear and arm sensors. However, the authors claim that feasibility should be taken into consideration because chest or knee-worn sensors would interfere with daily activities. Therefore, consistent yet realistic position of sensor is recommended, depending on the target application.
Similarly, another work [20] investigated seven activities including walking, walking upstairs, walking downstairs, running, standing, sitting, and lying using data acquired from six sensors placed on wrist, ankle, thigh, chest, lower back, and hip. Experiment results show that the hip was the best single location for activity recognition using SVM, and combining two sensors from different locations can achieve a significant increase in the recognition accuracy. The authors stated that combinations of sensors worn on the upper half and another on the lower half of the body can cover simple daily activities. However, the paper leaves identification of finer-grained activities (e.g., sitting and sitting working at a computer) as a future work.

Mobile Sensing and Public Datasets
Many efforts have been made to collect mobile sensing data for HAR, and there are many datasets publicly available. The ExtraSensory dataset [15,16] aims to analyze behavioral context in-the-wild from mobile sensors. It uses sensors from everyday devices including smartphone and smartwatch, where users were engaged in regular natural behavior and self-reported their context with 100 labels. Among many sensors, motion-reactive accelerometer data was sampled with high frequency: 40 Hz from a smartphone and 25 Hz from a smartwatch for 20 s every one minute. However, as the dataset adopts heterogeneous sensors on limited position of the human body, it is not sufficient to evaluate the recognition performance according to the sensor position and sampling rate.
The UCI smartphone dataset [17] continuously records IMU sensor data from a waist-mounted smartphone with 50 Hz sampling rate. In a lab environment, test subjects performed six activities including walking, walking upstairs and walking downstairs, sitting, standing, and lying. Likewise, the UniMiB-SHAR dataset (University of Milano Bicocca Smartphone-based HAR) [18] takes accelerometer samples at a maximum frequency of 50 Hz from a smartphone located in the front trouser pockets of the subjects. Both UCI and UniMiB-SHAR datasets leverage only one sensor mounted on a single point of the body, which limits the scope of experimental scenario.
Adopting high sampling rate may provide sufficient information for data analysis, but imposes burden on the system with large data size and computation load. On the other hand, applying low sampling rate may fail to capture intrinsic attributes of each activity, which can differentiate one from another. To the best of our knowledge, there has been no consensus made regarding the sampling frequency of motion-reactive sensors, and previous work accepts rates commonly used in literature for activity recognition from data acquired through smartphones. Therefore, we aim to investigate the optimal sampling frequency of body-worn sensors attached on different positions of the human body to collect activity data with low burden.

Deep Learning in Ubiquitous Computing
Deep learning model is renowned for its feature learning ability through neural networks, which accepts a large amount of raw data for training and identify unseen data through knowledge transfer methods [31]. As a derivative of Recurrent Neural Network (RNN), LSTM network model is known to perform well on extracting signal patterns in the input feature space. By accepting input data that spans over long sequences, it is specialized in time series problems, especially in time-series prediction domain. The gated architecture of LSTM manipulates the memory states that model the local dependencies on features. HAR is based on the assumption that motion-reactive sensor signals present discriminative patterns corresponding to different activities. As such activities usually last for a certain period of time when they occur, HAR is regarded as a time-series classification problem with temporal dependency, where input data that are close in space may be dependent while distant sequence of samples in time are assumed as independent. Thus, LSTM network model is renowned for its performance on HAR domains [32]. Papers [33,34] also compared deep-learning methods for HAR, and verified that variants of the LSTM architecture lead to flexible and robust results in several sensor-based HAR datasets.

Multimodal Sensor Fusion
Many efforts have been made to combine multimodal sensor data to ensure the classification accuracy by taking the advantage of diverse data sources. Sensor fusion in HAR field is usually divided into two categories, feature-level and classifier-level fusion. Authors of [16] referred to the feature-level fusion as early fusion, which merges feature data extracted from each modality prior to the classification procedure. Using a single feature vector for the training model is straightforward, but feature compatibility issues regarding heterogeneous sampling frequencies and configuration parameters can deteriorate the classifier performance [35]. Contrarily, classifier-level fusion aims to combine the classification results given by the base-learners trained on different sensor modalities. As class-probability predictions are uniform in their expression, classifier-level ensemble learning techniques are widely adopted to integrate multiple machine learning models [34][35][36]. In this paper, we apply two-level stacking and voting ensembles, which fall into the classifier-level sensor fusion technique.

Testbed and Data Acquisition
This section first describes the system architecture of the proposed body-worn sensor testbed that acquires human activity data. Then, the experimental protocols for data collection are explained in detail.

Testbed System
We build a testbed system that acquires continuous motion data from on-body IMU devices. We adopt a MPU-9250 IMU sensor that samples tri-axial accelerometer, gyroscope, and magnetometer data at maximum 100 Hz frequency. As shown in Figure 1a, each IMU device is implemented on a ESP8266 Micro Controller Unit (MCU) using a Wi-Fi networking and nRF24L01 RF module, which is powered by a 3.7 V-1200 mAh micro-USB rechargeable Li-Po battery. The corresponding battery capacity is empirically verified to guarantee at least 12 h of continuous data collection with maximum sampling frequency. IMU devices are designed to be attached on both wrists, upper arms, ankles in addition to chest and waist of a test subject. Medical Velcro straps are used to fasten devices on the body to reduce signal noise.
We implement a relay module that operates as a data collection hub between IMU devices. As a large amount of activity data is accumulated at each IMU device, we choose Wi-Fi protocol to assure fast and stable data transmission through high bandwidth network channels. Therefore, the relay module first transceives operation commands between IMU devices using the RF protocol to initiate time-synchronized one-to-many communication channels between IMU devices. After establishing a Wi-Fi ad-hoc network through the soft-AP setting, the relay module then can retrieve activity data from each IMU device. When data gathering is completed, the relay module transmits the collated data to an Android mobile device over the Wi-Fi network. Overall communication protocol is illustrated in Figure 2.

ADLs
Statistics Canada annually publishes Time-use survey of the General Social Survey (GSS) [37], which includes respondents' time spent on a wide variety of day-to-day activities in Canada. In addition, Statistics Korea issues Koreans' average time-use [38] categorized by mandatory, personal, and leisure time activities. Authors of [18] surveyed ADLs and their occurrence in public HAR datasets, and we find that six primitive activity labels including walking, standing, sitting, lying, walking upstairs, and walking downstairs can cover daily activities listed in time-use survey results.
In addition to the primitive activities, we take three combinatorial activities into account: eating, driving a car, and moving in a car as a passenger. For example, eating activity is an integration of sitting activity (i.e., primitive) and actions using cutlery during a meal (i.e., secondary activity). Similarly, driving a car can be interpreted as a combinations of sitting activity and actions related to handling the steering wheel. In addition, moving in a car as a passenger can include additional actions such as using smartphones while sitting. Table 1 specifies nine activity labels used in our experiments.

Experiment Protocols
As a pilot experiment, we collected activity data from five researchers working at a research institute. Test subjects were composed of two males and three females with ages ranging from 34 to 48 years, heights ranging from 160 cm to 175 cm. All subjects were right-handed with no bodily impairments. Test subjects participated with informed consent and instructed to perform daily activities freely with no restrictions on their body movements. A total of three experiment sessions were carried out on three different workdays, accumulating about three hours of activity data per each subject.
Two different experiment protocols were designed for data collection: a real-world scenario and a controlled environment scenario. Experiments in a real-world scenario were conducted twice during the lunch time of ones' daily routine. No instructions were given to the test subject with regard to the experiment protocol, so test subjects behaved as usual in a natural manner. An instructor monitored the whole experiment session and manually tagged the corresponding activity labels. The sequential activities labeled during a round trip between the workplace and a restaurant is summarized in Table 2. On the other hand, an instructor informed each target activity and test subjects followed the instruction during the controlled environment scenario. One iteration of experiment is performed to collect activity labels presented in Table 3. Driving or Moving 5 7 Walking 10 8 Sitting 5  9 Walking 20 10 Walking Upstairs 3 Total 81 Table 3. Experiment protocol in the controlled environment (one iteration).

Final Testbed Configuration
At the initial stage of the system design, we developed eight body-worn sensors as shown in Figure 1b. During the primitive study, we observed that applying an excessive number of sensors interrupts identifying the unique characteristic in each activity data. In addition, sensors attached on both ankles did not show significant difference in distinguishing the activities. Therefore, we adopted the sensor on right ankle, and for the same reason, we excluded sensors on upper arms. Participants experienced significant discomfort incurred by the sensor on the chest; we dropped the chest sensor in the final testbed configuration for its feasibility.

Data Analysis
We extract activity data and their corresponding labels, and segment them using a fixed-width sliding window with 50% overlap. Sliding window size is a base unit to be applied to the data analysis, which can be interpreted as the resolution of activity. By leveraging a long window size, we can handle a stream of activity data as one phase to observe long-term activity patterns or trends, but may fail to extract distinguishing characteristics of each activity if discrete activities are contained in one segment of activity data. On the other hand, a short window size may be suitable when targeting finer primitive actions that occur in a very short period of time, especially with regard to low-level hand gestures such as reach, grasp, and release.
We aim to perceive coarser-grained activities related to modes of locomotion, which usually lasts for several seconds. Previous work in the HAR literature empirically adopted 5.12 s-long window size [35,39]. Based on the average free-living cadence (i.e., 76 ± 6 steps per minute [40]), we leverage the duration of one gait cycle that ranges from 1.46 to 1.71 s as a baseline. Therefore, we evaluate the recognition accuracy by applying a various sliding window that starts from 1.4 to 5 s. The result in Figure 3 indicates that the 2.5 s-long sliding window outperforms other parameter settings in our dataset. Deep learning classification models are capable to give better inference results using raw sensor data than the feature data, where pre-processing transformations usually decrease the classification accuracy [21]. Thus, raw sensor data is directly fed into LSTM cells for offline activity classification with no additional feature engineering [22,23].

Deep Learning Architecture
As a variation of RNN architecture, LSTM [25] network is renowned for its capability to learn long-term dependencies in time-series data by the selective forget-and-update procedure. The gated architecture operates the memory states that model the temporal dependencies. Through a sigmoid layer called forget gate, the LSTM cell decides how much information should be saved based on the time dependency of the previous data. Then, the input gate layer creates new candidate values from the input and updates the old cell state to a new one. These LSTM recurrent layers automatically learn feature representations and model the temporal dynamics between cell activation.
Our deep learning framework accepts time-series of sensor data as an input vector and convert them to a probability vector at the output for activity classification. A single data vector is assigned per time step, where a pre-defined number of time steps compose a unit (e.g., 2.5 s) of the time series as shown in Figure 4. The input vector consists of 3-axis accelerometer values sampled from four sensors for each cell, where we apply gyroscope and magnetometer readings for sensor fusion. We implemented a LSTM network using Google's TensorFlow framework by manually tuning the following parameters: batch size = 150, number of stacked LSTM layer = 1, hidden neuron size = 32, number of steps = 25, keep probability = 0.8.

Sensor Fusion
In addition to the accelerometer data, we also consider values collected from gyroscope and magnetometer sensors for multimodal sensor fusion to ensure the activity recognition performance. Leveraging the LSTM network as a base learner, we employ classifier-level fusion to combine results from multimodal data as shown in Figure 5. For each sensor modality including accelerometer, gyroscope, and magnetometer, prediction results from LSTM network indicating the classification probabilities of each activity class are accepted as meta-features. As a meta-learner, we first adopt stacking ensemble with 10-fold cross-validation, which trains the aggregator model on the stacked class-probabilities. We also apply simple soft-voting ensemble based on the weighted average of class-probabilities, where the weights are determined according to the prediction accuracy of the base-learner. Finally, the ensemble model makes combined predictions for the new multimodal input data.

Evaluation
In this section, we first evaluate the activity recognition performance with different configuration parameters applied using the accelerometer data. We empirically investigate the optimal sampling rate of body-worn IMU sensors and their position on the body. Then, we perform evaluation on the sensor fusion technique with multimodal sensor data.
In the classification problem, class imbalance is common because most datasets do not contain exactly equal number of instances in each class, and this phenomenon also applies to human activities during daily lives. F1-score is a commonly used measure in the class imbalanced settings [41]. In the multi-class classification domain, however, micro-averaging F1-score is equivalent to accuracy which measures a ratio of correctly predicted observation to the total instances. Micro-F1 reflects the unbalanced distributions of samples across classes, placing more weights on the dominant classes. On the other hand, macro-averaging F1-score is a harmonic mean of precision and recall measures, and takes both false positives and false negatives into account based on a per-class average [42].
We basically apply 10-fold cross-validation on the data to generalize the result to an independent dataset. By randomly partitioning the dataset into 10 equal-sized subsamples, it leaves a single subsample (i.e., 10% of total activity data) out as test data, while the remaining nine subsamples (i.e., 90% of data) are used as training data. In this way, the system is tested with unseen dataset. We also perform leave-one-subject-out (LOSO) cross-validation for reference by leveraging data augmentation [43][44][45][46], which is widely adopted in the deep learning research that requires large dataset for model training. We apply the jittering method to simulate additive sensor noise using random parameters as described in [47] for label-preserving augmentation. In this way, we handle limited data availability and secure a large amount of training data to increase robustness of our deep learning model. The mean performance of the cross-validation results is used in the remainder of this paper.

Impact of Sampling Rate
We first tested the recognition accuracy with different sampling rates applied. Starting from 100 Hz, we sub-sampled the sensor data and reduce the frequency to 10 Hz. The number of data was 20,132 when sampling at 100 Hz, and the number slightly varied with different sampling rates applied. Figure 6 shows the test results using four sensors attached on the left and right wrist, right ankle, and waist. As the sampling rate varies, the number of sample entries included in the 2.5 s-long sliding window vary accordingly. From the result, we found that a continuous time-series data with only 25 sensor samples that cover 2.5 s of body movement can successfully recognize the activity without significant loss in the accuracy. Figure 6 also represents the execution time to perform a 10-fold cross-validation on a desktop computer equipped with Intel Core i7 3.47 Hz CPU. Reduced sampling rate decreases the time to retrieve results from the neural network as the resolution of data becomes coarse. It is an encouraging result for mobile platforms, as devices usually do not have enough resources such as computation power, battery capacity and storage capability. By leveraging lower data sampling frequency, a mobile-based HAR system can operate in a cost-effective manner. Our findings are in line with the policy of many health trackers in the market that usually support low-frequency sampling rate. For example, a motion and physiological signal tracker Empatica E4 [48] supports the accelerometer data sampling at 32 Hz. In addition, a motion sensor manufacturer [49] provides a calculation tool that estimates maximum data acquisition time depending on the sampling rate of a sensor. According to the calculation, a MetaMotionC IMU sensor can operate up to 11.6 h at 12.5 Hz, while the battery lasts only 5.8 h at 25 Hz sampling frequency, respectively. As HAR applications in-the-wild require continuous and seamless data collection with the day-long battery life, it is worthwhile to identify that low-frequency activity data is still effective in recognizing daily activities.

Impact of Sensor Position and Combinations
To find an optimal sensor combination for the HAR systems, we investigated the recognition accuracy with different numbers of sensors attached on various locations of the body. The experiment with four sensors on left and right wrists, right ankle, and waist is used as the baseline, which is represented as experiment no.1 in Figure 7. The baseline experiment shows 93.0% of recognition accuracy as shown in Figure 7a. Starting from this baseline, we excluded each sensor in turn and extensively evaluated all sensor combinations on various body positions. Experiment no.3 that leverages right wrist, right ankle, and waist sensors presents the most competitive results among experiments using three sensors, even when considering inter-subject learning model as shown in Figure 7b. With only two sensors, positioning one on the upper half (e.g., wrists) and the other on the lower half (e.g., ankle or waist) of the body as experiment no.6 performs better than placing both sensors on the upper half (experiment no.11) or the lower half (experiment no.9) of the body. We can expect performance improvement when appending an additional sensor on the lower half of the body to the sensor combinations on the upper half of the body by comparing experiment results of no.2, no.5 to no.11. Similarly, we can observe performance gain by attaching a sensor to the upper half of the body, especially on the right wrist (experiment no.3) to the sensor combinations on the lower half of the body (experiment no.9). From the aspect of the wrist sensor, we found that an additional single sensor on the ankle, waist, and the opposite side of the wrist in sequence can influence the performance when adopting two sensor locations as shown in experiment results [no.6, no.8, no.11] and [no.7, no.10, no.11].  6, no.7], and [no.8, no.10] pairs, we observed that the sensor attached on right wrist tends to perform better than the left wrist. This is due to the fact that all five test subjects were right-handed persons, and eating activity is included in the experiment. Next, we explored the impact of sensors on the lower half of the body by comparing results of experiment [no.2, no.5], [no.6, no.8], and [no.7, no.10] pairs, and noticed that the sensor on the right ankle plays a more critical role than the sensor on the waist in the activity recognition.
By taking the trade-off between accuracy and feasibility into account, we concluded that placing sensors on the right wrist and right ankle can give reasonable performance when adopting two sampling positions, and an additional sensor may be placed on the waist to achieve better recognition accuracy.

Impact of Multimodal Sensor Fusion
For multimodal sensor fusion, we adopt stacking ensemble with various combinations of sensor modality. The meta-learner accepts stacked class-probabilities of each sensor modality to train the classification model. As shown in Figure 8, we compare the performance of Random Forest, kNN, and SVM when applied as a meta-learner for stacking ensemble, and adopt Random Forest algorithm. We also apply voting ensemble as the meta-learner that combines the accelerometer, gyroscope, and magnetometer sensor data. In hard-voting, the final class label is determined as the class label that has been predicted most frequently by the classification models. On the other hand, soft-voting predicts the class labels with the highest class probability, averaged over all the individual classifiers. Weights are determined based on the accuracy of each sensor modality, while each sensor may show different performance depending on the type of activity. Therefore, we elaborate weights according to the characteristic of each activity. First, we evaluate the performance of single sensor modality in activity recognition. Figure 9 illustrates the recognition accuracy of each sensor data applied to the LSTM base-learner model. The accelerometer shows the best performance (92.18%) among IMU sensors, followed by magnetometer (85.22%) and gyroscope (78.33%). As each sensor measures different physical aspects, contribution in recognizing ADLs differs according to the characteristic of activities. For example, an accelerometer is useful in measuring changes in velocity and position, and in consequence detecting small movements. By leveraging multiple axes, an accelerometer is also widely used as an absolute orientation sensor in the up-down plane. A gyroscope measures changes in orientation and rotational velocity, but requires calibration from a known orientation to achieve accuracy due to high amount of drift. A magnetometer is useful to determine absolute orientation from magnetic north with low drift over time, but shows poor accuracy for fast movements. Therefore, these multimodal sensors are usually combined to compensate the performance of each other. To investigate the contribution of each sensor on recognizing different activities, we evaluate the recognition accuracy of each sensor modality on individual activities as shown in Figure 10. Activities illustrated in Figure 10a are a group of combinatorial activities of sitting, including driving, moving in a car as a passenger, and eating. As a vehicle in operation generates consistent vibration and periodic acceleration, both moving in a car as a passenger and driving activities are well distinguished from sitting activities in different environment. Additionally, eating activity includes repetitive hand movements while sitting, which incurs various yaw, pitch, and roll rotations that can be effectively captured through the gyroscope sensor. Although driving, moving in a car, and eating activities can be considered as semi-static because activities occur during sitting, individual sensor modality shows high recognition accuracy independent from each other due to additional gestures. Both combinations of sensor fusion (i.e., accelerometer + gyroscope and accelerometer + magnetometer) show a similar contribution on the performance improvement. Figure 10b presents the recognition accuracy of a group of static activities including sitting, standing, and lying. In static position, gravity of the earth may be the major reading of the accelerometer, which accordingly reflects the tilt angle of the body posture in the signal. No remarkable measurements take place at the gyroscope sensor, while the magnetometer measures heading of the sensor. Therefore, the accelerometer shows outstanding performance in detecting different static postures. Furthermore, fusing additional sensors to the accelerometer data even degrades the accuracy in detecting sitting and standing activities. Figure 10c summarizes results of a group of activities that involves walking, including walking upstairs and downstairs. As walking activity involves highly dynamic physical movements, it is detected with high accuracy using each sensor modality, and sensor fusion also contributes to the recognition performance. However, experiments in making distinctions between stair walking and walking on the plain produced equivocal results. As shown in Figure 11, both walking upstairs and walking downstairs are incorrectly predicted as walking with high probabilities, where the x-axis of the confusion matrix stands for the ground truth and the y-axis denotes predicted labels. According to our experiments, fusing the accelerometer and gyroscope sensor data through stacking ensemble led to substantial improvement in detecting stair walking activities. (c) Figure 10. Recognition accuracy of a group of activities using multiple sensor modalities. (a) A group of combinatorial activities of sitting, including driving, moving in a car (as a passenger), and eating; (b) A group of static activities including sitting, standing, and lying; (c) A group of activities that involves walking, including walking upstairs and walking downstairs.
In soft-voting ensemble, sensor modalities are weighted by their demonstrated accuracy shown in Figure 9, which covers all activity types. Weight proportion of the accelerometer, gyroscope, and magnetometer is applied as 0.36:0.31:0.33 by default. To consider different characteristics of each sensor on detecting various activities, we choose weights by the activity basis. Figure 12 summarizes the proportion of weights computed using the recognition accuracy of sensor modalities. For a group of combinatorial activities of sitting, including driving, moving in a car as a passenger, and eating activities, both the accelerometer and magnetometer had similar weights (e.g., 0.35∼0.38). Among static activities, weight proportion of sitting showed analogous result to standing: 0.49:0.2:0.31 for the accelerometer, gyroscope, and magnetometer, respectively. As weights for walking upstairs and walking downstairs are noticeably different from other activities, applying activity-specific weight attributes substantial increment to the recognition accuracy. As presented in Figure 13, customized weight of each activity contributed to accuracy increase of 7.65% and 36.97% in the case of walking upstairs and downstairs, respectively.   Finally, Figure 14 summarizes the sensor fusion results using both stacking and voting ensemble techniques. Overall, the contribution of gyroscope to the recognition accuracy is higher than the magnetometer when combined with the accelerometer data. Although the gyroscope sensor data is used without calibration procedure, the results indicate that applying sensor data from additional modalities can contribute to improving the recognition accuracy of nine activities up to 93.07%. In addition, soft-voting achieves higher overall performance (94.47%) than hard-voting (93.48%), because it can give more weight to highly confident votes. By applying multimodal sensor fusion, we achieved 2.48% of performance improvement compared to the single modal-based activity recognition. Even considering the balance between precision and recall measures, our model performs well and seems to be sufficiently competitive with uneven class distribution.

Limitations and Future Work
Test subjects: As this work is initially designed to investigate the testbed system before actual distribution on large-scale subjects in real-world environment, experiments took place within the institution on a small number of test subjects. Based on our findings, we plan to adopt off-the-shelf wearable sensors in a lightweight form to improve the practicality of our testbed system. By tuning the systematic parameters, our refined testbed will enable a large-scale and long-term experiment that captures human activity in-the-wild.
Meta-learner: We have evaluated the meta-learner using classical classification algorithms including RF, kNN, and SVM, which are widely adopted in the literature. We plan to implement more classification algorithms with enhanced performance [50], and also extend our meta-learner using deep learning models to improve the recognition accuracy.
Future work: We further expect to gain a deep understanding on individuals' daily routines or long-term habits through a longitudinal study, which can contribute to the anomaly and pattern detection that may change over time. We also aim to apply additional sensing modalities such as environmental (e.g., background audio recordings, GPS or network locations) and physiological (e.g., electrocardiogram, electro-dermal activity) sensors. For sensor fusion, learning weights of each sensor modality on-line and in a real-time fashion may help to achieve better recognition accuracy and fast adaptation on new experimental settings.

Conclusions
In this paper, we perform an empirical study about the on-body sensor positioning and data acquisition details for HAR systems. We train the raw activity data collected in both real-world and controlled environments using a deep learning framework, which automatically learns features through neural networks without heuristic domain knowledge applied.
By leveraging the LSTM network that manipulates temporal dependencies on the time-series activity data, we identify that low sampling rate-as low as 10 Hz-is sufficient for the activity recognition. As low sampling frequency reduces the system burden by preserving battery and storage capacity, it consequentially allows prolonged data collection for HAR applications that usually operate on resource-hungry mobile devices. Our experiment result indicates that only two sensors attached on the right wrist and right ankle can present reasonable performance in recognizing ADLs including eating and driving (e.g., either as a driver or a passenger) activity. In other words, placing one sensor on the upper half and the other on the lower half of the body is recommended as the practical parametric setting of sensor position for further HAR research.
Additionally, we investigate the impact of sensor fusion by applying a two-level ensemble technique including stacking and voting on the multimodal sensor data, and demonstrate the performance improvement. By analyzing the recognition accuracy of each sensor on different types of activity, we elaborate custom weights for sensor modalities which can reflect the characteristic of individual activities.


A COMPLETE FRAMEWORK OPERATING SPATIALLY-ORIENTED RTI IN A 3D/2D CULTURAL HERITAGE DOCUMENTATION AND ANALYSIS TOOL

Close-Range Photogrammetry (CRP) and Reflectance Transformation Imaging (RTI) are two of the most used image-based techniques when documenting and analyzing Cultural Heritage (CH) objects. Nevertheless, their potential impact in supporting study and analysis of conservation status of CH assets is reduced as they remain mostly applied and analyzed separately. This is mostly because we miss easy-to-use tools for of a spatial registration of multimodal data and features for joint visualisation gaps. The aim of this paper is to describe a complete framework for an effective data fusion and to present a user friendly viewer enabling the joint visual analysis of 2D/3D data and RTI images. This contribution is framed by the on-going implementation of automatic multimodal registration (3D, 2D RGB and RTI) into a collaborative web platform (AIOLI) enabling the management of hybrid representations through an intuitive visualization framework and also supporting semantic enrichment through spatialized 2D/3D annotations.

INTRODUCTION: FILLING THE GAP BETWEEN CRP AND RTI
Multimodal acquisitions in Digital Heritage studies are becoming commonplace while requiring efficient frameworks to fully enable the analytic potentials of diverse spatially-registered results produced with imaging techniques. The photographic-based techniques, Close Range Photogrammetry (CRP) and Reflectance Transformation Imaging (RTI) (Ciortan et al., 2016) are widely used and diffused among CH community, also because several open-source packages are available to process and visualize those data types. Nevertheless, they are usually not used together since there is a lack of integrated solutions able to immerse different data types in the same working space and to visually analyze them jointly. The so-called multi-view RTI have been presented as a viable solution (Gunawardane et al., 2009), but both a proper methodology for acquisition, and accessible tools for processing and visualization are missing. Therefore this paper proposes a methodological shift, since the two media are not treated in isolation or side-by-side; conversely, we present the first results of a combined approach to link CRP and RTI into a common framework. Starting from some open source solutions, we designed a computational workflow allowing to merge and visualize RTI images with 3D models (as well as the raw 2D images used to produce the 3D models). Therefore, the main goal of this paper is to demonstrate the potential of an open platform integrating spatially-oriented and re-ligthable images into a 3D environment able to visualize data produced with the multi-view reconstruction approach.

RELATED WORK
CRP and RTI, are considered the dual low-cost techniques for 3D and 2D documentation in ICT and CH literature, broadly describ- * Corresponding author ing either the basic technologies and practical guidelines (Remondino et al., 2014, Ciortan et al., 2016 or the many applications tested (coins, rock painting/inscriptions, painted canvas, etc). Several works have also proposed comparisons to highlight advantages and weaknesses among competitive technologies (Mathys et al., 2013). On one hand, Multi-view and imagebased modeling (CRP and SfM) are periodically discussed regarding their gain of completeness and/or accuracy in necessary critical overviews (Remondino et al., 2014), whereas they can be exploited nowadays as a stable and reliable technique. On the other hand, Multi-light and RTI's techniques are massively developing from the original open-source code (Malzbender et al., 2001, Mudge et al., 2006, with numerous enhancements concerning the fitting-viewing side (Palma et al., 2010, Giachetti et al., 2017a and going to multispectral (Hanneken, 2014, Giachetti et al., 2017b and new devices (Schuster et al., 2014) support.
A third part of the literature focuses on the technical fusion of both approaches aiming or to enhance a 3D model with religthable reflectance texture mapping (Berrier et al., 2015) either to improve 3D reconstruction (Wu et al., 2011). Finally, a recent work tries to link RTI and 3D model from a methodological side (Shi et al., 2016), but a common and reliable spatial registration to fully benefit of this type of data fusion and endorse human-driven/datadriven analysis is missing. Meanwhile, the issue of multi-view RTI has been discussed for a decade (Gunawardane et al., 2009)  lenses and light sources) is a key feature for the quality of the final results. Despite their major difference is that one is based on the motion of a light source and the other on the motion of a camera, a common data acquisition protocol could easily converge to a combined approach, making easier the further data merging phase. Moreover, the camera parameters for optimal image quality for the two techniques are well-known as very similar even if the two processes are independent (e.g. lossless encoding, lowest ISO, sweet spot aperture, sharpness and fixed white balance).
To combine efficiently CRP and RTI, data acquisition protocol have been evaluated in wide diversity of acquisition scenarios, recommending to include or anticipate the RTI views as a component of the photogrammetric sequence; at the same time, we consider important to minimize as far as possible the contextual deviations of shooting conditions, parameters and equipment. In this context, the spatial registration of RTI will be conditioned by the robustness of incremental approaches relying on: i) features extraction and matching score; ii) the evaluation of intrinsic parameters; iii) the spatial overlapping allowing an accurate estimation of camera poses. The Spatial Orientation of an RTI towards the respective portion of 3D model -hereafter referred as SO-RTI -is maximized by the general overlapping and consistency of the RTI picture with its closer neighbours in the set of images acquired and processed in the CRP pipeline. From the photogrammetric meaning, this notion of overlap is extended to all elements composing the picture, including all elements below hereinafter cited as causes of matching issues : • The spatial resolution (i.e. the minimal distance on the object between two pixels) scaled by modifying the focal length and/or the distance to the object.
• The pixel resolution (i.e. the number of pixels count) set by the image dimension and related to the sensor pixel density.
• The radiometric values, more especially the radiance (hence the luminosity of pixel in the image) that will affect features detection and extraction.
• The temporal or contextual gaps, (i.e. differences on the object itself or on its close environment of acquisition) have also to be considered (e.g. highlight detection ball of H-RTI).
For this purpose a benchmark composed of tens data sets (see Table 1) from real CH case-studies and experiments led by the MAP laboratory for the past years have been re-processed and evaluated using the photogrammetric-based registration described in the next section.

TOWARDS FULLY AND VERSATILE AUTOMATED CO-REGISTRATION
Our co-registration method is based on incremental processing as image sets had to be merged among several iterations, with multiple users and/or at different time. As tie-points extraction and matching is the core and the main issue of photogrammetric based 2D/3D registration, our pipeline is initialized by an enhancement of 2D matching based on best neighbours pre-selection (see Figure 15) aiming to optimize the numbers and the quality features extracted. It is followed by a self-adaptive method for internal and external calibration fitting to the complexity of the dataset. A complete description of the "Totally Automated Co-Registration and Orientation" photogrammetric based workflow, called TACO is given in a complementary publication (Pamart et al., 2019).
The global co-registration processing dedicated to RTI is structured as follow (see Figure 1), a first iteration (Iter0) aims to achieve a first photogrammetric result. At the current stage of development, this initial 3D scene composed of calibrated/oriented pictures and a 3D point cloud is set as a master acquisition and stay frozen. Hence, it implies that the RTI must overlap with the area covered by CRP dataset initially processed (completion and/or densification of the scene is a planned future work). On next iterations (Iter1 and next), RTI master views are extracted, matched and spatialized, sequentially trying two methods. Matching is performed on the albedo image, produced by exploiting the inverse lighting capacities of RTI. It has been chosen, among other master images derived from RTI such as Mean or Median composites, since our tests produced better results in feature extraction and matching, even if we noted a potential issue with some radiometric deviation. Figure 1. Diagram of the iterative co-registration processing for SO-RTI scene The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W9, 2019 8th Intl. Workshop 3D-ARCH "3D Virtual Reconstruction and Visualization of Complex Architectures", 6-8 February 2019, Bergamo, Italy The spatial registration is attempted firstly using a conventional Bundle-Block Adjustment (BBA) on the RTI image from which metadata have been copied from raw images so as to re-use and/or refine according to EXIF, the internal calibration from previous iterations. However, performance and accuracy of BBA aiming to retrieve both internal and external calibration on a single picture could be sub-optimal or fail, according to aforementioned matching issues. In the case of rebellious RTI (i.e. the ones with lower overlapping) we switch to a second approach, Spatial Resection (SR) using Direct Linear Transformation (Shih and Faig, 1988), which requires fewer matching points. Unfortunately, SR needs Ground Control Points (i.e. 2D and 3D coordinates) not always available and not suitable for automatic processing. Thus, an approximate but automatic solution has been implemented. Our GCP's generation is based on 2D tie-points and corresponding 3D coordinates extracted from depth-maps computed in the initial iteration. Of course the performance of the results is limited by the precision of tie-points and the relative accuracy of the initial dense cloud. Nevertheless, this surrogated method shows robust results in term of finding an homographic relation which is sufficient for our visualization purposes (as shown in Figure 14).
Nonetheless, the merging of multimodal CH dataset composed of photogrammetric and RTI acquisitions can be achieved through variable scenarios according to the type of data acquisition and set-up used for the respective technique. The first and optimal scenario consist of an ad-hoc combined data acquisition where the RTI have been acquired with an automated dome, thus using the same equipment (sensor and focal) in a fixed context (see data sets, results and illustrations of Tegulae, Manuscript and Vasarely). However, the RTI-to-CRP registration produced accurate results, even when slight changes in illumination occurred. As shown in Figure 2 Table 2. Results of orientations on multimodal dataset legend • Partially or oddly oriented • Not oriented 5. AN HYBRID 2D/3D VISUALIZATION FRAMEWORK Any modern collaborative system for visual analysis shall be designed as a web platform so as to allow multiple collaborators to work at the same time, and where data should not be replicated locally. This has been shown to be feasible even in cases where the typology and amount of data was massive, and the timeframe very short (Apollonio et al., 2017).
In this paper the data were integrated in the context of AIOLI, a web platform that allows to navigate, visualize, interact and annotate all those oriented images and 3D data in a viewer composed of multilayered cameras, light sources and 3D points. The 3D reconstruction produced with CRP is stored as a point cloud and it acts as a medium to transfer the annotations made on a single image to all the images which are depicting the same portion of surface. To do this, once the registration and images processing are complete, an indexation step aim to store the 2D/3D projective relationship that links all the contents together in a single coherent 3D representation space (Manuel et al., 2016). Thus, for each image an array is created to store the indexes of the 3D points corresponding to the projected pixels. This important step is more generally used to allow the automatic replication of annotations made on a specific resource (whether 2D or 3D) on all the others, but also to re-project on images some geometric descriptors (Manuel et al., 2018). All these contents are associated with layers and semantic attributes that can be freely controlled by the user. AIOLI's viewer distinguishes original oriented images, tagged as Master, and different auxiliary images, derived from calculation methods (normal maps, depth maps, curvature maps, uncorrelated images, etc.). Master images are used to simulate virtual cameras copying their intrinsic and extrinsic parameters, all images being visualized by texture projection over the near plane of its camera's frustum. According to his needs, the user can then switch from one texture to another and use it to make his annotations. The integration of RTI for the navigation and interaction is nearly straightforward. The image can be visualized as a master image together with all the others, possibly choosing a clear representative image (see Figure 3).
When an image is selected for tracing and annotation purposes, the interface is essentially the same as the one used to visualize RGB images, with the addition of a button that enables the light direction variation, so that the best possible illumination can be chosen to draw or select a new annotation area (see Figure 4).

PROCESSING AND ENCODING RTI IMAGES FOR WEB VISUALIZATION
While the potentials of collaborative web tools are clear, additional challenges rise especially when complex data sets have to be handled by several people at the same time. An annotation project, possibly including multimodal data, can easily reach hundreds of items to handle. A main issue is dealing with geometric data, for which several solutions has been devised (Ponchio and , Mwalongo et al., 2016. For the handling of complex point clouds, AIOLI makes use of an octreebased rendering library (Schutz, 2014) which is able to deal even with billions of points.
RTI images encode a quantity of information which is much bigger than RGB images, but even high resolution images can be an issue: in the past some solutions were devised (Palma et al., 2014 to provide easy web publication and visualization of RTI and high resolution images. For the integration of RTI in Aioli, we have endorsed a very recent encoding of RTI data (Ponchio et al., 2018b) which is able to fulfill the usually hard task of offering high compression without sacrificing sampling accuracy. The adopted RTI representation combines Principal Component Analysis (PCA), adopted for data reduction, with an interpolation scheme. PCA is heavily used for the compression and the transmission of Bidirectional Texture Func-tions (BTF) (Schwartz et al., 2013). In this context, the idea is to apply it to produce a compact representation of all the images acquired (N ) in different lighting conditions. Then, these images are recombined using Radial Basis Function (RBF) to obtain high visual fidelity. These two operations can be easily integrated into a weighted linear summation, making the visualization client simple and efficient to implement. The quality/size ratio obtained in this way is always higher than the one of standard RTI techniques, such as PTM or HSH. Additionally, the appearance of highly specular material is well preserved at a cost of moderate file size (300-400 Kbytes).
All the above presented results have been implemented in the Relight library (http://vcg.isti.cnr.it/relight/), that includes not only tools for creating a RTI image in any format, but also a Javascript library for web visualization. This initial code was modified to be integrated within the AIOLI framework, so that images are encoded in a tiled format, and refined progressively as soon as bandwidth is available.

CURRENT STATUS AND RESULTS OF THE IMPLEMENTATION IN THE AIOLI PLATFORM
In addition to the methodological explanations above, this section intends to describe and clarify the technical aspects of SO-RTI stages and framework currently implemented in AIOLI. Thanks to its flexible architecture, we have been able to easily implement this RTI dedicated pipeline (presented in Figure 5). All of the processes are containerized thanks to Docker (Merkel, 2014) and are triggered by the server as every other AIOLI process pipeline. This allows us to benefit from the scalability, while using regular AIOLI communications to speak with the database and the project storing unit, thus not requiring to build a specific communication layer.  Figure 5. Diagram of the iterative co-registration process for SO-RTI scene The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W9, 2019 8th Intl. Workshop 3D-ARCH "3D Virtual Reconstruction and Visualization of Complex Architectures", 6-8 February 2019, Bergamo, Italy The framework currently requires the upload of raw RTI images and the corresponding LP files (which include calibrated light direction data), to allow to control and automatize the entire process. This includes the RTI encoding, which is implemented by avoiding a user-driven or explicit ball and highlight detection. To simplify the first implementation of the workflow, we focused our experiments on this kind of scenario. In the future, a metadata management system (i.e. group by EXIF) will help to define the best strategy for incremental processing using the optimum or combined BBA/SR approach. The main difficulty is the integration of unprocessed Highlight-Based RTI but thanks to the flexibility of the Relight library in term of supported fitters (Ponchio et al., 2018a), importing processed data in RTI or PTM format will be an easy task. For further information on updates, next features releases or conditions to join AIOLI's the beta-testing program visit the website (http:www.aioli.cloud/) The first results show a robust and versatile co-registration process but remains perfectible and still has to be tested on: i) reproducibility, complex dataset (scaling, spectral or temporal gaps and numbers of iterations); and ii) accuracy, checking by a necessary shift from relative to absolute coordinates with known or measured positions. Not surprisingly, the best results are obtained with ad-hoc combined acquisition or at least with fixed set-up contexts, diminishing global overlapping and matching issues. The 2D matching approach seems robust enough even considering some spatial or resolution gaps (e.g. distance, focal or sensor) solved by an automatic increase of image resolution to obtain optimize tie-points extraction used for both methods of spatial registration. The BBA (Figure 2) shows decent residual error (below 1px) but it is still perfectible in computation time. As the automatic dense image matching is performed in medium density (i.e. 1pt for 16px) the alignment is therefore accurate enough to satisfy visualization and annotation requirements. The results of the second SR method, applied in case of failure of BBA, shows variable and perfectible results. The latter shows a "depth inconsistency", meaning that focal and/or object distance is incorrectly approximated while the orientation point of view correspond to a consistent homographic projection (see Figure  6). As for BBA the correctness of SR registration is not related to number of keypoints used for resection but more their redundancy and spatial distribution. Obviously, fewer are the key-points and their multiplicity, higher will be the uncertainty of the calibration and orientation computed. Solutions to overpass this issue will focus our future works as other approaches (see Figure 7) and algorithms (Ke and Roumeliotis, 2017) could improve the spatial registration of rebellious resources (i.e. RTI or others).

CONCLUSION AND PERSPECTIVES
Multimodal acquisition (intended as the production of not only images and 3D data, but also of other advanced types of data like multispectral images, HDR, RTI, etc) is becoming a standard in the inspection and diagnostic analysis of works of art (being them small objects or an architectural environments) (Pamart et al., 2016). Nevertheless, there is a lack of tools that provide visualization and interaction modalities for this plurality of data. The presented work is a first step towards this objective, since an advanced type of acquisition (RTI images) is fully integrated with usual 2D and 3D data. The current implementation of SO-RTI already participate in blurring the lines between acquisition modalities. This paper presented a first effort in the direction of integration of different types of media in the context of a collaborative web platform. The integration has been tested on some real cases, showing how RTI can be easily integrated and help the analysis and annotation work by enabling the transfer of quantitative and qualitative properties. From this first outcome, there are two main research directions toward a so called high-level data fusion (Ramos and Remondino, 2015) that could start from the promising results presented in this paper: • The seamless integration of other multimodal data. In principle, any advanced image that could be encoded with an interpolation (i.e. time lapse images) could be directly included in the Relight library. Another challenge is represented by hyper-spectral images, which are becoming more and more important in multimodal acquisitions. In this case, similar issues to RTI have to be taken into account: coregistration with RGB images, data compression, interaction in the context of Aioli, and seamless remote navigation.
• A further direct contribution of RTI images in the registration and geometry calculation process. As already shown in previous works (Dellepiane et al., 2006), the estimated normals that can be extracted from RTI are coherent with the real shape of the object, but not accurate enough. If accurate calibration of acquisition device is performed (with a better estimation of light direction, uneven light coverage, white balance), moving more towards a photometric stereo approach, the extracted normals might be more accurate that the ones extracted by CRP. Hence, the point cloud can be "refined" in a second photogrammetric iteration (Giang et al., 2017), helping in having a more accurate representation and re-projection of annotations. Reversely, the RTI could also benefit of the photogrammetric fusion, from a metric point of view, with the foreseen orthorectified or orthomosaic RTI with our on-going experimentations (Fig. 13). The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W9, 2019 8th Intl. Workshop 3D-ARCH "3D Virtual Reconstruction and Visualization of Complex Architectures", 6-8 February 2019, Bergamo, Italy Convergent CRP with 50mm of Iter0, Close-up CRP with 60mm macro of Iter1 and two RTI with 24mm and 60mm registered with BBA Figure 9. Multi-focal CRP blocks and RTI views oriented with BBA on Uffizi dataset.
CRP and RTI both with 35mm lens Figure 10. Successful low spatial overlapping co-registration on low-textured Vasarely dataset.
CRP and RTI both with 60mm macro lens Figure 11. Correct SO-RTI achieved with ad-hoc acquisition of Manuscript dataset.
CRP with iPad and 2 RTI with Nikon D800E+60mm macro lens Figure 12. Consistent BBA based registration of two RTI acquisitions on low-resolution overlapping of iSkull dataset. Figure 13. Experiments on ortorectified RTI made on Uffizi case-study.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W9, 2019 8th Intl. Workshop 3D-ARCH "3D Virtual Reconstruction and Visualization of Complex Architectures", 6-8 February 2019, Bergamo, Italy Convergent CRP with 35mm of Iter0, Orthomosaic CRP with 60mm macro of Iter1, 5 next iterations for BBA and RS based RTI's registration (highlight areas