C @(#)smosub1.for 17.1.1.1 (ES0-DMD) 01/25/02 17:40:01 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) 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 C the 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 C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C SUBROUTINE AVER1M(METH,A,B,NX,NY,RADIUS,TR,ALL,V,CUTS) C C ---------------------------------------------------------- C C handle median filtering with all different possibilities C C ---------------------------------------------------------- C IMPLICIT NONE C INTEGER NX,NY,ALL INTEGER ICOUNT,KCOUNT,INDX INTEGER NH,NL,NN,NO,NOZ INTEGER NP,NXMRX,KK,K INTEGER OFF,RX,RX1,RXM1 INTEGER RADIUS(2),NW1,NW2,BINDX,NAUX C REAL A(*),B(*),TR(*),CUTS(*) REAL AA,TEST,COUNT,W,THRESH,THRB,RR C REAL V(*) C CHARACTER*(*) METH C C initialize C BINDX = 1 RX = RADIUS(1) NW1 = RX*2 + 1 !no. of values on row of filter NW2 = 1 !no. of values on column of filter ICOUNT = NW1 COUNT = FLOAT(ICOUNT) !no. of pixels in neighbourh. with center NXMRX = NX - RX RX1 = RX + 1 RXM1 = RX - 1 C OFF = 0 !offset for center row THRESH = TR(1) C IF (ALL.EQ.0) THEN NH = ICOUNT/2 !index of median NOZ = RX*NW2 + 1 !index of center value KCOUNT = ICOUNT - 1 ELSE NH = (ICOUNT+1) / 2 !index of median KCOUNT = ICOUNT ENDIF NAUX = NH/2 + 1 C IF (METH(2:2).EQ.'Q') GOTO 5000 IF (METH(2:2).EQ.'Z') GOTO 7000 C C ------------------------- C IF ( (METH(2:2).NE.'F') .AND. + (THRESH.LE.0.) .AND. + (ALL.NE.0.) ) THEN !the fastest we can do C DO 300,NL=1,NY C C build up 1st V buffer INDX = OFF + RX1 !pixel to work on NO = 1 !init work array index... C DO 80, NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 80 CONTINUE CALL SORTIT(V,KCOUNT,NH,NAUX) C B(BINDX) = V(NH) BINDX = BINDX + 1 C C for following pixels check against V buffer before sorting it all C DO 200,NP=RX1+1,NXMRX INDX = OFF + NP !pixel to work on C IF (A(INDX-RX-1) .GT. V(NH)) THEN !outgoing pixel > median? AA = A(INDX+RX) IF (AA .GE. V(NH)) GOTO 180 !new pixel > median => O.K. C DO 150, K=NH-1,1,-1 IF (V(K) .LE. AA) THEN DO 100, KK=NH,K+1,-1 V(KK+1) = V(KK) 100 CONTINUE V(K+1) = AA !insert new pixel in sorted V GOTO 180 ENDIF 150 CONTINUE C DO 160, K=NH-1,1,-1 !AA is the smallest pixel in V V(K+1) = V(K) 160 CONTINUE V(1) = AA C ELSE !we have to fill + sort the V buffer NO = 1 DO 170, NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 170 CONTINUE CALL SORTIT(V,KCOUNT,NH,NAUX) ENDIF C 180 B(BINDX) = V(NH) BINDX = BINDX + 1 C 200 CONTINUE C OFF = OFF + NX 300 CONTINUE C CALL CHCUTS(B,BINDX,CUTS) RETURN ENDIF C C DO 4000 NL=1,NY C C fork according to threshold option (just for speed reasons) IF (METH(2:2).EQ.'F') GOTO 3000 IF (THRESH.GT.0.) GOTO 1000 C C extract values in the inner part DO 600,NP=RX1,NXMRX INDX = OFF + NP !pixel to work on NO = 1 !init work array index... C DO 400,NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 400 CONTINUE C IF (ALL.EQ.0) V(NOZ) = V(ICOUNT) !replace center value by last one CALL SORTIT(V,KCOUNT,NH,NAUX) C B(BINDX) = V(NH) BINDX = BINDX + 1 600 CONTINUE GOTO 3700 !continue as other part C C extract values in the inner part 1000 IF (METH(2:2).NE.'A') GOTO 2000 C C handle absolute threshold DO 1900,NP=RX1,NXMRX INDX = OFF + NP !pixel to work on W = A(INDX) !save original pixel value NO = 1 !init work array index... C DO 1400,NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 1400 CONTINUE C IF (ALL.EQ.0) V(NOZ) = V(ICOUNT) !replace center value by last one CALL SORTIT(V,KCOUNT,NH,NAUX) TEST = V(NH) C C compare meadian of neighbourhood with threshold IF (ABS(TEST-W).LE.THRESH) TEST = W B(BINDX) = TEST BINDX = BINDX + 1 1900 CONTINUE GOTO 3700 C C handle relative threshold 2000 DO 2900,NP=RX1,NXMRX INDX = OFF + NP !pixel to work on W = A(INDX) !save original pixel value NO = 1 !init work array index... C DO 2400,NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 2400 CONTINUE C IF (ALL.EQ.0) V(NOZ) = V(ICOUNT) !replace center value by last one CALL SORTIT(V,KCOUNT,NH,NAUX) TEST = V(NH) C C compare median of neighbourhood with pixel in question IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W B(BINDX) = TEST BINDX = BINDX + 1 2900 CONTINUE GOTO 3700 C C handle fixed threshold-interval 3000 DO 3600,NP=RX1,NXMRX INDX = OFF + NP !pixel to work on TEST = A(INDX) !save original pixel value IF ( (TEST.LT.TR(1)) .OR. (TEST.GT.TR(2)) ) THEN B(BINDX) = TEST ELSE NO = 1 !init work array index... DO 3400,NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 3400 CONTINUE IF (ALL.EQ.0) V(NOZ) = V(ICOUNT) CALL SORTIT(V,KCOUNT,NH,NAUX) B(BINDX) = V(NH) ENDIF C BINDX = BINDX + 1 3600 CONTINUE C 3700 OFF = OFF + NX !update offset... C 4000 CONTINUE !outer loop over all lines C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C C handle relative threshold C 5000 DO 6000,NL=1,NY DO 5900,NP=RX1,NXMRX INDX = OFF + NP !pixel to work on W = A(INDX) !save original pixel value NO = 1 !init work array index... C DO 5400,NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 5400 CONTINUE C IF (ALL.EQ.0) V(NOZ) = V(ICOUNT) !replace center value by last one CALL SORTIT(V,KCOUNT,NH,NAUX) TEST = V(NH) C C compare median of neighbourhood with pixel in question IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) TEST = W B(BINDX) = TEST BINDX = BINDX + 1 5900 CONTINUE C C check progress OFF = OFF + NX !update offset... 6000 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C C handle threshold = a + b*SQRT(median) C 7000 THRB = TR(2) DO 8000,NL=1,NY DO 7900,NP=RX1,NXMRX INDX = OFF + NP !pixel to work on W = A(INDX) !save original pixel value NO = 1 !init work array index... C DO 7400,NN=INDX-RX,INDX+RX V(NO) = A(NN) NO = NO + 1 7400 CONTINUE C IF (ALL.EQ.0) V(NOZ) = V(ICOUNT) !replace center value by last one CALL SORTIT(V,KCOUNT,NH,NAUX) TEST = V(NH) C C compare a + b*SQRT(median) of neighbourhood with pixel in question IF (TEST.LT.0.0) THEN RR = SQRT(-TEST) ELSE RR = SQRT(TEST) ENDIF IF (ABS(TEST-W).LE.(THRESH+(THRB*RR))) TEST = W B(BINDX) = TEST BINDX = BINDX + 1 7900 CONTINUE C C check progress OFF = OFF + NX !update offset... 8000 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C END SUBROUTINE AVER1A(METH,A,B,NX,NY,RADIUS,TR,ALL,V,CUTS) C C -------------------------- C C we handle boxcar filtering (smoothing, averaging) C C -------------------------- C IMPLICIT NONE C INTEGER NX,NY,ALL INTEGER ICOUNT,INDX INTEGER NL,NO,NOZ INTEGER NP,NXMRX,NXMRX1 INTEGER OFF,RX,RX1,RXM1,TOP INTEGER RADIUS(2),NW1,NW2,BINDX C REAL A(*),B(*),TR(*),CUTS(*) REAL TEST,COUNT,FACTOR,W,THRESH C DOUBLE PRECISION V(*) C CHARACTER*(*) METH C C initialize C BINDX = 1 RX = RADIUS(1) NW1 = RX*2 + 1 !no. of values on row of filter NW2 = 1 !no. of values on column of filter ICOUNT = NW1 COUNT = FLOAT(ICOUNT) !no. of pixels in neighbourh. with center NXMRX = NX - RX NXMRX1 = NXMRX - 1 RX1 = RX + 1 RXM1 = RX - 1 OFF = 0 !offset for center row THRESH = TR(1) C IF (ALL.EQ.0) THEN FACTOR = 1. / (COUNT-1.) ELSE FACTOR = 1. / COUNT ENDIF C NOZ = 0 !offset for lower left NO = 1 !offset for lower left + 1 C C fork according to test criterion (just for speed reasons) IF (METH(2:2).EQ.'F') GOTO 4000 IF (METH(2:2).EQ.'Q') GOTO 5000 IF (THRESH.GT.0.) THEN !check, if absolute or relative thresholding IF (METH(2:2).NE.'A') THEN GOTO 3000 ELSE GOTO 2000 ENDIF ENDIF C C ........... C C here average without thresholding for each line C C ........... C C DO 1500, NL=1,NY !for row-average in 2-dim frame V(1) = 0.D0 !build new high row sum DO 100, NP=NOZ+1,NOZ+NW1 V(1) = V(1) + A(NP) 100 CONTINUE NOZ = NOZ + NX !follow with offset... C C average in the inner part TOP = NO !init to left corner of interval C IF (ALL.EQ.0) THEN DO 400, INDX=OFF+RX1,OFF+NXMRX1 B(BINDX) = (V(1) - A(INDX)) * FACTOR BINDX = BINDX + 1 C C sum over neighbourhood + update the row-sum on the fly V(1) = V(1) - A(TOP) + A(TOP+NW1) TOP = TOP + 1 400 CONTINUE B(BINDX) = (V(1) - A(OFF+NXMRX)) * FACTOR OFF = OFF + NX !update offset of center line ELSE DO 600, INDX=RX1,NXMRX1 B(BINDX) = V(1) * FACTOR BINDX = BINDX + 1 V(1) = V(1) - A(TOP) + A(TOP+NW1) TOP = TOP + 1 600 CONTINUE B(BINDX) = V(1) * FACTOR !avoid out_of_bounds in TOP+NW1 ENDIF C BINDX = BINDX + 1 NO = NO + NX !follow with initial offset... 1500 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C C C ........... C C here average with absolute thresholding for each line C C ........... C C 2000 DO 2900, NL=1,NY V(1) = 0.D0 DO 2050,NP=NOZ+1,NOZ+NW1 V(1) = V(1) + A(NP) 2050 CONTINUE NOZ = NOZ + NX !follow with offset... C C average in the inner part TOP = NO !init to left corner of interval C IF (ALL.EQ.0) THEN DO 2200 INDX=OFF+RX1,OFF+NXMRX1 W = A(INDX) !save original pixel value C C sum over neighbourhood + update the row-sum on the fly TEST = (V(1) - W) * FACTOR V(1) = V(1) - A(TOP) + A(TOP+NW1) C C compare average of neighbourhood with threshold IF (ABS(TEST-W).LE.THRESH) THEN B(BINDX) = W ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 TOP = TOP + 1 2200 CONTINUE W = A(OFF+NXMRX) !for last pixel no A(TOP+NW1) TEST = (V(1) - W) * FACTOR ELSE DO 2300 INDX=OFF+RX1,OFF+NXMRX1 W = A(INDX) !save original pixel value TEST = V(1) * FACTOR V(1) = V(1) - A(TOP) + A(TOP+NW1) C C compare average of neighbourhood with threshold IF (ABS(TEST-W).LE.THRESH) THEN B(BINDX) = W ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 TOP = TOP + 1 2300 CONTINUE W = A(OFF+NXMRX) !for last pixel no A(TOP+NW1) TEST = V(1) * FACTOR ENDIF C IF (ABS(TEST-W).LE.THRESH) TEST = W B(BINDX) = TEST BINDX = BINDX + 1 C NO = NO + NX !follow with initial offset... OFF = OFF + NX !update offset of center line 2900 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C C ........... C C here average with relative thresholding for each line C C ........... C C 3000 DO 3900,NL=1,NY V(1) = 0.D0 DO 3050,NP=NOZ+1,NOZ+NW1 V(1) = V(1) + A(NP) 3050 CONTINUE NOZ = NOZ + NX !follow with offset... C C average in the inner part TOP = NO !init to lower left corner of matrix C IF (ALL.EQ.0) THEN DO 3200,INDX=OFF+RX1,OFF+NXMRX1 W = A(INDX) !save original pixel value C C sum over neighbourhood + update the row sum on the fly TEST = (V(1) - W) * FACTOR V(1) = V(1) - A(TOP) + A(TOP+NW1) C C compare average of neighbourhood with threshold IF (ABS(TEST-W).LE.ABS(THRESH*W)) THEN B(BINDX) = W ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 TOP = TOP + 1 3200 CONTINUE W = A(OFF+NXMRX) !for last pixel no A(TOP+NW1) TEST = (V(1) - W) * FACTOR ELSE DO 3300,INDX=OFF+RX1,OFF+NXMRX1 W = A(INDX) !save original pixel value TEST = V(1) * FACTOR V(1) = V(1) - A(TOP) + A(TOP+NW1) IF (ABS(TEST-W).LE.ABS(THRESH*W)) THEN B(BINDX) = W ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 TOP = TOP + 1 3300 CONTINUE W = A(OFF+NXMRX) !for last pixel no A(TOP+NW1) TEST = V(1) * FACTOR ENDIF C IF (ABS(TEST-W).LE.ABS(THRESH*W)) TEST = W B(BINDX) = TEST BINDX = BINDX + 1 C NO = NO + NX !follow with initial offset... OFF = OFF + NX !update offset of center line 3900 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C C C ........... C C here average with fixed threshold-interval for each line C C ........... C C 4000 DO 4900,NL=1,NY V(1) = 0.D0 DO 4050,NP=NOZ+1,NOZ+NW1 V(1) = V(1) + A(NP) 4050 CONTINUE NOZ = NOZ + NX !follow with offset... C C average in the inner part TOP = NO !init to lower left corner of matrix C IF (ALL.EQ.0) THEN DO 4200,INDX=OFF+RX1,OFF+NXMRX1 TEST = A(INDX) !save original pixel value C C compare average of neighbourhood with threshold IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) ) THEN B(BINDX) = (V(1) - TEST) * FACTOR ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 C C sum over neighbourhood + update the row-sum on the fly V(1) = V(1) - A(TOP) + A(TOP+NW1) TOP = TOP + 1 4200 CONTINUE TEST = A(OFF+NXMRX) !for last pixel no A(TOP+NW1) IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) ) + TEST = (V(1) - TEST) * FACTOR ELSE DO 4300,INDX=OFF+RX1,OFF+NXMRX1 TEST = A(INDX) !save original pixel value IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) ) THEN B(BINDX) = V(1) * FACTOR ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 V(1) = V(1) - A(TOP) + A(TOP+NW1) TOP = TOP + 1 4300 CONTINUE TEST = A(OFF+NXMRX) !for last pixel no A(TOP+NW1) IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) ) + TEST = V(1) * FACTOR ENDIF C B(BINDX) = TEST BINDX = BINDX + 1 C NO = NO + NX !follow with initial offset... OFF = OFF + NX !update offset of center line 4900 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C C ........... C C here average with relative thresholding for each line C C ........... C C 5000 IF (ALL.EQ.0) THEN FACTOR = 1. / (COUNT-1.) ELSE FACTOR = 1. / COUNT ENDIF C DO 5900,NL=1,NY V(1) = 0.D0 DO 5050,NP=NOZ+1,NOZ+NW1 V(1) = V(1) + A(NP) 5050 CONTINUE NOZ = NOZ + NX !follow with offset... C C average in the inner part TOP = NO !init to lower left corner of matrix C IF (ALL.EQ.0) THEN DO 5200,INDX=OFF+RX1,OFF+NXMRX1 W = A(INDX) !save original pixel value C C sum over neighbourhood + update the row sum on the fly TEST = (V(1) - W) * FACTOR V(1) = V(1) - A(TOP) + A(TOP+NW1) C C compare average of neighbourhood with threshold IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) THEN B(BINDX) = W ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 TOP = TOP + 1 5200 CONTINUE C ELSE DO 5300,INDX=OFF+RX1,OFF+NXMRX1 W = A(INDX) !save original pixel value TEST = V(1) * FACTOR V(1) = V(1) - A(TOP) + A(TOP+NW1) IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) THEN B(BINDX) = W ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 TOP = TOP + 1 5300 CONTINUE W = A(OFF+NXMRX) !for last pixel no A(TOP+NW1) TEST = V(1) * FACTOR ENDIF C IF (ABS(TEST-W).LE.ABS(THRESH*TEST)) TEST = W B(BINDX) = TEST BINDX = BINDX + 1 C NO = NO + NX !follow with initial offset... OFF = OFF + NX !update offset of center line 5900 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C END SUBROUTINE AVER2M(METH,A,B,NX,NY,RADIUS,TR,ALL,O,VB,CUTS) C C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine AVER2M version 1.00 861215 C K. Banse ESO - Garching C 1.10 870303 1.20 880707 1.30 910529 C 1.40 920210 C C.KEYWORDS C median, neighbourhood C C.PURPOSE C remove spot noise C C.ALGORITHM C a) smooth image by taking the average A(x,y)of all points in neighbourhood C of x,y with given radius, and replace original value I(x,y) by A(x,y), C if average differs more than THRESH*I(x,y) from original value C b) as a) but take median instead of average C C.INPUT/OUTPUT C call as AVER2M(METH,A,B,NX,NY,RADIUS,TR,ALL,O,VB,CUTS) C C input par: C METH: char.exp. (1:1) = method to use: always M for median C (2:2) = A for absolute, C R for relative threshold C Q is handled in BVER2M ... C F for fixed threshold-interval C A: R*4 array input array C B: R*4 array output array C NX: I*4 no. of pixels per line C NY: I*4 no. of lines in input array C RADIUS: I*4 array "radius" of neighbourhood in x and y, C i.e. no. of pixels in horizontal and vertical C direction to include C TR: R*4 array threshold(s) for center point to determine C replacement C ALL: I*4 flag for inclusion of center point in C averaging process: C 1 = include center point, 0 = no inclusion C O: R*4 working buffer C VB: R*4 buffer for columns if x-radius = 0 C temp buffer otherwise C C.VERSIONS C 1.00 has been built on 2.80 of old module AVERAGE C 1.10 add fixed threshold-interval option C 1.20 fix bugs with fixed-threshold option (jump into inner loop at 8000) C 1.30 also filter borders (input image is larger now) and use C heapsort and binary insertion for median determination C 1.40 fix bug with METH(2:2)=F (fixed interval) C 1.50 work on faster algorithm for 1-dim column filters C C -------------------------------------------------------------------- C IMPLICIT NONE C INTEGER NX,NY,ALL INTEGER ICOUNT,KCOUNT INTEGER INDX,AUX1,AUX2,AUX3,AUX4,AUX5 INTEGER N,NH,NINDX,NL,NN,NO,NOZ INTEGER NP,NXMRX,MFLAG INTEGER OFF,RX,RX1,RX2,RY,RYNX INTEGER RADIUS(2),NW1,NW11,NW2,NAUX INTEGER JWOFF,KWOFF,KWOF1,KINDX,BINDX,WINDX C REAL A(*),B(*),TR(*) REAL O(*),VB(*),CUTS(*) REAL TEST,W,THRESH,THRB,QQ,RR REAL COLO(200) C CHARACTER*(*) METH CHARACTER*80 CBUF C C initialize C BINDX = 1 RX = RADIUS(1) RY = RADIUS(2) NW1 = RX*2 + 1 !no. of values on row of filter NW2 = RY*2 + 1 !no. of values on column of filter ICOUNT = NW1 * NW2 RYNX = RY * NX NXMRX = NX - RX RX1 = RX + 1 AUX1 = RYNX + RX AUX2 = NX - NW1 C OFF = RYNX !offset for center row THRESH = TR(1) C C IF (ALL.EQ.0) THEN !center pixel not included NH = ICOUNT/2 !index of median NOZ = (RX * NW2) + RY + 1 !index of center value KCOUNT = ICOUNT - 1 !one pixel less NAUX = NH/2 + 1 !needed for sorting algorithm ELSE KCOUNT = ICOUNT NH = (ICOUNT+1) / 2 !index of median NAUX = NH/2 + 1 !needed for sorting algorithm ENDIF C !we need a size check (ugly...) IF ((KCOUNT.GT.2500) .OR. (NW1.GT.200)) THEN WRITE (CBUF,99778) KCOUNT CALL STTPUT(CBUF,N) RETURN ENDIF C C branch according to different options IF (METH(2:2).EQ.'F') GOTO 8000 !fixed interval elsewhere C MFLAG = 0 IF (METH(2:2).EQ.'Z') THEN THRB = TR(2) MFLAG = 3 !a+b*SQRT(median) elsewhere ELSE IF (METH(2:2).EQ.'Q') THEN MFLAG = 4 !relative threshold with median ELSE IF (THRESH.GT.0.) THEN IF (METH(2:2).NE.'A') THEN !relative thr > 0. MFLAG = 2 ELSE !absolute thr > 0. MFLAG = 1 ENDIF ENDIF ENDIF C IF (ALL.EQ.0) GOTO 6000 !exclude the center pixel C C --------------------------------------------------------------- C order the complete filter kernel completeley only for first C center pixel + first row C then just update the ordered set with the new pixels which are C added to the already ordered set C --------------------------------------------------------------- C OFF = OFF - AUX1 !saves 1 addition... NW11 = NW1 - 1 !additional constants AUX3 = RX1 + (NW2 - 1) * (NW1 + AUX2) !for this section AUX4 = KCOUNT - NW1 + 1 AUX5 = NXMRX + NW11 RX2 = RX1 + NW1 C C ......................... C C include the center pixel C C ......................... C NO = 1 !first line + first kernel NINDX = OFF + RX1 DO 150, NN=1,NW2 DO 100, N=1,NW1 O(NO) = A(NINDX) VB(NO) = O(NO) !save 1. neighbourhood NO = NO + 1 NINDX = NINDX + 1 100 CONTINUE NINDX = NINDX + AUX2 != NINDX + NX - NW1 150 CONTINUE C CALL XSORT1(O,KCOUNT,NH,KINDX) CALL XSAVX(1,KCOUNT) !copy internal sort arrays out C C compare median of neighbourhood with threshold IF (MFLAG.GT.0) THEN !threshold > 0 WINDX = OFF + AUX1 + RX1 !index for input pixel we work on W = A(WINDX) !save original pixel value TEST = O(KINDX) RR = ABS(TEST-W) IF (MFLAG.EQ.1) THEN !absolute threshold IF (RR.LE.THRESH) TEST = W ELSE IF (MFLAG.EQ.2) THEN !relative threshold IF (RR.LE.ABS(THRESH*W)) TEST = W ELSE IF (MFLAG.EQ.3) THEN IF (TEST.LT.0.0) THEN QQ = SQRT(-TEST) ELSE QQ = SQRT(TEST) ENDIF IF (RR.LE.(THRESH+(THRB*QQ))) TEST = W ELSE IF (RR.LE.ABS(THRESH*TEST)) TEST = W ENDIF B(BINDX) = TEST ELSE B(BINDX) = O(KINDX) ENDIF BINDX = BINDX + 1 C C now move along first line JWOFF = 1 DO 250, NP=RX2,AUX5 NINDX = OFF + NP DO 200, N=1,NW2 COLO(N) = A(NINDX) !new values NINDX = NINDX + NX 200 CONTINUE C CALL XSORT2(O,KCOUNT,NH,JWOFF,COLO,NW1,NW2,KINDX) JWOFF = JWOFF + 1 !point to next column start IF (JWOFF.GT.NW1) JWOFF = 1 C C compare median of neighbourhood with threshold IF (MFLAG.GT.0) THEN !threshold > 0 WINDX = WINDX + 1 W = A(WINDX) !save original pixel value TEST = O(KINDX) RR = ABS(TEST-W) IF (MFLAG.EQ.1) THEN !absolute threshold IF (RR.LE.THRESH) TEST = W ELSE IF (MFLAG.EQ.2) THEN !relative threshold IF (RR.LE.ABS(THRESH*W)) TEST = W ELSE IF (MFLAG.EQ.3) THEN IF (TEST.LT.0.0) THEN QQ = SQRT(-TEST) ELSE QQ = SQRT(TEST) ENDIF IF (RR.LE.(THRESH+(THRB*QQ))) TEST = W ELSE IF (RR.LE.ABS(THRESH*TEST)) TEST = W ENDIF B(BINDX) = TEST ELSE B(BINDX) = O(KINDX) ENDIF BINDX = BINDX + 1 250 CONTINUE OFF = OFF + NX !update offset... KWOFF = 1 KWOF1 = 0 C C use last 1. buffer for beginning of all following lines IF ((RY+2).GT.(NY-RY)) GOTO 3000 C DO 1000, NL=RY+2,NY-RY DO 330, N=1,KCOUNT !copy 1. neighbourhood O(N) = VB(N) 330 CONTINUE CALL XSAVX(-1,KCOUNT) !copy internal sort arrays back in C NINDX = OFF + AUX3 DO 400, N=1,NW1 COLO(N) = A(NINDX) VB(N+KWOF1) = COLO(N) NINDX = NINDX + 1 400 CONTINUE C CALL XSORT2(O,KCOUNT,NH,KWOFF,COLO,1,NW1,KINDX) CALL XSAVX(1,KCOUNT) !copy internal sort arrays out KWOFF = KWOFF + NW1 !point to next row start IF (KWOFF.GT.AUX4) KWOFF = 1 KWOF1 = KWOFF - 1 C C compare median of neighbourhood with threshold IF (MFLAG.GT.0) THEN !threshold > 0 WINDX = OFF + AUX1 + RX1 W = A(WINDX) !save original pixel value TEST = O(KINDX) RR = ABS(TEST-W) IF (MFLAG.EQ.1) THEN !absolute threshold IF (RR.LE.THRESH) TEST = W ELSE IF (MFLAG.EQ.2) THEN !relative threshold IF (RR.LE.ABS(THRESH*W)) TEST = W ELSE IF (MFLAG.EQ.3) THEN IF (TEST.LT.0.0) THEN QQ = SQRT(-TEST) ELSE QQ = SQRT(TEST) ENDIF IF (RR.LE.(THRESH+(THRB*QQ))) TEST = W ELSE IF (RR.LE.ABS(THRESH*TEST)) TEST = W ENDIF B(BINDX) = TEST ELSE B(BINDX) = O(KINDX) ENDIF BINDX = BINDX + 1 C C now move along each line JWOFF = 1 DO 800, NP=RX2,AUX5 !AUX5 = NXMRX + NW11 NINDX = OFF + NP DO 700, N=1,NW2 COLO(N) = A(NINDX) !new values NINDX = NINDX + NX 700 CONTINUE C CALL XSORT2(O,KCOUNT,NH,JWOFF,COLO,NW1,NW2,KINDX) JWOFF = JWOFF + 1 IF (JWOFF.GT.NW1) JWOFF = 1 C C compare median of neighbourhood with threshold IF (MFLAG.GT.0) THEN !threshold > 0 WINDX = WINDX + 1 W = A(WINDX) !save original pixel value TEST = O(KINDX) RR = ABS(TEST-W) IF (MFLAG.EQ.1) THEN !absolute threshold IF (RR.LE.THRESH) TEST = W ELSE IF (MFLAG.EQ.2) THEN !relative threshold IF (RR.LE.ABS(THRESH*W)) TEST = W ELSE IF (MFLAG.EQ.3) THEN IF (TEST.LT.0.0) THEN QQ = SQRT(-TEST) ELSE QQ = SQRT(TEST) ENDIF IF (RR.LE.(THRESH+(THRB*QQ))) TEST = W ELSE IF (RR.LE.ABS(THRESH*TEST)) TEST = W ENDIF B(BINDX) = TEST ELSE B(BINDX) = O(KINDX) ENDIF BINDX = BINDX + 1 800 CONTINUE C OFF = OFF + NX !update offset... 1000 CONTINUE C C finally update the cuts 3000 CALL CHCUTS(B,BINDX,CUTS) RETURN C C ......................... C C exclude the center pixel C C ......................... C 6000 DO 6900,NL=RY+1,NY-RY DO 6800,NP=RX1,NXMRX NO = 1 !extract neighbourhood INDX = OFF + NP !pixel to work on NINDX = INDX - AUX1 != INDX - RYNX + RX DO 6400,NN=1,NW2 DO 6350,N=1,NW1 O(NO) = A(NINDX) NO = NO + 1 NINDX = NINDX + 1 6350 CONTINUE NINDX = NINDX + AUX2 != NINDX + NX - NW1 6400 CONTINUE C O(NOZ) = O(ICOUNT) !replace center value by last one CALL SORTIT(O,KCOUNT,NH,NAUX) IF (MFLAG.GT.0) THEN W = A(INDX) !save original pixel value TEST = O(NH) RR = ABS(TEST-W) IF (MFLAG.EQ.1) THEN !absolute threshold IF (RR.LE.THRESH) TEST = W ELSE IF (MFLAG.EQ.2) THEN !relative threshold IF (RR.LE.ABS(THRESH*W)) TEST = W ELSE IF (MFLAG.EQ.3) THEN IF (TEST.LT.0.0) THEN QQ = SQRT(-TEST) ELSE QQ = SQRT(TEST) ENDIF IF (RR.LE.(THRESH+(THRB*QQ))) TEST = W ELSE IF (RR.LE.ABS(THRESH*TEST)) TEST = W ENDIF B(BINDX) = TEST ELSE B(BINDX) = O(NH) ENDIF BINDX = BINDX + 1 6800 CONTINUE C OFF = OFF + NX !update offset... 6900 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C C C --------------------------------------------------------------- C handle fixed threshold-interval C --------------------------------------------------------------- C C 8000 DO 8900,NL=RY+1,NY-RY DO 8800,NP=RX1,NXMRX INDX = OFF + NP !pixel to work on TEST = A(INDX) !save original pixel value C C check if pixel should be replaced IF ( (TEST.GE.TR(1)) .AND. (TEST.LE.TR(2)) ) THEN NO = 1 !extract neighbourhood NINDX = INDX - AUX1 != INDX - RYNX - RX DO 8400,NN=1,NW2 DO 8350,N=1,NW1 O(NO) = A(NINDX) NO = NO + 1 NINDX = NINDX + 1 8350 CONTINUE NINDX = NINDX + AUX2 != NINDX + NX - NW1 8400 CONTINUE C IF (ALL.EQ.0) + O(NOZ) = O(ICOUNT) !replace cent_val by last one CALL SORTIT(O,KCOUNT,NH,NAUX) B(BINDX) = O(NH) ELSE B(BINDX) = TEST ENDIF BINDX = BINDX + 1 8800 CONTINUE C OFF = OFF + NX !update offset... 8900 CONTINUE C C finally update the cuts CALL CHCUTS(B,BINDX,CUTS) RETURN C 99778 FORMAT('Filter size = ',I5,' > 2500 - we abort...') C END