C     Fx compiler version 2.1
CC    ALIGN workspace(I,J,K) WITH t(I)
CC    TEMPLATE t(50)
CC    DISTRIBUTE t(BLOCK_CYCLIC(2,32))
C     CALL getarg(1,filename)
C     blank = index(filename,' ')
C     WRITE (unit=6,fmt=*) 'Opening ',filename,' for input'
C     OPEN (unit=10,file=filename,form='unformatted',status='old')
C     CALL getarg(2,filename)
C     blank = index(filename,' ')
C     WRITE (unit=6,fmt=*) 'Opening ',filename,' for output'
C     INQUIRE (file=filename,exist=created)
C     IF (created) THEN
C     OPEN (unit=11,file=filename,form='unformatted',status='old')
C     ELSE
C     OPEN (unit=11,file=filename,form='unformatted',status='new')
C     ENDIF
C     READ (unit=10) dummy
C     nimgs = dummy
C     READ (unit=10) minx
C     READ (unit=10) miny
C     READ (unit=10) minz
C     READ (unit=10) maxx
C     READ (unit=10) maxy
C     READ (unit=10) maxz
C     PDO (d = 1:50)
C     POUT workspace(d,:,:)
C     LOCAL
C     IF (fxcellid.LT.32) THEN
C     fxloopstart0 = (MAX(((fxcellid * 2) + 1),1))
C     fxlmidx0 = ifxlm(fxworkspacedesc,0,fxloopstart0)
C     DO d = fxloopstart0, MIN(((fxcellid * 2) + 2),50), 1
C     DO i = 1, 50, 1
C     DO j = 1, 50, 1
C     workspace(d,j,i) = 0.
C     ENDDO
C     ENDDO
C     fxlmidx0 = (fxlmidx0 + 1)
C     ENDDO
C     ENDIF
C     ENDLOCAL
C     ENDPDO
C     total_time = 0.
C101  CONTINUE
C     READ (unit=10) cont
C     IF (cont.EQ.1.) THEN
C     DO i = 1, nimgs, 1
C     DO j = 1, 3, 1
C     READ (unit=10) (proj_xform(i,j,k),k=1,4)
C     ENDDO
C     ENDDO
C     READ (unit=10) (((img(i,k,j),k=1,528),j=1,486),i=1,nimgs)
C     CALL timer_start()
C     PDO (d = 0:49)
C     POUT workspace(d,:,:)
C     LOCAL
C     IF (fxcellid.LT.32) THEN
C     fxloopstart0 = (MAX(((fxcellid * 2) + 1),0))
C     fxlmidx0 = ifxlm(fxworkspacedesc,0,fxloopstart0)
C     DO d = fxloopstart0, MIN(((fxcellid * 2) + 2),49), 1
C     depth = (minz + (((d + (-1)) * (maxz - minz)) / 50))
C     DO j = 1, 50, 1
C     DO k = 1, 50, 1
C     sum(k,j) = 0.
C     sumsq(k,j) = 0.
C     ENDDO
C     ENDDO
C     DO i = 1, nimgs, 1
C     CALL warp(img,tempimg,i,minx,miny,maxx,maxy,depth)
C     DO j = 1, 50, 1
C     DO k = 1, 50, 1
C     sum(k,j) = (sum(k,j) + tempimg(k,j))
C     sumsq(k,j) = (sumsq(k,j) + (tempimg(k,j) ** 2))
C     ENDDO
C     ENDDO
C     ENDDO
C     CALL variance(sum,sumsq)
C     DO i = 1, 50, 1
C     DO j = 1, 50, 1
C     workspace(d,j,i) = (workspace(d,j,i) + (1 / sumsq(j,i)))
C     ENDDO
C     ENDDO
C     fxlmidx0 = (fxlmidx0 + 1)
C     ENDDO
C     ENDIF
C     ENDLOCAL
C     ENDPDO
C     total_time = (total_time + ((timer_stop() * 8) / 2e7))
C     GOTO 101 
C     ELSE
C     IF (cont.EQ.2.) THEN
C     WRITE (unit=6,fmt=*) 'Processors = ',32,'; Elapsed time =',
C    $total_time
C     DO d = 1, 50, 1
C     DO j = 1, 50, 1
C     DO k = 1, 50, 1
C     sumsq(k,j) = workspace(d,k,j)
C     ENDDO
C     ENDDO
C     WRITE (unit=11) ((sumsq(k,j),k=1,50),j=1,50)
C     ENDDO
C     PDO (d = 1:50)
C     POUT workspace(d,:,:)
C     LOCAL
C     IF (fxcellid.LT.32) THEN
C     fxloopstart0 = (MAX(((fxcellid * 2) + 1),1))
C     fxlmidx0 = ifxlm(fxworkspacedesc,0,fxloopstart0)
C     DO d = fxloopstart0, MIN(((fxcellid * 2) + 2),50), 1
C     DO i = 1, 50, 1
C     DO j = 1, 50, 1
C     workspace(d,j,i) = 0.
C     ENDDO
C     ENDDO
C     fxlmidx0 = (fxlmidx0 + 1)
C     ENDDO
C     ENDIF
C     ENDLOCAL
C     ENDPDO
C     total_time = 0.
C     GOTO 101 
C     ENDIF
C     ENDIF
C     CLOSE (unit=10)
C     CLOSE (unit=11)
C     STOP 
C     END
      PROGRAM main
      INTEGER fxlmidx0
      REAL rfxtemp0
      CHARACTER*100 filename
      REAL cont
      EXTERNAL variance
      EXTERNAL warp
      INTEGER fxcellid
      REAL minx
      REAL miny
      REAL minz
      REAL maxx
      REAL maxy
      REAL maxz
      INTEGER fxloopstart0
      INTEGER d
      INTEGER i
      INTEGER j
      INTEGER k
      EXTERNAL fx_procstart
      INTEGER blank
      REAL depth
      INTRINSIC float
      REAL float
      INTRINSIC index
      INTEGER index
      INTEGER nimgs
      EXTERNAL ifxlm
      INTEGER ifxlm
      EXTERNAL fx_deinit
      REAL dummy
      INTEGER fxdistarg
      EXTERNAL ifxownsol
      INTEGER ifxownsol
      EXTERNAL getarg
      EXTERNAL msg_reduce
      EXTERNAL fxinitdesc
      REAL total_time
      EXTERNAL timer_stop
      INTEGER timer_stop
      LOGICAL created
      EXTERNAL ifxceil
      INTEGER ifxceil
      EXTERNAL fx_init
      EXTERNAL fx_set_module_geometry
      EXTERNAL fxbroadcast
      EXTERNAL fxreceive_bcast
      EXTERNAL ifxowns
      INTEGER ifxowns
      EXTERNAL timer_start
      INTEGER fxworkspacedesc(71,3)
      REAL sumsq(50,50)
      REAL workspace(2,50,50)
      REAL tarray(2)
      REAL proj_xform(4,3,4)
      REAL img(4,528,486)
      REAL sum(50,50)
      REAL tempimg(50,50)
      INTEGER fxdistarray(6,20)
      COMMON proj_xform,nimgs
      COMMON /fxcommon/ fxcellid,fxdistarg
      CALL fx_init(fxcellid,1,1,0,0)
      CALL fx_set_module_geometry(0,0,0,8,8)
      fxdistarg = 0
      fxdistarray(1,1) = 1
      fxdistarray(2,1) = 0
      fxdistarray(3,1) = 2
      fxdistarray(4,1) = 32
      fxdistarray(5,1) = 0
      fxdistarray(6,1) = 0
      fxdistarray(3,2) = 50
      fxdistarray(3,3) = 50
      CALL fxinitdesc(1,3,fxdistarray,fxworkspacedesc)
C     ALIGN workspace(I,J,K) WITH t(I)
C     TEMPLATE t(50)
C     DISTRIBUTE t(BLOCK_CYCLIC(2,32))
      CALL getarg(1,filename)
      blank = index(filename,' ')
      IF (fxcellid.EQ.0) THEN
        WRITE (unit=6,fmt=*) 'Opening ',filename,' for input'
      ENDIF
      IF (fxcellid.EQ.0) THEN
        OPEN (unit=10,file=filename,form='unformatted',status='old')
      ENDIF
      CALL getarg(2,filename)
      blank = index(filename,' ')
      IF (fxcellid.EQ.0) THEN
        WRITE (unit=6,fmt=*) 'Opening ',filename,' for output'
      ENDIF
      IF (fxcellid.EQ.0) THEN
        INQUIRE (file=filename,exist=created)
      ENDIF
      IF (created) THEN
        IF (fxcellid.EQ.0) THEN
          OPEN (unit=11,file=filename,form='unformatted',status='old'
     $)
        ENDIF
      ELSE
        IF (fxcellid.EQ.0) THEN
          OPEN (unit=11,file=filename,form='unformatted',status='new'
     $)
        ENDIF
      ENDIF
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) dummy
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(dummy,1)
      ELSE
        CALL fxreceive_bcast(dummy,1)
      ENDIF
      nimgs = dummy
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) minx
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(minx,1)
      ELSE
        CALL fxreceive_bcast(minx,1)
      ENDIF
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) miny
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(miny,1)
      ELSE
        CALL fxreceive_bcast(miny,1)
      ENDIF
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) minz
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(minz,1)
      ELSE
        CALL fxreceive_bcast(minz,1)
      ENDIF
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) maxx
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(maxx,1)
      ELSE
        CALL fxreceive_bcast(maxx,1)
      ENDIF
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) maxy
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(maxy,1)
      ELSE
        CALL fxreceive_bcast(maxy,1)
      ENDIF
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) maxz
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(maxz,1)
      ELSE
        CALL fxreceive_bcast(maxz,1)
      ENDIF
C     PDO (d = 1:50)
C     POUT workspace(d,:,:)
C     LOCAL
      IF (fxcellid.LT.32) THEN
        fxloopstart0 = (MAX(((fxcellid * 2) + 1),1))
        fxlmidx0 = ifxlm(fxworkspacedesc,0,fxloopstart0)
        DO d = fxloopstart0, MIN(((fxcellid * 2) + 2),50), 1
          DO i = 1, 50, 1
            DO j = 1, 50, 1
              workspace(fxlmidx0,j,i) = 0.
            ENDDO
          ENDDO
          fxlmidx0 = (fxlmidx0 + 1)
        ENDDO
      ENDIF
C     ENDLOCAL
C     ENDPDO
      total_time = 0.
101   CONTINUE
      IF (fxcellid.EQ.0) THEN
        READ (unit=10) cont
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CALL fxbroadcast(cont,1)
      ELSE
        CALL fxreceive_bcast(cont,1)
      ENDIF
      IF (cont.EQ.1.) THEN
        DO i = 1, nimgs, 1
          DO j = 1, 3, 1
            IF (fxcellid.EQ.0) THEN
              READ (unit=10) (proj_xform(i,j,k),k=1,4)
            ENDIF
            IF (fxcellid.EQ.0) THEN
              CALL fxbroadcast(proj_xform,48)
            ELSE
              CALL fxreceive_bcast(proj_xform,48)
            ENDIF
          ENDDO
        ENDDO
        IF (fxcellid.EQ.0) THEN
          READ (unit=10) (((img(i,k,j),k=1,528),j=1,486),i=1,nimgs)
        ENDIF
        IF (fxcellid.EQ.0) THEN
          CALL fxbroadcast(img,1026432)
        ELSE
          CALL fxreceive_bcast(img,1026432)
        ENDIF
        CALL timer_start()
C     PDO (d = 0:49)
C     POUT workspace(d,:,:)
C     LOCAL
        IF (fxcellid.LT.32) THEN
          fxloopstart0 = (MAX(((fxcellid * 2) + 1),0))
          fxlmidx0 = ifxlm(fxworkspacedesc,0,fxloopstart0)
          DO d = fxloopstart0, MIN(((fxcellid * 2) + 2),49), 1
            depth = (minz + (((d + (-1)) * (maxz - minz)) / 50))
            DO j = 1, 50, 1
              DO k = 1, 50, 1
                sum(k,j) = 0.
                sumsq(k,j) = 0.
              ENDDO
            ENDDO
            DO i = 1, nimgs, 1
              CALL warp(img,tempimg,i,minx,miny,maxx,maxy,depth)
              DO j = 1, 50, 1
                DO k = 1, 50, 1
                  sum(k,j) = (sum(k,j) + tempimg(k,j))
                  sumsq(k,j) = (sumsq(k,j) + (tempimg(k,j) ** 2))
                ENDDO
              ENDDO
            ENDDO
            CALL variance(sum,sumsq)
            DO i = 1, 50, 1
              DO j = 1, 50, 1
                workspace(fxlmidx0,j,i) = (workspace(fxlmidx0,j,i) + 
     $(1 / sumsq(j,i)))
              ENDDO
            ENDDO
            fxlmidx0 = (fxlmidx0 + 1)
          ENDDO
        ENDIF
C     ENDLOCAL
C     ENDPDO
        total_time = (total_time + ((timer_stop() * 8) / 2e7))
        GOTO 101 
      ELSE
        IF (cont.EQ.2.) THEN
          IF (fxcellid.EQ.0) THEN
            WRITE (unit=6,fmt=*) 'Processors = ',32,
     $'; Elapsed time =',total_time
          ENDIF
          DO d = 1, 50, 1
            DO j = 1, 50, 1
              DO k = 1, 50, 1
                IF ((ifxowns(fxworkspacedesc,0,d,fxcellid).NE.0).AND.
     $(fxworkspacedesc((fxcellid + 8),1).NE.(-1000000))) THEN
                  rfxtemp0 = workspace(ifxlm(fxworkspacedesc,0,d),k,j
     $)
                  CALL fxbroadcast(rfxtemp0,1)
                ELSE
                  CALL fxreceive_bcast(rfxtemp0,1)
                ENDIF
                sumsq(k,j) = rfxtemp0
              ENDDO
            ENDDO
            IF (fxcellid.EQ.0) THEN
              WRITE (unit=11) ((sumsq(k,j),k=1,50),j=1,50)
            ENDIF
          ENDDO
C     PDO (d = 1:50)
C     POUT workspace(d,:,:)
C     LOCAL
          IF (fxcellid.LT.32) THEN
            fxloopstart0 = (MAX(((fxcellid * 2) + 1),1))
            fxlmidx0 = ifxlm(fxworkspacedesc,0,fxloopstart0)
            DO d = fxloopstart0, MIN(((fxcellid * 2) + 2),50), 1
              DO i = 1, 50, 1
                DO j = 1, 50, 1
                  workspace(fxlmidx0,j,i) = 0.
                ENDDO
              ENDDO
              fxlmidx0 = (fxlmidx0 + 1)
            ENDDO
          ENDIF
C     ENDLOCAL
C     ENDPDO
          total_time = 0.
          GOTO 101 
        ENDIF
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CLOSE (unit=10)
      ENDIF
      IF (fxcellid.EQ.0) THEN
        CLOSE (unit=11)
      ENDIF
      CALL fx_deinit()
      STOP 
      END
C     DO i = 1, 50, 1
C     DO j = 1, 50, 1
C     sumsq(j,i) = (((nimgs * sumsq(j,i)) - (sum(j,i) ** 2)) / (nimgs
C    $ * (nimgs + (-1))))
C     ENDDO
C     ENDDO
C     END
      SUBROUTINE variance(sum,sumsq)
      INTEGER fxcellid
      INTEGER i
      INTEGER j
      EXTERNAL fx_procstart
      INTEGER nimgs
      INTEGER fxdistarg
      EXTERNAL fxdblcheck
      EXTERNAL msg_reduce
      REAL sumsq(50,50)
      REAL proj_xform(4,3,4)
      REAL sum(50,50)
      COMMON proj_xform,nimgs
      COMMON /fxcommon/ fxcellid,fxdistarg
      CALL fxdblcheck(fxdistarg,0,'variance')
      DO i = 1, 50, 1
        DO j = 1, 50, 1
          sumsq(j,i) = (((nimgs * sumsq(j,i)) - (sum(j,i) ** 2)) / (
     $nimgs * (nimgs + (-1))))
        ENDDO
      ENDDO
      END
C     DO i = 1, 3, 1
C     DO j = 1, 4, 1
C     trans(i,j) = proj_xform(imgnum,i,j)
C     ENDDO
C     ENDDO
C     randint = irand(0)
C     randint = (randint - ((randint / 256) * 256))
C     DO i = 1, 50, 1
C     u = ((((trans(1,1) * minx) + (trans(1,2) * (miny + (((i + (-1))
C    $ * (maxy - miny)) / 50)))) + (trans(1,3) * z)) + trans(1,4))
C     v = ((((trans(2,1) * minx) + (trans(2,2) * (miny + (((i + (-1))
C    $ * (maxy - miny)) / 50)))) + (trans(2,3) * z)) + trans(2,4))
C     w = ((((trans(3,1) * minx) + (trans(3,2) * (miny + (((i + (-1))
C    $ * (maxy - miny)) / 50)))) + (trans(3,3) * z)) + trans(3,4))
C     DO j = 1, 50, 1
C     pu = (((u / w) + 264.) + 1)
C     pv = (((v / w) + 243.) + 1)
C     iu = pu
C     iv = pv
C     IF ((((iu.GE.1).AND.(iu.LT.528)).AND.(iv.GE.1)).AND.(iv.LT.486)
C    $) THEN
C     du = (pu - iu)
C     dv = (pv - iv)
C     a = inimg(imgnum,iu,iv)
C     b = inimg(imgnum,iu,(iv + 1))
C     c = inimg(imgnum,(iu + 1),iv)
C     d = inimg(imgnum,(iu + 1),(iv + 1))
C     mid1val = (a + ((b - a) * dv))
C     mid2val = (c + ((d - c) * dv))
C     outimg(j,i) = (mid1val + ((mid2val - mid1val) * du))
C     ELSE
C     outimg(j,i) = (randint)
C     ENDIF
C     u = (u + ((trans(1,1) * (maxx - minx)) / 50))
C     v = (v + ((trans(2,1) * (maxx - minx)) / 50))
C     w = (w + ((trans(3,1) * (maxx - minx)) / 50))
C     ENDDO
C     ENDDO
C     END
      SUBROUTINE warp(inimg,outimg,imgnum,minx,miny,maxx,maxy,z)
      INTEGER fxcellid
      REAL minx
      REAL miny
      REAL maxx
      REAL maxy
      REAL a
      REAL b
      REAL c
      REAL d
      INTEGER i
      INTEGER j
      EXTERNAL fx_procstart
      REAL u
      REAL v
      REAL w
      REAL z
      EXTERNAL irand
      INTEGER irand
      INTEGER nimgs
      INTEGER fxdistarg
      REAL du
      REAL dv
      INTEGER iu
      INTEGER iv
      REAL pu
      REAL pv
      EXTERNAL fxdblcheck
      INTEGER imgnum
      EXTERNAL msg_reduce
      REAL mid1val
      REAL mid2val
      INTEGER randint
      REAL inimg(4,528,486)
      REAL trans(3,4)
      REAL outimg(50,50)
      REAL proj_xform(4,3,4)
      COMMON proj_xform,nimgs
      COMMON /fxcommon/ fxcellid,fxdistarg
      CALL fxdblcheck(fxdistarg,0,'warp')
      DO i = 1, 3, 1
        DO j = 1, 4, 1
          trans(i,j) = proj_xform(imgnum,i,j)
        ENDDO
      ENDDO
      randint = irand(0)
      randint = (randint - ((randint / 256) * 256))
      DO i = 1, 50, 1
        u = ((((trans(1,1) * minx) + (trans(1,2) * (miny + (((i + 
     $(-1)) * (maxy - miny)) / 50)))) + (trans(1,3) * z)) + trans(1,4
     $))
        v = ((((trans(2,1) * minx) + (trans(2,2) * (miny + (((i + 
     $(-1)) * (maxy - miny)) / 50)))) + (trans(2,3) * z)) + trans(2,4
     $))
        w = ((((trans(3,1) * minx) + (trans(3,2) * (miny + (((i + 
     $(-1)) * (maxy - miny)) / 50)))) + (trans(3,3) * z)) + trans(3,4
     $))
        DO j = 1, 50, 1
          pu = (((u / w) + 264.) + 1)
          pv = (((v / w) + 243.) + 1)
          iu = pu
          iv = pv
          IF ((((iu.GE.1).AND.(iu.LT.528)).AND.(iv.GE.1)).AND.(iv.LT.
     $486)) THEN
            du = (pu - iu)
            dv = (pv - iv)
            a = inimg(imgnum,iu,iv)
            b = inimg(imgnum,iu,(iv + 1))
            c = inimg(imgnum,(iu + 1),iv)
            d = inimg(imgnum,(iu + 1),(iv + 1))
            mid1val = (a + ((b - a) * dv))
            mid2val = (c + ((d - c) * dv))
            outimg(j,i) = (mid1val + ((mid2val - mid1val) * du))
          ELSE
            outimg(j,i) = (randint)
          ENDIF
          u = (u + ((trans(1,1) * (maxx - minx)) / 50))
          v = (v + ((trans(2,1) * (maxx - minx)) / 50))
          w = (w + ((trans(3,1) * (maxx - minx)) / 50))
        ENDDO
      ENDDO
      END
