C @(#)necback.for 17.1.1.1 (ES0-DMD) 01/25/02 17:51:28 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 C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C. IDENTIFICATION C C Program ECHBACK version 1.00 : 27.06.90 C C Author : P. BALLESTER ESO - Garching C C. KEYWORDS C Echelle spectroscopy background rebinning C C. PURPOSE C Fit the background and generates a frame C C. ALGORITHM C Adapted from W. Verschueren, H. Hensberge, Astr. Astroph. 240,216(1990) C "Order extraction and background subtraction on CCD Caspec Echelle Spectra" C C. INPUT/OUTPUT C The following keywords are used : C IN_A/C/1/60 IN: input table (type ORDER) C IN_B/C/1/60 IN: reference image C INPUTI/I/1/1 IN: half-size of the measure window C INPUTI/I/2/1 IN: degree of the spline polynomials C INPUTR/R/1/1 IN: smoothing factor of the spline interpolation C OUT_A/C/1/60 OUT: output image (background estimate) C C. VERSIONS C 0.00 27.06.90 from version 1.20 of 19.02.90 of REBIN.FOR C 1.00 17.07.90 Reference version C 1.10 19.07.90 Speed improvements : from 9.26 Mmbbe to 2.06 Mmbbe C 22.06.99 Measurement mode (MMODE) added (swolf@eso.org) C. BUGS C Since the STDWRC standard interface does not work properly for felem=-1, C the first line of the HISTORY descriptor of the output frame is C overwritten (17.07.90) C C----------------------------------------------------------------------------- C IMPLICIT NONE C LOGICAL FLAG,NULL C LOGICAL DBG C INTEGER NAXIS,ROW,ABSROW,NPOINT,DEGREE INTEGER WIDTH,TYP,DUN,MMODE C INTEGER COL ! DBG INTEGER IMNOA,IMNOB INTEGER*8 PNTRA,PNTRB INTEGER NPIX(2),ERROR INTEGER IAV,UNI(1),NULO,STAT INTEGER MADRID(1) INTEGER TID,NCINT,NRINT,NACINT,NARINT,NSINT INTEGER POSX,POSY,POSBK INTEGER PSXSTA,PSXEND,PSYSTA,PSYEND REAL XSTA,XEND,YSTA,YEND C CHARACTER*60 ORDER,REFIMA,BAKIMA CHARACTER*80 STRING,HISTO CHARACTER CUNIT*48,IDENT*72 CHARACTER*16 LABX,LABY,LABBK,LBXSTA,LBXEND CHARACTER*16 LBYSTA,LBYEND,FORM,UNIT C PARAMETER (NPOINT=20000) REAL BACK(3,NPOINT),SMOOTH,MEASUR,CUTS(4) C DOUBLE PRECISION STEP(2),START(2) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA LABX,LABY,LABBK /'X' , 'YBKG' , 'BKG'/ DATA LBXSTA,LBXEND,LBYSTA,LBYEND /'XSTA','XEND','YSTA','YEND'/ C C --- 0. MIDAS environment set up and enables automatic error abort C CALL STSPRO ('ECHBAK') C C --- Set DEBUG flag or not C C DBG = (.FALSE.) C C --- 1. Read input/output descriptors C CALL STKRDC ('IN_A',1,1,60,IAV,ORDER,UNI,NULO,STAT) CALL STKRDC ('IN_B',1,1,60,IAV,REFIMA,UNI,NULO,STAT) CALL STKRDC ('OUT_A',1,1,60,IAV,BAKIMA,UNI,NULO,STAT) C C --- Translate names C CALL CLNTAB (ORDER,ORDER,0) CALL CLNFRA (REFIMA,REFIMA,0) CALL CLNFRA (BAKIMA,BAKIMA,0) C CALL STKRDI ('INPUTI',1,1,IAV,WIDTH,UNI,NULO,STAT) CALL STKRDR ('INPUTR',1,1,IAV,SMOOTH,UNI,NULO,STAT) CALL STKRDI ('INPUTI',2,1,IAV,DEGREE,UNI,NULO,STAT) CALL STKRDI ('INPUTI',3,1,IAV,MMODE,UNI,NULO,STAT) C C --- Display read descriptors C C IF (DBG) THEN C STRING = '*** Read descriptors : ' C CALL STTPUT (STRING,STAT) C STRING = ORDER(1:25)//REFIMA(1:25)//BAKIMA(1:25) C CALL STTPUT (STRING,STAT) C ENDIF C C --- 2. Get reference frame and map it. C CALL STIGET (REFIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIX,START,STEP,IDENT,CUNIT, + PNTRA,IMNOA,STAT) C C --- A call to STETER will abort the program... C IF (NAXIS.NE.2) + CALL STETER (1,'Input frame must have two dimensions...') C C --- 3.1. Get input table and open it C CALL TBTOPN (ORDER,F_IO_MODE,TID,STAT) C C --- 3.2. Read table descriptors C CALL TBIGET (TID,NCINT,NRINT,NSINT,NACINT,NARINT,STAT) IF (NRINT.EQ.0) THEN STRING = '*** FATAL : There are no data in the table' CALL STETER (9,STRING) ENDIF C C --- 4. Table columns processing C C --- 4.1. Search for columns defining background positions C CALL TBLSER (TID,LABX ,POSX, STAT) IF (POSX.EQ.-1) THEN STRING = 'Column label '//LABX//'not found in table '//ORDER CALL STTPUT (STRING,STAT) ENDIF C CALL TBLSER (TID,LABY ,POSY, STAT) IF (POSY.EQ.-1) THEN STRING = 'Column label '//LABY//'not found in table '//ORDER CALL STTPUT (STRING,STAT) ENDIF CALL TBLSER (TID,LBXSTA ,PSXSTA, STAT) CALL TBLSER (TID,LBXEND ,PSXEND, STAT) CALL TBLSER (TID,LBYSTA ,PSYSTA, STAT) CALL TBLSER (TID,LBYEND ,PSYEND, STAT) C C --- 4.2. Create column for measured background C Check that the column do not already exist and if not, creates C a new one. C CALL TBLSER (TID,LABBK,POSBK,STAT) IF (POSBK.GT.0) THEN STRING = 'Warning : Column label '//LABBK//'already exists' CALL STTPUT (STRING,STAT) STRING = ' Not created' CALL STTPUT (STRING,STAT) ELSE TYP = D_R4_FORMAT FORM = 'E12.3' UNIT = 'PIXELS' CALL TBCINI (TID,TYP,1,FORM,UNIT,LABBK,POSBK,STAT) ENDIF C C --- Displays columns positions C C IF (DBG) THEN C STRING = LABX(1:16)//LABY(1:16)//LABBK(1:16) C CALL STTPUT (STRING,STAT) C WRITE (STRING,200) POSX,POSY,POSBK C 200 FORMAT (3(I3,13X)) C CALL STTPUT (STRING,STAT) C ENDIF C C --- 5. Prepares array BACK by reading selected position values in the table C and measures local background level on these positions. C ABSROW=0 DO 100 ROW = 1,NRINT C C --- 5.1. Read selection flag C CALL TBSGET (TID,ROW,FLAG,STAT) IF (STAT.NE.0) THEN WRITE (STRING,350) ROW,ORDER 350 FORMAT ('*** Fatal : Cannot access row ',I5,' of table ',A40) CALL STETER (9,STRING) ENDIF C IF (FLAG) THEN ABSROW = ABSROW + 1 IF (ABSROW.GT.NPOINT) THEN STRING = 'Do you really want to process such a big table ???' CALL STTPUT (STRING,STAT) STRING= '*** Warning : Buffer overflow in module FITBACK.FOR' CALL STETER (9,STRING) ENDIF C C --- 5.2. Read positions C CALL TBERDR (TID,ROW,POSX, BACK(1,ABSROW),NULL,STAT) CALL TBERDR (TID,ROW,POSY, BACK(2,ABSROW),NULL,STAT) CALL TBERDR (TID,ROW,PSXSTA, XSTA,NULL,STAT) CALL TBERDR (TID,ROW,PSXEND, XEND,NULL,STAT) CALL TBERDR (TID,ROW,PSYSTA, YSTA,NULL,STAT) CALL TBERDR (TID,ROW,PSYEND, YEND,NULL,STAT) C C --- 5.3. Measures background C BACK (3,ABSROW) = MEASUR (MADRID(PNTRA),NPIX(1),NPIX(2), + XSTA,XEND,YSTA,YEND,WIDTH,MMODE) C C --- 5.4. Write measured background on the table C CALL TBEWRR (TID,ROW,POSBK,BACK(3,ABSROW),STAT) C ENDIF C IF (DBG) THEN C WRITE (STRING,300) ROW,FLAG,ABSROW,(BACK(COL,ABSROW),COL=1,3) C300 FORMAT(3I5,3F12.3) C CALL STTPUT (STRING,STAT) C ENDIF 100 CONTINUE C C --- 5.5 If DEGREE<0 EXIT C IF (DEGREE.LT.0) GOTO 1000 C C --- 6. Get output frame C IF (REFIMA.EQ.BAKIMA) THEN CALL STETER (2,'Input and Output frames must be different...') ENDIF C CALL STIPUT (BAKIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START,STEP,IDENT,CUNIT, + PNTRB,IMNOB,STAT) C C --- 7. Update descriptors of the output frame C C --- 7.1 Copy all descriptors from the input frame to the output frame C CALL STDCOP (IMNOA,IMNOB,1,' ',STAT) C C --- 7.2. Write specific IDENT on the output frame C WRITE (HISTO,800) WIDTH,SMOOTH,DEGREE 800 FORMAT ('Background estimate - Half-width:',I3,' Smooth:',E13.4, + ' Degree:',I3) CALL STTPUT (HISTO,STAT) C C CALL STDWRC (IMNOB,'ORIGIN',1,REFIMA,1,60,DUN,STAT) CALL STDWRC (IMNOB,'HISTORY',1,HISTO,1,80,DUN,STAT) C C CALL STDWRC (IMNOB,'HISTORY',1,HISTO,-1,80,DUN,STAT) C. ---> STDWRC does not work properly for felem = -1. (17.07.90) C C C --- 8. Call the routine which actually process the data C C - Check initial conditions C - Do the smoothing spline interpolation along the columns. C - Interpolates background row by row for residual undefined positions C - Writes interpolated values on the frame C - Returns an error code (0 = O.K.) C !!! The subroutine modifies the value of SMOOTH !!! C CALL FITBAK (NPOINT,BACK,ABSROW,MADRID(PNTRB),NPIX(1),NPIX(2), + DEGREE,SMOOTH,ERROR,CUTS) C IF (ERROR.NE.0) THEN STRING = '*** An error occured during execution' CALL STETER (9,STRING) ENDIF C C --- 8.2 Finish descriptor updating C CALL STDWRR (IMNOB,'LHCUTS',CUTS,1,4,DUN,STAT) C C --- 9. Release files, update keywords and exit C 1000 CONTINUE CALL STSEPI C STOP END SUBROUTINE FITBAK (DIMTAB,TABLE,NBVAL,IMAGE,NX,NY,DEGREE,SMOOTH, + ERR,CUTS) C C C C DESCRIPTION : C C This subroutine do the actual work of inter- or extrapolation C using a smoothing spline algorithm to fit the data in mainly C two steps : C C 1. First, assuming that background positions are defined column C by column, one fits the data for all defined columns. This operation C involves the following steps. C C - Search for a defined column in the table and load X,Y and W arrays C which represents respectively the row position, the background C measured value and the relative weight of the value (here W=1.0). C - Updates the array of defined columns (TABFLG) C - If necessary, forces the values on the edges of the column C - Computes knots and b-spline coefficients C - Computes interpolated values along the whole column and writes C them on the frame. C C 2. Residual undefined pixels are interpolated using a row by row C spline interpolation. This time, one forces the spline interpolation C to pass exactly by the already defined values (smoothing factor s=0.). C This is done in a similar way of the previous step : C C - In order to accelerate computation : C One loads an array of residual undefined columns positions (UNDCOL) C Since this array is constant, so are the knots positions of the C spline interpolation, which are computed only once, for the C first row. C - Computes the ratio of the job already done and displays it. C - Computes b-spline coefficients for each row C - Computes interpolated values along the whole row and writes C them on the frame. C C If something is wrong during the execution, an error message is C displayed, trying to fix the problem. An error code is returned to C the driver program and the execution immediately stopped. C C C INPUT/OUTPUT : C C ERR is the only output variable C IMAGE is an input/output variable C All others variables are input data C C DIMTAB INTEGER Second dimension of the array TABLE C (First dimension is supposed to be 3) C TABLE REAL array Col. 1 : X posiitons of background C Col. 2 : Y posiitons of background C Col. 3 : measured values of background C NBVAL INTEGER Number of significant data in array TABLE C IMAGE REAL array Adress of the output image C NX INTEGER X-dimension of the output image C NY INTEGER Y-dimension of the output image C DEGREE INTEGER Degree of the spline polynomials C SMOOTH REAL Smoothing factor of the spline interpolation C ERR INTEGER Error return code (0 = O.K.) C C VERSIONS : C C V. 0.00 27.06.90 Creation C V. 1.00 17.06.90 Reference version C V. 1.10 19.07.90 Remove MASK, introduce CUTS computation C Speed improv. : (use PREPOS, remove IMAGE init. and W init) C C --- 1. Variables declarations C IMPLICIT NONE C C --- 1.1. Input/Output variables C INTEGER DIMTAB REAL TABLE(3,DIMTAB),CUTS(4) INTEGER NBVAL,NX,NY,DEGREE,ERR REAL IMAGE(1),SMOOTH C C --- 1.2. Internal variables C INTEGER M,NEST,LEVEL PARAMETER (M=10000,NEST=4096) REAL X(M),Y(M),W(M),Q(M,6),XB,XE,KNOTS(NEST),COEF(NEST),FP REAL NRDATA(NEST),FPO,FPOLD,NPLUS,THRES INTEGER IOPT,NBKNOT,IER,NBUND,PREPOS REAL ARG,DERIV,UNDCOL(M),MEAN,RATIO INTEGER PTIMA,PTTAB,STAT,FLAG,NK1,NU,PASS1 INTEGER ROW,COL,POSARG,POSCOL,POSTAB,NPTS,TABFLG(M,2) INTEGER INTCOL,INTROW REAL BIGSMO CHARACTER STRING*80 C C --- 2. Initialisations C IF (NX.GT.M.OR.NY.GT.M) THEN STRING = 'Do you really want to process such a big image ???' CALL STTPUT (STRING,STAT) STRING = '*** Fatal : Parameter M in module FITBAK is too small' CALL STTPUT (STRING,STAT) IF (NY.GT.NX) NX = NY WRITE (STRING,5100) NX,M 5100 FORMAT ('Minimum value is : ',I8,' instead of :',I8) CALL STTPUT (STRING,STAT) ERR = 10 ! No "science" error RETURN ENDIF C C DO 50 COL = 1 , M TABFLG (COL,1) = 0 ! Number of data for interpolation TABFLG (COL,2) = 0 ! Interpolation done (=1) W (COL) = 1.0 ! Weights of values 50 CONTINUE ERR = 0 ! Normal return LEVEL = 2 ! Level of messages when calling IERSPL C SMOOTH = 0. ! Smoothing factor for this part of the algorithm INTCOL = 0 INTROW = 0 BIGSMO = 1.0E+20 C C --- 3. Building the table of values C 100 CONTINUE ! Beginning of the loop on table values C C --- 3.1. Identification of a new column position not yet interpolated C FLAG = 0 DO 110 PTTAB = 1 , NBVAL COL = NINT (TABLE(1,PTTAB)) IF (TABFLG(COL,1).EQ.0) THEN FLAG = 1 POSTAB = PTTAB GOTO 150 ENDIF 110 CONTINUE 150 CONTINUE C C --- 3.2. If FLAG#0, loading of table data for the given column number C IF (FLAG.EQ.0) THEN cd STRING = '*** No new columns' cd CALL STTPUT (STRING,STAT) GOTO 1000 ELSE NPTS = 0 ! Number of ORDER table values for column COL DO 200 PTTAB = POSTAB , NBVAL C POSCOL = NINT (TABLE(1,PTTAB)) IF (COL.EQ.POSCOL) THEN C NPTS = NPTS + 1 IF (NPTS.GT.M) THEN STRING = '*** Fatal : Parameter M in module FITBAK is too small' CALL STTPUT (STRING,STAT) WRITE (STRING,5100) NPTS,M ERR = 10 ! No "science" error RETURN ENDIF X(NPTS) = TABLE (2,PTTAB) ! Y position, in pixels Y(NPTS) = TABLE (3,PTTAB) ! Background value cd W(NPTS) = 1.0 ! Weight ENDIF 200 CONTINUE ENDIF C C --- 3.3. Update the array of interpolated columns and check number of data C TABFLG(COL,1) = NPTS cd WRITE (STRING,5000) COL,TABFLG(COL,1) cd CALL STTPUT (STRING,STAT) cd5000 FORMAT ('Process column : ',I5,' (',I4,' pixels)') C IF (TABFLG(COL,1).LE.DEGREE) THEN STRING = '*** Warning : Number insufficient for interpolation' CALL STTPUT (STRING,STAT) GOTO 100 ! Search for another column number ENDIF C C --- 4. Smoothing spline interpolation along columns C IOPT = 0 ! Restart all computations C C --- 4.1. Evaluating minimum abscisse XB and maximum XE and eventually, C forces the values on the edges of the column C C CALL SORSPL (X,Y,W,NPTS) ! If it is necessary to sort the data C IF (NINT(X(1)).NE.1) THEN DO 350 PTTAB = NPTS , 1 , -1 X(PTTAB+1) = X(PTTAB) Y(PTTAB+1) = Y(PTTAB) cd W(PTTAB+1) = W(PTTAB) ! Speed optim. 350 CONTINUE X(1) = 1.0 Y(1) = Y(2) cd W(1) = W(2) ! Speed optim. NPTS = NPTS + 1 ENDIF IF (NINT(X(NPTS)).NE.NY) THEN X(NPTS+1) = FLOAT(NY) Y(NPTS+1) = Y(NPTS) cd W(NPTS+1) = W(NPTS) ! Speed optim. NPTS = NPTS + 1 ENDIF XB = 1.0 XE = FLOAT(NY) C C --- 4.2. Computation of knots sequence and b-spline coefficients C CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,SMOOTH,NBKNOT, + KNOTS,COEF,FP,IOPT,IER,NRDATA,FPO,FPOLD,NPLUS) C C --- 4.3. Error management C IF (IER.EQ.3) THEN ! No convergence achieved CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,BIGSMO,NBKNOT, + KNOTS,COEF,FP,0,IER,NRDATA,FPO,FPOLD,NPLUS) PASS1 = 0 INTCOL = INTCOL + 1 ENDIF CALL IERSPL (IER,LEVEL) ! Displays error messages IF (IER.GT.0) THEN ERR = 1 ! Error during spline interpolation 5200 FORMAT ('Error while processing column ',I6) WRITE(STRING,5200) POSCOL CALL STTPUT(STRING,STAT) RETURN ELSE TABFLG (COL,2) = 1 ! Interpolation correctly done ENDIF C C --- 5. Computation of interpolated values C NU = 0 ! Order of the derivative PREPOS = DEGREE + 1 DO 400 ARG = XB , XE , 1.0 C C --- 5.1. Find the position of the argument within the knots sequence C CALL ARGSPL (ARG,KNOTS,NBKNOT,DEGREE,PREPOS,POSARG,NK1) C C --- 5.2. Computes interpolated values C C COL is already defined within the main loop ROW = NINT(ARG) PTIMA = (ROW-1)*NX+COL IMAGE(PTIMA) = DERIV (KNOTS,NBKNOT,COEF,NK1,NU,ARG,POSARG) cd WRITE (STRING,390) POSARG,ARG,COL,ROW,IMAGE(PTIMA) cd390 FORMAT (I3,F5.1,2I5,F15.2) cd CALL STTPUT (STRING,STAT) 400 CONTINUE C C --- End of loop on table values C GOTO 100 ! Search for another column number 1000 CONTINUE ! Continue the process IF (INTCOL.GT.0) THEN WRITE (STRING,395) INTCOL 395 FORMAT(1X,I5,' columns interpolated by a polynomial') CALL STTPUT(STRING,STAT) ENDIF C C --- 6. Interpolating along the rows C C --- 6.0. Initialisations C THRES = 20. PASS1 = 0 ! Indicates the first pass in the loop on ROW XB = 1.0 ! Boundaries of the interpolation XE = FLOAT(NX) C SMOOTH = 0. ! Smoothing factor for this part of the algorithm NBUND = 0 ! Initialises number of undefined columns DO 1050 COL = 1 , NX IF (TABFLG(COL,2).EQ.0) THEN NBUND = NBUND + 1 UNDCOL (NBUND) = FLOAT(COL) ENDIF 1050 CONTINUE IF (NBUND.EQ.0) RETURN ! All columns are defined. C That's no more work to do... DO 1100 ROW = 1 , NY C C --- C RATIO = 100.*FLOAT(ROW)/FLOAT(NY) IF (RATIO.GE.THRES) THEN WRITE (STRING,1150) THRES 1150 FORMAT ('Background interpolation:',F5.1,' % done') CALL STTPUT (STRING,STAT) THRES = THRES + 20.0 ENDIF NPTS = 0 DO 1200 COL = 1 , NX C C --- 6.1. Load values for the defined columns in arrays X,Y,W C IF (TABFLG(COL,2).NE.0) THEN NPTS = NPTS + 1 X(NPTS) = FLOAT (COL) PTIMA = (ROW-1)*NX+COL Y(NPTS) = IMAGE(PTIMA) cd W(NPTS) = 1.0 ! Speed optim. cd TYPE*,NPTS,X(NPTS),Y(NPTS) ENDIF 1200 CONTINUE C C --- 6.2. Computes b-splines coefficients. C IF (PASS1.EQ.0) THEN ! On the first pass, ... IOPT = 0 ! Restart all computations, ... PASS1 = 1 ELSE ! then ... IOPT = 1 ! starts with the knots found at the ENDIF ! last call. CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,SMOOTH,NBKNOT, + KNOTS,COEF,FP,IOPT,IER,NRDATA,FPO,FPOLD,NPLUS) C C --- 6.3. Error management C IF (IER.EQ.3) THEN ! No convergence achieved CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,BIGSMO,NBKNOT, + KNOTS,COEF,FP,0,IER,NRDATA,FPO,FPOLD,NPLUS) PASS1 = 0 INTROW = INTROW + 1 ENDIF CALL IERSPL (IER,LEVEL) ! Displays error messages IF (IER.GT.0) THEN ERR = 1 ! Error during spline interpolation C5300 FORMAT ('Error while processing row ',I6) WRITE(STRING,5200) ROW CALL STTPUT(STRING,STAT) RETURN ENDIF C C --- 6.4. Computes interpolated values C (Only where it is not already defined) C NU = 0 ! Order of the derivative PREPOS = DEGREE + 1 C DO 1400 NPTS = 1 , NBUND C ARG = UNDCOL (NPTS) ! The argument is an undefined column position DO 1400 ARG = 1.0 , NX C C --- 6.4.1. Find the position of the argument within the knots sequence C CALL ARGSPL (ARG,KNOTS,NBKNOT,DEGREE,PREPOS,POSARG,NK1) C C --- 6.4.2. Computes interpolated values C C ROW is already defined within the main loop COL = NINT(ARG) PTIMA = (ROW-1)*NX+COL IMAGE(PTIMA) = DERIV (KNOTS,NBKNOT,COEF,NK1,NU,ARG,POSARG) 1400 CONTINUE ! End of the loop on ARG 1100 CONTINUE ! End of the loop on ROW IF (INTROW.GT.0) THEN WRITE (STRING,396) INTROW 396 FORMAT(1X,I5,' rows interpolated by a polynomial') CALL STTPUT(STRING,STAT) ENDIF C C --- 7. CUTS computation C CUTS(3) = IMAGE(1) ! Minimum value CUTS(4) = IMAGE(1) ! Maximum value DO 2000 PTIMA = 2 , NX*NY IF (IMAGE(PTIMA).LT.CUTS(3)) CUTS(3)=IMAGE(PTIMA) IF (IMAGE(PTIMA).GT.CUTS(4)) CUTS(4)=IMAGE(PTIMA) 2000 CONTINUE MEAN = (CUTS(3)+CUTS(4))/2. CUTS(3) = MEAN - (MEAN-CUTS(3))*1.0 CUTS(4) = MEAN + (CUTS(4)-MEAN)*1.0 CUTS(1) = CUTS(3) CUTS(2) = CUTS(4) C C --- 8. That's all folks C RETURN END FUNCTION MEASUR (IMAGE,NX,NY,XSTA,XEND,YSTA,YEND,WIDTH,MMODE) C C ---------------------------------------------------------------------------- C C. IDENTIFICATION C C Author : P. Ballester ESO - Garching C C Versions : 1.00 30.06.90 C C. PURPOSE C C This function measures on the frame the median value of the signal C on the neighbourhood of the position (POSX,POY), inside C a sub-frame of dimension (2*WIDX+1,2*WIDY+1) C C. RESTRICTIONS : C C (2*WIDX+1)*(2*WIDY+1) < 225 (Parameter MAXPIX) C Executable only within MIDAS environment (STTPUT,STETER) C C. INPUTS : C C IMAGE Frames on which the signal is measured C NX,NY Dimensions of the frame C POSX,POSY Central position of the sub-image C WIDX X Half-size of the window C WIDY Y Half-size of the window C C. OUTPUTS : C C MEASUR Mean value of the signal in the sub-image C C. ALGORITHM : C C One checks that the central position is defined in the frame. Then C the pixels of the neighbourhood of the central position are C summed up (excepted the ones which are not defined in the C frame) and the median value is returned. C C ---------------------------------------------------------------------------- C C --- 1. Variables definition C IMPLICIT NONE C LOGICAL FLAG INTEGER MAXPIX,NPIX,MEDMIN,MEDMAX,KC,MMODE PARAMETER (MAXPIX=1024) REAL MEASUR,IMAGE(1),HISTO(MAXPIX),BUFFER INTEGER NX,NY,WIDTH INTEGER PTIMA,COL,ROW,POSY,KR, STP, NXSTA, NXEND REAL XSTA,XEND,YSTA,YEND,SLOPE,INTER CHARACTER STRING*80 C C --- 3. Computes mean value C SLOPE = (YEND-YSTA)/(XEND-XSTA) INTER = ((YEND+YSTA)-SLOPE*(XEND+XSTA))/2. STP = 1 IF (XSTA.GT.XEND) THEN NXSTA = NINT(XEND) NXEND = NINT(XSTA) ELSE NXSTA = NINT(XSTA) NXEND = NINT(XEND) ENDIF MEASUR = 0. NPIX = 0 DO 100 COL = NXSTA, NXEND POSY = NINT(COL*SLOPE+INTER) DO 50 KR = -WIDTH,WIDTH ROW = POSY + KR IF (COL.LT.1.OR.COL.GT.NX.OR.ROW.LT.1.OR.ROW.GT.NY) GOTO 50 NPIX = NPIX + 1 IF (NPIX.GT.MAXPIX) THEN STRING = '*** Fatal : Buffer HISTO is too small in MEASUR' CALL STETER (9,STRING) ENDIF PTIMA = (ROW-1)*NX + COL HISTO(NPIX) = IMAGE (PTIMA) 50 CONTINUE 100 CONTINUE C C --- 4. Sort the values C IF (NPIX.EQ.1) THEN MEASUR = HISTO(1) RETURN ENDIF 150 CONTINUE FLAG = (.FALSE.) DO 200 KC = 1, NPIX-1 IF (HISTO(KC).GT.HISTO(KC+1)) THEN BUFFER = HISTO (KC+1) HISTO(KC+1) = HISTO(KC) HISTO(KC) = BUFFER FLAG = (.TRUE.) ENDIF 200 CONTINUE IF (FLAG) GOTO 150 C C --- 5.a) Return the MINIMUM if desired C IF (MMODE.EQ.1) THEN MEASUR = HISTO(1) RETURN ENDIF C C --- 5.b) Define the median value C IF (2*NINT(NPIX/2.).EQ.NPIX) THEN MEDMIN = NPIX/2 MEDMAX = NPIX/2 + 1 MEASUR = (HISTO(MEDMIN)+HISTO(MEDMAX))/2. ELSE MEDMIN = (NPIX-1)/2+1 MEASUR = HISTO(MEDMIN) ENDIF C write(6,*) 'made it till here and measured ', MEASUR RETURN END