/* @(#)medfilt.c	17.1.1.1 (ESO-DMD) 01/25/02 17:11:16 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of
  the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public
  License along with this program; if not, write to the Free
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
  MA 02139, USA.

  Correspondence concerning ESO-MIDAS should be addressed as follows:
        Internet e-mail: midas@eso.org
        Postal address: European Southern Observatory
                        Data Management Division
                        Karl-Schwarzschild-Strasse 2
                        D 85748 Garching bei Muenchen
                        GERMANY
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT:  Copyright (c) 1987 European Southern Observatory,
                                         all rights reserved
.IDENTIFIER  module medfilt
.LANGUAGE    C
.AUTHOR      K. Banse                   IPG-ESO Garching
.KEYWORDS
.COMMENTS    

MEDFILT: Median filtering of single frame, general, v.2.1, NOV 97
  See: T.B.Georgiev, 1996, Bull.SAO (Russia) 39, 124
  MEDFILT  inp[.fts]  smt[.fts]  rsd[.fts]  W  C  Z 
  inpfr - input frame; 
  outfr - result (smoothed) frame;
  rsdfr - residual frame,  rsdfr = inpfr - outfr;
  W: window diameter, W=2/3/4/5..,99          DEFAULT W=25
  C=3-4: coeff. for deviation threshold -- the median M
    changes the current pixel value V if abs(V-M) > C*sigma
    the case of simple median filtering is:  DEFAULT C=0  
  T.B.Georgiev, IA Sofia, tsgeorg@phys.acad.bg

.VERSIONS    1.00       980708  implemented in Midas.  KB
 
------------------------------------------------------------------*/

#include <stdio.h>
#include <math.h>
#include <midas_def.h>

#define   MAXNC  8192   /* max image width in pixels */
#define   MAXFW  99  /*max window diameter in pixels */
#define   NRB    40  /* number of output buffer rows */
#define   MAXHI  65535         /* max histogram size */
#define   MAXIN  32767

/*

*/

main()

{
int  finp, fout, frsd, iav, uni, nul;
int  nr, nc, nim, fw, fw0, hw, sum, med, hlf, npix, nb, ib, maxhi;
int  me1, me2, su1, su2, q1, q2, ncor1, ncor2, iii;
int  niout, i, iinp, iout, j, k, l, ll, lk, n, np, npc, rr, ls, lr;

static short int m[MAXFW][MAXNC+MAXFW-1];             /* median frame band */
static short int ms[NRB][MAXNC];                      /* output-smt buffer */
static short int mr[NRB][MAXNC];                      /* output-rsd buffer */
static short int h[MAXHI];                      /* histogram in the window */
static short int r[MAXNC+MAXFW-1];            /* input or output image row */
static short int pn[MAXFW];               /*row number in the memory m[][] */
static short int lw[MAXFW/2+1];             /* limits of the circle window */

double	d, dls, dlr, dp, sm, sd, ssm, ssd, sig, csig, ssig, sssig;

float  value;

char   cbuf[88], infile[82], outfile[82], rsdfile[82];




(void) SCSPRO("medfilt");

(void) SCKGETC("IN_A",1,80,&iav,infile);
(void) SCKGETC("OUT_A",1,80,&iav,outfile);
(void) SCKGETC("OUT_B",1,80,&iav,rsdfile);
(void) SCKRDI("NUM_OUT",1,1,&iav,&nim,&uni,&nul);
(void) SCKRDI("FW",1,1,&iav,&fw,&uni,&nul);
(void) SCKRDD("CSIG",1,1,&iav,&csig,&uni,&nul);

if (fw < 2) fw=2; 
if (fw > MAXFW) fw=MAXFW;
if (csig < 0.) csig=0.;


/* --------------------- PRELIMINARY PART ------------------ */

maxhi = MAXHI - 2;          /* hist.numbers: from 1 to maxhi */
fw0 = fw; 
hw = fw / 2;         
if (2 * hw == fw) 
   rr = hw * hw;  
else 
   rr = (int)((hw + .5) * (hw + .5));

npix = 1; 
for (k=0; k<=hw; k++)
   { 
   lw[k] = 0;    /* the circle window limits  */
   for (l=0; l<=hw; l++) 
      {
      if (k*k + l*l <= rr)
         {
         lw[k] = l; 
         if (k > 0) npix += 4;
         }
      }
   }

fw = 2 * hw + 1;   
hlf = npix / 2 + 1;                    /* half of the histogram */
q1 = (int)((double)npix * 0.17 + .5);   /* left  quantil of hist */
q2 = (int)((double)npix * 0.83 + .5);   /* right quantil of hist */


/* ----------------------- MAIN PROCESS ----------------------- */

finp = fits_open(infile, "r", &nr, &nc);  

if (nc > MAXNC) 
   {
   sprintf(cbuf,"Too wide frame, ncol=%d",nc);
   SCETER(11,cbuf);
   }

if (nim != 1) fout = fits_open(outfile, "w", &nr, &nc);
if (nim != 0) frsd = fits_open(rsdfile, "w", &nr, &nc);

dp = (double)nr * (double)nc;
sprintf(cbuf,"%dx%d windim=%d (%dpix) thr=%2.1f",nr, nc, fw0, npix, csig);
SCTPUT(cbuf);

me1 = me2 = npix = ncor1 = ncor2 = 0;
sm = sd = ssm = ssd = ssig = sssig = 0.;


/* ---------------------- MAIN LOOP BY IMAGE ROWS ------------------- */

for (i= -hw; i<nr+hw; i++)
   {
   iinp = i;        /* inp.row number for the inner part of the image */
   if (i < 0) 
      iinp = nr + i;   /* input row number for upper margin */
   if (i > nr - 1)
      iinp = i - nr;   /* input row number for lower margin*/

   for (k=1; k<fw; k++) pn[k-1] = pn[k];   /* rotation of numbers */

   k = i + hw; 
   np = k - k / fw * fw; 
   pn[fw-1] = np;                            /* pos.number of inp.row */
   npc = pn[hw];       /* pos.num of curr.out.row for eventual output */
   fits_get_data (finp, iinp, 0, r, nc);     /* READING THE IMAGE ROW */

   for (j=0; j<nc; j++)
      {
      ll=r[j];  
      if (ll<-MAXIN+1) ll=-MAXIN+1; 
      if (ll> MAXIN-1) ll= MAXIN-1;
      m[np][hw+j] = ll;                   /* filling the buffer m[][] */
      }

   for (j=0; j<hw; j++)
      {                                         /* left-right margins */
      m[np][j] = m[np][nc+j];  
      m[np][nc+hw+j] = m[np][hw+j]; 
      }

   iout = i - hw;                                /* output row number */
   nb = iout / NRB;                       /* num of the output buffer */
   ib = iout - nb * NRB;       /* num of the row in the output buffer */
   if (iout < 0) goto Next;


   /* ---------------- LOOP IN THE CURRENT IMAGE ROW ---------------- */

   for (j=hw; j<nc+hw; j++)
      {
      if (j == hw)
         {                         /* first pixel in the current row */
         for (n=0; n<MAXHI; n++) h[n] = 0;  /* histogram inital. */
         
         for (k=0; k<fw; k++)    
            { 
            np = pn[k]; 
            l = abs(k - hw); 
            lk = lw[l];
            for (l=hw-lk; l<=hw+lk; l++) 
               { 
               n = m[np][l] + MAXIN;  
               h[n]++;  
               } 
            }
         med = 0; 
         sum = 0; 
         while (sum < hlf)
            { 
            med++; 
            sum += h[med]; 
            }                      /* MED */

         if (csig > 0.)
            {
            me1 = med - 1; 
            su1 = sum - h[med]; 

            while (su1 > q1)
               {
               su1 -= h[me1]; 
               me1--; 
               }    /* left q-limit of the histogram body */

            me2 = med + 1; 
            su2 = sum + h[me2]; 
            while (su2 < q2) 
               {
               me2++; 
               su2 += h[me2]; 
               } 
            }		/* right q-limit of the histogram body */
         }			/* end of j==hw, i.e. the first pix */

      else 
         {                  /* other pixels in the current row */
         for (k=0; k<fw; k++) 
            {			/* along the window edge columns only */
            np = pn[k]; 
            l = abs(k - hw); 
            lk = lw[l];
            n = m[np][j-lk-1] + MAXIN;  
            h[n]--;  
            if (n <= med)  sum--;        /* removing */
            n = m[np][j+lk] + MAXIN;    
            h[n]++;  
            if (n <= med)  sum++; 
            }                            /* adding */
         while (sum >= hlf)
            { 
            sum -= h[med]; 
            med--;
            }                   /* search for median */

         while (sum < hlf)
            { 
            med++; 
            sum += h[med];
            }                             /* MEDIAN */
         if (csig > 0.)
            {
            me1 = med - 1; 
            su1 = sum - h[med]; 
            while (su1 > q1) 
               {
               su1 -= h[me1]; 
               me1--; 
               }		/* left limit of the histogram body */

            me2 = med + 1; 
            su2 = sum + h[me2]; 
            while (su2 < q2)
               {
               me2++; 
               su2 += h[me2]; 
               } 
            }			/* right limit of the histogram body */
         }				/* end of  j -- "other pixels" */


/* ------------ DIFFERENT KINDS OF OUTPUT IMAGES ------------- */

      ls = ms[ib][j-hw] = med-MAXIN;  
      lr = mr[ib][j-hw] = m[npc][j] - ls;
      dls=(double)ls;      
      dlr=(double)lr;  

      if (csig==0.)
         { 
         sm += dls;  ssm += dls * dls; 
         sd += dlr;  ssd += dlr * dlr;
         }
      if (csig > 0.) 
         {
         sig = (double)(me2 - me1) / 2.;         /* hist.sigma */
         npix++; ssig += sig;  sssig += sig * sig;  
         ll = (int)(csig * sig + 0.5);
         if (lr < -ll) ncor1++;  
         if (lr > ll) ncor2++;  
         if (-ll < lr && lr < ll) 
            {			           /* no replacing */
            ms[ib][j-hw] = m[npc][j]; 
            mr[ib][j-hw] = 0;
            dls = (double)m[npc][j]; 
            sm += dls; 
            ssm += dls * dls; 
            } 
         else 
            {				/* save replacing made earlier */ 
            sm += dls;  ssm += dls * dls;
            sd += dlr;  ssd += dlr * dlr; 
            }
         }
      }                                       /* END OF j-LOOP */

   if (ib == NRB - 1 || iout == nr - 1)
      {
      if (nim != 1) 
         {
         for (l = 0; l <= ib; l++)
            { 
            for (j = 0; j < nc; j++) r[j] = ms[l][j];
            niout = nb * NRB + l;  
            fits_put_data(fout, niout, 0, r, nc); 
            }
         }
      if (nim != 0) 
         {
         for (l = 0; l <= ib; l++) 
            { 
            for (j = 0; j < nc; j++) r[j] = mr[l][j];
            niout = nb * NRB + l;  
            fits_put_data(frsd, niout, 0, r, nc); 
            } 
         }
      }
Next:  
   ;
   }                                         /* END OF i-LOOP */

sd /= dp; 
ssd /= dp; 
ssd = sqrt(ssd - sd * sd);
sm /= dp; 
ssm /= dp; 
ssm = sqrt(ssm - sm * sm);
sprintf(cbuf,"out %3.1f+-%3.1f",sm,ssm);
SCTPUT(cbuf);
sprintf(cbuf,"resid.out %3.1f+-%3.1f",sd,ssd); 
SCTPUT(cbuf);

if (csig > 0.)
   { 
   sm = (double)ncor1 / dp; 
   sd = (double)ncor2 / dp;
   d = (double)npix; 
   ssig /= d; 
   sssig /= d; 
   sssig = sqrt(sssig - ssig * ssig); 
   sprintf(cbuf,"number of replacings: -%d (%4.3f)  +%d (%4.3f)",
                          ncor1, sm, ncor2, sd); 
   SCTPUT(cbuf);
   }  

if (nim != 1) 
   {
   (void) SCDCOP(finp,fout,1,cbuf);	/* copy all descriptors */
   CGN_DSCUPD(finp,fout," ");		/* update history */
   }

if (nim != 0) 
   {
   (void) SCDCOP(finp,frsd,1,cbuf);	/* copy all descriptors */
   CGN_DSCUPD(finp,frsd," ");		/* update history */
   }

(void) SCSEPI();		/* close all files + exit */
}