C @(#)smooth.for 13.1.1.1 (ESO-IPG) 06/02/98 18:11:36 PROGRAM SMOOTH C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program SMOOTH version 3.30 850729 C K. Banse ESO - Garching C C.KEYWORDS C smoothing, average, median, template matching C C.PURPOSE C (1) filter an image with a Gaussian filter C (2) filter an image with a digital filter (given values) C (3) filter an image with a digital filter (kernel = another frame) C (4) smooth an image by a) averaging in the neighborhood C b) taking the median in the neighbourhood C (5) expand an image C (6) filter an image with a maximum filter C (7) filter an image with a minimum filter FILACT set to these values... C C.ALGORITHM C (1) get the radius in x and y + threshold for the averaging/median algorithm C from INPUTR and calculate new points via: C a) B(n) = SUM A(j), j in neighborhood of given radius around n, if C abs(SUM - A(n)) > abs(threshold*A(n)), else B(n) = A(n) C or if C abs(SUM - A(n)) > threshold , else B(n) = A(n) C b) B(n) = median of A(j), j in neighborhood of given radius around n, if C abs(median - A(n)) > abs(threshold*A(n)), else B(n) = A(n) C or if C abs(median - A(n)) > threshold , else B(n) = A(n) C C (2) get the weights for the template from INPUTR or IN_B and calculate new C points via: C B(n) = SUM A(j)*W(l), j in neighborhood around n, l according to template C numbering: 1 ................ nw (nw = 2*radius + 1) C nw+1 ........... 2*nw C ..................... C (nw-1)*nw+1 ... nw*nw 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 INPUTC/C/1/2 (1:1) = Y, include center pixel in calculation C = N, do not include center pixel C (2:2) = A, use absolute threshold C = R, use relative threshold with pixel value C = Q, use rel. threshold with average/median value C = F, use fixed threshold-interval C = Z, use a+b*SQRT(median) as threshold C ACTION/C/1/5 (1:1) = A for averaging, M for median, C D for digital, G for Gaussian C K for maximum filter C L for minimum filter C E for just expanding the frame C (2:2) = L for looping, N for not looping C (3:3) = D for immediate display, N for no display C (4:4) = E for expansion of input frame C (5:5) = W for auxiliary window, N for not C (5:5) not implemented yet! C INPUTR/R/1/4 radius ( >0 )in x + y, C threshold(s) for averaging algorithm C INPUTC/C/3/60 area for filtering, C if set to +, complete frame is filtered C [...:...] or CURSOR or table_name C INPUTC/C/63/2 :x or :y to indicate that 1-dim digital filter C IN_B/C/1/60 name of filter image or template C INPUTR/R/1/40 for (2) weights of template C for (3) xmean,xsigma,ymean,ysigma C INPUTI/I/1/2 x-, y-radius for Gaussian filter image (must be odd) C INPUTI/I/1/4 xlef,xri,ylef,yri expansion pixels for EXPANSION C C C.VERSIONS C 4.80 920430: add support for 3-dim frames C 4.90 921204: add flag for no-expansion of input frame C 5.00 940317: include "Z" method for median C C-------------------------------------------------- C IMPLICIT NONE C INTEGER AREAFL,IDUM(1),IAV,IOMODE INTEGER MINSZX,MINSZY,FF,ZFLAG INTEGER HIOFF,HOOFF,HISIZ,HOSIZ,HY,LINCO,REMAIN INTEGER*8 PNTRHI,PNTRHO,VPNTR INTEGER HACK,RETBUF(3) INTEGER N,NAXIS,NCOL,NMAL,NOPNTS,NROW,NWORK INTEGER CHUNKI,CHUNKO,K1,K2 INTEGER*8 PNTRA,PNTRC,PNTRI,PNTRO,PNTRX,PNTRY,FPNTR INTEGER QPIX(2) INTEGER STAT,STAT1,STAT2,SUBDIM,COOFF INTEGER SUBPIX(2,2),SUBOFA(2),SUBOFB(2),RADIUS(2) INTEGER SUBLO(3),SUBHI(3) INTEGER ICUR1(2),ICUR2(2),LOADFL(2),CUTIDX,LOW(3) INTEGER NPIX(3),TPIX(3),RESPIX(3),OPIX(2),KCOOR(4) INTEGER IMNOA,IMNOC,IMNOI,IMNOO,IMNOX,FIMNO,IMNOEX INTEGER TID,TBCLNM(4),ALL,TABNUL(4) INTEGER FPIX(2),FILACT,KLOOP,DAZO,TRUNAX,HCOUNT INTEGER FNAXIS,FSIZE,WW(4),OFFPIX,RESPLZ,INPLZ,COUNT3 INTEGER EXNPIX(3),EXBUF(4),OLDMAP,DEBUFL,EXPAZ,OFFEX INTEGER UNIT(1),NULLO,MADRID(1) C CHARACTER*60 FRAMEA,FRAMEC,FRAME,TABLE,ERMES3 CHARACTER CUNIT*64,IDENT*72,CBUF*80,ACTION*5,DEFAUL*2 CHARACTER TLABL(4)*16,CAREA*60,XYCH(2)*1,MIDUNI*2 CHARACTER TMPFIL*14,CACT*8 C DOUBLE PRECISION STEP(3),START(3),EXSTA(3) DOUBLE PRECISION FSTART(2),FSTEP(2) C REAL INPUTR(40),RMIN,RMAX REAL PCUR1(6),PCUR2(6) REAL RDUM(1),KUTS(4),CUTS(4),RBUF(17) REAL DGFIL1(9),DGFIL2(9),DGFIL3(9),DGFIL4(9) REAL DGFIL5(9),DGFIL6(9),DGFIL7(25),DGFIL8(25) C LOGICAL SELFLG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C INCLUDE 'MID_INCLUDE:IDIDEV.INC' INCLUDE 'MID_INCLUDE:IDIMEM.INC' C COMMON /VMR/ MADRID C DATA SUBLO /3*1/, SUBHI /3*1/ DATA TPIX /3*1/, LOW /3*1/ DATA NPIX /3*1/, KCOOR /-1,-1,-1,-1/ DATA RESPIX /3*1/, KUTS /4*0.0/ DATA IDENT /' '/, CUNIT /' '/ DATA KLOOP /0/, LOADFL /0,0/ DATA DAZO /0/, FIMNO /-1/ DATA ACTION /' '/, COOFF /0/ DATA EXNPIX /1,1,1/, EXSTA /0.,0.,0./ DATA IMNOI /-1/, IMNOO /-1/ DATA OLDMAP /0/, EXPAZ /1/ DATA START /3*0.D0/, STEP /3*1.D0/ C DATA ERMES3 /' both cursors are off, set one on ... '/ DATA TLABL /'XSTART ','YSTART ','XEND ','YEND '/ DATA XYCH /'x','y'/ DATA CACT /'NYYY?C0 '/ C C data for filters CONE, LAPLACE, LOW_PASS, POINT, CONE5: DATA DGFIL1 /1.,-6.,1.,-6.,36.,-6.,1.,-6.,1./ !CONE DATA DGFIL2 /0.,-1.,0.,-1.,5.,-1.,0.,-1.,0./ !LAPLACE DATA DGFIL3 /1.,2.,1.,2.,5.,2.,1.,2.,1./ !LOW_PASS DATA DGFIL4 /-1.,-2.,-1.,-2.,13.,-2.,-1.,-2.,-1./ !POINT DATA DGFIL5 /1.,1.,1.,1.,8.,1.,1.,1.,1./ !TENT DATA DGFIL6 /-1.,-2.,-1.,-2.,13.,-2.,-1.,-2.,-1./ !SHARP DATA DGFIL7 /0.,0.,1.,0.,0.,0.,2.,2.,2.,0., !CONE5 + 1.,2.,5.,2.,1.,0.,2.,2.,2.,0.,0.,0.,1.,0.,0./ DATA DGFIL8 /1.,2.,3.,2.,1.,2.,4.,6.,4.,2., !PYRAMID5 + 3.,6.,9.,6.,3.,2.,4.,6.,4.,2.,1.,2.,3.,2.,1./ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up MIDAS environment + get current MIDAS unit CALL STSPRO('SMOOTH') CALL STKRDC('MID$SESS',1,11,2,IAV,MIDUNI,UNIT,NULLO,STAT) TMPFIL(1:) = 'filt'//MIDUNI//'dum.bdf ' C C get method + special options and branch according to filter type CALL STKRDC('ACTION',1,1,5,IAV,ACTION,UNIT,NULLO,STAT) CALL UPCAS(ACTION,ACTION) IF (ACTION(2:2).EQ.'L') KLOOP = 1 !update loop flag IF (ACTION(3:3).EQ.'D') LOADFL(1) = -99 !set display flag IF ((ACTION(4:4).EQ.'N') .AND. !clear expansion flag + (LOADFL(1).NE.-99)) EXPAZ = 0 !only if no display... CALL STKRDC('MID$SPEC',1,1,5,IAV,CAREA,UNIT,NULLO,STAT) IF (CAREA(1:5).EQ.'DEBUG') THEN DEBUFL = 1 ELSE DEBUFL = 0 ENDIF C FF = 0 PNTRHI = -1 !to know if it's first time IF (ACTION(1:1).EQ.'G') THEN !Gaussian CALL STKRDI('INPUTI',1,2,IAV,RADIUS,UNIT,NULLO,STAT) CALL STKRDR('INPUTR',1,4,IAV,INPUTR,UNIT,NULLO,STAT) CALL STKRDC('IN_B',1,1,60,IAV,FRAME,UNIT,NULLO,STAT) ZFLAG = 1 !flag for routine TEMPLB ELSE IF (ACTION(1:1).EQ.'D') THEN !digital CALL STKRDC('IN_B',1,1,60,IAV,FRAME,UNIT,NULLO,STAT) ZFLAG = 2 ELSE IF ((ACTION(1:1).EQ.'K') .OR. (ACTION(1:1).EQ.'L')) THEN CALL STKRDI('INPUTI',1,2,IAV,RADIUS,UNIT,NULLO,STAT) FILACT = 6 IF (ACTION(1:1).EQ.'K') THEN ZFLAG = 3 !maximum filter ELSE ZFLAG = 4 !minimum filter ENDIF ELSE CALL STKRDC('INPUTC',1,1,2,IAV,DEFAUL,UNIT,NULLO,STAT) CALL UPCAS(DEFAUL,DEFAUL) IF ((DEFAUL(2:2).EQ.'F') .OR. (DEFAUL(2:2).EQ.'Z')) THEN CALL STKRDR('INPUTR',1,2,IAV,INPUTR,UNIT,NULLO,STAT) CALL STKRDC('P3',1,1,60,IAV,CBUF,UNIT,NULLO,STAT) N = INDEX(CBUF,',') CBUF(N:N) = ' ' N = INDEX(CBUF,',') + 1 !skip over 2. comma CBUF(1:) = CBUF(N:) !shift to front C IF (DEFAUL(2:2).EQ.'Z') THEN CALL GENCNV(CBUF,2,2,IDUM,INPUTR(3),STEP,IAV) IF (IAV.LT.2) CALL STETER + (18,'Invalid syntax for threshold factors') ELSE IF (CBUF(1:1).EQ.'<') THEN FF = 1 CALL GENCNV(CBUF(3:),2,1,IDUM,RDUM,STEP,IAV) INPUTR(4) = RDUM(1) ELSE CALL GENCNV(CBUF,2,1,IDUM,RDUM,STEP,IAV) IF (IAV.LT.1) CALL STETER + (18,'Invalid syntax for interval limits') INPUTR(3) = RDUM(1) N = INDEX(CBUF,',') + 1 IF (CBUF(N:N) .EQ. '>') THEN FF = 2 ELSE CALL GENCNV(CBUF(N:),2,1,IDUM,RDUM,STEP,IAV) INPUTR(4) = RDUM(1) ENDIF ENDIF IF (IAV.LT.1) + CALL STETER(18,'Invalid syntax for interval limits') ENDIF ELSE CALL STKRDR('INPUTR',1,3,IAV,INPUTR,UNIT,NULLO,STAT) ENDIF ENDIF C C C get area option, input frame + output frame CAREA(1:) = ' ' CALL STKRDC('INPUTC',1,3,60,IAV,CAREA,UNIT,NULLO,STAT) CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,UNIT,NULLO,STAT) CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,UNIT,NULLO,STAT) CALL GENEQF(FRAMEA,FRAMEC,IAV) !check that in, out_frame are different IF (IAV.NE.0) + CALL STETER + (2,'input and result frame have to be different...') C CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT) CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXIS,UNIT,NULLO,STAT) CALL STDRDI(IMNOA,'NPIX',1,NAXIS,IAV,NPIX,UNIT,NULLO,STAT) CALL STDRDD(IMNOA,'START',1,NAXIS,IAV,START,UNIT,NULLO,STAT) CALL STDRDD(IMNOA,'STEP',1,NAXIS,IAV,STEP,UNIT,NULLO,STAT) CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENT,UNIT,NULLO,STAT) IF (NAXIS.GT.3) CALL STETER(1,'max. 3-dim frames supported...') N = (NAXIS+1) * 16 CALL STDRDC(IMNOA,'CUNIT',1,1,N,IAV,CUNIT,UNIT,NULLO,STAT) C CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNIT,NULLO,STAT) IF (CUTS(1).GE.CUTS(2)) THEN CUTIDX = 3 ELSE CUTIDX = 1 ENDIF C C check for < or > in fixed interval option (only average+median) IF (FF.EQ.1) THEN INPUTR(3) = CUTS(3) - 0.1 ELSE IF (FF.EQ.2) THEN INPUTR(4) = CUTS(4) + 0.1 ENDIF C C determine "real" NAXIS TRUNAX = NAXIS IF (NAXIS.EQ.3) THEN IF (NPIX(3).EQ.1) THEN IF (NPIX(2).EQ.1) THEN TRUNAX = 1 ELSE TRUNAX = 2 ENDIF ENDIF ELSE IF (NAXIS.EQ.2) THEN IF (NPIX(2).EQ.1) TRUNAX = 1 ENDIF C C only for complete frames do we support 3-dim frames IF (TRUNAX.EQ.3) THEN IF (CAREA(1:1).NE.'+') THEN CBUF(1:) = + 'max 2-dim frames supported for chosen area option... ' CALL STETER(23,CBUF) ENDIF C CBUF(1:) = '3-dim input frame - all planes are processed ... ' CALL STTPUT(CBUF,STAT) ENDIF C COUNT3 = 1 !counter for planes in 3-dim case OFFPIX = 0 LOADFL(2) = D_R4_FORMAT !use real data format C C ------------------------------------------------------------------- C C construct Gaussian filter C C ------------------------------------------------------------------- C 500 IF (ACTION(1:1).EQ.'G') THEN IF (NAXIS.EQ.1) RADIUS(2) = 0 !force 1-dim filtering... IF ((RADIUS(1).LE.0).OR.(RADIUS(2).LT.0)) + CALL STETER(4,'invalid radius chosen...') IF (RADIUS(2).GT.0) THEN FNAXIS = 2 ELSE FNAXIS = 1 ENDIF C FPIX(1) = 2*RADIUS(1) + 1 !get no. of points FPIX(2) = 2*RADIUS(2) + 1 !and no. of lines... FSIZE = FPIX(1) * FPIX(2) FSTART(1) = 0.D0 FSTART(2) = 0.D0 FSTEP(1) = 1.D0 FSTEP(2) = 1.D0 C C check, if we should save the filter frame IF (FRAME(1:1).EQ.'+') THEN !No, don't keep... FRAME = 'dummy_frame ' IOMODE = F_X_MODE ELSE !Yes, keep filter frame... IOMODE = F_O_MODE ENDIF IF (FIMNO.GT.-1) CALL STFCLO(FIMNO,STAT) !because of loop... CALL STIPUT(FRAME,D_R4_FORMAT,IOMODE,F_IMA_TYPE, + FNAXIS,FPIX,FSTART,FSTEP, + 'Gaussian filter','none ',FPNTR,FIMNO,STAT) C C here we really build the filter (and normalize it) CALL GAUSS(MADRID(FPNTR),FNAXIS,FPIX,FSTART,FSTEP,INPUTR, + 4,RMIN,RMAX) CALL NORMIT(MADRID(FPNTR),FSIZE) IF (IOMODE.EQ.F_O_MODE) THEN !if we keep file, get cuts RBUF(3) = CUTS(3) RBUF(4) = CUTS(4) CALL MYMX(MADRID(FPNTR),FSIZE,RBUF(3)) RBUF(1) = 0. RBUF(2) = 0. CALL STDWRR(FIMNO,'LHCUTS',RBUF,1,4,UNIT,STAT) ENDIF FILACT = 1 C C ------------------------------------------------------------------- C C digital filter C C ------------------------------------------------------------------- C ELSE IF (ACTION(1:1).EQ.'D') THEN CALL UPCAS(FRAME,CBUF) RADIUS(1) = 1 !default to 3*3 matrix with weights RADIUS(2) = 1 FILACT = 9 !and fixed filter name K1 = 9 C C predefined filters RDUM(1) = 0. !we want to normalize filter kernels IF (CBUF(1:4).EQ.'CONE') THEN IF (CBUF(5:5).EQ.'5') THEN DO 800, N=1,25 INPUTR(N) = DGFIL7(N) RDUM(1) = RDUM(1) + INPUTR(N) 800 CONTINUE RADIUS(1) = 2 RADIUS(2) = 2 K1 = 25 ELSE DO 805, N=1,9 INPUTR(N) = DGFIL1(N) RDUM(1) = RDUM(1) + INPUTR(N) 805 CONTINUE ENDIF ELSE IF (CBUF(1:7).EQ.'LAPLACE') THEN DO 810, N=1,9 INPUTR(N) = DGFIL2(N) RDUM(1) = RDUM(1) + INPUTR(N) 810 CONTINUE ELSE IF (CBUF(1:8).EQ.'LOW_PASS') THEN DO 820, N=1,9 INPUTR(N) = DGFIL3(N) RDUM(1) = RDUM(1) + INPUTR(N) 820 CONTINUE ELSE IF (CBUF(1:5).EQ.'POINT') THEN DO 830, N=1,9 INPUTR(N) = DGFIL4(N) RDUM(1) = RDUM(1) + INPUTR(N) 830 CONTINUE ELSE IF (CBUF(1:4).EQ.'TENT') THEN DO 840, N=1,9 INPUTR(N) = DGFIL5(N) RDUM(1) = RDUM(1) + INPUTR(N) 840 CONTINUE ELSE IF (CBUF(1:5).EQ.'SHARP') THEN DO 850, N=1,9 INPUTR(N) = DGFIL6(N) RDUM(1) = RDUM(1) + INPUTR(N) 850 CONTINUE ELSE IF (CBUF(1:5).EQ.'PYRAM') THEN DO 860, N=1,25 INPUTR(N) = DGFIL8(N) RDUM(1) = RDUM(1) + INPUTR(N) 860 CONTINUE RADIUS(1) = 2 RADIUS(2) = 2 K1 = 25 C C any weights file ELSE IF (CBUF(1:2).NE.'++') THEN IF (FIMNO.GT.-1) CALL STFCLO(FIMNO,STAT) CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,FNAXIS,FPIX,FSTART,FSTEP, + CBUF,CBUF,FPNTR,FIMNO,STAT) C C make sure, we have odd NPIX(i) IF (FPIX(1).EQ.1) THEN !kernel is a column RADIUS(1) = 0 ELSE IF (MOD(FPIX(1),2).EQ.0) THEN FPIX(1) = FPIX(1) - 1 WRITE(CBUF,30001) XYCH(1),FPIX(1) CALL STTPUT(CBUF,STAT) ENDIF RADIUS(1) = FPIX(1)/2 ENDIF IF (FNAXIS.EQ.1) THEN !kernel is a row RADIUS(2) = 0 ELSE IF (MOD(FPIX(2),2).EQ.0) THEN FPIX(2) = FPIX(2) - 1 WRITE(CBUF,30001) XYCH(2),FPIX(2) CALL STTPUT(CBUF,STAT) ENDIF RADIUS(2) = FPIX(2)/2 ENDIF FILACT = 3 C C the weights are entered directly ELSE CALL STKRDC('P3',1,1,80,IAV,CBUF,UNIT,NULLO,STAT) K1 = 1 DO 880, N=1,80 !get no. of values entered IF (CBUF(N:N).EQ.',') K1 = K1 + 1 IF (CBUF(N:N).EQ.' ') GOTO 888 880 CONTINUE 888 IF (K1.GT.39) K1 = 39 IF (MOD(K1,2).EQ.0) THEN K1 = K1 - 1 WRITE(CBUF,30006) K1 CALL STTPUT(CBUF,STAT) ENDIF C CALL STKRDC('INPUTC',1,63,2,IAV,CBUF,UNIT,NULLO,STAT) IF (CBUF(1:1).EQ.':') THEN IF ((CBUF(2:2).EQ.'Y').OR.(CBUF(2:2).EQ.'y')) THEN RADIUS(1) = 0 RADIUS(2) = K1/2 ELSE RADIUS(1) = K1/2 RADIUS(2) = 0 ENDIF ELSE !2-dim square kernel C IF (K1.EQ.25) THEN RADIUS(1) = 2 RADIUS(2) = 2 ELSE IF (K1.NE.9) THEN IF (K1.LT.23) THEN K1 = 9 !3x3 or 5x5 are only possibilities ELSE K1 = 25 RADIUS(1) = 2 RADIUS(2) = 2 ENDIF WRITE(CBUF,30007) K1 CALL STTPUT(CBUF,STAT) ENDIF ENDIF C CALL STKRDR('INPUTR',1,K1,IAV,INPUTR,UNIT,NULLO,STAT) IF (IAV .LT. K1) THEN N = K1 - IAV WRITE(CBUF,30005) N DO 890, N=IAV+1,K1 INPUTR(N) = 0.0 890 CONTINUE ENDIF FILACT = 2 ENDIF C C final check of radii IF ( (RADIUS(1).LT.0) .OR. (RADIUS(2).LT.0) .OR. + ((RADIUS(1).EQ.0).AND.(RADIUS(2).EQ.0)) ) + CALL STETER(6,'invalid no. of weights...') C C for predefined filters divide by kernel sum IF (FILACT.EQ.9) THEN RDUM(1) = 1./RDUM(1) DO 895, N=1,K1 INPUTR(N) = INPUTR(N) * RDUM(1) 895 CONTINUE FILACT = 2 !we have numbers ENDIF C C ------------------------------------------------------------------- C C just expand input frame => result frame C C ------------------------------------------------------------------- C ELSE IF (ACTION(1:1).EQ.'E') THEN CALL STKRDI('INPUTI',1,4,IAV,EXBUF,UNIT,NULLO,STAT) CALL EXPAX(IMNOA,OFFPIX,FRAMEC,EXBUF,STAT) IF (STAT.NE.0) CALL STETER(17,'invalid radius...') CALL STFOPN(FRAMEC,D_OLD_FORMAT,0,F_IMA_TYPE,IMNOC,STAT) CALL DSCUPT(IMNOA,IMNOC,' ',STAT) GOTO 9900 !leave straight away ... C C ------------------------------------------------------------------- C C average (smoothing) or median filter C C ------------------------------------------------------------------- C ELSE IF (FILACT.NE.6) THEN FILACT = 4 C C test, if center pixel should be included IF (DEFAUL(1:1).EQ.'Y') THEN ALL = 1 ELSE ALL = 0 ENDIF ACTION(2:2) = DEFAUL(2:2) C C get radii + threshold(s) for the averaging algorithm RADIUS(1) = NINT(INPUTR(1)) IF ( (NAXIS.GT.1) .AND. (NPIX(2).GT.1) ) THEN RADIUS(2) = NINT(INPUTR(2)) ELSE RADIUS(2) = 0 ENDIF IF ( (RADIUS(1).LT.0) .OR. (RADIUS(2).LT.0) .OR. + ((RADIUS(1)+RADIUS(2)).EQ.0) ) + CALL STETER(5,'invalid x- or y-radius entered...') C NOPNTS = RADIUS(2)*2 + 1 IF (ACTION(1:1).EQ.'M') THEN N = (RADIUS(1)*2 + 1) * NOPNTS NWORK = 2 * N !2 buffers of filter size ELSE IF (RADIUS(2).GT.0) THEN NWORK = NOPNTS + NPIX(1) ELSE NOPNTS = 0 NWORK = 1025 !for max RADIUS(1) of 512 ENDIF NWORK = NWORK * 2 !we want double precision ENDIF C IF (KLOOP.GT.1) CALL STFCLO(IMNOX,STAT) CALL STFCRE('smoothdumx',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + NWORK,IMNOX,STAT) CALL STFMAP(IMNOX,F_X_MODE,1,NWORK,IAV,PNTRX,STAT) IF (ACTION(1:1).EQ.'M') THEN VPNTR = PNTRX + N ELSE PNTRY = PNTRX + (NOPNTS*2) !for double precision ENDIF ENDIF C C ++++++++++++++++++++++++++ C C build up the result frame C C ++++++++++++++++++++++++++ C MINSZX = RADIUS(1)*2 + 1 !minimum size of area for given filtersize MINSZY = RADIUS(2)*2 + 1 IF (EXPAZ.EQ.0) THEN !adapt size of result frame RESPIX(1) = NPIX(1) - RADIUS(1) - RADIUS(1) RESPIX(2) = NPIX(2) - RADIUS(2) - RADIUS(2) EXSTA(1) = START(1) + (RADIUS(1)*STEP(1)) EXSTA(2) = START(2) + (RADIUS(2)*STEP(2)) ELSE !we do expand... RESPIX(1) = NPIX(1) RESPIX(2) = NPIX(2) EXSTA(1) = START(1) EXSTA(2) = START(2) ENDIF RESPIX(3) = NPIX(3) EXSTA(3) = START(3) C C if it's complete frame + NPIX(2) > MONITPAR(20) - hack it up C but we do it only once in the beginning... C IF (KLOOP.LE.1) THEN NWORK = RESPIX(3) * RESPIX(2) * RESPIX(1) CALL STFCRE(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NWORK,IMNOC,STAT) CALL STDWRI(IMNOC,'NAXIS',NAXIS,1,1,UNIT,STAT) CALL STDWRI(IMNOC,'NPIX',RESPIX,1,NAXIS,UNIT,STAT) CALL STDWRD(IMNOC,'START',EXSTA,1,NAXIS,UNIT,STAT) CALL STDWRD(IMNOC,'STEP',STEP,1,NAXIS,UNIT,STAT) CALL STDWRC(IMNOC,'IDENT',1,IDENT,1,72,UNIT,STAT) N = (NAXIS+1) * 16 CALL STDWRC(IMNOC,'CUNIT',1,CUNIT,1,N,UNIT,STAT) CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNIT,STAT) C HACK = 0 !get max. lines for buffers CALL HACKUP(RESPIX,D_R4_FORMAT,RETBUF) NWORK = RETBUF(1) !no. of lines from input frame C IF (CAREA(1:1).EQ.'+') THEN !complete frame IF ( (NWORK.LT.RESPIX(2)) .OR. (TRUNAX.EQ.3) ) THEN IF (TRUNAX.EQ.3) THEN HACK = 2 !cube ... ELSE HACK = 1 N = RETBUF(2) IF (RETBUF(3).GT.0) N = N + 1 WRITE(CBUF,20010) N CALL STTPUT(CBUF,STAT) ENDIF CHUNKI = NWORK*RESPIX(1) !size of memory to allocate CALL STFXMP(CHUNKI,D_R4_FORMAT,PNTRHO,STAT) ELSE !map complete input frame CHUNKI = RESPIX(1) * RESPIX(2) CALL STFMAP(IMNOC,F_O_MODE,1,CHUNKI,IAV,PNTRC,STAT) ENDIF C C if not complete area, copy first ELSE CALL STTDIS + ('first we copy input frame to result frame...',0,STAT) C IF (EXPAZ.EQ.0) THEN !copy inner window to result frame LOW(1) = RADIUS(1) + 1 LOW(2) = RADIUS(2) + 1 CALL COPWI(IMNOA,NPIX,LOW,RESPIX, + IMNOC,RESPIX,SUBLO,NWORK) ELSE !straight file copy in chunks CHUNKI = NWORK * NPIX(1) CALL STFCRE('smoothdumi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + CHUNKI,IMNOI,STAT) CALL STFMAP(IMNOI,F_X_MODE,1,CHUNKI,IAV,PNTRI,STAT) CALL COPYFX(MADRID(PNTRI),IMNOA,IMNOC,NPIX,NWORK) CALL STFCLO(IMNOI,STAT) ENDIF IMNOI = -1 !so that we know later on ENDIF C IF (EXPAZ.EQ.0) THEN WRITE(CBUF,30003) RESPIX(1),RESPIX(2) CALL STTPUT(CBUF,STAT) INPLZ = NPIX(1) * NPIX(2) ELSE INPLZ = 0 ENDIF OFFEX = 0 RESPLZ = RESPIX(1) * RESPIX(2) ENDIF C C now expand frame to include also 2*RADIUS columns + lines EXBUF(1) = RADIUS(1) !expand input frame EXBUF(2) = RADIUS(1) EXBUF(3) = RADIUS(2) EXBUF(4) = RADIUS(2) C C ******************************************************* C ******************************************************* C C expand frame - this is a loop, if 3-dim frames involved C C ******************************************************* C ******************************************************* C 1200 HCOUNT = 1 FRAMEA = ' ' CALL GENTIM(FRAMEA) !get current time C IF (EXPAZ.EQ.1) THEN IF (TRUNAX.EQ.3)THEN WRITE(CBUF,20000) COUNT3,FRAMEA(1:40) ELSE WRITE(CBUF,20001) FRAMEA(1:40) ENDIF CALL STTPUT(CBUF,STAT) C CALL EXPAX(IMNOA,OFFPIX,TMPFIL,EXBUF,STAT) IF (STAT.NE.0) CALL STETER(17,'invalid radius...') C IF (CAREA(1:1).EQ.'+') THEN !complete frame IF (HACK.GE.1) THEN !if HACK >= 1, use STFOPN for input file CALL STFOPN(TMPFIL,D_R4_FORMAT,0,F_IMA_TYPE, + IMNOEX,STAT) CALL STDRDI(IMNOEX,'NPIX',1,2,IAV,EXNPIX,UNIT,NULLO,STAT) CALL STDRDD(IMNOEX,'START',1,2,IAV,EXSTA,UNIT,NULLO,STAT) IF (PNTRHI.EQ.-1) THEN NWORK = (RETBUF(1) + 2*RADIUS(2)) * EXNPIX(1) CALL STFXMP(NWORK,D_R4_FORMAT,PNTRHI,STAT) ENDIF ELSE CALL STIGET(TMPFIL,D_R4_FORMAT,F_I_MODE, + F_IMA_TYPE,2, !we keep NAXIS of FRAMEA + IAV,EXNPIX,EXSTA,STEP, + CBUF,CBUF,PNTRA,IMNOEX,STAT) ENDIF ELSE CALL STFOPN(TMPFIL,D_R4_FORMAT,0,F_IMA_TYPE,IMNOEX,STAT) CALL STDRDI(IMNOEX,'NPIX',1,2,IAV,EXNPIX,UNIT,NULLO,STAT) CALL STDRDD(IMNOEX,'START',1,2,IAV,EXSTA,UNIT,NULLO,STAT) ENDIF C ELSE IF (TRUNAX.EQ.3) THEN WRITE(CBUF,22000) COUNT3,FRAMEA(1:40) CALL STTDIS(CBUF,0,STAT) ENDIF IMNOEX = IMNOA EXNPIX(1) = NPIX(1) EXSTA(1) = START(1) EXNPIX(2) = NPIX(2) EXSTA(2) = START(2) C IF (CAREA(1:1).EQ.'+') THEN !complete frame IF (HACK.GE.1) THEN IF (PNTRHI.EQ.-1) THEN NWORK = (RETBUF(1) + 2*RADIUS(2)) * EXNPIX(1) CALL STFXMP(NWORK,D_R4_FORMAT,PNTRHI,STAT) ENDIF ELSE N = EXNPIX(1) * EXNPIX(2) !map input frame CALL STFMAP(IMNOEX,F_I_MODE,1,N,IAV,PNTRA,STAT) ENDIF ELSE CONTINUE !nothing to map now for subframes ENDIF ENDIF C C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C get the different area options for all filters C C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C if we are here from looping, we already have the area option IF (KLOOP.GT.1) THEN NMAL = 0 IF (AREAFL.EQ.2) THEN !it's the 2. time we're up here... COOFF = 0 GOTO 1500 !continue with cursor ELSE IF (AREAFL.EQ.3) THEN GOTO 2100 !continue with table ELSE NMAL = 1 IF (AREAFL.EQ.3) THEN GOTO 3000 !single subframe (fixed coords) ELSE OPIX(1) = EXNPIX(1) OPIX(2) = EXNPIX(2) GOTO 5000 !full image ENDIF ENDIF ENDIF C C a) complete input image - only here do we support 3-dim files IF (CAREA(1:1).EQ.'+') THEN AREAFL = 0 !use complete input frame... PNTRO = PNTRC !filter image is output image PNTRI = PNTRA !use full input image OPIX(1) = EXNPIX(1) !dimensions of expanded frame OPIX(2) = EXNPIX(2) GOTO 5000 ENDIF C C b) subimage defined by coordinates IF (CAREA(1:1).EQ.'[') THEN CALL EXTCOO(IMNOA,CAREA,2,SUBDIM,SUBLO,SUBHI,STAT) IF (STAT.NE.0) CALL STETER(4,' ') !stop, if invalid coordinates... SUBPIX(1,1) = SUBLO(1) SUBPIX(2,1) = SUBLO(2) SUBPIX(1,2) = SUBHI(1) SUBPIX(2,2) = SUBHI(2) AREAFL = 1 NMAL = 1 GOTO 3000 ENDIF C C c) cursor defined subimage(s) CALL UPCAS(CAREA(1:6),CUNIT) IF (CUNIT(1:6).NE.'CURSOR') GOTO 2000 !branch to table option C C get main control block for ImageDisplay + attach IF (DAZO.EQ.0) THEN CALL DTOPEN(1,STAT) DAZO = 1 ENDIF NMAL = 0 !reset counter CBUF(1:) = ' ' CALL CONCHA(QDSPNO,QOVCH,1,0) !clear overlay channel CALL SETCUR(QDSPNO,2,0,0,KCOOR,STAT) !set up cursor rectangle C C input cursor positions and check status 1500 CALL GETCUR(CACT,CBUF, + ICUR1,PCUR1(3),PCUR1(5),RDUM(1),STAT1, + ICUR2,PCUR2(3),PCUR2(5),RDUM(1),STAT2) C C check cursor status IF ( (STAT1.EQ.0) .AND. (STAT2.EQ.0) ) THEN IF (COOFF.EQ.0) THEN COOFF = 1 CALL STTDIS(ERMES3,0,STAT) CBUF(1:) = ' ' GOTO 1500 ELSE GOTO 8000 !cursor loop terminated ... ENDIF ELSE NMAL = NMAL + 1 !update counter COOFF = 1 ENDIF C C and get frame pixel no's SUBPIX(1,1) = PCUR1(3) SUBPIX(2,1) = PCUR1(4) SUBPIX(1,2) = PCUR2(3) SUBPIX(2,2) = PCUR2(4) AREAFL = 2 GOTO 3000 C C d) table defined subimage(s) - for this option we loop... 2000 CALL CLNTAB(CAREA,TABLE,0) CALL TBTOPN(TABLE,F_I_MODE,TID,STAT) CALL TBIGET(TID,NCOL,NROW,N,N,N,STAT) !get no. of rows + columns DO 2050, N=1,4 CALL TBLSER(TID,TLABL(N),TBCLNM(N),STAT) IF (TBCLNM(N).LE.0) THEN CBUF(1:) = 'column :'//TLABL(N)//'is missing ... ' CALL STETER(7,CBUF) ENDIF 2050 CONTINUE C AREAFL = 3 NMAL = 0 !init counter C C if we also want to display the filtered windows, we have to attach IF (LOADFL(1).NE.0) THEN IF (DAZO.EQ.0) THEN CALL DTOPEN(1,STAT) DAZO = 1 ENDIF CALL PIXXCV('INIT',IMNOA,RBUF,STAT) ENDIF C C here we loop 2100 IF (NMAL.GE.NROW) THEN WRITE(CBUF,30010) NMAL CALL STTPUT(CBUF,STAT) GOTO 8000 ENDIF C C get next row of values NMAL = NMAL + 1 CALL TBSGET(TID,NMAL,SELFLG,STAT) IF (.NOT.SELFLG) GOTO 2100 C CALL TBRGET(TID,'R*4 ',NMAL,4,TBCLNM,RBUF,TABNUL,STAT) PCUR1(3) = (RBUF(1)-START(1))/STEP(1) + 1. PCUR1(4) = (RBUF(2)-START(2))/STEP(2) + 1. PCUR2(3) = (RBUF(3)-START(1))/STEP(1) + 1. PCUR2(4) = (RBUF(4)-START(2))/STEP(2) + 1. IF (LOADFL(1).NE.0) THEN !we have to get related screen pixels RBUF(11) = PCUR1(3) RBUF(12) = PCUR1(4) CALL PIXXCV('_RS',0,RBUF(11),STAT) ICUR1(1) = NINT(RBUF(15)) ICUR1(2) = NINT(RBUF(16)) ENDIF C C get frame pixel no's SUBPIX(1,1) = PCUR1(3) SUBPIX(2,1) = PCUR1(4) SUBPIX(1,2) = PCUR2(3) SUBPIX(2,2) = PCUR2(4) C C ------------------------------------------------- C C common section for option b), c) and d) C C ------------------------------------------------- C C get no. of points involved 3000 TPIX(1) = SUBPIX(1,2) - SUBPIX(1,1) + 1 TPIX(2) = SUBPIX(2,2) - SUBPIX(2,1) + 1 C IF (TPIX(1).LT.MINSZX) THEN WRITE(CBUF,30004) NMAL,XYCH(1) CALL STTPUT(CBUF,STAT) GOTO 7000 ENDIF IF (TPIX(2).LT.MINSZY) THEN WRITE(CBUF,30004) NMAL,XYCH(2) CALL STTPUT(CBUF,STAT) GOTO 7000 ENDIF C C the memory allocation scheme wastes space for options c and d), but so what.. C OPIX(1) = TPIX(1) + 2*RADIUS(1) !size we need as input frame OPIX(2) = TPIX(2) + 2*RADIUS(2) N = OPIX(1) * OPIX(2) C IF (N.GT.OLDMAP) THEN OLDMAP = N !save size IF (IMNOI.GT.-1) THEN CALL STFCLO(IMNOI,STAT) CALL STFCLO(IMNOO,STAT) ENDIF CALL STFCRE('smoothdumi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + N,IMNOI,STAT) CALL STFMAP(IMNOI,F_X_MODE,1,N,IAV,PNTRI,STAT) N = TPIX(1) * TPIX(2) CALL STFCRE('smoothdumo',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + N,IMNOO,STAT) CALL STFMAP(IMNOO,F_X_MODE,1,N,IAV,PNTRO,STAT) ENDIF C IF (EXPAZ.EQ.0) THEN !get expanded subframe from original input SUBOFA(1) = SUBPIX(1,1) - RADIUS(1) SUBOFA(2) = SUBPIX(2,1) - RADIUS(2) LOW(1) = SUBOFA(1) + OPIX(1) - 1 LOW(2) = SUBOFA(2) + OPIX(2) - 1 IF ( (SUBOFA(1).LT.1) .OR. + (SUBOFA(2).LT.1) .OR. + (LOW(1).GT.NPIX(1)) .OR. + (LOW(2).GT.NPIX(2))) THEN WRITE(CBUF,30002) NMAL CALL STTPUT(CBUF,STAT) GOTO 7000 ENDIF ELSE !get expanded subframe from expanded input SUBOFA(1) = SUBPIX(1,1) !x-pixel to start (-radius+radius) SUBOFA(2) = SUBPIX(2,1) !y-pixel to start (-radius+radius) ENDIF SUBOFB(1) = 1 SUBOFB(2) = 1 CHUNKI = EXNPIX(1) * OPIX(2) K1 = (SUBOFA(2)-1)*EXNPIX(1) + 1 CALL STFMAP(IMNOEX,F_I_MODE,K1,CHUNKI,IAV,PNTRA,STAT) QPIX(1) = EXNPIX(1) QPIX(2) = OPIX(2) SUBOFA(2) = 1 CALL COPYF1(MADRID(PNTRA),QPIX,SUBOFA,OPIX,MADRID(PNTRI), + OPIX,SUBOFB) CALL STFUNM(IMNOEX,STAT) C 5000 IF (AREAFL.GE.2) THEN N = 1 CALL STKWRI('LOG',N,4,1,UNIT,STAT) !set LOG(4) = 1 ENDIF C C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C and finally do the filtering C C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ C IF (FILACT.EQ.4) THEN IF (HACK.EQ.0) THEN PNTRO = PNTRO + OFFPIX IF (RADIUS(2).EQ.0) THEN !average/median filter IF (ACTION(1:1).EQ.'M') THEN CALL AVER1M(ACTION,MADRID(PNTRI),MADRID(PNTRO), + OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),KUTS(3)) ELSE CALL AVER1A(ACTION,MADRID(PNTRI),MADRID(PNTRO), + OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),KUTS(3)) ENDIF ELSE IF (ACTION(1:1).EQ.'M') THEN CALL AVER2M(ACTION,MADRID(PNTRI),MADRID(PNTRO), + OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),MADRID(VPNTR),KUTS(3)) ELSE IF (ACTION(2:2).EQ.'Q') THEN CALL BVER2A(ACTION,MADRID(PNTRI),MADRID(PNTRO), + OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),MADRID(PNTRY),KUTS(3)) ELSE CALL AVER2A(ACTION,MADRID(PNTRI),MADRID(PNTRO), + OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),MADRID(PNTRY),KUTS(3)) ENDIF ENDIF ENDIF ELSE HIOFF = 1 + OFFEX HOOFF = 1 + OFFPIX LINCO = 0 REMAIN = RETBUF(1) !we know, that NPIX(2) >= RETBUF(1) C 6000 HY = REMAIN + 2*RADIUS(2) HISIZ = HY * EXNPIX(1) HOSIZ = REMAIN * RESPIX(1) CALL STFGET(IMNOEX,HIOFF,HISIZ,IAV,MADRID(PNTRHI),STAT) IF (HACK.EQ.1) THEN FRAMEA = ' ' CALL GENTIM(FRAMEA) !get current time WRITE(CBUF,20020) HCOUNT,FRAMEA(1:40) CALL STTDIS(CBUF,0,STAT) ENDIF C IF (RADIUS(2).EQ.0) THEN !average/median filter IF (ACTION(1:1).EQ.'M') THEN CALL AVER1M(ACTION,MADRID(PNTRHI),MADRID(PNTRHO), + OPIX(1),HY,RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),KUTS(3)) ELSE CALL AVER1A(ACTION,MADRID(PNTRHI),MADRID(PNTRHO), + OPIX(1),HY,RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),KUTS(3)) ENDIF ELSE IF (ACTION(1:1).EQ.'M') THEN CALL AVER2M(ACTION,MADRID(PNTRHI),MADRID(PNTRHO), + OPIX(1),HY,RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),MADRID(VPNTR),KUTS(3)) ELSE IF (ACTION(2:2).EQ.'Q') THEN CALL BVER2A(ACTION,MADRID(PNTRHI),MADRID(PNTRHO), + OPIX(1),HY,RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),MADRID(PNTRY),KUTS(3)) ELSE CALL AVER2A(ACTION,MADRID(PNTRHI),MADRID(PNTRHO), + OPIX(1),HY,RADIUS,INPUTR(3),ALL, + MADRID(PNTRX),MADRID(PNTRY),KUTS(3)) ENDIF ENDIF ENDIF C CALL STFPUT(IMNOC,HOOFF,HOSIZ,MADRID(PNTRHO),STAT) LINCO = LINCO + REMAIN !increase line count + check REMAIN = NPIX(2) - LINCO IF (REMAIN.GT.0) THEN IF (REMAIN.GT.RETBUF(1)) + REMAIN = RETBUF(1) !hack up into RETBUF(1)-chunks HIOFF = HIOFF + HISIZ - (2*RADIUS(2)*EXNPIX(1)) HOOFF = HOOFF + HOSIZ HCOUNT = HCOUNT + 1 GOTO 6000 !loop again ENDIF ENDIF C ELSE IF (HACK.EQ.0) THEN PNTRO = PNTRO + OFFPIX IF ((FILACT.EQ.1) .OR. (FILACT.EQ.3)) THEN CALL TEMPLB(ZFLAG,MADRID(PNTRI),MADRID(PNTRO),OPIX, + MADRID(FPNTR),RADIUS,KUTS(3),STAT) ELSE CALL TEMPLB(ZFLAG,MADRID(PNTRI),MADRID(PNTRO),OPIX, + INPUTR,RADIUS,KUTS(3),STAT) ENDIF IF (STAT.NE.0) + CALL STETER(5,'x- or y-radius of filter too large...') C ELSE HIOFF = 1 + OFFEX HOOFF = 1 + OFFPIX LINCO = 0 REMAIN = RETBUF(1) !we know, that NPIX(2) >= RETBUF(1) C 6600 HY = REMAIN + 2*RADIUS(2) HISIZ = HY * EXNPIX(1) HOSIZ = REMAIN * RESPIX(1) OPIX(2) = HY CALL STFGET(IMNOEX,HIOFF,HISIZ,IAV,MADRID(PNTRHI),STAT) IF (HACK.EQ.1) THEN FRAMEA = ' ' CALL GENTIM(FRAMEA) !get current time WRITE(CBUF,20020) HCOUNT,FRAMEA(1:40) CALL STTDIS(CBUF,0,STAT) ENDIF C IF ((FILACT.EQ.1) .OR. (FILACT.EQ.3)) THEN CALL TEMPLB(ZFLAG,MADRID(PNTRHI),MADRID(PNTRHO),OPIX, + MADRID(FPNTR),RADIUS,KUTS(3),STAT) ELSE CALL TEMPLB(ZFLAG,MADRID(PNTRHI),MADRID(PNTRHO),OPIX, + INPUTR,RADIUS,KUTS(3),STAT) ENDIF IF (STAT.NE.0) + CALL STETER(5,'x- or y-radius of filter too large...') C CALL STFPUT(IMNOC,HOOFF,HOSIZ,MADRID(PNTRHO),STAT) LINCO = LINCO + REMAIN !increase line count + check REMAIN = NPIX(2) - LINCO IF (REMAIN.GT.0) THEN IF (REMAIN.GT.RETBUF(1)) + REMAIN = RETBUF(1) !hack up into RETBUF(1)-chunks HIOFF = HIOFF + HISIZ - (2*RADIUS(2)*EXNPIX(1)) HOOFF = HOOFF + HOSIZ HCOUNT = HCOUNT + 1 GOTO 6600 !loop again ENDIF ENDIF C ENDIF C C if area option used, merge filtered frame back into output frame IF (AREAFL.GT.0) THEN SUBOFA(1) = 1 SUBOFA(2) = 1 IF (EXPAZ.EQ.0) THEN SUBOFB(1) = SUBPIX(1,1) - RADIUS(1) SUBOFB(2) = SUBPIX(2,1) - RADIUS(2) ELSE SUBOFB(1) = SUBPIX(1,1) SUBOFB(2) = SUBPIX(2,1) ENDIF CHUNKO = RESPIX(1) * TPIX(2) QPIX(1) = RESPIX(1) QPIX(2) = TPIX(2) K2 = (SUBOFB(2)-1)*RESPIX(1) + 1 CALL STFMAP(IMNOC,F_IO_MODE,K2,CHUNKO,IAV,PNTRC,STAT) SUBOFB(2) = 1 CALL COPYF1(MADRID(PNTRO),TPIX,SUBOFA,TPIX,MADRID(PNTRC), + QPIX,SUBOFB) CALL STFPUT(IMNOC,K2,CHUNKO,MADRID(PNTRC),STAT) IF (LOADFL(1).NE.0) THEN !for display option load subwindows LOW(1) = PCUR1(3) LOW(2) = PCUR1(4) WW(1) = ICUR1(1) WW(2) = ICUR1(2) WW(3) = SCALX WW(4) = SCALY CALL LOADWN(LOADFL,IMNOC,RESPIX,LOW,TPIX,WW,CUTS(CUTIDX)) ENDIF CALL STFUNM(IMNOC,STAT) ENDIF C C for cursor or table input loop 7000 IF (AREAFL.GE.2) THEN N = 0 CALL STKWRI('LOG',N,4,1,UNIT,STAT) !reset LOG(4) = 0 GOTO (1500,2100),AREAFL-1 C C also loop for 3-dim frames (with complete frame option) ELSE IF (TRUNAX.EQ.3) THEN COUNT3 = COUNT3 + 1 IF (COUNT3.LE.NPIX(3)) THEN OFFPIX = OFFPIX + RESPLZ OFFEX = OFFEX + INPLZ IF (EXPAZ.NE.0) CALL STFCLO(IMNOEX,STAT) GOTO 1200 ENDIF ENDIF C C update descriptors + add history 8000 IF (KLOOP.EQ.0) THEN CALL DSCUPT(IMNOA,IMNOC,' ',STAT) IF (AREAfl.EQ.0) + CALL STDWRR(IMNOC,'LHCUTS',KUTS,1,4,UNIT,STAT) GOTO 9000 ENDIF C C +++++++++++++++++++++++++++++++++ C C here we loop inside the module C C +++++++++++++++++++++++++++++++++ C IF (KLOOP.EQ.1) THEN KLOOP = 99 IF (ACTION(1:1).EQ.'D') THEN WRITE(*,10000) ELSE WRITE(*,10001) ENDIF ENDIF C 8800 CBUF(1:) = ' ' CALL STTDIS('> ',99,STAT) READ(*,30000) CBUF IF ((CBUF(1:1).EQ.'X').OR.(CBUF(1:1).EQ.'x')) THEN KLOOP = 0 GOTO 8000 ENDIF N = INDEX(CBUF,' ') C IF (ACTION(1:1).EQ.'G') THEN IF (CBUF(1:1).NE.'?') THEN CALL GENCNV(CBUF(1:N),1,2,RADIUS,RDUM,STEP,IAV) IF (IAV.LT.1) THEN WRITE(*,10011) GOTO 8800 ELSE IF (IAV.LT.2) THEN RADIUS(2) = RADIUS(1) ENDIF ENDIF N = N + 1 !move to next parameter INPUTR(5) = 1.0 IF ((CBUF(N:N).NE.'?') .AND. + (CBUF(1:1).NE.' ')) THEN CALL GENCNV(CBUF(1:N),2,4,IDUM,INPUTR,STEP,IAV) IF (IAV.LT.2) THEN INPUTR(5) = -1.0 ELSE IF (IAV.LT.4) THEN INPUTR(3) = INPUTR(1) INPUTR(4) = INPUTR(2) ENDIF ENDIF C ELSE IF (ACTION(1:1).EQ.'D') THEN CALL GENCNV(CBUF,2,9,IDUM,INPUTR,STEP,IAV) IF (IAV.LT.9) THEN FRAME(1:) = CBUF(1:) !it's a frame ELSE FRAME(1:2) = '++' ENDIF C ELSE IF ((CBUF(N:N).NE.'?') .AND. + (CBUF(1:1).NE.' ')) THEN CALL GENCNV(CBUF(1:N),2,3,IDUM,INPUTR,STEP,IAV) IF (IAV.LT.3) THEN WRITE(*,10010) GOTO 8800 ENDIF ENDIF N = N + 1 !move to next parameter IF (CBUF(N:N).NE.'?') + CALL UPCAS(CBUF(N:N+1),DEFAUL) ENDIF GOTO 500 C C free data 9000 IF (DAZO.NE.0) CALL DTCLOS(QDSPNO) CALL STFCLO(IMNOC,STAT) !close result frame CALL STFCLO(IMNOA,STAT) !close input frame CALL STFCLO(IMNOEX,STAT) IF (DEBUFL.EQ.0) THEN IF (EXPAZ.EQ.1) CALL STFDEL(TMPFIL,STAT) ENDIF 9900 CALL STSEPI C 10000 FORMAT(' Enter parameter P3 like in the command syntax,', + ' type X to exit') 10001 FORMAT(' Enter parameters P3 and P4 like in the command syntax,', + ' type X to exit') 10010 FORMAT(' Invalid input, we need 3 real numbers, try again...') 10011 FORMAT(' Invalid input, we need 2 integer numbers, try again...') 20000 FORMAT('expanding plane ',I3,' - ',A) 20001 FORMAT('expanding frame - ',A) 20010 FORMAT('we process frame in ',I2,' chunks ...') 20020 FORMAT('working on chunk ',I2,' - ',A) 22000 FORMAT('plane ',I3,' - ',A) 30000 FORMAT(A) 30001 FORMAT('Warning: No. of ',A,'-weights must be odd => changed to ', + I2,' ...') 30002 FORMAT('invalid subwindow no.',I4,' - so we skip...') 30003 FORMAT('no expansion - NPIX(1,2) of result frame = ',2I6) 30004 FORMAT + ('subimage no. ',I4,' too small in ',A,' for given filter size') 30005 FORMAT('missing kernel values (',I2,') are set to 0.0') 30006 FORMAT('Warning: No. of weights must be odd => changed to ', + I2,' ...') 30007 FORMAT + ('Warning: No. of weights must be 3x3 or 5x5 => changed to ', + I2,' ...') 30010 FORMAT(I4,' table entries processed ...') C END SUBROUTINE COPWI(IMNOA,NPIXA,BEGA,NPIXW,IMNOB,NPIXB,BEGB,STRIP) C IMPLICIT NONE C INTEGER IMNOA,NPIXA(2),BEGA(2),NPIXW(2) INTEGER IMNOB,NPIXB(2),BEGB(2),STRIP INTEGER*8 PNTRI,PNTRO INTEGER IMNOI,IMNOO INTEGER CHUNKI,CHUNKO,LOOPLM,RMAIND,STAT,IAV INTEGER IPIX(2),OPIX(2),WPIX(2),ILOW(2),OLOW(2) INTEGER K1,K2,N INTEGER MADRID(1) C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C LOOPLM = NPIXW(2) / STRIP RMAIND = NPIXW(2) - (LOOPLM*STRIP) C CHUNKI = STRIP * NPIXA(1) CALL STFCRE('copywork',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + CHUNKI,IMNOI,STAT) CALL STFMAP(IMNOI,F_X_MODE,1,CHUNKI,IAV,PNTRI,STAT) C CHUNKO = STRIP * NPIXB(1) CALL STFCRE('copzwork',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + CHUNKO,IMNOO,STAT) CALL STFMAP(IMNOO,F_X_MODE,1,CHUNKO,IAV,PNTRO,STAT) C IPIX(1) = NPIXA(1) IPIX(2) = STRIP OPIX(1) = NPIXB(1) OPIX(2) = STRIP ILOW(1) = BEGA(1) ILOW(2) = 1 OLOW(1) = BEGB(1) OLOW(2) = 1 WPIX(1) = NPIXW(1) WPIX(2) = STRIP C K1 = (BEGA(2)-1)*NPIXA(1) + 1 !full line offsets K2 = (BEGB(2)-1)*NPIXB(1) + 1 !full line offsets DO 1000, N=1,LOOPLM CALL STFGET(IMNOA,K1,CHUNKI,IAV,MADRID(PNTRI),STAT) !data in CALL COPYF1(MADRID(PNTRI),IPIX,ILOW,WPIX, + MADRID(PNTRO),OPIX,OLOW) CALL STFPUT(IMNOB,K2,CHUNKO,MADRID(PNTRO),STAT) !data out K1 = K1 + CHUNKI K2 = K2 + CHUNKO 1000 CONTINUE C IF (RMAIND.NE.0) THEN CHUNKI = RMAIND * NPIXA(1) CHUNKO = RMAIND * NPIXB(1) WPIX(2) = RMAIND CALL STFGET(IMNOA,K1,CHUNKI,IAV,MADRID(PNTRI),STAT) !data in CALL COPYF1(MADRID(PNTRI),IPIX,ILOW,WPIX, + MADRID(PNTRO),OPIX,OLOW) CALL STFPUT(IMNOB,K2,CHUNKO,MADRID(PNTRO),STAT) !data out ENDIF C CALL STFCLO(IMNOI,STAT) CALL STFCLO(IMNOO,STAT) C RETURN END