/* @(#)getdat.c	17.1.1.1 (ESO-DMD) 01/25/02 17:49:05 */
/*===========================================================================
  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 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding 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   (c) 1993 European Southern Observatory
.IDENTifer   GETDAT
.AUTHOR      R.M. van Hees IPG-ESO Garching
.KEYWORDS    low level plot routine
.LANGUAGE    C
.PURPOSE     extract a subimage from a FRAME and smoothes it (option)
        input: int   imf      file id of related MIDAS frame
               int   maxsize  maximum number of pixels stored
               int   *npix    standart descriptor of the frame
               float *image   contains in pixel units
               int   ism      smooting parameter
    in/output: float *p_img   pointer to the output frame

.COMMENTS    only two dimensional frames
             uses (little bit more than): 4 * 2 * maxsize Mbyte memory
             
.ENVIRONment MIDAS
             #include <plot_def.h>      Symbols used by the PLT interfaces
             #include <midas_def.h>     Prototypes for MIDAS interfaces

.VERSION     1.1     14-Dec-1993   rewritten: reads in chunks, RvH
             1.0     26-Nov-1993   Created by R.M. van Hees
------------------------------------------------------------*/
/*
 * Define _POSIX_SOURCE to indicate
 * that this is a POSIX program
 */
#define  _POSIX_SOURCE 1

/*
 * definition of the used functions in this module
 */
#include <math.h>

/*
 * define some macros and constants
 */
#include <midas_def.h>
#include <plot_def.h>

/*
 * here start the code of the function
 */
void GETDAT( imf, maxsize, npix, image, ism, p_img )
int   imf, maxsize, *npix, ism;
float *image, *p_img;

{
register int ic, ix, iy;

int    actvals, chunks, felem, imdum, ism1, ism2, navrg, nc, nr, nrcol, nrline,
       nrrow, size, ifram[4];
float  avrg, *pntr;
char   *cpntr;

if ( ism < 0 ) ism = 0;

/*
 * get sub-frame dimensions
 */
ifram[0] = NINT( MYMIN( image[0], image[1] ));
ifram[1] = NINT( MYMAX( image[0], image[1] ));
ifram[2] = NINT( MYMIN( image[2], image[3] ));
ifram[3] = NINT( MYMAX( image[2], image[3] ));

/*
 * determine  size of the sub image
 */
nrcol  = ifram[1] - ifram[0] + 1;
nrrow  = ifram[3] - ifram[2] + 1;

/*
 * we can extract NRLINes of data from the original
 */
nrline = (int) floor( (double) maxsize / (double) *npix );
nrline = MYMIN( nrrow + 2 * ism, nrline ); 

/*
 * this means effectively that the number of rows and chunks equals
 */
if ( (nrrow = nrline - 2 * ism) < 1 ) 
   SCETER( 1, "*** FATAL: GETDAT, maxsize too small given the smooth factor" );
chunks = (int) ceil( (double) (ifram[3] - ifram[2] + 1) / nrrow );

/*
 * allocate virual memory and scratch space
 */
size  = nrline * *npix;
(void) SCFCRE( "DUMMY", D_R4_FORMAT, F_X_MODE, F_IMA_TYPE, size, &imdum );
(void) SCFMAP( imdum, F_X_MODE, 1, size, &actvals, &cpntr );

if ( ism == 0 )
   { felem = *npix * (ifram[2] - 1) + 1;
     for ( ic = 0; ic < chunks; ic++ )
         { (void) SCFGET( imf, felem, size, &actvals, cpntr ); 
          
           pntr = (float *) cpntr;
           pntr += ifram[0] - 1;
           for ( nr = 0; nr < nrrow; nr++ )
               { for ( nc = 0; nc < nrcol; nc++ ) *p_img++ = *(pntr + nc);

                 pntr += *npix;
               }

/*
 get next first pixel to read, number of rows and siz eof the data block
 */
           felem += size;
           nrrow = MYMIN((ifram[3]-ifram[2]+1) - (ic+1) * nrrow, nrline );
           size  = nrrow * *npix;
         }
   }
else
   { for ( ic = 0; ic < chunks; ic++ )
         { felem = *npix * MYMAX( ifram[2] + ic * nrrow - ism - 1, 0 ) + 1;
           (void) SCFGET( imf, felem, size, &actvals, cpntr ); 

           ism1 = MYMIN( ism, ifram[2] - 1 + ic * nrrow );
           ism2 = MYMIN( ism, npix[1] - ifram[2] + 1 - (ic+1) * nrrow );
           nrline = nrrow + ism1 + ism2;
           nrrow = MYMIN((ifram[3] - ifram[2]+1) - ic * nrrow, nrrow );
          
           pntr = (float *) cpntr;
           pntr += ism1 * *npix + ifram[0] - 1;
           for ( nr = 0; nr < nrrow; nr++ )
               { for ( nc = 0; nc < nrcol; nc++ ) 
                     { navrg = 0;
                       avrg  = 0.0;
                       for ( iy = -ism; iy <= ism; iy++ )
                           { if ( nr+iy+ism1 >= 0 && nr+iy < nrline )
                                { for ( ix = -ism; ix <= ism; ix++ ) 
                                      { if ( nc+ix >= 0 && nc+ix < *npix )
                                           { register int k = nc+ix+iy* *npix;
                                             avrg += *(pntr + k);
                                             navrg++;
                                           }
                                      }
                                 }
                           }
                       if ( navrg > 0 )
                          *p_img++ = avrg / navrg;
                       else
                          *p_img++ = 0.0;
                     }
                 pntr += *npix;
               }
         }
   }
/*
 * release the allocated memory
 */
(void) SCFCLO( imdum );

return;
}