Proceedings 21.Wissenschaftlich-Technische Jahrestagung der DGPF, Konstanz, 4-7 September 2001, pp. 289-301
Automatic Reconstruction and Visualization of a Complex Buddha Tower of Bayon, Angkor, Cambodia
A. Gruen1, L. Zhang1, J. Visnovcova1
Abstract: The Hindu and Buddhist monuments of Angkor in Cambodia are the legacy of highly developed Khmer empires. Well-known for their structural and surface complexity they constitute a great challenge to any attempt towards precise and detailed 3D measurements and modeling. This paper reports about a pilot project using modern techniques of analytical and digital photogrammetry to derive a photorealistic 3D model of one of the very complex towers of the famous Bayon temple of the ancient city of Angkor Thom. This high quality model will then be subject to visualization and animation. On occasion of a balloon photogrammetry mission over the Bayon temple of Angkor , the first author took a number of tourist-type terrestrial images with a Minolta Dynax 500si analogue SRL camera of one of the many Buddha-faced towers of Bayon. We aim at deriving automatically, after scanning of the images, a texture mapped 3D model of the very complex object. In a first step we have generated already such a model with a mixture of manual (phototriangulation) and automated procedures (image matching for surface reconstruction, editing for blunder removal, texture mapping, visualization and animation). This result has been presented in Visnovcova et al., 2001 and includes already some novel methods for point cloud blunder editing and view-dependent texture mapping. Now we will include more images than before, intend to automate the phototriangulation procedure, develop a matching approach based on MPGC (Multiphoto Geometrically Constrained) Matching and also improve the procedure of color texture mapping. With these combined new approaches we expect a higher level of completion of the object, an improvement of the modeling quality and a shortening of the processing time to generate the 3D model. We streamline our new method of point cloud editing and improve our new technique of view-dependent texture mapping. Finally we show how different forms of high quality visualization techniques can be used on low-cost platforms. The ultimate goal of this activity and investigation is to develop a system that is capable of producing high quality photorealistic 3D models in fully automated mode.
Between the 9 and 15 centuries A.D. a succession of 42 kings ruled over one of the most remarkable empires in human history: The Khmer empire of Angkor in Cambodia.What used to be before and after an endless stretch of tropical jungle was turned into 75 square miles of fertile plains, with 72 major monuments, temples, palaces, canals with dikes, moats, and reservoirs. A sophisticated irrigation system, controlling equally abundant waters from monsoon rains and yearly droughts, provided for two to three rice harvests per year, feeding at times up to one million inhabitants. Today it is one of the most spectacular sites listed in the UNESCO World Heritage List (DAGENS, 1995, MALRAUX, 1930, MOUHOT, 1864, MOORE, 1960, PRESTON, MCCURRY, 2000). While statutes and other treasures were shipped off from the site from the very beginning of European interference, extensive excavations and renovation works have been under way
Institute of Geodesy & Photogrammetry, ETH-Hoenggerberg, CH-8093 Zurich, Switzerland, {agruen,zhangl,jana}@geod.baug.ethz.ch

1 Introduction th th

only since the early years of the 20th century. Renovation and conservation of the laterite, brick and sandstone monuments is extremely difficult under the local conditions of heat, heavy rains, humidity and an aggressive tropical vegetation. Today, after a decade of forced neglect during the rule of the Khmer Rouge, the precious Angkor monuments receive much help again from various countries. Among those is the JSA project (Japanese government team for Safeguarding Angkor), which aims at the conservation and restoration of the Northern Library of the Bayon temple in the ancient city of Angkor Thom (JSA NEWS, 1999). Bayon was the last great temple built in Angkor and remains one of the most enigmatic. Today Bayon is considered an interesting mixture of both Hindu and Buddhist elements of style. Bayon is situated right in the center of Angkor Thom. It is an extremely complex structure, which is best revealed by aerial photographs (compare Fig. 1). It is said to represent the Mount Meru of an old Hindu legend. It consists of three levels of platforms and two galleries with spectacular and well-preserved basreliefs, showing the live of the ordinary people and stories from the Hindu mythology. A special feature of Bayon are its 54 towers with four large faces each, pointing in all four geographical directions (such adding up to 216 faces in total, not counting the many more on other objects). These faces on the towers and gateways of Angkor Thom represent Lokesvara, a Bodhisattva from the Mahayana Buddhism, a holy one who stayed at earth to do good work and help people (Fig. 2). Some experts still think that they are stylized portraits of Jayavarman VII, the most prominent of all Khmer kings.
Fig. 1: Aerial image of Bayon, taken from a Helium balloon
While on sabbatical leave at the Asian Institute of Technology (AIT) in Bangkok in spring 1999 the first author was invited to join the Japanese team of geographers and surveyors from the Keio University, Tokyo as photogrammetrist for a balloon photogrammetry mission over Bayon. The main purpose of this mission was to generate a 3D-computer model of the very complex structure of Bayon from small format balloon images (Fig. 1). The results of this work will be published elsewhere. During this one-week mission the author took a number of film-based amateur photographs of one of the many Bodhisattva-faced towers. This paper reports about the processing of these images and shows the resulting texture mapped 3D model of the tower. In a previous paper (VISNOVCOVA ET AL., 2001) we described our semi-automated approach, which consisted of manual interior orientation, phototriangulation and some manual input into the joint registration of the individual surface models derived from single stereopairs. Since the aim of our recent work was on fully automated processing of the image sequence we had to develop some novel methods for phototriangulation and point cloud editing. We will report on all stages of the processing chain: Interior orientation, phototriangulation, image matching for surface model generation, point cloud editing and texture mapping.

2 Data acquisition

The main goal of the Bayon field campaign was the taking of small format-camera balloon images over Bayon for the 3D reconstruction of the complete and utterly complex temple (Fig. 1). During the mission a sequence of terrestrial small format images was taken with a Minolta Dynax 500si SRL camera (c = 35 mm) of one of the Bayon towers on the third platform in the southeast corner of the temple (Fig. 2). The sequence of 16 images, covering the full horizon, was rather meant as a lab exercise in photogrammetry than a serious project aiming at a complete recording of the object. Fig. 3 shows the arrangement of images, as they were used for bundle triangulation around the tower. Since a 360 degrees azimuth coverage was necessary and the light conditions were very extreme, the production of good, evenly illuminated pictures was practically impossible without artificial lighting. Therefore, the images suffer under strong variations of the Fig. 2: The smiling faces of Bayon. illuminated and shadow areas. Also, the shadow and One of the Minolta images used for light parts vary from image to image, depending on processing the time of the day the images were taken. This will cause problems with texture mapping from multiple images. Thus a good deal of our work went into a modification of the standard procedure of texture mapping.
Fig 3. Arrangement of images

3 Data processing

The data processing chain consists of the components interior orientation, phototriangulation, and surface reconstruction, point cloud editing and texture mapping. Fig. 4 shows a flow chart of the processing modules. We will report on these modules in detail, but focus on the novel approaches that we developed. In a first step, the images were scanned on an Agfa Horizon with a resolution of 21 m.

3.1 Interior orientation

Since the amateur images do not include fiducial marks we determine "virtual fiducials" by intersecting the straight lines of two adjacent and mutually perpendicular image edges. This procedure is applied to all images of the sequence and allows us to establish individual image coordinate systems, which are used for the transformation of the pixel coordinates. Since all pictures have been taken at constant focus, the camera constant is set to the same nominal value for all images. Since the lens distortion is not known at this point it is assumed to be zero. We will later estimate the systematic errors by using self-calibration in bundle triangulation.

Digitized Images data Additional Information
Interior Orientation and Automatic Relative Orientation

Epipolar Geometry

Disparity Surface
Semi-automatic Tie Point Extraction by ALSM Bundle Adjustment with Self-Calibration Fully Automatic 3D Surface Reconstruction by MPGCM
Exterior Orientation Parameters
3D Model Surface Data Editing and Triangulation View-dependent Texture Mapping
Visualization of Textured 3D Model
Fig. 4: Processing modules

3.2 Phototriangulation

The object coordinate system is established with minimum datum by measuring an object distance on a vertical leveling bar, and defining another point as third depth control point. While the distance provides for scale, the vertical bar also defines the direction of the Y-axis. The orientation of the X- and Zaxes is arbitrary in our kind of problem (Fig. 5). The phototriangulation is set up as a semi-automated procedure. It is

Y 2 [1.00; 5.05; 10.00]

1 [1.00; 1.80; 10.00] X 3 [-; -; 10.00] Z
Fig. 5: Definition of the object coordinate system via minimum datum. The coordinates are given in meters
based on manual tie point selection and Least Squares based Geometrically Constrained Multi-Image Matching. Since Least Squares Matching is a nonlinear procedure the success rate of matching depends largely on the reduction of the matching search spaces. This is achieved in two steps, executed prior to tie point matching: - Automated relative orientation of adjacent images for epipolar geometry recovery - Determination of disparity surfaces in all image pairs The first step allows for an epipolar-guided search and the second step reduces the search range along the epipolar lines. 3.2.1 Automated relative orientation Our approach for relative orientation considers the particular nature and orientation of these close-range images and can be characterized by the following detailed steps: (a) Feature point extraction by Moravec operator. This is done independently for each image of an adjacent stereopair. A good point distribution is assured by image subdivision in small patches and only keeping the points with the highest interest values in those patches. (b) Computation of correlation coefficients for all point combinations using 7x7 image windows, a search window size of 1/16 of the overall image size and a correlation threshold of r = 0.75. (c) Conjugate point determination by relaxation, following the approach of ZHANG, ET AL., 1995,1996, XU, ET AL., 1996. The energy function used in this relaxation procedure is based on continuity and uniqueness constraints. In our approach we apply the "looser-takes-nothing" iterative strategy for the minimization of the energy function. This procedure is essentially a combination of area-based and feature-based matching. (d) Derivation of the epipolar geometry by a linear method. Since we have to deal at least in parts with highly convergent images we apply the Direct Relative Transform (DRT) for relative orientation (CHANG, 1986). With the image coordinates ( x, y ) and ( x, y ) of the left and right images respectively and the joint camera constant c the DRT formulation results in L1 xy + L2 y y L3 y c + L4 xc L5c 2 + L6 xx + L7 xy L8 xc + cy = 0 (1) Where L1 ,., L8 are the eight unknown parameters, requiring at least eight conjugate image points. The orientation parameters are derived from these parameters by algebraic decomposition. Since we have to expect many blunders in the conjugate point pairs we use redundant data and the robust estimation method of minimization of the median of squared residuals (LMedS) min(med ri 2 ). Equation (1) is used as observation equation.

The chosen estimator must yield the smallest value for the median of squared residuals computed for the entire set of conjugate points. We are not aware of the existence of a direct solution for the LMedS estimator. Therefore a Monte Carlo technique is used to solve this problem, based on the choice of m = 1024 random sub-samples of p = 8 conjugate point pairs each. For each sub-sample j we determine the parameter vector r r L j = ( L1 ,., L8 ) j =1,.,m. For each L j we obtain the median of the squared residuals according to M j = med ri2 = med ( L1xy + L2 yy L3 yc + L4 xc L5c2 + L6 xx + L7 xy L8 xc + cy)2 (2) i =1,.,n i =1,.,n r We retain the solution L j for which Mj is minimal. Algorithmic details of this method will be published in a follow-up paper. Since the performance of LMedS is poor in case of presence of Gaussian noise we finally conclude our estimation with a Weighted Least Squares adjustment, using the weighting scheme 1 if ri 2 (2.5 ) 2 wi = (3) 0 otherwise
5 ) min ( M j ) (4) j n p Observations with wi = 0 are regarded as blunders and as such effectively excluded from the Least Squares solution min( wi ri 2 ) (5) with = 1.48(1 +
It turns out that this method is very robust to data sets containing blunders. In our experiments, the epipolar geometry can be robustly recovered from the matching results containing 5 - 40% blunders using the LMedS method. 3.2.2 Derivation of disparity surfaces In our application with 16 images, spread over the horizon, we have 16 image pairs, with some convergence angles larger than 45 degrees. In 11 models we get about 230 to 500 conjugate points, which give us a nice disparity surface. Another 4 models contain less than 20 conjugate points because of large convergence angles and as such big differences in local image scales. One model has such an unfavorable geometry (large base component perpendicular to the disparity direction, a so-called "vertical base"), that automated image matching did not work because of big image scale differences. This model had to be measured manually.
(a) The reference patch with marked point
(b) The approximate matches (arrows) and epipolar lines
(c) The final matching result of ALSM Fig. 6: An example of our tie point measurement procedure

3.2.3 Tie point measurement for phototriangulation The image measurement for triangulation was performed with Adaptive Least Squares Matching (GRUEN, 1985) in a semi-automated mode. The operator selects feature points in one image and ALS matching is done in an unconstrained mode since the current exterior orientation values can only be regarded as approximations. Since ALS matching is nonlinear and has a relatively small convergence radius we need fairly good approximate values. These are obtained by epipolar constrained cross-correlation. The epipolar geometry is taken from the previous relative orientation of the image pairs and also the disparity surfaces are considered as constraining conditions. This procedure delivers for the complete image arrangement 181 tie points, with (18,63,81,15,4) points connecting (2,3,4,5,6) images, thus providing for sufficient strength in image connection. Fig. 6 shows an example for tie point measurement. Here the respective approximate matches, the epipolar lines and the final ALS matching results are displayed. 3.2.4 Bundle adjustment Prior to the bundle adjustment, a strip adjustment with 16 models is performed to obtain the approximate values for the exterior orientation elements and the object coordinates of the tie points. This strip adjustment is done in model coordinate space and in a subsequent step the model coordinates are transformed to the object coordinate space by a 7 parameter similarity transformation using the 7 control point coordinates. For bundle adjustment we use our own software SGAP. In a first step, we apply SGAP to obtain the absolute orientation parameters for all images and the 3D object coordinates of all tie points without self-calibration, to guarantee global convergence. After the establishment of a stable adjusted image block, ten additional parameters (BROWN, 1971) are used to model systematic errors: the camera constant correction, principal point coordinate offsets, five parameters modelling the radial and tangential lens distortion and two parameters for a affine scale factor and shear (BEYER, 1992). In our case, the principal point coordinate offset in x direction, the parameter of the correction of camera constant, the affine scale factor and the first term of the radial lens distortion turned out to be significant. The SGAP bundle adjustment reported the following values:
Total number of tie points: 181 Total number of observations n = 1250 Global internal reliability r/n = 0.52 Standard deviation of image coordinates a posteriori: 0 = 20.3 micron Theoretical precision from 181 tie points (in mm):

X = 5.8, Y = 3.7, Z = 8.8
3.3 Surface reconstruction
Our surface reconstruction approach works in fully automated mode according to the following procedure: (1) Selection of 5 images to serve as master images. (2) Use of Moravec operator in each master image to extract a very dense pattern of feature points. In our procedure, the master image is subdivided into 77 small image patches and within each patch only one feature point which gains the highest interest value is selected.
(3) For each feature point we use the improved epipolar geometry, determined in phototriangulation, and the disparity surfaces between adjacent images to get good approximations for the following MPGC (Multi-Photo Geometrically Constrained) matching procedure by standard cross-correlation technique. (4) In a final step MPGC is applied for fine matching, including patch reshaping. MPGC exploits a priori known geometric information to constrain the solution and simultaneous use of more than two images (GRUEN, 1985, GRUEN, BALTSAVIAS, 1988, BALTSAVIAS, 1991). This algorithm can use some well-known features such as blunder and occlusion detection, quality control and self-diagnosis. In our procedure, for each feature point in the master image, we combine 6 adjacent images (totally 7 images) for the MPGC matching algorithm with automatic blunder and occlusion detection. With MPGC we can get sub-pixel accuracy matching results and their 3D object coordinates simultaneously. This procedure results in fairly reliable and very precise matching results. We obtain 5 point clouds with a total of 88140 points, which however are already jointly registered. Therefore, the matching of adjacent surfaces, usually a very errorprone procedure, is not necessary any more. Although we use the automatic blunder and occlusion detection, there are some blunders left in the final 3D point cloud. We apply the following surface data editing procedure to remove these blunders.
3.4 Surface data editing and triangulation
Determining surfaces from a set of 3D points containing outliers is a complex task. Many approaches have been designed to treat this problem, such as 3D deformable surfaces (COHEN ET AL., 1991) and iterative local surface fitting (FUA and SANDER, 1992). The last algorithm fits second order patches around each 3D point and groups points into surfaces. In this iterative process errors are eliminated without smoothing out relevant features. Our procedure is an adapted version of this algorithm without performing resampling and surface clustering. Because all the points should belong to one unique surface, only those outliers are deleted whose derivation from the fitted surface exceeds after 3-5 iterations the predefined threshold. Finally the errors are eliminated while preserving essential surface features (Fig. 7).

Fig. 7: The shaded triangulated point cloud before (a) and after editing (b).
For the conversion of the point cloud to a triangular surface mesh the 2.5D Delaunay triangulation was applied. Without losing its topology, the 3D surface model of the Bayon Tower was expanded to a plane by transforming the Cartesian coordinate system to a cylinder coordinate frame. In the defined z cylinder frame z is the vertical cylinder axis crossing the model center and parallel to the original Y-axis of the Cartesian object coordinate system. is the Euclidean distance from the surface point to the z-axis and is the angle around the zaxis. The 2.5D Delaunay triangulation was done in the z plane. The final shaded model of the triangulated mesh is shown in Fig. 7b.
3.5 View-dependent texture mapping
With the technique of texture mapping, gray-scale or true color imagery is mapped onto the 3D geometric surface in order to achieve photorealistic virtual models. Knowing the parameters of interior and exterior orientation, to each triangular face of the 3D surface the corresponding image coordinates are calculated. The gray-scale or color RGB values within the projected triangle are then attached to the face. A common approach of texture mapping is to use one frontal image for a related part of the object. In close-range applications this is often not satisfactory, because not enough image information is available for fully or partially occluded object parts. In the 3D model the texture appears as stretched (Fig. 8). Moreover, often-varying lighting conditions during image acquisition do not allow regular light distribution allover the object. This causes sharp transitions between neighbouring object parts, which are texture-mapped from different frontal images. To overcome these problems, a new method of texture mapping was developed a view-dependent texture mapping. Our method of view-dependent texture mapping is based on the selection of a combination of optimal image patches for each triangle of a 3D model. According to the best possible geometric and radiometric conditions a combination of image content is calculated from all images where a particular triangle appears. The locations of the image triangles are computed from object faces via collinearity equations. The procedure consists of three steps: preprocessing, selection of a geometrically optimal image and texture weighted averaging.

triangle t image plane p2

projection centre2

image plane p1

projection centre1

Fig. 8: Texture mapping with one frontal image sloped faces A and B appear with stretched pixels
Fig. 9: View-dependent texture mapping: for each object triangle the optimal texture image is selected (image 1 with the bigger area of projected triangle t)
(1) Pre-processing. In order to achieve a better local contrast and nearly the same color balance in all images high pass filtering and histogram equalization are used for each RGB channel separately. Similarly to Wallis filtering, the high pass filter enables a strong enhancement of the local contrast by removing low-frequency information in an image. It retains edge details in the specified radius where sharp color transitions occur. We apply high pass filtering with a radius of 50 pixels. Additionally, the equalization is performed
for all three channels in all images. With these procedures we avoid color shifts when merging the separate channels together. (2) Selection of a geometrically optimal image (Fig. 9). In all image planes where a particular triangle t is projected, the area of the triangle is calculated. The image where the triangle t appears largest contains the most texture information. Thus, it is considered as a geometrically optimal image candidate for texture mapping. However, if we use only one image patch in case of large differences in image radiometry and in case of strong local variability of the surface patch normal, the result will be a checkerboard-type 3D texture map. Therefore, in our procedure the radiometry of adjacent images related to the selected optimal image is also considered. (3) Texture weighted averaging. We consider all images where a particular triangle t appears simultaneously. The gray values of these images are averaged according to equation (6)

j =1 n

g i. gray value of the new image i g ij. gray values of old images j n. number of image patches used for texture weighted averaging wij. weight factors calculated from the areas of projected triangles wij = area j area max (7)
Obviously, the weights are chosen proportional to the area of the image patches, which gives priority to the image patch with the better geometry. Before this weighted averaging can take place the different image patches are transformed by an affine transformation to the geometry of the optimal master image. The weighted averaging reduces the effects of radiometric differences in adjacent images.

4 Visualization

The final model can be viewed with our own graphics program disp3D in wireframe, shaded or textured mode. Stereo display in form of anaglyphs and polarization is also feasible. Further, the models can be converted to VRML2 for viewing them with standard visualization packages. Some results are depicted in Fig.10.

5 Conclusions

We have presented a combination of algorithms that will finally allow us to generate a photorealistic high quality texture-mapped 3D model of a complex object in fully automated mode. Starting with a sequence of amateur images the various processing steps as scanning, interior orientation, relative orientation, tie point measurement, self-calibrating phototriangulation, multi-image matching for surface model generation, point cloud editing, view-dependent texture mapping and visualization are addressed. Novel approaches are specifically used in relative orientation, surface model generation, point cloud editing and texture mapping. The current system does not yet work fully automatically. The selection of tie point image patterns is still done in manual mode. However, we do not consider the automation of this step through interest point selection a particular problem and we will come up with a respective solution in the near future.
With an average coordinate accuracy of 6 mm the triangulation results are well within the usual specifications for such kind of project. It is difficult to give an accuracy figure for the surface model, but it can be expected to be roughly of the same accuracy. The overall results can be improved by the following measures: (a) Use of digital cameras. This saves scanning time and efforts for interior orientation and leads to a faster throughput. In our project such camera was not available on-site. (b) Generation of a denser sequence of images. This gives more network strength to triangulation, improves the reliability and precision of matching in relative orientation, tie point measurement and surface model generation and leads to better texture mapping results. (c) For image generation choose a better timing and/or artificial lighting. This will result in better texture maps. (d) Apply an efficient 3D Delauney triangulator. This will reduce problems in surface modeling (which was not a problem with our test object but could be with others). Our upcoming work will include some, if not all of these improvements and will hopefully lead to a robust and accurate system to solve such kind of problems within a minimal timespan and with very little operator interference.
(c) (d) Figure 10: Visualization of the Bayon Tower 3D model: (a) wireframe mode, (b) shaded mode, (c) textured mode of northern view, (d) textured mode of western view. Displays produced by disp3D.

6 Acknowledgement.

This is to acknowledge the support given by the JSA (Japanese government team for Safegarding Angkor) project. Above all we would like to thank Prof. Sachio Kubo, Keio University, Tokyo for giving the first author the opportunity to participate in this campaign. Also appreciated is the great help that we received in the field from Prof. Kubos young and ambitious team of surveying, climatology and economy students. Finally, our thanks go to Yoshito Miyatsuka for navigating and controlling the balloon even under the greatest difficulties, such that no harm was done to the precious instrument.

7 References

BALTSAVIAS, E., Multiphoto Geometrically Constrained Matching. Dissertation, IGP, ETH Zrich, Mitteilungen No. 49, Dec. 1991, 221p. BEYER, H. A., Geometric and Radiometric Analysis of a CCD-Camera Based Photogrammetric Close-Range System. IGP ETH Zrich, 1992. Mitteilungen Nr. 51, 186p. BROWN, D.C., Close-range Camera Calibration. Photogrammetric Engineering, 1971. 37(8), pp. 855 866. CHANG, B., The Formulae of the Relative Orientation for Non-metric Camera. International Archives of Photogrammetry and Remote Sensing, 1986, 26(5), pp.14-22 COHEN, I., COHEN, L., AYACHE, N., Introducing Deformable Surfaces to Segement 3D Images and Infer Differential Structure. Technical Report of INRIA, 1991. DAGENS, B., Angkor. Heart of an Asian Empire. Thames & Hudson Ltd., London, 1995, 191 pages. FUA, P., SANDER, P., Reconstructing Surfaces from Unstructured 3D Points. Proceedings of the Image Understanding Workshop, San Diego, California, 1991. GRUEN, A., Adaptive Least Squares Correlation: A powerful Image Matching Technique. South Africa Journal of Photogrammetry, Remote Sensing and Cartography, 14 (3), 1985, pp 175-187. GRUEN, A., BALTSAVIAS, E., Geometrically Constrained Multiphoto Matching. Photogrammetric Engineering and Remote Sensing, Vol. 54, No. 5, May 1988, pp. 633641. JSA News, JSA Project No 4, Commemorating the Completion of the Conservation and Renovation of the Northern Library of Bayon. JICE, Tokyo, Japan, 1999. MALRAUX, A., La voie royale. Editions Grasset & Fasquelle, Paris, 1930. MOUHOT, H., Travels in the Central Parts of Indo-China (Siam), Cambodia and Laos. 2 Vols., John Murray, London, 1864. MOORE, R., Angkor, Jewel of the Jungle. National Geographic Magazine, 1960, pp. 517569. PRESTON, D., MCCURRY, The Temples of Angkor Still Under Attack. National Geographic Magazine, August 2000, pp. 84 103. WHITE, P.T., The Temples of Angkor Ancient Glory in Stone. National Geographic Magazine, 1982. Vol. 161, No. 5, pp. 552 589.
VISNOVCOVA, J, GRUEN, A., ZHANG, L., Image-based Object Reconstruction and Visualization for Inventory of Cultural Heritage. Proceedings of EVA 2001 Conference Electronic Imaging and Visual Arts, 2001. Florence, Italy, 26.3.-30.3.2001, pp. 118-122. VISNOVCOVA, J, ZHANG, L., GRUEN, A., Generating a 3D Model of a Bayon Tower Using Non-metric Imagery. Proceedings of the International Workshop Recreating the Past-Visualization and Animation of Cultural Heritage. Ayutthaya, Thailand, 26.2.-01.3, 2001, pp. 30 39. XU, G., ZHANG, Z., Epipolar Geometry in Stereo, Motion and Object Recognition, A Unified Approach, Published by Kluwer Academic Publishers, 1996, ISBN: 0-79234199-6. 310p. ZHANG, Z., DERICHE, R., FAUGERAS, O., A Robust Technique for Matching Two Uncalibrated Image Through the Recovery of the Unknown Epipolar Geometry, Artificial Intelligence Journal, Vol 78, pp87-119, Oct. 1995. ZHANG, Z., Determining the Epipolar Geometry and Its Uncertainty: A Review. Research Report, No. 2927, INRIA Sophia-Antipolis, July, 1996.