C @(#)filtfl.for 17.1.1.1 (ESO-DMD) 01/25/02 17:19:20 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 Massachusetts 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 PROGRAM FILTFL C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program FILTFL version 1.00 881010 C D. Baade ESO - Garching C 1.10 890518 C C.KEYWORDS C filtering, noise reduction C C.PURPOSE C (Partly) fill up pixels having a flux below a user-supplied threshold with C the flux above the threshold in neighbouring pixels. C C.ALGORITHM C Get the names of input + output frames from IN_A + OUT_A (both 2-D). C Search for pixels with flux below threshold. Search box around each C such pixel for pixels with flux above threshold. Search box around each C such high-flux pixel for additional low-flux pixels. Now redistribute C flux from high pixels into central low pixel in such a fashion that C a) all high pixels contribute the same amount to the central low pixel, C b) a given high pixel gives the same flux to all its neighboring low C pixels, and C c) no high pixel drops below the threshold. C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/60 input frame C OUT_A/C/1/60 output frame C INPUTI/I/1/2 x and y dimensions of search box C INPUTR/R/1/1 threshold C C.VERSIONS C C 010201 last modif C C -------------------------------------------------------------------- C IMPLICIT NONE C INTEGER STAT,NAXIS,NPIX(2),IAV INTEGER INNO,OUTNO INTEGER RADIUS(2),THRESH INTEGER UNI(1),NULO,MADRID(1) INTEGER*8 INPT,OUTPT C CHARACTER INFRM*60,OUTFRM*60,CUNIT*64,IDENT*72 C DOUBLE PRECISION START(2),STEP(2) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up MIDAS environment + enable automatic error abort CALL STSPRO('FILTFL') C C get name of input image + map frame CALL STKRDC('IN_A',1,1,60,IAV,INFRM,UNI,NULO,STAT) CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIX,START,STEP, + IDENT,CUNIT,INPT,INNO,STAT) IF (NAXIS.NE.2) + CALL STETER(1,'Currently only 2-D frames supported!') C C Get name of output frame, map it for writing CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,STAT) IF (OUTFRM.EQ.INFRM) ! in- and output frame must be different + CALL STETER + (2,'Names of input + output frames must be different!') CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START, + STEP,IDENT,CUNIT,OUTPT,OUTNO,STAT) C C copy input frame to output IAV = NPIX(1) * NPIX(2) CALL COPYF(MADRID(INPT),MADRID(OUTPT),IAV) C C Get dimensions of filter box and the filtering threshold CALL STKRDI('INPUTI',1,2,IAV,RADIUS,UNI,NULO,STAT) CALL STKRDR('INPUTR',1,1,IAV,THRESH,UNI,NULO,STAT) C C Now comes the real thing: CALL FILLPIT(MADRID(INPT),MADRID(OUTPT),NPIX(1),NPIX(2), + RADIUS(1),RADIUS(2),THRESH) C C Update descriptor HISTORY CALL DSCUPT(INNO,OUTNO,' ',STAT) C C Free data CALL STSEPI END SUBROUTINE FILLPIT(IN,OUT,NX,NY,RX,RY,TR) C IMPLICIT NONE C INTEGER NX,NY,RX,RY,DX,DY,I,J,K,L,M,N,IK,JL,IM,JN INTEGER XSTART,XEND,YSTART,YEND,NHI1,NLO2,N20,NPERC INTEGER NCHECK C REAL IN(NX,NY),OUT(NX,NY),TR,SURPLUS,TAKE,GIVE C DX=2*RX+1 DY=2*RY+1 XSTART=DX+1 XEND=NX-DX YSTART=DY+1 YEND=NY-DY C C Prepare for progress messages C N20=NX/5 IF (N20.LT.10) THEN NCHECK=999999 ! not worth the effort ELSE NCHECK=N20 NPERC=20 ENDIF C C First treat all pixels farther than DX (DY) from frame boundaries C in X (Y). C DO 2000 I=XSTART,XEND DO 1800 J=YSTART,YEND IF (IN(I,J).LT.TR) THEN NHI1=0 GIVE=1.E+12 C C Consider all pixels within NX x NY box centered on the pixel to be filled C DO 1600 K=1,DX IK=I+K-RX-1 DO 1400 L=1,DY JL=J+L-RY-1 IF (IN(IK,JL).GT.TR) THEN NHI1=NHI1+1 NLO2=0 C C For all pixels with a flux above the threshold, check if there are C other low-value pixels in their respective NX x NY neighborhood boxes: C DO 1200 M=1,DX IM=IK+M-RX-1 DO 1000 N=1,DY JN=JL+N-RY-1 IF (IN(IM,JN).LT.TR) NLO2=NLO2+1 1000 CONTINUE ! Counted all secondary 1200 CONTINUE ! low-flux pixels ... C C Pixel (IK,JL) can donate at most SURPLUS. However, it may have to divide C its donation potential among NLO2 low-flux pixels. Therefore, earmark only C an accordingly smaller fraction for donation to current low pixel (I,J). C Furthermore, no pixel shall give more than it can afford. C SURPLUS=IN(IK,JL)-TR IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2 IF (GIVE.GT.SURPLUS) GIVE=SURPLUS ! always give maximimum ENDIF ! ... around this high-flux pixel 1400 CONTINUE 1600 CONTINUE C C If no high-flux pixel has been found in the box, skip to next pixel in C input image. C IF (NHI1.EQ.0) GOTO 1800 C C All flux donations available for pixel (I,J) have been summed up, now C work out if they really have to give all they can. C TAKE=NHI1*GIVE IF ((TAKE+IN(I,J)).GT.TR) THEN TAKE=TR-IN(I,J) GIVE=TAKE/NHI1 ENDIF C C Do the donation: C DO 1700 K=1,DX IK=I+K-RX-1 DO 1660 L=1,DY JL=J+L-RY-1 IF (IN(IK,JL).GT.TR) + OUT(IK,JL)=OUT(IK,JL)-GIVE 1660 CONTINUE ! all donations 1700 CONTINUE ! made OUT(I,J)=OUT(I,J)+TAKE ! and received ENDIF ! this low pixel filled C 1800 CONTINUE ! all inner pixels C C Issue progress messages C IF (I.GE.NCHECK) CALL PROGRS (20,N20,NPERC,NCHECK) C 2000 CONTINUE ! dealt with C C Next, consider outer pixels close to frame boundaries. Begin at left: C XSTART=1 XEND=DX C DO 6000 I=XSTART,XEND DO 4000 J=YSTART,YEND IF (IN(I,J).LT.TR) THEN NHI1=0 GIVE=1.E+12 DO 3600 K=1,DX IK=I+K-RX-1 IF (IK.LT.1) GOTO 3600 C DO 3400 L=1,DY JL=J+L-RY-1 IF (IN(IK,JL).GT.TR) THEN NHI1=NHI1+1 NLO2=0 DO 3200 M=1,DX IM=IK+M-RX-1 IF (IM.LT.1) GOTO 3200 DO 3000 N=1,DY JN=JL+N-RY-1 IF (IN(IM,JN).LT.TR) NLO2=NLO2+1 3000 CONTINUE ! Counted all secondary 3200 CONTINUE ! low-flux pixels ... SURPLUS=IN(IK,JL)-TR IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2 IF (SURPLUS.LT.GIVE) GIVE=SURPLUS ENDIF ! ... around this high-flux pixel 3400 CONTINUE 3600 CONTINUE C IF (NHI1.EQ.0) GOTO 4000 C TAKE=NHI1*GIVE IF ((TAKE+IN(I,J)).GT.TR) THEN TAKE=TR-IN(I,J) GIVE=TAKE/NHI1 ENDIF C DO 3800 K=1,DX IK=I+K-RX-1 IF (IK.LT.1) GOTO 3800 DO 3660 L=1,DY JL=J+L-RY-1 IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE 3660 CONTINUE ! all donations 3800 CONTINUE ! made OUT(I,J)=OUT(I,J)+TAKE ! and received ENDIF ! this low pixel filled C 4000 CONTINUE 6000 CONTINUE C C Continue at right edge: C XSTART=NX-DX+1 XEND=NX C DO 8000 I=XSTART,XEND DO 7800 J=YSTART,YEND IF (IN(I,J).LT.TR) THEN NHI1=0 GIVE=1.E+12 DO 7600 K=1,DX IK=I+K-RX-1 IF (IK.GT.NX) GOTO 7600 C DO 7400 L=1,DY JL=J+L-RY-1 C IF (IN(IK,JL).GT.TR) THEN NHI1=NHI1+1 NLO2=0 DO 7200 M=1,DX IM=IK+M-RX-1 IF (IM.GT.NX) GOTO 7200 C DO 7000 N=1,DY JN=JL+N-RY-1 IF (IN(IM,JN).LT.TR) NLO2=NLO2+1 7000 CONTINUE ! Counted all secondary 7200 CONTINUE ! low-flux pixels ... SURPLUS=IN(IK,JL)-TR IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2 IF (SURPLUS.LT.GIVE) GIVE=SURPLUS ENDIF ! ... around this high-flux pixel 7400 CONTINUE 7600 CONTINUE IF (NHI1.EQ.0) GOTO 7800 C TAKE=NHI1*GIVE IF ((TAKE+IN(I,J)).GT.TR) THEN TAKE=TR-IN(I,J) GIVE=TAKE/NHI1 ENDIF DO 7700 K=1,DX IK=I+K-RX-1 IF (IK.GT.NX) GOTO 7700 C DO 7660 L=1,DY JL=J+L-RY-1 IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE 7660 CONTINUE ! all donations 7700 CONTINUE ! made OUT(I,J)=OUT(I,J)+TAKE ! and received ENDIF ! this low pixel filled C 7800 CONTINUE 8000 CONTINUE C C Switch to bottom: C XSTART=DX+1 XEND=NX-DX YSTART=1 YEND=DY C DO 9000 I=XSTART,XEND DO 8900 J=YSTART,YEND C IF (IN(I,J).LT.TR) THEN NHI1=0 GIVE=1.E+12 C DO 8500 K=1,DX IK=I+K-RX-1 DO 8400 L=1,DY JL=J+L-RY-1 IF (JL.LT.1) GOTO 8400 C IF (IN(IK,JL).GT.TR) THEN NHI1=NHI1+1 NLO2=0 DO 8300 M=1,DX IM=IK+M-RX-1 DO 8200 N=1,DY JN=JL+N-RY-1 IF (JN.LT.1) GOTO 8200 IF (IN(IM,JN).LT.TR) NLO2=NLO2+1 8200 CONTINUE ! Counted all secondary 8300 CONTINUE ! low-flux pixels ... SURPLUS=IN(IK,JL)-TR IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2 IF (SURPLUS.LT.GIVE) GIVE=SURPLUS ENDIF ! ... around this high-flux pixel 8400 CONTINUE C 8500 CONTINUE C IF (NHI1.EQ.0) GOTO 8900 TAKE=NHI1*GIVE IF ((TAKE+IN(I,J)).GT.TR) THEN TAKE=TR-IN(I,J) GIVE=TAKE/NHI1 ENDIF C DO 8700 K=1,DX IK=I+K-RX-1 DO 8600 L=1,DY JL=J+L-RY-1 IF (JL.LT.1) GOTO 8600 IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE 8600 CONTINUE ! all donations 8700 CONTINUE ! made OUT(I,J)=OUT(I,J)+TAKE ! and received ENDIF ! this low pixel filled 8900 CONTINUE C 9000 CONTINUE C C Finish at top: C YSTART=NY-DY+1 YEND=NY C DO 10000 I=XSTART,XEND DO 9900 J=YSTART,YEND IF (IN(I,J).LT.TR) THEN NHI1=0 GIVE=1.E+12 DO 9500 K=1,DX IK=I+K-RX-1 DO 9400 L=1,DY JL=J+L-RY-1 IF (JL.GT.NY) GOTO 9400 IF (IN(IK,JL).GT.TR) THEN NHI1=NHI1+1 NLO2=0 DO 9300 M=1,DX IM=IK+M-RX-1 DO 9200 N=1,DY JN=JL+N-RY-1 IF (JN.GT.NY) GOTO 9200 IF (IN(IM,JN).LT.TR) NLO2=NLO2+1 9200 CONTINUE ! Counted all secondary 9300 CONTINUE ! low-flux pixels ... SURPLUS=IN(IK,JL)-TR IF (NLO2.GT.1) SURPLUS=SURPLUS/NLO2 IF (SURPLUS.LT.GIVE) GIVE=SURPLUS ENDIF ! ... around this high-flux pixel 9400 CONTINUE 9500 CONTINUE C IF (NHI1.EQ.0) GOTO 9900 TAKE=NHI1*GIVE IF ((TAKE+IN(I,J)).GT.TR) THEN TAKE=TR-IN(I,J) GIVE=TAKE/NHI1 ENDIF C DO 9700 K=1,DX IK=I+K-RX-1 DO 9600 L=1,DY JL=J+L-RY-1 IF (JL.GT.NY) GOTO 9600 IF (IN(IK,JL).GT.TR) OUT(IK,JL)=OUT(IK,JL)-GIVE 9600 CONTINUE ! all donations 9700 CONTINUE ! made OUT(I,J)=OUT(I,J)+TAKE ! and received ENDIF ! this low pixel filled 9900 CONTINUE C 10000 CONTINUE C RETURN END