C ++********************************************************************
C                                                                      *
C   BCKCQ.F              SPEEDED UP  FEB 2000 ARDEAN LEITH             *
C                                                                      *
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
C=*                                                                    *
C=* This program is free software; you can redistribute it and/or      *
C=* modify it under the terms of the GNU General Public License as     *
C=* published by the Free Software Foundation; either version 2 of the *
C=* License, or (at your option) any later version.                    *
C=*                                                                    *
C=* This program is distributed in the hope that it will be useful,    *
C=* but WITHOUT ANY WARRANTY; without even the implied warranty of     *
C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C                                                                      *
C  BCKCQ(CUBE,LTC,DM,B,N,IPCUBE,NN)                                    *
C                                                                      *
C  PURPOSE:         BACKPROJECTS B INTO CUBE                           *
C                                                                      *
C  PARAMETERS:      CUBE                                     SENT/RET. *
C                   LTC      DIMENSIONS OF CUBE ARRAY        SENT      *
C                   DM                                                 *
C                   B                                        SENT      *
C                   N        DIMENSIONS OF B                 SENT      *
C                   IPCUBE                                   SENT      *
C                   NN       2ND DIMENSION OF IPCUBE         SENT      *
C                                                                      *
C IMAGE_PROCESSING_ROUTINE                                             *
C                                                                      *
C        0         2         3         4         5         6         7 *
C23456789012345678901234567890123456789012345678901234567890123456789012
C***********************************************************************

        SUBROUTINE  BCKCQ(CUBE,LTC,DM,B,N,IPCUBE,NN)

        DIMENSION     DM(9),CUBE(LTC),B(4,N,N)
        INTEGER       IPCUBE(5,NN)
        COMMON /PAR/  LDP,NM,LDPNM

#ifdef SP_MP
c$omp   parallel do private(i,j,xb,yb,xbb,ybb,iqx,iqy,jt)
        DO   I=1,NN
           XB = ((IPCUBE(3,I)-LDP)*DM(1)+(IPCUBE(4,I)-LDP)*DM(2)+
     &          (IPCUBE(5,I)-LDP)*DM(3)) + LDPNM
           YB = ((IPCUBE(3,I)-LDP)*DM(4)+(IPCUBE(4,I)-LDP)*DM(5)+
     &          (IPCUBE(5,I)-LDP)*DM(6)) + LDPNM

           DO    J=IPCUBE(1,I),IPCUBE(2,I)
	      JT=J-IPCUBE(1,I)
	      XBB=XB+JT*DM(1)
	      YBB=YB+JT*DM(4)
              IQX  = IFIX(XBB)
              IQY  = IFIX(YBB)

              CUBE(J) = CUBE(J) + 
     &                      B(1,IQX,IQY) + (YBB-IQY) * B(2,IQX,IQY) + 
     &          (XBB-IQX) * (B(3,IQX,IQY) + (YBB-IQY) * B(4,IQX,IQY) )

           ENDDO
        ENDDO
#else
        DO   I=1,NN
           XB = ((IPCUBE(3,I)-LDP)*DM(1)+(IPCUBE(4,I)-LDP)*DM(2)+
     &          (IPCUBE(5,I)-LDP)*DM(3)) + LDPNM
           YB = ((IPCUBE(3,I)-LDP)*DM(4)+(IPCUBE(4,I)-LDP)*DM(5)+
     &          (IPCUBE(5,I)-LDP)*DM(6)) + LDPNM

           DO    J=IPCUBE(1,I),IPCUBE(2,I)

              IQX  = IFIX(XB)
              IQY  = IFIX(YB)

              CUBE(J) = CUBE(J) + 
     &                      B(1,IQX,IQY) + (YB-IQY) * B(2,IQX,IQY) + 
     &          (XB-IQX) * (B(3,IQX,IQY) + (YB-IQY) * B(4,IQX,IQY) )

              XB = XB + DM(1)
              YB = YB + DM(4)
           ENDDO
        ENDDO
#endif
        END



#ifdef NEVER
              XBB  = (J-IPCUBE(1,I)) * DM1 + XB
              YBB  = (J-IPCUBE(1,I)) * DM4 + YB
              IQX  = IFIX(XBB+FLOAT(LDPNM))
              IQY  = IFIX(YBB+FLOAT(LDPNM))
              DIPX = XBB+LDPNM-IQX
              DIPY = YBB+LDPNM-IQY
              IF (iqx .ne. iqx1 .or. iqy .ne. iqy1 .or.
     &            dipx1 .ne. dipx .or. dipy .ne. dipy1) then
                  write(6,*) 'simailar test failed'
                  write(6,*) iqx,iqy, ' != ',iqx1,iqy1
                  write(6,*) dipx,dipy, ' != ',dipx1,dipy1
                  stop 'simalar test failed'
              endif


C             REFORMATTED FASTER VERSION :
              CUBE(J) = CUBE(J) +
     &               B(IQX,  IQY) + 
     &         DIPY*(B(IQX,  IQY+1) - B(IQX,  IQY)) + 
     &         DIPX*(B(IQX+1,IQY)   - B(IQX,  IQY)  +
     &         DIPY*(B(IQX+1,IQY+1) - B(IQX+1,IQY)  -
     &              (B(IQX,  IQY+1) - B(IQX,  IQY))  ))

C             FASTER VERSION :
              CUBE(J)=CUBE(J)
     &        +B(IQX,IQY)+DIPY*(B(IQX,IQY+1)-B(IQX,IQY))
     &        +DIPX*(B(IQX+1,IQY)-B(IQX,IQY)
     &        +DIPY*(B(IQX+1,IQY+1)-B(IQX+1,IQY)
     &        -B(IQX,IQY+1)+B(IQX,IQY)))

C             SLOWER VERSION :
C     &                 +(1.0-DIPX)*(1.0-DIPY)*B(MAP(IQX,IQY))
C     &                 +     DIPX *(1.0-DIPY)*B(MAP(IQX+1,IQY))
C     &                 +(1.0-DIPX)*     DIPY *B(MAP(IQX,IQY+1))
C     &                 +     DIPX *     DIPY *B(MAP(IQX+1,IQY+1))

#endif