CC BY
12
The methodology for obtaining nonlinear and continuous three-dimensional topographic data using inertial and optical measuring instruments of unmanned
ground systems
P. Musa1, I. Purwanto1, D.A. Christie1, E.P. Wibowo1, R. Irawan 2 1 Faculty of Computer Science and Information Technology, Gunadarma University;
2 Faculty of Industrial Technology, Gunadarma University, Indonesia, Jl. Margonda Raya No.100, Pondok Cina, Kecamatan Beji, Kota Depok, Jawa Barat 16424
Abstract
Topography is the study of an area on the earth's surface. This term relates to the land's slope or contour, which is the interval of elevation differences between two adjacent and parallel contour lines. Topography generally presents a three-dimensional model of object surface relief and an identification of land or hilly areas based on horizontal coordinates such as latitude and longitude, and vertical position, namely elevation. The topography is essential information that must be provided in the execution of building or road construction based on the ground contour. The problem which is the ground contour which can provide visualization topography as a three-dimensional (3D) condition of the ground contour is not normal (non-linear). Another problem is that the traditional measurement techniques with wheel rotation only measure distances and cannot represent the trajectory of the ground contour in 3D. The proposed in-depth evaluation of orientation estimation results in the topography accuracy level. This methodology consists of several processes; Inertia and orientation of an object, Distance measurement, Terrestrial topocentric - Euclidean transformation, and Topography visualization. This research designed a prototype and proposed a new visualization method of the ground contours to reconstruct a topography map between three algorithms; Direct Cosine Matrix-3D Coordinate, Madgwick-3D Coordinate, and Complementary Filter. The methodology was tested and evaluated intensively by direct observation at three measurement locations with different difficulty levels. As a result, the Direct Cosine Matrix-3D Coordinate is able to visualize the ground contours by reconstructing a topography map much better than other methods.
Keywords: topography maps, visualization of a three-dimensional (3D), ground contour, inertia, distance measurement, direct cosine matrix, madgwick algorithms.
Citation: Musa P, Purwanto I, Christie DA, Wibowo EP, Irawan R. The methodology for obtaining nonlinear and continuous three-dimensional topographic data using inertial and optical measuring instruments of unmanned ground systems. Computer Optics 2022; 46(2): 280-297. DOI: 10.18287/2412-6179-CO-915.
Acknowledgments: The work was fully funded and supported by Gunadarma University, Indonesia.
Introduction
Computer science has been widely applied and implemented with existing technological advances in computer vision [1]. Computer science enables the function to visualize or see, identify the object, process images using data in the likewise process human vision does, and then deliver the appropriate output. One example of computer science image processses is creating and reconstructing 3D images for topographic maps.
A topographic map is a map that contains general information about the state of the land surface and information on its elevation using contour lines, namely the dividing lines, which are the locations of points with the same height to a particular surface reference of ground contour [2]. Contour is an imaginary line that is not straight but turns on the earth's surface, connecting points equal in height from sea level on the topographic map.
Meanwhile, the contour interval is the vertical distance between two height lines (the elevation differences between two adjacent and parallel contour lines) determined based on the scale [3]. A contour interval helps to build construction design planning [4]. The contour lines are sparse lines based on the linear estimated lines of the elevation changes, which are generally performed to fill in the gaps between two elevation points. In effect, the elevation information between the two contour lines is unknown. The process of identifying contour lines in individual subjects is a statistical analysis approach; adding straight points (linear lines) cannot be done.
Fig. 1 shows contour lines and contour intervals (elevation lines) on the map (the two-dimensional plane in the bottom image) and a land surface that has different ground contour (three-dimensional space in the top image). The results of the two contours form a winding and closed line based on the dots' results.
Contour lines
Fig. 1. The illustration of contour lines and contour interval on the topographic maps
A three-dimensional model makes it easier to read the contours of a place above the earth's surface because each height of line's height is immediately visible, rather than reading the two-dimensional model as shown in fig. 1 (bottom row). The research conducted by [5] to reconstruct a 3D requires data acquisition of multiple images from polarization experiments based on degrees and angles. However, the reconstruction of a 3D visual for a topographic map requires data acquisition of inertia and distance measurements. The input data in the form of a topographic map are analyzed and processed into a three-dimensional object model output to achieve this. The research proposed a solution to overcome these limitations by making a topographic map visually utilizing inertia data and distance measurements.
Designing ground contour is interval estimation (in the future, it will be referred to as ground contour estimation) which can be viewed from another perspective, namely the three-dimensional (3D) trajectory of a moving system estimation. Visualization in three dimensions (3D) by utilizing the orientation's attitude as the object's points from the scene will make to understand the algorithm from analyzing ground contour estimating two methods: DCM and Madgwick [6]. Toward visualization of the ground contour estimation in 3D, Dead reckoning and odometry are scanning techniques utilizing the counter of wheel rotation [7] and inertia sensor modules [8 - 10].
Measurement is an activity comparing a measured-quantity with a measuring instrument; everything in shape must have a size, length, height, weight, volume, or dimensions. The distance measuring tool is one of the civil works measuring tools in recording a road length, as shown in fig. 2.
The measurement technique used is to enumerate distance using a manual tool based on wheel rotation. Some distance measuring instruments have been designed to measure and record distances automatically by utilizing a rotary encoder sensor. A rotary encoder is an electromechanical device that can monitor and control movement and position [12]. A rotary encoder uses an optical sensor that produces a serial pulse that can be interpreted as information on movement, position, and direction.
The disadvantage of traditional measuring wheels is their inability to estimate 3D trajectories because the result of measurements as odometry is entirely dependent on the measurement of wheel rotation, which in turn the resulting trajectory is only in 2D.
Fig. 2. The tools of measuring the distance with rotation of the wheels [11]
Representing a topographic map is a complicated process, but this process can be simplified with the help of a computer graphics application. For example [12], suggests to make a topographic map (example, a vector format, raster format, or printed on paper) based on data sources INSPIRE (Infrastructure for Spatial Information in the European Community) [13]. Indonesia developed the National Geospatial Information Network (NIGN) intending to apply the principles of one standard, one reference, one database, and one geoportal in the management of Geospatial Information so that it can be used as a reference for spatial data, spatial planning, sustainable natural resource management, and policy formulation and decision making. There are 7 (seven) main aspects of NIGN [14 - 15]: regional and inter-state boundaries, forestry areas, spatial planning, infrastructure, licensing and land, natural resources and the environment, particular areas, and transmigration. For example, aspects related to our research are spatial planning, infrastructure related to civil works, and the need for topographic maps. The engineers in space geodesy in the MIIGAiK approve solution of the fundamental geodetic problem using tools and methods of space geodesy [16 - 17] and another researcher in the field of aerial photography about the terrain and survey of the object which uses drone [18].
Numerous researchers have explored and extracted geographic information from scanned topographic maps (STM). Their approaches are grouped into two categories which are segmentation and vectorization [19]. The other method is based on the appearance on topography maps from Google Earth which is converted into X, Y, Z coordinates. The coordinates can be mapped in the form of a three-dimensional map showing the hills' slope [20]. However, Google Maps being obtained from satellite image information in making topographic maps have weaknesses. One of the weaknesses is to predict the slope or the surface height of the land contour, especially in hilly areas where there are a lot of trees that make satellite images difficult to obtain valid values. The research [21] to reductions the cost and time of experimental investigation of the topography modeling using Gaussian and non-Gaussian of the highest probability distribution.
On the other hand, Dead-Reckoning can produce 3D trajectories, but since it depends entirely on accelerome-ter sensor measurements, the error accumulation occurs significantly [22, 23]. Odometry could perform a position measurement more accurate than dead reckoning because it combines an optical-based rotary encoder attached to the wheel, which has higher degree of accuracy in measuring distance. This research focuses on visualizing topographic maps which can be the accurate three-dimensional ground-contour measure that will get a land contour plot based on non-linear and continuous data.
1. Literature study 1.1. Topographic may
The topographic map is a view of land containing general information about the state of the land surface and information on its height using contour lines, namely the dividing line, which is the location of points with the same height to a specific reference plane (guideline / reference).
The altitude line aims to find out how high a place is from the sea level. The height line has the following characteristics:
a) The lower altitude line continuously circles the higher altitude line.
b) The height lines must not intersect each other and must not branch.
c) On sloping areas, the elevation lines are far apart. On the contrary, in steep areas, the elevation lines are close each other. For special regional conditions (such as cliffs, craters, ravines), the elevation line is mainly drawn.
d) The elevation line that juts out is a ridge and is always in the shape of the letter 'U'.
e) The inward indentation of the elevation line is a valley and is always like the shape of letter 'V'.
f) The height difference between two successive elevation lines (interval) is half thousand scales (example: 1 / 2000 x 50.000 = 25 meters). Unless otherwise stated.
g) The helper height line represents the height between two consecutive height lines.
h) The color of the elevation lines on the map is depicted in a specific color (usually brown).
1.2. Orientation representation
The orientation of a system can be described with several representations, including matrix representation and quaternion representation. Each has its own operating rules.
1.3. Rotation matrix
Rotation is a basic of motion of a particular space from geometry in physics. A rotation is the motion of a rigid body around a fixed point and preserves one point differently. Based on Euler's rotation theorem from geometry in 3D, a rigid body's displacement is equivalent
to a single rotation about some axis that runs through the fixed point [24, 25]. The theoretical and applied rotation matrices are discussed in rigid body mechanics, robotics, spacecraft attitude dynamics, navigation, three-dimensional, and computer graphics. The proven use of rotation matrices to realize rotations in R3 by [26 - 28] arranges the Euler parameters in quaternion form. They compute the double quaternion representation of rotations in four dimensions from 4 x 4 rotation matrices.
A rotation matrix is a matrix used in linear algebra to perform a rotation (roll, pitch, dan yaw) in Euclidean space. The rotation matrix used to represent an orientation in 3D space is a 3 x 3 matrix in matrix form as follows (1), where the notation describes the orientation of the b coordinate frame (body) to the navigation coordinates n.
R = Rb =
'11
'12
'13
r21 r22 r23
r31 r32 r33
(1)
The elements represente by rij, in which variable i (the row) and j (the column) are the angles between the i-axis in the navigation coordinates and the j-axis in the body coordinates. The relationship between vectors in body coordinates, rotation matrix, and vectors in navigation coordinates can be represented below.
V = RVh
(2)
where R is the rotation matrix, V„ is the vector of the navigation coordinates, and V„ is the vector of the body coordinates.
In Euler's method, the transformation from one coordinate form into another form can construct three successive rotations with different axes. The transformation can be from one coordinate form into a new coordinate form which can be expressed as follows:
Rotation through the x-axis of the body is called the roll angle (®). The rotation angle of the roll on the x-axis uses the Ri matrix with equation (3).
R1 =
1 0 0 0 cos ^ - sin^ 0 sin ^ cos ^
(3)
The rotation through the y-axis of the body is called the pitch angle (0). The representation for the rotation of the pitch angle on the y-axis uses the R2 matrix with equation (4).
R2 =
cos 8 0 sin 8 0 1 0 - sin 8 0 cos 8
(4)
Rotation through the body's z-axis is called the yaw angle (y). Moreover, the rotation of the yaw angle on the z-axis uses the R3 matrix with equation (5).
R3 =
cos y - sin y 0 sin y cos y 0 0 0 1
(5)
From equation (3) to (5), the matrix does multiplication - to produce a 3 x 3 rotation matrix (R) as in equation (6) below.
R =
cos 9 cos y sin ^ sin 9 cos y-cos ^ sin y cos ^ sin 9 cos y + sin ^ sin y cos 9 sin y sin ^ sin 9 sin y + cos ^ cos y cos ^ sin 9 sin y-sin ^ scos y - sin y sin ^ cos 9 cos ^ cos 9
(6)
Numerical errors stemming from the quantization of values caused by sensors and insufficient time resolution cause the rotation matrix to lose its properties over time. As a result, the rotational transformation performed by the rotation matrix is no longer rigid body transformation. The rotation matrix has unique properties, which are orthogonal and normal.
c = a x b
e/2 d/2
Fig. 3. Illustration of orthogonal and normalization process
1.4. Quaternions
A quaternion concept is a complex number with four components: one real part and three imaginary parts. Those components are a four-dimensional (4D) number system that extends the definition of complex two-dimensional (2D) numbers. According to [6], the quaternions perform a rotation (likewise rotation matrix) in Euclidean space, but they have different rules of operations and properties. Therefore, a 4D vector should represent quaternions with a length of 1. When creating a number system capable of representing rotation in three-dimensional space, William Hamilton proposed a new number system, namely a four-dimensional system; one actual number (R) and three imaginary numbers i, j, and k, called Quaternions [29 - 31]. Due to the property of quaternions, each quaternion unit represents the rotation of an object in three dimensions (3D) [24], [32 - 36]. The quaternion formula has the following equation (7).
Or
q = qi + qii + q^i + q4 k
q = w + xi + yj + zk,
(7)
where #1, q2, #3, #4 e M are real numbers, and i, j, k are imaginary numbers that satisfy i2 =j2 = k2 = ijk=-1. Also, the real numbers are w, x, y, z or as a 4D vector q = [w, v] which is called scalar part and v = (x, y, z) which is called vector part. It is assumed that all real numbers are quater-
nions with q2 = q3 = q4 = 0 and all complex numbers are quaternions with q3 = q4 = 0. It is impossible to represent four-dimensional space because perception is only limited in three-dimensional space, so quaternions cannot be illustrated in a Cartesian diagram. The properties of the quaternions possess are [37]:
Rotations: The method of quaternions has proved rotations in 3-dimensional space-the different methods of traditional rotations around the axes of the Cartesian coordinates by angles y, 9, and 0. However, Euler in the general rotation of a rigid object can describe as a single rotation about some fixed vector. Given v = [l, m, n] over V3, then a rotation by an angle 0 about v is given by.
Lq(v) = qvq *, where
cos9 , sin9
q =
(8)
m-
sin9 sin9
Computer graphics: The application of quaternions can be realistic animation. This technique is called spherical linear interpolation (SLERP) and uses all unit quaternions to form a unit sphere. By representing the quaternions of key frames as points on the unit sphere, a SLERP defines the intermediate sequence of rotations as a path along the great circle between the two points on the sphere. Physics: In physics, the application of quaternions has been used in a wide variety of research, such as the making work on Special and General Relativity, Newtonian Mechanics, scattering experiments (crystallography and quantum mechanics).
The related research applies the quaternion method to spacecraft, such as the research of relative navigation, position, and attitude in a unified form of a spacecraft using dual quaternion-based relative with Kalman filter [38 -40]. The simulated spacecraft to identify and target asteroid is Kleopatra with a dual quaternion relative navigation filter [41]. In four-dimensional, it proposes the quaternion equations, the spatial inertial navigation systems for an azimuthally stabilized and a gyrostabilized keep orientation invariant in inertial space [42].
The advantage of the quaternion that other number systems do not have is its ability to represent rotations in three-dimensional space. However, the drawback is that the quaternion loses its commutative property, #1 #2 ^ #2 #1. This property can logically arise if associated with a three-dimensional rotation property, producing a different output of the order of rotation is also different.
1.5. IMU sensor and calibration
Inertial Measurement Unit (IMU) is an electronic device that is able to measure and report the direction, angular velocity, and magnetic field imposed on objects, using a combination of several sensors, including an accel-erometer, gyroscope, and magnetometer. Each sensor takes measurements in three axes (x, y, z) in the body frame. This three-axis measurement can be used to find the orientation of the objects represented in Euler angles, rotation matrix, and quaternions [43]. The IMU is part of the navigation system known as the Inertial Navigation System (INS). The most commonly used representation is the Euler angles in the three principal axes; x for roll, y for pitch, z for yaw. Pitch is a change in position (e.g., an airplane) in the direction of the flight path but causes the aircraft to climb or descend. Roll is a change in position (e.g., an airplane) in the direction of the flight path but causes the rotation of the aircraft body or the aircraft to be tilted. Yaw is a change in position (e.g., an airplane) in the direction of flight and a change in position with a combination of roll or pitch. The IMU does use on aircraft, Unmanned Aerial Vehicles (UAV), spacecraft, satellites, etc. The accelerometer is a sensor used to measure acceleration, detect and measure vibration, and measure the acceleration due to gravity. Gyroscope sensor is to measure the angular velocity of an object. In comparison, the magnetometer sensor does use to measure the direction that comes from the magnet's strength or geomagnetic. The advantage of the IMU system is that it calculates continuously from the orientation of an object, known as the Attitude and Heading Reference System (AHRS) [44].
IMU sensor calibration is performed with two approaches; offset calibration [45] for accelerometer and gyro and hard-soft iron calibration [46] for the magnetometer. Offset calibration is conducted by finding the offset value, xo, which is the mean difference between the measured value, x, and the value of the sensor should be measured, x, expressed by
Zn K
i=X - X .
Hard-soft iron calibration is performed by estimating the transformation matrix M and translation, t, where the value measured by sensor, p, is transformed to actual value p; expressed by p = M (p -1 )p = M (p -1)
1.6. Orientation measurement
Naive measurement
Naive measurement is meant to be an orientation measurement without control / correction. This measurement involves only two sensors, an accelerometer, and a gyroscope. Both sensors measure orientation in Euler, ea = ((a, 9a), eg = ((g, 9g, ^g), according to their respective measurement. ea and eg are then subjected to low pass and high pass filters, respectively [47], resulting in a system orientation estimation.
Direct cosine matrix
The details of the DCM orientation estimation algorithm can be found in the original publication [48, 49]. Generally, the DCM consists of 4 stages: (a) Rotation matrix update, (b) Ortho-normalization, (c) accelerometer and magnetometer drift detection, (d) drift correction with PI-controller.
Madgwick
The details of the Madgwick's orientation estimation algorithm can be found in the original publication [50]. Generally, the Madgwick's algorithm consists of 4 stages: orientation in quaternion calculation, gradient error from accelerometer and magnetometer calculation, gradient descent, integration and normalization to the quaternion.
1.7. Distance measurement
This system may include a rotary encoder, a management that measures the distance of contour lines and accumulates it with distance measured from the rotary encoder by the robot. The rotary response conversion to the distance units is highly dependent on the characteristics of the wheels used by the system, such that,
s = n X
2%r
(9)
The measured distance in units is used by r, n is the number of pulses measured, c is the number of pulses given by the sensor in one revolution, and r is the system's wheel radius. [51] expanded the odometry method with PID in order to correct a position and determine position.
1.8. Terrestrial topocentric to euclidean coordinate system transformation
Terrestrial topocentric systems define a coordinate system in the plane tangential to the earth's physical surface, as follows [52]: (a) at least the point of origin is at a position near the surface of the earth, (b) the primary of a ground contour is the plane related to the earth's surface at the point, (c) the elemental axis is from north point, (d ) the explain topocentric systems are usually left-handed. Planar A
Fig. 4. Geometric illustration for terrestrial topocentric-to-euclidean coordinate system transformation calculations
Consider a point p e V3 at the origin (fig. 4). A vector v (V3) describes the motion from p to p' in 3D space. Vector v can be projected into two planar planes A and B, with the origin at p. Consider plane A, which is represented by a normal vector nA, placed in V3, such that v • nA = 0 and nA • nB = 0, where plane B is parallel to the ground and all notations are represented in the world frame reference.
In-plane A coordinate system, vector v can be referred to as ae V2 and a • nA = 0. Vector v projected into plane B can also be referred to as beV2 and b • nB = 0 . Vector a is obtained using pitch measurement, and vector b is obtained using yaw measurement.
Elevation transformation to Z-Axis
Consider vector a as a linear change in height over a given time interval in the plane A coordinate system. The current height is obtained by projecting vector a onto y A in plane A in V2, and because y is colinear with Z in V3, then.
yA = Z„ = J|a,| sin e,. (10)
t=1
Where Z„ is the current height, | at |= | vt | = st, is the current distance obtained from the rotary encoder, et is the pitch value obtained from the orientation estimation algorithm, and t indicates a specific time interval. yA can be viewed as the Zn axis of the system's position. Azimuth tranformation to XY-axis
Consider vector be V2on plane B, which is the projection of vector beV3. Vector b can be projected onto the x and y axes, on the plane B in V2. Since x and y are colinear with the X and Y axes in V3 and it is assumed the linearity of the function when t ^ 0 (as used in Elevation Transformation to Z-Axis), then.
xBn = Xn =£ |at| cos (et )cos (yt), (11)
t=1
n
yH = Yn = at | sin (et )cos (y t). (12)
t=1
Where | at |= | vt | = st is the current distance obtained from the rotary encoder, et and y are the pitch, yaw values are obtained from the orientation estimation algorithm, and t indicates a specific time interval. We can view xnB and ynB as the Xn and Yn axis of the system's position.
1.9. Geodesy and cartography
Geodesy is a branch of applied mathematics that aims to determine the shape and size of the Earth, the coordinate position, length, and direction of the Earth's surface lines on a particular scale. Along with technological advances, geodesy science, which describes the position (global navigation satellite system (GNSS) by [53]) on the Earth using latitude, longitude, and altitude, is called Cartography [54]. Cartography is the art, science, and technology of making maps related to
geography, spatial information, and topography. The result of research by [55] represents a city guide tour application utilizing the concept of augmented reality in the area of geodesy and cartography.
1.10. Degree of _freedom (DOF)
The degree of freedom (DOF) is a formal description of the state parameters of the system concerning the degree of independence that states the position in the form of coordinates of the system. There are two degrees of freedom models: single degrees of freedom (SDOF) and multiple degrees of freedom (MDOF). According to
[56], SDOF is an object movement system that only has one coordinate axis to move in one direction. The MDOF is a system that has n-degrees of freedom of an object
[57]. For example, the MDOF case with three rotations (x, y, z) provides a displacement moment on the x, y, and z directions.
1.11. Visualization
Odometry visualization can be viewed as a 3D positional trajectory. Mathematically, the trajectory T is a set of 3D coordinate points such that T = {%1,..., Xn} where X e V3 and n represent time step.
2. Methodology 2.1. General system description
The methodology to obtain information on an area in making a topographic map by direct observations makes the study more accurate. It provides the result in three-dimensional form, distance, height, and angle by utilizing various tools or instruments.
The visual requirement of the three-dimensional modeling is the input required on the map topography in general, shown in fig. 5, which describes the system design in detail so that the required data domains, functions, processes, or procedures can be determined along with their performance and interface.
The block diagram describes the system as a whole, which the system is divided into two parts:
a) The device, hardware consists of a microcontroller and sensor modules such as an Inertial Measurement Unit (IMU) and Rotary Encoder.
b) Software, using tools the numerical computing to proposed methodology is implemented.
2.2. Proposed methodology
The objectives can be formulated into the following points:
a) Creating a new instrument that generates nonlinear and continuous plot contour charts to obtain information on an area in making a topographic map.
b) Proposing a new method adopting odometry and dead reckoning by combining the orientation measurement from the inertia sensor and the dis-
Fig. 5. Overall System Description.
Accelerometer
Gyroscope
Magnetometer
Acceleration (x, y. 7) axis
Angular velocity (x, y, z) axis
Magnetic Field (x, y, z) axis
C
s '-5
es
a
ii 'l 0
S
£
®
M
Rotary Fncoder
Distance Measurement
Topoeentric-to-
Euclidean
Transformations
Transformation Elevation ~ '
Pitch Length
Transformation Azimuth I í
Yaw
Length '
Roll
Odometry
z Visualization
X
Visualization
Y
Fig. 6. Proposed methodology
The proposed methodology, namely DCM-3DC (Direction Cosine Matrix - 3D Coordinate), is described by the block diagram in fig. 6. There are four processes (blocks), namely:
a) Orientation Estimation Algorithm processes IMU data into an orientation using an inertia module sensor.
b) Distance Measurement processes encoder data into a distance using a rotary encoder sensor.
c) Terrestrial Topocentric-to-Euclidean Coordinate System Transformation processes orientation and distance range data into a 3D coordinate system.
d) Odometry visualization is the process of interpreting the 3D trajectory of the system into a Topography Map.
2.3. Proposed methodology stages Orientation estimation algorithm
Direct Cosine Matrix is used as the orientation estimation algorithm in this study. Comparison between two popular algorithms (DCM and Madgwick's) has been intensively conducted in [6]. The results show that DCM
can estimate orientation more accurately in the situation with a lot of noise.
Ground contour estimation algorithm
To perform the ground contour estimation algorithm, the conditions that must be met are the data availability relative to the ground in Euler angles (9, 0, and y) and the displacement of the system, ds.
Algorithm 1 Topocentric to Euclidean Coord. System Transformation
Input:
Collection of orientation in each time step in euler angle E = {ei,..., e„}\ ei = {(, 9, y} e SO (3), Collection of distance changes in each time step s = {dsi,..., ds„} \ ds„ e R
Output:
Collection of 3D coordinate, X = {Xi,...,Xn} \X, = {x,y,z} e R3
Proses:
1. x ^ 0, y ^ 0, z ^ 0
2. For ei e E do
z ^ equation (10) x ^ equation (11) y ^ equation (12)
Ground plane estimation algorithm
Assuming the slope of the ground, which is parallel to the X-axis, is a linear plane where t ^ 0, then the plane at that specific point can be represented by two points at 8 distance from the origin of the coordinates of the system (system's center point) in both X-axes direction.
Algorithm 2 Ground Plane Estimation around System Estimation
Input:
Collection of point pairs in both direction of the X-axis of the system P = {. ,,p„} |pt = {(8,0), (- 8,0)}, 8 << 1, Collection of 3D coordinates,
X = {X1,., Xn} |X- = {x, y, z} e M3 Output:
Collection of point pairs that have been transformed according to the orientation of the X-axis of the system P' = {p1',.,p„'} |p' = {(8, 0), (- 8, 0)}, 8 << 1 Proses:
1. Rxz Extract the rotation of the X-axis and Z-axis of the system
2. Forp,X eP, X do
pi ^ Rxz x pi + Xi (rigid body transformation)
2.4. Testing and evaluation _ framework Testing with direct observation
Several locations were selected to perform direct observation. Three locations were selected, namely;
a) Simple ramps (at house terrace),
b) Ramps and junction (indoor area of office at Gunadarma University),
c) Double ramps and turns (outdoor area of office at Gunadarma University).
The locations represent three different difficulty levels, from the least challenging, simple ramps, to the most challenging, double ramps and turns (fig. 7, upper rows).
The three-dimensional trajectory is created based on the sketch numerical computing tools (see fig. 7, bottom row). The contour lines based on the topography are shown in fig. 7 (bottom row) for each observation at three locations, while the contours of the intervals (for example, observe the double ramps and turns) are shown in fig. 8.
Three methods are being compared;
a) The proposed methodology uses DCM as its orientation estimation, which will then be referred to as DCM-3DC (DCM 3D Coordinate),
b) The comparison method uses Madgwick's as its orientation estimation, which will then be referred to as Madgwick-3DC (Madgwick 3D Coordinate),
c) The baseline method uses fundamental 3D estimation, Complementary Filter. This third method is the primary method used as a reference for comparison with other algorithms.
a) b) c)
Fig. 7. Observation in three different locations. The upper row is the experiment place; the bottom row is the contour lines
(ground truth) built based on the actual topography
Qualitative testing
Qualitative testing is conducted by measuring the deviation of the estimated ground contour with the ground truth. The deviation can be viewed as a Euclidean distance of each measurement sample to the ground truth. The exploration of two sample points of the nearest neighbor which can estimate the ground truth [58] is performed (see fig. 9).
Four metrics are used to evaluate the performance of the associated algorithm, including:
Fig. 8. The contour interval on maps of topography Average euclidean distance (AED)
AED
- n SIx - x
(13)
Fig. 9. Pair of estimated sample points (blue dot with green outline) with the closest point to its ground truth (red dot)
Median euclidean distance (MED)
MED - medln-1
x - X
Mean error growth
Mean Error Growth = ^ ",=2 ( - e!-1).
Median error growth
Median Error Growth = medn= 2 ( - et-1).
(14)
(15)
(16)
Qualitative evaluation
Qualitative evaluation is conducted by analyzing the visualization results of ground contour and topographic map. The estimated visualization of the ground contours and the topographic map is collocated with the ground truth visualization in the same graph. Deviation will be shown if the algorithm is not well capable of estimating the contour. The measurement of error contains a quantitative analysis that reflects the advantages of the developed algorithm upon others algorithms using such metrics as mean difference, standard deviation, standard error, and measurement error at some (90 %, 95 %, or 99) confidence level. The formula of metrics is used for quantitative analysis, namely:
Mean difference
The measurement error of distance in the unit of
measure cm
x - — V Averageu de D t + Median, de ri.an .
/-k ^^^ <J Euclidean Dis tan ce Euclidean Dis tan ce
(17)
The measurement error of distance in the unit of measure cm / sampe
x - — V Mean „ G h + Median
/-k ^^^ ErrorGrowth
Standard deviation
s -a-
n-1 V (x-x )2.
(18)
(19)
Standard error
a
sx =-T • Vn
Standard error
Number of Parts Percentage =---x 100%.
I (x, )
(20)
(21)
Confidence level
The confidence level is equal to (1-alpha)* 100 %. If it is assumed that alpha equals to 0.05, it indicates a confidence level of 95 percent and the result of calculating the area on the standard normal curve is equal to (1 - alpha), or 95 percent is ± 1.96.
C
Confidence -1.96 (('yjnumber of samples V (22)
A confidence interval (or interval estimate) is an interval of values used to estimate the precise value of a population parameter. Therefore, the confidence interval is as follows:
x ± 1.96(aj'yjnumber of samples).
3. Results analysis and discussion
3.1. Developed device
(23)
Our research project built a prototype that could generate continuous and non-linear ground contour plot data to create topographic maps in 3D space between the two orientation estimation algorithms, DCM and Madgwick. fig. 10 is the measurement device that had been built in this study. It contains a microcontroller, two sensors (a rotary encoder and IMU), and an SD-Card. The reason is to make the robot mobile (next it is called unmanned ground systems (UGS) to carry out testing and survey on simple ground contours such as Simple Ramps, Ramps, and Junction and Double Ramps and Turns.
This research developed an instrument in the UGS as a measurement of wheel rotation (as shown in fig. 2) using an electrical encoder rotating component (fig. 10a) packaged in a single pack on the wheel drive, namely a DC motor. The rotary encoder sensor is combined with an inertia sensor (fig. 10b) to obtain the estimated orientations on 9 DoF from the Accelerometers, Gyroscopes, and Magnetometer.
To optimize the result, when the data are collected from the inertia sensor, this sensor is placed in the center of the UGS. The results from both sensors are recorded on the memory card during the scanning process against a surface area.
3.2. Sensor acquisition and calibration results
The parameter configured on the accelerometer is a known measurement range of ± 2 g, which is in units of gravity (g) as in the following equation:
i-1
1g = 9.8 m
(24)
The system will not experience extreme acceleration (above 2 g) because of the implications for measurements with the highest sensitivity. In contrast, the resolution is full-resolution, 10 bits. The two configured parameters (measurement range and resolution) will have implications for measurement performance, including sensitivity and scale. The smaller the measuring range, the higher the sensitivity. Furthermore, the ability to measure acceleration is accurate, as much as 256 Least Significant Bit (LSB) to describe 1 g or as much as 4.3 mg (milli-gravity) for every LSB.
a)
b)
Fig. 10. The developed ground contour measurement device
a)
The parameters of the gyroscope are similar to the accelerometer; and the configured parameters are measuring range and resolution. However, the size of the selected measurement range is ± 2000 dps (degrees per second), so the system is able to perform measurements at an angular velocity up to quite extreme conditions, namely 2000 dps. In addition, in contrast to the accelerometer measurement with a resolution of 10 bits, the resolution of standard configuration gyroscope measurements performed have 16 bits. The two configured parameters will have implications on measurement performance, including sensitivity and scale factor. For example, the ability of the system to measure angular velocity is relatively accurate, which is an average of 14 LSB to describe a change in angular velocity of 1 dps, or 69.5 mdps (milli-degree per second) for every increase in LSB.
Fig. 11 is the accelerometer and gyroscope offset calibration results. The red, green, blue, yellow, cyan, and magenta are the X, Y, Z, and calibrated X-Y-Z, axes, respectively.
On the magnetometer, the parameters that are configured are the measurement range and the resolution as well. The chosen measurement range configuration is the standard configuration, which is ± 1.3 G (Gauss - unit of magnetic field, 1 G = 10 -4 T, which T denotes teslas). That value of standard configuration was chosen because the average value of the Earth's magnetic field above the surface falls within the measurement range, around 0.5 Gauss and the configuration parameter at the resolution is 12 bits. The two configured parameters have implications on measurement performance, including sensitivity and scale factor. For example, the system's ability to measure the Earth's magnetic field is relatively accurate, i.e., 1090 LSB to describe a 1 Gauss or 0.92 mGauss (milli-Gauss) change for every LSB increase.
2 5
15
Calibrated gyroscope data acquisitions
> 0.5 -
0.5
_JLAJLT\*nW YVi YIT^W—T~
X
Y Z
X Calibrated
Y Calibrated Z Calibrated
*-\JWYJ—ITVTttV-ir
**j> UL.,vv««.'>w\m'YWY'/v'» " xnrnmrtiw \/-
» rWYV V~~
"■onn1"
b)
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (second)
Fig. 11. Accelerometer and gyroscopes offset calibration results
The plot of calibrated collection data to create a 3D scatter calibration is performed by rotating the IMU along
the X-axis, Y-axis, and Z-axis. Magnetometer captures magnetic field value in all possible rotation axis [43].
2
a)
0.4 0.3 -0.2 0.1 -0 --0.1 -02 --0 3 --0 4 --0.5 -
b)
V ll .-4 M »
t-•■■ ■■ . . r■
-05 -0.4 3D view
Fig. 12. Magnetometer measurements before (red) and after calibration (green)
Fig. 12 is a validation of the magnetometer calibration results. Before calibration (red), it shows an ellipsoid shape, while after calibration (green), it shows a spherical shape, which means that other magnetic fields that do not originate from the earth have compensated.
The purpose of the calibrated magnetometer is to be able to distinguish the translation results of the two ellipsoids in the Centroid to Origin Distance line. Before calibration, the distance from the midpoint of the ellipsoid to the origin is 0.077 G, and after calibration, the value decreases to 0.045 G. It means that the ellipsoid is closer to the origin. In addition, the difference in the shape of the two ellipsoids can be seen in the Standard Deviation of the Distance Points to the Centroid; the standard deviation before calibration is 0.02 G, while the standard deviation after calibration is 0.007 G; The closer to zero the standard deviation value is, the more spherical it is. It indicates that the ellipsoid after calibration is more spherical than before calibration.
3.3. Orientation estimation results
The comparison between two popular algorithms has been conducted intensively before this study [6], which is part of this comprehensive research. It shows that DCM could estimate orientation more accurately in a situation with lots of noise.
3.4. Ground plane and contour estimation results to making a topographic map
Data preparation
Sensor calibration and orientation estimation have provided a satisfactory performance; hence the ground plane and contour estimation can be conducted. When the ground plane with a contour estimation algorithm is evaluated, the system gets observational data as data collected.
The observation data acquisition format consists of 11 fields; time interval (ms), change in the distance (cm), three-axis accelerometer measurements, three-axis gyro-
scope measurements, and three-axis magnetometer measurements. For the last nine fields, the IMU sensor data, the raw data format (ADC sensor) are used. This method is useful toward reducing the process in the microcontroller, hence increasing system responsiveness.
1D smooth filter with time window applied to the 9 DoF IMU data is used to reduce the noise. Filtration is proven to be effective in increasing the performance of subsequent processes. Analysis
All data are used as the inputs for the ground contour estimation algorithm (Algorithm 1) and the ground plane estimation algorithm (Algorithm 2); we call it DCM-3DC. As a comparison, we built another algorithm using Madgwick's and Complementary filter as its orientation estimation.
Tab. 1 shows the metric results of three methods at three different locations. quantitatively, DCM-3DC provides significantly better performance than two others. Italways closes to its ground truth, in contrast to the other two deviations from its ground truth for all of the performance measures used in this study.
The quantitative measurement results shown in Tab. 2 provide evidence that the DCM-3D method is excellent. Based on the quantitative analysis results (used equations from equation 17 until equation 20), experiments of simple ramps with the DCM-3DC method result in a mean of 12.82 cm. The mean using the Madgwick-3DC method was 33.93 cm; and the Complementary Filter method was 29.72 cm in the simple ramps experiments. Therefore, the two experimental results provide significantly higher mean measurement result than the measurement result of the DCM-3DC method. It can be concluded that the DCM-3DC method is recommended for this experiment.
The quantitative analysis in the experiment of simple ramps by measuring the standard deviation shows that the measurements on the Madgwick method are 1.28 cm, and the Complementary Filter is 0.58 cm. In contrast to DCM-3DC, where the measurement result of the standard
deviation is smaller than the two previous methods, the standard deviation is 0.06 cm. Accordingly, the standard
deviation results in the DCM-3DC method are recommended to measure the simple ramps cases.
Tab. 1. Quantitative Evaluation of three algorithms at three different direct observation locations (the smaller, the better)
Number of Samples Distance (cm) Average Median Mean Error Median Error
Observation Method Distance Error Distance Error Growth Growth
(cm) (cm) (cm/sample) (cm/sample)
DCM-3DC 1.2861 1.2770 0.0005 0.0002
Simple Ramps 912 300 Madgwick-3DC 34.834 33.018 0.0028 0.0026
Comp. Filter 30.126 29.306 0.0024 0.0021
Ramps and Junction DCM-3DC 4.9086 4.3029 0.0012 0.0006
1312 1400 Madgwick-3DC 15.511 8.7419 0.0067 0.0048
Comp. Filter 7.5670 6.0253 0.0027 0.0031
Double DCM-3DC 13.3776 14.0543 0.0016 0.0014
Ramps and 2863 2200 Madgwick-3DC 25.2415 17.1909 0.0076 0.0067
Turns Comp. Filter 35.8518 33.0592 0.0077 0.0079
The subsequent quantitative analysis measures the standard error in the case of the simple ramp against the three methods that have been carried out in the previous measurement. The standard error measurement results for quantitative analysis are shown in table 2, and the DCM-3DC method gets the minor standard error of 0.05 cm.
Meanwhile, the results of other standard error measurements show that the standard error is above 0.05 cm; namely, the Madgwick method is 0.91 cm, and the Complementary Filter method is 0.41 cm. Furthermore, the quantitative analysis recommends that DCM-3DC is a better method for the simple ramp.
Tab. 2. Quantitative analysis of the mean, standard deviation, standard error
Mean Standard Deviation Standard Error Average
cm cm/sample % cm cm / sample % cm cm/sample % (%)
DCM-3DC 12.82 0.0004 16.76 0.06 0.0002 3.34 0.05 6.38 3.34 7.81
Simple Ramps Madgwick-3DC 33.93 0.0027 44.37 1.28 0.0001 66.59 0.91 16.51 66.59 59.19
Comp. Filter 29.72 0.0023 38.87 0.58 0.0002 30.07 0.41 14.65 30.07 33.00
Ramps DCM-3DC 46.06 0.0009 3.71 0.43 0.0003 0.69 0.30 2.15 0.69 1.70
and Madgwick-3DC 51.47 0.0058 41.49 50.85 0.0010 81.78 35.95 43.71 81.78 68.35
Junction Comp. Filter 67.96 0.0029 54.79 10.90 0.0002 17.53 7.71 30.13 17.53 29.95
Double DCM-3DC 137.16 0.0015 19.77 4.78 0.0001 5.87 3.38 70.27 5.87 10.51
Ramps and Madgwick-3DC 212.16 0.0072 30.58 56.93 0.0005 69.88 40.25 85.95 69.88 56.78
Turns Comp. Filter 344.55 0.0078 49.66 19.75 0.0001 24.24 13.96 165.29 24.24 32.71
The research process for quantitative analysis in the ramps and junction case, three measurement results, namely the mean, standard deviation, and standard error, are shown in tab. 2. The measurement results using DCM-3DC show that the mean value is 46.06 cm, which is smaller than Madgwick's measurement's mean value (51.47 cm) and Complementary Filter's (67.96 cm). In the standard deviation measurement in the ramps and junction, the DCM-3DC method still got the smallest value based on the results of the quantitative analysis. The values obtained from the standard deviation measurements are as follows: DCM-3DC is 0.43 cm; Madgwick is 50.85 cm; Complementary Filter is 10.90 cm. The standard error measurement result representing DCM-3DC is 0.30 cm. The smallest values based on the results of quantitative analysis of Madgwick and Complementary Filter are 35.95 cm, and are 7.71 cm, respectively. Thus, from the result of the experiment with three methods, the DCM-3DC produces the smallest value compared to the other two experiments; DCM-3DC is the best method for ramps and junctions.
The experiment at the third location is on double ramps and turns, where the quantitative analysis uses the DCM-3DC, Madgwick, and Complementary Filter methods. The three methods have performed a quantitative analysis to measure the mean, standard deviation, and standard error. The mean to quantitative analysis measurement results for double ramps and turns in the DCM-3DC, Madgwick and Complementary Filter methods are 137.16 cm 212.16 cm, and 344.55 cm, respectively. The smallest mean value was obtained using the DCM-3DC method. The quantitative analysis to obtain the standard deviation results in that the DCM-3DC method is 4.78 cm; Madgwick is 56.93 cm; the Complementary Filter is 19.75 cm. While the results of the quantitative analysis to calculate the standard error obtains that the DCM-3DC method is 3.38 cm; Madgwick is 40.25 cm; the Complementary Filter is 13.96 cm. Therefore, in the quantitative analysis process for testing the case of double ramps and turns, the mean calculation, standard deviation, and standard error of each measurement provides the smallest value by using the DCM-3DC method com-
pared to the Madgwick method and Complementary Filter method. Accordingly, the DCM-3DC can be recommended as the best method for double ramps and turns.
The average results by equation (21) in percentage unit's show that the average value of DCM-3DC is the
lowest; and it is better algorithm than that of the others. The average value of error of DCM-3DC at the simple ramps location is 7.81 %, while at the ramps and junction location, the average of error is 1.70 %, and at the double ramps and turns location the average of error is 10.51 %.
Tab. 3. Quantitative analysis of the observational / measurement error
Observation
Samples
Simple Ramps
Ramps and Junction
Double Ramps and Turns
912
1312
2863
idencelevel (%) DCM-3DC Madgwick-3DC Complementary Filter
90 0.004 0.070 0.032
95 0.004 0.083 0.038
99 0.005 0.110 0.050
90 0.019 2.311 0.495
95 0.023 2.754 0.590
99 0.031 3.621 0.776
90 0.147 1.751 0.607
95 0.175 2.086 0.724
99 0.231 2.742 0.951
The confidence level is the probability where the confidence interval formed. It contains the population parameters if a repeated estimation process is carried out. In Table 3, the measurement results in the quantitative analysis of the confidence level using three measurement values of error, namely 90 % (the alpha of 90 % is 0.1), 95 % (the alpha of 95 % is 0.05), and 99 % (the alpha of 99 % is 0.01). Therefore, the equation (22) for measuring the confidence level of each method of observation and three algorithms needs data from each standard deviation in tab. 2. However, from three options, the confidence level used in our research is 95 % because it can provide a balance between precision and reliability.
Tabl. 3 shows the measurement of confidence level with DCM-3DC algorithm has a 0.004 cm in the experiment in case simple ramps. However, the confidence level with the Madgwick-3DC algorithm is 0.083 cm, and the filter complementarity algorithm is 0.038 cm. In the experiment in the cases of ramps and junctions, the result of the confidence level of the DCM-3DC is 0.023 cm, Madgwick-3DC is 2.754 cm, and the complementary filter is 0.590 cm. For the experiment in the cases of double ramps and turns, the result of the confidence level of the DCM-3DC is 0.175 cm, Madgwick-3DC is 2.086 cm, and the complementary filter is 0.724 cm.
Equation (23) is a measure of the confidence interval in our research to provide an estimate of the interval of values used to estimate the true value of the error measurement. Based on each mean value in the DCM-3DC algorithm that has been obtained previously (in tabl. 2), the results of the 0.95 (95 %) confidence interval estimate that the proportion of distance of error for the case of the simple ramp is in the range of 33.92 cm < distance of error < 33.93 cm, for the case of ramps and junctions in the interval is 51.44 cm < distance of error < 51.49 cm, and for the case of double ramps and turns ramps in the interval is 211.99 cm < distance of error < 212.34 cm.
Fig. 13 shows the description of the ground plane visualized in 3D from each case of land contours at three
different locations. Location one is called simple ramps, location two is called ramps and junction, and location three is called double ramps and turns. These images will be compared to the test results of the three different algorithms in fig. 15.
a)
b)
c)
Fig. 13. The ground plane in the case: (a) Simple Ramps, (b) Ramps and Junction, and (c) Double Ramps and Turns
Before showing the results of the three different algorithms, fig. 14. shows the use of the IMU sensor, which shows the orientation of the attitude on the UGS that will experience roll, pitch, yaw rotation. It is also possible to combine rotation simultaneously. Fig. 14. using the DCM-3DC method can show the visualization process from time to time while the UGS is moving. For example, it can be seen from the figure that the process undergoes a
roll condition on the x-axis rotation so that the contour line seems not really straight.
Fig. 15a shows the 3D visualization results obtained from collecting data from 9 DoF from the IMU that has filtered time intervals. The rotary encoder distance range data in the direct observation process 1 (simple ramps) are a graph of the contour estimation of the plot of land using the DCM-3DC, Madgwick-3DC and filters complementarity.
The main feature of orientation estimation from observations on simple ramps, where the ramp starts from +/- 60 cm to +/- 120 cm (see fig. 13a), can be estimated well by the system. Moreover, the results of the contour estimation of the land plot using the proposed methodology (DCM-3DC) have a very satisfactory performance by displaying the results on the estimated plot of land (red line) that almost resembles the actual plot of land (blue line). The qualitative evaluation of performance is in Table 1, where DCM-3DC has a mean error of 1.2861 cm and a MED of 1.2770 cm.
Meanwhile, the Madgwick-3DC methodology (magenta line) displays the results of land plots estimated at +/-120 cm, experiencing actual plane differences. The mean error is 3.4834 cm, and the MED is 3.3018 cm (tab. 1). The Complementary Filter methodology (green line) displays the results of the land plot, which is estimated to be +/-120 cm experiencing a difference in the actual plane so that it has a mean error of 3.0126 cm and a MED of 2.9306 cm (tab. 1).
The data collected consisting of 11 fields: time intervals, distance changes (rotary encoder), 3-axis accel-erometer, 3-axis gyroscope, 3-axis magnetometer will produce a 3D visualization (as seen in fig. 13b) of the measurement data on the observation case of ramps and junctions. Qualitative test with this level of observation difficulty, which is more complicated than the observation in the previous case, is simple ramps.
The estimation of the contour of the plot of land using the DCM-3DC methodology showed satisfactory performance in this observation (ramps and junctions). The main features of this observation are two junctions and one ramp, which can be estimated well by the system. As a result, the estimated plot (red line) nearly resembles the actual plot (blue line).
The summary of the comparison evaluation is shown in tab. 1 on the observation of ramps and junctions. It can be concluded that the DCM-3DC method provides significantly better performance and the results of 3D graphics (see fig. 15b) are always closer to ground truth, even though it has passed through two junctions and one ramp. DCM-3DC has a mean error of 4.9086 cm and MED of 4.3029 cm.
Madgwick-3DC has a mean error of 15.511 cm and MED of 8.7419 cm, so that in fig. 15b, the Madwick-3DC method (magenta line) shows poor performance after passing through two junctions. The estimated ground contour deviates far from ground truth, especially after the second junction. The same result for the Complementary Filter method is almost identical to the Madgwick-3DC method, which has a mean error of 7.5670 cm and a MED of 6.0253 cm. However, the Complementary Filter method (green line) produces a 3D graph of the ground contour, estimated to deviate slightly from the ground truth (see fig. 15b).
As in two previous observations, the data visualization of the measurement results on observations of double ramps and turns (as initially shown in fig. 13c) can represent visualization by collecting data as many as 11 items consisting of: time intervals, distance changes, 3-axis ac-celerometer, 3-axis gyroscope, and 3-axis magnetometer. In addition, qualitative test was carried out to observe double ramps and turns, showing the highest level of difficulty between the two previous observations. Based on tab. 1, the DCM-3DC method has a mean error of 13.3776 cm and MED of 14.0543 cm.
Roll Condition
i Ground Truth Estimation DCM
110 9
I 5 g
>
4 H 3
I'
Fig. 14. Illustration of the rotation on roll condition
a)
b)
c)
Fig. 15. The ground contour and topographic map result with DCM-3DC, Madgwick-3DC, and Complementary Filter in the case (a)
simple ramps, (b) ramps and junction, and (a) double ramps and turns
In estimating the plot contour of land using the DCM-3DC methodology, it also has a satisfactory performance in observing double ramps and turns. In fig. 15c, it can be seen from the DCM-3DC trajectory, which is always closer to ground truth, even though it has passed three turns and two ramps. In the initial conditions, the contour estimation of the plot of land showed perfect results, but over time and distance, there was a deviation between the estimated results of the DCM-3DC (red line) and the actual plot (blue line). This phenomenon is caused by minor errors that accumulate over time. However, the results of the DCM-3DC methodology are able to compensate of the error very well amidst the challenges of a reasonably
extensive distance range (+/- 20 meters) with three turns at an angle of 90o.
Representation of fig. 15c, the Madgwick-3DC method (magenta line) does have a deviation from ground truth, but the deviation is much better than the complementary filter method. On the other hand, the complementary filter method (green line) gives poor performance. The elevation deviation (Z-axis) between the estimate and the ground truth looks very clear and unstable and it is the deviation of the X Y-axis. The two results of these errors can be seen in Table 1, where Madgwick-3DC has a mean error of 25.2415 cm and MED of 17.1909 cm. Moreover, the Complementary Filter has a mean error of 35.8518 cm and a MED of 33.0592 cm.
Conclusions
The new methodology described in this paper is designed for unmanned ground systems obtaining spatial data from ground contours. The data acquisition process utilizes measurement instruments with optical and inertia sensors as digital elevation models to represent the topographic map of one part of the cartographic product. The results of the topographic maps in the three-dimensional model provide a picture of a surface or contour of the land that is easier to understand and easier to adjust to the real world's actual situation than the two-dimensional shape depicted on the map. The result of three-dimensional object model has a dense distance and a tenuous point distance using a three-dimensional representation model in contour lines and contour intervals.
An analysis of three different algorithms at three different locations for ground plane and contour visualization was conducted. Qualitatively, a significant deviation begins to appear in observation double ramps and turns. It can be seen that Madgwick-3DC (magenta line) and Complementary Filter (green line) have a pretty significant deviation relative to the ground truth. In contrast, DCM-3DC (red line) is closer to the ground truth even though it has passed two ramps and two turns. The results show that DCM-3DC provides much higher performance than the other two methods. This performance is seen based on standard measurements (metrics) of ground truth, including; (a) AED, (b) MED, (c) Mean Error Growth, and (d) Median Error Growth. Likewise, for the quantitative analysis to measure the mean, standard deviation, standard error, and confidence level, the DCM-3DC algorithm provides significantly lower error values across metrics than the other methods across all metrics and all observed locations.
The sustainability of the 3D visualization results can be displayed in more detail, effectively, at zero error rate, and complex by utilizing spatial data collection such as sea level height, air pressure, base location systems, and other instruments. In addition, it needs to integrate with a connected system in a global network.
References
[1] Subaweh MB, Wibowo EP. Implementation of Pixel Based Adaptive Segmenter method for tracking and counting
vehicles in visual surveillance. Int Conf on Informatics and Computing (ICIC) 2016: 1-5. DOI: 10.1109/IAC.2016.7905679.
[2] Rostianingsih S, Gunadi K. Pemodelan Peta Topografi ke Objek Tiga Dimensi. Jurnal Informatika University Petra Kristian 2004; 5(1): 14-21. DOI: 10.9744/informatika.5.1 .pp.14-21.
[3] Jenks GF, Steinke T, Buchert B, Armstrong L. Illustrating the concepts of the contour symbol, interval, and spacing via 3-D maps. J Geog 1971; 70(5): 280-288. DOI: 10.1080/00221347108981640.
[4] Yu C, Lee J, Munro-Stasiuk MJ. Extensions to least-cost path algorithms for roadway planning. Int J Geogr Inf Sci 2003; 17(4): 361-376. DOI: 10.1080/1365881031000072645.
[5] Wibowo EP, Talita AS, Iqbal M, Mutiara AB, Lu CK, Meriaudeau F. An improved calibration technique for polarization images. IEEE Access 2019; 7: 28651-28662. DOI: 10.1109/ACCESS.2019.2900538.
[6] Purwanto I, Wibowo EP, Christie DA, Musa P, Harahap RK, Irawan R. Comparative performance evaluation of direction cosine matrix and Madgwick's as 3D orientation estimation algorithm. ICIC Express Letter 2021; 15(4): 409-420.
[7] Doh N, Choset H, Chung W. Accurate relative localization using odometry. Conf on Robotics and Automation 2003; 2: 1606-1612. DOI: 10.1109/robot.2003.1241824.
[8] Jin Y, Toh HS, Soh WS, Wong WC. A robust dead-reckoning pedestrian tracking system with low cost sensors. Int Conf on Pervasive Computing and Communications (PerCom) 2011: 222-230. DOI: 10.1109/PERC0M.2011.5767590.
[9] Fathurrochman D, Musa P, Wimananda DD, Lestari OB. Remote sensing system of odometry and telemetry data in real-time. 3rd Int Conf on Informatics and Computing (ICIC) 2018: 1-6. DOI: 10.1109/IAC.2018.8780484.
[10] Fadlillah HM, Musa P, Wimananda DD. Heading correction in rocket flight system using odometry trajectory information. 3rd International Conference on Informatics and Computing (ICIC) 2018. DOI: 10.1109/IAC.2018.8780428.
[11] 11 Best measuring wheels of 2021 reviewed. Source: (https://cdn.architecturelab.net/wp-content/uploads/2020/03/Best-Measuring-Wheels-3.jpg).
[12] Iafolla L, Filipozzi M, Freund S, Zam A, Rauter G, Cattin PC. Proof of concept of a novel absolute rotary encoder. Sens Actuator A Phys 2020; 312: 112100. DOI: 10.1016/j.sna.2020.112100.
[13] Kljucanin S. INSPIRE specifications in the service of making a topographic map. Geodetski List 2020; 74(4): 373-387.
[14] Rusmanto A. Pedoman Geoparsial Kebijakan Satu Peta [In Indonesian]. Badan Informasi Geospasial 2018. Source: (https://portalksp.ina-
sdi.or.id/Buku_Pedoman_Geoportal%20KSP.pdf).
[15] Nurwadjedi N, Rosalina L, Wibisono Y. Membangun Satu Peta Untuk Penataan Ruang [In Indonesian]. Seminar Nasional Geomatika 2019; 3: 157-166. DOI: 10.24895/sng.2018.3-0.946.
[16] Kluykov AA, KrylovVI. Space geodesy: Past, present and future. To the 50th anniversary of the first set of students MIIGAiK on specialty "Space Geodesy" [In Russian]. Geodezy and Cartography 2019; 80(3): 48-56. DOI: 10.22389/0016-7126-2019-945-3-48-56.
[17] Savinikh VP, Krasnorylov II, Shlapak VV. Training engineers of space geodesy in the MIIGAiK [In Russian]. Geodezy and Cartography 2014; (5): 30-34. DOI: 10.22389/0016-7126-2014-887-5-30-34.
[18] Makarov AP, Chibunichev AG, Poliakova EV. Application of robotic system for obtaining information about the area. ISPRS Archives 2020; XLIII-B2-2020: 611-616. DOI: 10.5194/isprs-archives-XLIII-B2-2020-611-2020.
[19] Liu T, Xu P, Zhang S. A review of recent advances in scanned topographic map processing. Neurocomputing 2019; 328: 75-87. DOI: 10.1016/j.neucom.2018.02.102.
[20] Islami N. The use of online media and topography map in the topic of landslide natural disaster, advanced earth physics course. International Journal of Educational Best Practices 2019; 3(2): 1-9. DOI: 10.31258/ijebp.v3n2.p1-9.
[21] Pawlus P, Reizer R, Wieczorowski M. A review of methods of random surface topography modeling. Tribol Int 2020; 152: 106530. DOI: 10.1016/j.triboint.2020.106530.
[22] Arraigada M, Parti M. Calculation of displacements of measured accelerations, analysis of two accelerometers and application in road engineering. Swiss Transport Research Conference 2006: I-II, 3-30.
[23] Doumiati M, Victorino A, Charara A, Lechner D. Estimation of road profile for vehicle dynamics motion: Experimental validation. American Control Conference 2011: 5237-5242. DOI: 10.1109/acc.2011.5991595.
[24] Cariow A, Cariowa G, Majorkowska-Mech D. An algorithm for quaternion-based 3D rotation. Int J Appl Math Comput Sci 2020; 30(1): 149-160. DOI: 10.34768/amcs-2020-0012.
[25] Corrales-Rodrigáñez C. Rotations and units in quaternion algebras. J Number Theory 2012; 132(5): 888-895. DOI: 10.1016/j.jnt.2011.12.009.
[26] Sarabandi S, Thomas F. A survey on the computation of quaternions from rotation matrices. J Mech Robot 2019; 11(2): JMR-18-1163. DOI: 10.1115/1.4041889.
[27] Chudá H. Universal approach to derivation of quaternion rotation formulas. MATEC Web of Conferences 2019; 292: 01060. DOI: 10.1051/matecconf/201929201060.
[28] Terzakis G, Lourakis M, Ait-Boudaoud D. Modified Rodrigues parameters: An efficient representation of orientation in 3D vision and graphics. J Math Imaging Vis 2018; 60(3): 422-442. DOI: 10.1007/s10851-017-0765-x.
[29] Hamilton WR. Theory of quaternions. Proceedings of the Royal Irish Academy (1836-1869) 1844; 13(48): 1-16. Source: <https://www.jstor.org/stable/20489494).
[30] Griffin S. Quaternions: Theory and applications (Mathematics research developments). Nova Science Publishers Inc; 2017.
[31] Xu FN, Wang B, Deng ZC, Li QJ, Wei Y. Attitude control of targets captured by tethered space robots based on the quaternion theory. Appl Math Mech 2017; 38(12): 13091318. DOI: 10.21656/1000-0887.380168.
[32] Ramasubramani V, Glotzer SC. rowan: A Python package for working with quaternions. J open source softw 2018; 3(27): 787. DOI: 10.21105/joss.00787.
[33] Vince J. Calculus for computer graphics. Springer International Publishing; 2019.
[34] Bekta§ DemIrcI B, Aghayev N. On geometric applications of quaternions. Turkish J Math 2020; 44(4): 1289-1303. DOI: 10.3906/MAT-1907-120.
[35] Baek J, Jeon H, Kim G, Han S. Visualizing quaternion multiplication. IEEE Access 2017; 5: 8948-8955. DOI: 10.1109/ACCESS.2017.2705196.
[36] Novelia A, O'Reilly OM. On geodesics of the rotation group SO(3). Regul Chaotic Dyn 2015; 20(6): 729-738. DOI: 10.1134/S1560354715060088.
[37] Bennett A, Kindratenko V. Quaternion C++ class. Source: <www.ncsa.illinois.edu/People/kindr/emtc/quaternions/).
[38] Lee D, Pernicka H. Vision-based relative state estimation using the unscented kalman filter. Int J Aeronaut Space Sci 2011; 12(1): 24-36. DOI: 10.5139/IJASS.2011.12.1.24.
[39] Zivan Y, Choukrouny D. Dual quaternion Kalman filters for spacecraft relative navigation. AIAA Guidance, Navigation, and Control Conference 2018: 1-29. DOI: 10.2514/6.2018-1347.
[40] Na Y, Bang H, Mok SH. Vision-based relative navigation using dual quaternion for spacecraft proximity operations. Int J Aeronaut Space Sci 2019; 20(4): 1010-1023. DOI: 10.1007/s42405-019-00171-8.
[41] Razgus B, Mooij E, Choukroun D. Relative navigation in asteroid missions using dual quaternion filtering. J Guid
Control Dyn 2017; 40(9): 2151-2166. DOI: 10.2514/1.G002805.
[42] Chelnokov YN. Inertial navigation in space using the regular quaternion equations of astrodynamics. Mechanics of Solids 2019; 54(2): 157-168. DOI: 10.3103/S0025654419030063.
[43] Musa P, Christie DA, Wibowo EP. An implementation of direction cosine matrix in rocket payload dynamics attitude monitoring. Int Conf on Informatics and Computing (ICIC) 2016: 271-276. DOI: 10.1109/IAC.2016.7905728.
[44] Fernández JRM, Junco AH. Design and implementation of an Attitude and Heading Reference System (AHRS) using Direction Cosine Matrix. Revista Cubana de Ciencias Informáticas 2017; 11(1): 15-28.
[45] Tuck K. Implementing auto-zero calibration technique for accelerometers. Freescale Semiconductor Inc; 2007. Source: (https://forum.pololu.com/uploads/default/original/2X/7/7f 6a6201387635b3597e5b19cdb692250e0e3260.pdf).
[46] Li Q, Griffiths JG. Least squares ellipsoid specific fittin. Proc Geometric Modeling and Processing 2004: 335-340. DOI: 10.1109/gmap.2004.1290055.
[47] Shenoi BA. Introduction to digital signal processing and filter design. Hoboken, New Jersey: John Wiley and Sons Inc; 2005.
[48] Premerlani W, Bizard P. Direction cosine matrix imu: Theory. Diy Drone Usa 2009; 13-15. Source: (https://owenson.me/build-your-own-quadcopter-autopilot/DCMDraft2.pdf).
[49] Wang Y. Rajamani R. Direction cosine matrix estimation with an inertial measurement unit. Mech Syst Signal Process 2018; 109: 268-284. DOI: 10.1016/j.ymssp.2018.02.038.
[50] Madgwick SOH. An efficient orientation filter for inertial and inertial/magnetic sensor arrays. 2010. Source: (https://www.samba.org/tridge/UAV/madgwick_internal_r eport.pdf).
[51] Albab AN, Rahmawati E, Yantidewi M, Sucahyo I, Firmansyah RR. Control position of mobile robot based on odometry method and PID controller. J Phys: Conf Ser 2020; 1491(1): 1-7. DOI: 10.1088/1742-6596/1491/1/012039.
[52] Krakiwsky AJ, Wells DE. Coordinate systems in geodesy. Department of Geodesy and Geomatics Engineering, University of New Brunswick 1971. Source: (http://www2.unb. ca/gge/Pubs/LN16.pdf).
[53] Altamimi Z, Gross R. Geodesy. In Book: Springer handbook of global navigation satellite systems. New York: Springer; 2017: 1039-1061. DOI: 10.1007/978-3-319-42928-1_36.
[54] Ormeling F. Cartography. The ideal and its history. Int J Cartogr 2021: 1-4. DOI: 10.1080/23729333.2021.1890343.
[55] Wçgrzyn M, Moscicka A. CityGuideTour Torun - tourist application using augmented reality. Geodesy Cartogr 2017; 66(2): 317-331. DOI: 10.1515/geocart-2017-0018.
[56] Marino L, Cicirello A. Experimental investigation of a single-degree-of-freedom system with Coulomb friction. Nonlinear Dyn 2020; 99: 1781-1799. DOI: 10.1007/s11071-019-05443-2.
[57] Jiang J, Newman ST, Zhong RY. A review of multiple degrees of freedom for additive manufacturing machines. Int J Comput Integr Manuf 2021; 34(2): 195-211. DOI: 10.1080/0951192X.2020.1858510.
[58] Friedman JH, Bentley JL, Finkel RA. An algorithm for finding best matches in logarithmic expected time. ACM Trans Math Softw 1977; 3(3): 209-226. DOI: 10.1145/355744.355745.
Authors' information
Purnawarman Musa (b. 1975), he received his B.CS. Degree in Computer Science from Gunadarma University, Indonesia, in 1999, his M.T. degree in Electronics and Telecommunication from Gunadarma University, Indonesia, in 2006, and his Ph.D. degree in Electronics Informatics from Burgundy University, France, in 2013. He is a member of the Professional Organization in the field of Information and Computer Technology (IPKIN) and Indonesian Association of Higher Education in Informatics and Computing (APTIKOM). He has worked as a Head of the Research Center Aeronautics and Robotics at Gunadarma University. His current interests include the Internet of Things, Robotics, Artificial Intelligent, Sensing, Embedded System, and System on Chips Design. E-mail: p_musa@staff.gunadarma.ac.id .
Imam Purwanto, (b. 1972), he received his B.S. degree in Computer System from Gunadarma University, Indonesia, in 1997; his M.S. degree in Information Management System from Gunadarma University, Indonesia, in 2000; and his doctoral degree in Computer Science from Gunadarma University, Indonesia, in 2020. He worked as an academic service administration staff and as a lecturer in Information System Engineering at Gunadarma University since 1997. His current research interest includes robotic navigation system. E-mail: imampur@staff.gunadarma. ac. id .
Dennis Aprilla Chistie (b. 1993), he received his B.S. degree in Informatics Engineering from Gunadarma University, Indonesia, in 2014, his M.S. degree in Information Management System from Gunadarma University, Indonesia, in 2015, and his M.S. degree in Computer Science from Université de Bourgogne Franche Comté, France, in 2016. He is a lecturer in computer science in Gunadarma University and currently conducting his doctoral degree research in Gunadarma University. His current research interests include 3D computer vision, robotics, machine learning and deep learning. E-mail: dennis_christie@staff.gunadarma.ac.id .
Eri Prasetyo Wibowo, (b. 1966), is a Professor at Gunadarma University, Indonesia. He was born in Kendal, Indonesia in 1966. He received the B.S degree in Electronics and Instrumentation at Gadjah Mada University, Indonesia, in 1991, the M.S. degree in Information System at Gunadarma University, Indonesia in 1994 and the Ph.D. degree in Electronics Informatics from Universite de Bourgogne, France in 2005. He is member of IEEE and EACOVIROE Project to promote Master VIBOT which founded by Erasmus Mundus. His research interests are focused on sensor CMOS design for face tracking and recognition, designed and developed ADC for video application, VLSI design and real time image processing. E-mail: eri@staff.gunadarma.ac.id .
Rudi Irawan (b. 1964). In 1988, he received B.Sc degree in Physics-Electronics & Instrumentation from Gadjah Mada University, Indonesia. His M.Sc. degree is in Physics from the University of Tennessee, Knoxville, USA in 1994 and Ph.D in Physics & Electronics Engineering (Optical Sensing) from the University of New England, Armidale, Australia in 2001. His research interests are in the areas of instrumentation, sensor, photonics, energy efficiency and renewable energy. He worked at Nanyang Technological University, Republic Polytechnic and Nitto Denko Asia Technical Centre. Currently, he works as a lecturer at Mechanical Engineering Department, Gunadarma University, Indonesia. E-mail: rudi_irawan@staff.gunadarma. ac.id .
Code of State Categories Scientific and Technical Information (in Russian - GRNTI)): 29.31.15, 29.33.43, 20.53.23.
Received April 24, 2021. The final version - October 14, 2021.
