C ++******************************************************************** C * C * 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 * C * C PURPOSE: * C * C PARAMETERS: * C * C 0 2 3 4 5 6 7 * C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** SUBROUTINE REPR2_S(BCKE,BCN,NSAM,LCYL,NSLICE, & NROWL,NROWH,NANG,IPCUBE,NN,PROJ,PRN,IRI,LTB,LTBN,ABA,INPIC) INCLUDE 'CMBLOCK.INC' DIMENSION BCKE(Nsam*nslice*3),BCN(Nsam*nslice) dimension PROJ(NSAM,NANG),PRN(NSAM,NANG) dimension ipcube(5,nn) double precision aba C logical*1 active_min,active_max,endit,dodo,dod c data IOFF/6/ CALL RDPRM2(ALA,aim,NOT_USED,'Lambda, Correction limit') CALL RDPRMI(maxit,inter,NOT_USED, & 'Iteration limit, Number of internal steps') CALL RDPRMI(MODE,idum,NOT_USED,'Mode') CALL RDPRM2(tmin,tmax,NOT_USED,'Minimum, maximum') tmin=tmin-aba tmax=tmax-aba write(nout,2059) tmin,tmax 2059 format(' Minimum and maximum after average subtraction',/, & 2(5x,1pe10.3)) CALL RDPRM(t,NOT_USED,'Smoothing constant') t=amin1(amax1(t,0.0),1.0) c if(.not. & (mode.eq.1.or.mode.eq.3.or.mode.eq.6.or.mode.eq.8)) then if(maxit.gt.1) then inter=inter*maxit maxit=1 dodo=.false. write(nout,6066) inter 6066 format(/ & ' Since no smoothing was requested iteration limit was set to 1' & ,/,' and number of internal sets was set to =',i5) else dodo=.true. endif endif c nmat=nsam*nslice active_min=.false. active_max=.false. sqold=1.e23 endit=.false. ITER = 0 DO WHILE(.not.endit) ITER = ITER + 1 WRITE(NOUT,971) ITER 971 FORMAT(/' ITERATION ',I5) c sq=0.0 do ls=nrowl,nrowh if(ITER.gt.1) then do irc=1,nslice irb=1+(irc-1)*nsam call redlin(inpic,bcke(irb),nsam,(irc-1)*lcyl+ls-nrowl+1) enddo else do i=1,nmat bcke(i)=0.0 enddo endif c c get the projection data c do j=1,nang call redlin(IOFF+j,proj(1,j),nsam,ls) enddo c c remove average c do j=1,nang do i=1,nsam proj(i,j)=proj(i,j)-aba enddo enddo sqint=1.0e23 do itr=1,inter if(.not.endit) then if(ITER.eq.1.and.itr.eq.1) then c$omp parallel do private(i,j) do j=1,nang do i=1,nsam prn(i,j)=proj(i,j) enddo enddo c goto 981 else c bcke -> prn call prjs2 if(dodo) then dod=itr.eq.inter else dod=.true. endif c if(dod) then call fmin_s(prn,ltbn,gmin) call fmax_s(prn,ltbn,gmax) write(nout,2062) ls,gmin,gmax 2062 format(' Min Max in slice #',i3,' =',2(1x,1pe10.3) ) endif if((mode.eq.2.or.mode.eq.3.or.mode.eq.7.or.mode.eq.8) & .and. .not.active_min .and. dod) then if(gmin.lt.tmin) then call bmin_s(bcke,nmat,ipcube,nn,bmin) write(nout,2051) bmin 2051 format(' Min constraint activated, value in 3D =',1pe10.3) active_min=.true. endif endif c if((mode.eq.5.or.mode.eq.6.or.mode.eq.7.or.mode.eq.8) & .and. .not.active_max .and. dod) then if(gmax.gt.tmax) then call bmax_s(bcke,nmat,ipcube,nn,bmax) write(nout,2052) bmax 2052 format(' Max constraint activated, value in 3D =',1pe10.3) active_max=.true. endif endif c c$omp parallel do private(i,j) do j=1,nang do i=1,nsam prn(i,j)=proj(i,j)-prn(i,j) enddo enddo c 981 endif sqt=0.0 c$omp parallel do private(i,j),reduction(+:sqt) do j=1,nang do i=1,nsam sqt=sqt+prn(i,j)*prn(i,j) enddo enddo c if(sqt.le.sqint) then sqint=sqt if(dod) sq=sq+sqt c c prn -> bcn call bprj2 c call docorr3_s(bcke,bcn,nmat,ipcube,nn,ala,sbq) write(nout,*) ' Squared correction of the structure ',ls,itr,sbq if(mode.ne.0) then if(active_min) call domin3_s(bcke,nmat,ipcube,nn,bmin) if(active_max) call domax3_s(bcke,nmat,ipcube,nn,bmax) endif else c break it up, the error increased. endit=.true. endif c end of do-loop over itr endif enddo c c5078 continue do irc=1,nslice irb=1+(irc-1)*nsam call wrtlin(inpic,bcke(irb),nsam,(irc-1)*lcyl+ls-nrowl+1) enddo c enddo c write(nout,*) ' Squared discrepancies between projections ',sq if(mode.eq.1.or.mode.eq.3.or.mode.eq.6.or.mode.eq.8) then do ls=nrowl+1,nrowh-1 if(ls.eq.nrowl+1) then do irc=1,nslice irb=1+(irc-1)*nsam call redlin(inpic,bcke(irb),nsam,(irc-1)*lcyl+ls-1-nrowl+1) enddo do irc=1,nslice irb=1+(irc-1)*nsam+nsam*nslice call redlin(inpic,bcke(irb),nsam,(irc-1)*lcyl+ls-nrowl+1) enddo endif do irc=1,nslice irb=1+(irc-1)*nsam+2*nsam*nslice call redlin(inpic,bcke(irb),nsam,(irc-1)*lcyl+ls+1-nrowl+1) enddo call smt3(t,bcke,bcn,nsam,nslice,ipcube,nn) do irc=1,nslice irb=1+(irc-1)*nsam call wrtlin(inpic,bcn(irb),nsam,(irc-1)*lcyl+ls-nrowl+1) enddo call cppb(bcke(nsam*nslice+1),bcke,2*nsam*nslice) enddo endif if(.not.endit) then if(sq.gt.aim .and. iter.lt.maxit) then if(iter.eq.1) then c continue iterations elseif(sq.lt.sqold) then sqold=sq c continue iterations else endit=.true. c stop iterations endif endif endif C ENDDO END