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 solution yet. Break 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* - fix optimization problems with qinp 45************************************************************************ 46 47CCC Process single object descriptor during input phase 48CC 49C 50 function qinp( buf, detail, shadow, sdtail ) 51* Use MODULE LISTS not COMMON /LISTS/ for dynamic allocation 52 USE LISTS 53* 54 IMPLICIT NONE 55 logical qinp 56 real buf(100) 57 real detail(17), sdtail(14) 58 logical shadow 59* 60 real QQ(4,4), QP(4,4), QT(4,4) 61 real xq, yq, zq, radlim, red, grn, blu 62 real xc, yc, zc, rc 63 real xr, yr, zr, xs, ys, zs, rs 64 real pfac 65* 66 integer ix,iy,ixlo,ixhi,iylo,iyhi 67c VOLATILE ix,iy,ixlo,ixhi,iylo,iyhi 68* 69* Array sizes 70 INCLUDE 'parameters.incl' 71* 72 EXTERNAL PERSP 73 REAL PERSP 74* 75 COMMON /RASTER/ NTX,NTY,NPX,NPY 76 INTEGER NTX,NTY,NPX,NPY 77* 78 COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT, 79 & TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT 80 & ,RAFTER, TAFTER 81 REAL XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT 82 REAL TMAT(4,4), TINV(4,4), TINVT(4,4) 83 REAL SROT(4,4), SRTINV(4,4), SRTINVT(4,4) 84 REAL RAFTER(4,4), TAFTER(3) 85* 86C Replace COMMON with MODULE 87C COMMON /LISTS/ KOUNT, MOUNT, TTRANS, ISTRANS 88C INTEGER KOUNT(MAXNTX,MAXNTY), MOUNT(NSX,NSY) 89C INTEGER TTRANS(MAXNTX,MAXNTY), ISTRANS 90* 91 COMMON /NICETIES/ TRULIM, ZLIM, FRONTCLIP, BACKCLIP 92 & , ISOLATION 93 REAL TRULIM(3,2), ZLIM(2), FRONTCLIP, BACKCLIP 94 LOGICAL ISOLATION 95* 96* Assume this is legitimate 97 qinp = .TRUE. 98* 99* Update limits (have to trust the center coords) 100* 101 xq = buf(1) 102 yq = buf(2) 103 zq = buf(3) 104 radlim = buf(4) 105 call assert(radlim.ge.0,'limiting radius < 0 in quadric') 106 if (radlim.gt.0) then 107 trulim(1,1) = MIN( trulim(1,1), xq) 108 trulim(1,2) = MAX( trulim(1,2), xq) 109 trulim(2,1) = MIN( trulim(2,1), yq) 110 trulim(2,2) = MAX( trulim(2,2), yq) 111 trulim(3,1) = MIN( trulim(3,1), zq) 112 trulim(3,2) = MAX( trulim(3,2), zq) 113 endif 114* 115* Standard color checks 116* 117 red = buf(5) 118 grn = buf(6) 119 blu = buf(7) 120 call assert(red.ge.0,'red < 0 in quadric') 121 call assert(red.le.1,'red > 1 in quadric') 122 call assert(grn.ge.0,'grn < 0 in quadric') 123 call assert(grn.le.1,'grn > 1 in quadric') 124 call assert(blu.ge.0,'blu < 0 in quadric') 125 call assert(blu.le.1,'blu > 1 in quadric') 126* 127* Transform center before saving 128* (But can we deal with perspective????) 129* 130 call transf (xq, yq, zq) 131 radlim = radlim / TMAT(4,4) 132 if (eyepos.gt.0) then 133c pfac = 1./(1.-zq/eyepos) 134 pfac = persp(zq) 135 xq = xq * pfac 136 yq = yq * pfac 137 zq = zq * pfac 138 radlim = radlim * pfac 139 endif 140 xc = xq * scale + xcent 141 yc = yq * scale + ycent 142 zc = zq * scale 143 rc = radlim * scale 144* save transformed Z limits 145 zlim(1) = min( zlim(1), zc ) 146 zlim(2) = max( zlim(2), zc ) 147* 148* check for Z-clipping 149 if (zc.gt.FRONTCLIP .or. zc.lt.BACKCLIP) then 150 qinp = .FALSE. 151 return 152 endif 153* 154 detail(1) = xc 155 detail(2) = yc 156 detail(3) = zc 157 detail(4) = rc 158 detail(5) = red 159 detail(6) = grn 160 detail(7) = blu 161* 162* This is a terrible kludge, but necessary if called from normal3d -size BIGxBIG 163* Thes test should really be if we are called from normal3d but no flag for that 164 if (ntx.gt.size(kount,1) .or. nty.gt.size(kount,2)) goto 101 165* 166* Tally for tiles the object might impinge on 167* Again we are relying on the correctness of the center coordinates 168* 169 ixlo = INT(xc-rc) / npx + 1 170 ixhi = INT(xc+rc) / npx + 1 171 iylo = INT(yc-rc) / npy + 1 172 iyhi = INT(yc+rc) / npy + 1 173 if (ixlo.lt.1) ixlo = 1 174 if (ixlo.gt.NTX) goto 101 175 if (ixhi.lt.1) goto 101 176 if (ixhi.gt.NTX) ixhi = NTX 177 if (iylo.lt.1) iylo = 1 178 if (iylo.gt.NTY) goto 101 179 if (iyhi.lt.1) goto 101 180 if (iyhi.gt.NTY) iyhi = NTY 181 do iy = iylo,iyhi 182 do ix = ixlo,ixhi 183 KOUNT(ix,iy) = KOUNT(ix,iy) + 1 184 TTRANS(ix,iy) = TTRANS(ix,iy) + istrans 185 enddo 186 enddo 187101 continue 188* 189* build matrix from coeffients describing quadric surface in standard form 190* 191 qq(1,1) = buf(8) 192 qq(2,2) = buf(9) 193 qq(3,3) = buf(10) 194 qq(1,2) = buf(11) 195 qq(2,1) = buf(11) 196 qq(2,3) = buf(12) 197 qq(3,2) = buf(12) 198 qq(1,3) = buf(13) 199 qq(3,1) = buf(13) 200 qq(1,4) = buf(14) 201 qq(4,1) = buf(14) 202 qq(2,4) = buf(15) 203 qq(4,2) = buf(15) 204 qq(3,4) = buf(16) 205 qq(4,3) = buf(16) 206 qq(4,4) = buf(17) 207* 208* Transformed matrix QP = TINV(Transpose) * QQ * TINV 209* where TINV is the inverse of TMAT 210* 211 call tmul4( qt, qq, tinv ) 212 call tmul4( qp, tinvt, qt ) 213CD noise = 0 214CD write (noise,191) 'QT ',((QT(i,j),j=1,4),i=1,4) 215CD 191 format(a,4(/,4f8.4)) 216CD write (noise,191) 'QP ',((QP(i,j),j=1,4),i=1,4) 217* 218* Save components of quadric surface built from transformed matrix 219* for use during rendering 220* 221 detail(8) = qp(1,1) 222 detail(9) = qp(2,2) 223 detail(10) = qp(3,3) 224 detail(11) = qp(1,2) 225 detail(12) = qp(2,3) 226 detail(13) = qp(1,3) 227 detail(14) = qp(1,4) 228 detail(15) = qp(2,4) 229 detail(16) = qp(3,4) 230 detail(17) = qp(4,4) 231* 232* Do it all over again for the shadow buffers NOT TESTED YET! 233* (since I'm a little confused about what transformation 234* I need to apply to the QQ matrix in shadow space) 235* 236 if (.not.shadow) return 237* first transform center and limiting sphere 238 xr = srot(1,1)*xq + srot(1,2)*yq + srot(1,3)*zq 239 yr = srot(2,1)*xq + srot(2,2)*yq + srot(2,3)*zq 240 zr = srot(3,1)*xq + srot(3,2)*yq + srot(3,3)*zq 241 xs = xr * scale + sxcent 242 ys = yr * scale + sycent 243 zs = zr * scale 244 rs = radlim * scale 245 sdtail(1) = xs 246 sdtail(2) = ys 247 sdtail(3) = zs 248 sdtail(4) = rs 249* tally shadow tiles the object might impinge on 250 ixlo = INT(xs-rs) / npx + 1 251 ixhi = INT(xs+rs) / npx + 1 252 iylo = INT(ys-rs) / npy + 1 253 iyhi = INT(ys+rs) / npy + 1 254 if (ixlo.lt.1) ixlo = 1 255 if (ixlo.gt.NSX) goto 209 256 if (ixhi.lt.1) goto 209 257 if (ixhi.gt.NSX) ixhi = NSX 258 if (iylo.lt.1) iylo = 1 259 if (iylo.gt.NSY) goto 209 260 if (iyhi.lt.1) goto 209 261 if (iyhi.gt.NSY) iyhi = NSY 262 do iy = iylo,iyhi 263 do ix = ixlo,ixhi 264 MOUNT(ix,iy) = MOUNT(ix,iy) + 1 265 enddo 266 enddo 267* transform QQ into shadow space as well 268 call tmul4( qt, qp, srtinv ) 269 call tmul4( qp, srtinvt, qt ) 270 sdtail(5) = qp(1,1) 271 sdtail(6) = qp(2,2) 272 sdtail(7) = qp(3,3) 273 sdtail(8) = qp(1,2) 274 sdtail(9) = qp(2,3) 275 sdtail(10) = qp(1,3) 276 sdtail(11) = qp(1,4) 277 sdtail(12) = qp(2,4) 278 sdtail(13) = qp(3,4) 279 sdtail(14) = qp(4,4) 280* 281 209 continue 282* 283* Done with input 284* 285 return 286 end 287