DLTDSP README is a companion dataset for file DLTDSP FORTRAN which contains the following photogrammetric and rigid-body routines: (1) DLTC: camera calibration subroutine based on the Direct Linear Transformation for a spatial or planar calibration point distribution, with various constraints on the equivalent, internal camera parameters (image skew factor, principal point, principal distances); (2) CNVDLT and DLTCNV: 1-to-1 mapping between DLT-parameters and `conventional' (internal/external) camera parameters; (3) DLTPOS and DISP3D: iterative, minimum-variance 3-D point reconstruction from DLT-parameters and measured image data, and rigid-body reconstruction from 3-D point data, with a number of supplementary subroutines. Two testroutines are provided in this file (DLTDSP README): (A) TSTDLTC: simulated generation of calibration data for DLTC (B) TSTDLT: test for CNVDLT/DLTCNV and simulated generation of noisy position data for DLTPOS It is recommended that the routines in file DLTDSP FORTRAN are compiled and placed in a library for easy access. The routines have been tested with MS-Fortran 5.0 under MS-DOS 4.01, and were found to compile properly under VAX/VMS FORTRAN. Some minor editing may be necessary for running the test programs under other compilers than MS-FORTRAN, especially w.r.t. the backslash (\) edit descriptor for CR/LF suppression in the two test programs and w.r.t. random noise generation in subroutine GAUSS. COMMENTS Public/private: while DLTDSP README is a public file, DLTDSP FORTRAN has been classified as private, i.e., for Biomch-L subscribers only who can retrieve it by sending the SEND DLTDSP FORTRAN request to the LISTSERVer from their subscription address (NB: the LISTSERVer may have matching problems if the request is sent via a different mailer than the one from which subscription took place, e.g., EARN/BITNET's jnet% mailer v. Internet mailers). If the file cannot be retrieved, the Biomch-L moderators can send a GIVE DLTDSP FORTRAN TO YOUR_ID@YOUR_SITE request to the LISTSERVer. Least-squares: in order to maximize numerical stability, the iterative procedures in these routines utilize orthogonal triangularization by Givens rotations rather than normal equations [1]; convergence assessment takes place by comparing the residual standard deviation with the norm of the right-hand-side of the linearized equations, at each iteration. DLTC: for a truely spatial calibration point distribution (i.e., at least 6 known points, no three of which should be collinear), DLTC can be run in unconstrained mode unless the camera is very far from the landmark distrib- ution. If the camara image axes are known to be orthogonal, this constraint can be imposed [3]. If a camera's principal point is known, also this value may be imposed. Furthermore, equality of the principal distances or known values for the principal distances may be imposed [3]. These constraints are par- ticularly useful if the landmark distribution is shallow or planar; however, it is recommended that a calibration distribution encompass the full volume in which 3-D reconstruction is to occur in order to avoid extrapolation errors. DLTCNV and CNVDLT: there is a 1-to-1 relationship between the conventi- onal, internal/external camera parameters (position, attitude angles; principal point, principal distances, image skew) and the DLT-parameters [1,3]. DLTPOS: the procedure iterates the fundamental DLT-equations in order to minimize nonlinear image noise biases, and it accomodates weight factors for different noise levels per camera. DISP3D: given two sets of corresponding 3-D co-ordinates, DISP3D assesses the least-squares optimal rotation matrix and translation vector between the two sets, and an optional scaling factor. The points in each set must be non-collinear. If the two sets are truely spatial, the routine will also determine whether the displacement is a pure or an impure rotation (i.e., a reflection). DISP3D is a pre- cursor to the procedure described in [2]. These routines have been tested to a limited extent; they are provided `as is' for non-commercial purposes, and the author would be grateful for comments on experiences obtained therewith in order to improve the package. See the guidelines in file BIOMCH-L INFO at LISTSERV@HEARN. REFERENCES: [1] H.J. Woltring (1980), Planar control in multi-camera calibration for 3-D gait studies. Journal of Biomechanics 13(1), 39-48. [2] F.E. Veldpaus, H.J. Woltring & L.J.G.M. Dortmans (1988), A least- squares algorithm for the equiform transformation from spatial marker co-ordinates. Journal of Biomechanics 21(1), 45-54. [3] H.J. Woltring (1991), One hundred years of photogrammetry in biolocomotion. In: A. Cappozzo, M. Marchetti & V. Tosi (Eds), Biolocomotion: A century of Research using moving pictures, (Formia, Italy 1989) [in press, 1991]. Please communicate any comments or suggestions to: Herman J. Woltring 1991-08-08 Brussellaan 29, NL-5628 TB EINDHOVEN, The Netherlands Tel. +31.40.480 869, voice/fax/modem +31.40.413 744 1991-08-15 (rev.) - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TSTDLTC.FOR, Woltring 1991-08-15 C C*********************************************************************** C C TSTDLTC (revised, 1991-08-15) C C Purpose: C ******* C C Testprogramme for subroutine DLTC (DLT camera calibration with C constraints, initialized with DISP3D) C C Interactive parameters: C ********************** C C SD Image noise s.d. per coordinate (in image units) C IC Constraint switch: C <= 0: no constraints C >= 1: ALPHA == 0.0 C >= 2: prior X0 and Y0 (2048,2048)' C = 3: Cx == Cy C >= 4: prior Cx and Cy (8192,8192)' C NP Number of control points C (points to be entered manually, in XYZ order) C SDx,SDy,SDz Object point s.d. per coordinate C C*********************************************************************** C PROGRAM TSTDLTC C IMPLICIT REAL*8 (Q) PARAMETER (PI=3.1415926536, RADDEG=180.0/PI, DEGRAD=PI/180.0, 1 MP=32, X0=2048.0, Y0=X0, CX=8192.0, CY=CX) DIMENSION XYZP(3,MP), XYI(2,MP), VARP(3), VARI(2), DLT(11), 1 DU(11,11), CNV(11),CNV2(11),TMP(3),XYZ(3),ROT(3,3),XYZI(3,MP) C C*** Assess object coordinates C 10 WRITE(*,700) 700 FORMAT('0NP = '\) READ(*,710) NP 710 FORMAT(I5) NP = MIN0(MP,MAX0(3,NP)) WRITE(*,720) NP 720 FORMAT('0NP =',I3,'; Enter control point co-ordinates:') DO 20 I=1,NP WRITE(*,730) I 730 FORMAT(I5,': '\) READ(*,740) (XYZP(J,I),J=1,3) 740 FORMAT(3E15.0) 20 CONTINUE C C*** Get object and noise parameters C 30 WRITE(*,750) 750 FORMAT('0SD, IC, SDxyz = '\) READ(*,760,END=100) SD,IC,VARP 760 FORMAT(E15.0,I5,3E15.0) IF (SD.LT.0.0) GO TO 100 IF ((VARP(1).eq.0.).and.(VARP(2).eq.0.).and.(VARP(3).eq.0.)) THEN VARI(1) = 1.0 ELSE VARI(1) = SD * SD ENDIF VARI(2) = VARI(1) C*** Defaults testing IF (IC.lt.0) IC = 0 C C*** Change object space s.d.'s into variances C DO 40 J=1,3 VARP(J) = VARP(J)**2 40 CONTINUE C C*** Get camera parameters C 50 WRITE(*,770) 770 FORMAT('0Camera pos, phi: '\) READ(*,780) (CNV(J),J=1,6) 780 FORMAT(6E15.0) CALL MATSCL(3,1,CNV(4),3,DEGRAD,CNV(4),3) CALL PHIROT(CNV(4),ROT,1) DEN = DINPRD(1,3,ROT(3,1),3,CNV,1,0.0) IF (DEN.EQ.0.0) STOP 'DLT: singular camera station' CNV(7) = X0 CNV(8) = Y0 CNV(9) = CX CNV(10) = CY CNV(11) = 0.0 C C*** Assess image data and copy into real*8 for DISP3D C DO 70 IP=1,NP CALL MATDIF(3,1,XYZP(1,IP),3,CNV,3,TMP,3) CALL PRAB(3,3,1,ROT,3,TMP,3,XYZ,3) XYI(1,IP) = X0 - CX * XYZ(1) / XYZ(3) + GAUSS(SD) XYI(2,IP) = Y0 - CY * XYZ(2) / XYZ(3) + GAUSS(SD) XYZI(1,IP) = (XYI(1,IP) - X0) / CX XYZI(2,IP) = (XYI(2,IP) - Y0) / CY XYZI(3,IP) = 1.0 60 CONTINUE 70 CONTINUE C C*** Estimate position and attitude via DISP3D C S = 0.0 CALL DISP3D(XYZP,XYZI,NP,CNV2,ROT,S) IF (S.eq.0.D0) THEN WRITE(*,*) 'DISP3D failure' GO TO 30 ENDIF CALL ROTPHI(ROT,CNV2(4)) CALL MATCOP(5,1,CNV(7),5,CNV2(7),5) CALL CNVDLT(CNV2,DLT) CALL DLTCNV(DLT,CNV2,.FALSE.) WRITE(*,799) S 799 FORMAT('0DISP3D: S =',1PE15.6) WRITE(*,800) (CNV2(J),J=1,3),(RADDEG*CNV2(J),J=4,6), 1 (CNV (J),J=7,11) C C*** Estimate DLT parameters and convert vice-versa C C*** Estimate DLT parameters IT = 4 !initialize CALL DLTC(XYZP,XYI,NP,VARP,VARI,IC,CNV(7),DLT,DU,SDRES,IT) WRITE(*,790) IT,SDRES,DLT 790 FORMAT('0DLTCAL: IT =',I3,', SDres =',1PE15.6/ 1 ' x:',4E15.6/' y:',4E15.6/' z:',3E15.6/) IF (IT.lt.0) GO TO 30 C*** Convert DLT to standard parameters CALL DLTCNV(DLT,CNV2,.FALSE.) WRITE(*,800) (CNV2(J),J=1,3),(RADDEG*CNV2(J),J=4,6), 1 (CNV2(J),J=7,11) 800 FORMAT(' Standard parameters:'/ 1 ' position: ',3F11.6/' angles: ',3F11.3/ 2 ' pr. point:',2F11.3/' pr. dist.:',2F11.3/ 3 ' alpha: ',1PE15.6) C GO TO 30 C 100 STOP 'DLTCTST ready' END - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TSTDLT.FOR, Woltring 1989-04-07 (1991-08-08: MS-DOS) C C************************************************************************** C C Testprogramme for subroutines C C (1) CNVDLT, conversion from conventional to DLT parameters C (2) DLTCNV, conversion from DLT to conventional parameters C (3) DLTPOS, position reconstruction from DLT parameters and C noisy image data in two cameras. C C Useful camera configurations are: #1 5,0,5,0,-45,0 C #2 -5,0,5,0,+45,0 C C (Positions X, Y, Z; Cardanic angles in X-Y-Z sequence) C C The principal point is at (2048,2048)', and the principal distances C are (8192,8192), while the image obliqueness (skew) factor ALPHA = 0.0. C C The noise is uncorrelated, Gaussian, and additive to the image data, C with standard deviation SIGMA. C C************************************************************************** C PROGRAM TSTDLT PARAMETER (PI=3.1415926536, DEGRAD=PI/180.0, RADDEG=180.0/PI) DIMENSION CNVPAR(11,2), DLTPAR(11,2), XY0CXY(5), W(2), 1 XYZP(3), XYI(2,2), XYZR(3), XYZD(3) DATA XY0CXY/2*2048.0,2*8192.0,0.0/ !Internal parameters C C*** Get camera parameters and type out after CNV --> DLT --> CNV C 10 WRITE(*,700) 700 FORMAT('0Camera #: X, Y, Z, chi, eta, zeta: ') DO 20 I=1,2 WRITE(*,710) I 710 FORMAT(' ',I1,': '\) READ(*,720) (CNVPAR(J,I),J=1,6) 720 FORMAT(6E15.0) CALL MATSCL(3,1,CNVPAR(4,I),3,DEGRAD,CNVPAR(4,I),3) CALL MATCOP(5,1,XY0CXY,5,CNVPAR(7,I),5) CALL CNVDLT(CNVPAR(1,I),DLTPAR(1,I)) CALL DLTCNV(DLTPAR(1,I),CNVPAR(1,I),.FALSE.) WRITE(*,725) (CNVPAR(J,I),J=1,3), 1 (RADDEG*CNVPAR(J,I),J=4,6) 725 FORMAT(1P3E15.5, 0P3F10.4) W(I) = 1.0 20 CONTINUE C C*** Get noise level C 30 WRITE(*,730) 730 FORMAT('0Sigma: '\) READ(*,720) SIGMA C C*** Get landmark position and convert to image data C 40 WRITE(*,740) 740 FORMAT('0Xp, Yp, Zp: '\) READ(*,720) XYZP DO 50 I=1,2 CALL POSDLT(XYZP,DLTPAR(1,I),XYI(1,I),SIGMA) 50 CONTINUE C C*** Reconstruct position and type results C CALL DLTPOS(XYI,2,W,DLTPAR,11,2,XYZR,SDRES) !reconstruct XYZr CALL MATDIF(3,1,XYZR,3,XYZP,3,XYZD,3) !assess differences XYZd IF (SDRES.GE.0.0) THEN WRITE(*,750) XYZR, XYZD, SDRES 750 FORMAT('0XYZr =',1P3E12.4/' XYZd =',3E12.4,', SDres =',E12.4) ELSE WRITE(*,760) IFIX(-SDRES) 760 FORMAT('0DLTPOS error code',I4) ENDIF WRITE(*,770) 770 FORMAT('0go to (0=point, 1=noise, 2=cameras, else=exit): '\) READ(*,780) I 780 FORMAT(I2) IF ((I.lt.0).or.(I.gt.2)) STOP 'TSTDLT ready' GO TO (40,30,10), I + 1 END C SUBROUTINE POSDLT(XYZP,DLT,XYI,SIGMA) DIMENSION XYZP(3), DLT(11), XYI(2) C DEN = DINPRD(1,3,XYZP,1,DLT(9),1,1.0) DO 10 I=1,2 J = (I-1)*4 + 1 XYI(I) = DINPRD(1,3,XYZP,1,DLT(J),1,DLT(J+3)) / DEN + 1 GAUSS(SIGMA) 10 CONTINUE C RETURN END