C ++******************************************************************** C * C WATERSHED K WAS WRONG MAR 01 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 MEDIAN(LUN1,LUN2,NSAM,NROW,NSLICE) * C * C PURPOSE: MEDIAN FILTRATION USING BOX OR CROSS PATTERN KERNAL * C * C PARAMETERS: * C * C IMAGE_PROCESSING_ROUTINE * C * C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** SUBROUTINE WATERSHED(LUN1,LUN2,NSAM,NROW,NSLICE,FMINT) INCLUDE 'CMBLOCK.INC' REAL, ALLOCATABLE, DIMENSION(:) :: VOL,VOLOUT NPIX = NSAM * NROW * NSLICE ALLOCATE (VOL(NPIX),VOLOUT(NPIX),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NPIX * 2 CALL ERRT(46,'CE WA, VOL ',MWANT) RETURN ENDIF IF (NSLICE .EQ. 1) THEN C FOR IMAGE C READ INPUT VOLUME CALL REDVOL(LUN1,NSAM,NROW,1,1,VOL,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 CALL WATER2(VOL,VOLOUT,NSAM,NROW,FMINT,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 C CREATE OUTPUT VOLUME CALL WRTVOL(LUN2,NSAM,NROW,1,1,VOLOUT,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ELSE C FOR VOLUME CALL ERRT(101,'NOT IMPLEMENT FOR VOLUMES',NE) GOTO 9999 ENDIF 9999 IF (ALLOCATED(VOL)) DEALLOCATE(VOL) IF (ALLOCATED(VOLOUT)) DEALLOCATE(VOLOUT) END C ----------------- WATER2 ------------------------------------- SUBROUTINE WATER2(XIMG,XOUT,NSAM,NROW,FMINT,IRTFLG) REAL, DIMENSION(NSAM*NROW) :: XIMG,XOUT REAL, ALLOCATABLE, DIMENSION(:) :: XSORT,RLOC REAL, ALLOCATABLE, DIMENSION(:) :: NEXTPIX INCLUDE 'CMBLOCK.INC' IRTFLG = 1 NPIX = NSAM * NROW NQUESIZE = NPIX ALLOCATE (XSORT(NPIX),RLOC(NPIX),NEXTPIX(NQUESIZE), & STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NPIX * 3 + NQUESIZE CALL ERRT(46,'CE WA, RLOC,... ',MWANT) RETURN ENDIF DO IPIX = 1,NPIX XSORT(IPIX) = XIMG(IPIX) RLOC(IPIX) = IPIX ENDDO C SORT THE PIXELS (SORT IS A SINGLETON SORT), RLOC IS REORDERD CALL SORT(XSORT,XOUT,RLOC,NPIX) C SET CURRENT SEGMENTATION LEVEL RLEVEL = 1.0 NLEVEL = 0 C INITIALIZE ARRAY XOUT XOUT = -1.0 DO IPIX = 1,NPIX ! LOOP OVER ALL SORTED PIXELS ILOC = RLOC(IPIX) ! STARTING PIXEL LOCATION C UPDATE LEVEL IF (NLEVEL .GT. 0) THEN RLEVEL = RLEVEL + 1.0 NLEVEL = 0 DO I=1,NPIX IF(XOUT(I) .EQ. 0) XOUT(I) = -1 ENDDO ENDIF C PUT STARTING PIXEL IN QUE IQUE = 0 NQUE = 1 NEXTPIX(NQUE) = ILOC C write(6,*)' ----------- stack( ',nstack,'): ',ILOC DO WHILE(IQUE .LT. NQUE) C GET NEXT PIXEL FROM QUE IQUE = IQUE + 1 ISLOC = NEXTPIX(IQUE) c write(6,*)' from stack( ',istack,') got: ',ISLOC IF (XOUT(ISLOC) .LE. 0.0) THEN C PIXEL NOT VISITED OR SET YET C GIVE UNVISITED PIXEL CURRENT WATERSHED LEVEL XOUT(ISLOC) = RLEVEL c write(6,*)' XOUT(',ISLOC,'): ',RLEVEL C INCREMENT NUMBER OF PIXELS IN AT THIS LEVEL NLEVEL = NLEVEL + 1 C PLACE PIXEL'S UNVISITED HIGHER NEIGHBORS IN QUE C UPPER NEIGH = ISLOC - NSAM IF (NEIGH .GT. 0 .AND. & XOUT(NEIGH) .LT. 0.0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,' QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF C LOWER NEIGH = ISLOC + NSAM IF (NEIGH .LE. NPIX .AND. & XOUT(NEIGH) .LT. 0.0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,'QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF C LEFT NEIGH = ISLOC - 1 IF (MOD(NEIGH,NSAM) .GT. 0 .AND. & XOUT(NEIGH) .LT. 0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,'QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF C RIGHT NEIGH = ISLOC + 1 IF (MOD(ISLOC,NSAM) .GT. 0 .AND. & XOUT(NEIGH) .LT. 0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,'QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF C UPPER-LEFT NEIGH = ISLOC - NSAM - 1 IF (NEIGH .GT. 0 .AND. & MOD(NEIGH,NSAM) .GT. 0 .AND. & XOUT(NEIGH) .LT. 0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,'QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF C LOWER-LEFT NEIGH = ISLOC + NSAM - 1 IF (NEIGH .LE. NPIX .AND. & MOD(NEIGH,NSAM) .GT. 0 .AND. & XOUT(NEIGH) .LT. 0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,'QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF C UPPER-RIGHT NEIGH = ISLOC - NSAM + 1 IF (NEIGH .GT. 0 .AND. & MOD(ISLOC,NSAM) .GT. 0 .AND. & XOUT(NEIGH) .LT. 0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,'QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF C LOWER-RIGHT NEIGH = ISLOC + NSAM + 1 IF (NEIGH .LE. NPIX .AND. & MOD(ISLOC,NSAM) .GT. 0 .AND. & XOUT(NEIGH) .LT. 0 .AND. & (XIMG(NEIGH) .GE. XIMG(ISLOC))) THEN NQUE = NQUE + 1 IF (NQUE .GT. NQUESIZE) THEN CALL ERRT(102,'QUE OVERFLOW',NQUE) RETURN ENDIF NEXTPIX(NQUE) = NEIGH XOUT(NEIGH) = 0.0 ENDIF ENDIF ! END OF: IF (XOUT(ISLOC) .LE. 0 .AND...... ENDDO ! END OF: DO WHILE (IQUE .LE. NQUE) NSHEDS = RLEVEL IF (NLEVEL .GT. 0) THEN WRITE(NOUT,90)NSHEDS,NLEVEL 90 FORMAT(' WATERSHED: ',I8,' HAS: ',I8,' PIXELS') ENDIF ENDDO ! END OF: IPIX = 1,NPIX LOOP OVER SORTED PIXELS NSHEDS = RLEVEL IF (NLEVEL .LE. 0) NSHEDS = NSHEDS -1 WRITE(NOUT,*)' ' WRITE(NOUT,*)' WATERSHEDS: ',NSHEDS IRTFLG = 0 END