/* $Id: cpl_image_stats_body.h,v 1.21 2007/11/07 13:09:16 llundin Exp $
 *
 * This file is part of the ESO Common Pipeline Library
 * Copyright (C) 2001-2004 European Southern Observatory
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/* Type dependent macros */
#if CPL_CLASS == CPL_CLASS_DOUBLE
#define CPL_TYPE double
#define CPL_TYPE_T CPL_TYPE_DOUBLE

#elif CPL_CLASS == CPL_CLASS_FLOAT
#define CPL_TYPE float
#define CPL_TYPE_T CPL_TYPE_FLOAT

#elif CPL_CLASS == CPL_CLASS_INT
#define CPL_TYPE int
#define CPL_TYPE_T CPL_TYPE_INT

#else
#undef CPL_TYPE
#undef CPL_TYPE_T
#endif

#define CPL_TYPE_ADD(a) CPL_CONCAT2X(a, CPL_TYPE)

#if CPL_OPERATION == CPL_IMAGE_STATS_MEDIAN_SUBW

    case CPL_TYPE_T:
    {
        /* Point to first pixel in first row to read */
        const CPL_TYPE   * pi = (const CPL_TYPE*)image->pixels
            + (lly-1)*image->nx;
        /* Number of pixels inside window */
        const int          npix = (urx-llx+1) * (ury-lly+1);
        /* Duplicate the pixels inside the window :-( */
        CPL_TYPE         * copybuf = cpl_malloc(npix * sizeof(CPL_TYPE));

        if (image->bpm == NULL ||
            cpl_mask_is_empty_window(image->bpm, llx, lly, urx, ury)) {

            /* All pixels are good */
            ngood = npix;

            /* Cannot fail at here */
            (void)cpl_tools_copy_window((void*)copybuf,
                                        image->pixels,
                                        sizeof(CPL_TYPE),
                                        image->nx, image->ny,
                                        llx, lly, urx, ury);

        } else {
            /* Point to first pixel in first row to read */
            const cpl_binary * pbpm
                = cpl_mask_get_data_const(image->bpm) + (lly-1)*image->nx;
            int i, j;

            for (j = lly - 1; j < ury; j++, pi += image->nx, pbpm += image->nx) {
                for (i = llx - 1; i < urx; i++) {
                    /* Take only good pixels */
                    if (pbpm[i] == CPL_BINARY_0) {
                        copybuf[ngood++] = pi[i];
                    }
                }
            }
        }
        if (ngood > 0) {
            /* There are good pixels */
            /* Compute the median */
            median = CPL_TYPE_ADD(cpl_tools_get_median)(copybuf, ngood);
        }
        cpl_free(copybuf);
    
        break;
    }

#elif CPL_OPERATION == CPL_IMAGE_STATS_MEDIAN_STAT

  case CPL_TYPE_T:
  {
      CPL_TYPE * pi;

      pi = (CPL_TYPE*)image->pixels ;
      for (i=0 ; i<npix ; i++)
          if (!badmap[i]) *sigma+= fabs((double)(pi[i]-median_val));
      break ;
  }


#endif

#undef CPL_TYPE
#undef CPL_TYPE_T