1************************************************************************ 2* Support routines for quadric surfaces * 3************************************************************************ 4* EAM Jun 1997 - initial version, supports version 2.4(alpha) of render 5* EAM May 1998 - additional error checking to go with Parvati/rastep 6* EAM Jan 1999 - version 2.4i 7* EAM Mar 2008 - Gfortran optimization breaks the object accounting. 8* No real solution yet. Move qinp.f into separate file 9************************************************************************ 10* 11* Quadric surfaces include spheres, cones, ellipsoids, paraboloids, and 12* hyperboloids. The motivation for this code was to allow rendering 13* thermal ellipsoids for atoms, so the other shapes have not been 14* extensively tested. 15* A quadric surface is described by 10 parameters (A ... J). 16* For efficiency during rendering it is also useful to know the center and 17* a bounding sphere. So a QUADRIC descriptor to render has 17 parameters: 18* 14 (object type QUADRIC) 19* X Y Z RADLIM RED GRN BLU 20* A B C D E F G H I J 21* 22* The surface itself is the set of points for which Q(x,y,z) = 0 23* where 24* Q(x,y,z) = A*x^2 + B*y^2 + C*z^2 25* + 2D*x*y + 2E*y*z + 2F*z*x 26* + 2G*x + 2H*y + 2I*z 27* + J 28* 29* It is convenient to store this information in a matrix QQ 30* | QA QD QF QG | QA = A QB = B QC = C 31* QQ = | QD QB QE QH | QD = D QE = E QF = F 32* | QF QE QC QI | QG = G QH = H QI = I 33* | QG QH QI QJ | QJ = J 34* 35* Then Q(x,y,x) = XT*QQ*X where X = (x,y,z,1) 36* The point of this is that a 4x4 homogeneous transformation T can be 37* applied to QQ by matrix multiplication: QQ' = TinvT * QQ * Tinv 38* 39* The surface normal is easily found by taking the partial derivatives 40* of Q(x,y,z) at the point of interest. 41************************************************************************ 42* TO DO: 43* - can we distinguish an ellipsoid from other quadrics? Do we care? 44************************************************************************ 45 46 47CCC Calculate the matrices for quadric transformation 48CC QQ' = TINVT * QQ * TINV (TINV is inverse of transposed TMAT) 49C 50 subroutine qsetup 51* 52 COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT, 53 & TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT 54 & ,RAFTER, TAFTER 55 REAL XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT 56 REAL TMAT(4,4), TINV(4,4), TINVT(4,4) 57 REAL SROT(4,4), SRTINV(4,4), SRTINVT(4,4) 58 REAL RAFTER(4,4), TAFTER(3) 59* 60 COMMON /ASSCOM/ noise, verbose 61 integer noise 62 logical verbose 63* 64 real TRAN(4,4), POST(4,4), det 65* 66* TMAT is a post-multiplier matrix, but unfortunately the 67* quadric surface math was worked out for pre-multipliers 68* So the quadric math uses the transpose of TMAT 69 call trnsp4( TRAN, TMAT ) 70* 71* Remove translation component from TMAT. 72* I had to make a choice whether the X,Y,Z "center" of the quadric 73* is implicitly coded into the coefficients, or whether it is to be 74* maintained explicitly. Since all other object types do the latter, 75* I have chosen to keep all translation components out of the quadric 76* coefficients 77 tran(1,4) = 0.0 78 tran(2,4) = 0.0 79 tran(3,4) = 0.0 80* 81* July 1999 - Allow post-hoc rotation matrix also 82 call tmul4( POST, RAFTER, TRAN ) 83* 84 det = tinv4( TINV, POST ) / TMAT(4,4) 85 call trnsp4( TINVT, TINV ) 86* 87* While we're at it, check for legality of rotation matrix 88 if (abs(1. - abs(det)) .gt. 0.02) then 89 write (noise,901) abs(det) 90 901 format('>>> Warning: Determinant of rotation matrix =', 91 & F7.3,' <<<') 92 endif 93 if (TMAT(1,4).ne.0.or.TMAT(2,4).ne.0.or.TMAT(3,4).ne.0) then 94 write (noise,902) 95 902 format('>>> Warning: Non-zero cross terms in 4th column of', 96 & ' TMAT <<<') 97 endif 98* 99* Same thing again for shadow rotation matrix 100 call trnsp4( TRAN, SROT ) 101 det = tinv4( SRTINV, TRAN ) 102 call trnsp4( SRTINVT, SRTINV ) 103 if (abs(1. - abs(det)) .gt. 0.02) then 104 write (noise,903) abs(det) 105 903 format('>>> Warning: Determinant of shadow matrix =', 106 & F7.3,' <<<') 107 endif 108* 109 return 110 end 111 112 113************************************************************************ 114* This routine does the exact test for pixel impingement of a quadric * 115* in both pixel and shadow space. * 116* Find Z coordinate of point on quadric surface with given X, Y * 117* Also return surface normal at that point * 118************************************************************************ 119CCC 120CC 121C 122 function qtest( center, coeffs, xp, yp, zp, qnorm, shadowspace, 123 & backside ) 124 IMPLICIT NONE 125 logical qtest, shadowspace 126 logical backside 127 real center(3), coeffs(10), xp, yp, zp, qnorm(3) 128* 129 COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT, 130 & TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT 131 & ,RAFTER, TAFTER 132 REAL XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT 133 REAL TMAT(4,4), TINV(4,4), TINVT(4,4) 134 REAL SROT(4,4), SRTINV(4,4), SRTINVT(4,4) 135 REAL RAFTER(4,4), TAFTER(3) 136* 137 real QA,QB,QC,QD,QE,QF,QG,QH,QI,QJ 138 real AA,BB,CC,DD,EE 139 real x,y,z 140* 141 QA = coeffs(1) 142 QB = coeffs(2) 143 QC = coeffs(3) 144 QD = coeffs(4) 145 QE = coeffs(5) 146 QF = coeffs(6) 147 QG = coeffs(7) 148 QH = coeffs(8) 149 QI = coeffs(9) 150 QJ = coeffs(10) 151* 152* xp and yp are in pixel coordinates, rather than in normalized ones 153* First displace center of quadric surface from the implicit origin 154* to the point specified by center(3), then convert to pixel coords 155* 156 x = (xp - center(1)) / scale 157 y = (yp - center(2)) / scale 158* 159 AA = QC 160 BB = 2.0 * (QF*x + QE*y + QI) 161 CC = QA*x*x + QB*y*y + 2.0*(QD*x*y + QG*x + QH*y) + QJ 162 DD = BB*BB - 4*AA*CC 163 if (DD.LT.0) then 164 qtest = .false. 165 return 166 else 167 qtest = .true. 168 EE = sqrt(DD) 169 endif 170 if (AA .eq. 0.) then 171CDEBUG Can this happen???? 172 z = 9999. 173 else if (backside) then 174 if (AA .le. 0.) then 175 z = (-BB + EE) / (2*AA) 176 else 177 z = (-BB - EE) / (2*AA) 178 endif 179 else 180 if (AA .gt. 0.) then 181 z = (-BB + EE) / (2*AA) 182 else 183 z = (-BB - EE) / (2*AA) 184 endif 185 endif 186 zp = z * scale 187 zp = zp + center(3) 188* 189* Surface normal comes from partial derivatives at this point 190 if (.not. shadowspace) then 191 qnorm(1) = QA*x + QD*y + QF*z + QG 192 qnorm(2) = QB*y + QD*x + QE*z + QH 193 qnorm(3) = QC*z + QF*x + QE*y + QI 194 endif 195* 196 return 197 end 198 199CCC Convert ANISOU description of anisotropic displacement parameters 200CC into quadric surface enclosing a given probability volume 201C Returns -1 if the Uij are non-positive definite, 1 otherwise 202C 203 function anitoquad( ANISOU, PROB, QUADRIC, EIGENS, EVECS) 204 implicit NONE 205 integer anitoquad 206 real ANISOU(6), PROB, QUADRIC(10) 207 real EIGENS(3), EVECS(4,4) 208c 209 COMMON /ASSCOM/ noise, verbose 210 integer noise 211 logical verbose 212c 213 real UU(4,4), UINV(4,4) 214 integer i 215c 216 real tinv4, det 217c 218c Save FPU traps later by checking that ANISOU is non-zero 219 do i = 1,3 220 EIGENS(i) = 0.0 221 end do 222 if (ANISOU(1).eq.0.and.ANISOU(2).eq.0.and.ANISOU(3).eq.0 .and. 223 & ANISOU(4).eq.0.and.ANISOU(5).eq.0.and.ANISOU(6).eq.0) then 224 anitoquad = -2 225 return 226 end if 227c 228c Build matrix from Uij coefficients and invert it 229 UU(1,1) = ANISOU(1) 230 UU(2,2) = ANISOU(2) 231 UU(3,3) = ANISOU(3) 232 UU(4,4) = -(1/PROB**2) 233 UU(1,2) = ANISOU(4) 234 UU(2,1) = ANISOU(4) 235 UU(2,3) = ANISOU(6) 236 UU(3,2) = ANISOU(6) 237 UU(1,3) = ANISOU(5) 238 UU(3,1) = ANISOU(5) 239 UU(1,4) = 0.0 240 UU(4,1) = 0.0 241 UU(2,4) = 0.0 242 UU(4,2) = 0.0 243 UU(3,4) = 0.0 244 UU(4,3) = 0.0 245 det = tinv4( UINV, UU ) 246c 247c and from that we extract the coefficients of the surface 248 QUADRIC(1) = UINV(1,1) 249 QUADRIC(2) = UINV(2,2) 250 QUADRIC(3) = UINV(3,3) 251 QUADRIC(4) = UINV(1,2) 252 QUADRIC(5) = UINV(2,3) 253 QUADRIC(6) = UINV(1,3) 254 QUADRIC(7) = UINV(1,4) 255 QUADRIC(8) = UINV(2,4) 256 QUADRIC(9) = UINV(3,4) 257 QUADRIC(10) = UINV(4,4) 258c 259c Find eigenvalues of the ellipsoid 260 call jacobi( UU, 3, 4, EIGENS, EVECS ) 261c 262c Units of this matrix are A**2; we want values in A 263 anitoquad = 1 264 do i = 1,3 265 if (EIGENS(i).gt.0.) then 266 EIGENS(i) = sqrt(EIGENS(i)) 267 else 268C write (noise,*) 'Non-positive definite ellipsoid!' 269 EIGENS(i) = 0. 270 anitoquad = -1 271 endif 272 enddo 273c 274 return 275 end 276 277CCC Matrix manipulation routines for 4x4 homogeneous transformations 278CC 279C 280 subroutine tmul4( C, A, B ) 281 real C(4,4), A(4,4), B(4,4) 282 integer i,j 283 do i = 1,4 284 do j = 1,4 285 C(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) 286 & + A(i,3)*B(3,j) + A(i,4)*B(4,j) 287 enddo 288 enddo 289 return 290 end 291 292 subroutine trnsp4( B, A ) 293 real B(4,4), A(4,4) 294 integer i,j 295 do i=1,4 296 do j=1,4 297 B(i,j) = A(j,i) 298 enddo 299 enddo 300 return 301 end 302 303CCC Matrix inversion for a 4x4 matrix A 304CC 305C 306 function tinv4( AI, A ) 307 real tinv4 308 real AI(4,4), A(4,4) 309c 310 real TMP(4,4), D 311 integer index(4) 312c 313 do i=1,4 314 do j=1,4 315 TMP(i,j) = A(i,j) 316 AI(i,j) = 0. 317 enddo 318 AI(i,i) = 1. 319 enddo 320 call ludcmp( TMP, 4, index, D ) 321 tinv4 = D * TMP(1,1)*TMP(2,2)*TMP(3,3)*TMP(4,4) 322 do j=1,4 323 call lubksb( TMP, 4, index, AI(1,j) ) 324 enddo 325 return 326 end 327 328************************************************************************ 329* Matrix inversion via LU decomposition * 330* adapted from Numerical Recipes in Fortran (1986) * 331************************************************************************ 332* 333CCC input NxN matrix A is replaced by its LU decomposition 334CC output index(N) records row permutation due to pivoting 335C output D is parity of row permutations 336c 337 subroutine ludcmp( A, n, index, D ) 338 implicit NONE 339 integer n,index(n) 340 real A(n,n) 341 real D 342c 343 integer i,imax,j,k 344 real aamax,dum,sum 345 integer NMAX 346 parameter (NMAX=10) 347 real vv(NMAX) 348c 349 d = 1. 350 do i=1,n 351 aamax = 0. 352 do j=1,n 353 if (abs(A(i,j)).gt.aamax) aamax = abs(A(i,j)) 354 enddo 355 call assert(aamax.ne.0.,'Singular matrix') 356 vv(i) = 1. / aamax 357 enddo 358c 359 imax = 1 360 do j=1,n 361 if (j.gt.1) then 362 do i=1,j-1 363 sum = A(i,j) 364 if (i.gt.1) then 365 do k=1,i-1 366 sum = sum - A(i,k)*A(k,j) 367 enddo 368 A(i,j) = sum 369 endif 370 enddo 371 endif 372 aamax = 0. 373 do i=j,n 374 sum = A(i,j) 375 if (j.gt.1) then 376 do k=1,j-1 377 sum = sum - A(i,k)*A(k,j) 378 enddo 379 A(i,j) = sum 380 endif 381 dum = vv(i) * abs(sum) 382 if (dum.ge.aamax) then 383 imax = i 384 aamax = dum 385 endif 386 enddo 387 if (j.ne.imax) then 388 do k=1,n 389 dum = A(imax,k) 390 A(imax,k) = A(j,k) 391 A(j,k) = dum 392 enddo 393 d = -d 394 vv(imax) = vv(j) 395 endif 396 index(j) = imax 397 call assert(A(j,j).ne.0.,'Singular matrix') 398 if (j.ne.n) then 399 dum = 1. / A(j,j) 400 do i=j+1,n 401 A(i,j) = A(i,j) * dum 402 enddo 403 endif 404 enddo 405 return 406 end 407 408CCC corresponding back-substitution routine 409CC 410C 411 subroutine lubksb( A, N, index, B ) 412 implicit NONE 413 integer n, index(n) 414 real A(n,n), B(n) 415c 416 integer i,ii,j,ll 417 real sum 418c 419 ii = 0 420 do i=1,n 421 ll = index(i) 422 sum = B(ll) 423 B(ll) = B(i) 424 if (ii.ne.0) then 425 do j=ii,i-1 426 sum = sum - A(i,j)*B(j) 427 enddo 428 else if (sum.ne.0.) then 429 ii = i 430 endif 431 B(i) = sum 432 enddo 433c 434 do i=n,1,-1 435 sum = B(i) 436 if (i.lt.n) then 437 do j=i+1,n 438 sum = sum - A(i,j)*B(j) 439 enddo 440 endif 441 B(i) = sum / A(i,i) 442 enddo 443 return 444 end 445 446************************************************************************ 447* Find eigenvalues and eigenvectors of nxn symmetric matrix using * 448* cyclic Jacobi method. NP is (Fortran) physical storage dimension. * 449* On return A is destroyed, D contains eigenvalues,and each column of * 450* V is a normalized eigenvector * 451* Adapted from Numerical Recipes in Fortran (1986) * 452* We only need it for 3x3 symmetric matrices, so it's overkill. * 453************************************************************************ 454CCC 455CC 456C 457 subroutine jacobi( A, n, np, D, V ) 458 PARAMETER (NMAX=4) 459c Machine dependent! (converge when off-diagonal sum is less than this) 460 PARAMETER (TINY = 1.e-37) 461 real A(np,np), D(n), V(np,np) 462 real B(NMAX), Z(NMAX) 463c 464 call assert(n.le.NMAX,'Matrix too big for eigenvector routine') 465c 466c Initialize V to identity matrix, B and D to diagonal of A 467 do ip = 1,n 468 do iq = 1,n 469 V(ip,iq) = 0.0 470 enddo 471 V(ip,ip) = 1.0 472 B(ip) = A(ip,ip) 473 D(ip) = B(ip) 474 Z(ip) = 0.0 475 enddo 476 nrot = 0 477c 478c 50 sweeps is never expected to happen; Press et al (1986) claim 479c 6 to 10 are typical for moderate matrix sizes 480c Empirical trials of rastep show 6 sweeps, 8-10 rotations 481 do i = 1,50 482 sm = 0.0 483 do ip = 1,n-1 484 do iq = ip+1,n 485 sm = sm+abs(A(ip,iq)) 486 enddo 487 enddo 488 if (sm .lt. TINY) return 489c After 4 sweeps skip rotation if the off-diagonal is small 490 if (i .lt. 4) then 491 thresh = 0.2*sm/(n*n) 492 else 493 thresh = 0.0 494 endif 495 do ip = 1,n-1 496 do iq = ip+1,n 497 g = 100. * abs(A(ip,iq)) 498 if ((i.gt.4) .and. 499 & (abs(D(ip)) + g .eq. abs(D(ip))) .and. 500 & (abs(D(iq)) + g .eq. abs(D(iq)))) then 501 A(ip,iq) = 0.0 502 else if (abs(A(ip,iq)).gt.thresh) then 503 h = D(iq) - D(ip) 504 if (abs(h) + g .eq. abs(h)) then 505 t = A(ip,iq) / h 506 else 507 theta = 0.5 * h / A(ip,iq) 508 t = 1. / (abs(theta) + sqrt(1.+theta*theta)) 509 if (theta.lt.0.) t = -t 510 endif 511 c = 1./sqrt(1.+t*t) 512 s = t * c 513 tau = s/(1.+c) 514 h = t * A(ip,iq) 515 Z(ip) = Z(ip) - h 516 Z(iq) = Z(iq) + h 517 D(ip) = D(ip) - h 518 D(iq) = D(iq) + h 519 A(ip,iq) = 0.0 520 do j = 1,ip-1 521 g = A(j,ip) 522 h = A(j,iq) 523 A(j,ip) = g - s*(h+g*tau) 524 A(j,iq) = h + s*(g-h*tau) 525 enddo 526 do j = ip+1,iq-1 527 g = A(ip,j) 528 h = A(j, iq) 529 A(ip,j) = g - s*(h+g*tau) 530 A(j,iq) = h + s*(g-h*tau) 531 enddo 532 do j = iq+1,n 533 g = A(ip,j) 534 h = A(iq,j) 535 A(ip,j) = g - s*(h+g*tau) 536 A(iq,j) = h + s*(g-h*tau) 537 enddo 538 do j = 1,n 539 g = V(j,ip) 540 h = V(j,iq) 541 V(j,ip) = g - s*(h+g*tau) 542 V(j,iq) = h + s*(g-h*tau) 543 enddo 544 nrot = nrot + 1 545 endif 546 enddo 547 enddo 548c 549c Update D with sum and reinitialize Z 550 do ip = 1,n 551 B(ip) = B(ip) + Z(ip) 552 D(ip) = B(ip) 553 Z(ip) = 0.0 554 enddo 555 enddo 556c 557 call assert(.false.,'Failed to find eigenvectors') 558 return 559 end 560c 561 562 563C This one really has nothing to do with quadrics per se, 564C but since it's invoked by qinp, both render and rastep 565C need to be able to see it. 566C Should really be in separate file of support routines 567C 568 FUNCTION PERSP( Z ) 569 REAL PERSP, Z 570 COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT, 571 & TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT 572 & ,RAFTER, TAFTER 573 REAL XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT 574 REAL TMAT(4,4), TINV(4,4), TINVT(4,4) 575 REAL SROT(4,4), SRTINV(4,4), SRTINVT(4,4) 576 REAL RAFTER(4,4), TAFTER(3) 577 IF (Z/EYEPOS .GT. 0.999) THEN 578 PERSP = 1000. 579 ELSE 580 PERSP = 1. / (1. - Z/EYEPOS) 581 ENDIF 582 RETURN 583 END 584 585