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 MAIN SUBROUTINE FOR MARKER PROGRAM. C C WARNING ! C PARAMETERS LS AND LV ARE HARD-WIRED IN ALL THE SUBROUTINES. C C ********************************************************** C ATTEMPTS TO ALIGN A DATA SET USING FIDUCIARY MARKERS. C C THE PROCEDURE IS TO FIRST ROUGHLY TRANSLATIONALLY C ALIGN THE IMAGES. USING THIS SET, RECONSTRUCT THE C COORDS OF THE POINTS IN THE 3D OBJECT. THEN PROJECT C THE POINTS BACK DOWN ONTO A PLANE AND COMPARE THEM C WITH THE POINTS FROM PREVIOUS TRIAL. CORRECT THE C POINTS AND ONCE AGAIN RECONSTRUCT THE 3D, ETC. C C THE CORRECTIONS ACCUMULATE IN A FEW REGISTERS. C ALL COMPARISONS EXCEPT FOR THE ROUGH ALIGNMENT ARE C BETWEEN THE ORIGINAL PROJECTIONS AND THEIR PSUEDO- C PROJECTIONS. C ********************************************************** SUBROUTINE MRK(MAXDIM) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' PARAMETER (LV=300) PARAMETER (LS=256) PARAMETER (MAXLOOP=19999) REAL MAXPER PARAMETER (MAXPER=1.0E-7) PARAMETER (NLIST=2) DIMENSION DLIST(NLIST) LOGICAL*1 FULL CHARACTER*1 ASKFULL,NULL CHARACTER(LEN=MAXNAM):: DOCNAM LOGICAL*1 PTACTIVE(LS,LV),PARAMQ(4) COMMON /GENERAL/PTACTIVE,NUMPTS(LV),NTVW,NTPT,CIR(2),PARAMQ DIMENSION XYPTS(2,LS,LV),PRJ(2,LS,LV) DIMENSION ROT(3,3),ANGLE(3,LV),TPSI(LV) DIMENSION P3D(3,LS),SCALE(LV) DIMENSION NUMBER(LV),SCALEI(LV) DIMENSION CVPT(3) DIMENSION SHFT(2),SHIFT(2,LV),TSHIFT(2,LV) DIMENSION ERVW(LV),ERPT(LS) LOGICAL*1 FIRST INTEGER NTIMES DOUBLE PRECISION TCH,ERTOT DATA PI/3.141592654/ NDOUTT = 70 NTIMES = 0 NULL = CHAR(0) C GET COORDS OF MARKERS IN EACH PROJECTION CALL MRGETINFO(XYPTS,IREF,NUMBER,ANGLE,SCALE,TSHIFT,FIRST,NTTVL) DO IVIEW=1,NTVW TPSI(IVIEW)=ANGLE(1,IVIEW) SHIFT(1,IVIEW)=TSHIFT(1,IVIEW) SHIFT(2,IVIEW)=TSHIFT(2,IVIEW) CALL MRALIGN(XYPTS(1,1,IVIEW),XYPTS(1,1,IVIEW), & ANGLE(1,IVIEW),SHIFT(1,IVIEW),SCALE(IVIEW), & PTACTIVE(1,IVIEW),NTPT) ENDDO C GET 3D VOLUME CALL MR2TO3D(P3D, XYPTS, ANGLE) CALL MRPROJ(P3D, PRJ, ANGLE) CALL MRCALERR(XYPTS,PRJ,ERTOT,ERVW,ERPT) WRITE(NOUT,*)' INITIAL INFO ' FULL=.FALSE. CALL MRPUTINFO(XYPTS,PRJ,ANGLE,P3D, & NUMBER,TSHIFT,SCALE,ERVW,ERPT,FULL) TCH=ERTOT WRITE(NOUT,701) NTIMES,ERTOT 701 FORMAT(' PASS #',I5,' TOTCHANGE= ',10('*'),' ERTOT= ',G11.4) 700 FORMAT(' PASS #',I5,' TOTCHANGE= ',G10.3,' ERTOT= ',G11.4) C CREATE PSUEDO-PROJECTIONS AND COMPARE TO LAST RUN 200 NTIMES=NTIMES+1 C DO 210 ONLY IF _IN PLANE_ ROTATION OR _SCALE_ CORRECTIONS REQUESTED. IF (PARAMQ(1).OR.PARAMQ(3)) THEN SHFT(1)=0.0 SHFT(2)=0.0 DO 210 JVIEW=0,NTVW C THIS STRANGE TRICK IS TO MAKE SURE THAT FIRST PASS IS OVER REFERENCE C IMAGE, IN SUBSEQUENT ITERATIONS REFERENCE IMAGE IS SKIPPED. IVIEW=JVIEW IF(JVIEW.EQ.IREF) GOTO 210 IF(JVIEW.EQ.0) IVIEW=IREF C ASSUMPTION IS THAT TILT OF PROJ AND TILT OF ORIG ARE EQUAL C THIS IS THE ONLY REASON WE CAN REALLY COMPARE THE TWO C CALL MRQUATER(RPT,VPT,IVIEW,ROT) C CALL MRROTATE(ROT,PHI2,THETA,PSI) C FIND PSI (IN PLANE) ANGLE IF (PARAMQ(1)) THEN CALL MRANG2(PRJ(1,1,IVIEW),XYPTS(1,1,IVIEW),IVIEW,PSI) ELSE PSI=0.0 ENDIF C FIND SCALE IF (PARAMQ(3)) THEN CALL MRSCALE(PRJ(1,1,IVIEW),XYPTS(1,1,IVIEW),IVIEW,SCALEI) SCAL=SCALEI(IVIEW) IF (IVIEW.EQ.IREF) SCAD=1/SCALEI(IVIEW) SCALEI(IVIEW)=SCALEI(IVIEW)*SCAD ELSE SCAL=1.0 SCALEI(IVIEW)=1.0 ENDIF C UPDATE CORRECTION INFO IN REGISTERS, NO SHIFTS AT THIS MOMENT C=COS(PSI) S=SIN(PSI) QT=TSHIFT(1,IVIEW) TSHIFT(1,IVIEW)=SCALEI(IVIEW)*(QT*C+TSHIFT(2,IVIEW)*S) TSHIFT(2,IVIEW)=SCALEI(IVIEW)*(-QT*S+TSHIFT(2,IVIEW)*C) TPSI(IVIEW)=PSI+TPSI(IVIEW) SCALE(IVIEW)=SCALE(IVIEW)*SCALEI(IVIEW) ANGLE(1,IVIEW)=PSI c apply corrections to points with SHFT set to zero ... CALL MRALIGN(XYPTS(1,1,IVIEW),XYPTS(1,1,IVIEW), & PSI,SHFT,SCALEI(IVIEW),PTACTIVE(1,IVIEW),NTPT) 210 CONTINUE C NOW GOTO 3D CALL MR2TO3D(P3D, XYPTS, ANGLE) ENDIF C FIND TRANSLATION (SHIFT) IF (PARAMQ(4)) THEN C FIND CENTER OF GRAVITY OF POINTS IN 3D C AND APPLY ANY SHIFTS FOUND TO THESE POINTS CALL MRCG3D(P3D) C GOTO 2D CALL MRPROJ(P3D, PRJ, ANGLE) C FIND THE SHIFT BETWEEN POINTS AND SHIFT XYPTS CALL MRSHIFT(PRJ,XYPTS,SHIFT) C UPDATE PARAMETERS DO IVIEW=1,NTVW TSHIFT(1,IVIEW)=TSHIFT(1,IVIEW)+SHIFT(1,IVIEW) TSHIFT(2,IVIEW)=TSHIFT(2,IVIEW)+SHIFT(2,IVIEW) ENDDO C GOTO 3D AGAIN CALL MR2TO3D(P3D,XYPTS,ANGLE) ENDIF C FIND THETA (TILT) ANGLE IF (PARAMQ(2)) THEN C GOTO 2D CALL MRPROJ(P3D, PRJ, ANGLE) DO JVIEW=0,NTVW C THIS STRANGE TRICK IS TO MAKE SURE THAT FIRST PASS IS OVER REFERENCE C IMAGE, IN SUBSEQUENT ITERATIONS REFERENCE IMAGE IS SKIPPED. C WON'T MODIFY XYPTS, ONLY THETA IVIEW=JVIEW IF ( JVIEW.NE.IREF) THEN IF (JVIEW.EQ.0) IVIEW=IREF THETA = ANGLE(2,IVIEW) C PRINT *,IVIEW,THETA C SET SCALE=1.0 AND PSI=0.0 CALL MRTHETA(PRJ(1,1,IVIEW),XYPTS(1,1,IVIEW),IVIEW,P3D,THETA) IF (IVIEW.EQ.IREF) THAD = THETA THETA=THETA-THAD CC CORRECTION OF THETA IS FOUND, SO CALCULATE NEW THETA C ANGLE(2,IVIEW)=ANGLE(2,IVIEW)+THETA ANGLE(2,IVIEW)=THETA ENDIF ENDDO C GOTO 3D AGAIN CALL MR2TO3D(P3D, XYPTS, ANGLE) ENDIF C GOTO 2D AND GET THE ERROR ... CALL MRPROJ(P3D,PRJ,ANGLE) CALL MRCALERR(XYPTS,PRJ,ERTOT,ERVW,ERPT) TTCH=SNGL((TCH-ERTOT)/TCH) C IF (MOD(NTIMES,10).EQ.0) THEN C WRITE(NOUT,700)NTIMES,BOXERR,TTCH,ERTOT C WRITE(NDAT,700)NTIMES,BOXERR,TTCH,ERTOT C ENDIF TCH=ERTOT IF (MOD(NTIMES,150).EQ.0) THEN WRITE(NOUT,700)NTIMES,TTCH,ERTOT FULL=.FALSE. CALL MRPUTINFO(XYPTS,PRJ,ANGLE,P3D, & NUMBER,TSHIFT,SCALE,ERVW,ERPT,FULL) ENDIF C DO THE NEXT ITERATION ? IF (ABS(TTCH) .GT. MAXPER & .AND. NTIMES.LT.MAXLOOP) GOTO 200 IF (NTIMES.GE.MAXLOOP) & WRITE(NOUT,*)' EXCEEDED LOOP LIMIT, ENDING CYCLE' C JUST DO SOME CLEANING UP BEFORE PUTTING OUT VALUES TILR=ANGLE(2,IREF) c XSH=TSHIFT(1,IREF) c YSH=TSHIFT(2,IREF) DO IVIEW=1,NTVW ANGLE(1,IVIEW)=TPSI(IVIEW) ANGLE(2,IVIEW)=ANGLE(2,IVIEW)-TILR C THIS IS DONE TO PUT TILT AXIS AT (0,0) OF REFERENCE VIEW C TSHIFT(1,IVIEW)=TSHIFT(1,IVIEW)-(XSH*COS(ANGLE(2,IVIEW))) C TSHIFT(2,IVIEW)=TSHIFT(2,IVIEW)-YSH ENDDO WRITE(NOUT,700)NTIMES,TTCH,ERTOT CALL RDPRMC(ASKFULL,NUMC,.TRUE.,'FULL OUTPUT (Y/N)',NULL,IRT) FULL = (ASKFULL.NE.'N') CALL MRPUTINFO(XYPTS,PRJ,ANGLE,P3D, & NUMBER,TSHIFT,SCALE,ERVW,ERPT,FULL) 987 CALL OPENDOC(DOCNAM,.TRUE.,NLET,NDOUTT,NDOUT,.TRUE., & 'ERROR PER VIEW OUTPUT DOC', & .FALSE.,.FALSE.,.TRUE.,NEWFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 51 IF (NDOUT .LE. 0) THEN CALL ERRT(102,'INCORE DOC FILE NOT ALLOWED HERE.',NE) GOTO 987 ENDIF DO IVIEW=1,NTVW C DLIST(1) = NUMBER(IVIEW) DLIST(2) = ERVW(IVIEW) CALL LUNDOCWRTDAT(NDOUT,NUMBER(IVIEW),DLIST(2),1,IRTFLG) ENDDO C DLIST(1) = -1 DLIST(2) = ERTOT CALL LUNDOCWRTDAT(NDOUT,-1,DLIST(2),1,IRTFLG) 51 CLOSE(NDOUT) 986 CALL OPENDOC(DOCNAM,.TRUE.,NLET,NDOUTT,NDOUT,.TRUE., & 'ERROR PER MARKER OUTPUT DOC', & .FALSE.,.FALSE.,.TRUE.,NEWFILE,IRTFLG) IF (IRTFLG .NE. 0) GOTO 52 IF (NDOUT .LE. 0) THEN CALL ERRT(102,'INCORE DOC FILE NOT ALLOWED HERE.',NE) GOTO 986 ENDIF DO IPT=1,NTPT C DLIST(1) = IPT DLIST(2) = ERPT(IPT) CALL LUNDOCWRTDAT(NDOUT,IPT,DLIST(2),1,IRTFLG) ENDDO C DLIST(1) = -1 DLIST(2) = ERTOT CALL LUNDOCWRTDAT(NDOUT,-1,DLIST(2),1,IRTFLG) 52 CLOSE(NDOUTT) CALL MRDOC(SCALE,TSHIFT,ANGLE,NUMBER,P3D) END