RSD-TR-18-83 SOFTiWARE SUPPORT FOR THE STEREO CAMERA: PHANTOM AND MEASUREMENT PROGRAMS1 Gideon Frieder and Myoung Lee Department of Computer and Communication Sciences The University of Michigan Ann Arbor, Michigan 48109 November 1983 CENTER FOR ROBOTICS AND INTEGRATED MANUFACTURING Robot Systems Division COLLEGE OF ENGINEERING yTHE UNIVERSITY OF MICHIGAN ANN ARBOR, MICHIGAN 48109,', CSupported in part by USAF/SAM via SCEEE Contract No. F-33616-78-D-0317 Task 81 and by CRIM via.XFOSR Contract No. F49620-82-C-0089.

^ I F "-, a

TABLE OF CONTENTS 1. Introductio.................................................................................................... 2. The Phantom Program.................................................................................. 2 2.1. The Coordinate System............................................................................ 2 2.2. Program Data............................................................... 2.3. Set Up of the Projector and the Camera.................................................. 2.4. The Phantom Program Algorithm.............................................................. 5 2.5. Calibration...................6....................................................................... 6 2.6. The Phantom Program..................................... 7 2.6.1. ReadParameters.................8.................................... 8 2.6.2. Read_Cal_Para.......................................................................... 8 2.6... Iuitialize............................................................................................ 2.6.4. Find_ Tmat................................................. 8 2.6.5. Find Lmat......................9...................... 9 2.6.6. Find BeamEqs....................................9..................... 9 2.6.7. Find Ray Eqs................................................................................... 10 2.6.8. Find_BeamPlane ~ ~...............~............ 10 2.6.9. MaptoCamera.................................................................................. 11 2.6.10. Map.to_ Cam_Tmat................................. 6.. 11 2.6.11. Map..toCam_,Geo................. ~ ~ ~... 12 2.6.12. MakeImagePattern............................ 12 2.6.13. OutputPattern....................................................................... 13 2.6.14. Find CalFace..................................................................................... 13 2.6.15. Output_CubeImage............................................................... 14 2.6.16. Main Program....................................... 14 3. Surface Mapping Program............................................................................. 15 3.1. Init.................................................................................................................. 15 3.3. Imgin..........................~............................................... 16 3.4. Tmat............................................................................................................. 16 3.5. Lmat........................... 6........... 6 3.6. Map.............................................................. 17

3.7. Column............................................................................................... 17 3.8. Setcof...........17..............................................................................*~ 17 3.9. Space Coding................................................................................................. 18 4. The Analysis Program............................................................................................... 19 5. General Description of the Program Files................................................................ 20 6. Simulations..................................................................... 21 7. References.....~ ~......................................................................... 25 8. Program Listings.................................................................................................... 26 8.1. Phantom Program........................................................................................ 26 8.2. Eye.Cal Program Listing.............................................................. 27 8.3. Eye.Tlmat Program Listing.............................................................................. 28 8.4. Support Program Listing........................................................... 29 8.5. Analysis Program Listing................................................................................. 30 9. Sample Run................................................................................................... 31

1 1. Introduction This set of programs is designed to provide a test bed for the verification and testing of the stereo camera system[l,21. The stereo camera system can rapidly acquire the location of a large number of sample points on a surface. The system involves a coded-laser-beam-array projector with a programmable shutter, and a camera. The projector illuminates the surface with predetermined patterns. Points which are illuminated by the projector are imaged on the camera plane. The stereo camera testing program consists of three parts: a phantom program, a surface mapping program, and an analysis program. The phantom program simulates a measurement of a given surface. It calculates the points on the surface, called "surface points", which are illuminated by the projector, and tne points on the camera plane,cailed "image points", which are the images of the surface points. The surface mapping program reads the image points which were produced by the phantom program and calculates the surface points. The analysis program compares the surface points from the phantom program and the computed surface points from the surface mapping program, calculates the distance between the two and provides some error statistics. image point \ z ^^^ /~ ^ ~/ camera projector plane shutter plane / / / slface - - > /. point / Y x Figure 1. Stereo Camera System

2 2. The Phantom Program The phantom program simulates a measurement process and thus provides data to be used by the surface mapping program. It produces one file of surface points and another of image points. The phantom program also produces a file of calibration data for the calibration option of the surface mapping program. For a surface mapping program which does not calibrate, it provides T and L matrices calculated directly from geometrical formulas. 2.1. The Coordinate System The testbed programs use three coordinate systems. The world coordinate system is a three dimensional Cartesian coordinate system which is used to measure the position of the camera, the shutter and the surface to be mapped. Points on the surface illuminated by the projector and surface points which are computed by the surface mapping program are measured in this coordinate system. x, y, and z designate the three axes of this system. The shutter coordinate system is a two dimensional Cartesian coordinate system on the shutter plane of the projector. u and v are the two axes of this system. The camera coordinate system is a two dimensional Cartesian coordinate system on the camera plane. x, and y, are the two axes of this system. When positions and parameters of the shutter ( camera ) are known, a coordinate of a point on the shutter ( camera ) plane can be converted into a coordinate in the world coordinate system, and vice versa.

3 2.2. Program Data In the phantom program, the position and parameters of the shutter and the camera can be changed. In the current version, the surface to be mapped is always a plane, the equation of which can be changed. The size, the position, and the orientation of the calibration cube can be changed for the surface mapping program option which calibrates the T and L matrices. When the phantom program is invoked, it requests the following parameters from the user: For the projector Distsh, theta.sh, phish r 8,0 of the focal point of the projector D_sh Distance from the focal point of the projector to the center of the shutter plane Su, Sv Half width and half length of the shutter plane Nw, Nv Number of sample points in the whole shutter plane in u and v axis at present,N N, N, 64 For the camera Dist_ca, theta_ca, phica r 8,0 of the focal point of the camera D_ca Distance from the focal point of the camera to the center of the camera plane Sx, Sy Half width and half length of the camera plane N, Ny Number of sample points on the camera in x, and y, axis At present, N, Ny _ 128 and are powers of 2 For the calibration cube Length Length of the side of the cube Step Step in u and v axis. To control the number of the calibration data points, the calibration cube is illuminated by a projector at each shutter point ( u, v ) with "step" apart. Axis and Angle Rotation of the cube around the world "axis" xt, yt, zt Translation of the cube along x, y, z axis. For the surface to be mapped A, B, C, D Equation of the surface, Ax + By + Cz +D - 0

Figure 2. Parameters of the Phantom Program j focal \S Cpo D npoint 7/ I ______ ____ x N @~. 1 ~ X. sample Y. __ I _g. _____l __ points l I 1. * A__,_____ 1 2 N x

5 2.3. Set Up of the Projector and the Camera The projector and the camera are set up to aim at the same point, which is chosen as the origin of the world coordinate system. That means that the line perpendicular to the plane of the shutter (camera) passing through the center of the shutter (camera) surface and through the focal point of the shutter (camera), is also passing through the origin of the world coordinate system. Thus, for example, the shutter plane of the projector is perpendicular to the line which passes through the focal point of the projector and the origin of the world coordinate system.The orientation of the plane is such that two unit vectors a, a, are the same as the unit vectors in the spherical world coordinate system a, - a#. where 6, and i are the parameters of the projector. The camera is positioned in the same e way as the projector. 2.4. The Phantom Program Algorithm The phantom program requests parameters from the user. When all the parameters are given, computation begins. The phantom program first finds the world coordinates of each sample point on the shutter plane (called "shutter points") of the projector. For each shutter point, it finds the line equation of each beam which passes through the focal point of the projector and the shutter point. The phantom program then computes the intersection of the beam and the plane to be mapped, which provides a point which is imaged on the camera plane. The image points are calculated by three different methods. The first is the geometrical method. Here the equation of the line which passes through the focal point of the camera and the actual surface point is computed. Then the phantom program calculates the intersection of this line and the camera plane. That intersection point on the camera plane is rounded to the nearest sample point, which is the image point we want. The second method is the one employed in graphics. We first find the transformation matrix Te which transforms the world coordinate into the clipping coordinatel[3]: (xyx 1) = (xy zl)T,

6 where xc, Ye, and z, are the clip coordinates, and x, y, z are the world coordinates of the surface point. After it clips, the clipped coordinate is converted to the camera coordinate, which is the image point we want. And finally, we find the T matrix which transforms the world coordinate directly into the camera coordinate. The T matrix is computed by multiplying Te by the perspective transformation matrix Tp. The three methods are basically equivalent. The image points obtained by these methods are compared to ensure that computational errors due to roundoffs and catastrofic cancellations in floating point arithmetic do not provides erroneous results. For each shutter point (u, v) on the shutter plane, the phantom program provides the shutter points (u, v), image points (x,,y.), and the intersection (x, y, z) of a beam passing through the shutter point and the surface plane, on file unit 2. The phantom outputs image patterns on file unit 4. These patterns are determined by the space coding. When the surface mapping program does not calibrate the T and L matrices, the phantom program outputs the T and L matrices on file unit 3. When the surface mapping program calibrates the T and L matrices, the phantom program outputs the calibration data on file unit 3. 2.5. Calibration Calibration is the process of determining the T and L matrices[2J. These matrices are calculated by the least square method41], using data points which are measured on the surface of a known geometrical objects, in our case a cube. To determine the T and L matrices for various positions of the projector and the camera, the calibration cube should be properly positioned. Two or three faces of the calibration cube should be visible to both the projector and the camera, ( three is the maximum ), and there should be at least seven image points on the calibration cube. At present, the maximum number of calibration data points is 100. In the phantom program, it is possible to vary the size and to move the calibration cube. Positioning of the calibration cube involves rotation around the world coordinate x, y, z axes, which can be

7 followed by a translation along the x, y, z axes. Since the rotation operation is not commutative, it is possible to specify rotation in any order, e.g. rotating 81 around x axis, followed by 82 around z axis, followed by 83 around x axis, followed by 84 around y axis, etc. After completing the rotation, translation of the cube is optionally performed. Translation is specified by xt, yt, zt, which are the translation along the x, y, z axes respectively. By rotating and translating the cube, it is possible to position the cube in any location, in any orientation. The position of the calibration cube is important for the accurate calibration of the T and L matrices. When the size of the cube is entered, the calibration cube is located with its center at the origin of the world coordinate system, and its sides parallel to the x, y, z axes. With rotation and translation specified, a transformation matrix H and its inverse H-1 are calculated. The new position of the eight vertices and the new equation of the six faces are computed by multiplying them with the translation matrices, H and H', respectively[51. Among the six faces of the calibration cube, only the faces which are visible to both the projector and the camera are chosen. If the number of the faces visible to both the projector and the camera is one or zero, a different position of the calibration cube should be tried, since in this case the least square method does not provide a meaningful solution. When the number of the faces is two or three, intersections of the faces with the beams of the projector are calculated at each shutter point (u, v) with step apart, i.e. (u + m step, v + n step), where m, n s,, 2,... The intersections are imaged on the camera. The phantom program outputs the world coordinates of the intersections, (x, y, z), shutter column u, and the image points on the camera plane (x,, y8) on file unit 3. The number of the calibration data points is controlled by the parameter "step". 2.6. The Phantom Program The phantom program consists of one main program and the following procedures. Read_Parameters, Read_Cal_Para, Initialize, FindTmat, FindLmat, Find_Ray_Eqs, FindBeam_Eqs, Find_BeamPlanes, MapTo_Camera, Map_To_Cam_Tmat, MapTo_Cam_Geo, Make_Image_Pattern, Output_Pattern, Find_Cal_Face, OutputCubelmage. Some procedures are redundant and are written to ensure the correctness

8 of the program. 2.6.1. Read_Parameters This procedure reads the parameters of the projector, the camera, and the equation of the plane to be mapped. It also writes some parameters on the calibration file (unit 3) for use by the surface mapping programs. 2.6.2. Read_CalPara This procedure reads the calibration parameters, the size of the calibration cube, step in u, and v, and the rotation and the translation of the calibration cube. Any combination of the axes and angles can be entered. While reading the axis of the rotation and the angle of the rotation, this routine calculates the rotation matrix R. When it reads a character which is not an'x','y', or'z', it stops reading the rotation parameters and request translation parameters xt, yt, zt, along x, y, and z axis. It then computes the translation T and its inverse. The H matrix which accounts for both rotations and translations is found by multiplying the rotation matrix R by the translation matrix T. H-~ can be found by multiplying T' by R-'. R- is equal to Rt since the R matrix is orthogonal, and T' is equal to a translation matrix which translates - xt, - yt, - zt. 2.6.3. Initialize This procedure initializes the phantom program. 2.6.4. FindTmat This procedure computes the T matrix from the camera parameters given by the user. Te is the transformation matrix which transforms the world coordinates into the clip coordinates[61.

-Fsin(6) -Fyos(#)cos(O) -cos(#)cos(O) 0 F~cos(O) -Fysin(G)cos(O) -sin()sin(#) 0., 0 Fzsin(O) -cos(O) 0 FA-zssin(() -cos(4') +ysin(O)cos(e)) +ysin(G)sin(q) where Fx D / S, and Fy = D / S. Tp is the perspective transformation matrix. [v oo 1o, - 0 0 0 T= Vex Vq 0 1 0 000 where V,= - Nx / 2, V - Ny / 2, and Vcx-( Nx + 1 )/2, V ( Ny+ 1 ) /. Te matrix is computed by T T. Tp 2.6.5. FindLmat This procedure computes the L matrix from the parameters of the projector given by the user. All the equations are the same as for the T matrix except that the projector parameters are used instead of the camera parameters. 2.86.6 Find_Beam_Eqs This procedure computes the line equation of every beam that starts at the focal point of the projector and passes through the shutter point. The unit vectors on the shutter plane of the projector a^ and a, are a = a = -i sin() + j cos(8) a =~ -a = -i cos(O)cos(k) - j sin(O)cos(k) + k sin(X) where i, j, and k are the unit vectors of the Cartesian world coordinate system, and a#, and al

10 are the unit vectors of the spherical world coordinate system. The world coordinates of the shutter points r,, are ruv = rh + (u- ) - 0.51(2 ) + [(v — )- 0.5(2 )a, 2 ra+(u —.) 2 - N where r-h is the world coordinate of the center of the shutter plane of the projector. The direction numbers a, b, and c of the lines passing through the focal point ( xyfzf ) and the shutter point rv are defined as x-xf y-yf Z-Zf a b c where a = x coordinate of rv-Xf b = y coordinate of rv-yf c = z coordinate of rv,-Zf 2.6.7. Flnd_RayEqs This procedure computes the line equation of every ray which starts from the focal point of the camera and passes through the sample points of the camera plane. All the equations are same as Find_Beam_Eqs except that camera parameters used instead of the projector parameters. This routine is used to ensure the correctness of the T matrix. 2.8.8. FindBeamPlane This procedure computes for every column u the equation of the plane which is defined by the focal point (xyrz)f the oetoa of the projector and two points (u, ), (u, N) on thesame column u of the shutter plane of the projector. This is used to ensure the correctness of the L matrix. When al,bl,cl and a,b,,, are the direction numbers of the beams passing through the points (u, 1) and (u, N,) respectively, the equation of the plane is Ax + By + Cz + D = 0, where

II A - blc —bncl B - c0a-c0al C - alb,-anbl D - -xfA-yB-zfC 2.6.9. Map_o_GCamera This procedure computes the image of the surface point that is illuminated by the projector. The homogeneous coordinate of the surface point in the world coordinate system is multiplied by transformation matrix Te. (Xe YXcX 1) (x y z 1)Te where x~, Ye, and Zc are the clip coordinates, and x, y, z are the world coordinate of the surface point[31. Those points with -zc<xc<z and -yc<Xe<Ye are imaged on the camera plane. Finally image points are found by x, = round((-X )V,+V4) Ze,= round((ee )V,+Vy) ze where V,, V, Ve, and Vy are defined and computed as in the procedure Find_Tmat. 2.6.10. Map_to_Cam_Tmat This procedure performs the same calculations as MapTo_Camera. It is written to compare the image points obtained by this procedure with that obtained by MapToCamera. This procedure computes image points by multiplying the homogeneous coordinate of the surface point by the T matrix and then rounding it.

12 2.6.11. Mapto_Cam_Geo This procedure performs again the same computation as the two former ones, and is written to compare the image points obtained by this routine with those obtained by the former methods. The equation of the line which passes through the focal point of the camera (xf,yfrz) and the surface point (Xp,yp,Zp ) is x-Xf y-yf z-zf a b c where a - Xp-x, b = yp-Y c - Zp-Zf. The equation Ax + By + Cz + D = 0 of the camera plane is easily computed. A, B, C, and D are the direction numbers of the line which passes through the origin of the world coordinate and the center of the camera plane (xc,y~,z~). A = cos(O)sin(o) B = sin(6)sin(q) C = cos(,) D = -(Ax+Byc+Czc) Finally the intersection of the line and the plane is ( xr+ta,y+tb,zrttc ), where Ax^+Byr+Czr^D Aa+Bb+Cc The intersection point is rounded to the nearest sample point on the camera plane. 2.6.12. Make_ImagePattern This procedure computes the image points (xs, y,) that correspond to the shutter points (u, v). For each shutter point, the line equation of the beam that passes through the focal point and the shutter point is already computed by the procedure Find_BeamEqs. In Make_Image_Pattern the intersection of the beam and the plane to be mapped are computed, and are imaged on the camera plane by the procedures Map_to_Camera, Map_to_Cam_Tmat, and Map_to_Cam_Geo. The computed image points are then checked for equality.

13 If the two different surface points are mapped unto the same image point on the camera, the space coding methodology will fail to identify a proper column value for the image points and all the surface acqusition methodology will therefore fail. To avoid this situation, only the first shutter point which generated that image is recorded, and other shutter points are fagged on the terminal. Such a situation essentially signals a failure of the measurement. 2.6.13. Output Pattern This procedure creates a file (on unit 4) of all the patterns as dictated by the space coding method. Each entry consists of pairs of x~, y, Each pattern is terminated by the pair -1, -1. The number of patterns is a function of the shutter parameters. 2.6.14. Find_CalFace This procedure computes the new position of the eight vertices of the calibration cube and the new equations of the six faces of the calibration cube. At first, the calibration cube is located with its center at the origin of the world coordinate system. Denoting the length of the edge of the cube by 2a, the homogeneous coordinates of the eight vertices of the cube are V1 (a,a,a,l) V = (a,a,-al) V- (a,-a,-a,l) V4 = (a,-a,a,l) Vs - (-a,a,a,l) V = (-a,a,-a,l) V7 = (-a,-a,-a,l) V8 = (-a,-a,a,l) and the homogeneous equations of the six faces of the cube are F1 = (0,0,1,-a) F2 = (-1,0,0,-a) Fs - (0,0,-l,-a) F, - (1,0,0,-a) F6 - (0,-1,0,-a) F, = (0,1,0,-a) The new homogeneous coordinate V~ of the moved vertices are found by

14 V, - Vi-H where H is the transformation matrix computed at the procedure Read_Cal_Para The new homogeneous equations of the moved faces F,i are found by[5] F, — H-'F The number of the faces visible by both the projector and the camera are then computed. For the face F1 to be visible to both the projector and the camera, the focal points of both the projector and the camera should be at the same side of the face and should be at the opposite side of the face of the vertex V+i1. Since the sign of V-F tells at which side of the plane F lies the point V, the faces visible to both can be determined. 2.5.15. Output_Cube_Image This procedure finds the intersection of the beams of the projector and the faces of the calibration cube for each shutter point, where u and v are separated by "step". This routine outputs the world coordinate of the surface point, the column of the shutter which illuminated the surface point, and the image points on the camera of that surface point. 2.8.16. Main Program The main program opens file units 2, 3, and 4, and then calls the necessary procedures.

15 3. Surface Mapping Program This program calculates the surface point from the patterns of image points provided by either the phantom program or by an actual measurement. There are two versions of the surface mapping program, one which reads the T and L matrices provided by the phantom program, the other which calibrates the T and L matrices by the least square method. The surface mapping program reads the patterns of image points from a file (unit 4). From this data, the procedure computes the column of the shutter point which generated the image point, utilizing the space coding method[21. Using the T and L matrices, the program then solves a set of linear equations to find the surface point. The surface mapping program consists of a main program and the fol* lowing routines: Init, Calin, Imgin, Surout, Tmat, Lmat, Map, Column, and Setcof. 3.1. Init This routine calls the Calin routine, and Performs initialization. 3.2. Calin This routine reads the calibration data from file unit 3. The format of the file on unit 3 is line contents 1 number of sample points along u, v axes on the shutter plane 2 number of sample points along x,, y, axes on the camera plane 3 a point (x, y, z) on the face of the calibration cube which is illuminated by the projector 4 column u on the shutter plane, and image point (x,y,) lines 3 - 4 are repeated for each illuminated point

16 3.3. Imgin This routine reads the patterns of images on the camera plane from file unit 4. File unit 4 contains lines of x, and y, for each pattern. End of a pattern is indicated by a line of -1 -1. 3.4. Tmat This routine calculates the T matrix by the least square method utilizing the L2 routine. The routine calculates the normalized T matrix by setting T4 = 1 and dividing every element of the T matrix by To. The routine sets up the matrices B, C that are passed to the L2 routine to solve the equation Bx = C by the least square method, where x is the normalized T matrix, xt (TlT12,Ti,T 14,T21,T22,TT24,T4iT4T) and B, C are the least square matrices[7], x1 o -xs y' 0 -ys, z 0 -z, 1 0 o xl. -xl, 1 0 y' -yy, O z' -zyo 0 1 B= x 0y -x y - zY 0 - zy ZI - 1 0 o xi -x. 0 yo -yyl 0 z1 -z, 0 1 C= (-,-y1._ _ y ) where i = 1, 2,... N, N is the number of calibration data points, (x',y',z') are the world coordinates of a point on the face of the calibration cube, and (x,,yg) is the image point on the camera plane. When the surface mapping program does not calibrate the T matrix, this routine reads the T matrix calculated by the phantom program. 3.5. Lmat This routine finds the L matrix by the least square method. By setting L44= 1, and dividing every element of L matrix by L44, one finds the normalized L matrix. Following this,

17 the routine calculates the matrices B, C to be passed to the L2 routine to solve the equation Bx - C by the least square method, where x is the normalized L matrix, x- (L1,L14,L21,L^,L3,Ls,LL1) and B, C are the least square matrices, -x ul yl -y lul -_zlul ~B = x - Xslul y -y'u' z' -zI I C -I where ul is the column of the shutter. When the surface mapping program does not calibrate the L matrix, this routine reads the L matrix calculated by the phantom program. 3.6. Map This routine computes the surface points using the T and L matrices. It first finds the column of the shutter which generated the image by calling the routine Column. Next the routine sets up a linear equation by calling the routine Setcof, and then solves the linear equation by calling L2. 3.7. Column This routine uses the space coding method to infer the shutter column value for each image point. 3.8. Setcof This routine sets up the linear equation Ax = B where x is the world coordinate (x, y, z) of the surface point, and[lj

18 [Tu-Tl4x. T21-T2x, Ts1-Ts4x A = Ti-T1y, T=-T2y, T3-T3y, Lll-Ll4U Li2-L2u Lsr-L3u -T41+T44 B =-T4+T4Y. -L^i+LLu 3.9. Space Coding Space coding is used to find the column of the shutter which generated the image. In this coding method, the number of patterns necessary to encode and decode the column is log2( number of columns ) + 1. The space coding method is explained in detail elsewhere[l]. As an example, for the case of N. - 8, the 4 patterns of the space coding are: columns patterns 1 2 3 4 5 6 7 8 pattern 1 O O O O O O O O pattern 2 X X X O O O O pattern 3 X X O O X X O O pattern 4 X X O X 0 X where "O" indicates that the column is open, and "X" indicates that the column is closed.

19 4. The Analysis Program This program reads the file of surface points from the phantom program. This file is found on unit 2, and it contains, for all points, the shutter coordinate (u, v), the intersection (xl,yl,zl) of the beam which passes through the shutter point (u, v) and the surface plane, and the image point (x,, y,) which is the image of the intersection on the camera plane. The analysis program reads another file of computed surface points provided by the surface mapping program. This file consists of lines of a shutter column u, an image point (x,, y,); and an computed surface point (x2,y2z2). By comparing u, x,, and y,, the analysis program matches each surface point from the two files. It then calculates the distance between the two, the average distance, and the standard deviation of the distances. When the two files does not match, the program flags the surface points which do not match. This indicates something amiss in the phantom program or in the surface mapping program.

20 5. General Description of the Program Files There are four sets of files: Phan.geo and Eye.geo, Phan.tmat and Eye.tmat, Phan.tlmat and Eye.tlmat, and finally Phan.cal and Eye.cal. Programs beginning with "phan" are phantom programs, and programs beginning with "eye" are surface mapping programs. Phan.geo and Eye.geo do not use the T and L matrices. They use purely geometrical computation. Instead of T matrix, Eye.geo uses equations of rays which starts at the focal point of the camera and passes through the sample points of the camera plane. Instead of L matrix, Eye.geo uses equations of planes which pass through the focal point of the projector and two points on the same column of the shutter. Phan.tmat and Eye.tmat use the T matrix as computed from the parameters of the camera, and generate the same output as the.geo programs. This provides for an independent calculation which assures the correctness of the T matrix Phan.tlmat and Eye.tlmat use both the T matrix and the L matrix which is computed from the parameters of the projector. Again the same output as.geo and.tmat programs is generated, which ensuring the correctness of the L matrix. All the previous programs do not calibrate the T and L matrices. All the previous set of programs are basically the same, and Phan.tlmat and Eye.tlmat are chosen to represent the geometrical approach. Phan.cal and Eye.cal use T and L matrices, and calibrate the T and L matrices. The T and L matrices calibrated here should be close to the T and L matrices calculated above, which ensures the correctness of the calibration. The Support.2 file includes among other routines the routine L2 which solves the least square problem. The solution of the L2 routine was compared with the International Mathematical and Statistical Library routine LLBQF[IMSL] and was found accurate. The analysis program is found in the file Analysis.

21 8. Simulations With the above set of programs, it is possible to simulate various situations and analyze the output obtained, e.g. the effect of the precision of the phantom output on the accuracy of the surface mapping program can be tested. The precision of the phantom output is determined by the number of the sample points on the camera plane, the number of the sample points on the shutter plane of the projector, and the precision of the intersection points of the beam of the projector with the calibration cube. Since phan.tlmat and eye.tlmat calculate the T and L matrices by geometrical formulas, they represent the geometrical approach and the data from them is the bound of the accuracy which can be attained by the surface mapping program. Table 1 illustrates the effect of the positions of the shutter and the camera on the accuracy of the stereo camera system. With the shutter fixed, various positions of the camera has been tried. The numbers are the average distance between the actual surface point and the computed surface point, and its standard deviation. Table 2 illustrates the effect of the number of sample points on the accuracy of the stereo camera system. Here, the column headed by "geometry" is the lower bound which can be obtained by the stereo camera system. In Table 3, the number of decimal places of the surface point on the face of the calibration cube has been changed. In most situations, the calibration of the T and L matrices was satisfactory. However, sometimes calibration does not provide accurate T and L matrices even though the condition number from the L2 routine is not high. Table 4 illustrates the effect of the position of the calibration cube on the accuracy of the stereo camera system. This set of programs can be used to find the satisfactory position of the calibration cube: First, find the average distance by the set of programs which use geometry. This is the lower bound. Next, try various positions of the calibration cube until you obtain the average distance close to the lower bound.

22 camera position standard average 0 deviation 45 60 1.26 1.33 60 60.63.60 90 60.32.22 135 60.29.12 210 60 1.40 2.03 135 5 6.91 21.1 135 15.91.73 135 30.73.17 135 60.29.12 135 75.27.12 135 85.25.12 Table 1. Effect of Camera Position shutter: r - 20, 0 - 30, -= 45, D - 4, Su, Sv = 2 N,,N, = 16 camera: r = 20, D - 4, Sx, Sy = 2 N,,Ny = 32 equation of the plane: z - 1

23 geometry calibration standard standard NNy average average deviation deviation 32 1.26 1.33 2.42 2.92 N,,Nv 8 128.29.31.49.53 32 1.23 1.42 2.04 2.62 N,NN, 16 128.31.32.35.37 32 1.44 1.26 2.58 4.78,,Nu, 32 128.30.31.45.55 Table 2. Effect of Number of Sampling Points shutter: r - 20, 0 - 30, ~ = 45, D - 4, Su, Sv = 2 camera: r - 20, 0 = 60, 4 = 60, D = 4, Sx, Sy = 2 calibration: size = 8, rotation = x 10, y 45, z 30, translation = 0, 0, 0 equation of the plane: z = - 1

24 number of standard average decimal places deviation 7 1.08 1.37 3 1.08 1.37 2 1.06 1.35 1 1.06 1.32 0 4.80 6.23 Table 3. Effect of Number of Decimal Places shutter: r = 20, 8 = 30, H = 45, D - 4, Su, Sv = 2 Nu,N, = 16 camera: r -20, 8 = 60, - 60, D -4, Sx, Sy = 2 NN, Ny 32 calibration: size = 8, rotation = x 10, y 45, z 30, translation = 0, 0, 0 equation of the plane: z = - 1

25 rotations standard average x y z deviation 30 30 60 10.7 16.0 10 45 30 1.07.98 N,,N, -8 geometry.70.52 30 30 60.99.71 10 45 30.95.70 N,,N, = 16 geometry.77.49 Table 4. Effect of Rotation of Calibration Cube shutter: r - 20, 0 = 30, b = 45, D - 4, Su, Sv = 2 camera: r - 20, 0 = 60, = - 60, D - 4, Sx, Sy = 2 calibration: size = 8, translation - 0, 0, 0 equation of the plane: z = - 1 7. References Ill M. D. Altschuler, J. L. Posdamer, G. Frieder, "The Numerical Stereo Camera," SPIE, vol. 283, 1982 12) J. L. Posdamer, M. D. Altschuler, "Surface Measurement by Space-encoded Projected Beam System," Computer Graphics and Image Processing, vol. 18 pp 1-17, 1982 [3] W. M. Newman, Principles of Interactive Computer Graphics, 2nd ed. pp 339-344, McGraw-Hill, 1979 14) G. E. Forsythe, M. Malcolm, C. B. Moler, Computer Methods for Mathematical Computations, pp 192-235, Prentice-Hall, 1977 [5] W. M. Newman, et. al., op. cit. pp 492-501 16) W. M. Newman, et. al., op. cit., pp348-351 71 D. V. Roger, J. A. Adams, Mathematical Elements for Computer Graphics, pp 78-82, McGraw-Hill, 1976

26 8. Program Listings 8.1. Phantom Program

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:9 PAGE LINE NUMBER B P C I STMT # S 0 U R C E P R GR A M V —-+ 1 + 2 + — -- 3- + -- 4 --— + — 5 --— + -— 6- - +-.7. +-8 9..+ —- -IV 1.000 PROGRAM PhanCal ( INPUT, OUTPUT ); 2.000 3.000 / 4.000 5.000 6.000 7.000 8.000 *Phantom Program 9.000 * for Testing Numerical Stereo Camera 10.000 11.000 * Programmed by Myoung Lee 12.000 13.000 * July 1983 14.000 15.000 * 16.000 * The Phantom Program provides data for the Surface Mapping 17.000 * Program. This program can be used for the testing and 18.000 * experiment of Surface Mapping Program. 19.000 * This also generates actual surface points to be mapped 20.000 * Into camera points. 21.000 * This calculates T and L matrix by geometrical consideration. 22.000 23.000 * For a detailed description of the phantom program and the 24.000 * related discussion, refer to the separate documentation. 25.000 26.000 / 27.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 2 LINE NUMBER B P C I STMT N SOUR C E PROGRAM V -— + —— I —— + —-— 2 —— + —— 3 —----— 4 —-----— 5 —----— 6 —----— 7 —----— 8 28.000 29.000 CONST 30.000 31.000 u_max = 64: maximum number of sample points along u axis ) 32.000 v max = 64; maximum number of sample points along v axis) 33.000 xs_max = 128; 4 maximum number of smaple points along Xs axis) 34.000 ysmax a 128 4 maximum number of sample points along Ys axis) 35.000 36.000 unobserved = -999; 37.000 pre - 7; number of decimal places of real number ) 38.000 wid 14; 4 length of the real field ) 39.000 widf wid - 1; 4 for FORTRAN interface ) 40.000 len = 6; length of the Integer field ) 41.000 lenf = len - 1 4 for FORTRAN interface ) 42.000 43.000 44.000 TYPE 45.000 46.000 Vector a ARRAY (. 1..3.) of REAL 47.000 Mat_4 = ARRAY (. 1..4, 1..4.) of REAL; 48.000 49.000 VAR 50.000 51.000 Calf, 4 Calibration file ) 52.000 lmg_f. 4 Image file ) 53.000 Ana_f: TEXT;( File for Anlysis program ) 54.000 55.000 Tmat, ( T matrix ) 56.000 Lmat, 4 L matrix ) 57.000 Viewxform: ARRAY (. 1..4, 1..4.) of REAL; 58.000 4 Transforms world coordinates into clip coordinates ) 59.000 60.000 Ray: ARRAY (. 1..xs max, 1..ys-max, 1..3.) of REAL; 61.000 4 direction numbers of each ray passing through the focal point ) 62.000 4 of the camera and sample points on the camera plane ) 63.000 64.000 Beam: ARRAY (. 1..umax, i..vr max, 1..3.) of REAL 65.000 4 direction numbers of each beam passing through the focal point ) 66.000 4 of the shutter and sample points on the shutter plane ) 67.000 68.000 Beam_plane: ARRAY (. 1..umax, 1..4.) of REAL 69.000 4 equation of the planes passing through the focal point of the shutter ) 70.000 4 and two sample points on the same column of the shutter plane i 71.000 Screen: ARRAY (. 1..xs_max, 1..ysmax.) of BOOLEAN; 72.000 4 If two different shutter points Imaged on the same camera point ) 73.000 4 Space Coding gives a spurious column i 74.000 To avoid this, the camera point is recorded true when it is imaged ) 75.000 76.000 Pattern: ARRAY (. I..umax...vmax, 1..2.) of INTEGER; 77.000 4 image point for each sutter point ) 78.000 79.000 H mat, ( calibration cube transformation matrix } a f fNr\ U 4- A...-~x A

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19 PAGE 3 LINE NUMBER B P C I STMT # S 0 U R C E P R 0 G R A M V —-+ --- -— + —-- 2 —— + —-- 3 —----- 4 —----- + 5 + —----- 6 —-- + —-- 7 — + - + —---— IV 81.000 82.000 83.000 { Every name ending with "sh" is associated with shutter ) 84.000 "u" and "v" are coordinate system in the shutter plane ) 85.000 86.000 ( Every name ending with "ca" is associated with camera ) 87.000 ( "x" and "y" are coordinate system In the camera plane ) 88.000 89.000 ( shutter parameters ) 90.000 91.000 92.000 Thetash, 93.000 Phi sh. 94.000 Distsh: REAL; { position of the center of the shutter ) 95.000 96.000 Center sh, { Coordinate of the center of the shutter ) 97.000 FP sh, { coordinate of the Focal Point of the shutter) 98.000 Dir no sh:Vector; { direction number of the center of the shutter ) 99.000 100.000 D_sh, ( distance from center of the shutter to focal point of the shutter ) 101.000 Su, ( width of the shutter plane ) 102.000 Sv: REAL; 4 height of the shutter plane ) 103.000 104.000 Nu, { Number of sample points in u direction ) 105.000 Nv: INTEGER; ( Number of sample points In v direction ) 106.000 107.000 108.000 { camera parameters ) 109.000 110.000 111.000 Theta ca, 112.000 Phi ca, 113.000 01st_ca: REAL; 4 position of the center of the camera } 114.000 115.000 Centerca, 116.000 FP_ca, 117.000 Dir_noca: Vector 118.000 names are ) 119.000 Dca, ( same as ) 120.000 Sx, ( above } 121.000 Sy: REAL; ( with camera ) 122.000 ( Instead of ) 123.000 Nx, ( shutter ) 124.000 Ny: INTEGER 125.000 126.000 Map_plane: ARRAY (. 1..4.) of REAL; ( Plane to be mapped ) 127.000 128.000 129.000 { calibration cube ) 130.000 131.000 Face; ARRAY (. 1..3, 1..4.) of REAL; 132.000 ( equation of the visible faces of the calibration cube) 133.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 4 LINE NUMBER B P C I STMT # S 0 U R C E P R O G R AM V —-+ + —--— + ----- 4 2 —-- + —. 3 —----— 4 —------ 5 —------ -6 —----— 7 —------ 1V 134.000 Num_face: INTEGER; { number of faces visible by shutter and camera ) 135.000 136.000 a: REAL; { length of the side of a cube = 2 * a ) 137.000 138.000 xt, yt, zt: REAL; { translation of the cube ) 139.000 140.000 Step: INTEGER; { in calibration ) 141.000

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19 PAGE 5 LINE NUMBER B P C I STMT SOUR C E PRO R AM V —-7 —— I —-----— 2 —— 9 —-— 3 —-- -t —142.000 143.000 144.000 / 145.000 * Make_ldentity 146.000 * 147.000 * This makes a 4 by 4 matrix an identity matrix 148.000 149.000 150.000 1 PROCEDURE Make_ldentity ( VAR Matrix: mat_4 151.000 1 152.000 1 VAR 153.000 1. j: INTEGER; ( indices ) 154.000 1 155.000 1 BEGIN 156.000 1 157.000 1 Initialize 4 by 4 Matrix as Identity matrix ) 158.000 1 159.000 1 1 1 FOR I:- I TO 4 00DO 160.000 1 2 2 FOR j:- 1 TO 4 00DO 161.000 3 Matrix (.1. J.):= 0 162.000 1 163.000 1 1 4 FOR 1:= 1 TO 4 DO 164.000 1 5 Matrix (.1, 1.): I 165.000 1 166.000. END: 167.000 168.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 6 LINE NUMBER B P C I STMT i SOUR C E PROGRAM v —+ —— ^ —-------— 2 —+ —-— 3 —------— 4 —-----— 1V 169.000 170.000 171.000 1/ 172.000 MultMat4 173.000 174.000 175.000 * This multiplies two 4 by 4 matrix 176.000 */ 177.000 178.000 1 PROCEDURE MultMat4 ( VAR Mat m Mat_4; 179.000 1 Mat_a, Mat_b: mat_4 ); 180.000 1 181.000 i VAR 182.000 1 1, j, k: INTEGER { indices ) 183.000 1 184.000 1 BEGIN 185.000 1 186.000 1 ( multiply 4 by 4 Matrix Mat_a and Matb ) 187.000 1 188.000 1 1 1 FOR 1:= 1 TO 4 DO 189.000 1 2 2 FOR j:= 1 TO 4 DO 190.000 1 2 BEGIN 191.000 1 1 2 3 Mat (.1, j.):0 192.000 1 1 3 4 FOR k: i TO 4 DO 193.000 1 1 2 5 Mat (.1, j.):- Mat (.1. j.) + Mat_a (.1, k.) * Mat_b (.k, j.); 194.000 1 2 END 195.000 END; 196.000

pASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19 PAGE 7 LINE NUMBER B P C I STMT N S 0 U R C E P R O G R AM V —+ —— I —— + —— 32 —---— 3 —----— 4 — ----------- 6-7 197.000 198.000 199.000 * 200.000 * ReadParameters 201.000 * 202.000 * This reads parameters for phantom program 203.000 / 204.000 205.000 1 PROCEDURE ReadParameters; 206.000 1 207.000 208.000 1 BEGIN 209.000 I 210.000 1 1 WRITELN; 211.000 1 2 WRITELN (' Enter Parameters' ); 212.000 1 3 WRi TELN; 213.000 1 214.000 I 4 WRITELN (' —-- shutter parameters ----' ); 215.000 1 5 WRITELN; 216.000 1 217.000 1 6 WRITELN (' focal point -- r' ) 218.000 1 7 READLN ( 0stsh; 219.000 1 220.000 1 8 WRITELN (' focal point -- thata' ) 221.000 1 9 READLN ( Theta_sh; 222.000 1 223.000 1 10 WRITELN (' focal point — phi'; 224.000 1 11 READLN Phish; 225. 000 1 226.000 112 WRITELN (' shutter plane -- 0' ) 227.000 113 READLN ( Dsh ); 228.000 229.000 1 14 WRITELN ( * shutter plane -- Su; 230.000 1 15 READLN ( Su ) 231.000 1 232.000 1 16 WRITELN (' shutter plane -- Sv' ) 233.000 1 17 READLN ( Sv ); 234.000 235.000 1 li WRITELN (' shutter plane -- Nu' 236.000 19 READLN ( Nu ); 237.000 1 20 WRITE ( calf, Nu:lenf 238.000 1 239.000 1 21 WRITELN (' shutter plane -- Nv' 240.000 1 22 READLN ( Nv ); 241.000 1 23 WRITELN ( calf, Nv:len 242.000 1 243.000 1 24 WRITELN i' ---- camera parameters --' 244.000. 25 WRITELN 245.000 1 246.000 1 26 WRITELN (' focal point -- r' ); 247.000 1 27 READLN ( Distca ) 248.000 249.000 1 28 WRITELN (' focal point -- theta' );

PASCAL/VS RELEASE 2.1 PHAN_CAL:READ_PARAMETERS 10/10/83 11:28:19 PAGE 8 LINE NUMBER B P C I STMT S 0 U R C E P R OG R A M V —-+ — 1 — + 2 -— + -— 3 -— + --— 4 -— + --— 5 + --- 6..+..7..+.-8..+ -9.+ —- 1V 250.000 1 29 READLN ( Theta_ca ); 251.000 1 252.000 1 30 WRITELN (' focal point -- phi' ); 253.000 1 31 READLN ( Ph_ca ); 254.000 1 255.000 1 32 WRITELN (' camera plane -- D'); 256.000 1 33 READLN ( Dca ); 257.000 1 258.006 1 34 WRITELN (' camera plane -- Sx' ); 259.000 1 35 READLN ( Sx ); 260.000 t 261.000 1 36 WRITELN (' camera plane -- Sy' ); 262.000 1 37 READLN ( Sy ): 263.000 1 264.000 1 38 WRITELN (' camera plane -- Nx' ); 265.000 1 39 READLN ( Nx ); 266.000 1 40 WRITE ( cal f, Nx:lenf ); 267.000 1 268.000 1 41 WRITELN (' camera plane -- Ny' ); 269.000 1 42 READLN ( Ny ); 270.000 1 43 WRITELN ( calf, Ny:len ); 271.000 1 272.000 1 44 WRITELN (' Equatton of the plane to be mapped -- A, B, C, D' ); 273.000 1 45 READLN ( Map plane(.l.), Mapplane(.2.), Map_plane(.3.), Map_plane(.4.) ); 274.000 1 46 WRITELN; 275.000 1 276.000 END

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19PAGE 9 LINE NUMBER B P C I STMT N S O. U R C E P R 0 G R A M V —------- I ------— 2 -------— 3 —-— 1 —-- 277.000 278.000 279.000/* 280.000 * Intlalize 281.000 * 282.000 * 283.000 / 284.000 285.000 1 PROCEDURE Initlalize; 286.000 1 287.000 1 VAR 288.000 1 ang_rad: REAL; ( angle to radian conversion ) 289.000 1 1, j: INTEGER; 4 indices ) 290.000 291.000 1 BEGIN 292.000 1 293.000 1 convert angles to radian ) 294.000 1 295.000 1 1 ang_rad:~ 3.141592 / 180; 296.000 1 297.000 12 Theta_sh:= Theta_sh * ang_rad; 298.000 1 3 Phi_sh:a Phish * ang_rad; 299.000 1 300.000 1 4 Thetaca:- Theta_ca * ang_rad; 301.000 1 5 Phl_ca:- Phi ca * ang_rad; 302.000 1 303.000 1 calculate direction numbers ) 304.000 1 305.000 1 6 DIr_no_sh (.1.): COS ( Theta_sh )* SIN ( phi_sh; 306.000 1 7 Dirnosh (.2.):~ SIN ( Theta_sh )* SIN ( phish 307.000 1 8 Dir_no_sh (.3.):w COS phi_sh ); 308.000 1 309.000 1 9 Dir_no_ca (.1.):~ COS ( Thetaca )* SIN ( phica ) 310.000 1 10 Dir_no_ca (.2.):~ SIN ( Thetaca )* SIN ( phica ) 311.000 1 11 Dir_no_ca (.3.):u COS ( phlca ); 312.000 1 313.000 I calcuate Cartesian coordlate of focal points and center of shutter and camera ) 314.000 315.000 1 1 12 FOR I:1I TO 3 DO 316.000 1 1 BEGIN 317.000 1 1 1 13 FPsh (:1.):= Distsh * Dirnosh (.1.) 318.000 1 1 1 14 FPca (.1.): Dist_ca * Dirno_ca (.1.); 319.000 1 1 1 320.000 1 1 1 15 Centersh (.1.): ( Dlst_sh - D_sh ) * Dir_nosh (.1.); 321.000 1 1 I 16 Ceniteica (.1.):~ ( Distca - Dca ) * Dirnoca (.1.) 322.000 1 END 323.000 1 324.000 1 I nitialize screen ) 325.000 1 1 17 FOR I:= 1 TO Ny DO 326.000 1 2 18 FOR j:= 1 TO Nx D00'327.000 1 2 19 Screen (. 1. j.): FALSE 328.000 I.2 329.000 1 2

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 10 LINE NUMBER B P C I STMT # S 0 U R C E P R 0 G R A M V —— + - --- +-+ 2 -+ —3 - -+ — — 4 + --— 5 — — + -— 6 + 7 + 8 + 9 330.000 END 331.000 332.000

PASCAL/VS RELEASE 2.1'PHAN_CAL 10/10/83 11:28:19 PAGE 11 LINE NUMBER B P C I STMT S O U R C E P R O G R AM V —-+ + --- 2 + --— + —-3 —— + 4 + ----— 5 —------ 6 —— + --— 7 —— + 8 + 9 -— 1V 333.000 334.000 335.000 /* 336.000 * Find Tmat 337.000 * 338.000 * This calculates T matrix by geometrical consideration 339.000 340.000 341.000 1 PROCEDURE Find Tmat 342.000 1 343.000 1 VAR 344.000 1 Pers_xform: ARRAY (. 1..4, 1..4.) of REAL; 345.000 1 { Perspective transformation matrix ) 346.000 1 SFx, SFy: REAL; ( Scale Factors') 347.000 1 x, y, z: REAL; ( focal point of the camera ) 348.000 1 1, J k: INTEGER; { Indices ) 349.000 1 temp: REAL; { temporary ) 350 000 1 351.000 1 BEGIN 352.000 1 353.000 1 1 SFx:- D_ca / Sx; 354.000 1 2.SFy:= Dca / Sy; 355.000 1 356.000 1 3 x: FP ca (.1.) 357.000 1 4 y:" FP_ca (.2.); 358.000 5 z:= FP ca (.3.) 359.000 1 360.000 1 6 Viewxform (. I, 1.): - SFx * SIN(Theta ca) 361.000 1 7 Vlewxform (. 2, 1.):- SFx * COS(Thetaca) 362.000 1 8 1Vew_xform (. 3, 1.): 0; 363.000 9 Viewxform (. 4, 1.): SFx * ( x * SIN ( Theta_ca ) - y * COS ( Theta_ca ) ); 364.000 1 365.000 1 10 Vlew_xform (. 1, 2.): - SFy * COS ( Theta_ca ) * COS ( Phlca ); 366.000 1 II Vlew_xform (. 2, 2.) - SFy SIN ( Theta ca ) * COS ( Phlca ); 367.000 1 12 Vlew_xform (. 3, 2.):2 SFy * SIN ( Phi ca; 368.000 1 13 View_xform (. 4, 2.): SFy * ( - z SIN(Phl_ca) + x * COS(Theta_ca) * COS(Phl_ca) 369.000 1 + y * SIN(Theta_ca) * COS(Phi_ca) ); 370.000 1 371.000 1 14 View_xform (. 1, 3.): -COS ( Theta_ca ) * SIN ( Phica ); 372.000 1 15 Viewxform (. 2, 3.):- SIN ( Theta_ca ) * SIN ( Ph_ca ); 373.000 1 16 Vew_xform (. 3, 3.):- COS ( Phi ca ) 374.000 1 17 Vlew_xform (. 4. 3.): z * COS(Phl ca) + x * COS(Thetaca) * SIN(Phl_ca) 375.000 1 + y * SIN(Theta_ca) * SIN(Ph_ca); 376.000 I 377.000 1 18 View_xform (., 4.):0; 378.000 1 19 Viewxform (. 2, 4.):0; 379.000 1 20 Viewxform (. 3, 4.):; 380.000 1 21 Viewxform (. 4, 4.):; 381.000 1 382.000 O 383.000 1 ( Perspective transform ) 384.000 1 385.000 I 1 22 FOR 1:= I TO 4 DO

PASCAL/VS RELEASE 2.1 PHAN CAL:FIND TMAT 10/10/83 11:28:19 PAGE 12 LINE NUMBER B P C I STMT SOUR C E PROGRAM V --------- --------- 2 —------ 3 —------------- 386.000 i 2 23 FOR J:- 1 TO 4 DO 387.000 1 24 Peraxform (. i, j.): 0 388.000 389.000 25 Pers-xform (. 1, 1.): Nx / 2 390.000 1 26 Pers_xform (. 2, 2.): Ny/ 2 391.000 1 27 Pers_xform (. 3, 1.): ( i + Nx ) 2; 392.000 1 28 Pers_xform (. 3. 2.): ( i + Ny ) 2; 393.000 1 29 Pers_xform (. 3, 4.):I; 394.000 1 395.000 1 396.000 1 multiply Viewxform by Persxforn ) 397.000 1 398.000 1 30 WRITELN; WRITELN 399.000 1 32 WRITELN (' --— T matrix ----' ); 400.000 1 1 33 FOR 1:~ 1 TO 4 00 401.000 1 1 BEGIN $02.000 1 1 403.000 1 1 2 34 FOR J:a 1 TO 4 D00 404.000 1 2 BEGIN 405.000 2 1 2 35 temp:; 406.000 2 1 3 36 FOR K:I TO 4 DO 407.000 2 1 2 37 temp:= temp + Viewxform (. 1, k ) * Pers_xform (. k, j.); 408.000 2 1 2 38 Tmat (. I, j.):= temp; 409.000 2 1 2 39 WRITE ( Tmat (.1, j.):w1d:pre ) 410.000 1 1 1 END; (for j) 411.000 1 1 1 40 WRITELN 412.000 1 I 1 413.000 1 END ( for 1) 414.000 i 415.000 1 41 WRITELN 416.000 1 1 42 FOR I: 1 TO 4 DO 417.000 1 1 BEGIN 418.000 1 1 2 43 FOR j:- 1 TO 4 DO 419.000 i i 1 44 WRITE ( Tmat (.1, j.) / Tmat (.4, 4.):w1d:pre ); 420.000 1 I WRITELN 421.000 1 45 END 422.000 1 1 423.000 END; 424.000 425.000

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19 PAGE 13 LINE NUMBER B P C I STMT # SOU R C E P R O G R AM V —-+.t.- -+ 2 + -— 3 -----— 4 --— + - -- + -— 6 ---- --— 7 + + -------------- 426.000 427.000 428.000 /* 429.000 * FindLmat 430.000 * 431.000 * This calculates L matrix from the geometrical consideration 432.000 */ 433.000 434.000 1 PROCEDURE FindLmat 435.000 1 436.000 1 VAR 437.000 1 View mat: ARRAY (. 1..4, 1..4.) of REAL; 438.000 1 ( L matrix equivalent of View_xform of T matrix ) 439.000 1 Pers mat: ARRAY (. 1..4, 1..4.) of REAL; 440.000 1 { Perspective transformation matrix of shutter ) 441.000 1 SFx, SFy: REAL; ( Scale Factors ) 442.000 1 x, y, z: REAL; ( focal point ) 443.000 1, J, k: INTEGER; { Indices ) 444.000 1 temp: REAL; I temporary ) 445.000 t 446.000 t ~ BEGIN 447.000 1 448.000 1 1 SFx:' Dsh / Su; 449.000 1 2 SFy: D_sh / Sv; 450.000 1 451.000 1 3 x: FP sh (.1.); 452.000 1 4 y: FPsh (.2.); 453.000 1 5 z: FP sh (.3.); 454.000 1 455.000 1 6 View mat (. 1, I.): - SFx * SIN(Thetash); 456.000 1 7 Viewmat (. 2, 1.):- SFx * COS(Theta sh); 457.000 1 8 View mat (. 3,1.):0; 458.000 1 9 View_mat (. 4, 1.):= SFx * ( x * SIN ( Theta_sh ) - y * COS ( Theta_sh ); 459.000 1 460.000 1 10 View_mat (. 1, 2.):- SFy * COS ( Theta_sh ) * COS ( Phl_sh ); 461.000 1 11 View_mat (. 2, 2.): - SFy * SIN ( Thetash ) * COS ( Phl_sh ); 462.000 1 12 Viewmat (. 3, 2.): SFy * SIN ( Phi sh ); 463.000 1 13 Vlew_mat (. 4, 2.):~ SFy * ( - z * SIN(Phish) + x + COS(Theta_sh) * COS(Phi_sh) 464.000 1 + y * SIN(Theta_sh) * COS(Phl_sh) ); 465.000 1 466.000 1 14 Viewmat (. 1, 3.):- COS ( Theta sh ) * SIN ( Phi sh ); 467.000 1 15 View_mat (. 2, 3.): - SIN ( Thetash ) * SIN ( Phi_sh );.468.000 1 16 View_mat (. 3, 3.):- COS ( Phi_sh ); 469.000 1 17 Viewmat (.4, 3.): z * COS(Phi sh) + x * COS(Theta sh) * SIN(Phi_sh) 470.000 1 + y * SIN(Theta_sh) * SIN(Phi_sh); 471.000 1 472.000 1 18 Viewmat (. 1, 4.): 0; 473.000 19 View mat (. 2. 4.):0; 474.000 1 20 View mat (. 3, 4.):= 0 475.000 1 21 View mat (. 4, 4.) 1; 476.000 1 477.000 1 478.000 1

PASCAL/VS RELEASE 2.1 PHAN_CAL:FIND_LMAT 10/10/83 11:28:19 PAGE 14 LINE NUMBER B P C I STMT N SO U R C E PROGORAM V ------- I —-----— 2 —--— + —^-3 —— + —-— 4 —-----— 5 —— + —-— 6 —----— 7 —— + —-— 8 —----— 9 — 479.000 1 ( Perspective transform ) 480.000 1 481.000 1 1 22 FOR 1:= I TO 4 DO 482.000 2 23 FOR j:= i TO 4 00DO 483.000 1 24 Pers_mat (. 1, j.): 0; 484.000 1 485.000 i 25 Persmat (.1, 1.): Nu 2 486.000 1 26 Pers_mat (. 2, 2.): Nv /2 487.000 1 27 Pers_mat (. 3. i.): ( 1 + Nu )/2; 488.000 1 28 Persmat (. 3, 2.): I ( i + Nv ) /2; 489.000 1 29 Persmat (. 3, 4.):1; 490.000 1 491.000 492.000 1 multiply Viewmat by Pero_mat I 493.000 494.000 1 30 WRITELN; WRITELN; 495.000 1 32 WRITELN (' ---- L matrix ----' ); 496.000 1 1 33 FOR 1:~ 1 TO 4 00DO 497.000 1 1 BEGIN 498.000 1 1 499.000 1 1 2 34 FOR j:= I TO 4 DO 500.000 1 1 2 BEGIN 501.000 2 1 2 35 temp: 0; 502.000 2 1 3 36 FOR K: TO 4 D00 503.000 2 1 2 37 temp:- temp + View mat 4. 1, k.) * Pers_mat (. k, j.); 504.000 2 1 2 38 Lmat (. 1, j.):u temp; 505.000 2 1 2 39 WRITE( Lmat (.1, J.):w1ld:pre 506.000 1 1 1 END (; for j) 507.000 1 1 508.000 1 1 40 WRITELN 509.000 1 1 510.000 1 END;(for 1) 511.000 1 512.000 1 513.000 1 41 WRITELN 514.000 1 1 42 FOR 1:= 1 TO 4 DO 515.000 1 1 BEGIN 516.000 1 1 2 43 FOR j:= i TO 4 DO 517.000 1 1 44 WRITE ( Lmat (.1, J.) / Lmat (.4, 4.):w1d:pre; 518.000 1 i1 WRITELN 519.000 1 1 45 END 520.000 1 1 521.000 1 1 522.000 END 523.000

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/3 11:28:9 PAGE 15 LINE NUMBER B P C I STMT # S 0 U R C E P R O G R AM V —-+ -— + —— 2 —— + —— 3 —-- ---- 4 —----— + 5 --— + -6 —----— + 7 --— 8 —----— + —- 1V 524.000 525.000 526.000 / 527.000 FindRay_Eqs 528.000 529.000 * This calculates equations of the lines from the focal point 530.000 * of the camera to the sample points on the camera plane. 531.000 * This is used for verifing the T matrix calculated by 532.000 * geometrical consideration. 533.000 / 534.000 535.000 1 PROCEDURE FindRay_Eqs; 536.000 1 537.000 1 VAR 538.000 1 axe, 539.000 1 aye: Vector; { unit vectors In the camera plane ) 540.000 1 541.000 1 Camera: ARRAY (. 1..xsmax,..ys_max, 1..3.) of REAL; ( coordinates of camera points ) 542.000 1 543.000 1 xs, ys, k: INTEGER; ( indices } 544.000 1 545.000 1 BEGIN 546.000 1 547.000 1 ( find unit vectors In the camera plane ) 548.000 1 549.000 1 axe (.1.):- - SIN ( Theta_ca ); 550.000 1 2 axe (.2.):- COS ( Theta_ca ); 551.000 1 3 axe (.3.):- 0; 552.000 1 553.000 1 4 aye (.1.):= - COS ( Phi_ca ) * COS ( Theta_ca ); 554.000 1 5 aye (.2.):= - COS ( Phl_ca ) * SIN ( Theta_ca ); 555.000 1 6 aye (.3.):~ SIN ( Phi_ca ); 556.000 1 557.000' 558.000 1 I find the line eq of each ray ) 559.000 1 560.000 1 561.000 1 1 7 FOR xs:= 1 TO Nx DO 562.000 1 2 8 FOR ys: I TO Ny DO 563.000 1 2 BEGIN 564.000 1 2 565.000 1 1 3 9 FOR k: I TO 3 DO 566.000 1 1 3 10 Camera (.xs,ysk.):- Center_ca (.k.) + ( 2 * xs - Nx - i ) * (Sx / Nx) * axe(.k.) 567.000 1 1 2 + ( 2 * ys - Ny - 1 ) * (Sy / Ny) * aye(.k.); 568.000 1 I 2 569.000 1 1 3 I1 FOR k:= 1 TO 3 DO 570.000 1 1 3 BEGIN 571.000 2 1 3 12 Ray (.xs, ys. k.):= Camera (.xs, ys, k.) - FP_ca (.k.); 572.000 1 2 END; 573.000 t 1 2 574.000 1 2 END 575.000 1 2 576.000 END; ( find_ray_eq )

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 16 LINE NUMBER B P C I STMT N SO U R C E PROGRAM V —+ —---- -— + —— 2 —— + —— 3 —— + —— 4 —----— 5 —----— 6 —----— 7 —----— 8 577.000 578.000 579.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 17 LINE NUMBER B P C I STMT # SOU R C E P R O G R AM V —-+ —-- -— + —-— a2 —-----— 3 —— + —— 4 —----— 5 —----— 6 —-----— 7 —------ 8-V 580.000 581.000 582.000 /* 583.000 * Find_Beam_Eqs 584.000 585.000 586.000 * This calculates equations of the lines from the focal point 587.000 * of the shutter to the sample points on the shutter plane. 588.000 * This is used for finding the Intersection of each beam 589.000 * with the plane to be mapped. 590.000 / 591.000 592.000 1 PROCEDURE Find_BeamEqs; 593.000 1 594.000 1 VAR 595.000 1 axe. 596.000 1 aye: Vector; { unit vectors in the shutter plane ) 597.000 1 598.000 t Shutter: ARRAY (. 1..u_max, 1..v_ max, 1..3.) of REAL: { coordinates of shutter points ) 599.000 1 600.000 1 u, v, k: INTEGER; i ndices ) 601.000 1 602.000 1 BEGIN 603.000 1 604.000 1 ( find unit vectors in the shutter plane ) 605.000 1 606.000 1 axe (.1.):= - SIN ( Thetash ); 607.000 1 2 axe (.2.):- COS ( Theta_sh ); 608.000 1 3 axe (.3.):= 0 609.000 1 610.000 1 4 aye (.1.):= - COS ( Phi_sh ) * COS I Theta_sh ); 611.000 1 5 aye (.2.):= - COS ( Phi_sh ) * SIN ( Theta_sh ); 612.000 1 6 aye (.3.): SIN ( Phl_sh ); 613.000 1 614.000 1 615.000 1 { find the line eq of each beam ) 616.000 1 617.000 1 618.000 1 1 7 FOR v: 1 TO Nv DO 619.000 1 2 8 FOR u:= 1 TO Nu DO 620.000 1 2 BEGIN 621.000 1 2 622.000 I 3 9 FOR k:- I TO 3 DO 623.000 1 I 3 10 Shutter (.u,v,k.):= Center_sh (.k.) + ( 2 * u - Nu - 1 ) * (Su / Nu) * axe(.k.) 624.000 1 2 + ( 2 v - Nv - ) * (Sv / Nv) aye(.k.); 625.000 1 1 2 626.000 1 1 3 11 FOR k:= 1 TO 3 DO 627.000 1 1 3 12 Beam (.u, v, k.):~ Shutter (.u, v, k.) - FP_sh (.k.) 628.000 1 1 3 629.000 1 2 END 630.000 1 2 631.000 END; ( find_beam_eq ) 632.000

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19 PAGE 18 LINE NUMBER B P C I STMT # S 0 U R C E P R OG R AM V —-+ -- - +- 2 —---—. -3 —-- + —4 + 5 + 6 + 7 8 + 633.000 1

PASCAL/VS RELEASE 2.1 PIAN_CAL 10/10/83 11:28:19 PAGE 19 LINE NUMBER B P C I STMT $ 0 U R C E P R O G R A M V —-+- - -- ------ 2 --— + 3- — + --— 4 -+ 5 + —-6....+ -— 7....+ 8...+ 9..+ -— + 1V 634.000 635.000 636.000 /* 637.000 * Find_Beam Planes 638.000 639.000 * This calculates equations of the lines from the focal point 640.000 * of the shutter to the sample points on the shutter plane. 641.000 This is used for verif ng the L matrix calculated by 642.000 * geometrical consideration. 643.000 / 644 000 645.000 1 PROCEDURE FindBeam Planes; 646 000 1 647.000 1 VAR 648.000 1 at, b1, cl, { direction numbers on column 1 ) 649.000 1 an, bn, cn:'REAL; { direction numbers on column n ) 650.000 1 u: INTEGER; { Index ) 651.000 1 652.000 1 BEGIN 653.000 1 654.000 1 1 1 FOR u:- I TO Nu DO 655.000 1 1 BEGIN 656.000 11 1 2 at:= Beam (.u, 1, 1.); 657.000 1 1 1 b: Beam (.u, 1, 2.); 658.000 1 1 1 4 c:= Beam (.u, 1, 3.); 659.000 1 1 660.000 1 1 1 5 an:~ Beam (.u, Nv, 1.); 661.000 1 1 1 6 bn:~ Beam (.u, Nv, 2.); 662.000 1 1 1 7 cn:- Beam (.u, Nv, 3.); 663.000 1 1 664.000 i 1 1 8 Beam Plane (.u, 1.):~ bi * cn - bn * cl; { A ) 665.000 1 1 1 9 Beam Plane (.u, 2.):- c * an - cn * at; B ) 666.000 1 1 1 10 Beam Plane (.u, 3.):~ at * bn - an * bi; { C ) 667.000 1 I I 11 Beam_Plane (.u, 4.): - FP_sh(.1.) * Beam_plane (.u, 1.) - FPsh(.2.) * Beam_plane (.u, 2.) 668.000 1 1 1 - FPsh(.3.) * Beam_plane (.u, 3.); { D ) 669.000 1 1 END 670.000 1 1 671.000,END; 672.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 20 LINE NUMBER B P C I STMT # S O U R C E P R O G R A M V —-+ + 2 + 3 + 4 + 5 + 6 -— + 7 + 8 + 9 + —-1V 673.000 674.000 675.000 / 676.000 * Map_to_Camera 677.000 678.000 * This calcuates the sample point on the camera plane 679.000 * using the clip coordinates. This is the standard 680.000 * method used In graphics. 681.000 * The output Is checked against the output obtained 682.000 * using T matrix and the output obtained by geometrical 683.000 * consideration. 684.000 / 685.000 686.000 1 FUNCTION Map_to_Camera ( Pt: Vector; 687.000 1 VAR Xs, Ys: INTEGER ): BOOLEAN; 688.000 1 VAR 689.000 1 Clip: Vector; { clip coordinates ) 690.000 1 691.000 1 sum: REAL; { sum ) 692.000 1 i, k: INTEGER; { ndices ) 693.000 1 694.000' 1 BEGIN 695.000 1 696.000 1 find clip coordinates by multiplying Vlew_xform matrix ) 697.000 1 698.000 1 1 1 FOR i:= I TO 3 DO 699.000 1 1 BEGIN 700.000 1 1 2 sum:=;.701.000 1 1 2 3 FOR k:= I TO 3 DO 702.000 1 1 4 sum:= sum + Pt(.k.) * Viewxform (.k, i.); 703.000 1 1 5 Clip (.i.):= sum + View_xform (.4, 1.) 704.00 10 END; 705.QOO 1 706.000 1 { clip and find camera coordinates ) 707.000 1 1 6 IF ( - Cllp.(.3.) <= Clip(.1.) ) & ( Clip(.1.) <= Clip(.3.) ) & 708.000 1 1 ( - Cllp(.3.) <= Clip(.2.) ) & ( Clip(.2.) <= Clip(.3.) ) 709.000 1 1 THEN 710.000 1 1 BEGIN 711.000 I 1 1 7 Xs:= ROUND 4 ( Cllp(.1.) / Cl1p(.3.) ) * ( Nx / 2 ) + ( Nx / 2 + 0.5 ) ); 712.000 1 1 1 8 Ys: ROUND ( ( Clip(.2.) / Clip(.3.) ) * ( Ny / 2 ) + ( Ny / 2 + 0.5 ) ); 713.000 1 1 1 9 Map_to Camera:= TRUE 714.000 1 1 END 715.000 1 1 716.000 1 1 ELSE BEGIN 717.000 1 I 1 10 Map_to_camera:= FALSE; 718.000 I 1 11 Xs:r unobserved 719.000 1 1 1 12 s: unobserved 720.000 I I END 721.000 1 1 722.000 END; 723.000 724.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 21 LINE NUMBER B P C I STMT S 0 U R C E P R O G R A M ~V —--— 4- -I__~~ - -------- --— 8 — 3 —-----— ~~-4- 4 —— 9 —---- - 5 —- 6tV 7 725.000 726.000 727.000 /* 728.000 * Map_to_Cam Tmat 729.000 * 730.000 * This finds the sample points on the camera plane 731.000 * using the T matrix. 732.000 */ 733.000 734.000 1 FUNCTION Map_to_Cam_Tmat ( Pt: Vector; 735.000 1 VAR Xs, Vs: INTEGER ): BOOLEAN; 736.000 1 VAR 737.000 1 Homo: ARRAY(.1..4.) of REAL; ( homogeneous coordinates ) 738.000 I 739.000 1 sum: REAL; ( sum ) 740.000 1 1. k: INTEGER; ( indices ) 741.000 1 742.000 1 BEGIN 743.000 1 744.000 1 ( find homogeneous coordinates by multiplying Tmat matrix ) 745.000 1 746.000 1 1 1 FOR I: 1 TO 4 DO 747.000 1 1 BEGIN 748.000 1 1 1 2 sum:-; 749.000 1 I 2 3 FOR k:= I TO 3 DO 750.000 I 1 1 4 sum:= sum + Pt(.k.) * Tmat (.k, 1.); 751.000 1 1 1 5 Homo (.1.):- sum + Tmat (.4, 1.) 752.000 1 END; 753.000 1 754.000 1 { find screen coordinates ) 755.000 1 756.000 1 6 Xs:= ROUND ( Homo(.1.) / Homo(.4.) ); 757.000 7 Ys: ROUND ( Homo(.2.) / Homo(.4.) ); 758.000 1 759.000 1 { clip ) 760.000 1 761.000 1 1 8 IF ( 1 < Xs ) & ( Xs <- Nx ) 8 762.000 1 1 ( < Ys ) & ( Ys < Ny ) 763.000 1 1 9 THEN Map_to_Cam_Tmat:~ TRUE 764.000 1 1 ELSE BEGIN 765.000 1 1 1 10 Map_to_cam_Tmat:t FALSE 766.000 1 1 1 II Xs:- unobserved 767.000 1 I 12 Ys:= unobserved 768.000 i I END 769.000 1 1 770.000 END; 771.000 772.000 773.000

PASCAL/VS. RELEASE 2.1 PHAN_CAL 10/10/83 11128: 19 PAGE 22 LINE NUMBER B P C I STMT # S O U R C E P R O G R A M V —-+ -4.0 + 2003 4+ 5 + 6 7 + + - 774.000 775.000 776.000 1* 777.000 * Map_to_Cam_Geo 778.000 * 779.000 * This calculates sample points on the camera plane 780.000 * by the geometrical consideration. 781.000 */ 782.000 783.000 1 FUNCTION Map_to_Cam_Geo ( Pt: Vector 784.000 1 VAR Xs, Ys: INTEGER ): BOOLEAN 785.000 1 VAR 786.000 1 Ap, Bp, Cp, Dp, ( equation of the plane } 787.000 1 al, bl, cl: REAL; { equalton of line } 788.000 1 789.000 1 axe, aye: Vector; { unit vectors on the camera plane } 790.000 1 791.000 1 Cam_pt: Vector; i{ ntersection of the camera plane and the line } 792.000 1 793.000 1 S1, S2, t: REAL; { sums, temporary } 794.000 1 1: INTEGER; ({ ndex ) 795.000 I 796.000 1 BEGIN 797.000 1 798.000 1 { equation of camera plane: Ap x + Bp y + Cp z + Dp = 0 ) 799.000 1 800.000 1 1 Ap: Dlr_no_ca (.1.); 801.000 1 2 p:= Dir_no_ca (.2.); 802.000 3 Cp:= Dir_noca (.3.); 803.000 4 Dp:= - ( Ap * Center_ca (.1.) + Bp * Center_ca (.2.) + Cp * Centerca (.3.) ); 804.000 1 805.000 1 { direction numbers of a line from the point to the focal point of a camera ) 806.000 1 807.000 1 5 al: Pt (.1.) - FPca (.1.); 808.000 1 6 bl: Pt (.2.) - FPca (.2.); 809.000 1 7 cl: Pt (.3.) - FP_ca (.3.); 810.000 1 811.000 1 { Intersection of a camera plane and the line ) 812.000 1 813.000 1 8 t: - ( Ap * FP_ca(.1.) + Bp * FPca(.2.) + Cp * FP_ca(.3.) + Dp ) / 814.000 1 ( Ap * al + Bp * bl + Cp * cl); 815.000 1 816.000 1 9 Cam_pt (.1.):= FP_ca (.1.) + t * al; 817.000 1 10 Cam_pt (.2.):= FP_ca (.2.) + t * bi; 818.000 1 11 Cam_pt (.3.):= FP_ca (.3.) + t * cl; 819.000 1 820.000 1 ( find unit vectors In the camera plane ) 821.000 1 822.000 1 12 axe (.1.):= - SIN ( Thetaca ); 823.000 1 13 axe (.2.):= COS ( Thetaca ); 824.000 1 4 axe (.3.):; 825.000 1 826.000 t 15 aye (.1.):= - COS ( Phica ) * COS ( Thetaca );

PASCAL/VS RELEASE 2.1 PHAN CAL:MAP_TOCAM_GEO 10/10/83 11:28:19 PAGE 23 LINE NUMBER B P C I STMT # SOUR C E P.R O G R AM V —- --— + 1 + --— 2 --— + 3 + ---- 4 + 5 +- 6 + +- 8 + + V 827.000 1 16 aye (.2.):= - COS ( Ph1_ca ) + SIN ( Theta_ca ); 828.000 1 17 aye (.3.):= SIN ( Phi_ca ); 829.000 1 830.000 1 { find camera coordintes for the world coordinate point Campt ) 831.000 I 832.000 1 18 S1:; S2:; 833.000 1 1 20 FOR 1: 1I TO 3 DO 834.000 1 1 BEGIN 835.000 1 1 1 21 S:= S + ( Campt(.1.) - Center_ca(.t.) ) * axe (.1.); 836.000 1 1 1 22 S2: S2 + ( Cam_pt(..) - Center_ca(.l.) ) * aye (.1.) 837.000 1 END; 838.000 1 839.000 1 23 Xs:= ROUND ( ( S1 * ( Nx / Sx ) + I + Nx ) /2 ); 840.000 1 24 Ys: - ROUND ( ( - S2 * ( Ny / Sy ) + + Ny ) / 2 ) + Ny + 1; 841.000 1 842.000 1 ( Clip ) 843.000 1 844.000 1 25 IF ( 1 <= Xs )& ( Xs <- Nx ) & 845.000 1 ( 1 <- YS ) & ( Ys < Ny ) 846.000 1 1 26 THEN Map to_Cam_Geo: TRUE 847.000 1 ELSE BEGIN 848.000 1 1 1 27 Mapto_CamGeo: FALSE; 849.000 1 1 1 28 Xs:= unobserved 850.000 1 1 29 Ys:= unobserved 851.000. 1 END 852.000 1 853.000 END; 854.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 24 LINE NUMBER B P C I STMT S O U R C E P R O G R A M V —-+ 1 -- — + 2 -- + 3. -+ -— 4 + ---— 5 ---— +6...-+- -— 78 -+- — 9 + 1V 855.000 856.000 857.000 / 858.000 * MakeImage_pattern 859.000 860.000 * This makes a pattern of Images to be passed to the 861.000 * Surface Mapping Program. 862.000 * If two different points mapped Into the same point 863.000 * on the camera, Surface Mapping Program generates a 864.000 * spurious output. This routine delete all such points 865.000 * from the pattern and flags that point. 866.000 * 867.000 868.000 1 PROCEDURE Make_ImagePattern; 869.000 t 870.000 1 VAR 871.000 1 Intersect: Vector; { Intersection of the plane and the beam of the shutter ) 872.000 1 873.000 1 A, B, C, D: REAL; { Ax + By + Cz + D 0 0 ) 874.000 1 875.000 1 Xsl, Ysl, 876.000 1 Xs2, Ys2, 877.000 1 Xs3, Ys3: INTEGER; ( Screen coordinates } 878.000 1 879.000 1 t: REAL; ( temporary ) 880.000 1 dummy: BOOLEAN; ( dummy ) 881.000 1 u, v, k: INTEGER; { Indices ) 882.000 1 883.000 1 BEGIN 884.000 1 885.000 1 A:= Map_plane (.1.); 886.000 1 2 B: Mapplane (.2.); 887.000 1 3 C:= Mapplane (.3.); 888.000 1 4 D: Map_plane (.4.) 889.000 1 890.000 1 1 5 FOR u:=.1 TO Nu DO 891.000 1 2 6 FOR v:= i TO Nv DO 892.000 1 2 893.000 1 2 BEGIN 894.000 1 1 2 7 t:= - ( A * FP_sh(.1.) + B * FP_sh(.2.) + C * FP sh(.3.) + D ) / 895.000 1 1 2 ( A * Beam(.u, v, 1.) + B * Beam(.u, v, 2.) + C * Beam(.u, v, 3.) ); 896.000 1 1 2 897.000 1 1 3 8 FOR k: 1 TO 3 DO 898.000 1 1 2 9 Intersect (.k.):= FP sh (.k.) + t * Beam (.u, v, k.); 899.000 1 1 2 900.000 1 1 2 10 IF Mapto Camera ( Intersect, Xsl, Ysl ) = TRUE 901.000 1 1 2 THEN BEGIN 902.000 2 1 2 2 11 IF Screen (.Xsl, Ysl.) = FALSE 903.000 2 1 2 THEN BEGIN 904.000 3 1 2 2 12 WRITELN ( anaf, u:lenf, v:len, Xsl:len, Ysl:len, Intersect(.1.):w1d 905.000 3 1 2 2:pre, Intersect(.2.):w1d:pre, Intersect(.3.):w1d:pre ); 906.000 3 1 2 907.000 3 1 2 2 13 Screen (.Xsl, Ysl.):= TRUE;

PASCAL/VS RELEASE 2.1 PHAN_CAL:MAKE_IMAGE_PATTE 10/10/83 11:28:19 PAGE 25 LINE NUMBER B P C I STMT # S 0 U R C E P R O G R AM V — -+....1.. +.-..2 --— + --— 3 —-- + —-- 4 - + 5 +.6.+ 7 +.8 + 9...+ —- 1V 908.000 3 1 2 2 14 Pattern (.u, v, 1.): Xsl; 909.000 3 1 2 2 15 Pattern (.u, v, 2.):~ Ysl 910.000 2 1 2 2 END 911.000 2 2 2 ELSE BEGIN 912.000 2 1 2 2 Different beams map Into same camera point ) 913.000 3 1 22 16 Pattern (.u, v, 1.):= unobserved; 914.000 3 1 2 2 17 WRITELN (' **** Duplicate points —-', u:len, v:len ) 915.000 2 1 2 2 END 916.000 1 1 1 2 END 917.000 1 1 2 918.000 1 1 1 2 ELSE BEGIN 919.000 2 1 1 2 18 Pattern (.u, v, 1.):= unobserved; 920.000 2 1 1 2 19 WRITELN (' VMAT *** not observed ***', u, v ) 921.000 1 2 END; 922.000 i 1 2 923.000 1 1 2 ( comparison of 3 methods } 924.000 1 1 2 925.000 1 1 2 20 dummy:- Map_toCam_Geo ( Intersect, Xs2, YS2 ); 92G.000 1 1 2 21 dummy. Map_to_Cam_Tmat ( Intersect, Xs3, YS3 ); 927.000 1 I 2 928.000 1 1 1 2 22 If (Xsi <> Xs2) (Xs2 <> Xs3) I (Ysi <> Vs2) I (Ys2 <> Ys3) 929.000 1 1 1 2 THEN BEGIN 930.000 2 1 2 23 WRITE ('N#### Inconsistency', u:len, v:len ); 931.000 2 i 1 2 24 WRITELN (' VMAT', Xs1:len, YsI:len ); 932.000 2 1 1 2 25 WRITELN (' GEO', Xs2:len, Ys2:len ); 933.000 2 1 1 2 26 WRITELN (' Tmat', Xs3:len, Ys3:len ) 934.000 1 1 1 2 END 935.000.1 1 2 936.000 1 2 END {for u,v) 937.000 1 2 938.000 END; 939.000

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19 PAGE 26 LINE NUMBER B P C I STMT # S O U R C E P R O G R A M V —-+ -— i + —-- 2 — +. —-3 -- 6 +..7..8.+ 9 + —-1V 940.000 941.000 942.000 /* 943.000 * Output_pattern 944.000 945.000 * ++ HERE ++ Nu is a power of 2 946.000 947.000 * This generates patterns of images to be used 948.000 * by the Surface Mapping Program. 949.000 / 950.000 951.000 1 PROCEDURE Output_Pattern; 952.000 1 953.000 1 VAR 954.000 1 No_pattern, { Number of patterns ) 955.000 1 Column, { column on the shutter plane } 956.000 1 Run, ( run run run.........) 957.000 1 Pat: INTEGER; This goes 1, 2, 4, 8, 16,... ) 958.000 1 959.000 1 Mask val: BOOLEAN; shutter is open or close? ) 960.000 i 1, J. k, v: INTEGER; { Indices } 961.000 1 962.000 1 BEGIN 963.000 1 964.000 1 1 No_pattern':= ROUND ( LN ( Nu ) / LN ( 2 ) ) + 1; 965.000 1 966.000 1 2 Pat:= 1; This goes 1, 2, 4, 8, 16,.... 967.000 1 968.000 1 1 3 FOR I:= 1 TO Nopattern DO 969.000 1 1 BEGIN 970.000 1 1 971.000 1 1 1 4 Column:= Nu + 1; 972.000 1 1 1 5 Maskval:= TRUE; ( FALSE --- Closed shutter ) 973.000 1 1 1 ( TRUE --- Open shutter } 974.000 1 t 1 6 Run:= Nu DIV Pat; 975.000 1 1 976.000 1 1 2 7 FOR j:= 1 TO Pat DO 977.000 I1 2 BEGIN 978.000 1 1 2 979.000 2 1 3 8 FOR k:= 1 TO Run DO 980.000 2 1 3 BEGIN 981.000 3 1 3 9 Column:= Column -; 982.000 3 1 3 983.000 3 1 1 3 10 IF Mask_val = True ( open } 984.000 3 1 1 3 THEN 985.000 3 1 1 4 11 FOR v:= 1 TO Nv DO 986.000 3 I 2 4 12 IF Pattern (.Column, v, 1.)'= unobserved 987.000 3 1 2 4 13 THEN WRITELN ( Img_f, Pattern (.Column, v, 1.):lenf, 988.000 3 1 2 4 Pattern (.Column, v, 2.):len ) 989.000 2 1 2 END; (for k) 990.000 2 1 2 991.000 2 1 2 14 Mask_val:= NOT Mask_val 0QQ Onnn I I o

PASCAL/VS RELEASE 2.1 PHAN CAL:OUTPUT PATTERN 10/10/83 11:28:19 PAGE 27 LINE NUMBER B P C I STMT S O U R C E P R O G R A M V —-+. + 2 + + 4 + 5 — + 6 + —- 1V7+ 8 + 993.000 1 1 1 END; (for j) 994.000 1 1 995.000 1 1 1 15 Pat:= 2 * Pat; 996.000 1 1 1 997.000 1 1 1. 16 WRITELN ( Img_f, -1:lenf, -1:len ) 998.000 t I 1 999.000 1 1 END {for 1) 1000.000 1 1 1001.000 END; 1002.000 1003.000 1004.000

PASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:28:19 PAGE 28 LINE NUMBER B P C I STMT N SOURCE PROGRAM V —-+ —-7 i —------- 2 —------ 3 —-------- 1005.000 1006.000 1007.000 /* 1008.000 * ReadCalPara 1009.000 * 1010.000 * This reads parameters about the calibration cube. 1011.000 1012.000 1013.000 1 PROCEDURE ReadCal Para 1014.000 1 1015.000 1 VAR 1016.000 1 1017.000 1 R_mat, 4 Rotation matrix ) 1018.000 1 Rot, 4 temporary rotation matrix ) 1019.000 1 T_inv,( inverse translation matrix ) 1020.000 1 Tempmat: mat_4; temporary matrix 1021.000 1 1022.000 1 Axis: CHAR; ( axis of rotation of the cube ) 1023.000 i Angle, { amount of rotation of the cube In degree j 1024.000 1 Length: REAL; length of the side of the cube ) 1025.000 1 1026.000 1 1, J: INTEGER 4 Indices i 1027.000 1 1028.000 1029.000 1 BEGIN 1030.000 1 1031.000 1- 1 WRITELN 1032.000 1 2 WRITELN (' ------ Enter Calibration Parameters'; 1033.000 1 3 WRITELN 1034.000 1 1035.000 1 4 WRITELN (' length of the side of a cube -- Length' ) 1036.000 I 5 READLN ( Length; 1037.000 1 1038.000 1 6 WRITELN (' Step In u and v -- Step' ); 1039.000 1 7 READLN ( Step; 1040.000 1 1041.000 1 8 a:~ Length / 2; 1042.000 1043.000 1 Rotation of the cube about x, y, and z axis ) 1044.000 1 1045.000 1 9 WRITELN (' ---- rotation of the cube ----' 1046.000 1 10 WRITELN; 1047.000 1048.000 II WRITELN (' Enter axis ( x, y, z, or / ) and angle' 1049.000 1 12 READLN ( Axis, Angle ); 1050.000 113 Angle:= Angle * 3.141592 / 180; 1051.000 1052.000 1 Find rotation matrix Rmat ) 1053.000 i 1054.000 1 14 Make_Identity ( Rmat; 1055.000 I 1056.000 1 15 WHILE ( Axis ='x' )I ( Axis'y' i ( yAxis ='' ) DO 1057.000 1 1 RFGTN

PASCAL/VS RELEASE 2.1 PHANCAL:READ_CAL_PARA 10/10/83 11:28:19 PAGE 29 LINE NUMBER B P C I STMT # SOU R C E P R O G R AM V —-+ —- 1 — -+ —— 2 —— + —— 3 —— + —— 4 + 5 + 6 + 7 S-+ 8...+ 1V 1058.000 1 1 1059.000 1 1 1 16 Make_Identity ( Rot); 1060.000 1 I 1 1061.000 11 11 17 CASE Axis OF 1062.000 1 1 1 1063.000 1 1 I 1'x': BEGIN 1064.000 2 1 1 1 18 Rot(.2.2.): COS (Angle); 1065 000 2 1 1 1 19 Rot(.2, 3.): - SIN (Angle); 1066.000 2 1 1 1 20 Rot(.3, 2.): SIN (Angle); 1067.000 2 1 1 1 21 Rot(.3, 3.):- COS (Angle) 1068.000 1 1 1 END 1069.000 1 1 1 1 1070.000 1 1 1'y': BEGIN 1071.000 2 1 1 1 22 Rot(.1, 1.):'COS (Angle); 1072.000 2 1 1 1 23 Rot(.1, 3.): SIN (Angle); 1073.000 2 1 1 24 Rot(.3, 1.): - SIN (Angle); 1074.000 2 1 1 1 25 Rot(.3, 3.):- COS (Angle) 1075.000 1 1 1 END; 1076.000 1 1 1 1 1077.000 1 1 1'z': BEGIN 1078.000 2 1 1 1 26 Rot(.1, 1.):' COS (Angle); 1079.000 2 1 1 1 27 Rot(.l, 2.):- - SIN (Angle); 1080.000 2 1 1 28 Rot(.2, 1.):~ SIN (Angle); 1081.000 2 1 1 1 29 Rot(.2, 2.):- COS (Angle) 1082.000 1 1 1 END 1083.000 1 1 1 1 1084.000 1 I 1 END; case) 1085.000 1 1 1 1086.000 1 1 1 30 Tempmat:= R_mat; 1087.000 1 1 1 31 Mult_Mat4 ( R_mat, Tempmat. Rot ); 1088.000 1 1 1 1089.000 1 I 1 32 WRITELN (' Enter axis ( x, y, z, or / ) and angle' ); 1090.000 1 1 1 33 READLN ( Axis, Angle ); 1092.000 1 I 34 Angle: Angle * 3.141592 /180 1093.000 1 1 1 1094.000 1 END; (while) 1095.000 1 1096.000 1 1097.000 1 ( Translation of the cube ) 1098.000 1 1099.000 1 35 WRITELN (' - Translation of the cube ----' ) 1100.000 1 36 WRITELN; 1101.000 1 1102.000 1 37 WRITELN (' Enter translation along x, y, and z axs' ); 1103.000 1 38 WRITELN (' e.g. 3.8 10 8.5'); 1104.000 1 39 REAULN ( xt, yt, zt ); 1105.000 1 1106.000 1 1107.000 1 40 Make_ldentity ( Tnv ); 1108.000 1109.000 I 41 T_1nv (.4, 1.): - xt 1110.000 1 42 T nv (.4, 2.): - yt 1111.000 1 43 T nv (.4, 3.):' - zt;

PASCAL/VS RELEASE 2.1 PHAN_CAL:READ_CAL_PARA 10/10/83 11:28:19 PAGE 30 LINE NUMBER B P C I STMT # S OU R C E P R O G R AM V —-+ --—....1.. + —-- 2 ---- + 3 ——.. + ---- 4 —-- ---- 5 —— 6 —--. 7. --- 8 -+ 9 -+ —--- 78V 1112.000 1 1113.000 1 { find H mat and Inverse H_1nv ) 1114.000 1 1115.000 1 44 H mat:= R mat 1116.000 1 45 Hmat (.4, 1.):= xt 1117.000 1 46 H_liat (.4, 2.):= yt 1118.000 1 47 H_mat (.4, 3.):= zt 1119.000' 1 1120.000 1 48 FOR I:= 1 TO 4 DO 1121.000 1 2 49 FOR J:= 1 TO 4 DO 1122.000 1 50 Temp_mat (.1, j.):= R_mat (.J, 1.); { transpose ) 1123.000.1 1124.000 1 51 Mult_Mat4 ( H_lnv, T_tnv, Temp_mat ) 1125.000 1 1126.000 END

pASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 31 LINE NUMBER 8 P C I STMT # S 0 U R C E P R O G R A M V --— 1 + 2 + -— 3 --— + ------— 35 ------------ --— 6 + ---— 7 + --— 8 + 9 -+ —-1V 1127.000 1128.000 1129.000 / 1130.000 * Find Cal_Face 1131.000 1132.000 * This finds the equation the faces of the cube. 1133.000 + Then this test If the face is visible by shutter and 1134.000 * camera, and output the number of faces visible by shutter 1135.000 * and camera. 1136.000 / 1137.000 1138.000 1 PROCEDURE F1ndCalFace; 1139.000 1 1140.000 1 VAR 1141.000 1 1142.000 1 V. { vertices of the cube ) 1143.000 1 Vt: ARRAY (.1..8, 1..4.) of REAL; moved vertices of the cube) 1144.000 1 F, faces of the cube ) 1145.000 1 Ft: ARRAY (.1..6, 1..4.) of REAL; ( moved faces of the cube ) 1146.000 1 11471000 1 signl, sign2, sign3: REAL; ( signs ) 1148.000 1 1149.000 1, J, k, n: INTEGER; ( indices ) 1150.000 1151.000 1 BEGIN 1152.000 1153.000 1 position of the vertices before the cube is moved ) 1154.000 1 1155.000 1 V(., 1.) a;; V(.1, 2.) 3.): a; V(.1, 4.): 1156.000 1 5 V(.2,.): a; V(.2, 2.): a; V(.2. 3.): - a; V(.2, 4.):; 1157.000 1 9 V(.3,.): a; V(.3, 2.): - a; V(.3, 3.): - a; V(.3, 4.):; 1158.000 1 13 V(.4, 1.): a; V(.4, 2.):- - a; V(.4, 3.):= a; V(.4. 4.):; 1159.000 1 17 V(.5. 1.): -a; V(.5, 2.): a; V(.5, 3.): a; V(.5, 4.): 1; 1160.000 1 21 V(.6,.): -a; V(.6, 2.): a; V(.6, 3.): - a; V(.6, 4.):; 1161.000 1 25 V(.7 1.): - a V(.7, 2.):- - a; V(.7, 3.): - a; V(.7, 4.):; 1162.000 1 29 V. 1.): - a V(.8, 2.): - a;V(.8. 3.): a; V(.8, 4.): 1163.000 1 1164.000 1 equation of faces of a cube before the cube is moved ) 1165.000 1 1166.000 33 F(.1, 1.): 0; F(.1, 2.):0 F(.. 3.); F(.1, 4.): - a 1167.000 1 37 F(.2,.): - 1 (.2,.2 2.):= 0; F(.2, 3.): 0; F(.2. 4.): -; 1168.000 1 41 F(.3,.):; F(.3, 2.): 0; F(. 3.): - 1;F(.3. 4.): -; 1169.000 145 F(.4.: 1; (.4 2.): 0; F(.4, 3.):= 0; F(.4, 4.): -a; 1170.000 1 49 F(.5 1.):; F(.5, 2.): - I; F(.5, 3.): 0; F(.5, 4.): - a; 1171.000 53 F(.6. 1.): 0; F(.6, 2.): 1; F(.6, 3.): 0; F(.6, 4.): - a; 1172.000 1 1173.000 1 Move Vertices to new position -- Vt(.i, 4.) is I 1174.000 1 1175.000 1 1 57 FOR I:= I TO 8 DO 1176.000 1 2 58 FOR J: 1I TO 4 DO 1177.000 1 2 BEGIN 1178.000 1 t 2 59 Vt(.t, j.):0; 1179.000 1 1 3 60 FOR K:= I TO 4 DO

PASCAL/VS RELEASE 2.1 PHAN_CAL:FIND_CAL_FACE 10/10/83 11:28:19 PAGE 32 LINE NUMBER B P C I STMT # S OUR C E P R O G R A M V —-+ ---— + 2 -— +- — 3 - -+ 4 + --- 5 --— + --— 6 +.7 + 8 -— + 9 + 1V 1180.000 1 3 61 Vt(.1. j.):= Vt(.1, J.) + V(.1, k.) * H_mat(.k, j.) 1181.000 1 END; 1182.000 1 1183.000 1 { Move Faces to new position ) 1184.000 1 1185.000 1 1 62 FOR I:= TO 6 DO 1186.000 1 2 63 FOR J:= 1 TO 4 DO 1187.000 1 2 BEGIN 1188.000 1 1 2 64 Ft(.i, j.):= 0 1189.000 1 I 3 65 FOR K:= I TO 4 DO 1190.000 1 1 3 66 Ft(.i, j.):= Ft(.1, J.) + H_tnv(.J, k.) * F(.1, k.) 1191.000 1 END; 1192.000 1 1193.000 1 1194.000 1 { find faces which are visible by both shutter and camera ) 1195.000 1 1196.000 1 67 n:= 0 1197.000 1 1198.000 1 1 68 FOR 1: 1 TO 6 DO 1199.000 1 1 BEGIN 11200.000 I I 1201.000 i 1 By the sign of signi, we know at which side of the face focal point of shutter s 1202.00Q 1 1 1203.000 1 1 1 69 signl:= 0; 1204.000 1 1 2 70 FOR k: I TO 3 DO 1205.000 1 1 1 71 signl:= signl + FP_ca(.k.) * Ft(.i, k.); 1206.000 1 1 72 signl:= signl + Ft(.,1 4.); 1207.000 1 1 1 1208.000 1 1 1 73 sign2:= 0; 1209.000 1 1 2 74 FOR k:= 1 TO 3 DO 1210.000 1 1 1 75 slgn2:= sign2 + FP_sh(.k.) * Ft(.i, k.); 1211.000 1 1 1 76 sign2: sign2 + Ft(.1, 4.); 1212.000 I 1 1 1213.000 1 1 1 77 slgn3:= 0; 1214.000 1 i 2 78 FOR K: 1 TO 4 DO 1215.000 1 t 1 79 slgn3:- sign3 + Vt(.l+1, k.) * Ft(.1. k.); 1216.000 1 1 1 1217.000 1 I 1 { focal points of shutter and camera should be at the same side of the face 1218.000 1 i I { and they should be at the opposite side of the face to a vertex selected ) 1219.000 1 1 1 1220.000 1 1 1 1 80 IF ( ( sign < 0 ) & ( sign2 < 0 )& ( sign3 > 0 ) ) I 1221.000 1 1 1 1 ( ( signl > 0 ) & ( sign2 > 0 ) & ( sign3 < 0 ) ) THEN 1222.000 1 i 1 1 BEGIN 1223.000 2 1 1 I 81 n:= n +; 1224.000 2 1 12 82 FOR k:= I TO 4 DO 1225.000 2 1 1 2 83 Face (.n. k.):= Ft(.l, k.) 1226.000 1 1 1 END; 1227.000 1 1 1 1228.000 1 END;(for) 1229.000 1 84 Num face:= n; 1230.000 1 1231.000 1 85 WRITELN; 1232.000 1 86 WRITELN (' Num face', Num face )

PASCAL/VS RELEASE 2.1 PHAN CAL:FIND CAL FACE 10/10/83 11:28:19 PAGE 33 LINE NUMBER B P C I STMT n S 0 U R C E P R 0 G R A M V —-+ —-----— + —— 2 —— + —— g —----— 4 —— + —— 5 —-----— 6 —— 0 —-— 7 —------- -t 1233.000 1 1234.000 END; 1235.000 1236.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 34 LINE NUMBER B P C I STMT # SOUR C E PROGRAM V-~-+ —-- I —— + —-- 2 —— + —— 3 —— + —— 4 —----— 5 —----— 6 —----— 7 —------- 1237.000 1238.000 1239.000/* 1240.000 * Output_CubeImage 1241.000 * 1242.000 * This outputs sample points of the cube on the 1243.000 * camera plane. 1244.000 1245.000 1246.000 1 PROCEDURE Output_Cube_Image; 1247.000 1248.000 1 VAR 1249.000 1250.000 1 A. B, C, D. { equation of the faces ) 1251.000 1 bound, ( square of length from the center of the cube to the vertex ) 1252.000 1 dist, ( distance ) 1253.000 1 t: REAL; i temporary ) 1254.000 1 1255.000 1 Intersect: Vector; i ntersection of the face and the beam of the shutter ) 1256.000 1 1257.000 1 u, v, 1, k, 1258.000 1 Xs, Vs: INTEGER; 4 indices ) 1259.000 1 1260.000 1 BEGIN 1261.000 1 1262.000 1 1 IF Numface <= 0 THEN 1263.000 1 1 BEGIN 1264.000 1 1 1 2 WRITELN (' ** Calibration cube Is not visible *****' ); 1265.000 1 1 3 RETURN 1266.000 1 1 END; 1267.000 1 1268.000 1 4 bound: 3 * a * a 1269.000'1 1270.000 1 5 u: i 1271.000 1 1 6 WHILE ( u <- Nu ) DO 1272.000 1 I BEGIN 1273.000 I 1 1 7 v:= 1 1274.000 1 1 2 8 WHILE v <= Nv DO 1275.000 1 1 2 BEGIN 1276.000 1 1 2 1277.000 2 1 3 9 FOR 1: i TO Num face DO 1278.000 2 1 3 BEGIN 1279.000 2 1 3 1280.000 2 1 3 4 find Intersection of beam(.u, v, *.) and Face ) 1281.000 2 1 3 1282.000 3 1 3 10 A Face (.1. 1.) 1283.000 3 1 3 ii B: Face (.1, 2.); 1284.000 3 1 3 12 C: Face (.1. 3.) 1285.0003 1 3 13 0: Face (.1. 4.); 1286.000 3 1 3 1287.000 3 1 3 14 t: - ( A * FPsh(.i.) + B * FPsh(.2.) + C * FPsh(:3.) + D ) / 1288.000 3 1 3 (A * Beam(.u, v, 1.) + B * Beam(.u, v, 2.) + C * Beam(.u, v, 3.) ); 1289.000 3 1 3

PASCAL/VS RELEASE 2.1 PHAN CAL:OUTPUT CUBE IMAG 10/10/83 11:28:19 PAGE 35 LINE NUMBER B P C I STMT S 0 U R C E P R 0 R A M V -----— + —- ---— 8 —— 2 -------— 3 —-1 —--- - 1290.000 3 1 4 15 FOR k:- 1 TO 3 DO 1291.000 3 1 3 16 Intersect (.k.):~ FP_8h (.k.) + t * Beam (.u, v, k.) 1292.000 3 1 3 1293.000 3 1 3 17 diet:= SOR( Intersect(.1.) - xt ) + SQR( Intersect(.2.) - yt ) + 1294.000 3 1 3 SOR( Intersect(.3.) - zt ); 1295.000 3 1 3 1296.000 3 1 1 3 18 IF ( dist <- bound ) THEN ( Inside cube ) 1297.000 4 2 3 19 BEGIN IF Map_to Camera ( Intersect, Xs, Ys ) Y TRUE 1298.000 4 1 2 3 THEN BEGIN 1299.000 5 1 2 3 20 WRITELN ( calf, Intersect(.1.):wldf:3, Intersect(.2.):wld:3. 1303.000 5 1 2 3 Intersect(.3.):wld:3 ) 1301.000 5 1 2 3 1302.000 5 1 2 3 21 WRITELN ( calf, u:lenf, Xs:len, Ys:len ); 1303.000 4 1 2 3 END 1304.000 3 1 1 3 END 1305.000 3 1 1 3 1306.000 2 2 END;(for) 1307.000 2 1 2 1308.000 2 1 2 22 v:w v + step 1309.000 1 1 1 END;(while v) 1310.000 1 1 1 1311.000 1 1. 1 23 u:= u + step 1312.000 1 END (while u) 1313.000 1 1 1314.000 END; 1315.000 1316.000 1317.000

PASCAL/VS RELEASE 2.1 PHANCAL 10/10/83 11:28:19 PAGE 36 LINE NUMBER B P C I STMT i SOUR C E PROGRAM V —-+ —-7 --------— 2 —-+ —--- 3 —------ 1 1318.000 1319.000 1320.000 /****+*** ****** **** 1321.000 * 1322.000 * MAIN 1323.000 * 1324.000 1325.000 * 1326.000 1327.000 1328.000 1329.000 BEGIN 1330.000 1331.000 4 Open flies ) 1332.000 1333.000 1 RESET ( INPUT,'UNIT=SCARDS, INTERACTIVE ); 1334.000 2 REWRITE ( Calf,'UNIT3' ); 1335.000 3 REWRITE ( lmg_f,'UNIT~4' ); 1336.000 4 REWRITE ( Anaf,'UNIT=2' ); 1337.000 1338.000 1339.000 5 ReadParameters 1340.000 6 ReadCal Para; 1341.000 7 Initialize; 1342.000 1343.000 8 FindTmat 1344.000 9 FInd_-Imat 1345.000 1346.000 10 Find Ray_Eqs; 1347.000 11 Find_Beam_Eqs; 1348.000 12 FindBeamPlanes; 1349.000 1350.000 13 Make IlmagePattern; 1351.000 14 OutputPattern 1352.000 1353.000 15 FindCal_Face; 1354.000 Output_CubeImage 1355.000 1356.000 END. 1357.000 NO COMPILER DETECTED ERRORS OPTIONS IN EFFECT: MARGINS( 1,100), NOSEQUENCE, LANGLVL(EXTENDEO), LINECOUNT(60), PAGEWIDTH(132), CHECK, DEBUG, GOSTMT, SOURCE, WARNING SOURCE LINES: 1356; COMPILE TIME: 1.16 SECONDS; COMPILE RATE: 70004 LPM SOURCE LINES: 1356; TRANSLATE TIME: 0.71 SECONDS; TRANSLATE RATE: 115277 LPM TOTAL TIME: 2.33 SECONDS; TOTAL RATE: 34963 LPM

27 8.2. Eye.Cal Program

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:37:46 PAGE P001 C 1.000 C Mapping Program: Eye.Cal 2.000 C 3.000 C 4.000 C 5.000 C Version: This is based on the program of Dr. Martin Altschuler 6.000 C Calibration of T matrix and L matrix are done 7.000 C by least square method 8.000 C 9.000 C 10.000 C 11.000 C 12.000 C 13.000 C DIRECTORY 14.000 C SROW # OF ROWS OF CONTOUR GRID 15.000 C SCOL # OF COLUMNS OF CONTOUR GRID 16.000 C CROW N OF ROWS IN ENTIRE RASTER SCREEN 17.000 C CCOL # OF COLUMNS IN ENTIRE RASTER SCREEN 18.000 C NIMAGE # OF INAGES NEEDED TO CODE SCENE 19.000 C 20.000 C LUCAL LOGICAL UNIT FOR CALIBRATION INFO 21.000 C LUDAT LOGICAL UNIT FOR INPUT DATA 22.000 C LUOUT LOGICAL UNIT FOR FILE OUTPUT 23.000 C LUTERM LOGICAL UNIT FOR TERMINAL OUTPUT 24.000 C 25.000 C PATI MATRIX OF RASTER IMAGE # 1 26.000 C 27.000 C: 28.000 C PAT7 MATRIX OF RASTER IMAGE # 7 29.000 C 30.000 C 31.000 C LOGICAL DEVICES 32.000 C 33.000 C 34.000 C LOGICAL UNIT #3 CALIBRATION INFO 35.000 C LOGICAL UNIT #4 INPUT FILE 36.000 C LOGICAL UNIT #5 OUTPUT FILE 37.000 C LOGICAL UNIT #6 TERMINAL 38.000 C 39.000 C 40.000 C —-----------------— r —--------------- 41.000 C 42.000 C Main program 43.000 C 44.000 C —----- ------------- ------------------- 45.000 C 46.000 C Local variables 47.000 C 48.000 C Transformation matrices 49.000 C 50.000 0001 REAL T(4,4). L(4,4) 51.000 C 52.000 C Calibration cube data 53.000 C.9 A

MICHIGAN TERMINAL SYSItM FUNHIAN UG21.8I MAIN 10-10-83 11:37:46 PAGE P002 0003 INTEGER CUBIMG(100,2), CUBCOL(100) 56.000 C Begin 58.000 C 59.000 C Acquire and check the initial parameters. 60.000 0004 CALL INIT ( CUBSUR, CUBIMG, CUBCOL ) 61.000 C 62.000 C Calibrate the system. 63.000 0005 CALL TMAT( CUBSUR. CUBIMG, T ) 64.000 0006 CALL LMAT( CUBSUR. CUBCOL, L ) 65.000 C 66.000 C Solve for the unknown surface. 67.000 0007 CALL IMGIN 68.000 0008 CALL MAP( T. L ) 69.000 C 70.000 0009 STOP 71.000 0010 END 72.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK.LOADNOMAP *OPTIONS IN EFFECT* NAME = MAIN LINECNT a 57 *STATISTICS* SOURCE STATEMENTS = 10,.PROGRAM SIZE = 2916 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:37:46 PAGE P001 C 73.000 C —------------------------------------------ ~.~74.000 C 75.000 0001 SUBROUTINE INIT ( CUBSR, CUBIM, CUBCO ) 76.000 C 77.000 C —---------------------------------------------------- 78.000 C 79.000 C include GLOBAL and 10O common 80.000 C 81.000 C "GLOBAL",common (100,120) 100.000 C 101.000 0002 INTEGER SROW, SCOL 102.000 C Dimensions of the shutter grid 103.000 C 104.000 0003 INTEGER CROW, CCOL 105.000 C Dimensions of the camera grid 106.000 C 107.000 0004 INTEGER NIMAGE, NUMCAL 108.000 C Number of patterns needed 109.000 C Number of calibration data 109.500 C 110.000 0005 COMMON/GLOBAL/ 111.000 1 SROW, SCOL, CROW, CCOL, NIMAGE, NUMCAL 112.000 C 113.000 C "IO" common (200, 210 ) 200.000 0006 INTEGER LUCAL, LUDAT., LUOUT, LUTERM 201.000 C Logical units for the calibration file, data, output, 202.000 C and the terminal. 203.000 C 204.000 0007 COMMON/IO/ 205.000 1 LUCAL. LUDAT, LUOUT, LUTERM 206.000 C 207.000 C 208.000 C 83.000 C parameters -- calibration cube data 84.000 C 86.000 0008 REAL CUBSR(100, 3) 86.000 0009 INTEGER CUBIM(100, 2), CUBCO(100) 87.000 C 88.000 C Begin initialization 89.000 C 90.000 0010 LUCAL = 3 91.000 0011 LUDAT - 4 92.000 0012 LUOUT a 5 93.000 0013 LUTERM a 6 94.000 C 95.000 0014 CALL CALIN( CUBSR, CUBIM, CUBCO ) 96.000 C 97.000 0015 NIMAGE = INT(1.001 + ALOG(FLOAT(SCOL))/ALOG(2.)) 98.000 C 99.000 0016 RETURN 100.000 0017 END 101.000 *OPTIONS IN EFFECT* ID.EBCDICSOURCE.NOLIST,NODECKLOADNOMAP *OPTIONS IN EFFECT* NAMF = INIT I TmJJrMT a r.

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-1,O-83 11:37:47 PAGE POO0 C 102.000 C************************************************************** O103.000 C 104.000 C File I/O Module 105.000 C 106.000 C This module consists of procedures to handle the input/output 107.000 C activities for the surface acquisition program. 10I.000 C Include procedures are 109.000 C 110.000 C Callbin- reads all the data from the calibration file 111.000 C Imagein- reads the Image data 112.000 C Surfout- writes the surface points 113.000 C 114.000 C**+**********+************** ****** ********************************* 115.000 C 116.000 C —-------------------------------------------- 117.000 C 118.000 0001 SUBROUTINE CALIN (CUBS, CUBI, CUBC) 119.000 C * 120.000 C —-------------------------------------—.. —------------ 121.000 C 122.000 C INPUT ROUTINE FOR the CALIBRATION file 123.000 C 124.000 C INCLUDE Global and 10 commons 125.000 C "GLOBAL" common (100,120) 100.000 C 101.000 0002 INTEGER SROW, SCOL 102.000 C Dimensions of the shutter grid 103.000 C 104.000 0003 INTEGER CROW, CCOL 105.000 C Dimensions of the camera grid 106.000 C 107.000 0004 INTEGER NIMAGE, NUMCAL 108.000 C Number of patterns needed 109.000 C Number of calibration data 109.500 C 110.000 0005 COMMON/GLOBAL/ 111.000 1 SROW, SCOL, CROW, CCOL, NIMAGE, NUMCAL 112.000 C 113.000 C "10" common (200, 210 ) 200.000 0006 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 201.000 C Logical units for the calibration file, data, output, 202.000 C and the terminal. 203.000 C 204.000 0007 COMMON/10/ 205.000 1 LUCAL, LUDAT, LUOUT, LUTERM 206.000 C 207.000 C 208.000 C 127.000 C 128.000 C Declare parameters 129.000 0008 REAL CUBS(100,3) 130.000 0009 INTEGER CUBI(00,2), CUBC(100) 131.000 C'o% nnn

*STATISTICS* NO DIAGNOSTICS GENERATED

0010 READ(LUCAL, 2000) SROW, SCOL, CROW, CCOL 134.000 0011 2000 FORMAT ( 216 / 216 ) 135.000 C 136.000 C Now read the data for the calibration cube 137.000 C 138.000 0012 CALL CLRTBL( CUBS, 100, 3 ) 139.000 0013 CALL CLRTBL( CUBI. 100, 2 ) 140.000 0014 CALL CLRVC( CUBC, 100 ) 141.000 C 142.000 0015 NUMCAL = 0 143.000 0016 DO 100 I = 1,100 144.000 0017 READ(LUCAL, 300, END = 2100) ( CUBS(I,J), J = 1, 3 ) 145.000 0018 READ(LUCAL, 400) ( CUBC(I), CUBI(I,1), CUBI(I,2) ) 146.000 0019 NUMCAL = NUMCAL + 1 147.000 0020 100 CONTINUE 148.000 0021 2100 CONTINUE 148.500 C 148.550 C If there are too many data, stop 148.600 C 148.700 0022 IF (( NUMCAL.LE. 100 ).AND. ( NUMCAL.GE. 7 ) ) GO TO 4000 148.800 0023 WRITE(-LUTERM, 3015) NUMCAL 148.900 0024 STOP 148.950 0025 4000 CONTINUE 148.970 C 149.000 C WRITE(LUTERM, 3000) ((CUBS(I,d),Jd1,3),I=1,100) 150.000 C WRITE(LUTERM, 3005) ((CUBI(I,J),J-1,2), I=1,100) 151.000 C WRITE(LUTERM, 3010) (CUBC(I),I-,.100) 152.000 C 153.000 0026 300 FORMAT( 3F14.3 ) 154.000 0027 400 FORMAT( 316 ) 155.000 C 156.000 0028 3000 FORMAT(99(3F10.2,/),3F10.2) 157.000 0029 3005 FORMAT(99(2I10,/),2110) 158.000 0030 3010 FORMAT(100(I10,/)) 159.000 0031 3015 FORMAT(IX,'Number of calibration data is not adequate', 159.500 1' NumCal =', 110 ) 159.700 C 160.000 0032 RETURN 161.000 0033 END 162.000 *OPTIONS IN EFFECT* IDEBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = CALIN, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 33,PROGRAM SIZE ~ 966 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:37:47 PAGE POO1 C 163.000 C ----------— 1 —-----— 64.000 C 165.000 0001 SUBROUTINE IMGIN 166.000 C 167.000 C —---------------------------------------------------------- 68.000 C 169.000 C INPUT ROUTINE FOR RASTER DATA 170.000 C 171.000 C COMMONS Included here 172.000 C "GLOBAL* common (100.120) 100.000 C 101.000 0002 INTEGER SROW, SCOL 102.000 C Dimensions of the shutter grid 103.000 C 104.000 0003 INTEGER CROW, CCOL 105.000 C Dimensions of the camera grid 106.000 C 107.000 0004 INTEGER NIMAGE, NUMCAL 108.000 C Number of patterns needed 109.000 C Number of calibration data 109.500 C 110.000 0005 COMMON/GLOBAL/ 111.000 SROW, SCOL, CROW, CCOL, NIMAGE, NUMCAL 112.000 C 113.000 C "10" common (200, 210 ) 200.000. 0006 INTEGER LUCAL, LUDAT., LUOUT. LUTERM 201.000 C Logical units for the calibration file, data, output, 202.000 C and the terminal. 203.000 C 204.000 0007 COMMON/10/ 205.000 1 LUCAL, LUDAT, LUOUT, LUTERM 206.000 C 207.000 C 208.000 C "Patterns" common (300.310) 300.000 0008 INTEGER PAT1(128,128), PAT2(128,128), PAT3(128, 128). 301.000 1 PAT4(128.128). PAT5(128,128). PAT6(128,128), 302.000 2 PAT7(128.128), PAT8(128,128) 303.000 C Contain the Image data. 304.000 C 305.000 0009 COMMON/PATS/ 306.000 1 PATI, PAT2, PAT3, PAT4, PAT5, PAT6. 307.000 2 PAT7, PAT8 308.000 C 174.000 C 175.000 0010 CALL CLRTBL( PATI, 128. 128 ) 176.000 0011 CALL CLRTBL( PAT2. 128, 128 ) 177.000 0012 CALL CLRTBL( PAT3, 128. 128 ) 178.000 0013 CALL CLRTBL( PAT4, 128, 128 ) 179.000 0014 CALL CLRTBL( PAT5. 128. 128 ) 180.000 0015 CALL CLRTBL( PAT6. 128. 128 ) 181.000 0016 CALL CLRTBL( PAT7. 128, 128 ) 182.000 0017 CALL CLRTBL( PATS, 128, 128 ) 183.000 C A - A

C TO SIGNAL END OF PATTERN 186.000 C 187.000 0018 1 READ (LUDAT, 2000, END=100) I,J 188.000 0019 IF ((I.EQ.-1).AND. (J.EQ.-1)) GO TO 2 189.000 0020 PATI(I,J) = 1 190.000 0021 GO TO 1 191.000 C 192.000 0022 2 READ (LUDAT, 2000,END=100) I,J 193.000 0023 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 3 194.000 0024 PAT2(I,J) = 1 195.000 0025 GO TO 2 196.000 C 197.000 0026 3 READ (LUDAT, 2000, END=100) I,d 198.000 0027 IF ((I.EQ.-1).AND. (J.EQ.-1)) GO TO 4 199.000 0028 PAT3(I,J) = 1 200.000 0029 GO TO 3.201.000 C 202.000 0030 4 READ (LUDAT, 2000, END=100) 1,d 203.000 0031 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 5 204.000 0032 PAT4(I,J) = 1 205.000 0033 GO TO 4 206.000 C 207.000 0034 5 READ (LUDAT, 2000, END"100) I,J 208.000 0035 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 6 209.000 0036 PAT5(I,d) = 1 210.000 0037 GO TO 5 211.000 C 212.000 0038 6 READ (LUDAT, 2000, END=100) I., 213.000 0039 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 7 214.000 0040 PAT6(I,.) - 1 215.000 0041 GO TO 6 216.000 C 217.000 0042 7 READ (LUDAT, 2000, END=100) I,d 218.000 0043 IF ((I.EQ. -1).AND. (d.EQ.-1)) GO TO 8 219.000 0044 PAT7(I.d) = 1 220.000 0045 GO TO 7 221.000 C 222.000 0046 8 READ (LUDAT, 2000, ENDi100) I,J 223.000 0047 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 100 224.000 0048 PAT8(1,J) = 1 225.000 0049 GO TO 8 226.000 C 227.000 C1io WRITE( LUTERM, 3000 ) (( PATI(I,d), d1,CCOL), I11,CROW ) 228.000 C 229.000 0050 100 RETURN 230.000 0051 2000 FORMAT( 216 ) 231.000 C3000 FORMAT( 1612 ) 232.000 C 233.000 0052 END 234.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = IMGIN, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 52,PROGRAM SIZE s 1702 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:37:48 PAGE PiO1 C 235.000 C 236.000 C- - - --------------------------- 237.000 0001 SUBROUTINE SUROUT( SURPNT ) 239.000 C 241.000 C OUTPUT ROUTINE FOR ESTIMATED LOCATION 242.000 C OF POINTS 243.000 ----— C —------------------------------------------------ 244.000 C 245.000 C 10 Commons Included here 246.000 C "1O" common (200, 210 ) 200.000 0002 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 201.000 C Logical units for the calibration file, data, output, 202.000 C and the terminal. 203.000 C 204.000 0003 COMMON/IO/ 205.000 LUCAL, LUDAT, LUOUT, LUTERM 206.000 C 207.000 C 208.000 C 248.000 C 249.000 C PARAMETERS 250.000 0004 REAL SURPNT(3) 251.000 C 252.000 0005 WRITE (LUOUT, 2000) SURPNT(1), SURPNT(2), SURPNT(3) 253.000 0006 2000 FORMAT(3F 10.3) 254.000 C 255.000 0007 RETURN 256.000 0008 END 257.000 *OPTIONS IN EFFECT* ID.EBCDIC,SOURCE,NOLIST.NODECKLOAD.NOMAP *OPTIONS IN EFFECT* NAME = SUROUT, LINECNT 57 *STATISTICS* SOURCE STATEMENTS = 8,PROGRAM SIZE ~ 332 *STATISTICS* NO DIAGNOSTICS GENERATED

C 258.000 C 259.000 C 260.000 C —----------------------------------—..... ---- ----- 261.000 C 262.000 0001 SUBROUTINE TMAT ( CUB, DTA, TM ) 263.000 C 264.000 C calibration of T matrix by least square method 265.000 C 266.000 C Here T matrix is normalized by T44 267.000 C NOTE T matrix is used only In (3-49) of Rogers in Eye. 268.000 C But In Phantom, T matrix should not be mormalized. 269.000 C 270.000 C —--------------—, —. —---------—. 271.000 C. 272.000 C Include global and 10 common here 273.000 C "GLOBAL" common (100,120) 100.000 C 101.000 0002 INTEGER SROW, SCOL 102.000 C Dimensions of the shutter grid 103.000 C 104.000 0003 INTEGER CROW, CCOL 105.000 C Dimensions of the camera grid 106.000 C 107.000 0004 INTEGER NIMAGE, NUMCAL 108.000 C Number of patterns needed 109.000 C Number of calibration data 109.500 C 110.000 0005 COMMON/GLOBAL/ 111.000 1 SROW, SCOL, CROW, CCOL, NIMAGE, NUMCAL 112.000 C 113.000 C "10" common (200, 210 ) 200.000 0006 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 201.000 C Logical units for the calibration file, data, output, 202.000 C and the terminal. 203.000 C 204.000 0007 COMMON/10/ 205.000 1 LUCAL, LUDAT. LUOUT, LUTERM 206.000 C 207.000 C 208.000 C 275.000 C Solve BX * C 276.000 C 277.000 0008 REAL CUB(100,3),TM(4,4) 278.000 0009 INTEGER DTA(iOO,2), IWK(11),IER 279.000 0010 REAL A(200,12),B(200,11),C(200),X(li),COMMU(4),WK(3000) 280.000 C 280.100 C More least square parms 280.200 0011 REAL EPS/4.768E-07/,QL, QLMAX/1.E07/ 280.400 0012 INTEGER ITS/10/. PRT1, PRT2 280.600 C 281.000 C begin 282.000 C 283.000 0013 CALL CLRTBL(A,200,12) 284.000 0014 CALL CLRTBL(TM,4,4) 285.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) TMAT 10-10-83 11:37:49 PAGE P002 C 286.000 0015 NDATA = 2 * NUMCAL 287.000 0016 J=O 288.000 0017 DO 1005 1=1, NDATA, 2 289.000 0018 d=J+1 290.000 0019 A(I,1) a CUB(J,1) 291.000 0020 A(1+1,2) = CUB(J,1) 292.000 0021 A(I,4) - CUB(J.2) 293.000 0022 A(I+1,5) * CUB(J,2) 294.000 0023 A(I,7) - CUB(J,3) 295.000 0024 A(I+1,8) = CUB(J,3) 296.000 0025 A(I,10) = 1. 297.000 0026 A(I+1,11) = 1. 298.000 C 299.000 0027 A(I,12) = - FLOAT( DTA(J,1) 300.000 0028 A(I+1,12) " - FLOAT( DTA(J,2) ) 301.000 0029 A(1,3) a - CUB(J.I) * FLOAT( DTA(J,1) ) 302.000 0030 A(I+1,3L ~ - CUB(J,1) * FLOAT( DTA(J,2) ) 303.000 0031 A(I,6) w - CUB(J,2) * FLOAT( DTA(J,I) ) 304.000 0032 A(I+1,6) ~ - CUB(J.2) * FLOAT( OTA(J,2) ) 305.000 0033 A(I.,9) * - CUB(J,3) * FLOAT( DTA(J.1) ) 306.000 0034 A(I+1,9)' - CUB(J,3) * FLOAT( DTA(J,2) ) 307.000 0035 1005 CONTINUE 308.000 C 309.000 0036. DO 1006 I1, 200 310.000 0037 DO 1007 J=1,11 311.000 0038 B(I,J) a A(I,J) 312.000 0039 1007 CONTINUE 313.000 0040 C(I) = - A(I,12) 314.000 0041 1006 CONTINUE 315.000 C 316.000 C%%%% This Is a call to IMSL package 317.000 C%%%3 318.000 C%%%% CALL LLBOF ( B,200,NOATA,11,C,200, 1.0COMMU.X,11.IWK,WK.IER ) 319.000 C%%%% 320.000 C%%%% IF (IER.EQ. O ) GO TO 1010 320.100 C 320. 150 C This Is a call to L2 320.200 C 320.300 0042 PRT1 a 1 320.360 0043 PRT2 = 0 320.420 0044 CALL L2(NDATA, 1 B.,C,X, EPS ITS,PRTI.PRT2,200,QL) 320.500 C 320.700 0045 IF ( L.LT. OLMAX ) GO TO 1010 321.000 C 322.000 0046 WRITE(LUTERM,2001) OL 323.000 0047.STOP 324.000 C 325.000 0048 1010 CONTINUE 326.000 0049 TM(I.1) * X(1) 327.000 0050 TM(1.2) a X(2) 328.000 0051 TM(1,4) = X(3) 329.000 0052 TM(2.1) a X(4) 330.000 0053 TM(2.2) - X(5) n

0055 TM(3,1) = X(7) 333.000 0056 TM(3,2) = X(8) 334.000 0057 TM(3,4) = X(9) 335.000 0058 TM(4,1) = X(10) 336.000 0059 TM(4.2) = X(11) 337.000 0060 TM(4,4) = 1. 338.000 C 339.000 0061 WRITE(LUTERM,2000) ((TM(I,J),J1=,4),I=1,4) 340.000 C 341.000 0062 RETURN 342.000 C 343.000 0063 2000 FORMAT(/,1X,'T Matrix determined',/,4(4F14.7,/)) 344.000 C 345.000 0064 2001 FORMAT(IX,'Abnormal termination — MESSAGE FROM SUBROUTINE TMAT',/, 346.000 1 1X,'CONDITION NUMBER OF T MATRIX =', E20.5 ) 347.000 C 348.000 0065 END 349.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD.NOMAP *OPTIONS IN EFFECT* NAME = TMAT, LINECNT 57 *STATISTICS* SOURCE STATEMENTS = 65,PROGRAM SIZE = 33300 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:37:50 PAGE P001 C 350.000 C —--- - ------------------ ---------------- --------- 351.000 C 352.000 0001 SUBROUTINE LMAT ( SURPTS, COLPTS, LM ) 353.000 C 354.000 C calibration of L matrix by least square method 355.000 C 356.000 C Here L matrix is normalized by L44 357.000 C NOTE L matrix is used only In (3-49) of Rogers In Eye. 358.000 C But In Phantom, L matrix should not be mormalized. 359.000 C 360.000 C —-------------- ---------- ---------- -----------— 361t.000 C 362.000 C Include global and J0 common here 363.000 C "GLOBAL" common (100,120) 100.000 C 101.000 0002 INTEGER SROW, SCOL 102.000 C Dimensions of the shutter grid 103.000 C 104.000 0003 INTEGER CROW., CCOL 105.000 C Dimensions of the camera grid 106.000 C 107.000 0004 INTEGER NIMAGE, NUMCAL 108.000 C Number of patterns needed 109.000 C Number of calibration data 109.500 C 110.000 0005 COMMON/GLOBAL/ 111.000 1 SROW, SCOL, CROW, CCOL, NIMAGE, NUMCAL 112.000 C 113.000 C "10" common (200, 210 ) 200.000 0006 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 201.000 C Logical units fob the calibration file, data, output, 202.000 C and the terminal. 203.000 C 204.000 0007 COMMON/IO/ 205.000 1 LUCAL, LUDAT, LUOUT. LUTERM 206.000 C 207.000 C 208.000 C 365.000 C Solve BX = C 366.000 C 367.000 0008 REAL SURPTS(100,3),LM(4,.4) 368.000 0009 INTEGER COLPTS(I00), IWK(11), IER 369.000 0010 REAL A(100,8),B(100,7).C(100),.X(7),COMMU(4),WK(3000) 370.000 C 370.100 C More least square parms 370.200 0011 REAL EPS/4.768E-07/.,QL, QLMAX/I1.E07/,WKAREA(200) 370.300 0012 INTEGER ITS/10/, PRT1, PRT2 370.400 C 370.500 C 371.000 C begin 372.000 C 373.000 0013 CALL CLRTBL(A,100,8) 374.000 0014~ CALL CLRTBL(LM,4,4) 375000

0015 DO 1005 1=1, NUMCAL 377.000 0016 A(I,1) = SURPTS(I,1) 378.000 0017 A(l,2) = - SURPTS(I,1) * FLOAT( COLPTS(I) ) 379.000 0018 A(I,3) = SURPTS(I.2) 380.000 0019 A(I,4) = - SURPTS(I,2) * FLOAT( COLPTS(I) ) 381.000 0020 A(I,5) = SURPTS(I,3) 382.000 0021 A(I,6) = - SURPTS(1,3) * FLOAT( COLPTS(I) ) 383.000 0022 A(1,7) = 1. 384.000 0023 A(I,8) = - FLOAT ( COLPTS(I) ) 385.000 0024 1005 CONTINUE 386.000 C 387.000 0925 DO 1006 1=1,100 388.000 0026 DO 1007 d-1,7 389.000 0027 B(I.J) = A(I,J) 390.000 0028 1007 CONTINUE 391.000 0029 C(I) = - A(I,8) 392.000 0030 1006 CONTINUE 393.000 C 394.000 C%%%% This is a call to IMSL package 395.000 C%%%% 396.000 C%%%% CALL LLBQF ( B,100,NUMCAL,7,C,100,1,0, COMMU,X,7,IWK,WK,IER ) 397.000 C%%%% 398.000 C%%%% IF( IER.EQ. 0 ) GO TO 1010 399.000 C 400.000 C This is a call to L2 400.100 C 400.150 0031 PRTi = I 400.200 0032 PRT2 = 0 400.300 0033 CALL L2(NUMCAL,7,B,C,X,EPS,ITS,PRT1,PRT2.100,QL) 400.400 C 400.500 0034 IF ( QL.LT. QLMAX ) GO TO 1010 400.600 C 400.700 0035 WRITE(LUTERM,2001) QL 401.000 0036 STOP 402.000 C 403.000 0037 1010 CONTINUE 404.000 0038 LM(1,1) = X(1) 405.000 0039 LM(1,4) = X(2) 406.000 0040 LM(2,1) = X(3) 407.000 0041 LM(2,4) = X(4) 408.000 0042 LM(3,1) = X(5) 409.000 0043 LM(3,4) = X(6) 410.000 0044 LM(4,1) X(7) 411.000 0045 LM(4,4) = 1. 412.000 C 413.000 0046 WRITE(LUTERM,2000) ((LM(I,J),.J1,4),1N1,4) 414.000 C 415.000 0047 RETURN 416.000 C 417.000 0048 2000 FORMAT(/,IX,'L Matrix determined'./,4(4F14.7,/)) 418.000 C 419.000 0049 2001 FORMAT( IX,'Abnormal termination — MESSAGE FROM SUBROUTINE LMAT', 420.000 1 /, IX,'CONDITION NUMBER OF L MATRIX "', E20.5 ) 421.000 C 422.000 0050 END 423.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) LMAT 10-10-83 11:37:50 PAGE E003 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE.NOLISTNODECK,LOAD.NOMAP *OPTIONS IN EFFECT* NAME - LMAT. LINECNT * 57 *STATISTICS* SOURCE STATEMENTS a 50.PROGRAM SIZE a 20628 *STATISTICS* NO DIAGNOSTICS GENERATED

C 424.000 C 425.000 C*** f*** ***************** S********* ******************************* 426.000 C 427.000 C MAP Module 428.000 C 429.000 C This module does the mapping of the unknown surface. 430.000 C 431.000 C Procedures: 432.000 C Map - Drives module 433.000 C Column - Determines the shutter column, u, for a given 434.000 C Camera grid point. 435.000 C Setcoef - Sets the coefficients for least sqares equations 436.000 C 437.000 C* *** ********************************* 438.000 C 439.000 C —------------------------------------- 440.000 C 441.000 0001 SUBROUTINE MAP (TM,LM) 442.000 C 443.000,C This is the driver routine for the module. 444.000 C For each point Imaged, it determines the 445.000 C shut'ter column which created the Image, 446.000 C sets the coefficients for the equations 447.000 C which determine the surface point, and does 448.000 C the least squares procedure which determines 449.000 C the surface point. 450.000 C 451.000 C Parameters: 452.000 C TM - Transformation matrix for camera geometry 453.000 C LM - Equation of the plane for each column of the shutter 454.000 C 455.000 C —------------------------------------- 456.000 C 457.000 C Include all commons here 458.000 C "GLOBAL" common (100,120) 100.000 C 101.000 0002 INTEGER SROW. SCOL 102.000 C Dimensions of the shutter grid 103.000 C 104.000 0003 INTEGER CROW, CCOL 105.000 C Dimensions of the camera grid 106.000 C 107.000 0004 INTEGER NIMAGE, NUMCAL 108.000 C Number of patterns needed 109.000 C Number of calibration data 109.500 C 110.000 0005 COMMON/GLOBAL/ 111.000 1 SROW, SCOL, CROW, CCOL, NIMAGE, NUMCAL 112.000 C 113.000 C "IO" common (200, 210 ) 200.000 0006 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 201.000 C Logical units for the calibration file, data, output, 202.000 C and the terminal. 203.000 C 204.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAP 10-10-83 11:37:51 PAGE P002 0007 COMMON/O10/ 205.000 1 LUCAL. LUDAT, LUOUT, LUTERM 206.000 C 207.000 C 208.000 C "Patterns" common (300,310) 300.000 0008 INTEGER PAT1(128,128). PAT2(128,128), PAT3(128,128), 301.000 1 PAT4(128,128). PAT5(128,128). PAT6(128,.128), 302.000 2 PAT7( 128.128), PATB( 128,128) 303.000 C Contain the Image data. 304.000 C 305.000 0009 COMMON/PATS/ 306.000 1 PATi1, PAT2. PAT3, PAT4, PATS, PAT6, 307.000 2 PAT7. PATS 308.000 C 460.000 0010 REAL TM(4,4), LM(4.4) 461.000 C 462.000 C Local variables: 463.000 C 464.000 C Matrice's for least squares AX=B 465.000 0011 REAL LSA(3.3), LSB(3), LSX(3) 466.000 C 467.000 C Known values In least squares equations 468.000 0012 REAL COL, XSTAR, YSTAR 469.000 0013 INTEGER U 470.000 C 471.000 C More least square parms 472.000 0014 REAL EPS/4.768E-07/,QL, QLMAX/1.E07/,WKAREA(200) 473.000 0015 INTEGER ITS/10/, PRT1, PRT2 474.000 C 475.000 C Begin 476.000 C 477.000 0016 DO 20 J, 1. CCOL 478.000 C 479.000 0017 DO 10 I = 1, CROW 480.000 C 481.000 0018 IF (PAT1(I,J).EQ. 0) GO TO 10 482.000 C 483.000 0019 CALL CLRVC (LSB.3) 484.000 0020 CALL CLRVC (LSX.3) 485.000 0021 CALL CLRTBL (LSA.3,3) 486.000 C 487.000 0022 CALL COLYMN(I.JU) 488.000 C 489.000 0023 XSTAR. FLOAT(I) 490.000 0024 YSTAR a FLOAT(J) 491.000 0025 COL ~ FLOAT(U) 492.000 C 493.000 0026'CALL SETCOF(LSALSB.XSTAR.YSTARCOL,TM.LM) 494.000 C 495.000 0027 PRT1 = 0 496.000 0028 PRT2 = 0 497.000 0029 CALL L2(3,3,LSA.LSBLSX,.EPS, ITS.PRT1,PRT2,3.QL) 498.000 C 499.000 C%%%%% CALL LEQT2F(LSA.1.3,3. LSB.O.WKAREAQL) 500.000

0030 IF (QL.GT. QLMAX) GO TO 10 502.000 C 503.000 C%%%%% CALL SUROUT( LSB ) 504.000 C 505.000 0031 WRITE (LUOUT, 100) U. I, J, LSX(1), LSX(2), LSX(3) 506.000 I C 507.000 0032 10 CONTINUE 508.000 0033 20 CONTINUE 509.000 C 510.000 0034 100 FORMAT (IX, 316, 3F14.7 ) 511.000 C 512.000 0035 RETURN 513.000 0036 END 514.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOADNOMAP *OPTIONS IN EFFECT* NAME = MAP, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 36,PROGRAM SIZE ~ 1832 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:37:52 PAGE P001 C 515.000 C —--------------- ------- ----------------- -----------— 516.000 C 517.00 0001 SUBROUTINE COLUMN (RWCL.U) 518.000 C 519.000 C Column finds the shutter column (u) for the given point by looking 520.000 C at the Images the point appeared In. It is assumed that pattern two 521.000 C illuminates the right half of the shutter, etc.. 522.000 C 523.000 C Input: RW - camera row for point. 524.000 C CL - camera column for point. 525.000 C Output: U - shutter column number for point 526.000 C 527.000 C —------ ---------------------------------------- 580528.000 C 529.000 C Include GLOBAL and PATTERN commons here 530.000 C "GLOBAL" common (100.120) 100.000 C 101.000 0002 INTEGER SROW, SCOL 102.000 C Dimensions of the shutter grid 103.000 C 104.000 0003 INTEGER CROW, CCOL 105.000 C Dimensions of the camera grid 106.000 C 107.000 0004 INTEGER NIMAGE, NUMCAL 108.000 C Number of patterns needed 109.000 C Number of calibration data 109.500 C 110.000 0005 COMMON/GLOBAL/ I11.000 1 SROW, SCOL, CROW, CCOL, NIMAGE, NUMCAL 112.000 C 113.000 C "Patterns" common (300,310) 300.000 0006 INTEGER PAT1(128,128), PAT2(128,128), PAT3(128,128), 301.000 1 PAT4(128,128), PAT5(128,128), PAT6(128,128), 302.000 2 PAT7(128,128), PAT8(128,128) 303.000 C Contain the image data. 304.000 C 305.000 0007 COMMON/PATS/ 306.000 1 PAT1, PAT2, PAT3, PAT4, PAT5, PAT6, 307.000 2 PAT7, PAT8 308.000 C 533.000 0008 INTEGER RWCL,U 534.000 C 535.000 C begin 536.000 C 537.000 0009 U = PAT2(RW.CL) * ( 2 ** (NIMAGE-2)) 538.000 1 + PAT3(RW.CL) * ( 2 ** (NIMAGE-3)) 539.000 2 + PAT4(RW.CL) * ( 2 ** (NIMAGE-4)) 540.000 3 + PAT5(RW.CL) * ( 2 ** (NIMAGE-5)) 541.000 4 + PAT6(RW,CL) * ( 2 ** (NIMAGE-6)) 542.000 5 + PAT7(RW,CL) * ( 2 ** (NIMAGE-7)) 543.000 6 + PAT8(RW,CL) * ( 2 ** (NIMAGE-8)) 544.000 7 + 1 545.000 0010 RETURN 546.000

0011 END 548.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME w COLUMN. LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 11,PROGRAM SIZE a 986 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:37:52 PAGE POO1 C 549.000 C —------------------------------- ------ ---------- 550.000 C 551.000 0001 SUBROUTINE SETCOF(A. B, XS, VS, U, T, L) 552.000 C 553.000 C Setcoef sets the coefficients on the following three 554.000 C equations which are used to solve for the point (x,y,z): 555.000 C 556.000 C (T11-T14x*)x + (T21-T24x*)y + (T31-T34x*)z + (T41-T44x*) a 0 557.000 C (T12-T14y*)x + (T22-T24y*)y + (T32-T34y*)z + (T42-T44y*) - 0 558.000 C (L11-L14u )x + (L21-L24u )y + (L31-L34u )z + (L41-L44u ) a 0 559.000 C 560.000 C Parameters: 561.000 C 562.000 C Input: XS - x* in the equations 563.000 C YS - y* in the equations 564.000 C U - column of the shutter 565.000 C T - Transformation matrix for camera geometry 566.000 C L - Equation of the plane for each shutter column 567.000 C 568.000 C Output: A, B: Least square matrices 569.000 C 570.000 C —-------------------------------- 571.000 C 572.000 0002 REAL XS, VS. U. T(4.4), L(4.4), A(3,3), B(3) 573.000 C 574.000 C Begin... 575.000 0003 A(1,1) - T(1.1) - T(1,4) * XS 576.000 0004 A(1,2) = T(2,1) - T(2.4) * XS 577.000 0005 A(1,3) - T(3,1) - T(3,4) * XS 578.000 0006 B(1) a - ( T(4,1) - T(4,4) * XS ) 579.000 C 580.000 0007 A(2,1) = T(1,2) - T(1,4) * YS 581.000 0008 A(2,2) = T(2,2) - T(2,4) * YS 582.000 0009 A(2,3) = T(3,2) - T(3.4) * YS 583.000 0010 B(2) = - ( T(4,2) - T(4,4) * YS ) 584.000 C 585.000 0011 A(3, 1) - L(1,1) - L(1,4) * U 586.000 0012 A(3,2) = L(2,1) - L(2,4) * U 587.000 0013 A(3,3) L(3.1) - L(3.4) * U 588.000 0014 8(3) = - ( L(4.1) - L(4,4) * U ) 589.000 C 590.000 C 591.000 0015 RETURN 592.000 0016 END 593.000 *OPTIONS IN EFFECT* ID.EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP tOPTIONS IN EFFECT* NAME = SETCOF, LINECNT = 57 SSTATISTICS* SOURCE STATEMENTS = 16,PROGRAM SIZE ~ 614 *STATISTICS* NO DIAGNOSTICS GENERATED NO STATEMENTS FLAGGED IN THE ABOVE COMPILATIONS.

C 1.000 C Mapping Program: Eye.TLmat 2.000 C 3.000 C 4.000 C 5.000 C Version:This is based on the program of Dr.. Martin AltschulerC 6.000 C Calibration of T matrix and L matrix are not 7.000 C performed. 8.000 C 9.000 C Instead, T and L matrix are calculated by geometric 10.000 C consideration. 11.000 C 12.000 C 13.000 C 14.000 C 15.000 C DIRECTORY 16.000 C SROW # OF ROWS OF CONTOUR GRID 17.000 C SCOL # OF COLUMNS OF CONTOUR GRID 18.000 C CROW # OF ROWS IN ENTIRE RASTER SCREEN 19.000 C CCOL # OF COLUMNS IN ENTIRE RASTER SCREEN 20.000 C NIMAGE # OF INAGES NEEDED TO CODE SCENE 21.000 C 22.000 C LUCAL LOGICAL UNIT FOR CALIBRATION INFO 23.000 C LUDAT LOGICAL UNIT FOR INPUT DATA 24.000 C LUOUT LOGICAL UNIT FOR FILE OUTPUT 25.000 C LUTERM LOGICAL UNIT FOR TERMINAL OUTPUT 26.000 C 27.000 C PATi MATRIX OF RASTER IMAGE # 1 28.000 C: 29.000 C: 30.000 C PAT7 MATRIX OF RASTER IMAGE # 7 31.000 C 32.000 C 33.000 C LOGICAL DEVICES 34.000 C 35.000 C 36.000 C LOGICAL UNIT #3 CALIBRATION INFO 37.000 C LOGICAL UNIT #4 INPUT FILE 38.000 C LOGICAL UNIT #5 OUTPUT FILE 39.000 C LOGICAL UNIT #6 TERMINAL 40.000 C 41.000 C 42.000 C —, —--------------------- ----------------- 443.000 C 44.000 C Main program 45.000 C 46.000 C —----------------------------------------------- 47.000 C 48.000 C Local variables 49.000 C 50.000 C Transformation matrices 51.000 C 52.000 0001 REAL T(4,4), L(4,4) 53.000 C 54.000 C Begin 55.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 1t:39: 12 PAGE P002 C Acquire and check the Initial parameters. 56.000 0002 CALL INIT 57.000 C 58.000 C Calibrate the system. 59.000 0003 CALL TMAT( T ) 60.000 0004 CALL LMAT( L 61.000 C 6.2.000 C Solve for the unknown surface. 63.000 0005 CALL IMGIN 64.000 0006 CALL MAP( T. L ) 65.000 C 66.000 0007 STOP 67.000 0008 END 68.000 *OPTIONS IN EFFECT* ID.EBCDICSOURCE.NOLIST.NODECK,.LOADNOMAP *OPTIONS IN EFFECT* NAME = MAIN LINECNT a 57 *STATISTICS* SOURCE STATEMENTS a 8,PROGRAM SIZE ~ 484 *STATISTICS* NO DIAGNOSTICS GENERATED

~~~~~~~~~~~~~~C ~~~~69.000 C —------------------------------------------- 70.000 C 71.000 0001 SUBROUTINE INIT 72.000 C 73.000 C ----------------------------------- --- -- 74.000 C 75.000 C include GLOBAL and 10 common 76.000 C 77.000 C 1.000 C "GLOBAL" common (1,15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW. SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST, STHETA, SPHI, SSEP, SWIDTH. SHIGH, SROW, SCOL, 14. 000. 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 16.000 C "IO" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/IO/ 23.000 1 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 79.000 C 80.000 C Begin initialization 81.000 C 82.000 0012 LUCAL = 3 83.000 0013 LUDAT = 4 84.000 0014 LUOUT a 5 85.000 0015 LUTERM = 6 86.000 C 87.000 C read global variables 88.000 C. 89.000 0016 READ (LUCAL, 2000) SDIST, STHETA, SPHI, SSEP, SWIDTH, SHIGH, 90.000 1 SCOL, SROW 91.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) INIT 10-10-83 11:39:12 PAGE P002 0017 READ (LUCAL, 2000) CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, 92.000 1 CCOL, CROW 93.000 0018 2000 FORMAT ( F9.4 5F8.4, 216 ) 94.000 C 95.000 0019 NIMAGE = INT( 1.001 + ALOG( FLOAT( SCOL ) ) / ALOG(2.) ) 96.000 C 97.000 C convert to radian 98.000 C 99.000 0020 ANGRAD = 3:141592 / 180 100.000 0021 CTHETA x CTHETA * ANGRAD 101.000 0022 CPHI a CPHI * ANGRAD 102.000 C 103.000 C 104.000 0023 RETURN 105.000 0024 END 106.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE.NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME a INIT, LINECNT 5 57 *STATISTICS* SOURCE STATEMENTS = 24,PROGRAM SIZE = 692 *STATISTICS* NO DIAGNOSTICS GENERATED

C 107.000 C 108.000 C -------- ----------------- 109.000 C 110.000 0001 SUBROUTINE IMGIN 111.000 C 112.000 C t113.000 C — ---------- ------------------- - - - ---.113.000 C 114.000 C INPUT ROUTINE FOR RASTER DATA 115.000 C 116.000 C COMMONS Included here 117.000 C 1.000 C "GLOBAL" common (1,15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW, SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST, STHETA, SPHI, SSEP, SWIDTH, SHIGH, SROW, SCOL. 14.000 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH. CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 16.000 C "IO" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/10/ 23.000 1 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 26.000 C "Patterns" common (26,37) 27.000 C 28.000 0012 INTEGER PAT1(128,128), PAT2(128,128), PAT3(128,128), 29.000 1 PAT4(128,128), PAT5(128,128), PAT6(128,128), 30.000 2 PAT7(128.128), PAT8(128,128) 31.000 C Contain the Image data. 32.000 C 33.000 0013 COMMON/PATS/ 34.000 1 PAT1, PAT2, PAT3, PAT4, PATS, PAT6, 35.000 2 PAT7, PAT8 36.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) IMGIN 10-10-83 11:39:13 PAGE P002 C 119.000 C 120.000 0014 CALL CLRTBL( PATI. 128. 128 121.000 0015 CALL CLRTBL( PAT2. 128. 128 122.000 0016 CALL CLRTBL( PAT3, 128, 128 123.000 0017 CALL CLRTBL( PAT4, 128. 128 124.000 0018 CALL CLRTBL( PAT5, 128. 128 125.000 0019 CALL CLRTBL( PAT6. 128, 128 126.000 0020 CALL CLRTBL( PAT7, 128, 128 127.000 0021 CALL CLRTBL( PAT8, 128, 128 128.000 C 129.000 C LOOP THRU DATA UNTIL YOU GET A (-1,-I) ENTRY 130.000 C TO SIGNAL END OF PATTERN 131.000 C 132.000 0022 1 READ (LUDAT. 2000, END=100) I.J 133.000 0023 IF ((I.EQ.-I).AND. (J.EQ.-O) GO TO 2 134.000 0024 PATI(IJ) = 1 135.000 0025 GO TO 1 136.000 C 137.000 0026 2 READ (LUDAT, 2000,END"100) I.J 138.000 0027 IF ((I.EQ. -I) AND. (J.EQ.-1)) GO TO 3 139.000 0028 PAT2(IJ) 1 140.000 0029 GO. TO 2 141.000 C 142.000 0030 3 READ (LUDAT. 2000, END0100) IJ 143.000 0031 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 4 144.000 0032 PAT3(I.J) = 1 145.000 0033 GO TO 3 146.000 C 147.000 0034 4 READ (LUDAT, 2000, ENDI100) I.J 148.000 0035 IF ((I.EQ. -1).AND. (J.EQ.-I)) GO TO 5 149.000 0036 PAT4(I.J) = 1 150.000 0037 GO TO 4 151.000 C 152.000 0038 5 READ (LUDAT, 2000, END100) 1. 153.000 0039 IF ((I EQ. -1) AND. (J.EQ.-1)) GO TO 6 154.000 0040 PAT5(IJ) = 1 155.000 0041 GO 10 5 156.000 C 157.000 0042 6 READ (LUDAT. 2000, END=100) I. J 158.000 0043 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 7 159.000 0044 PAT6(1.J) - 1 160.000 0045 GO TO 6 161.000 C 162.000 0046 7 READ (LUDAT, 2000, END~100) I.J 163.000 0047 IF ((I EQ. -1) AND. (J.EQ.-0) GO TO 8 164.000 0048 PAT7(I.J) = 1 165.000 0049 GO TO 7 166.000 C 167.000 0050 8 READ (LUDAT, 2000. END~100) IJ 168.000 0051 IF ((I EQ. -1) AND. (tJ.EQ.-1)) GO TO 100 169.000 0052 PAT8(I.J) - 1 170.000 0053 GO TO 8 171.000

I'!!t11 II I CKMIItL 311ICM rUcm rn ui " I.ol 1M U1 I)U- IU-OJ l I:u:3 I J rAU rVVUU C 174.000 0054 100 RETURN 175.000 0055 2000 FORMAT( 216 ) 176.000 0056 3000 FORMAT( 1614 ) 177.000 C 178.000 0057 END 179. 000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = IMGIN, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 57,PROGRAM SIZE 1710 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:39:13 PAGE P001 C 180.000 C 181.000 C —-----------------------— ~ —--- ------ 182.000 0001 SUBROUTINE SUROUT( SURPNT ) 184.000 C 186.000 C OUTPUT ROUTINE FOR ESTIMATED LOCATION 187.000 C OF POINTS 188.000 C —------------------------------------------------------- 189.000 C 190.000 C Commons Included here 191.000 C 16.000 C "10D" COMMON (16,25) 17.000 C 18.000 0002 INTEGER LUCAL, LUDAT, LUOUT. LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0003 COMMON/10/ 23.000 1 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 193.000 C 194.000 C PARAMETERS 195.000 0004 REAL SURPNT(3) 196.000 C 197.000 0005 WRITE (LUOUT, 2000) SURPNT(1), SURPNT(2), SURPNT(3) 198.000 0006 2000 FORMAT(3F 10.3) 199.000 C 200.000 0007 RETURN 201.000 0008 END 202.000 *OPTIONS IN EFFECT* ID,EBCDICSOURCE,NOLIST.,NODECK.LOAD.NOMAP *OPTIONS IN EFFECT* NAME - SUROUT, LINECNT a 57 *STATISTICS* SOURCE STATEMENTS 8,PROGRAM SIZE 332 *STATISTICS* NO DIAGNOSTICS GENERATED

C 203.000 C 204.000 C 205.000 C ------------------------------------------------------- 206.000 C 207.000 0001 SUBROUTINE TMAT ( TM ) 208.000 C 209.000 C' read in T matrix from phantom 210.000 C —------------------------------------ 211.000 C 212.000 C Include IO common here 213.000 C 16.000 C "IO" COMMON (16,25) 17.000 C 18.000 0002 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0003 COMMON/I0/ 23.000 1 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 215.000 C 216.000 0004 REAL TM(4, 4) 217.000 C 218.000 0005 CALL CLRTBL(TM, 4, 4) 219.000 C 220.000 0006 READ (LUCAL,2000) (( TM(I, J), J = 1, 4), I = 1, 4) 221.000 C WRITE (LUTERM,3000) (( I, J, TM(I, J), dJ 1, 4), I - 1, 4) 222.000 C 223.000 0007 2000 FORMAT( F15.7, 3F14.7 ) 224.000 C3000 FORMAT( 4( 214, F14.7 )) 225.000 C 226.000 0008 RETURN 227.000 0009 END 228.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = TMAT, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 9,PROGRAM SIZE ~ 460 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:39:14 PAGE POO0 C 229.000 C 230.000 C --------------------------------------------------------- 3231.000 C 232.000 0001 SUBROUTINE LMAT ( LM ) 233.000 C 234.000 C read In L matrix from phantom 235.000 C —------------------------------------------------------ 236.000 C 237.000 C Include GLOBAL and 10 Common here 238.000 C 1.000 C "GLOBAL" common (,1.5) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW, SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6. 600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST. STHETA. SPHI, SSEP. SWIDTH. SHIGH, SROW. SCOL. 14.000 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE' 15.000 C 16.000 C "IO" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/IO/ 23.000 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 240.000 C 241.000 0012 REAL LM(4,4) 242.000 C 243.000 0013 CALL CLRTBL(LM. 4, 4) 244.000 C 245.000 0014 READ (LUCAL, 2000) ( ( LM(I, J), J l1. 4 ), I " 1, 4 ) 246.000 C WRITE (LUTERM, 3000) (( I, J. LM(I, J), J ~ 1, 4 ), I = 1, 4 ) 247.000 C 1 248.000 0015 2000 FORMAT( F15.7, 3F14.7 ) 249.000 C3000 FORMAT( 4( 214, F14.7 ) ) 250.000

0016 RETURN 252.000 0017 END 253.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLISTNODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = LMAT, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 17,PROGRAM SIZE = 468 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 11:39:14 PAGE POO1 C 254.000 ---------------------------------------------------- 255.000 C 256.000 C MAP Module 257.000 c 258.000 C This module does the mapping of the unknown surface. 259.000 C 260.000 C Procedures: 261.000 C Map - Drives module 262.000 C Column - Determines the shutter column, u, for a given 263.000 C Camera grid point. 264.000 C Setcoef - Sets the coefficients for least sqares equations 265.000 C 266.000 C —------------------------------------- 267.000 C- -------------------------------------------------------------- 269.000 C 270.000 0001 SUBROUTINE MAP (TM,LM) 271.000 C 272.000 C This Is the driver routine for the module. 273.000 C For each point Imaged, It determines the 274.000 C shutter column which created the Image, 275.000 C sets the coefficients for the equations 276.000 C which determine the surface point, and does 277.000 C the least squares procedure which determines 278.000 C the surface point. 279.000 C 280.000 C Parameters: 281.000 C TM - Transformation matrix for camera geometry 282.000 C LM - Equation of the plane for each column of the shutter 283.000 C 284.000 C —-------------------------------------------------------- 285.000 C 286.000 C Include all commons here 287.000 C 1.000 C "GLOBAL" common (1,15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW, SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP. CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 rnfMun /ni nn A /

2 CDIST, CTHETA. CPHI, CSEP, CWIDTH, CHIGH, CROW, dCOL, 14.500 3 NIMAGE 15.000 C 16.000 C "1O" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT. LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/10O/ 23.000 1 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 26.000 C "Patterns" common (26,37) 27.000 C 28.000 0012 INTEGER PAT1(128,.128), PAT2(128,128), PAT3(128,128), 29.000 1 PAT4(128.,128), PAT5(128,.128), PAT6(128,128), 30.000 2 PAT7(128,128), PAT8(128,128) 31.000 C Contain the Image data. 32.000 C 33.000 0013 COMMON/PATS/ 34.000 1 PATI, PAT2, PAT3, PAT4, PATS, PAT6, 35.000 2 PAT7, PAT8 36.000 C 289.000 0014 REAL TM(4,4), LM(4,4) 290.000 C 291.000 C Local variables: 292.000 C 293.000 C Matrices for least squares AX=B 294.000 0015 REAL LSA(3,3). LSB(3), LSX(3) 295.000 C 296.000 C Known values In least squares equations 297.000 0016 REAL COL, XSTAR. YSTAR 298.000 0017 INTEGER U 299.000 C 300.000 C More least square parms 301.000 0018 REAL EPS/4.768E-07/,QL, QLMAX/1.E07/,WKAREA(200) 302.000 0019 INTEGER ITS/10/, PRTi, PRT2 303.000 C 304.000 C Begin... 305.000 C 306.000 0020 DO 20 J = 1, CCOL 307.000 C 308.000 0021 DO 10 1 = 1i, CROW 309.000 C 310.000 0022 IF (PATi(I.J).EQ. 0) GO TO 10 311.000 C 312.000 0023 CALL CLRVC (LSB,3) 313.000 0024 CALL CLRVC (LSX,3) 314.000 0025 CALL CLRTBL'(LSA,3,3) 315.000 C 316.000 0026 CALL COLUMN(I.J,U) 317.000 C 318.000 0027 XSTAR = FLOAT(I) 319.000 0028 YSTAR = FLOAT(J) 320.000 0029 COL = FLOAT(U) 321.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAP 10-10-83 11:39:14 PAGE P003 C 322.000 0030 CALL SETCOF(LSA,LSB,XSTAR,YSTAR,COLTM,LM) 323.000 C 324.000 0031 PRT1 = 0 325.000 0032 PRT2 = 0 326.000 0033 CALL L2(3,3,LSALSB,LSX,EPS,ITS,PRT1,PRT2,3,QL) 327.000 C 328.000 C%%%%% CALL LEQT2F(LSA,1,3,3, LSB,O,WKAREA,QL) 329.000 C 330.000 0034. IF (QL.GT. QLMAX) GO TO 10 331.000 q' 332.000 C%%%%% CALL SUROUT( LSB ) 333.000 C 334.000 0035 WRITE (LUOUT, 100) U. 1, J, LSX(1), LSX(2), LSX(3) 336.000 C 337.000 0036 10 CONTINUE 338.000 0037 20 CONTINUE 339.000 C 340.000 0038 100 FORMAT (IX, 316, 3F14.7 ) 342.000 C 343.000 0039 RETURN 344.000 0040 END 345.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE.NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = MAP, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS - 40,PROGRAM SIZE ~ 1832 *STATISTICS* NO DIAGNOSTICS GENERATED

C 346.000 C —----------------------------------------- ---------- 347.000 C 348.000 0001 SUBROUTINE COLUMN (RW,CL,U) 349.000 C 350.000 C Column finds the shutter column (u) for the given point by looking 351.000 C at the Images the point appeared In. It is assumed that pattern two 352.000 C illuminates the right half of the shutter, etc.. 353.000 C 354.000 C Input: fW - camera row for point. 355.000 C CL - camera column for point. 356.000 C Output: U - shutter column number for point 357.000 C 358.000 C —--- ------------------ --------------------------- 359.000 C 360.000 C Include GLOBAL and PATTERN commons here 361.000 C 1.000 C "GLOBAL" common (1,15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW, SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL'CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST, STHETA, SPHI, SSEP, SWIDTH, SHIGH, SROW, SCOL, 14.000 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 26.000 C "Patterns" common (26,37) 27.000 C 28.000 0010 INTEGER PATI(112,128), PAT2(128,128), PAT3(128,128), 29.000 1 PAT4(128,128), PAT5(128,128), PAT6(128,128), 30.000 2 PAT7(128,128), PAT8(128,128) 31.000 C Contain the Image data. 32.000 C 33.000 0011 COMMON/PATS/ 34.000 1 PATI, PAT2, PAT3, PAT4, PATS, PAT6, 35.000 2 PAT7, PAT8 36.000 C 364.000 0012 INTEGER RW,CL,U 365.000 C 366.000 C begin 367.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) COLUMN 10-10-83 11:39:16 PAGE PO02 C 368.000 0013 U = PAT2(RW,CL) * ( 2 ** (NIMAGE-2)) 369.000 1 + PAT3(RW,CL) * ( 2** (NIMAGE-3)) 370.000 2 + PAT4(RWCL) * ( 2 ** (NIMAGE-4)) 371.000 3 + PAT5(RW,CL) * ( 2 ** (NIMAGE-5)) 372.000 4 + PA16(RWCL) * ( 2 ** (NIMAGE-6)) 373.000 5 + PAT7(RWCL) * ( 2 ** (NIMAGE-7)) 374.000 6 + PAT8(RW.CL) * ( 2 ** (NIMAGE-8)) 375.000 7 + 1 376.000 0014 RETURN 377.000 C 378.000 0015 END 379.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK.LOAD,NOMAP *OPTIONS IN EFFECT* NAME a COLUMN, LINECNT a 57 *STATISTICS* SOURCE STATEMENTS ~ 15.PROGRAM SIZE - 986 *STATISTICS* NO DIAGNOSTICS GENERATED

C 380.000 C —------------------------------------------ ------- 381.000 C 382.000 0001 SUBROUTINE SETCOF(A, B, XS, YS, U, T, L) 383.000 C 384.000 C Setcoef sets the coefficients on the following three 385.000 C equations which are used to solve for the point (x,y,z): 386.000 C 387.000 C (T11-T14x*)x + (T21-T24x*)y + (T31-T34x*)z + (T41-T44x*) = 0 388.000 C (T12-T14y*)x + (T22-T24y*)y + (T32-T34y*)z + (T42-T44y*) = 0 389.000 C (L11-Li4u )x + (L21-L24u )y + (L31-L34u )z + (L41-L44u ) = 0 390.000 C 391.000 C Parameters: 392.000 C 393.000 C Input: XS - x* in the equations 394.000 C YS - y* In the equations 395.000 C U - column of the shutter 396.000 C T - Transformation matrix for camera geometry 397.000 C L - Equation of the plane for each shutter column 398.000 C 399.000 C Output: A, B: Least square matrices 400.000 C 401.000 C —---- -------------— 402.000 C 403.000 0002 REAL XS, YS, U, T(4,4), L(4,4), A(3,3), B(3) 404.000 C 405.000 C Begin... 406.000 0003 A(1,1) = T(1,1) - T(1,4) * XS 407.000 0004 A(1,2) = T(2,1) - T(2,4) * XS 408.000 0005 A(1,3) = T(3,1) - T(3,4) * XS 409.000 0006 B(1) = - ( T(4,1) - T(4,4) * XS ) 410.000 C 411.000 0007 A(2,1) = T(1,2) - T(1,4) * YS 412.000 0008 A(2,2) = T(2,2) - T(2,4) * YS 413.000 0009 A(2,3) = T(3,2) - T(3,4) * YS 414.000 0010 B(2) = - ( T(4,2) - T(4,4) * YS ) 415.000 C 416.000 0011 A(3,1) = L(1,1) - L(1,4) * U 417.000 0012 A(3,2) = L(2,1) - L(2,4) * U 418.000 0013 A(3,3) = L(3,1) - L(3,4) * U 419.000 0014 B(3) = - ( L(4,1) - L(4,4) * U ) 420.000 C 421.000 C 422.000 0015 RETURN 423.000 0016 END 424.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = SETCOF, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 16,PROGRAM SIZE = 614 *STATISTICS* NO DIAGNOSTICS GENERATED NO STATEMENTS FLAGGED IN THE ABOVE COMPILATIONS.

28 8.3. Eye.Tlmat Program

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 17:05:51 PAGE P001 C 1.000 C Mapping Program: Eye.TLmat 2.000 C 3.000 C 4.000 C 5.000 C Version: This is based on the program of Dr. Martin AltschulerC 6.000 C Calibration of T matrix and L matrix are not 7.000 C performed. 8.000 C 9.000 C Instead, T and L matrix are calculated by geometric 10.000 C consideration. 11.000 C 12.000 C 13.000 C 14,000 C 15.000 C DIRECTORY 16.000 C SROW # OF ROWS OF CONTOUR GRID 17.000 C SCOL N OF COLUMNS OF CONTOUR GRID 18.000 C CROW N OF ROWS IN ENTIRE RASTER SCREEN 19.000 C CCOL N OF COLUMNS IN ENTIRE RASTER SCREEN 20.000 C NIMAGE N OF INAGES NEEDED TO CODE SCENE 21.000 C I 22.000 C LUCAL LOGICAL UNIT FOR CALIBRATION INFO 23.000 C LUDAT LOGICAL UNIT FOR INPUT DATA 24.000 C LUOUT LOGICAL UNIT FOR FILE OUTPUT 25.000 C LUTERM LOGICAL UNIT FOR TERMINAL OUTPUT 26.000 C 27.000 C PAT1 MATRIX OF RASTER IMAGE N 1 28.000 C 29.000 C 30.000 C PAT7 MATRIX OF RASTER IMAGE N 7 31.000 C 32.000 C 33.000 C LOGICAL DEVICES 34.000 C 39.000 C 36.000 C LOGICAL UNIT #3 CALIBRATION INFO 37.000 C LOGICAL UNIT N4 INPUT FILE 38.000 C LOGICAL UNIT #5 OUTPUT FILE 39.000 C LOGICAL UNIT #6 TERMINAL 40.000 C 41.000 C 42.000 C- --- -- - -—..... —-- - - - - - - - - - -,-. —-—._ —-—..- 43.000 C 44.000 C Main program 45.000 C 46.000 C- -----------------------------------------------—., — 47.000 C 48.000 C Local variables 49.000 C 50.000 C Transformation matrices 51.000 C 52.000 0001 REAL T(4.4). L(4,4) 53.000

C Acquire and check the Initial parameters. 56.000 0002 CALL INIT 57.000 C 58.000 C Calibrate the system. 59.000 0003 CALL TMAT( T ) 60.000 0004 CALL LMAT( L ) 61.000 C 62.000 C Solve for the unknown surface. 63.000 0005 CALL IMGIN 64.000 0006 CALL MAP( T, L ) 65.000 C 66.000 0007 STOP 67.000 0008 END 68.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST.NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = MAIN, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 8,PROGRAM SIZE - 484 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 17:05:51 PAGE P001 C 69.000 C -------—. —...-.-......... ---------... - 70.000 C 71.000 0001 SUBROUTINE INIT 72.000 C 73.000 C —------------------- 74.000 C 75'. 000 C Include GLOBAL and 10 common 76.000 C 77.000 C 1.000 C "GLOBAL" common (1,15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP. SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW, SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 / C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST, STHETA, SPHI, SSEP, SWIDTH, SHIGH. SROW, SCOL, 14.000 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 16.000 C "10" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/IO/ 23.000 LUCAL. LUDAT, LUOUT, LUTERM 24.000 C 79.000 C 80.000 C Begin Initialization 81.000 C 82.000 0012 LUCAL = 3 83.000 0013 LUDAT = 4 84.000 0014 LUOUT = 5 85.000 0015 LUTERM - 6 86.000 C 87.000 C read global variablos 88.000 C 89.000 0016 READ (LUCAL. 2000) SDIST. STHETA, SPHI, SSEP, SWIDTH. SHIGH, 90.000

0017 READ (LUCAL, 2000) CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, 92.000 1 CCOL, CROW 93.000 0018 2000 FORMAT ( F9.4, 5F8.4, 216 ) 94.000 C 95.000 0019 NIMAGE = INT( 1.001 + ALOG( FLOAT( SCOL ) ) / ALOG(2.) ) 96.000 C 97.000 C convert to radian 98.000 C 99.000 0020 ANGRAD = 3. 41592 / 180 100.000 0021 CTHETA = CTHETA * ANGRAD 101.000 0022 CPHI = CPHI * ANGRAD 102.000 C 103.000 C 104.000 0023 RETURN 105.000 0024 END 106.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = INIT, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 24,PROGRAM SIZE = 692 *STATISTICS*'NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 17:05:52 PAGE P001 C 107.000 C 108.000 C ------------------------------------------------------------- 109.000 C 110.000 0001 SUBROUTINE IMGIN 111.000 C 112.000 C --------------- 113.000 C 114.000 C INPUT ROUTINE FOR RASTER DATA 115.000 C 116.000 C COMMONS Included here 117.000 c 1.000 C "GLOBAL" common (1,.15) 2.000 C 3.000 0002 REAL SDIST. STHETA. SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH. SHIGH 3.600 C 33.800 0004 INTEGER SROW. SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST, SIHETA, SPHI, SSEP, SWIDTH, SHIGH, SROW. SCOL, 14.000 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 16.000 C "10" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/IO/ 23.000 1 LUCAL, LUDAT, LUOUT. LUTERM 24.000 C 26.000 C "Patterns" common (26.37) 27.000 C 28.000 0012 INTEGER PATi(128.128), PAT2(128,128), PAT3(128,128), 29.000 ~1 PAT4(128.128). PAT5(128,128), PAT6(128,128). 30.000 2 PAT7(128.128), PAT8(128,128) 31.000 C Contain the Image data. 32.000 C 33.000 0013 COMMON/PATS/ 34.000 1 PATI, PAT2, PAT3, PAT4, PAT5, PAT6, 35.000

c 119.000 C 120.000 0014 CALL CLRTBL( PAT1, 128, 128 ) 121.000 0015 CALL CLRTBL( PAT2. 128, 128 ) 122.000 0016 CALL CLRTBL( PAT3, 128, 128 ) 123.000 00117 CALL CLRTBL( PAT4, 128, 128 ) 124.000 0018 CALL CLRTBL( PATS, 128, 128 ) 125.000 0019 CALL CLRTBL( PAT6, 128, 128 ) 126.000 0020 CALL CLRTBL( PAT7, 128, 128 ) 127.000 0021 CALL CLRTBL( PAT8, 128, 128 ) 128.000 C 129.000 C LOOP THRU DATA UNTIL YOU GET A (-1,-1) ENTRY 130.000 C TO SIGNAL END OF PATTERN 131.000 C 132.000 0022 1 READ (LUDAT, 2000, END=100) I,J 133.000 0023 IF ((I.EQ.-1).AND. (J.EQ.-1)) GO TO 2 134.000 0024 PAT1(I,J) = 1 135.000 0025 GO TO 1 136.000 C 137.000 0026 2 READ (LUDAT, 2000,END=100) I,J 138.000 0027 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 3 139.000 0028 PAT2(I,J) = 1 140.000 0029 GO TO 2 141.000 C 142.000 0030 3 READ (LUDAT, 2000, END=100) I, 143.000 0031 IF ((I.EQ.-1) AND. (J.EQ.-1)) GO TO 4 144.000 0032 PAT3(1,J) = 1 145.000 0033 GO TO 3 146.000 C' 147.000 0034 4 READ (LUDAT, 2000, END=100) I,J 148.000 0035 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 5 149.000 0036 PAT4(1,J) 1 150.000 0037 GO TO 4 151.000 C 152.000 0038 5 READ (LUDAT, 2000, END=100) I,J 153.000 0039 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 6 154.000 0040 PAT5(I,J) = 1 155.000 0041 GO TO 5 156.000 C 157.000 0042 6 READ (LUDAT, 2000, END=100) I,J 158.000 0043 IF ((I.EQ. -1).AND. (J.EQ.-t)) GO TO 7 159.000 0044 PAT6(1,J) = 1 160.000 0045.GO TO 6 161.000 C 162.000 0046 7 READ (LUDAT, 2000, END=100) I,J 163.000 0047 IF ((I.EO. -1).AND. (d.EQ.-1)) GO TO 8 164.000 0048 PAT7(I,J) = 1 165.000 0049 GO TO 7 166.000 C 167.000 0050 8 READ (LUDAT, 2000, END 100) I,J 168.000 0051 IF ((I.EQ. -1).AND. (J.EQ.-1)) GO TO 100 169.000 0052 PAT8(I,J) = 1 170.000 0053 GO TO 8 171.000 C 172.000 C100 WRITE( LUTERM, 3000 ) (( PATI(I.J), J=1,CCOL), I=1,CROW ) 173.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) IMGIN 10-10-83 17:05:52 PAGE P003 C 174.000 0054 100 RETURN 175.000 0055 2000 FORMAT( 216 ) 176.000 0056 3000 FORMAT( 1614 ) 177.000 C 178.000 0057 END 179.000 *OPTIONS IN EFFECT* ID.EBCDIC,SOURCE,NOLISTNODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = IMGIN, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 57,PROGRAM SIZE ~ 1710 *STATISTICS* NO DIAGNOSTICS GENERATED

C 180.000 C 181.000 C, _ _______ —------------------------------- ------------ 182.000 0001 SUBROUTINE SUROUT( SURPNT ) 184.000 C 186.000 C OUTPUT ROUTINE FOR ESTIMATED LOCATION 187.000 C OF POINTS 188.000 C —--------------------------------------------------------- 189.000 C 190.000 C.Commons Included here 191.000 C 16.000 C "I0" COMMON (16,25) 17.000 C 18.000 0002 INTEGER LUCAL, LUDAT, LUOUT, LUTERM. 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0003 COMMON/IO/ 23.000 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 193.000 C 194.000 C PARAMETERS 195.000 0004 REAL SURPNT(3) 196.000 C 197.000 0005 WRITE (LUOUT, 2000) SURPNT(1), SURPNT(2), SURPNT(3) 198.000 0006 2000 FORMAT(3F10.3) 199.000 C 200.000 0007 RETURN 201.000 0008 END 202.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCENOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = SUROUT, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 8,PROGRAM SIZE = 332 *STATISTICS* NO DIAGNOSTICS GENERATED

hutrCGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 17:05:53 PAGE POO1 C 203.000 C 204.000 C 205.000 C — - - - - - - --------------------------------------- 206.000 C 207.000 0001 SUBROUTINE TMAT ( TM ) 208.000 C 209.000 C read in T matrix from phantom 210.000 C-~- - -- - -- - -----— ~- — ~ —. -.~ ~-~-. —--- -. ~ ~..-~ —- - 211.000 C 212.000 C Include 10 common here 213.000 C 16.000 C "IO" COMMON (16,25) 17.000 C 18.000 0002 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical unite for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0003 COMMON/I0/ 23.000 1 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 215.000 C 216.000 0004 REAL TM(4, 4) 217.000 C 218.000 0005 CALL CLRTBL(TM, 4. 4) 219.000 C 220.000 0006 READ (LUCAL,2000) (( TM(I, 0), J x 1, 4). 1 I 1, 4) 221.000 C WRITE (LUTERM,3000) (( 1. J. TM(I, J), J 1, 4), I 1., 4) 222.000 C 223.000 0007 2000 FORMAT( F15.7, 3F14.7 ) 224.000 C3000 FORMAT( 4( 214, F14.7 )) 225.000 C 226.000 0008 RETURN 227.000 0009 END 228.000 *OPTIONS IN EFFECT* ID,.EBCDICSOURCENOLISTNODECKLOAD.NOMAP *OPTIONS IN EFFECT* NAME = TMAT, LINECNT * 57 *STATISTICS* SOURCE STATEMENTS 9.,PROGRAM SIZE ~ 460 *STATISTICS* NO DIAGNOSTICS GENERATED

C 229.000 C 230.000 C —----------- ------------------------- 231.000 C 232.000 0001 SUBROUTINE LMAT ( LM ) 233.000 C 234.000 C read In L matrix from phantom 235.000 C —------------------------------------------- 236.000 C 237.000 C Include GLOBAL and 10 Common here 238.000 C 1.000 C "GLOBAL" common (1.15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW, SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST, STHETA, SPHI, SSEP, SWIDTH, SHIGH, SROW, SCOL, 14.000 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 16.000 C "IO" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/IO/ 23.000 1 LUCAL, LUDAT, LUOUT, LUTERM 24.000 C 240.000 C 241.000 0012 REAL LM(4,4) 242.000 C 243.000 0013 CALL CLRTBL(LM. 4, 4) 244.000 C 245.000 0014 READ (LUCAL, 2000) ( ( LM(I, J), Jd 1, 4 ), I ~ 1, 4 ) 246.000 C WRITE (LUTERM, 3000) (( I, J, LM(I, J), J ~ 1, 4 ), I! 1, 4 ) 247.000 C 248.000 0015 2000 FORMAT( F15.7, 3F14.7 ) 249.000 C3000 FORMAT( 4( 214. F14.7 ) ) 250.000 C 251.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) LMAT 10-10-83 17:05:54 PAGE P002 0016 RETURN 252.000 0017 END 253.000 *OPTIONS IN EFFECT* ID.EBCDICSOUHCENOLISr,NODECK.LOAD.NOMAP *OPTIONS IN EFFECT* NAME * LMAT, LINECNT a 57 *STATISTICS* SOURCE STATEMENTS' 17,PROGRAM SIZE * 468 *STATISTICS* NO DIAGNOSTICS GENERATED

C 254.000 C —------------------------------------- 255.000 C -255.000 C 256.000 C MAP Module 257.000 C. 258.000 C This module does' the mapping of the unknown surface. 259.000 C 260.000 C Procedures: 261.000 C Map - Drives module 262.000 C Column - Determines the shutter column, u, for a given 263.000 C Camera grid point. 264.000 C Setcoef - Sets the coefficients for least sqares equations 265.000 C 266.000 C —------------------------------------- 267.000 C C —-269.000 C 270.000 0001 SUBROUTINE MAP (TM,LM) 271.000 C 272.000 C This Is the driver routine for the module. 273.000 C For each point Imaged, It determines the 274.000 C shutter column which created the Image, 275.000 C sets the coefficients for the equations 276.000 C which determine the surface point, and does 277.000 C the least squares procedure which determines 278.000 C the surface point. 279.000 C 280.000 C Parameters: 281.000 C TM - Transformation matrix for camera geometry 282.000 C LM - Equation of the plane for each column of the shutter 283.000 C 284.000 C —------------------ 285.000 C 286.000 C Include all commons here 287.000 C 1.000 C "GLOBAL" common (1.15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIDTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW, SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW, CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST, STHETA, SPHI, SSEP, SWIDTH, SHIGH, SROW, SCOL, 14.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAP 10-10-83 17:05:54 PAGE P002 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 16.000 C "10" COMMON (16,25) 17.000 C 18.000 0010 INTEGER LUCAL, LUDAT, LUOUT, LUTERM 19.000 C Logical units for the calibration file, data, output, 20.000 C and the terminal. 21.000 C 22.000 0011 COMMON/IO/ 23.000 1 LUCAL, LUDAT, LUOUT. LUTERM 24.000 C 26.000 C "Patterns" common (26.37) 27.000 C 28.000 0012 INTEGER PAT(128, 128), PAT2(128,128), PAT3(128,128). 29.000 1 PAT4(128, 128), PAT5(128,128), PAT6(128,128), 30.000 2 PAT7(128,128), PAT8(128,128) 31.000 C Contain the Image data. 32.000 C 33.000 0013 COMMON/PATS/ 34.000 1 PAT1, PAT2. PAT3, PAT4, PAT5. PAT6, 35.000 2 PAT7. PAT8 36.000 C 289.000 0014 REAL TM(4,4), LM(4,4) 290.000 C 291.000 C Local variables: 292.000 C 293.000 C Matrices for least squares AX=B 294.000 0015 REAL LSA(3,3), LSB(3), LSX(3) 295.000 C 296.000 C Known values In least squares equations 297.000 0016 REAL COL, XSTAR, YSTAR 298.000 0017 INTEGER U 299.000 C 300.000 C More least square parms 301.000 0018 REAL EPS/4.768E-07/,QL, QLMAX/1.E07/,WKAREA(200) 302.000 0019 INTEGER ITS/10/, PRT1, PRT2 303.000 C 304.000 C Begin... 305.000 C 306.000 0020 DO 20 J = 1, CCOL 307.000 C 308.000 0021 DO 10 I = 1, CROW 309.000 C 310.000 0022 IF (PATI(I.J).E. 0) GO TO 10 311.000 C 312.000 0023 CALL CLRVC (LSB,3) 313.000 0024 CALL CLRVC (LSX,3) 314.000 0025 CALL CLRTBL (LSA,3,3) 315.000 C 316.000 0026 CALL COLUMN(I.J.U) 317.000 C 318.000 0027 XSTAR = FLOAT(I) 319.000 0028 VTAR r= FlnA^. ^ i

C 322.000 0030 CALL SETCOF(LSA,LSB,XSTAR,YSTAR,COL,TM,LM) 323.000 C 324.000 0031 PRTI = 0 325.000 0032 PRT2 = 0 326.000 0033 CALL L2(3,3.LSA.LSB.LSX.EPS, ITS,PRT 1PRT2,3,QL) 327.000 C 328.000 C%%%%% CALL LEQT2F(LSA,1,3,3, LSB,O,WKAREA,QL) 329.000 C 330.000 0034 IF (QL.GT. QLMAX) GO TO 10 331.000 C 332.000 C%%%%% CALL SUROUT( LSB ) 333.000 C 334.000 0035 WRITE (LUOUT, 100) U, 1.,J LSX(1), LSX(2), LSX(3) 336.000 C 337.000 0036 10 CONTINUE 338.000 0037 20 CONTINUE 339.000 C 340.000 0038 100 FORMAT (IX, 316, 3F14.7 ) 342.000 C 343.000 0039 RETURN 344.000 0040 END 345.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME - MAP. LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 40,PROGRAM SIZE " 1832 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 17:05:55 PAGE P001 C 346.000 C —----- ------------------- --- 347.000 C 348.000 0001 SUBROUTINE COLUMN (RW,CL,U) 349.000 C 350.000 C Column finds the shutter column (u) for the given point by looking 351.000 C at the images the point appeared in. It is assumed that pattern two 352.000 C illuminates the right half of the shutter. etc.. 353.000 C 354.000 C Input: RW - camera row for point. 355.000 C CL - camera column for point. 356.000 C Output: U - shutter column number for point 357.000 C 358.000 C —---------------------------------------------------------------- 359.000 C 360.000 C Include GLOBAL and PATTERN commons here 361.000 C 1.000 C "GLOBAL" common (1.15) 2.000 C 3.000 0002 REAL SDIST, STHETA, SPHI 3.200 C 3.400 0003 REAL SSEP, SWIOTH, SHIGH 3.600 C 3.800 0004 INTEGER SROW. SCOL 4.000 C Dimensions of the shutter grid 5.000 C 6.000 0005 REAL CDIST, CTHETA, CPHI 6.300 C 6.600 0006 REAL CSEP, CWIDTH, CHIGH 6.700 C 6.800 0007 INTEGER CROW. CCOL 7.000 C Dimensions of the camera grid 8.000 C 9.000 0008 INTEGER NIMAGE 10.000 C Number of patterns needed. 11.000 C 12.000 0009 COMMON/GLOBAL/ 13.000 1 SDIST. STHETA, SPHI, SSEP, SWIDTH, SHIGH, SROW, SCOL, 14.000 2 CDIST, CTHETA, CPHI, CSEP, CWIDTH, CHIGH, CROW, CCOL, 14.500 3 NIMAGE 15.000 C 26.000 C "Patterns" commor (26,37) 27.000 C 28.000 0010 INTEGER PAT1(128,128). PAT2(128,128), PAT3(128,128), 29.000 1 PAT4(128,128), PAT5(128,128). PAT6(128,128), 30.000 2 PAT7(128.128). PAT8(128.128) 31.000 C Contain the image data. 32.000 C 33.000 0011 COMMON/PATS/ 34.000 1 PATI, PAT2, PAT3, PAT4, PAT5, PAT6, 35.000 2 PAT7, PAT8 36.000 C 364.000 0012 INTEGER RW.CL,U 365.000 C'aC. fvv%

C 368.000 0013 U = PAT2(RW,CL) * ( 2 ** (NIMAGE-2)) 369.000 1 + PAT3(RW,CL) * ( 2 ** (NIMAGE-3)) 370.000 2 + PAT4(RW,CL) * ( 2 ** (NIMAGE-4)) 371.000 3 + PAT5(RW,CL) t ( 2 ** (NIMAGE-5)) 372.000 4 + PAT6(RW,CL) * ( 2 ** (NIMAGE-6)) 373.000 5 + PAT7(RWCL) * ( 2 ** (NIMAGE-7)) 374.000 6 + PAT8(RW,CL) * ( 2 ** (NIMAGE-8)) 375.000 7 + 1 376.000 0014 RETURN 377.000 C 378.000 0015 END 379.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = COLUMN, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 15,PROGRAM SIZE = 986 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 17:05:56 PAGE P001 C 380.000 C —----------------------------------------- -------- 381.000 C 382.000 0001 SUBROUTINE SETCOF(A. B, XS. VS, U, T, 1) 383.000 C 384.000 C Setcoef sets the coefficients on the following three 385.000 C equations which are used to solve for the point (xy,z): 386.000 C 387.000 C (T11-Tl4x*)x + (T21-T24x*)y + (T31-T34x*)2 + (T41-T44x*) a 0 388.000 C (T12-Ti4y*)x + (T22-124y*)y + (T32-T34y*)z + (T42-T44y*) a 0 389.000 C (L1i1-Li4u )x + (L21-L24u )y + (L31-L34u )z + (L41-L44u ) a 0' 390.000 C 391.000 C Parameters: 392.000 C 393.000 C Input: XS - x* In the equations 394.000 C VS - y* in the equations 395.000 C U - columna of the shutter 396.000 C T - Transformation matrix for camera geometry 397.000 C L - Equation of the plane for each shutter column 398.000 C 399.000 C Output: A, B: Least square matrices 400.000 C 401.000 C —------------------- ----------—, - -- -.....-... ------------- 402.000 C 403.000 0002 REAL XS, VS. U, T(4,4), 1(4.4), A(3,3), 8(3) 404.000 C 405.000 C Begin... 406.000 0003 A(1,1) a T(1,1) - T(1,4) * XS 407.000 0004 A(1,2) a T(2.1) - T(2,4) * XS 408.000 0005 A(t,3) a T(3,1) - T(3,4) * XS 409.000 0006 B(l) ~ - ( T(4,1) - T(4,4) * XS ) 410.000 C 411.000 0007 A(2.1) a T(1,2) - T(1,4) * YS 412.000 0008 A(2.2) a T(2.2) - T(2,4) * YS 413.000 0009 A(2,3) a T(3,2) - T(3.4) * VS 414.000 0010 8(2) a - ( T(4.2) - T(4,4) * VS ) 415.000 C 416.000 0011 A(3,i) L 1(1,1) - L(1.4) * U 417.000 0012 A(3,2) " L(2.1) - 1(2,4) * U 418.000 0013 A(3.3) - L(3.1) - 1(3,4) * U 419.000 0014 B(3) ~ - ( 1(4,1) - 1(4.4) * U ) 420.000 C 421.000 C 422.000 0015 RETURN 423.000 0016 END 424.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLISTNODECK.LIOAD.NOMAP *OPTIONS IN EFFECT* NAME - SETCOF LINECNT = 57 *STATISTICS* SOURCE STATEMENTS * 16.PROGRAM SIZE a 614 *STATISTICS* NO DIAGNOSTICS GENERATED NO STATEMENTS FLAGGED IN THE ABOVE COMPILATIONS.

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 16:41:56 PAGE POO1 C ---------------------------------------------------------------- 1.000 C 3.000 0001 SUBROUTINE L2(M,N.A,B,X,EPS,11.12,13,M1,QL) 6.000 C 7.000 C ----------------------— 7.500 C 7.700 C THIS SUBROUTINE SOLVES THE LINEAR EQUATIONS AX = B. 8.000 C 8.500 C A -- M X N matrix 9.000 C EPS -- threshold for the singular vaules. 10.000 C Those S.V.'s (in ARRAY Q) smaller than EPS * QMAX 11.000 C ~ are set to zero. 11.500 C EPS1 -- precision desired for the solution. 12.000 C II -- maximum number of Iterations allowed to obtain the 13.000 C precision desired. 14.000 C 12, 13 - flags for various printouts. 15.000 C M1 -- row dimension of matrix A as declared In the calling 16.000 C routine. 17.000 C TOL -- Smallest REAL number in the particular machine used. 18.000 C 18.600 C 19.200 C REFERENCE: G.H. GOLUB & C. REINSCH 20.000 C NUMER. MATH. 14, 403-420, (1970). 21.000 C FORTRAN BY ALAN CLINE, NCAR, BOULDER, CO., (1970). 22.000 C 23.000 C SUBROUTINES CALLED: SVD Si 25.000 C 26.000 C 26.500 0002 COMMON/BLKA/LDTA,LCAL.LGCI.,LUNOUT,INFL,.IMIM2,LCOD 27.000 0003 COMMON/SUB/U(200,200U),V(11,11).Q(11)DX(11).R(200) 28.000 0004 REAL A(M1,N),B(M),X(N) 29.000 0005 DOUBLE PRECISION S 30.000 0006 DATA EPSI/I.E-05/ 31.000 C 32.000 C begin 33.000 C 34.000 C LUNOUT ALONE INITIALIZED HERE TO MINIMIZE INTERFACE 35.000 C WITH REST OF PROGRAM AND MINIMIZE CHANGES TO CANNED PROGRAM. 36.000 C 36.500 0007 LUNOUT-6 37.000 C 38.000 0008 CALL SVD(A.M,N,I,1.MI) 39.000 C 40.000 0009 QL=Q(1) 41.000 0010 QOHQ(1) 42 000 C 42.500 0011 DO 2 J=2,N 43.000 0012 QH=AMAX1(Q(J),QH) 44.000 0013 QL=AMINI(Q(d),QL) 45.000 0014 2 CONTINUE 46.000 C 46.500 0015 QL=QH/QL 47.000 C 48.000 0016 IF (12.EQ.1) WRITE(LUNOUT.104) OL.(O(I).I-iN o.N) v

0017 OH=QH*EPS 51.000 0018 DO 3 J=1,N 52.000 0019 IF(Q(J).LT.OH)Q(J)=0. 53.000 0020 3 CONTINUE 54.000 C 55.000 0021 IF (12.EQ.1) WRITE(LUNOUT, 104) QL,(Q(I),I1,N) 56.000 C 57.000 0022 CALL S1(X,B,MN) 58.000 C 59.000 0023 IF(13.NE.O)WRITE(LUNOUT, 100)(X(J),J=1,N) 60.000 0024 IF(I1.LE.O)RETURN 61.000 C 62.000 C Iterative Improvement 62.300 C 62.600 0025 XNORM=O. 63.000 0026 DO 5 J=1,N 64.000 0027 XNORM-AMAX1(XNORM,ABS(X(J))) 65.000 0028 5 CONTINUE 66.000 C 66.500 0029 IF(XNORM.EQ.O.)RETURN 67.000 C 68.000 0030 DO 10 K-1,11 69.000 C 69.500 0031 DO 7 I=1,M 70.000 0032 S=O. 71.000 0033 DO 6 J=1,N 72.000 0034 S=S+DBLE(A(I, J))*DBLE(X(J)) 73.000 0035 6 CONTINUE 74.000 0036 R(I) DBLE( B(I) ) - S 75.000 0037 7 CONTINUE 76.000 C' 76.500 0038 CALL S1(DX,R,M,N) 77.000 C 77.500 0039 DXNORM=O. 78.000 0040 DO 8 J=1.N 79.000 0041 T=X(J) 80.000 0042 X(J)=X(J)+DX(J) 81.000 0043 DXNORM=AMAX1(DXNORM,ABS(X(J)-T)) 82.000 0044 8 CONTINUE 83.000 C 83.500 0045 IF (13.NE. O) WRITE(LUNOUT, 100) (X(J),J=1,N) 84.000 0046 IF (.NOT. ( DXNORM.LE. EPSI*XNORM ) ) GO TO 11 85.000 0047 IF (13.NE. 0) WRITE(LUNOUT,102) K 86.000 C 86.500 0048 RETURN 87.000 C 87.500 0049 11 CONTINUE 88.000 0050 10 CONT I NUE 89.000 C 90.000 0051 WRITE(LUNOUT, 101) 1 91.000 0052 RETURN 92.000 C 93.000 0053 100 FORMAT(1HO,(5E12.3)) 94.000 0054 101 FORMAT(' CONVERGENCE IN ITERATIVE IMPROVEMENT DID NOT OCCUR IN', 95.000 1 12,' ITERATIONS.') 96.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) 12 10-10-83 16:41:56 PAGE P003 0055 102 FORMAT(13,' STEPS OF ITERATIVE REFINEMENT WERE TAKEN.') 97.000 0056 103 FORMAT(' ITERATION',13,' L2 NORM OF RESIDUAL ~',E12.3,/. 98.000 1 (5E12.3)) 99.000 0057 104 FORMAT(/,' CONDITION OF MATRIX a',E12.3./. 100.000 1' SINGULAR VALUES'/.(5E12.3)) 101.000 C 102.000 0058 END 103.000 *OPTIONS IN EFFECT* ID.EBCDIC.SOURCE,NOLIST.NODECK.LOAD,NOMAP *OPTIONS IN EFFECT* NAME a 2, LINECNT - 57 *STATISTICS* SOURCE STATEMENTS 58.,PROGRAM SIZE * 2284 *STATISTICS* NO DIAGNOSTICS GENERATED

C 104.000 C 105.000 C --------- -- ------- ---------------------------- - -00106.000 C 107.000 C 108.000 0001 SUBROUTINE SVD(A,M,N,IU,IV,M1) 113.000 C 114.000 C —----------------------------------------------------------- 114.500 C 114.700 C COMPUTATION OF THE SINGULAR VALUES AND COMPLETE ORTHOGONAL 115.000 C DECOMPOSITION OF A RECTANGULAR MATRIX A. 116.000 C 116.500 C A = U * DIAG(G) * V(TRANSPOSE) 117.000 C U(TRANSPOSE) * U = V(TRANSPOSE) * V = I 118.000 C 118.500 C.IF IU=1, THEN U IS COMPUTED. 119.000 C IF IV=1. THEN V IS COMPUTED. 120.000 C 121.000 C dimensions A(MN).U(M,N).V(N,N),G(N) 122.000 C 123.000 0002 REAL A(M1,N) 124.000 0003 COMMON/SUB/U(200,200),V(11,11),Q( 11),E(200) 125.000 0004 DATA TOL/1.102E-30/ 126.000 C 127.000 C begin 128.000 C 128.500 0005 EPS=4.768E-07 129.000 C 129.500 0006 DO 1 I=1,M 130.000 0007 DO 1001 J=1,N 131.000 0008 U(I,J)=A(I,d) 132.000 0009 1001 CONTINUE 133.000 0010 1 CONTINUE 134.000 C 134.500 0011 G=0. 135.000 0012. X=O. 136.000 C 136.500 0013 DO 11 I=1,N 137.000 0014 E(I )G 138.000 0015 S0=. 139.000 0016 L=I+1 140.000 C 140.500 0017 DO 2 J=I,M 141.000 0018 S=S+U(J, )*U(J.I) 142.000 0019 2 CONTINUE 143.000 C 143.500 C BEGIN HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM. 144.000 C 144.500 0020 G=O. 145.000 0021 IF(.NOT.(S.GE.TOL))GO TO 5 146.000 0022 F U( I,) 147.000 0023 G=SQRT(S) 148.000 0024 IF(F.GE.O.)G=-G 149.000 0025 H=F*G-S 150.000 0026 U(I,I)=F-G 151.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) SVD 10-10-83 16:41:57 PAGE P002 0027 IF(.NOT.(L.LE.N))GO TO 1005 152.000 0028 DO 4 J=LN 153.000 0029 S=O. 154.000 0030 DO 3 K=I,M 155.000 0031 S=S+U(K,I)*U(K,J) 156.000 0032 3 CONTINUE 157.000 0033 F'S/H 158.000 0034 DO 1004 KI1,M 159.000 0035 U(K,dJ)U(K,.)+F*U(K,I) 160.000 0036 1004 CONTINUE 161.000 0037 4 CONTINUE 162.000 0038 1005 CONTINUE 163.000 0039 5 CONTINUE 164.000 C 164.500 0040 Q(I)-G 165.000 0041 G=O. 166.000 0042 IF(.NOT.(L.LE.N))GO TO 10 167.000 0043 S0O. 168.000 0044 DO 6 J'L,N 169.000 0045 StS+U(I,J)*U(I,) 170.000 0046 6 CONTINUE 171.000 0047 IF(.NOT.(S.GE.TOL))GO TO 1010 172.000 0048 F=U(1,I+1) 173.000 0049 G=SQRT(S) 174.000 0050 IF(F.GE.O.)G=-G 175.000 0051 H=F*G-S 176.000 0052 U(, + 1 )=F-G 177.000 0053. DO 7 1JL,N 178.000 0054 E(J) U(IJ),/H 179.000 0055 7 CONTINUE 180.000 0056 DO 9 J=L,M 181.000 0057 S=0. 182.000 0058 DO 8 K=L.N 183.000 0059 S=S+U(JK)*U(I.K) 184.000 0060 8 CONTINUE 185.000 0061 DO 1009 K=L,N 186.000 0062 U(J,K)=U(J,K)+S*E(K) 187.000 0063 1009 CONTINUE 188.000 0064 9 CONTINUE 189.000 0065 1010 CONTINUE 190.000 0066 10 CONTINUE 191.000 0067 Y*ABS(Q( ) )ABS(E( )) 192.000 0068 IF(Y.GT.X)X=Y 193.000 0069 11 CONTINUE 194.000 C 194.500 C END HOUSEHOLDER'S REDUCTION 195.000 C 195.500 C BEGIN COMPUTATION OF VECTOR V 196.000 C 196.500 0070 IF(.NOT.(IV.EQ.1))GO TO 18 197.000 0071 00 17 I1 J,N 198.000 0072 IN+1-I1 199.000 0073 IF(.NOT.(L.LE.N))GO TO 165 200.000 0074 IF(.NOT.(G.NE.O.))GO TO 15 201.000

0076 DO 12 J=L,N 203.000 0077 V(J, I )=U( I,J)/H 204.000 0078 12 CONTINUE 205. 000 C 205.500 0079 DO 14 J=L,N 206.000 0080 S=O. 207.000 0081 DO 13 K=L,N 208.000 0082 S=S+U(I,K)*V(K.) 209.000 0083 13 CONTINUE 210.000 0084 DO 1014 K=L,N 211.000 0085 V(K, J)=V(K, )+S*V(K,I) 212.000 0086 1014 CONTINUE 213.000 0087 14 CONT INUE 214.000 C 214'.500 0088 15 CONTINUE 215.000 C 215.500 0089 DO 16 J=L,N 216.000 0090 V(I,J)=O. 217.000 0091 V(, I)=0. 218.000 0092 16 CONTINUE 219.000 C 219.500 0093 165 CONTINUE 220.000 C 220.500 0094 V(I, I)1. 221.000 0095 G=E(I) 222.000 0096 -L 223.000 0097 17 CONTINUE 224.000 0098 18 CONTINUE 225.000 C 225.500 C END COMPUTATION OF V 226.000 C 226.500 C BEGIN COMPUTATION OF VECTOR U 227.000 C 227.500 0099 IF(.NOT.(IU.EQ. 1))GO TO 26 228.000 0100 DO 25 I1=1,N 229.000 0101 I=N+1-I1 230.000 0102 L=I+1 231.000 0103 G=Q(I) 232.000 0104 IF(.NOT.(L.LE.N))GO TO 195 233.000 0105 DO 19 J=L,N 234.000 0106 U(I,J)=0. 235.000 0107 19 CONTINUE 236.000 0108 195 CONTINUE 237.000 C 237. 500 0109 IF(.NOT.(G.NE.O.))GO TO 23 238.000 0110 H=U(I, )*G 239.000 0111 IF(.NOT.(L.LE.N))GO TO 215 240.000 0112 00 21 J=L,N 241.000 0113 S=O. 242.000 0114 DO 20 K=L,M 243.000 0115 S=S+U(K. I )*U(K,J) 244.000 0116 20 CONTINUE 245.000 0117 F=S/H 246.000 0118 DO 1021 K=I,M 247.000 0119 U(K,J)=U(KJ)+F*U(K, ) 248.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) SVD 10-10-83 16:41:57 PAGE P004 0120 1021 CONTINUE 249.000 0121 21 CONTINUE 250.000 0122 215 CONTINUE 251.000 C 251.500 0123 DO 22 JI.,M 252.000 0124 U(J,I)-U(dI)/G 253.000 0125 22 CONTINUE 254.000 C t 254.500 0126 GO TO 1025 255.000 C ELSE 256.000 0127 23 CONTINUE 257.000 c 257.500 0128 DO 24 J=I.M 258.000 0129 U(J,1)=O. 259.000 0130 24 CONTINUE 260.000 C 260.500 0131 1025 CONTINUE 261.000 0132 U(1,I)=U(I.1)+1. 262.000 0133 25 CONTINUE 263.000 0134 26 CONTINUE 264.000 C 264.500 C END COMPUTATION OF U 265.000 C 265.500 C BEGIN DIAGONALIZATION OF THE BIDIAGONAL FORM Q USING QR ALGORITHM 266.000 C 266.500 0135 EPS=EPS*X 267.000 0136 DO 41 K1=1,N 268.000 0137 K=N+1-K1 269.000 C BEGIN BIG LOOP 270.000 0138 27 CONTINUE 271.000 C 271.500 0139 DO 28 LI=1.K 272.000 0140 L=K+1-L1 273.000 0141 IF(.NOT.(ABS(E(L)).GT.EPS))GO TO 32 274.000 0142 IF(.NOT.(ABS(Q(L-1)).GT.EPS))GO TO 29 275.000 0143 28 CONTINUE 276.000 C 276.500 0144 29 CONTINUE 277.000 c 277.500 0145 C-O. 278.000 0146 S=1. 279.000 0147 L1=L-1 280.000 0148 DO 31 IL,K 281.000 0149 F=S*E(I) 282.000 0150 E(I)=C*E(1) 283.000 0151 IF(ABS(F).LE.EPS)GO TO 32 284.000 0152 G-Q(I) 285.000 0153 H=SQRT(F*F+G*G) 286.000 0154 ( I)=H 287.000 0155 C=G/H 288.000 0156 S=-F/H 289.000 0157 IF(.NOT.(IU.EQ.1))GO TO 1031 290.000 0158 DO 30 Jd=,M 291.000 0159 VY=II1. I I^ ^^^

0161 U(J. L )=Y*C+Z*S 294.000 0162 U(J. I )=-Y*S+Z*C 295.000 0163 30 CONTINUE 296.000 0164 1031 CONTINUE 297.000 0165 31 CONTINUE 298.000 0166 32 CONTINUE 299.000 C 299.500 0167 Z Q(K) 300.000 0168 IF(L.EQ.K)GO TO 39 301.000 C 301.500 C ACCELERATION SCHEME FOR CONVERGENCE BY SHIFTING TO 302.000 C DISPLACED EIGENVALUES 303.000 C 303.500 0169 X=Q(L) 304.000 0170 Y=Q(K-1) 305.000 0171 G=E(K-1) 306.000 0172 H=E(K) 307.000 0173 F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.*H*Y) 308.000 01,74 G=SQRT(F*F+1.) 309.000 0175 IF(.NOT.(F.LT.O.))GO TO 33 310.000 0176 F=((X-Z)*(X+Z)+H*(Y/(F-G)-H))/X 311.000 0177 GO TO 34 312.000 0178 33 CONTINUE 313.000 0179- F=((X-Z)*(X+Z)+H*(Y/(F+G)-H))/X 314.000 0180 34 CONTINUE 315.000 C 315.500 0181 C=1. 316.000 0182 S 1. 317.000 0183 L2=L+1 318.000 0184 DO 38 I=L2,K 319.000 0185 G=E(I) 320.000 0186 Y=Q(I) 321.000 0187 H=S*G 322.000 0188 G=C*G 323.000 0189 Z=SQRT(F*F+H*H) 324.000 0190 E(I-1)=Z 325.000 0191 C=F/Z 326.000 0192 S=H/Z 327.000 0193 F=X*C+G*S 328.000 0194 G=-X*S+G*C 329.000 0195 H=Y*S 330.000 0196 Y=Y*C 331.000 0197 IF(.NOT.(IV.EQ.1))GO TO 36 332.000 0198 DO 35 dJ1,N 333.000 0199 X=V(J, -1) 334.000 0200 Z=V(J, I) 335.000 0201 V(J, I- )=X*C+Z*S 336.000 0202 V(J. I )=-X*S+Z*C 337.000 0203 35 CONTINUE 338.000 0204 36 CONTINUE 339.000 C 339.500 0205 Z SORT(F*F+H*H) 340.000 0206 0(I-1)=Z 341.000 0207 C=F/Z 342.000 0208 S=H/Z 343.000

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) SVD 10-10-83 16:41:57 PAGE P006 0209 F C*G+S*Y 344.000 0210 X=-S*G+C*Y 345.000 0211 IF(.NOT.(IU.EQ. 1))GO TO 1038 346.000 0212 DO 37 J-1,M 347.000 0213 Y=U(J. I-1) 348.000 0214 Z=U(J. ) 349.000 0215 U(J, I- )Y*C+Z*S 350.000 0216 U(., I )s-Y*S+Z*C 351.000 0217 37 CONTINUE 352.000 0218 1038 CONTINUE 353.000 0219 38 CONT INUE 354.000 C 354.500 0220 E(L)0-. 355.000 0221 E(K)=F 356.000 0222 Q(K)=X 357.000 0223 GO TO 27 358.000 C 358.. 500 C END BIG LOOP 359.000 C 359.500 0224 39 CONTINUE 360.000 C 360.500 C CHANGE SIGN OF NEGATIVE SQUARE ROOTS 361.000 C 361.500 0225 IF(Z.GE.O.)GO TO 41 362.000 0226 Q(K) —Z 363.000 0227 IF(IV.NE.1)GO TO 41 364.000 C 364.500 0228 DO 40 JI,N 365.000 0229. V(d,K)=-/(J,K) 366.000 0230 40 CONTINUE 367.000 C 367.500 0231 41 CONTINUE 368.000 C 368.500 0232 RETURN 369.000 0233 END 370.000 *OPTIONS IN EFFECT* IDEBCDICSOURCENOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME - SVD, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 233,PROGRAM SIZE * 5466 *STATISTICS*, NO DIAGNOSTICS GENERATED

C 371.000 C 372.000 C ----------------------------------------------------------------- 373.000 C *374.000 C 375.000 0001 SUBROUTINE S1(X.B,M,N) 380.000 C 381.000 c —---—, —---------------- ------------ 381.500 C 381.700 C COMPUTES SOLUTION VECTOR X. 382.000 C X = V * DIAG(G INVERSE) * U(TRANSPOSE) 383.000 C A = U * DIAG(G) * V(TRANSPOSE) 384.000 C 385.000 0002 REAL X(N),B(M),Y(11) 386.000 0003 COMMON/SUB/U(200,200) V(11,11),Q(11),OX(11),R(200) 387.000 C 388.000 0004 DO 3 J=1,N 389.000 0005 IF(.NOT.(Q(J).LE.O))GO TO 1 390.000 0006 Y(J)=0. 391.000 0007 GO TO 3 392.000 0008 1 CONTINUE 393.000 C ELSE 394.000 0009 S=0. 395.000 0010 DO 2 I1,M 396.000 0011 S=S+U(IJ)*B(I) 397.000 0012 2 CONTINUE 398.000 0013 Y(J)=S/Q(d) 399.000 0014 3 CONTINUE 400.000 C 400.500 0015 DO 6 I=1,N 401.000 0016 S=0. 402.000 0017 DO 5 J=1,N 403.000 0018 S=S+V(I,J)*Y(J) 404.000 0019 5 CONTINUE 405.000 0020 X(I)=S 406.000 0021 6 CONTINUE 407.000 C 407.500 0022 RETURN, 408.000 0023 END 409.000 *OPTIONS IN EFFECT* ID.EBCDIC.SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = SI, LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 23,PROGRAM SIZE = 804 *STATISTICS* NO DIAGNOSTICS GENERATED

MICHIGAN TERMINAL SYSTEM FORTRAN G(21.8) MAIN 10-10-83 16:41:59 PAGE P001 C 413.000 c —---------------------------------------------------------------- 414.000 C 415.000 C 416.000 0001 SUBROUTINE CLRTBL(MATRIX,NROW,NCOL) 417.000 C 418.000 c —--------------------------------------------------------------- 418.500 C 419.000 C * 420.000 C THIS SUBROUTINE CLEARS A 2-D TABLE (CALLED'MATRIX') OF 421.000 C DIMENSION NROW X NCOL. 422.000 C 423.000 0002 REAL MATRIX(NROW,NCOL) 424.000 C 425.000 C begin 426.000 C 426.500 0003 IF(.NOT.((NROW.LT.1).OR.(NCOL.LT.1)))GO TO 1001 427.000 0004 STOP 428.000 C 428.500 0005 1001 CONTINUE 429.000 0006 DO 1002 J=1,NCOL 430.000 0007 DO 1003 It1,NROW 431.000 0008 MATRIX(I,J)=O. 432.000 0009 1003 CONTINUE 433.000 0010 1002 CONTINUE 434.000 C 434.500 001 RETURN 435.000 C 436.000 0012 END 437.000'OPTIONS IN EFFECT* ID,EBCDIC.SOURCE.NOLIST,NODECK.LOAD,NOMAP *OPTIONS IN EFFECT* NAME = CLRTBL, LINECNT a 57 *STATISTICS* SOURCE STATEMENTS = 12,PROGRAM SIZE ~ 544 *STATISTICS* NO DIAGNOSTICS GENERATED

C 438.000 C 439.000 C ------------------- 440.000 C 442.000 C 443.000 0001 SUBROUTINE CLRVC(VECTOR,N) 447.000 C 448.000 C --------- ---------— 448.500 C 449.000 C 450.000 C THIS SUBROUTINE CLEARS AN ARRAY (CALLED'VECTOR') OF 451.000 C DIMENSION N. 452.000 C 453.000 0002 REAL VECTOR(N) 454.000 C 455.000 C begin 456.000 C 456.500 0003 IF(.NOT.(N.LT.1))GO TO 1001 457.000 0004 STOP 458.000 0005 1001 CONTINUE 459.000 C 459.500 0006 DO 1002 I1,N 460.000 0007 VECTOR(I)=O. 461.000 0008 1002 CONTINUE 462.000 C 462.500 0009 RETURN 463.000 C 464.000 0010 END 465.000 *OPTIONS IN EFFECT* ID,EBCDIC,SOURCE,NOLIST,NODECK,LOAD,NOMAP *OPTIONS IN EFFECT* NAME = CLRVC LINECNT = 57 *STATISTICS* SOURCE STATEMENTS = 10,PROGRAM SIZE = 430 *STATISTICS* NO DIAGNOSTICS GENERATED NO STATEMENTS FLAGGED IN THE ABOVE COMPILATIONS.

30 8.5. Analysis Program

pASCAL/VS RELEASE 2.1 PHAN_CAL 10/10/83 11:32:54 PAGE LINE NUMBER B P C I STMT # S 0 U R C E P R 0 G R AM v-~- ----— I --- -+ —-— 3 —— ^ —-— 3 — — + —-— 4-' —+ —--— 5 —— ^ —— 6- — ~~- -7 — ^ —- -~ —----- 1.000 PROGRAM Phan_Cal ( INPUT. OUTPUT ); 2.000 3.000 CONST 4.000 5.000 u max = 64 6.000 vmax * 64; 11.000 12.000 13.000 VAR 14.000 15.000 Phaf, 16.000 Eye_f: TEXT 17.000 18.000 Data-: ARRAY (. 1..umax, 1..v_max, 1..8.)of REAL; 19.000 20.000 Xs, Ys, U, V: INTEGER; 21.000 Xi, Y1, Z1, 22.000 X2, Y2, Z2: REAL 23.000 24.000 Dist, 25.000 Dlstsqr, 26.000 Sum, 27.000 Sumrsqr: REAL 28.000 29.000 Numdata: INTEGER 30.000 31.000 Mean, 32.000 Standev: REAL; 33.000 34.000 BEGIN 35.000 36.000 Open files ) 37.000 38.000 1 RESET ( Phaf,'UNIT2' ); 39.000 2 RESET ( Eyef,'UNIT5' ); 39.200 39.400 ( Initialize ) 39.600 39.640 1 3 FOR U:u I TO u max 00 39.680 2 4 FOR V:~ 1 TO v max DO 39.720 5 Data(.u, v, 1.): 0; 40.000 41.000 read intersection points from phantom file ) 42.000 43.000 1 6 WHILE NO) EOF ( Phaf ) DO 44.000 1 BEGIN 45.000 1 7 1 READLN ( Phaf. U, V. Xs, Ys, Xi, V1, Z1 ); 46.000 I a8 Data (.U, V, 1.): Xs 46.500 1 1 9 Data (.U, V. 2.): s; 47.000 1 10 Data (.U. V, 3.): Xi 48.000 1 1 11I Data (.U, V, 4.): Y1i 49.000 1 1 12 Data (.U, V. 5.): ZI 50.000 END;

PASCAL/VS RELEASE 2.1 PHANCAL:<MAIN> 10/10/83 11:32:54 PAGE 2 LINE NUMBER B P C I STMT N SO U R C E PROGRAM v —-4 —— I —— + —— 2 —— + —— 3 —— + —-— 4 —-----— 5 —----— 6 —-----— 7 —-----— 8 —---— 9 - -1V 51.000 52.000 { read calculated points from eye file ) 53.000 54.000 1 13 WHILE NOT EOF ( Eye_f ) DO 55.000 1 BEGIN 56.000 1 1 14 READLN ( Eye_f, U, Xs, Ys, Xi, Yi, Z1; 56.100 1i 56.200 1 1 find V ) 56.300 1 1 56.400 1 1 15 V:= 0 56.500 1 2 16 REPEAT V: V + I 56.600 1 1 17 UNTIL ( V > v_max ) (( Data(.U, V, 1.) = Xs ) &( Data(.U, V, 2..) = Ys )); 56.800 1 1 56.860 1 18 IF V <= v max 56.920 1 1 i THEN BEGIN 57.000 2 1 1 19 Data (.U, V. 6.):= Xi; 58.000 2 1 1 20 Data (.U, V, 7.): Y1 59.000 2 1 1 21 Data (.U, V, 8.): Zi1 59.500 1 1 i END 59.700 1 1 1 22 ELSE WRITELN ('* Not Found * ) 59.800 1 1 1 60.000 END 61.000 62.000 4 find average and standard deviation of distances between ) 63.000 intersection points from phantom file and calculated points ) 64.000 4 from eye file 65.000 65.100 23 WRITELN (' ----- intersection ----- ----- computed point ----' ); 65.300 24 WRITELN (' XI Y1 Z1 X2 V2 Z2 distance 65.600 25 WRITELN; 65.800 66.000 26 Numdata:= 0 67.000 27 Sum:; 68.000 28 Sum_sqr: 0 69.000 70.000 1 29 FOR U:= I TO umax DO 71.000 2 30 FOR V:= 1 TO vmax DO 72.000 1 2 31 IF Data (.U. V, 1.) > 0 THEN 73.000 1 2 BEGIN 74.000 1 2 75.000 1 1 2 32 XI:= Data (.U. V, 3.) X2: Data (.U, V, 6.); 76.000 1 1 2 34 VY1: Data (.U, V, 4.); Y2: Data (.U, V. 7.); 77.000 i I 2 36 Z1: Data (.U, V. 5.); Z2: Data (.U, V, 8.) 78.000 1 i 2 79.000 1 I 2 38 Distsqr: SQR ( X1 - X2 ) + SQR ( Yi - Y2 ) + SOR ( ZI - Z2 ); 80.000 1 1 2 39 Sum_sqr:= Sum_sqr + Dlstsqr; 81.000 i 1 2 82.000 1 1 2 40 D1st:= SORT ( Dist_sqr ); 83.000 1 i 2 41 Sum:= Sum + Dist; 84.000 1 1 2 42 WRITELN ( XI:1 0:3, Vi:10:3, Zi1:10:3, X2:12:3, Y2:10:3, Z2:10:3, Dist:12:3; 85.000 1 1 2 86.000 1 1 2 43 Numdata:= Num data + 1 87.0001 1 2

PASCAL/VS RELEASE 2.1 PHANCAL:<MAIN> 10/10/83 11:32:54 PAGE 3 LINE NUMBER B P C I STMT # S 0 U R C E P R O G R AM V —+ —------ ------— 2 —-----— 3 —-----— 4- 88.000 1 1 2 89.000 1 2 END 90.000 1 2 91.000 44 Mean:= Sum / Num data 92.000 45 Stan_dev:= SORT ( Sum_sqr / Numdata - Mean * Mean ) 93.000 93.500 46 WRITELN WRITELN 94.000 48 WRITELN (' — Average Distance:', mean:7:3 ) 94.500 49 WRITELN (' —- Standard Deviation:'. stanDev:7:3 ) 95.000 96.000 END. NO COMPILER DETECTED ERRORS OPTIONS IN EFFECT: MARGINS(1, 100). NOSEQUENCE. LANGLVL(EXTENDED), LINECOUNT(60), PAGEWIDTH(132), CHECK. DEBUG, GOSTMT, SOURCE, WARNING SOURCE LINES: 1i7; COMPILE TIME: 0.13 SECONDS; COMPILE RATE: 53867 LPM SOURCE LINES. 117; TRANSLATE TIME: 0.12 SECONDS; TRANSLATE RATE: 59338 LPM TOTAL TIME: 0.31 SECONDS; TOTAL RATE: 22308 LPM

PASCAL/VS RELEASE 2.1 PHANCAL.08/20/83 22:25:57 PAGE LINE NUMBER B P C I STMT n 0 U R C E P R O G R A M v — -+ —---- --- --- ------— 3 —------- ------— 5 -----— 6 7 8 —-+ 1.000 PROGRAM PhanCa) ( INPUT, OUTPUT ) 2.000 3.000 CONST 4.000 5.000 u max = 64 6.000 v max a 64 11.000 12.000 13.000 VAR 14.000 15.000 Phaf, 16.000 Eyef: TEXT 17.000 18.000 Data: ARRAY (. I..umax, 1..vmax, 1.. 8.) of REAL 19.000 20.000 Xs. Vs. U, V: INTEGER 21.000 XI. YI. Z1, 22.000 X2, Y2, Z2: REAL 23.000 24.000 Dist, 25.000 Distsqr, 26.000 Sum, 27.000 Sumsqr: REAL 28.000 29.000 Numrdata: INTEGER; 30.000 31.000 Mean, 32.000 Standev: REAL 33.000 34.000 BEGIN 35.000 36.000 ( Open files 37.000 38.000 1I RESET ( Phaf,'UNIT=2'; 39.000 2 RESET ( Eyef,'UNIT=5'; 39.200 39.400 ( Initialize ) 39.600 39.640 1 3 FOR U:= I TO u max DO 39.680 2 4 FOR V:= 1 TO vmax DO 39.720 5 Data(.u, v, 1.):= 0 40.000 41.000 ( read Intersection points from phantom file i 42.000 43.000 1 6 WHILE NOT EOF ( Phaf ) 00DO 44.000 1 BEGIN 45.0001 1 7 READLN ( Pha4f. U, V. Xs, Vs. Xi, VYI, Z1); 46.000 1 1 8 Data (.U, V, 1.): Xs; 46.500 1 1 9 Data (.U, V, 2.): s; 47.000 1 10 Data (.U, V, 3.):= X1 i An nfVnA 4 4* a U.- a. I. - ~

PASCAL/VS RELEASE 2.1 PHANCAL:<MAIN> 08/20/83 22:25:57 PAGE 2 LINE NUMBER B P C I STMT # S O U R C E P R O G R A M V -----— + —----- 2- + —- 3 + — 4 + ------ 5 —--- + 6 + 7 +8 + -— 9- -— + —-1V 51.000 52.000 ( read calculated points from eye file ) 53.000 54.000 1 13 WHILE NOT EOF ( Eye_f ) DO 55.000 1 BEGIN 56.000 1 1 14. READLN ( Eye_f, U, Xs, Ys, XI, Yi, Z1 ); 56.100 1 1 56.200 1 1 {find V 56.300 1 1 56.400 1 1 15 V:= 0 56.500 1 2 16 REPEAT V = V + 56.600 1 17 UNTIL ( V > v_max )| (( Data(.U, V, 1.) = Xs ) & ( Data(.U, V, 2.) = Ys )); 56.800 1 1 56.860 1 1 1 18 IF V <= v max 56.920 1 1 1 THEN BEGIN 57.000 2 1 1 19 Data (.U, V, 6.):= X1; 58.000 2 1 20 Data (.U, V, 7.): Y1; 59.000 2 1 1 21 Data (.U, V, 8.):= Z1 59.500 1 1 1 END 59.700 1 1 1 22 ELSE WRITELN ('*** Not Found ***' ) 59.800 1 1 60.000 END 61.000 62.000 { find average and standard deviation of distances between 63.000 { Intersection points from phantom file and calculated points ) 64.000 ( from eye file 65.000 65.100 23 WRITELN (' ----- intersection ----- ----- computed point --—'; 65.300 24 WRITELN (' X Y1 Z1 X2 Y2 Z2 distance' ) 65.600 25 WRITELN 65.800 66.000 26 Numdata:= 0; 67.000 27 Sum:= 0; 68.000 28 Sum_sqr: 0 69.000 70.000 1 29 FOR U:= 1 TO u max DO 71.000 2 30 FOR V:= 1 TO vmax DO 72.000 1 2 31 IF Data (.U, V, 1.) > 0 THEN 73.000 1 2 BEGIN 74.000 I 2 75.000 1 1 2 32 X1:= Data (.U, V, 3.); X2:= Data (.U, V, 6.); 76.000 i i 2 34 Y1:= Data (.U, V, 4.); Y2:= Data (.U, V, 7.); 77.000 1 1 2 36 Z1:= Data (.U, V, 5.); Z2:= Data (.U. V, 8.); 78.000 1 1 2 79.000 1 1 2 38 Dist_sqr:= SQR ( X1 - X2 ) + SQR ( Y1 - Y2 ) + SOR ( Z1 - Z2 ); 80.000 1 1 2 39 Sum_sqr:= Sum_sqr + Dlst_sqr; 81.000 1 1 2 82.000 1 1 2 40 Dist:= SORT ( Dist_sqr ) 83.000 1 1 2 41 Sum:= Sum + Dst; 84.000 1 1 2 42 WRITELN ( XI:10:3, Y1:10:3, Z1:10:3, X2:12:3, Y2:10:3, Z2:10:3. Dlst:12:3 ); 85.000 1 1 2 86.000 1 1 2 43 Numdata:= Num data + 1 87.000 1 1 2

PASCAL/VS RELEASE 2.1 PHANCAL:<MAIN> 08/20/83 22:25:57 PAGE 3 LINE NUMBER B PCI STMT SOUR C E PROGRAM V —-+ —— I1 —— + —— 2 —— + —— 3 —— + —— 4 —----— 5 —-----— 6 —----— 7 —------- 88.000 1 1 2 89.000 1 2 END; 90.000 1 2 91.000 44 Mean:= Sum / Numdata; 92.000 45 Standev:= SORT ( Sumsqr / Numdata - Mean * Mean ); 93.000 93.500 46.WRITELN; WRITELN 94.000 48 WRITELN (' --- Average Distance:', mean:7:3 ); 94.500 49 WRITELN (' --- Standard DeviatIon:' stanDev:7:3 ) 95.000 96.000 END. NO COMPILER DETECTED ERRORS OPTIONS IN EFFECT: MARGINS( 1,100), NOSEQUENCE, LANGLVL(EXTENDED), LINECOUNT(60), PAGEWIDTH(132), CHECK, DEBUG, GOSTMT, SOURCE, WARNING SOURCE LINES: 117; COMPILE TIME: 0.14 SECONDS; COMPILE RATE: 50289 LPM SOURCE LINES: 117; TRANSLATE TIME: 0.12 SECONDS; TRANSLATE RATE: 58792 LPM TOTAL TIME: 0.33 SECONDS; TOTAL RATE: 21596 LPM

Q,~~~9. SQampls~e~~ Run~~31 g* Sample Run

,,Sn exec.e $.00, $13.97T fenrpty -5 I/Done. $.00. $13.97T /,$rti )iamp:eyve.obj+lamnp:stup.obj.2+naas:imsI 3=-3 4=-4 5=-5 t=3 #Eyecut loi beg ris CONDITION OF MATRIX = 0.375E+03 SINGULAR VALUES 0.540EF03 0.331E+03 0.390E+03 0.217E+02 0.172E+02 0.159E+02 0.590E+01 0.302E+01 0.219E+01 0.144E+01 0. 546E+01 CONDITION OF MATRIX 0.375E+03 SINGULAR VALUES 0.54OE+03 0.331E403 0.390E+03 0.217E+02 0.172Er02 0.159E*02 0.590EI01 0.302E401 0.219E+01 0.144E+01 0.546E+01 r Matrix determined -1.7496386 -0.8302813 0.0 -0.0260207 0.0814062 -1.3786592 0.0 -0.0421989 -0.4219891 0.9700564 0.0 -0.0260956 16.3746338 16.4649811 0.0 1.0000000 CONDITION OF MAIRIX = 0.658E+02 SINGULAR VALUES 0. 100E403 0.644E+02 0.770E+02 0.737E-101 0.337E+401 0.261E+O1 0.152E+01 CONDITION OF MATRIX = 0.658E+02 SINGULAR VALUES 0.100E+03 0.644E+02 0.770E+02 0.737E+01 0.337E+01 0.261E+01 0.152E01 L Matrix determiiiined -0.3377837 0.0 0.0 -0.0306186 0.2668607 0.0 0.0 -0.0176777 -0.1590991 0.0 0.0 -0.0353554 4.5000000 0.0 0.0 1.0000000 hFyxrcut Ilon terminated $15. $14.12T

2 2' Tinat 3 12 2 3 Geome t 6 15 2 3'mat 6 15 2 4 Geome t 8 17 2 Tinmat 8 17 2 5 Geomet 11 19 2 5 Tmat 11 19 2 6 Geomet 13 21 2 6 Tmat 13 21 2 7 Geomet 16 24 2 7 Tmat 16 24 2 8 Geomet 18 26 2 8 Triat 18 26 3 1 Geomet 3 9 3 1 Tmat 3 9 3 2 Geomet 6 11 3 2 Tlmat 6 t 1 3 3 Geomet 8 14 3 3 Tmat 8 14 3 4 Geomet 11 16 3 4 Tmat 11 16 3 5 Geomet 14 18 3 5 Tmat 14 18 3 6 Geomet 16 21 3 6 Tmat 16 21 3 7 Geomet 19 23 3 7 Tmat 19 23 3 8 Geomet 21 25 3 8 Tmat 21 25 4 1 Geomet 6 7 4 1 Tmat 6 7 4 2 Geomet 9 10 4 2 Tmat 9 10 4 3 Geomet 11 13 4 3 Tmat 11 13 4 4 Geomet 14 15 4 4 Tmat 14 15 4 5 Geomet 17 18 4 5 Tmat 17 18 4 6 Geomet 20 20 4 6 Tmat 20 20 4 7 Geomet 22 22 4 7 Tmat 22 22 4 8 Geomet 25 25 4 8 Tmat 25 25 5 1 Geomet 9 6 5 1 TTmat 9 6 5 2 Geomet 12 9 5 2 lmat 12 9 5 3 Geomet 15 11 5 3 Tmat 15 11 5 4 Geomet 18 14 5 4 Tnmat 18 14 5 5 Geomet 20 17 5 5 Tmat 20 17 5 6 Geomet 23 19 5 f6 Tmat 23 19 H 7, r -tn n %

5 8 Tmat 29 21 6 1 Goome t 12.I 6 I mat f2 41 6 2 Goomet 15 7 6 2 ma t 15 7 G 3 Geone t 19 0 G:' ITm<at 19 10 6 4 Geomet 22 13 6 4 Tnat 22 13 6 5 Geomnet 24 16 G 5 Tmeat 24 16 6 6 Geomnet 27 18 6 G Tmat 27 18 6 7 Geomet 30 21 6 7 I mat 30 21 VMAT "'' not observed +* 6 8 GEO "'* not observed t+' 6 8 TMAT t*-' nrot observed *+ 6 8 7 I Geomet 16 2 7 1 Tmat 16 2 7 2 Geomet 20 5 7 2 na t 20 5 7 3 Geomet 23 8 7 3 1 emat 23 8 7 4 Geonmet 26 11 7. Tiat 26 11 7 5 Geomet 29 14, 7 5 Tlmat 29 14 7 6 Geomet 32 17 7 6 Tmat 32 17 VMAT ** not observed ** 7 7 GEO * i not observed * 7 7 TMAT **+ not observed +* 7 7 VMAT + * not observed *** 7 8 GEO + not observed + 7 8 TMAT + nlot observed tt* 7 8 VMAT +* not observed + 8 1 GEO * "- not observed t * 8 1 TMAT **+ not observed **+ 8 1 8 2 Geomet 24 3 8 2 Tmat 24 3 8 3 Geomet 28 6 8 3 Tmat 28 6 8 4 Geoie t 31 10 8 4 Imat 31 10 VMAT 1"' not observed *++ 8 5 GEO +' not observed * 8 5 TMAT + not observed *+ 8 5 VMAT *' not observed ** 8 6 GEO + not observed * 8 6 TMAT +*+ not observed +* 8 6 VMAT "' not observed $* 8 7 GEO +* not observed 1 * 8 7 TMAT +** not observed + 8 7 VMAT *' not observed * 8 8 GEO +11 not observed 8' 8 IMAT'* not observed * 8 8 Niin. face 3'i F-ecu t inn tet minl ate (

lengtl of the side of a cube --- Length 8 8.0000 Step In u aned v -- Step ---- rotatiol of the cube ---- Enter axis ( x, y, z, or / ) and angle x 10 x 10.0000 Enter axis ( x, y, z, or / ) and angle y 45 y 45.0000 Enter axis ( x, y. z, or / ) and angle z 30 z 30.0000 Enter axis I x, y. z, or / ) and angle / 0 0.0 ---- Tranlslation of the cube --- Enter translation along x, y, and z axis e.g. 3.8 10 8.5 0 0 0.0 0.0 0.0 ---- natri >: -34.8575208 -15.1447174 0.0 -0.4330128 3.6250092 -26.2314068 0.0 -0.7499998 -8.2500031 19.4628063 0.0 -0.5000002 330.000000 330.0000000 0.0 20.0000000 L matrix ---- -6.7556749 -7.6516563 0.0 -0.6123724 5.3372140 -4.4194168 0.0 -0.3535533 -3.1819810 2.4748723 0.0 -0.7071069 90. 0000000 90.0000000 0.0 20.0000000 VMAT * not observed *** 1 1 GEO * not observed *** 1 1 TMAT ~' not observed * 1 1 1 2 Geomet 1 14 1 2 Tmat 1 14 1 3 Geoie t 3 16 1 3 imat 3 16 1 4 Geomet 6 18 1 4 Tmat 6 18 1 5 Geomet 8 20 1 5 Tmat 8 20 6 Geome t 10 22 1 6 Tinat 10 22 1 7 Geomet 13 24 1 7 Tmat 13 24 1 8 Geomet 15 26 8 Tmat 15 26 39 I (rsncnt I' I*

$Log Output: LAMP, 01:32:11 Thu Aug 04/83 #so50 exec.f # $.00. $13.79T #empty -2 #Done. $.00. $13.79T Nenipty -3 #Dotne. $.00. $13.80T #empty -4 #Done. # $.00. $13.80T *r lamp:phan.obj4+pascalvsiib 2=-2 3=-3 4=-4 t=3 #Execution begins Enter Parameters ---- shutter parameters ---- focal point -- r 20 focal point -- tlieta 30 focal point -- phi 45 shutter plnne -- D 4 shutter plane -- Su 2 shutter plane -- Sv 2 shutter plan? -- Nu 8 shutter plane — N 8 ---- camera parameters focal point -- r 20 focal point -- theta 6O focal point -- phi 60 camera plane -- D 4 camera plane -- Sx 2 camera plane -- Sy 2 camera plane -- Nx 32 camniera plane -- Ny 32 Equation of the plane to he mapped -- A. B, C, D 0 0 1 0.0 0.0 1.0000000 1.0000000 Enter Calibration Parameters I

4 4F' 4 F 4 * 1 t ~ I' t +~ 2 + - * * 4- + 3 k f ~ f < t k 4 4 t 4 + 44 2 4 I 5 4 I- 4 7 4 4 7' I,,, 4 4'I3 0I. ll'III. 3 -------- --— 69230:3 —-- 33 - -------- 6923 —3 ---- 92303 3- - 9 3 -- 230 23 -- - 62303 2 ---- Job number' 692303 University of Michicint Terminal System (Model U1_223) 01:35-47 EDT Thu Atlg 01/83'(l / I MMMMM MMMMM TTTTTTTTTTTTTT fTTTI[TTTT SSSSSSSSS MMMMMM MMMMM M M TTTTTTTTT TTTTT Tl TTTTTTTT SSSSSSSSSSSSS MMMMMMM MMMMMMM TTrTTTTTTTTTTTTTTTTTTTTT SSSSSSSSSSSSSSSSS MMMMMMMM MMMMMMMM TTTTT SSSSSS SSSSSSS MMMMMMMMM MMMMMMMMM TTTTT SSSSS SSSSS MMMMMMMMMM MMMMMMMMMM TTT T SSSSS MMMMM MMMMM MMMMM MMMMM T T TT SSSSSS MMMMM MMMMM MMMMM MMMMM lTTTT SSSSSSSSSSSSS MMMMM MMMMM MMMMM MMMMM TTTTT SSSSSSSSSSSSS MMMMM MMMMMMMMMMM MMMMM TTTT SSSSSSSSSSSSS MMMMM MMMMMMMMM MMMMM TTTIT SSSSSSS MMMMM MMMMMMM MMMMM TTfT SSSSS MMMMM MMMMM MMMMM TTITT SSSSS MMMMM MMM MMMMM TTTTT SSSSS SSSSS MMMMM MMMMM TTITT SSSSSSS SSSSSSS MMMMM MMMMM TTTTT SSSSSSSSSSSSSSSSS MMMMM MMMMM TTTT1 SSSSSSSSSSSSS MMMMM MMMMM TTTTT SSSSSSSS LL AAAAAAAAA tiM MM FrPPPPPPPPPP IL. AAAAAAAAAAAA MMM MMM PPIPPPPPPPPP L AA AA A'MMMM MM!rr.IM PP PP LL. AA AA MM MM MM MM PP PP L AA AA MM MMMM MM P P PP LL AAAAAAAAAAA MM M1 M MM ['PPPPPPPPPP LL.AAAAAAAAAAA MM MM PPPPPPPPPP LL AA AA MM MM PP LL AA AA MM MM PP LI AA AA MM MM PP LLLLLLLI. LLI. AA AA MM MM PP LLLLLLLLLLLL AA AA MM MM PP

aso exec.a' $.01. $14.13T ftI lamp:analy.obj+tpascalvslib 2=-2 5=-5 t=3.Execut orion bgris intersection ----- ----- computed point ---- XI YI ZI X2 Y2 Z2 distance 8.948 -3.076 -1.000 8.828 -3.004 -0.538 0.183 7.220 -4.942 -1.000 7.522 -4.298 -0.270 1.019 5.086 -7.246 -1.000 5.144 -7.038 -0.667 0.397 2.382 -10.164 -1.000 2.966 -9.240 -0.326 1.284 -1.153 -13.980 -1.000 0.200 -12.039 0.107 2.613 -5.974 -19.184 -1.000 -6.096 -19.072 -0.437 0.587 -12.937 -26.701 -1.000 -12.493 -25.686 0.233 1.658 9.444 0.077 -1.000 9.174 -0.221 -1.175 0.438 7.928 -1.310."1.000 8.065 -1.195 -1.036 0.183 6.093 -2.989 -1.000 6.010 -3.000 -0.777 0.238 3.826 -5.064 -1.000 4.282 -4.517 -0.559 0.837 0.955 -7.691 -1.000 0.864 -7.818 -1.149 0.216 -2.800 11.127 -1.000 -2.309 -10.648 -0.899 0.693 -7.920 -15.812 -1.000 -9.297 -16.882 -0.348 1.861 -15.316 -22.580 -1.000 -16.639 -23.430 0.232 1.997 8.513 1.690 -1.000 8.351 1.611 -0.729 0.325 6.908 0.456 -1.000 6.775 0.324 -1.180 0.260 4.966 -1.037 -1.000 5.303 -0.677 -0.395 0.781 2.566 -2.882 -1.000 2.736 -2.736 -0.906 0.243 -0.473 -5.219 -1.000 -0.999 -5.730 -1.649 0.980 -4.448 -8.274 -1.000 -4.290 -8.059 -0.429 0.630 -9.867 -12.441 -1.000 -12.269 -14.350 -1.379 3.092 -17.696 -18.459 -1.000 -20.681 -20.783 -1.177 3.787 7.582 3.302 -1.000 7.358 3.154 -1.110 0.291 5.888 2.222 -1.000 5.620 2.048 -1.057 0.325 3.839 0.915 -1.000 4.130 1.133 -0.376 0.722 1.307 -0.700 -1.000 1.462 -0.597 -0.930 0.199 -1.901 -2.746 -1.000 -2.332 -3.011 -0.815 0.539 -6.095 -5.421 -1.000 -8.233 -6.826 -1.787 2.677 -11.814 -9.069 -1.000 -14.065 -10.548 -1.840 2.821 -20.075 -14.338 -1.000 -29.069 -20.126 -1.975 10.740 6.651 4.915 -1.000 6.404 4.771 -0.691 0.421 4.869 3.988 -1.000 4.678 3.873 -0.671 0.398 2.712 2.867 -1.000 2.493 2.766 -1.256 0.351 0.0417 1.482 -1.000 -0.533 1.197 -1.326 0.723 -3.328 -0.273 -1 000 -3.147 -0.200 -0.578 0.465 -7.742 -2.568 -1.000 -8.830 -3.108 -1.517 1.320 -13.761 -5.697 -1.000 -18.048 -7.889 -1.730 4.870 -22.45-1 10.217 — 1.000 -37.726 -17.993 -4.267 17.446 5.720 6.527 -1.000 5.717 6.452 -0.472 0.533 3.849 5.755 -1.000 4.055 5.770 -0.509 0.533 1.584 4.819 -1.000 1.107 41.641 -1.133 0.527 -1.213 3.663 -1.000 -2.072 3.351 -1.301 0.963 -4.756 2.200 -1.000 -4.8411 2.122 -0.703 0.318 -9.389 0.285 -1.000 -10.885 -0.224 -1.766 1.756 -15.707 -2.3293 -1.000 -20.841 -4.262 -2.292 5.637 4.788 8.140 -1.000 4.823 8.003 0.345 0.670 2.829 7.521 l.O0O 2.560 7.417 -0.919 0.300 0.457 6.771 -1.000 0.214 6.724 -1.130 0.280 -2.472 5.845 -1.000 -2.951 5.788 -1.414 0.636 -6.184 4.672 -1.000 -7.456 4.455 -1.819 1.529 -11.037 3.138;1.000 -14.380 2.408 -2.441 3.714 1.809 9.287 -1.000 2.062 9.167 -0.415 0.618

-0.670 8.723 -1.C00 -0.970 8.707 -1.171 0.34G -3.732 8.027 -1.000 -4.269 7.926 -1.069 0.551 -- Average Distance 1.597 -- Standard Deviation 2.780 FExecution terminated i $.04, $14171:mts $.13. $14.32T Fcopy -iog2'printit