Logo Search packages:      
Sourcecode: raster3d version File versions  Download package

qinp.f

************************************************************************
*              Support routines for quadric surfaces                   *
************************************************************************
* EAM Jun 1997    - initial version, supports version 2.4(alpha) of render 
* EAM May 1998    - additional error checking to go with Parvati/rastep
* EAM Jan 1999    - version 2.4i
* EAM Mar 2008    - Gfortran optimization breaks the object accounting.
*             No solution yet.  Break qinp.f into separate file
************************************************************************
*
* Quadric surfaces include spheres, cones, ellipsoids, paraboloids, and
* hyperboloids.  The motivation for this code was to allow rendering 
* thermal ellipsoids for atoms, so the other shapes have not been
* extensively tested.
* A quadric surface is described by 10 parameters (A ... J).
* For efficiency during rendering it is also useful to know the center and 
* a bounding sphere. So a QUADRIC descriptor to render has 17 parameters:
* 14 (object type QUADRIC)
*    X Y Z RADLIM RED GRN BLU
*    A B C D E F G H I J
*
* The surface itself is the set of points for which Q(x,y,z) = 0
* where 
*      Q(x,y,z)   =  A*x^2 +  B*y^2 +  C*z^2 
*           + 2D*x*y + 2E*y*z + 2F*z*x
*           + 2G*x   + 2H*y   + 2I*z   
*           +  J
* 
* It is convenient to store this information in a matrix QQ
*     | QA QD QF QG |         QA = A    QB = B    QC = C
* QQ =      | QD QB QE QH |         QD = D    QE = E    QF = F  
*     | QF QE QC QI |         QG = G    QH = H    QI = I  
*     | QG QH QI QJ |         QJ = J
*
* Then Q(x,y,x) = XT*QQ*X     where X = (x,y,z,1)
* The point of this is that a 4x4 homogeneous transformation T can be 
* applied to QQ by matrix multiplication:     QQ' = TinvT * QQ * Tinv
* 
* The surface normal is easily found by taking the partial derivatives
* of Q(x,y,z) at the point of interest.
************************************************************************
* TO DO:
*     - can we distinguish an ellipsoid from other quadrics? Do we care?
*     - fix optimization problems with qinp
************************************************************************

CCC     Process single object descriptor during input phase
CC
C
      function qinp( buf, detail, shadow, sdtail )
*       Use MODULE LISTS not COMMON /LISTS/ for dynamic allocation
        USE LISTS
*
      IMPLICIT NONE
      logical  qinp
      real    buf(100)
      real    detail(17), sdtail(14)
      logical shadow
*
      real    QQ(4,4), QP(4,4), QT(4,4)
      real        xq, yq, zq, radlim, red, grn, blu
      real    xc, yc, zc, rc
      real    xr, yr, zr, xs, ys, zs, rs
      real  pfac
*
      integer  ix,iy,ixlo,ixhi,iylo,iyhi
c     VOLATILE ix,iy,ixlo,ixhi,iylo,iyhi
*
*     Array sizes
      INCLUDE 'parameters.incl'
*
      EXTERNAL PERSP
      REAL     PERSP
*
      COMMON /RASTER/ NTX,NTY,NPX,NPY
      INTEGER         NTX,NTY,NPX,NPY
*
      COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT,
     &                  TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT
     &                 ,RAFTER, TAFTER
      REAL   XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT
      REAL   TMAT(4,4), TINV(4,4),   TINVT(4,4) 
      REAL   SROT(4,4), SRTINV(4,4), SRTINVT(4,4)
      REAL   RAFTER(4,4), TAFTER(3)
*
C     Replace COMMON with MODULE 
C     COMMON /LISTS/ KOUNT, MOUNT, TTRANS, ISTRANS
C     INTEGER KOUNT(MAXNTX,MAXNTY), MOUNT(NSX,NSY)
C     INTEGER TTRANS(MAXNTX,MAXNTY), ISTRANS
*
      COMMON /NICETIES/ TRULIM,      ZLIM,    FRONTCLIP, BACKCLIP
     &                , ISOLATION
      REAL              TRULIM(3,2), ZLIM(2), FRONTCLIP, BACKCLIP
      LOGICAL           ISOLATION
*
* Assume this is legitimate
      qinp = .TRUE.
*
* Update limits (have to trust the center coords)
*
      xq     = buf(1)
      yq     = buf(2)
      zq     = buf(3)
      radlim = buf(4)
      call assert(radlim.ge.0,'limiting radius < 0 in quadric')
      if (radlim.gt.0) then
          trulim(1,1) = MIN( trulim(1,1), xq)
          trulim(1,2) = MAX( trulim(1,2), xq)
          trulim(2,1) = MIN( trulim(2,1), yq)
          trulim(2,2) = MAX( trulim(2,2), yq)
          trulim(3,1) = MIN( trulim(3,1), zq)
          trulim(3,2) = MAX( trulim(3,2), zq)
      endif
*
* Standard color checks
*
      red    = buf(5)
      grn    = buf(6)
      blu    = buf(7)
      call assert(red.ge.0,'red < 0 in quadric')
      call assert(red.le.1,'red > 1 in quadric')
      call assert(grn.ge.0,'grn < 0 in quadric')
      call assert(grn.le.1,'grn > 1 in quadric')
      call assert(blu.ge.0,'blu < 0 in quadric')
      call assert(blu.le.1,'blu > 1 in quadric')
*
* Transform center before saving
* (But can we deal with perspective????)
*
      call transf (xq, yq, zq)
      radlim = radlim / TMAT(4,4)
      if (eyepos.gt.0) then
c         pfac = 1./(1.-zq/eyepos)
          pfac = persp(zq)
          xq = xq * pfac
          yq = yq * pfac
          zq = zq * pfac
          radlim = radlim * pfac
      endif
      xc = xq * scale + xcent
      yc = yq * scale + ycent
      zc = zq * scale
      rc = radlim * scale
*     save transformed Z limits
      zlim(1) = min( zlim(1), zc )
      zlim(2) = max( zlim(2), zc )
*
*     check for Z-clipping
      if (zc.gt.FRONTCLIP .or. zc.lt.BACKCLIP) then
          qinp = .FALSE.
          return
      endif
*
      detail(1) = xc
      detail(2) = yc
      detail(3) = zc
      detail(4) = rc
      detail(5) = red
      detail(6) = grn
      detail(7) = blu
*
* This is a terrible kludge, but necessary if called from normal3d -size BIGxBIG
* Thes test should really be if we are called from normal3d but no flag for that
      if (ntx.gt.size(kount,1) .or. nty.gt.size(kount,2)) goto 101
*
* Tally for tiles the object might impinge on
* Again we are relying on the correctness of the center coordinates
*
      ixlo = (xc-rc) / npx + 1
      ixhi = (xc+rc) / npx + 1
      iylo = (yc-rc) / npy + 1
      iyhi = (yc+rc) / npy + 1
      if (ixlo.lt.1)   ixlo = 1
      if (ixlo.gt.NTX) goto 101
      if (ixhi.lt.1)   goto 101
      if (ixhi.gt.NTX) ixhi = NTX
      if (iylo.lt.1)   iylo = 1
      if (iylo.gt.NTY) goto 101
      if (iyhi.lt.1)   goto 101
      if (iyhi.gt.NTY) iyhi = NTY
      do iy = iylo,iyhi
      do ix = ixlo,ixhi
          KOUNT(ix,iy) = KOUNT(ix,iy) + 1
          TTRANS(ix,iy) = TTRANS(ix,iy) + istrans
      enddo
      enddo
101   continue
*
* build matrix from coeffients describing quadric surface in standard form
*
      qq(1,1) = buf(8)
      qq(2,2) = buf(9)
      qq(3,3) = buf(10)
      qq(1,2) = buf(11)
      qq(2,1) = buf(11)
      qq(2,3) = buf(12)
      qq(3,2) = buf(12)
      qq(1,3) = buf(13)
      qq(3,1) = buf(13)
      qq(1,4) = buf(14)
      qq(4,1) = buf(14)
      qq(2,4) = buf(15)
      qq(4,2) = buf(15)
      qq(3,4) = buf(16)
      qq(4,3) = buf(16)
      qq(4,4) = buf(17)
*
* Transformed matrix QP = TINV(Transpose) * QQ * TINV
* where TINV is the inverse of TMAT
*
      call tmul4( qt, qq,  tinv )
      call tmul4( qp, tinvt, qt )
CD    noise = 0
CD    write (noise,191) 'QT  ',((QT(i,j),j=1,4),i=1,4)
CD 191      format(a,4(/,4f8.4))
CD    write (noise,191) 'QP  ',((QP(i,j),j=1,4),i=1,4)
*
* Save components of quadric surface built from transformed matrix 
* for use during rendering
*
      detail(8)  = qp(1,1)
      detail(9)  = qp(2,2)
      detail(10) = qp(3,3)
      detail(11) = qp(1,2)
      detail(12) = qp(2,3)
      detail(13) = qp(1,3)
      detail(14) = qp(1,4)
      detail(15) = qp(2,4)
      detail(16) = qp(3,4)
      detail(17) = qp(4,4)
*
* Do it all over again for the shadow buffers NOT TESTED YET!  
* (since I'm a little confused about what transformation
* I need to apply to the QQ matrix in shadow space)
*
      if (.not.shadow) return
*     first transform center and limiting sphere
      xr = srot(1,1)*xq + srot(1,2)*yq + srot(1,3)*zq
      yr = srot(2,1)*xq + srot(2,2)*yq + srot(2,3)*zq
      zr = srot(3,1)*xq + srot(3,2)*yq + srot(3,3)*zq
      xs = xr * scale + sxcent
      ys = yr * scale + sycent
      zs = zr * scale
      rs = radlim * scale
      sdtail(1) = xs
      sdtail(2) = ys
      sdtail(3) = zs
      sdtail(4) = rs
*     tally shadow tiles the object might impinge on
      ixlo = (xs-rs) / npx + 1
      ixhi = (xs+rs) / npx + 1
      iylo = (ys-rs) / npy + 1
      iyhi = (ys+rs) / npy + 1
      if (ixlo.lt.1)   ixlo = 1
      if (ixlo.gt.NSX) goto 209
      if (ixhi.lt.1)   goto 209
      if (ixhi.gt.NSX) ixhi = NSX
      if (iylo.lt.1)   iylo = 1
      if (iylo.gt.NSY) goto 209
      if (iyhi.lt.1)   goto 209
      if (iyhi.gt.NSY) iyhi = NSY
      do iy = iylo,iyhi
      do ix = ixlo,ixhi
          MOUNT(ix,iy) = MOUNT(ix,iy) + 1
      enddo
      enddo
*     transform QQ into shadow space as well
      call tmul4( qt, qp, srtinv )
      call tmul4( qp, srtinvt, qt )
      sdtail(5)  = qp(1,1)
      sdtail(6)  = qp(2,2)
      sdtail(7)  = qp(3,3)
      sdtail(8)  = qp(1,2)
      sdtail(9)  = qp(2,3)
      sdtail(10) = qp(1,3)
      sdtail(11) = qp(1,4)
      sdtail(12) = qp(2,4)
      sdtail(13) = qp(3,4)
      sdtail(14) = qp(4,4)
*
 209  continue
*
* Done with input
*
      return
      end


Generated by  Doxygen 1.6.0   Back to index