C++********************************************************************* C C WPRO.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 WPRO(B,NSAM,NROW,CUBE,LTB,IPCUBE,NN,PHI,THETA,PSI,RI) C C IPCUBE: 1 - BEGINNING C 2 - END C 3 - IX C 4 - IY C 5 - IZ C C IMAGE_PROCESSING_ROUTINE C C 1 2 3 4 5 6 7 C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE WPRO & (B,NSAM,NROW,CUBE,LTB,IPCUBE,NN,PHI,THETA,PSI,RI) DIMENSION B(NSAM*NROW),CUBE(LTB) INTEGER IPCUBE(5,NN) DIMENSION DM(6) DOUBLE PRECISION CPHI,SPHI,CTHE,STHE,CPSI,SPSI COMMON /PAR/ LDPX,LDPY,LDPZ DOUBLE PRECISION QUADPI,DGR_TO_RAD,RAD_TO_DGR PARAMETER (QUADPI = 3.141592653589793238462643383279502884197) PARAMETER (DGR_TO_RAD = (QUADPI/180)) R = RI * RI CPHI = DCOS(DBLE(PHI)*DGR_TO_RAD) SPHI = DSIN(DBLE(PHI)*DGR_TO_RAD) CTHE = DCOS(DBLE(THETA)*DGR_TO_RAD) STHE = DSIN(DBLE(THETA)*DGR_TO_RAD) CPSI = DCOS(DBLE(PSI)*DGR_TO_RAD) SPSI = DSIN(DBLE(PSI)*DGR_TO_RAD) DM(1) = CPHI*CTHE*CPSI-SPHI*SPSI DM(2) = SPHI*CTHE*CPSI+CPHI*SPSI DM(3) = -STHE*CPSI DM(4) = -CPHI*CTHE*SPSI-SPHI*CPSI DM(5) = -SPHI*CTHE*SPSI+CPHI*CPSI DM(6) = STHE*SPSI C DM(7) = STHE*CPHI C DM(8) = STHE*SPHI C DM(9) = CTHE DM1 = DM(1) DM4 = DM(4) C ZERO THE WHOLE B ARRAY B = 0.0 DIM = MIN(NSAM,NROW) IF ((2*(RI+1) + 1) .GT. DIM) THEN C MUST CHECK IQX & IQY BOUNDARIES DO I=1,NN XB = (IPCUBE(3,I)-LDPX)*DM(1) +(IPCUBE(4,I)-LDPY)*DM(2) + & (IPCUBE(5,I)-LDPZ)*DM(3) + LDPX YB = (IPCUBE(3,I)-LDPX)*DM(4) +(IPCUBE(4,I)-LDPY)*DM(5) + & (IPCUBE(5,I)-LDPZ)*DM(6) + LDPY DO J=IPCUBE(1,I),IPCUBE(2,I) C CHECK FOR PIXELS OUT OF BOUNDS IQX = IFIX(XB) IF (IQX.LT.1 .OR. IQX.GE.NSAM) GOTO 2 IQY = IFIX(YB) IF (IQY.LT.1 .OR. IQY.GE.NROW) GOTO 2 CT = CUBE(J) DIPX = (XB - IQX) DIPY = (YB - IQY) * CT DIPY1M = (CT - DIPY) DIPX1M = (1.0 - DIPX) ILOC = (IQY - 1) * NSAM + IQX B(ILOC) = B(ILOC) + DIPX1M * DIPY1M B(ILOC+1) = B(ILOC+1) + DIPX * DIPY1M B(ILOC+NSAM) = B(ILOC+NSAM) + DIPX1M * DIPY B(ILOC+NSAM+1)= B(ILOC+NSAM+1) + DIPX * DIPY 2 XB = XB + DM1 YB = YB + DM4 ENDDO ENDDO ELSE C NO NEED TO CHECK IQX & IQY BOUNDARIES DO I=1,NN XB = (IPCUBE(3,I)-LDPX)*DM(1) +(IPCUBE(4,I)-LDPY)*DM(2) + & (IPCUBE(5,I)-LDPZ)*DM(3) + LDPX YB = (IPCUBE(3,I)-LDPX)*DM(4) +(IPCUBE(4,I)-LDPY)*DM(5) + & (IPCUBE(5,I)-LDPZ)*DM(6) + LDPY DO J=IPCUBE(1,I),IPCUBE(2,I) IQX = IFIX(XB) IQY = IFIX(YB) CT = CUBE(J) DIPX = (XB - IQX) DIPY = (YB - IQY) * CT DIPY1M = (CT - DIPY) DIPX1M = (1.0 - DIPX) ILOC = (IQY - 1) * NSAM + IQX B(ILOC) = B(ILOC) + DIPX1M * DIPY1M B(ILOC+1) = B(ILOC+1) + DIPX * DIPY1M B(ILOC+NSAM) = B(ILOC+NSAM) + DIPX1M * DIPY B(ILOC+NSAM+1)= B(ILOC+NSAM+1) + DIPX * DIPY XB = XB + DM1 YB = YB + DM4 ENDDO ENDDO ENDIF END #ifdef NEVER XB = (IPCUBE(3,I)-LDPX)*DM1 +(IPCUBE(4,I)-LDPY)*DM2 + & (IPCUBE(5,I)-LDPZ)*DM3 + LDPX YB = (IPCUBE(3,I)-LDPX)*DM4 +(IPCUBE(4,I)-LDPY)*DM5 + & (IPCUBE(5,I)-LDPZ)*DM6 + LDPY CT = CUBE(J) DIPX = (XB - IQX) DIPY = (YB - IQY) B(IQX,IQY) =B(IQX,IQY) +(1.0-DIPX)*(1.0-DIPY)*CT B(IQX+1,IQY) =B(IQX+1,IQY) + DIPX *(1.0-DIPY)*CT B(IQX,IQY+1) =B(IQX,IQY+1) +(1.0-DIPX)* DIPY *CT B(IQX+1,IQY+1)=B(IQX+1,IQY+1)+ DIPX * DIPY *CT c IQXP1 = IQX + 1 c IQYP1 = IQY + 1 c T1MXC = (1.0 - DIPX) * CT c T1MY = (1.0 - DIPY) c B(IQX, IQY) = B(IQX, IQY) + T1MXC * T1MY c B(IQXP1,IQY) = B(IQXP1,IQY) + DIPX * T1MY * CT c B(IQX, IQYP1) = B(IQX, IQYP1) + T1MXC * DIPY c B(IQXP1,IQYP1) = B(IQXP1,IQYP1) + DIPX * DIPY * CT C ORIGINAL (SLOWER) C B(IQX,IQY) =B(IQX,IQY) +(1.0-DIPX)*(1.0-DIPY)*CT C B(IQX+1,IQY) =B(IQX+1,IQY) + DIPX *(1.0-DIPY)*CT C B(IQX,IQY+1) =B(IQX,IQY+1) +(1.0-DIPX)* DIPY *CT C B(IQX+1,IQY+1)=B(IQX+1,IQY+1)+ DIPX * DIPY *CT #endif #ifdef NEVER igo = IPCUBE(1,I) iend = IPCUBE(2,I) igo = ifix(xb) iend = ifix(IPCUBE(2,I)) igo = min(igo,iend) iend = max(igo,iend) xbmax = (IPCUBE(2,I) - IPCUBE(1,I)) * DM1 ixbmax = IFIX(xbmax) if (ixbmax .gt. ybmax = (IPCUBE(2,I) - IPCUBE(1,I)) * DM4 iybmax = IFIX(ybmax) #endif