C @(#)modif.for 17.1.1.1 (ESO-DMD) 01/25/02 17:40:00 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 MODIF C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program MODIF version 1.00 890630 C K. Banse ESO - Garching C C.KEYWORDS C modify pixels C C.PURPOSE C merge of MODPIX.FOR, MODAR.FOR, MODCUT.FOR C C.ALGORITHM C see the subroutines C C.INPUT/OUTPUT C the following keywords are used: C ACTION/C/1/2 action = MP for modify/pixel C = MA for modify/area C = MC for modify/cuts C C.VERSIONS C 1.00 created C see SCCS C C 010706 last modif C C------------------------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,STAT INTEGER UNI(1),NULO C CHARACTER ACTION*2 C C get into MIDAS environment CALL STSPRO('MODIF') C CALL STKRDC('ACTION',1,1,2,IAV,ACTION,UNI,NULO,STAT) CALL UPCAS(ACTION,ACTION) C C branch according to option IF (ACTION .EQ. 'MP') THEN CALL MODPIX(1) ELSE IF (ACTION .EQ. 'MA') THEN CALL MODPIX(2) ELSE IF (ACTION .EQ. 'MC') THEN CALL STKRDC('P2',1,1,1,IAV,ACTION,UNI,NULO,STAT) CALL MODCUT(ACTION) ENDIF C CALL STSEPI END SUBROUTINE MODPIX(MODMOD) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine MODPIX version 1.00 881028 C K. Banse ESO - Garching C 1.50 980923 S.Wolf C C.KEYWORDS C bulk data frame, bad pixels C C.PURPOSE C to replace a group of undefined or bad pixels values by surface C interpolation of neighboring area. The bad pixel area is defined via the C cursor. C C.ALGORITHM C a) (original input from H. Waldthausen, had to be cleaned a lot...) C a surface is fitted using the pixels surrounding the cursor area C the surrounding area is chosen as C x-pixels = x-area factor * x-size of cursor rectangle C y-pixels = y-area factor * y-size of cursor rectangle C the program forms an approximation to the weighted, least-squares Chebyshev C series surface fit; the weights are derived from the Poisson error. C C b) (original input from A. Lauberts) C a LSQ surface is fitted using the pixels around a circle but inside C the cursor rectangle C C.INPUT/OUTPUT C call as MODPIX(MODMOD) C C MODMOD: Integer 1 = Modify/Pixel C 2 = Modify/Area C C the following keywords are used: C a) C P1/C/1/60 input frame, table C or C CURSOR C IN_A/C/1/60 input frame, if P1 = CURSOR C OUT_A/C/1/60 result frame C INPUTI/I/1/6 area factors in x,y C x,y degree of fitting polynomial, C no. of iterations + overlay_draw_flag C NOISE/C/1/1 Y or N - whether to add noise C IMASEC/C/1/100 allowed image section C C b) C P1/C/1/80 input frame, table C or C CURSOR C IN_A/C/1/60 input frame, if P1 = CURSOR C OUT_A/C/1/60 result frame C INPUTI/I/1/2 degree of fitted SURFACs, draw flag C INPUTR/R/1/1 constant in case that INPUTI = 0 C C------------------------------------------------------------------- C IMPLICIT NONE C INTEGER MODMOD INTEGER ICX,ICY,IPX,IPY INTEGER ISX,ISY,IT,XDEG,YDEG INTEGER IPIXS(4),ISUR(4) INTEGER DEGREE,IAV,IDRAW INTEGER N,N2,NAXIS,NCOL,NROW,NMAL INTEGER*8 PNTRA,PNTRB INTEGER IMNOA,IMNOB,SIZE,LFLAGS(2) INTEGER STAT1,STAT2,STAT,TABFLG INTEGER NPIXA(3),KCOOR(4),WPIX(3),ISCAL(2) INTEGER SPIX(3),EPIX(3) INTEGER IXLO,IXUP,IYLO,IYUP,INPUTI(6),DAZOP INTEGER ICUR1(2),ICUR2(2),WW(4) INTEGER TID,TABCOLN(4) INTEGER LOW1(3),LOW2(3),HI(3),CUTIND,COOFF INTEGER UNI(1),NUL,MADRID(1) INTEGER XFIGU(6),YFIGU(8) C CHARACTER FRAMEA*80,TABLE*80,FRAMEB*80 CHARACTER*80 OUTPUT CHARACTER*100 CIMASEC CHARACTER CUNITA*64,IDENTA*72,CACT*8,CNOISE*1 CHARACTER TABLAB(4)*16 CHARACTER*40 ERMS12,ERMS13,ERMS15,ERMS16 CHARACTER INFO*44,TIP15*56,TIP12*28 CHARACTER CSYST C DOUBLE PRECISION STARTA(3),STEPA(3) DOUBLE PRECISION DD1(4),DD2(4) !for WCS conversions C REAL RBUF(6),CUTS(4),COOR(4),STPINV(2) REAL RCUTS(2),ZCUTS(2),PCUR1(6),PCUR2(6),CONST REAL OFFSET,OMIN,CRMS,PCEN(6) REAL INPUTR(8) C LOGICAL NOISE,TABNUL(4),SELFLG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C INCLUDE 'MID_INCLUDE:IDIDEV.INC' INCLUDE 'MID_INCLUDE:IDIMEM.INC' C COMMON /VMR/ MADRID C C messages DATA + INFO /' switch cursors on - next time we exit ... '/, + ERMS12 /' ERR: surface dimension too large '/, + ERMS13 /' ERR: in mathematical subroutines '/, + ERMS15 /' ERR: surface-column undersampled '/, + ERMS16 /' ERR: area not fully in image section '/ DATA + TIP12 /' TIP: reduce cursor area '/, + TIP15 /' TIP: enlarge cursor area, reduce degree and #iter.'/ C DATA OUTPUT /' '/, TABLE /' '/, FRAMEA /' '/, FRAMEB /' '/ DATA TABLAB /'XSTART ','YSTART ','XEND ','YEND '/ DATA KCOOR /-1,-1,-1,-1/ DATA LOW1 /1,1,1/, LOW2 /1,1,1/, HI /1,1,1/ DATA NPIXA /1,1,1/, WPIX /1,1,1/ DATA ISCAL /1,1/, DAZOP /1/, COOFF /0/ DATA CUNITA /' '/, IDENTA /' '/ DATA SPIX /1,1,0/, EPIX /1,1,0/ DATA PCUR1 /6*0.0/, PCUR2 /6*0.0/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get input option (CURSOR or table) + read keys CSYST = 'W' !coord.system default: world LFLAGS(1) = -99 CALL STKRDC('P1',1,1,80,IAV,OUTPUT,UNI,NUL,STAT) N = INDEX(OUTPUT,',') IF (N.GT.0) THEN !table specified TABFLG = 2 FRAMEA(1:80) = OUTPUT(1:N-1)//' ' CALL CLNFRA(FRAMEA,FRAMEA,0) TABLE(1:80) = OUTPUT(N+1:)//' ' N2 = INDEX(TABLE,',') IF (N2.GT.0) THEN !only for MODI/PIX implem. CSYST = TABLE(N2+1:N2+2) CALL UPCAS(CSYST,CSYST) TABLE(1:80) = TABLE(1:N2-1)//' ' IF (CSYST.EQ.'P') THEN !pixel coord. selected DO 90, N=1,4 N2 = INDEX(TABLAB(N),' ') !labels have to end TABLAB(N)(1:) = TABLAB(N)(1:N2-1)//'PIX ' 90 CONTINUE !with PIX ELSE IF (CSYST.NE.'W') THEN !Wrong flag used OUTPUT(1:) = 'Wrong table flag used. Only PIXEL or'// + ' WORLD is possible ' CALL STETER(3,OUTPUT) ENDIF ELSE CSYST = 'W' !world coordinates ENDIF CALL CLNTAB(TABLE,TABLE,0) CALL TBTOPN(TABLE,F_I_MODE,TID,STAT) CALL TBIGET(TID,NCOL,NROW,N,N,N,STAT) !get no. of columns + rows C DO 100, N=1,4 CALL TBLSER(TID,TABLAB(N),TABCOLN(N),STAT) !find columns... IF (TABCOLN(N).LE.0) THEN OUTPUT(1:) = 'column labelled '//TABLAB(N)// + 'missing... ' CALL STETER(1,OUTPUT) ENDIF 100 CONTINUE ELSE C TABFLG = 1 CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NUL,STAT) ENDIF C C get result frame + integer parameters CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEB,UNI,NUL,STAT) CALL STKRDI('INPUTI',1,6,IAV,INPUTI,UNI,NUL,STAT) CALL STKRDR('INPUTR',1,3,IAV,INPUTR,UNI,NUL,STAT) C C get remaining parameters depending on method IF (MODMOD.EQ.1) THEN !get the noise flag !(RNH) CALL STKRDC('NOISE',1,1,1,IAV,CNOISE,UNI,NUL,STAT) IF (CNOISE.EQ.'Y' .OR. CNOISE.EQ.'y') THEN NOISE=.TRUE. ELSE NOISE=.FALSE. ENDIF !get allowed image section CALL STKRDC('IMASEC',1,1,100,IAV,CIMASEC,UNI,NUL,STAT) C C verify area factors, degree of polynomial + no. of iterations INPUTR(1) = MAX(1.0,INPUTR(1)) !arfacts specified by INPUTR(2) = MAX(1.0,INPUTR(2)) !... MODI/PIX INPUTI(1) = MAX(1,INPUTI(1)) !only used by MODI/AREA INPUTI(2) = MAX(1,INPUTI(2)) !only used by MODI/AREA INPUTI(3) = MIN(10,MAX(0,INPUTI(3))) INPUTI(4) = MIN(10,MAX(0,INPUTI(4))) INPUTI(5) = MAX(0,INPUTI(5)) IDRAW = INPUTI(6) ELSE C C verify degree (only relevant for MODI/AREA) DEGREE = INPUTI(1) IF (DEGREE.LT.0) THEN DEGREE = 0 ELSE IF (DEGREE.GT.2) THEN DEGREE = 2 ENDIF IDRAW = INPUTI(2) IF (DEGREE.EQ.0) !for degree 0 get constant value + CONST = INPUTR(1) !changed by swolf@eso.org c + CALL STKRDR('INPUTR',1,1,IAV,CONST,UNI,NUL,STAT) ENDIF C C map input + result frame IF (FRAMEB(1:1).NE.'+') THEN CALL GENEQF(FRAMEA,FRAMEB,STAT) ELSE FRAMEB(1:) = FRAMEA(1:) !result frame = input frame STAT = 1 ENDIF IF (STAT.EQ.1) THEN CALL STIGET(FRAMEA,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE, + 2,NAXIS,NPIXA,STARTA,STEPA,IDENTA, + CUNITA,PNTRA,IMNOA,STAT) IMNOB = IMNOA PNTRB = PNTRA SIZE = NPIXA(1)*NPIXA(2) IF (NAXIS.GT.2) THEN NAXIS = 2 CALL STTPUT('Only 1. plane of frame used!',STAT) ENDIF CALL STTPUT + ('Warning: Input frame will be updated...',STAT) ELSE CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIXA,STARTA,STEPA,IDENTA, + CUNITA,PNTRA,IMNOA,STAT) IF (NAXIS.GT.2) THEN NAXIS = 2 CALL STTPUT('Only 1. plane of input frame used!',STAT) ENDIF CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIXA,STARTA,STEPA,IDENTA, + CUNITA,PNTRB,IMNOB,STAT) SIZE = NPIXA(1)*NPIXA(2) OUTPUT(1:) = '... copying input frame to result frame ... ' CALL STTDIS(OUTPUT,99,STAT) CALL COPYF(MADRID(PNTRA),MADRID(PNTRB),SIZE) IF (IDRAW .GE. 2) THEN !copy to disk IPIXS(1) = 1 !IPIXS will be reset later on... IPIXS(2) = 1 CALL MEMDSK(IMNOB,0,MADRID(PNTRB),NPIXA, + IPIXS,NPIXA,STAT) ENDIF ENDIF LFLAGS(2) = D_R4_FORMAT !indicate real data type C C get cuts of input image CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NUL,STAT) IF (CUTS(1).GE.CUTS(2)) THEN CUTIND = 3 ELSE CUTIND = 1 ENDIF ZCUTS(1) = CUTS(3) ZCUTS(2) = CUTS(4) STPINV(1) = 1.0 / STEPA(1) STPINV(2) = 1.0 / STEPA(2) C IF (MODMOD.EQ.1) THEN IF (CIMASEC(1:1).NE.'+') THEN CALL EXTCOO(IMNOA,CIMASEC,3,IAV,SPIX,EPIX,STAT) IF (STAT.NE.0) + CALL STETER(2,'Invalid coords for image section...') C ELSE SPIX(1) = 1 EPIX(1) = NPIXA(1) SPIX(2) = 1 EPIX(2) = NPIXA(2) ENDIF IXLO = SPIX(1) !lower x-coordinate IXUP = EPIX(1) !upper x-coordinate IYLO = SPIX(2) !lower y-coordinate IYUP = EPIX(2) !upper y-coordinate of allowed image section OMIN = CUTS(3) IF (OMIN.LE.0.) THEN !set offset value... OFFSET = ABS(OMIN) + 1. ELSE OFFSET = 0. ENDIF ENDIF C C fork according to input mode NMAL = 0 IF (TABFLG.EQ.2) THEN IF (IDRAW.GE.2) THEN CALL DTOPEN(1,STAT) !check, if image is loaded CALL DTGICH(QDSPNO,QIMCH,OUTPUT,INPUTR,STAT) CALL GENEQF(FRAMEA,OUTPUT,STAT) IF (STAT.NE.1) THEN CALL STTPUT + ('Input image not lodaded => draw flag ignored...', + STAT) IDRAW = 0 ELSE CALL PIXXCV('INIT',IMNOA,RBUF,STAT) ENDIF ELSE DAZOP = 0 !indicate, that we don't use display ENDIF IF (CSYST.EQ.'W') !init wc --> frame pixel conversion + CALL FPXWCO(0,IMNOA,DD1,DD2,STAT) GOTO 3000 ENDIF C IF ((IDRAW.EQ.1).OR.(IDRAW.EQ.3)) THEN CACT = 'NYYY?C0 ' ELSE CACT = 'NNYY?C0 ' ENDIF C C for interactive (cursor) mode - attach ImageDisplay CALL DTOPEN(1,STAT) FRAMEA(1:) = ' ' !to init GETCUR ... CALL SETCUR(QDSPNO,2,0,0,KCOOR,STAT) !set cursor rectangle CALL CONCHA(QDSPNO,QOVCH,1,0) !clear overlay channel CALL PIXXCV('INIT',IMNOA,RBUF,STAT) C C input cursor positions and check status 1000 WRITE(OUTPUT,10000) NMAL+1 CALL STTDIS(OUTPUT,99,STAT) CALL GETCUR(CACT,FRAMEA, + ICUR1,PCUR1,COOR(1),RBUF(1),STAT1, + ICUR2,PCUR2,COOR(3),RBUF(1),STAT2) C C check cursor status IF ((STAT1.EQ.0).AND.(STAT2.EQ.0)) THEN IF ((NMAL.EQ.0) .AND. (COOFF.EQ.0)) THEN COOFF = 1 FRAMEA(1:) = ' ' CALL STTPUT(INFO,STAT) GOTO 1000 ELSE GOTO 5000 ENDIF ENDIF C C now goto section, where we really do the work NMAL = NMAL + 1 !update counter IF (MODMOD.EQ.1) THEN PCEN(1) = (ICUR1(1)+ICUR2(1)) * .5 PCEN(2) = (ICUR1(2)+ICUR2(2)) * .5 CALL PIXXCV('IRW',0,PCEN,STAT) ENDIF GOTO 4000 C C here we get our input data from a table C 3000 NMAL = NMAL + 1 IF (NMAL.GT.NROW) THEN CALL TBTCLO(TID,STAT) WRITE(OUTPUT,10001) (NMAL-1) CALL STTPUT(OUTPUT,STAT) GOTO 5000 ENDIF C C read xstart, ystart, xend, yend from table CALL TBSGET(TID,NMAL,SELFLG,STAT) IF (.NOT.SELFLG) GOTO 3000 C CALL TBRRDR(TID,NMAL,4,TABCOLN,COOR,TABNUL,STAT) IF (CSYST.EQ.'W') THEN !world coordinates --> transfer to pixels DD1(1) = COOR(1) !convert wc -> fp DD1(2) = COOR(2) CALL FPXWCO(-1,0,DD1,DD2,STAT) IF (STAT .GT. 0) + CALL STETER(36,'Problems with WCS of frame...') PCUR1(1) = DD2(1) PCUR1(2) = DD2(2) DD1(1) = COOR(3) DD1(2) = COOR(4) CALL FPXWCO(-1,0,DD1,DD2,STAT) IF (STAT .GT. 0) + CALL STETER(36,'Problems with WCS of frame...') PCUR2(1) = DD2(1) PCUR2(2) = DD2(2) ELSE !table entries area already in pixels PCUR1(1) = COOR(1) PCUR1(2) = COOR(2) PCUR2(1) = COOR(3) PCUR2(2) = COOR(4) END IF IF (IDRAW.GE.2) CALL PIXXCV('_RS',0,PCUR1,STAT) ICUR1(1) = NINT(PCUR1(5)) ICUR1(2) = NINT(PCUR1(6)) C C do the replacement 4000 IF (MODMOD.EQ.1) THEN !define cursor area + surrounding surface IPIXS(1) = NINT(PCUR1(1)) !real pixel no.s IPIXS(2) = NINT(PCUR2(1)) IPIXS(3) = NINT(PCUR1(2)) IPIXS(4) = NINT(PCUR2(2)) IF ((IPIXS(2).LT.IXLO) .OR. (IPIXS(1).GT.IXUP) .OR. + (IPIXS(4).LT.IYLO) .OR. (IPIXS(3).GT.IYUP)) THEN CALL STTPUT(ERMS16,STAT) GOTO (1000,3000),TABFLG ENDIF IPX = IPIXS(2) - IPIXS(1) + 1 !no. of pixels in cursor area... IPY = IPIXS(4) - IPIXS(3) + 1 ICX = IPIXS(1) + IPX/2 !center of cursor area ICY = IPIXS(3) + IPY/2 C ISX = NINT(INPUTR(1)*(IPX-1)) + 1 !no. of pixels in surface ISY = NINT(INPUTR(2)*(IPY-1)) + 1 ISUR(1) = MAX(IXLO,ICX - ISX/2) !xlow ISUR(2) = MIN(IXUP,ISUR(1) + ISX - 1) !xhi ISUR(3) = MAX(IYLO,ICY - ISY/2) !ylow ISUR(4) = MIN(IYUP,ISUR(3) + ISY - 1) !yhi C XDEG = INPUTI(3) !erase spot YDEG = INPUTI(4) IT = INPUTI(5) CRMS = -1. C CALL ERASE(MADRID(PNTRB),NPIXA(1),NPIXA(2),IPIXS,ISUR,XDEG, + YDEG,OFFSET,IT,CRMS,NOISE,STAT) C IF (STAT.NE.0) THEN !check return status IF (STAT.GE.1 .AND. STAT.LE.5) THEN !NAGLIB error WRITE(OUTPUT,111) STAT CALL STTPUT(OUTPUT,STAT) CALL STTPUT(ERMS13,STAT) GOTO 5500 ELSE IF (STAT .EQ. 7) THEN !surface dimension too large CALL STTPUT(ERMS12,STAT) CALL STTPUT(TIP12,STAT) GOTO (1000,3000),TABFLG ELSE IF (STAT.EQ.8) THEN !subimage-row undersampled CALL STTPUT(ERMS15,STAT) CALL STTPUT(TIP15,STAT) GOTO (1000,3000),TABFLG ELSE IF (STAT.EQ.9) THEN !subimage-not well defined CALL STTPUT(ERMS16,STAT) GOTO (1000,3000),TABFLG ENDIF ENDIF C PCEN(5) = (COOR(1)+COOR(3)) * .5 PCEN(6) = (COOR(2)+COOR(4)) * .5 WRITE(OUTPUT,222) IPX,IPY,PCEN(5),PCEN(6) CALL STTPUT(OUTPUT,STAT) C ELSE CALL STAREM(MADRID(PNTRB),NPIXA(1),NPIXA(2), + STARTA,STEPA,COOR,CSYST,DEGREE,CONST,RCUTS) IF (ZCUTS(1).GT.RCUTS(1)) ZCUTS(1) = RCUTS(1) IF (ZCUTS(2).LT.RCUTS(2)) ZCUTS(2) = RCUTS(2) ENDIF C C if required, reload modified subwindow IF (IDRAW.GE.2) THEN LOW1(1) = NINT(PCUR1(1)) LOW1(2) = NINT(PCUR1(2)) HI(1) = NINT(PCUR2(1)) HI(2) = NINT(PCUR2(2)) WPIX(1) = HI(1) - LOW1(1) + 1 WPIX(2) = HI(2) - LOW1(2) + 1 CALL MEMDSK(IMNOB,0,MADRID(PNTRB),NPIXA,LOW1,WPIX,STAT) WW(1) = ICUR1(1) WW(2) = ICUR1(2) WW(3) = SCALX WW(4) = SCALY CALL LOADWN(LFLAGS,IMNOB,NPIXA,LOW1,WPIX,WW,CUTS(CUTIND)) ENDIF C IF ((IDRAW.EQ.1) .OR. (IDRAW.EQ.3)) THEN WW(1) = ICUR1(1) WW(2) = ICUR1(2) WPIX(1) = HI(1) - LOW1(1) + 1 WPIX(2) = HI(2) - LOW1(2) + 1 IF (SCALX.GT.1) THEN WPIX(1) = WPIX(1)/SCALX ELSE IF (SCALX.LT.-1) THEN N2 = -SCALX WPIX(1) = WPIX(1)*N2 ENDIF IF (SCALY.GT.1) THEN WPIX(2) = WPIX(2)/SCALY ELSE IF (SCALY.LT.-1) THEN N2 = -SCALY WPIX(2) = WPIX(2)*N2 ENDIF WW(3) = WW(1) + WPIX(1) - 1 WW(4) = WW(2)+ WPIX(2) - 1 CALL BLDGRA('REC',WW,INPUTR,XFIGU,YFIGU,6,N2) CALL IIGPLY(QDSPNO,QOVCH,XFIGU,YFIGU,N2,255,1,STAT) ENDIF C IF (TABFLG.EQ.1) THEN GOTO 1000 ELSE GOTO 3000 ENDIF C C get new dynamic range + write descriptor LHCUTS 5000 IF (CUTS(3).GT.ZCUTS(1)) CUTS(3) = ZCUTS(1) IF (CUTS(4).LT.ZCUTS(2)) CUTS(4) = ZCUTS(2) C C append HISTORY CALL DSCUPT(IMNOA,IMNOB,' ',STAT) CALL STDWRR(IMNOB,'LHCUTS',CUTS,1,4,UNI,STAT) !but use old cuts... C C That's it folks... 5500 IF (DAZOP.EQ.1) THEN CALL DTCLOS(QDSPNO) IF (IDRAW.GT.0) CALL REFOVR(STAT) ENDIF OUTPUT(1:) = ' ' CALL STTDIS(OUTPUT,99,STAT) RETURN C C format statements 10000 FORMAT(I3,', ready for cursor input...') 10001 FORMAT(I6,' rows processed.') 111 FORMAT('IFAIL = ',I3) 222 FORMAT('cleaned area:',I3,' x',I3,' pixels at xc = ',F10.3, + ', yc = ',F10.3) END SUBROUTINE + STAREM(A,NPIX1,NPIX2,STARTA,STEPA,COOR,CSYST,IDEG,CONST,RMAX) C IMPLICIT NONE C INTEGER NPIX1,NPIX2,IDEG INTEGER I,I1,I2,J,J1,J2 INTEGER M,N,INC INTEGER STAT C REAL A(NPIX1,NPIX2),COOR(4),RMAX(2) REAL X(200),Y(200),Z(200),B(7) REAL AIJ,B1,B2,B3,B4,B5,B6,D2 REAL Q2,Q22,QX,QX2,QY,QY2 REAL X1,X2,XC,XX,Y1,Y2,Y22,YC,YY REAL ZMAX,ZMIN REAL STEP1,STEP2,START1,START2 REAL CONST,STEP11,STEP22 C DOUBLE PRECISION STARTA(*),STEPA(*) C CHARACTER CSYST*1, OUTPUT*80 C 222 FORMAT('clean area:',I3,' x',I3,' pixels at xc = ',F10.3, + ', yc = ',F10.3) C C init STEP1=STEPA(1) STEP2=STEPA(2) START1=STARTA(1) START2=STARTA(2) C C compute X1=COOR(1) Y1=COOR(2) X2=COOR(3) Y2=COOR(4) C XC=(X1+X2)*0.5 YC=(Y1+Y2)*0.5 QX=X2-XC QY=Y2-YC QX2=QX**2 QY2=QY**2 Q2=QX*QY Q22=Q2**2 C C define outer rectangular boundary IF (CSYST.EQ.'W') THEN J1=(Y1-START2)/STEP2+1. J2=(Y2-START2)/STEP2+1. IF (J1.LT.1) THEN J1=1 ELSE IF (J2.GT.NPIX2) THEN J2=NPIX2 ENDIF C I1=(X1-START1)/STEP1+1. I2=(X2-START1)/STEP1+1. IF (I1.LT.1) THEN I1=1 ELSE IF (I2.GT.NPIX1) THEN I2=NPIX1 ENDIF Y1 = FLOAT(J1-1)*STEP2+START2-YC X1 = FLOAT(I1-1)*STEP1+START1-XC ELSE J1=NINT(Y1) J2=NINT(Y2) I1=NINT(X1) I2=NINT(X2) IF (J1.LT.1) THEN J1=1 ELSE IF (J2.GT.NPIX2) THEN J2=NPIX2 ENDIF IF (I1.LT.1) THEN I1=1 ELSE IF (I2.GT.NPIX1) THEN I2=NPIX1 ENDIF Y1 = Y1 - YC X1 = X1 - XC STEP1 = 1.0 STEP2 = 1.0 ENDIF C C assemble points between outer ellipse and circumscribed rectangle C XX and YY are coordinates relative to the centre (XC,YC) C YY = Y1 C WRITE(OUTPUT,222) I2-I1+1,J2-J1+1,XC,YC CALL STTPUT(OUTPUT,STAT) C C if degree = 0, replace intensities inside ellipse by constant IF (IDEG.EQ.0) THEN DO 400, J=J1,J2 Y22=YY*YY*QX2 XX=X1 DO 380, I=I1,I2 D2 = (XX*XX*QY2) + Y22 IF (D2.LE.Q22) THEN A(I,J)=CONST ENDIF XX=XX+STEP1 380 CONTINUE YY=YY+STEP2 400 CONTINUE C RMAX(1) = CONST RMAX(2) = CONST RETURN !we're done... ENDIF C C compute approx max no of points (the exact no depends on ellipse centre) C if necessary take coarser sample C M = (I2-I1+1)*(J2-J1+1) - NINT(3.14159*Q2/STEP1**2) !no. of points INC = 1 + (M/200) STEP11 = STEP1*FLOAT(INC) STEP22 = STEP2*FLOAT(INC) M = 0 ZMIN = 1.E10 ZMAX = -ZMIN DO 2000, J=J1,J2,INC Y22=YY*YY*QX2 XX=X1 DO 1000, I=I1,I2,INC D2 = (XX*XX*QY2) + Y22 IF(D2.LT.Q22) GOTO 990 M=M+1 X(M)=XX Y(M)=YY AIJ=A(I,J) Z(M)=AIJ ZMIN=AMIN1(ZMIN,AIJ) ZMAX=AMAX1(ZMAX,AIJ) IF(M.EQ.200) GOTO 3000 990 XX=XX+STEP11 1000 CONTINUE YY=YY+STEP22 2000 CONTINUE C C test M RMAX(1) = ZMIN RMAX(2) = ZMAX IF (M.LT.4) RETURN C C compute N IF ( (IDEG.EQ.1) .OR. + ((ZMAX-ZMIN).LT.(ZMAX+ZMIN)*0.2) ) THEN N=3 ELSE N=6 ENDIF C C fit max 2nd order SURFAC to M points C if M < 12, then use a plane only (degree = 1) C the 4 corner points should always be included C otherwise skip replacement C 3000 CALL SURFAC(X,Y,Z,B,M,N) B1 = B(1) IF (N.NE.3.AND.(B1.LT.ZMIN.OR.B1.GT.ZMAX)) THEN N=3 GOTO 3000 ENDIF C B2=B(2) B3=B(3) B4=B(4) B5=B(5) B6=B(6) C C replace intensities inside ellipse YY=Y1 DO 4000, J=J1,J2 Y22=YY*YY*QX2 XX=X1 DO 3800, I=I1,I2 D2 = (XX*XX*QY2) + Y22 IF (D2.LE.Q22) THEN A(I,J)=B1 + XX*(B2+XX*B4+YY*B5) + YY*(B3+YY*B6) IF (RMAX(1) .GT. A(I,J)) THEN RMAX(1) = A(I,J) ELSE IF (RMAX(2) .LT. A(I,J)) THEN RMAX(2) = A(I,J) ENDIF ENDIF XX=XX+STEP1 3800 CONTINUE YY=YY+STEP2 4000 CONTINUE C RETURN END SUBROUTINE SURFAC(X,Y,Z,B,M,N) C IMPLICIT NONE C INTEGER M,N INTEGER N1,K,I,J C DOUBLE PRECISION A(7,7),C(7) C REAL X(M),Y(M),Z(M),B(7) C C LSQ fitting of surface C Z = B1 + B2*X + B3*Y + B4*X*X + B5*X*Y + B6*Y*Y C to M points C IF (M.LT.12) N=3 N1 = N+1 DO 2, I=1,7 DO 1, J=1,7 A(I,J)=0.D0 1 CONTINUE B(I)=0. 2 CONTINUE C DO 5, K=1,M C(1)=1.D0 C(2)=DBLE(X(K)) C(3)=DBLE(Y(K)) C(N1)=DBLE(Z(K)) IF (N.EQ.6) THEN C(4)=C(2)**2 C(5)=C(2)*C(3) C(6)=C(3)**2 ENDIF C DO 4, I=1,N1 DO 3, J=1,N1 A(I,J)=A(I,J)+C(I)*C(J) 3 CONTINUE 4 CONTINUE 5 CONTINUE C CALL SIMUL(N,A,B) C RETURN END SUBROUTINE SIMUL(N,A,X) C C the GAUSS-JORDAN complete elimination method is employed C with the maximum pivot strategy C IMPLICIT NONE C INTEGER N,I,J,K,IROWI,IROWK,JCOLI,JCOLK INTEGER ISCAN,JSCAN,KM1,MAX INTEGER IROW(6),JCOL(6),STAT C DOUBLE PRECISION A(7,7),EPS,AIJCK,PIVOT,FPIVOT C REAL X(N) C DATA EPS /1.D-15/ DATA IROW /6*0/, JCOL /6*0/ C MAX=N+1 C C begin elimination procedure C DO 19,K=1,N KM1=K-1 PIVOT=0.D0 !search for the pivot element DO 10,I=1,N DO 11,J=1,N C C search irow and jcol arrays for invalid pivot subscripts IF (K.EQ.1) GOTO 9 DO 8,ISCAN=1,KM1 DO 7,JSCAN=1,KM1 IF (I.EQ.IROW(ISCAN)) GOTO 11 IF (J.EQ.JCOL(JSCAN)) GOTO 11 7 CONTINUE 8 CONTINUE 9 CONTINUE IF (DABS(A(I,J)).LE.DABS(PIVOT)) GOTO 11 PIVOT=A(I,J) IROW(K)=I JCOL(K)=J 11 CONTINUE 10 CONTINUE C C insure that selected pivot is larger than eps IF (DABS(PIVOT).LE.EPS) THEN CALL STTPUT(' Warning! Singular Matrix',STAT) RETURN ELSE FPIVOT = 1.D0 / PIVOT ENDIF C IROWK=IROW(K) JCOLK=JCOL(K) C C normalize pivot row elements DO 14,J=1,MAX A(IROWK,J)=A(IROWK,J)*FPIVOT 14 CONTINUE C C carry out elimination C A(IROWK,JCOLK)=FPIVOT DO 18,I=1,N AIJCK=A(I,JCOLK) IF (I.EQ.IROWK) GOTO 18 A(I,JCOLK) = -(AIJCK*FPIVOT) DO 17,J=1,MAX IF (J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J) 17 CONTINUE 18 CONTINUE 19 CONTINUE C C order solution values C DO 20,I=1,N IROWI=IROW(I) JCOLI=JCOL(I) X(JCOLI)= A(IROWI,MAX) 20 CONTINUE C RETURN END SUBROUTINE MODCUT(ACTIO) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine MODCUT version 1.00 890119 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frame, scaling C C.PURPOSE C to rescale subimages on the ImageDisplay C C.ALGORITM C extract subimage via cursor window + apply new cut values C C.INPUT/OUTPUT C ACTIO char. array 1 = C - cursor option, else full image C C the following keywords are used: C INPUTI/I/1/4 no_of_cursors,draw_flag C and xsiz,ysiz if single cursor C C last used cut values are written to C OUTPUTR/R/1/2 C C------------------------------------------------------------------- C IMPLICIT NONE C INTEGER IMNO,IAV,N,NMAL,NAXPIX(4),KEY INTEGER STAT1,STAT2,STAT,LFLAGS(2) INTEGER NPIXA(3),KCOOR(4),WPIX(3) INTEGER LOW1(3),HI(3),IFIRST INTEGER BOX(4),COOFF,UNI(1),NULO C CHARACTER ACTIO*(*) CHARACTER FRAMEA*100 CHARACTER INFO*44,AREA*80,FRMSTR*8,OPTION*5,CACT*8 C DOUBLE PRECISION STARTA(3),STEPA(3) C REAL RR,RVAL1,RVAL2,COOR(4),RCUTS(2) REAL PCUR(4),ZBINS(3),RBUF(8),ROLD(4) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:IDIDEV.INC' INCLUDE 'MID_INCLUDE:IDIMEM.INC' C DATA INFO + /' switch cursors on - next time we exit ... '/ C DATA FRAMEA /' '/ DATA KCOOR /-1,-1,-1,-1/ DATA LOW1 /1,1,1/, HI /1,1,1/ DATA NPIXA /1,1,1/, WPIX /1,1,1/ DATA COOFF /0/, IFIRST /0/ DATA FRMSTR /'E15.6'/ DATA ZBINS /256.0,0.0,0.0/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C attach Image Display + get loaded image CALL DTOPEN(1,STAT) CALL DTGICH(QDSPNO,QIMCH,FRAMEA,RBUF,STAT) C CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNO,STAT) CALL STDRDI(IMNO,'NAXIS',1,1,IAV,NAXPIX,UNI,NULO,STAT) CALL STDRDI(IMNO,'NPIX',1,NAXPIX(1),IAV,NPIXA,UNI,NULO,STAT) CALL STDRDD(IMNO,'START',1,NAXPIX(1),IAV,STARTA,UNI,NULO,STAT) CALL STDRDD(IMNO,'STEP',1,NAXPIX(1),IAV,STEPA,UNI,NULO,STAT) DO 100, N=1,3 NAXPIX(N+1) = NPIXA(N) 100 CONTINUE C LFLAGS(2) = D_R4_FORMAT !tell LOADWN that we opened as R4 data IF (ACTIO(1:1).EQ.'C') GOTO 490 C C --------------------------------------- C C here we handle the full image C C --------------------------------------- C OPTION = 'YYHNH' !full frame, histogram only BOX(1) = SSPX BOX(2) = SSPY LOW1(1) = SFPX LOW1(2) = SFPY WPIX(1) = NSX WPIX(2) = NSY LFLAGS(1) = -100 !whole frame in display window C 400 AREA(1:) = ' ' WRITE(AREA,10001) COOR(1),COOR(2),COOR(3),COOR(4) CALL BLANKO(AREA) !make statistics in displayed area CALL ZMSTAT(IMNO,AREA,NAXPIX,STARTA,STEPA,ZBINS,FRMSTR,OPTION) C 440 CALL PLOHI(IMNO) IF (IFIRST.EQ.0) THEN RVAL1 = 50 RVAL2 = 50 IFIRST = 1 ELSE CALL PTDATA(20,0,0,ROLD(1),ROLD(2),0.0,1) !show old cuts CALL PTDATA(20,0,0,ROLD(3),ROLD(4),0.0,1) ENDIF CALL PTGCUR(RVAL1,RVAL2,KEY,STAT) !graphCursor input IF (KEY.EQ.32) GOTO 8900 !gCursor exit C ROLD(1) = RVAL1 ROLD(2) = RVAL2 CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1) RCUTS(1) = RVAL1 CALL PTGCUR(RVAL1,RVAL2,KEY,STAT) IF (KEY.EQ.32) GOTO 8900 !gCursor exit C ROLD(3) = RVAL1 ROLD(4) = RVAL2 CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1) RCUTS(2) = RVAL1 IF (RCUTS(1).GT.RCUTS(2)) THEN RR = RCUTS(1) RCUTS(1) = RCUTS(2) RCUTS(2) = RR ENDIF WRITE(AREA,10010) RCUTS(1),RCUTS(2) CALL STTDIS(AREA,99,STAT) CALL STDWRR(IMNO,'LHCUTS',RCUTS,1,2,UNI,STAT) C C load image with new cuts C CALL LOADWN(LFLAGS,IMNO,NPIXA,LOW1,WPIX,BOX,RCUTS) CALL DAZVIS(QDSPNO,QIMCH,2,1) GOTO 440 C C --------------------------------------- C C here we handle cursor defined sub-frames C C --------------------------------------- C 490 LFLAGS(1) = -99 !subframes in display window CACT = 'YNYY?C0 ' CALL SETCUR(QDSPNO,2,0,2,KCOOR,STAT) !set cursor shape FRAMEA(1:) = ' ' !to init GETCUR ... NMAL = -1 OPTION = 'NYHNH' !subframe, histogram only C C input cursor positions and check status 500 NMAL = NMAL + 1 550 CALL GETCUR(CACT,FRAMEA, + BOX(1),PCUR(1),COOR(1),RVAL1,STAT1, + BOX(3),PCUR(3),COOR(3),RVAL2,STAT2) C C check cursor status IF ((STAT1.EQ.0).AND.(STAT2.EQ.0)) THEN IF ((NMAL.EQ.0) .AND. (COOFF.EQ.0)) THEN COOFF = 1 FRAMEA(1:) = ' ' CALL STTPUT(INFO,STAT) GOTO 550 ELSE GOTO 9000 !get out ENDIF ENDIF C AREA(1:) = ' ' WRITE(AREA,10001) COOR(1),COOR(2),COOR(3),COOR(4) CALL BLANKO(AREA) CALL ZMSTAT(IMNO,AREA,NAXPIX,STARTA,STEPA,ZBINS,FRMSTR,OPTION) C 700 CALL PLOHI(IMNO) IF (IFIRST.EQ.0) THEN RVAL1 = 50 RVAL2 = 50 IFIRST = 1 ELSE CALL PTDATA(20,0,0,ROLD(1),ROLD(2),0.0,1) !show old cuts CALL PTDATA(20,0,0,ROLD(3),ROLD(4),0.0,1) ENDIF CALL PTGCUR(RVAL1,RVAL2,KEY,STAT) IF (KEY.EQ.32) GOTO 1000 !move to next subframe C ROLD(1) = RVAL1 ROLD(2) = RVAL2 CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1) RCUTS(1) = RVAL1 CALL PTGCUR(RVAL1,RVAL2,KEY,STAT) IF (KEY.EQ.32) GOTO 1000 !move to next subframe C ROLD(3) = RVAL1 ROLD(4) = RVAL2 CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1) RCUTS(2) = RVAL1 IF (RCUTS(1).GT.RCUTS(2)) THEN RR = RCUTS(1) RCUTS(1) = RCUTS(2) RCUTS(2) = RR ENDIF WRITE(AREA,10010) RCUTS(1),RCUTS(2) CALL STTDIS(AREA,99,STAT) C C load sub-image with new cuts C LOW1(1) = NINT(PCUR(1)) LOW1(2) = NINT(PCUR(2)) HI(1) = NINT(PCUR(3)) HI(2) = NINT(PCUR(4)) WPIX(1) = HI(1) - LOW1(1) + 1 WPIX(2) = HI(2) - LOW1(2) + 1 CALL LOADWN(LFLAGS,IMNO,NPIXA,LOW1,WPIX,BOX,RCUTS) GOTO 700 C 1000 CALL STTPUT(AREA,STAT) CALL STTDIS + ('Move to next subframe ',99,STAT) IFIRST = 0 GOTO 500 C C That's it folks... 8900 CALL STTPUT(AREA,STAT) !write permanently CALL STDWRR(IMNO,'LHCUTS',RCUTS,1,2,UNI,STAT) !update descr. LHCUTS C 9000 CALL STTDIS(' ',-9,STAT) CALL DTCLOS(QDSPNO) CALL STKWRR('OUTPUTR',RCUTS,1,2,UNI,STAT) !save last cuts RETURN C 10001 FORMAT('[',G12.6,',',G12.6,':',G12.6,',',G12.6,'] ') 10010 FORMAT('Cut values used: ',2F12.5) END SUBROUTINE ERASE(CIMG,NPIX,NLIN,ICUR,IMAP, + K,L,OFS,IT,CRMS,NOISE,ISTAT) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine ERASE C (H.Waldthausen) C K. Banse version 1.00 880922 C R. Hook version 1.10 890621 C - Less restrictive area of fit C - Optional noise according to new logical flag (NOISE) C C.TYPE subroutine C C.KEYWORDS interpolation C C.PURPOSE C C This routine removes an unwanted detail in an image. C The surface is C an approximation to the weighted,least- squares Chebyshev C series surface fit.The weights are derived from the Poisson C error for the two first modes. C C.ALGORITHM 1) read cursor coordinates C 2) construct surface area around cursor C 3) replace bad pixels inside cursor C 4) optionally add noise C C.MODULES CALLS E02CAF,E02CBF from NAGLIB C C.INPUT/OUPUT parameters C C ->CIMG input image to be updated C ->NPIX x-pixels and C ->NLIN y-lines of CIMG C ->IMAP surface area indices C ->ICUR cursor area indices C ->K surface fit order in x C ->L surface fit order in y C ->OFS add this pixel offset to CIMG C ->IT number of iterations C ->CRMS rms of CIMG (if pos. take this value) C ->NOISE logical flag for whether patch should be filled with noise C <-ISTAT output status C : 0 O.K. C : POS 1-5 failure of E02CAF or E02CBF routines C : 6 fit order exeeded C : 7 surface-size larger than 102400 pixels C : 8 subimage-row undersampled C : 9 subimage-bounds are wrong C C CHANGED: C 28.4.98 swolf@eso.org - blow area up: 60x60 --> 320x320 C 24.9.98 swolf@eso.org - check areas (xlo