Augmented reality for art , design and cultural heritage-system design and evaluation

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.

Introduction
This paper describes the design of an optical see-through head-mounted system for Augmented Reality (AR) and its quantitative and qualitative performance. Augmented Reality is a technique that can be placed in the so-called mixed reality continuum [1], with at one far end the real world that dominates the perception (Reality) and the other end the virtual world that dominates the perception (Virtual Reality); see Figure 1.
In contrast with Virtual Reality (VR), where a complete virtual world must be created, in AR usually only virtual objects or avatars are added to the real world as the rest of the world is the real world. In this paper we focus on mobile immersive AR, which implies that a headset is worn in which the real world view is augmented with virtual objects.
Since in VR only the virtual world is shown, walking with a headset in this world is difficult because the user has little clue in which direction he walks. In Video-See-Through AR the user perceives the real and virtual world by looking at displays in front of his eyes, whereas the merging of both worlds is performed by the digital mixing of video data from the virtual content and the real world. The real world is perceived by two video cameras placed directly before the displays in front of the user's eyes. A problem in this setup is that the real world looks pixilated, that the entire field of view of a person must be covered by the displays, and that the displaying of the real world usually has a delay of one or more hundreds of milliseconds, which might cause motion sickness when walking (for some people), since there is a mismatch between visual information, the information from the inner ear and the information from the muscles [2][3][4].
In Optical-See-Through AR the real world information and the virtual world information is merged through optical mixing using half-translucent prisms. The benefit of this setup is that headsets can be made that are open, as we did in our project. As with normal glasses that people wear, one can also look underneath and left and right of the glasses, relaxing the "scuba-diving" feeling. Since the real world is not delayed at all and one can also look below the displays, walking is in general no problem.
In contrast with Video-See-Through, the real world can only be suppressed by increasing the illumination level of the virtual objects, which is of course limited. Creating dark virtual objects in a bright real world is hence cumbersome.
The biggest problem in AR is to exactly overlay the real and virtual world. This problem has some analogy with color printing, where the various inks must be exactly in overlay to obtain full color prints. However, in AR this is a 3D problem rather than a 2D problem and, worse, the human head can move rapidly with respect to the real world. A first solution was worked out in 1999 [5] after which we refined this in later phases [6,7]. We used one or more visual markers, with known size, position, and distances to each other, which can be found and tracked by a measurement camera on the headset. In order to cope with fast head movements that the camera cannot follow, the head pose data from the camera was merged with data from an inertia tracker. This setup is in analogy with the visual system-inner ear combination of humans. In 2004 HITLab published the AR-Toolkit [8] that used the same type of markers as well as a WebCam in which AR on the computer screen can be displayed. Recently it has been made fit for web-based and iPhone-3GS-based applications.
The ultimate goal of our research, which started in 1998, was to design an immersive, wearable light-weight AR system that is able to provide stereoscopic views of virtual objects exactly in overlay with the real world: a visual walkman, equivalent to the audio walkman. Note, however, that with an audio walkman the virtual music source (e.g., an orchestra) turns with the user when the user turns his head. Using visual anchor points like markers, both virtual visual and virtual audio data can be fixed to a specific location in the real world. Figure 2 shows our current system that evolved during the past decade and that we evaluated during the last three years in real applications.
We measured its accuracy and performance in our laboratory using an industrial robot and in order to get a feeling how the system performs in real life, we tested it with artists and designers in various art, design, and cultural heritage projects in museums and at exhibitions.
The possibilities of immersive AR for applications are plentiful. It can be fruitfully used in area development, architecture, interior design, product design, as it may diminish the number of mock-ups and design changes in too late stage of the process. It can be used for maintenance of complex machines, and possibly in future for medical interventions. A main benefit of AR is that new designs or repair procedures can be shown in an existing environment. Its future possibilities in online gaming and tele-presence are exiting. Our initial application idea was to provide a tool for guided tours and a narrative interface for museums. Hence, with the AR system, one must be able to easily roam through indoor environments with a head-tracking system that is largely independent of the environment.
Similar AR systems exist already, such as LifePLUS [9] and Tinmith [10] but they use video-see-through methods which makes registration easier but at the cost of loss of detail of the world. Other projects like BARS [11] and MARS [12] use optical-see-through methods but do not care for precise pose tracking or do not use a camera for tracking.
In the remainder of this paper we describe the technical setup of our system (Section 2) and its application in art, design, and cultural heritage projects (Section 3). Figure 3 shows the components of the system.  Note that Figure 2 shows a prototype AR headset that, in our project, was designed by Niels Mulder, student of the Postgraduate Course Industrial Design of the Royal Academy of Art with as basis the Visette 45SXGA.

Main System Setup.
Off-line virtual content is made using Cinema-4D [19]; its Open-GL output is online rendered on the laptop to generate the left and right-eye images for the stereo headset. The current user's viewpoint for the rendering is taken from a pose prediction algorithm, also online running on the laptop, which is based on the fusion of data from the inertia tracker and the camera, looking at one or more markers in the image. In case more markers are used, their absolute positions in the world are known. Note that also markers with no fixed relation to the real world can be used. They can be used to represent moveable virtual objects such as furniture.
For interaction with virtual objects a 5DT data glove [20] is used. A data-glove with RFID reader (not shown here) was made to make it possible to change/manipulate virtual objects when a tagged real object is touched.

Head Pose Tracking. The Xsens MTx inertia tracker
[15] contains three solid state accelerometers to measure acceleration in three orthogonal directions, three solid state gyroscopes to measure the angular velocity in three orthogonal directions, and three magnetic field sensors (magnetometers) that sense the earth's magnetic field in three orthogonal directions. The combination of magnetometers and accelerometers can be used to determine the absolute 3D orientation with respect to the earth. The inertia tracker makes it possible to follow changes in position and orientation with an update rate of 100 Hz. However, due to inaccuracies in the sensors, as we integrate the angular velocities to obtain angle changes and double integrate accelerations to obtain position changes, they can only track reliably for a short period. The error will grow above 10 to 100 meter within a minute. This largest error is due to errors in the orientation that leads to an incorrect correction for the earth's gravitational pull. This should be corrected by the partial, absolute measurements of the magnetometers, as over short distances the earth's magnetic field is continuous; but this field is very weak and can be distorted by metallic objects nearby. Therefore, although the magnetic field can be used to help "anchoring" the orientation to the real world, the systematic error can be large depending on the environment. We measured deviations of 50 • near office tables. Hence, in addition to the magnetometers, other positioning systems with lower drift are necessary to correct the accumulating errors of the inertia tracker. A useful technique for this is to use visual information acquired by video cameras. Visual markers are cheap to construct and easily mounted (and relocated) on walls, doors, and other objects. A marker has a set of easy detectable features such as corners or edges that enable recognition of the marker and provide positional information. Many different marker types exist, circular [21] or barcode like [22]. We chose a marker with a rectangular border to be able to easily detect and localize the marker and chose a 2D barcode as its identity is detectable even when the marker is very small (Figure 4).
If the marker is unique, then the detection of the marker itself restricts the possible camera positions already. From four coplanar points, the full 6D pose can be calculated with respect to the marker with an accuracy that depends on the distance to the marker and on the distance between the points. In case more markers are seen at the same time, and their geometric relation is known, our pose estimation will use all available detected points in a more precise estimation. In a demo situation with multiple markers, the marker positions are usually measured by hand.
Tracking is not restricted to markers, also pictures, doorposts, lamps, or all that is visible could be used. However, finding and tracking natural features, for example, using SIFT [23,24], GLOH [25], or SURF [26] comes at a cost of high process times (up to seconds as we use images of 1280 × 1024), which is undesirable in AR due to the possibility of a human to turn his head very quickly. To give an impression: in case of a visual event in the peripheral area of the human retina, after a reaction time of about 130 ms in which the eye makes a saccade to that periphery, the head starts to rotate accelerating with 3000 • /s 2 to a rotational speed of 150 • /s to get the object of interest in the fovea. When the eye is tracking a slow moving object (smooth pursuit) the head rotates with about 30 • /s [27,28].
Moreover, sets of natural features have to be found that later enable recognition from various positions and under various lighting conditions to provide position information. The biggest issue with natural features is that their 3D position is not known in advance and should be estimated using, for instance, known markers or odometry (Simultaneous Localization And Mapping [29,30]). Hence, we think that accurate marker localization will remain crucial for a while in mobile immersive AR.

Required Pose Accuracy.
The question rises what should be the accuracy of a tracking system if we want to have adequate alignment of virtual and real objects. For an eye with a visual acuity of about 0.01 • , looking through a headmounted display at 10 cm distance with an opening angle of 36 • × 27 • , we actually need a resolution of about 3000 × 2000 pixels. As our HMD has 1280 × 1024 pixels the maximum accuracy we can obtain is one pixel of our display, which translates to roughly 0.03 • or 0.5 mm at 1 meter distance of the eye. Hence, currently an AR user at rest will always perceive static misalignment due to the limitations of the HMD. Dynamically, we can present virtual objects on our HMD at a rate of 60 Hz. Assuming instantaneous head pose information from the pose measuring system, and assuming head movements in smooth pursuit we obtain a misalignment lag of 1/60 * 30 • /s = 0.5 • . If we assume head motions as reaction on attention drawing, we obtain a temporary misalignment lag due to head movements of 1/60 * 150 • /s = 2.5 • . Consequently, with the current headset technology the user will inevitably notice both static and dynamic misalignment due to head motion. Reasoning the other way around, the extra dynamic misalignment due to the current headset cannot be noticed (less than the static misalignment) if we rotate our head with less than 0.03 * 60 = 1.8 • /s. Concluding, the target accuracies for our pose measurement system are based on the accuracies for the pose of virtual objects that can be realized by the current HMD and we can distinguish three scenarios.
(i) A static misalignment of <0.03 • , that is, a position misalignment of <0.05 cm of a virtual object at 1 m.
(ii) A dynamic misalignment of <0.5 • when smoothly pursuing an object, that is, a temporal position error of < 0.9 cm of a virtual object at 1 m.
(iii) A dynamic misalignment of <2.5 • when another event in the image draws the attention and the head rotates quickly, that is, a position error of <4.3 cm of virtual object at 1 m.
These are theoretical values. Given the flexible and versatile human vision system users might not find these errors disturbing. We address this in Section 3.

Camera-Only
Tracking. Below we describe our methods to calculate the pose of a camera from an image of a known marker. Our aim was to use as few markers as possible, ultimately a single marker seen from quite a distance. Hence, we also use a lens with a very large opening angle of 108 • .
We investigated the influence of image noise and parameters such as line thickness and marker size on the accuracy of the estimated pose. We used a rectangular pattern with a big black border on a white field with inside a 2D barcode to identify the individual markers [7,8] (see Figure 4). Figure 5 shows the real-time image processing steps that we use to track the pose of the camera with respect to a marker.
To minimize latency we need fast methods. Therefore, we first detect candidate markers (single closed contours) using a Canny edge detector, with a fixed threshold on the gradient to suppress noise from the imaging system. While following the edges in the Canny algorithm we keep track of connected edge points and count the number of points that are not part of a line (end-points, T crossings, etc.). Only contours with no special points (single closed contour) are interesting.
Then we search for corners only along these contours and keep contours with four corners. The corners are found by using a modified Haralick-Shapiro corner detector [31,32]. As the gradients are high on the edge, we only need a threshold on their circularity measure and search for local maxima of that measure along the edge. After splitting the contour in the four segments, we find the accurate location of the edge points, correct for lens distortions, and fit a line through each segment. The intersections of the lines give an unbiased location of the four corners needed for pose estimation. Other corner detectors as [31][32][33] did not perform well as they need either a large patch around the corner (impairs speed and makes them less robust against nearby other edges) or have a bias in their estimate. To reach our unbiased estimate we had to correct the location of the edge points for lens distortion prior to fitting the lines.
Accurate edge-point locations are crucial to find accurate corner points; hence, we must eliminate systematic errors and noise as well as possible [34,35]. Using the step-edge model (Gaussian blurred edge) we can calculate the edge location accurately from three pixels centered on and perpendicular to the edge. To increase processing speed we evaluate three pixels along the horizontal or vertical direction, depending on which one is most perpendicular to the edge. Where usually the gradient magnitudes are used to find the location as the top of a parabola, we use the logarithm of the gradients. This makes sure that the parabolic profile assumption is valid for sharp images as well, and an unbiased estimate for the edge location of our model edge is obtained. In an experiment with a linearly moving edge the bias in location was measured to be up to 0.03 px without the logarithm, and 0.01 px with the logarithm.
We first investigated the influence of the thickness of the black border on our step-edge locator. We found that when the black border is thicker than 8 pixels in the image, the edge points on the outer contour of the border can be located with practically zero bias and an RMS error <0.01 pixel using integer Gaussian derivative operators with a scale of 1.0 px. We use integer approximations of the Gaussians because of Grab   their fast implementations using SIMD instructions. Using simpler derivatives, this bias will stay low even at a thickness of 3-5 pixels; however, this error is then symmetrically dependent on the subpixel location of the edge. If a large number of points are used for fitting a line through the edgepoints-usually 12-30 points are used-the bias error can be regarded as a zero mean noise source, but for short edges the fit will have an offset. We tried several edge detectors/locators and in the presence of noise, the most accurate and robust detector was using an integer Gaussian derivative filter with the three gradient magnitude values to calculate the edge position not from neighboring pixels but from pixels at a distance of two pixels, provided that the line thickness was big enough. We used this detector but with three neighboring pixels as we expect line thicknesses of near five pixels (markers at a few meters distance). The detector to use in other situations should be chosen on basis of the expected line thickness and noise, for example, marker distance, marker viewing angle, and illumination (indoor/outdoor) circumstances.
We then determined the size of the marker pattern that is needed when it should be detected at 5 m distance under an angle of 45 • . With a 5-pixel line thickness and leaving 2 × 2 pixels for the black and white blocks, the minimum size of a marker is 18.2 × 25 cm, fitting on A4. The bias per edge location will then be between 0.01 and 0.04 pixels, depending on the scale of the edge. When the camera is not moving, the scale is 0.8 pixels corresponding to a bias of 0.01 pixels. Because the edge location has only a small bias, the error of our algorithm is noise limited, and in the absence of noise, it is model limited.
We then verified our step-edge model and found that it fits well to experimental data. We still found a bias of around 0.004 pixel and an RMS error around 0.004 pixel as well. This bias we attribute to the small error we still make in assuming a Gaussian point spread function of the imaging system. When the Contrast to Noise Ratio-CNR = 2a/σ noise -is around 26 dB, the standard deviation of the edge location is 0.1 pixel. This is also the residual error of the saddle points after a lens calibration.
When the CNR is higher, the biggest source of error in our experimental setup seems to be the (model of the) lens. In order to be able to use a pinhole camera model, we tried to calibrate all distortions away, but even with an elaborate lens distortion model we obtained a residual calibration error of 0.37 pixel maximum (standard deviation 0.1 pixel). We found an increased blurring at the borders of the image, suggesting lens artifacts. In photography, these artifacts are minimized using more elaborate lens systems. More research is needed to investigate how to further reduce this systematic error, with a better lens (model) as a starting point. Our lens distortion model is given by with r u = − → p U ; D and U denote distorted/undistorted metric sensor plane coordinates. This model performs better in our case than the other models we tried [36][37][38][39]. The parameters were estimated using the Zhang calibration method [38]. We found that we can detect the contours of a marker robustly down to a CNR of 20 dB and now we only need to worry about the detection of the four corners along these contours. The Haralick-Shapiro corner detector [31,32] is the least sensitive to noise while it performs well along the Canny edge, and we found it can be used with CNR ratios higher than 20 dB. Along the edge we can reliably detect corners with an angle of less than 120 • . When the CNR is 25 dB, corners can be detected up to 150 • . Corner angles of 120 • and 150 • relate to marker pitch angles of 35 • and 65 • , respectively. To realize our target of detecting the marker up to pitch angles of 60 • , we need the CNR to be around 25 dB.
For online estimation of the pose from four corners we used a variation of the Zhang calibration algorithm; only the external parameters need to be estimated. Using static measurements to determine the accuracy of our pose estimation algorithm we determined that the position of a marker in camera coordinates is very accurate when the marker is on the optical axis at 5 m, that is, less than 0.5 mm in x and y, and less than 1 cm along the optical axis. The marker orientation accuracy, however, highly depends on that orientation. The angular error is less than 5.2 • (1.5 • due to noise) when the marker pitch is less than 20 • at 5 m. When we convert the marker pose in camera coordinates to the camera pose in marker coordinates, the stochastic orientation error results in an error in position of 2.7 cm/m. With a pitch larger than 20 • , the orientation accuracy is much better, that is, less than 1.4 • (0.5 • due to noise), resulting in a stochastic positional error of the camera of less than 0.9 cm/m. Hence, markers can best be viewed not frontally but under a camera pitch of at least 20 • . Finally, with this data, we can determine the range where virtual objects should be projected around a marker to achieve the required precision for our AR system. We found that with one marker of size 13 × 16.5 cm (at 1.5 m-6 m from the camera), a virtual object should not be projected at more than 60 cm from that marker in the depth direction, or within 1 m from that marker in the lateral direction to achieve the target accuracy of 0.5 • error in the perceived virtual object position.

Camera Data Fused with Inertia Data.
We need fast inertia data to keep up with fast head movements. However, cheap solid-state inertia trackers build up severe pose errors within a second. Consequently, these pose measurements should be corrected using the much slower but more stable camera pose data that is acquired by locking onto features of markers in the real world. We used an inertia tracker fixed onto a camera. Our sensor fusing Kalman filter [40,41] combines the absolute pose estimate from the camera with acceleration sensors, angular velocity sensors and magnetic sensors to get a better estimate of the HMD pose. The Kalman filter is also necessary to interpolate the pose inbetween the slow pose estimates from the camera. Figure 6 shows the problem we encounter when we fuse pose data from the camera with pose data from the inertia tracker. The inertia pose data has a frequency of 100 Hz. The camera with image processing has an update rate of about 15 Hz. Note that the online viewpoint-based rendering costs also time. The Kalman filter with inertia tracker data can be used to predict the head pose at the precise moment we display the virtual objects precisely aligned on the headset.
From now on, we refer to the pose of the camera with respect to a marker at a certain point in time as its state. This state does not only include the position and orientation of the camera at that point in time, but also its velocity and angular velocity, and where necessary their derivatives. The error state is the estimation of the error that we make with respect to the true state of the camera.
Our fusion method takes latencies explicitly into account to obtain the most accurate estimate; other work assumes synchronized sensors [42,43] or incorporates measurements only when they arrive [44] ignoring the ordering according to the time of measurement.
Our filter is event based, which means that we incorporate measurements when they arrive, but measurements might be incorporated multiple times as explained next.
We synchronize the camera data with the filter by rolling back the state updates to the point in time at which the camera has acquired its image. We then perform the state update using the camera pose data and use stored subsequent inertia data again to obtain a better estimate of the head pose for the current point in time, and to predict a point of time in the near future, as we need to predict the pose of the moving head at the moment in time that the image of the virtual objects are projected onto the LCD displays of the headset. In this way, we not only get a better estimate for the current time, but also for all estimates after the time of measurement; this was crucial in our case as camera pose calculations could have a delay of up to 80 ms, which translates to 8 inertia measurements.
A Kalman filter can only contribute to a limited extend to the total accuracy of the pose estimates. The estimate can only be made more accurate when the filter model is accurate enough; that is, that the acceleration/angular speed is predictable, and that the inertia sensors are accurate enough. A bias in the sensors-for instance caused by a systematic estimation error or an unknown delay in the time of measurement-will prevent the filter from giving a more accurate result than the camera alone. We minimized the errors introduced by the Kalman filter by using robust methods to represent the orientation and time update of the orientation, and decreased the nonlinearity be using a nonadditive error state Kalman filter in which the error state is combined with the real state using a nonlinear function (see the transfer of the orientation error in Figure 8). We used Quaternions [45] for a stable differentiable representation.
To make the orientation model more linear, we used an indirect Kalman filter setup where the error states are estimated instead of the actual state. Due to this choice the error-state update is independent of the real state. Effectively we created an extended kalman Filter for the error state. If the error state is kept at zero rotation by transferring the error-state estimate to the real state estimate immediately after each measurement update, the linearization process for the Extended Kalman Filter [46] becomes very simple and accurate. In addition, we convert all orientation measurements to error-quaternions: q e,k = q −1 k|k−1 ⊗ q m,k . This makes the measurement model linear (the state is also an errorquaternion) and stable in case of large errors, at the expense of a nonlinear calculation of the measurement and its noise.
In simulations we found that the position sensor accuracy has the largest influence on the total filter accuracy in absence of orientation errors. Changing the sampling rates or using more accurate acceleration measurements had less influence. We can argue that when the process noise in acceleration (or angular velocity for that matter) due to the user's motion is high compared to the measurement noise of the inertia sensors, it is of little use to filter the inertia sensor measurements, meaning that a computationally cheaper model can be used in which the inertia sensors are treated as an input during the time update. Figure 7 shows the process models of the two Kalman filters as we implemented them. The orientation-error Kalman filter at the top estimates errors in orientation and errors in gyroscope bias. The position-error filter estimates errors in position, speed, and accelerometer bias. When gyroscope and accelerometer data is received-they are transmitted simultaneously by the inertia tracker-all real states are x = (dp, dv, db) T Application σ da

Process model
Accelerometers σ q0 , σ dω Gyroscopes Figure 7: The prediction steps of the two implemented error-state Kalman filters and separately maintained position and orientation states when gyroscope and accelerometer data is processed.
updated. In addition, both filters perform a prediction step using their respective process models. In our current setup, we immediately transfer predicted errors to the real states, so the error states will always be zero-or more precisely, they indicate zero error. With zero error input, the output of the prediction step will also be zero. However, the uncertainty of this zero error will increase due the noisy measurements and the expected change in the acceleration and angular velocity. These expected changes should be provided by the application. In our demos we did not make special assumptions for the motions and used the same process noise values for all axes. For the position-error filter we could find a full solution for the process noise due to acceleration change and bias change. We could also find a full solution for the orientation-error filter's process noise. The resulting equation, however, was not practical for implementation. We further assumed the angular velocity to be zero and used the result presented in the figure. The process noise values can be increased a bit to account for the error in this extra assumption, but in practice these values are determined experimentally already. Position observation the best estimate, we reorder all measurements by their measurement time. Therefore, when a camera measurement is received, both error-state filters and the states themselves are rolled back synchronously to the closest state n to the time t, the capture time of the image for the camera pose measurement. All measurements taken after time t n will now be processed again, ordered in time. This reprocessing starts at state i = n. Gyroscope and accelerometer measurements are again processed using the process models, and they will advance the state i → i + 1. Position and orientation measurements will be used to update the a priori estimates at state i to a posteriori estimates in the observation update steps of the Kalman filters. First, these measurements need to be transformed into error observations. We do this using the nonlinear transformations, and thereby circumvent the linearization step of the measurement model for better accuracy. Then, these error measurements are incorporated using the standard Kalman observation update equations. The resulting estimates of the errors are transferred to the separately maintained states of position, orientation, bias and so forth. Hence, all pose estimates up to the present time will benefit from this update.

AR System
Accuracies. Finally, we measured our complete tracking system: camera, inertia tracker and Kalman filter, using an industrial robot as controllable motion source and a marker at 3.2 m. The robot motions are shown in Figure 9. The positional accuracy of the system is shown in Figure 10. The values along the x-axis were the most inaccurate. Without the filter to correct for the slow and delayed camera measurements, the positional error would be up to 20 cm depending on the speed of the robot ( Figure 10(a)). With the filter, the accuracy is generally just as good as the accuracy of the camera measurements. The camera pose shows a position dependent systematic error of up to 3 cm (Figure 10(b)). This proved to be due to a systematic error in the calculated orientation from the camera. When we correct for the orientation error, the positional error becomes less than 1 cm (Figure 10(c)). However, in normal situations the ground truth orientation will not be available. Using the orientation from the inertia tracker did not help in our experiments; the high accelerations are misinterpreted as orientation offsets, which introduces a systematic error in its output.
From our experiments we conclude that our data fusion does its task of interpolating the position in between camera measurements very well.
The tracking system has an update rate of 100 Hz. However, the pose estimates-albeit at 100 Hz-were less accurate than the estimates from the camera because of the high process noise (unknown jerk and angular acceleration from user movements).
We measured that the required orientation accuracy of <0.5 • when moving slowly can be met only when the encountered systematic error in camera pose estimation is ignored: 1 cm at 3 m translates to 0.2 • . Since the camera is the only absolute position sensor, the encountered error of up to 4 cm (0.9 • ) cannot be corrected by inertia tracker data.
Ways to diminish this static error are the following. (ii) Use multiple markers, spread out over the image; this will average the pose errors. (iii) Find ways to calibrate the lens better, especially at the corners. (iv) Use a better lens with less distortion.
A systematic static angular error leads to the fact that an acceleration measured by the inertia tracker is wrongly corrected. This is also visible in static situations due to the acceleration due to gravity. For example with a 1 • error, the Kalman filter will first output an acceleration of sin(1 • ) * 9.81 = 17 cm/s 2 , which is slowly adjusted by the filter since the camera indicates that there is no acceleration. When the camera faces the marker again with a zero error, the wrongly estimated accelerometer bias now generates the same error but then in the other direction and hence this forms jitter on the pose of the virtual object. We found that the bias of the accelerometer itself is very stable. When the process noise for this bias is set very small, the bias will not suffer much from this systematic error. To counter a systematic orientation error it seems more appropriate to estimate a bias in the orientation. However, when the user rotates, other markers will come into view at another location in the image, with another bias. The real effective solution is to minimize camera orientation errors. However, knowing that systematic errors occur we can adapt our demos such that these errors are not disturbing, by letting virtual objects fly for instance. Of all errors, jitter is the most worrying. This jitter is due to noise in the camera image in bad illumination conditions and due to the wrong correction of the earth gravitational field. Note that the first jitter also occurs in, for example, ARToolkit. Jitter in virtual objects makes that it draws the attention of the user, as the human eye cannot suppress saccades to moving objects. Finally, to make a working optical-see-through AR system, many extra calibrations are needed, such as the poses of the sensors, displays, and the user's eyes, all of them crucial for accurate results. Most of these calibrations were done by hand, verifying a correct overlay of the virtual world with the real world.

Application in Art, Design, and Cultural Heritage
In order to obtain insight in how the AR system performs also in qualitative sense, we tested it with artists and designers in various art, design, and cultural heritage projects. The application of artists and designers and curators is of course in no way a replacement for a full user study, but it did lead to some useful observations for the profitable use of the system. For this, within the context of the projects Visualization techniques for Art and Design (2006)(2007) and Interactive Visualization techniques for Art and Design (2007)(2008)(2009)) the Royal Academy of Art (KABK), the Delft University of Technology (TUD), and various SME founded an AR lab [47] in which two prototype AR systems had been developed and tested. The aim of the first project was to research the applicability of AR technique in art and design and to disseminate the technology to the creative industry. The aim of the second project was to combine AR with interaction tools and disseminate the technology to public institutes like museums. The basic idea behind this cooperative projects was that AR technology is new; hence designing with it has no precedent and most probably needs a new approach. Like the first iron bridge (1781); being the first of its kind and therefore its design was based on carpentry, for example, using dovetails [48]. A number of projects have been realized within the context of the ARlab, some of which are recalled below.
30/1/2007 Augmented Letter Soup. The 325th anniversary of the typography design institute of the KABK leads to a project where AR was used to combine typography with interior and interaction design. Wearing the AR headset, users can experience a virtual, typographic interior placed in a real, physical environment and write text in augmented space using 3D, animated letters attached to tangible optical markers; see Figure 11.
By juxtaposing the markers, representing letters of the alphabet, the visitors could write their own name or a short (ii) motion of the virtual objects covers up for jitter; the human attention is already drawn and the jitter is less noticed. The same is true if the human moves; (iii) virtual objects are not bound to the floor, ceiling, walls, or tables; they only need to be within some distance to their nearest marker(s). This means that also information display and interaction does not necessarily have to take place on a wall or table, but might also take place in the air; (iv) the image of the tracker camera can also be used to beam the augmented view of the user on a screen, by which a broad audience can see (almost) through the user's eye.
10-15/4/2007. Augmented Reality Theater. It was an interactive installation at the unDEAF/DEAF festival with virtual 3D animated puppets in AR, data gloves, and physical objects tagged with RFID. Using a data-glove the user could control the position and face expression of an animated virtual puppet. In various physical objects an RFID tag was hidden, which was used to trigger changes in the behavior and looks of the puppet, Figure 12. The following are the Lessons learned: (i) using design packages such as Cinema 4D enlarges the possibilities of the interaction designers; making interaction with animated figures possible; (ii) for real 3D animated films with large plots, game engines must be used; (iii) manipulation of real objects that influence (through RFID) the virtual world is "magic" for many people; (iv) more image processing on the tracker camera is useful, for example, to segment the user's hand and fingers making unhandy data gloves superfluous.   (i) the sounds that the ellipsoids made were coupled to their 3D position, which added to their pose recognition by the user and made it possible to draw his attention; (ii) by applying VR design techniques (i.e., normally in AR only objects are drawn; the walls and floors are taken from the real world) the virtual objects seem real and the real objects, that is, humans walking around, appear virtual or ghosts; (iii) the graphics rendering done on the laptop to generate the stereoscopic view does not show entirely geometric correct rendered images. Research is needed into rendering for AR headsets, taking the deformation of the presented images by the prisms into account; (iv) using image processing on the tracker, the camera can be used to segment walking persons, thus enabling virtual objects (e.g., birds) to encircle them realistically.
16-21/4/2008: At the Dutcheese exhibition at the Salone Internazionale del Mobile 2008 in Milan, apart from real furniture and textiles designs, a large space was augmented with animated furniture and interactive textile (RFID tagged). Two AR systems were simultaneously used, making it possible for the bearers to discuss the designs each from a     museum that was actually out of sight. See Figure 18. The following are the Lessons learned: (i) using a standard laptop is on the one hand rather heavy to wear but does enable fast connection of new interaction devices such as the Wii, but also webcams; (ii) webcams can be used to generate life video streams inside the virtual world.    objects from the period 1450-1550. Sgraffito is an ancient decorative technique in which patterns are scratched into the wet clay. This archaeological collection was made accessible for a broad public using 3D visualization and reconstruction techniques. The original objects were scanned in a CT system and after processing the data, the handmade plates, bowls and pots and their relief decorations were both rendered in virtual representations and rapid prototyped to provide 3D copies of the originals. In the exhibition, video projections show the actual CT scans; whereas the virtual renderings enable visitors to view the objects from all angles. The printed clones competed with the hundred original Sgraffito objects in the exhibition. In screen-based AR setups, the visitors could manipulate objects by manipulating special rapid prototyped pieces of the earthenware with markers or by browsing through a book of markers; at each page one object was visualized and explained in text and accompanied by music from the era of the pottery. Headset-based AR was used in a setup in which virtual plates and bowls could be arranged on a large   The following are the Lessons learned: (i) augmented reality can be fruitfully used to attract a broad public to displays of cultural heritage. Its narrative power is huge; (ii) screen-based AR is a low cost replacement of HMD based AR and can be fruitfully used to introduce the topic at hand and the AR technology itself; (iii) HMD-based AR is at its best when a full immersive experience is required and people can walk around larger objects. (i) for outdoor AR it is necessary that the ambient light intensity and the intensity of the LCD displays on the HMD are in balance. Hence also the real world light intensity needs to be controlled, for example, using self-coloring sunglass technology.

Conclusions
In this paper we described the design of an optical-seethrough head-mounted system for indoor and outdoor roaming Augmented Reality (AR) and its quantitative and qualitative evaluation. Our ultimate goal was that virtual world objects are indistinguishable from real world objects. Hence, for optical see-through AR, measuring the head movements with respect to the physical world is mandatory.  As it makes no sense to have the alignment better than the resolution of the current headset displays, this forms the theoretical limiting factor for the head-pose tracking system: a static misalignment of <0.03 • , a dynamic misalignment, when smoothly pursuing an object of <0.5 • and a dynamic misalignment of <2.5 • when an event in the image draws the attention. Based on these requirements we developed a head-mounted AR system, of which the hardest problem was to develop an accurate tracking system. We implemented a combination of camera and inertia tracking, alike the human visual/vestibular system. Although our ambition was to use natural features, we had to focus on a marker tracking camera system, as for now the processing of natural features is still too slow for this application. After realizing two prototypes, one of which incorporated a redesign of the head-mounted displays, making it more lightweight and open, we measured our system by mounting it on an industrial robot to verify if our requirements were met. To obtain qualitative conclusions, an ARlab was founded with the Royal Academy of Art (KABK), the Delft University of Technology (TUD), and various SME as partners, and we tested the system with artists, designers, and curators in art, design, and cultural heritage projects. This collaboration provided us with very useful observations for profitable use of the system. 14 EURASIP Journal on Image and Video Processing 4.1. Quantitative Conclusions. We can conclude that our tracker based on the fusion of data from the camera and the inertia tracker works well at 100 Hz, albeit that the required orientation accuracy of 0.5 • when moving the head slowly (smooth pursuit) is just met with one 13 × 16.5 cm marker at 5 m distance when the camera's systematic orientation error can be calibrated away. Because the camera is the only absolute position sensor to "anchor" to the real world, these errors cannot be corrected by the inertia sensors. In addition, to obtain this error one has to view the markers under an angle of more than 20 • , which restricts the user's movements a bit. However, the real improvement should come from a more accurate lens calibration or better lens, and/or higher resolution cameras and/or putting more markers, with known geometric relations, in the field of view of the camera and/or using natural features in combination with markers. The current systematic error, that is dependent on the location of the marker in the image, is compensated by the Kalman filter using the bias states, leading to over and undershoots upon user movements. This leads to visible jitter of the virtual objects on top of jitter from noisy camera measurements when the marker is far away or the illumination conditions are not within range.
Although, the jitter is visible for the user, it is not as bad as it seems as the human eye seems to cope with it; the fovea tracks the virtual objects especially when they move.

Qualitative Conclusions.
Lessons learned from experiments with audience on various events and exhibitions showed the following.
(i) The augmented view can be peeked from the tracker camera and used to let the public see through the user's eye.
(ii) Information display and interaction do not necessarily have to take place on a wall or table, but might also take place in the air.
(iii) Positioning virtual objects in the air covers up for static misalignment.
(iv) Motion of the virtual objects covers up for misalignment and jitter; the human visual attention is already drawn by the motion of the object. The same is true when the user moves.
(v) Design packages such as Cinema 4D make design with animated figures possible. For real 3D animated films with large plots, game engines must be incorporated.
(vi) Manipulation of real objects can influence (through RFID) the virtual world. This is "magic" for many people.
(vii) More image processing on the tracker camera is useful, for example, to segment the user's hand and fingers making unhandy data gloves superfluous. Segmenting walking people enables virtual objects to encircle them.
(viii) The sound that virtual objects make adds to their pose recognition and attention drawing.
(ix) By applying VR design techniques, virtual objects appear real and real objects virtual.
(x) More research is needed into the rendering of stereoscopic images for AR headsets, taking the deformation of the presented images by the prisms into account.
(xi) Design discussions are more vividly using HMD based AR as each user can now individually select his (the best) viewpoint.
(xii) Standard laptops are heavy to wear but enable easy connections to new interaction devices such as the Wii.
(xiii) Life video streams inside the virtual world give a telepresence awareness.
(xiv) Screen-based AR is a low cost replacement of HMD based AR and can be fruitfully used to introduce the topic at hand and the AR technology itself.
(xv) Headset-based AR is at its best when a full immersive experience is required and people can walk around larger objects.
(xvi) For outdoor AR it is necessary that the ambient light intensity and the intensity of the LCD displays on the HMD are in balance.
(xvii) Augmented reality can be fruitfully used to attract a broad public to displays of cultural heritage as a three-month exhibition in museum Boijmans van Beuningen in Rotterdam showed. Its narrative power is huge.
The collaboration between researchers in the area of image processing with artists, designers, and curators appeared to be very fruitful and has led to many amazing productions and exhibitions.


extfieldsofstudy:
Computer Science
Medicine

text:
On calculating the probability of a set of orthologous sequences.

Probabilistic DNA sequence models have been intensively applied to genome research. Within the evolutionary biology framework, this article investigates the feasibility for rigorously estimating the probability of a set of orthologous DNA sequences which evolve from a common progenitor. We propose Monte Carlo integration algorithms to sample the unknown ancestral and/or root sequences a posteriori conditional on a reference sequence and apply pairwise Needleman–Wunsch alignment between the sampled and nonreference species sequences to estimate the probability. We test our algorithms on both simulated and real sequences and compare calculated probabilities from Monte Carlo integration to those induced by single multiple alignment.

Introduction
Comparative genomics/proteomics research often focuses on a set of orthologous sequences arising from evolutionary speciation. For example, multiple related species (for example, human, mouse, and rat) can have a common gene as well as the corresponding promoters in the upstream region of such a gene, although these matched sequences may have minor difference across species. For simplicity the set of sequences studied in the sequel are assumed to have almost equal length in light of these examples. Sequence alignment algorithms 1 have substantially facilitated comparative genomics/proteomics research by showing conservation pattern along orthologous sequences, and biologically functional segments are likely to be those more conserved regions along the genome. For the vast body of related literature, we refer to Liu and colleagues, 2 Kellis and colleagues, 3 Moses and colleagues, 4 Xie and colleagues, 5 Wei and Jensen, 6 Sinha and He, 7 and many others. As another major tool, statistical modeling approaches are devoted to comprehensively describing the probabilistic uncertainties linked to those established biological evolution models which may include two topological structures: parallel and phylogenic models (see Figure 1).
The joint parallel evolution process probability Pr(Ancestor, Species 1, 2, and 3) is Pr(Ancestor) ∏ Jukes and Cantor 8 proposed the first probabilistic nucleotide evolution model which assumes substitution to take place randomly among four types of nucleotides "A[1]T [2]C [3]G [4]". The transition (from nucleotide i to nucleotide j) probability up to time t is derived as We assume that the substitution rate parameter (α) is constant for different species and the evolution duration (t) is represented by specific time period (t 0 , t 1 , t 2 , or t 3 ) for the associated divergence process (see Figure 1). Our question is how to effectively estimate the marginal probability for the given orthologous species sequence set without knowing the genotype of the ancestor and/or root.

Material and methods
For the unknown ancestor sequence, we simply assume that the nucleotide on any site follows a tetranomial distribution with categories {ATCG} and equal proportion, (1/4). We further assume that each nucleotide on the ancestor sequence evolves independently (under the probability law, Eq. (3)), so that each species sequence is a series of nucleotides which follow another tetranomial distribution identically and independently. The state space is {ATCG} and the state proportions are (P A , P T , P C , P G ) which can be calculated by Thus, each species' nucleotide follows the same tetranomial distribution as the ancestor nucleotide. Under the independence assumption, the probability for the species sequence is simply a product of all nucleotide marginal probabilities, (1/4). This formulation can also be used to sample the unknown ancestor state among {ATCG} given the reference species state j = (A, T, C or G), since the posterior distribution among {ATCG} for the unknown ancestor state can be easily derived to be We now briefly investigate the ambiguity extent to which different sources of sequence are aligned. For simplicity, we use the Jukes and Cantor 8 model and assume the ancestor vs species nucleotide identity ("ancestor = species") probability is which equals p t ii ( ), in Eq. (3), the substitution probability is p, which equals p t ij ( ) in Eq. (3) for i j ≠ . The identity probability between two species ("species = species") nucleotides with equal evolution duration is thus The statistical sequence evolution model works on probabilistic transition from the ancestor nucleotide to species nucleotide. Since the ancestor sequence is never known for a direct alignment, we may sample it a posteriori given the reference species nucleotide. The probability for event "X", nucleotide identity between such a posterior ancestor nucleotide and another species' nucleotide ("posterior ancestor = species") other than the reference species nucleotide, is derived as The three identity probabilities (Eqs. (6), (7), and (8)) are plotted in Figure 2 where the ancestor-to-species transition probability (p) varies. We find that, the identity probabilities for these three types of matched nucleotides follow the order Pr(ancestor = species)  Pr(species = species)  Pr(posterior ancestor = species), or (6)  (7)  (8). On calculating the probability of a set of orthologous sequences Dovepress submit your manuscript | www.dovepress.com

Dovepress
Thus, alignment between the ancestor and species sequences may be less ambiguous than between two species sequences. Figure 2 also indicates that, the difference among these three types of identity probability is most significant in the middle of interval (0, 1/4). The dominance of "Pr(ancestor = species) over Pr (reference species = another species)" is more significant than that of "Pr (reference species = another species) over Pr (posterior ancestor = another species)", and the difference between the lower two curves ( Figure 2) seems not to be relatively large. Thus pairwise alignment between the posterior ancestor and another species sequence may achieve similar unambiguity as alignment between species. Now we study a multispecies orthologous sequence set (say Human = Species 1, Mouse = Species 2, and Rat = Species 3 in Figure 1). We denote by B a , B 1 , B 2 , and B 3 the sequences of the unknown ancestor, Species 1, 2, and 3. Under the nucleotide substitution model and unambiguous matching, the probability for the set of sequences under parallel evolution is where π ( ) ⋅ is the identical and independent tetranomial distribution for the ancestor nucleotide with state space {ATCG} and equal (1/4) proportion. The result is obtained by integrating out four possible ancestor nucleotides on each site for a marginal nucleotide group (three members across species) probability and multiplying these individual marginal probabilities along the sequence. Similarly, the phylogenic evolution model requires integrating out both the ancestor and root nucleotides on each site to get the result. Note that multiple alignment is not needed under the substitution model since no gaps are allowed. For general nucleotide substitution-insertion-deletion, the probabilistic evolution model developed by Rivas 9 gives the overall "substitution, insertion, and deletion" probabilities from the ancestor to species given divergence time. Calculating the evolution probability from the ancestor (assumed to be known) to the observed species sequence using the Rivas 9 model may require multiple alignments up front in order to match those nucleotides between the ancestor and species. Aside from not knowing the ancestor sequence, unambiguous alignment may not exist due to moderate sequence divergence. 7 Thus, one can underestimate the sequence set probability which is induced in a similar way to Eq. (9) because it simply picks one alignment to calculate the sequence set probability without incorporating other possible alignments. Ignoring ambiguous alignment may also lead to incorrect phylogenic inference and/or misleading sequence taxa partition pattern. 10,11 Under a Dovepress moderate sequence length (∼100 nucleotides), a tetranomial distribution for each nucleotide along the ancestor sequence may be used to sample ancestor nucleotides independently to form a large number (n) of sequences which are further used to induce a probabilistically evolved set of species sequences based on the Rivas 9 model. The sequence set probability is simply estimated as the number of exact duplicates of the given sequence set divided by n. However, this is highly impractical under moderate sequence length due to the small chance of duplicate sequence sets. Another possible way is applying pairwise alignment between each sampled ancestor sequence and observed species-specific sequences, and the sequence set probability may be done by averaging these evolution probabilities over all sampled ancestor sequences. However, this may also be inefficient due to noninformative ancestor sampling and lack of reliable alignment between a random sequence and species sequences. Thus it becomes desirable to propose and investigate more efficient multiple-imputation-like approaches such as using posterior ancestor samples which may offer multiple representative alignment results conditional on a reference species sequences for sequence set probability elicitation. Instead of following the theme of Eq. (9), we turn to calculating P (B 1 , where [II] is obtained after multiplicity over all nucleotide marginal probabilities for Species 1 (see Eq. (4)). As for [I], since B a (posterior ancestor) is sampled from the reference sequence (offspring) B 1 and the integrand is the offspring (B 2 , B 3 ) probability derived from the representative ancestor B a which is already linked to offspring B 1 through posterior sampling. Monte Carlo integration introduced in Eq. (10) realistically implements the joint probability of multiple post-evolution sequences by working on pairwise alignments between the sampled ancestor sequence and observed species-specific sequences. Under phylogenetic (tree-structured) evolution (the right panel in Figure 1), the sequence set probability can be written as where B 1 is Species 1 (Human) sequence, B r is the root sequence, and B a is the ancestor sequence in Figure 1. However, if we use Species 2 (Mouse) sequence B 2 as the reference sequence, then we have the following decomposition Note that, As in Eq. (4), we assign 1/4 to the probability for each nucleotide along the reference sequence B 2 after applying Eq. (13). Only pairwise alignment between the posterior ancestor sequence and species sequence is used for Monte Carlo integration (Eqs. (10), (11), and (12)). Since the probability of a sequence evolving from an ancestor is obtained by multiplying over all individual nucleotide evolution probabilities along a sequence, a large sequence length (say 100) may result in an overly small probability and lead to numerical overflow. The log-probability (LogPr) for a species On calculating the probability of a set of orthologous sequences Dovepress submit your manuscript | www.dovepress.com Dovepress evolutionary sequence from the ancestor is the summation of individual nucleotide evolutionary log-probabilities, and the evolutionary probability expectation obtained from Monte Carlo integration (exp(LogPr) mean), can be implemented by using moment generating function with argument one. Normality of these randomly produced LogPrs leads to the simple result of exp (µ + σ 2 /2) where µ and σ 2 are the sample mean and variance for these LogPrs.

simulation and real data study
We first introduce in detail the extended Jukes and Cantor model by Rivas 9 which will be used for our simulation study. The transition probabilities among general states {-ATCG} ("-" is the gap or covalent bond between two nucleotides) until time t are

s t r t t s t s t r t s t t s t s t s t r t
where the {-ATCG} column to the left of the transition probability matrix represents the initial (ancestor) states and the {-ATCG} row on top of the matrix represents the final (species) states. Specifically,

Parallel evolution model
We refer to the left panel of Figure 1. 1. Simulate the ancestor sequence with length = L 0 ; 2. Simulate species "1, 2, 3" sequences from this simulated ancestor sequence; 3. Apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 2 and 3 conditional on Species 1 sequence, where the unknown ancestor sequences are sampled using Eq. (16) with corresponding divergence time; 4. As a numerical verification, we apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 1 and 3 conditional on Species 2 sequence, where the unknown ancestor sequences are sampled using Eq. (16) with corresponding divergence time; 5. As another numerical verification, we apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 1 and 2 conditional on Species 3 sequence, where the unknown ancestor sequences are sampled using Eq. (16) with corresponding divergence time; 6. We investigate the consistency among different references. 7. Various divergence time vector (t 1 , t 2 , t 3 ) in the left panel of Figure 1 and transition parameter (β and q 0 in Eq. (15)) configurations are given in Table 1, where transition parameter (α in Eq. (15)) is standardized into one unit.

Phylogenic evolution model
We refer to the right panel of Figure 1. 1. Simulate root sequence with length = L 0 and the evolved ancestor sequence for Species 2 and 3; 2. Simulate Species 1 sequence from this simulated root sequence, and simulate the Species 2 and 3 sequences from this simulated ancestor sequence; 3. Apply Monte Carlo integration to randomly produced log(evolution probabilities) for Species 2 and 3 conditional on Species 1 sequence, where the unknown root and ancestor sequences are sampled using Eq. (16) with corresponding divergence times; 4. As a numerical verification, we apply Monte Carlo integration to randomly produced log(evolution probabilities) for Carlo integration to randomly produced log(evolution probabilities) for Species 1 and 2 conditional on Species 3 sequence, where the unknown root and ancestor sequences are sampled using Eq. (16) with corresponding divergence times; 6. We investigate the consistency among different references. 7. Various divergence time vector (t 0 , t 1 , t 2 , t 3 ) in the right panel of Figure 1 and transition parameter (β and q 0 in Eq. (15)) configurations are given in Table 1, where transition parameter (α in Eq. (15)) is standardized into one unit. We use the same transition parameter (β and q 0 in Eq. (15)) from the parallel model simulation and make evolution divergence times comparable between the parallel and phylogenic models. We collect LogPrs from 1000 Monte Carlo simulations. The distribution of these LogPrs are plotted in Figures 3, 4, 5 and 6. A Kolmogorov-Smirnov normality test gives p-value (0.15) for all LogPr sets, which means that the difference between the produced LogPrs and a normally distributed random variable is not significant. The normality assumption for LogPrs holds and the probability approximation based on log-normal distribution is reasonable. For such an assumption, a heuristic justification without rigorous theoretical proof is as follows: Given each randomly produced ancestor sequence, each nucleotide (event) LogPr on the non-reference sequences acts as an independently and identically distributed random variable, and the summation of these LogPrs follows the central limit theorem for a large sample size (sequence length). By increasing the ancestor or root length from 100 to 500, we can see that the relationship between LogPr and the sequence length is approximately linear. Another observation from Tables 2 and 3 is that different reference species sequences may lead to inconsistent sequence set probabilities due to different evolution durations and/or topological locations within the phylogenic structure. The phylogenic evolution model (the right panel of Figure 1) seems to show more inconsistency than the parallel evolution model (the left panel of Figure 1) does due to the dual missing sequences (the root and ancestor) instead of ancestor only in the parallel evolution model. A reference sequence which is closer to the root and/or ancestor is preferable since the imputed multiple roots and/or ancestors tend to be more informative due to shorter divergence. The CLUSTAL W multiple alignment 12 -induced probability is obtained by moving along the sequence set which holds nucleotides {ATCG} and possible gaps (covalent bonds) and applying Rivas 9 model, where "one gap with two nucleotides" across three matched sequence sites stands for a deletion and "two gaps with one nucleotide" across three matched sequence sites stands for an insertion. For each simulated sequence set, the discrepancy between Monte Carlo integration (MCI) and single multiple alignment (MA) induced probabilities are clearly more significant than that among probabilities estimated from different reference species sequences.

CREB promoters study
From the ABS database, 13 we extracted the promoter regions of transcription factor CREB for three mammals (human,    On calculating the probability of a set of orthologous sequences Dovepress submit your manuscript | www.dovepress.com Dovepress mouse and rat). We used MEGA 4.1 package 14 to construct the phylogeny tree with corresponding divergence times under uniform transition rate 1 (see Figure 7). These are used for sampling the posterior ancestor and root.
Since no current packages give us β and q 0 maximum likelihood estimation for Rivas 9 model (Eq. (14)), we mainly investigate the sequence set probability sensitivity to β and q 0 input by trying different values. We report the means and variances as well as estimated Log(sequence set probabilities) under two parameter settings for β and q 0 under different reference species. We also report the p-values from Kolmogorov-Smirnov normality test ( Table 4). The normality test results are sensitive to parameter input and reference species selection, which may be due to the fact that conservation levels/transition probabilities are likely to be nonhomogeneous along the sequences. Two LogPr distributions under associated parameter inputs are plotted in Figure 8. As a verification, we apply multiple alignment to these three promoters and at each site we calculate the nucleotide identity proportion within the window (with size 23) starting from this site at the gene direction ( Figure 9). The conservation levels show that evolutionary transition rates  Dovepress are approximately constant piece-wisely, thus the central limit theorem discussed for the simulated data study may still apply to these LogPrs on each promoter segment with quasi-constant conservation level under certain nucleotide insertion-deletion parameter (β, q 0 ) values.

Discussion
We proposed and investigated some promising numerical algorithms for accurately estimating the probability of a set of orthologous sequences with equal length under certain assumptions. Our approach was to informatively shuffle the unknown ancestors and/or roots and to find the distributional characteristics of simulated log-probabilities in order to reasonably approximate the true probability. The merit of our approach depends on how well the ancestor and/or root is imputed based on certain pentanomial distribution proportions (p -, p A , p T , p C , p G ) in Eq. (16) using the evolution model 9 and how reliably the pairwise Needleman-Wunsch alignment is applied to cross-species matching of nucleotides which are supposed to come from the same ancestor entry {-, A, T, C or G}. The former depends on the divergence duration from the ancestor/root to the reference sequence and the latter may depend on the species-specific adjustment of pairwise alignments based on phylogenic information.
When this piece of information is not immediately available, the algorithms by Yang, 15 Redelings and Suchard, 11 and MEGA package 14 are useful. Recently, Wong and colleagues 16 demonstrated that various alignments may lead to quite inconsistent inference. Although distance estimation for multiple species from a common ancestor may lack some accuracy using only one sequence set (Figure 9), we used MEGA package for phylogenic structure information for real sequence set probability estimation. Note that we only use background sequences as examples to demonstrate our algorithms by assuming independent tetranomial distribution among {ATCG} along sequences. For the set of orthologous sequences involving many species (3), we follow the evolutionary process (described by a phylogenic tree) to sample the internal nodes within the phylogenic tree conditional on one selected reference sequence (a terminal node on the phylogenic tree) and apply Monte Carlo integration to these imputed internal nodes for obtaining LogPrs (we omit the details). As one referee points out, it may be unreliable to directly apply our algorithms to sequences with very irregular lengths, since the insertion-deletion events need to be identified by matching nucleotides across all involved species other than due to artificial sequence truncation. Thus a crude multiple alignment across such

Publish your work in this journal
Submit your manuscript here: http://www.dovepress.com/advances-and-applications-in-bioinformatics-and-chemistry-journal Advances and Applications in Bioinformatics and Chemistry is an international, peer-reviewed open-access journal that publishes articles in the following fields: Computational biomodelling; Bioinformatics; Computational genomics; Molecular modelling; Protein structure modelling and structural genomics; Systems Biology; Computational Biochemistry; Computational Biophysics; Chemoinformatics and Drug Design; In silico ADME/Tox prediction. The manuscript management system is completely online and includes a very quick and fair peerreview system, which is all easy to use. Visit http://www.dovepress.com/ testimonials.php to read real quotes from published authors.

48
Liu et al Dovepress submit your manuscript | www.dovepress.com Dovepress Dovepress sequences may overly produce insertion and/or deletions. A rough solution may involve first applying multiple alignment procedures to these sequences and then segmenting the aligned sequences into subsequences involving different numbers of species followed by segment-wise Monte Carlo integration. However, the internal edge-effects introduced by segmentation deserves further study. Lastly, we highlight that applying the proposed algorithms to real sequences is not so straightforward in view of heterogeneous conservation patterns along the orthologous sequences, which poses as an important future research topic.