C ++******************************************************************** C * C EROSION CREATED FEB 01 ARDEAN LEITH * C ADDED 'ER DOC' MAR 01 ARDEAN LEITH * C ADDED NPIXER MAR 01 ARDEAN LEITH * C INCORE LUNDOC JUL 03 ARDEAN LEITH * C ********************************************************************** C=* * C=* This file is part of: SPIDER - Modular Image Processing System. * C=* SPIDER System Authors: Joachim Frank & ArDean Leith * C=* Copyright 1985-2010 Health Research Inc., * C=* Riverview Center, 150 Broadway, Suite 560, Menands, NY 12204. * C=* Email: spider@wadsworth.org * C=* * C=* SPIDER 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=* SPIDER 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=* You should have received a copy of the GNU General Public License * C=* along with this program. If not, see * C=* * C ********************************************************************** C C EROSION(LUN1,LUN2,NSAM,NROW,NSLICE) C C PARAMETERS: C C PURPOSE: ERODE (SHRINK) AN OBJECT IN AN IMAGE OR VOLUME C * C ********************************************************************** SUBROUTINE EROSION(LUN1,LUN2,NSAM,NROW,NSLICE) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER(LEN=MAXNAM) :: DOCNAM REAL, ALLOCATABLE, DIMENSION(:,: ) :: VIN2 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: VIN3 CHARACTER(LEN=1) :: MODE LOGICAL :: NEWFILE LUNDOC = 0 IF (FCHAR(4:4) .EQ. 'D') LUNDOC = 77 CALL RDPRMC(MODE,NA,.TRUE.,'BOX OR CROSS NEIGHBORHOOD (B/C)', & CHAR(0),IRTFLG) IF (IRTFLG .NE. 0) RETURN LENGTH = 3 10 CALL RDPRI1S(LENGTH,NOT_USED,'LENGTH OF NEIGHBORHOOD',IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (LENGTH .LE. 1) THEN CALL ERRT(102,'LENGTH MUST BE GREATER THAN',2) GOTO 10 ELSEIF (MOD(LENGTH,2) .EQ. 0) THEN LENGTH = LENGTH + 1 WRITE(NOUT,90) LENGTH 90 FORMAT(' EFFECTIVE LENGTH OF NEIGHBORHOOD:',I5) ENDIF LH = LENGTH / 2 IF (NSLICE .LE. 1) THEN IF (MODE .EQ. 'B') THEN LAT = LENGTH * LENGTH - 1 ELSE LAT = 2 * LENGTH - 2 ENDIF ELSE IF (MODE .EQ. 'B') THEN LAT = LENGTH * LENGTH * LENGTH - 1 ELSE LAT = 3 * LENGTH - 3 ENDIF ENDIF WRITE(NOUT,91) LAT 91 FORMAT(' NUMBER OF NEIGHBORS:',I5) CALL RDPRI1S(LAT,NOT_USED, & 'ERODE IF NUMBER OF OCCUPIED NEIGHBORS IS < THAN',IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (LUNDOC .GT. 0) THEN CALL OPENDOC(DOCNAM,.TRUE.,NLET,LUNDOC,LUNDOCT,.TRUE., & 'DOCUMENT FILE', & .FALSE.,.FALSE.,.TRUE.,NEWFILE,IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (LUNDOCT .LE. 0) THEN CALL ERRT(101,'OPERATION CAN NOT USE INCORE DOC FILE.',NE) RETURN ENDIF ENDIF IF (NSLICE .GT. 1) THEN ALLOCATE(VIN3(NSAM,NROW,NSLICE),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'DI, VIN3',IER) RETURN ENDIF C LOAD INPUT VOLUME INDEX = 0 DO K = 1,NSLICE DO J = 1,NROW INDEX = INDEX + 1 CALL REDLIN(LUN1,VIN3(1,J,K),NSAM,INDEX) ENDDO ENDDO CALL EROSION3(VIN3,NSAM,NROW,NSLICE,LH,LAT, MODE,LUNDOC,LUN2, & NPIXER) DEALLOCATE(VIN3) WRITE(NOUT,*) ' VOXELS ERODED: ',NPIXER ELSE ALLOCATE(VIN2(NSAM,NROW),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'DI, VIN2',IER) RETURN ENDIF C LOAD INPUT VOLUME INDEX = 0 DO J = 1,NROW INDEX = INDEX + 1 CALL REDLIN(LUN1,VIN2(1,J),NSAM,INDEX) ENDDO CALL EROSION2(VIN2,NSAM,NROW,NSLICE,LH,LAT, MODE,LUNDOC,LUN2, & NPIXER) WRITE(NOUT,*) ' PIXELS ERODED: ',NPIXER DEALLOCATE(VIN2) ENDIF IF (LUNDOC .GT. 0) CLOSE(LUNDOC) END C ------------------------- EROSION2 ----------------------------- SUBROUTINE EROSION2(X,NSAM,NROW,NSLICE,LH,LAT,MODE,LUNDOC,LUN2, & NPIXER) REAL, DIMENSION(NSAM,NROW) :: X REAL, DIMENSION(NSAM) :: Y CHARACTER(LEN=1) :: MODE IKEY = 0 NPIXER = 0 DO J=1,NROW IF (MODE .EQ. 'C') THEN C "CROSS" DO I=1,NSAM C COPY UNERODED VOXEL VALUE Y(I) = X(I,J) IF (X(I,J) .GT. 0.0) THEN C CENTRAL VOXEL IS OCCUPIED (NON-ZERO) LB = 0 DO M=-LH,LH IF (M. NE. 0) THEN C COUNT NUMBER OF OCCUPIED NEIGHBORS IF (X(I,MOD(J+M+NROW-1,NROW)+1) .GT. 0.0) & LB = LB + 1 IF (X(MOD(I+M+NSAM-1,NSAM)+1,J) .GT. 0.0) & LB = LB + 1 IF (X(I,J).GT.0.0) LB = LB + 1 ENDIF ENDDO IF (LB .LT. LAT) THEN C ERODE CENTRAL VOXEL (SET IT TO ZERO) Y(I) = 0.0 IF (LUNDOC .NE. 0 .AND. LB .EQ. 0) THEN C RECORD ULTIMATE ERODED PIXEL LOCATION CALL EROSION_TODOC(LUNDOC,IKEY,I,J,1,IRT) ENDIF NPIXER = NPIXER + 1 ENDIF ENDIF ENDDO ELSE C "BOX" CONNECTIVITY DO I=1,NSAM Y(I) = X(I,J) IF (X(I,J) .GT. 0.0) THEN C CENTRAL VOXEL IS OCCUPIED (NON-ZERO) LB = 0 C COUNT NUMBER OF OCCUPIED NEIGHBORS DO MJ=-LH,LH MJM = MOD(J+MJ+NROW-1,NROW)+1 DO MI=-LH,LH IF (X(MOD(I+MI+NSAM-1,NSAM)+1,MJM) & .GT. 0.0) LB = LB + 1 ENDDO ENDDO C ADJUST FOR THE CENTRAL ELEMENT LB = LB - 1 IF (LB .LT. LAT) THEN C ERODE CENTRAL VOXEL (SET TO ZERO) Y(I) = 0.0 IF (LUNDOC .NE. 0 .AND. LB .EQ. 0) THEN C RECORD ULTIMATE ERODED PIXEL LOCATION CALL EROSION_TODOC(LUNDOC,IKEY,I,J,1,IRT) ENDIF NPIXER = NPIXER + 1 ENDIF ENDIF ENDDO ENDIF C OUTPUT VOLUME CALL WRTLIN(LUN2,Y,NSAM,J) ENDDO END C ------------------------- EROSION3 ----------------------------- SUBROUTINE EROSION3(X,NSAM,NROW,NSLICE,LH,LAT,MODE,LUNDOC,LUN2, & NPIXER) REAL, DIMENSION(NSAM,NROW,NSLICE) :: X REAL, DIMENSION(NSAM) :: Y CHARACTER(LEN=1) :: MODE IKEY = 0 NPIXER = 0 DO N=1,NSLICE DO J=1,NROW IF (MODE .EQ. 'C') THEN C "CROSS" DO I=1,NSAM C COPY UNERODED VOXEL VALUE Y(I) = X(I,J,N) IF (X(I,J,N) .GT. 0.0) THEN C CENTRAL VOXEL IS OCCUPIED (NON-ZERO) LB = 0 DO M=-LH,LH IF (M. NE. 0) THEN C COUNT NUMBER OF OCCUPIED NEIGHBORS IF (X(I,MOD(J+M+NROW-1,NROW)+1,N) .GT. 0.0) & LB = LB + 1 IF (X(MOD(I+M+NSAM-1,NSAM)+1,J,N) .GT. 0.0) & LB = LB + 1 IF(X(I,J,MOD(N+M+NSLICE-1,NSLICE)+1).GT.0.0) & LB = LB + 1 ENDIF ENDDO IF (LB .LT. LAT) THEN C ERODE CENTRAL VOXEL (SET IT TO ZERO) Y(I) = 0.0 IF (LUNDOC .NE. 0 .AND. LB .EQ. 0) THEN C RECORD ULTIMATE ERODED PIXEL LOCATION CALL EROSION_TODOC(LUNDOC,IKEY,I,J,N,IRT) ENDIF NPIXER = NPIXER + 1 ENDIF ENDIF ENDDO ELSE C "BOX" CONNECTIVITY DO I=1,NSAM Y(I) = X(I,J,N) IF (X(I,J,N) .GT. 0.0) THEN C CENTRAL VOXEL IS OCCUPIED (NON-ZERO) LB = 0 C COUNT NUMBER OF OCCUPIED NEIGHBORS DO MN=-LH,LH MNM = MOD(N+MN+NSLICE-1,NSLICE)+1 DO MJ=-LH,LH MJM = MOD(J+MJ+NROW-1,NROW)+1 DO MI=-LH,LH IF (X(MOD(I+MI+NSAM-1,NSAM)+1,MJM,MNM) & .GT. 0.0) LB = LB + 1 ENDDO ENDDO ENDDO C ADJUST FOR THE CENTRAL ELEMENT LB = LB - 1 IF (LB .LT. LAT) THEN C ERODE CENTRAL VOXEL (SET TO ZERO) Y(I) = 0.0 IF (LUNDOC .NE. 0 .AND. LB .EQ. 0) THEN C RECORD ULTIMATE ERODED PIXEL LOCATION CALL EROSION_TODOC(LUNDOC,IKEY,I,J,N,IRT) ENDIF NPIXER = NPIXER + 1 ENDIF ENDIF ENDDO ENDIF C OUTPUT VOLUME CALL WRTLIN(LUN2,Y,NSAM,J+(N-1)*NROW) ENDDO ENDDO END C ------------------------- EROSION_TODOC ---------------------- SUBROUTINE EROSION_TODOC(LUNDOC,IKEY,ISAM,IROW,ISLICE,IRTFLG) REAL, DIMENSION(3) :: DLIST IKEY = IKEY + 1 C PUSH DLIST INTO DOC. FILE NVAL = 3 DLIST(1) = ISAM DLIST(2) = IROW DLIST(3) = ISLICE CALL LUNDOCWRTDAT(LUNDOC,IKEY,DLIST,NVAL,IRTFLG) END