/* $Id: cpl_stats_body.h,v 1.6 2006/08/23 12:24:31 yjung 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


#if CPL_OPERATION == CPL_IMAGE_STATS_VARIANCE

  case CPL_TYPE_T:
  {
      CPL_TYPE * pi;

      double delta ;
      double ipix = 0; /* Counter of pixels used */

      /* Compute the variance of the elements using the recurrence relation
         var_n = var_{n-1} + ((x_n-mean)*(x_n-mean) - var_{n-1})/(n+1) */
      pi = (CPL_TYPE*)image->pixels ;
      for (j=lly-1 ; j<ury ; j++) {
          pos = (llx-1)+j*image->nx ;
          for (i=llx-1 ; i<urx ; i++) {
              if (badmap == NULL || !badmap[pos]) {
                  delta = pi[pos] - pix_sum / stats->npix;
                  variance += (delta * delta - variance) / ++ipix;
              }
              pos++;
          }
      }
      break ;
  }

#elif CPL_OPERATION == CPL_IMAGE_STATS_CENTROID

  case CPL_TYPE_T:
  {
      CPL_TYPE * pi;

      if (min_pix < 0)  min_pix_tmp = min_pix ;
      else              min_pix_tmp = 0.0 ;
      pi = (CPL_TYPE*)image->pixels ;
      for (j=lly-1 ; j<ury ; j++) {
          pos = (llx-1)+j*image->nx ;
          for (i=llx-1 ; i<urx ; i++) {
              if (badmap == NULL || !badmap[pos]) {
                  sum_xz += ((double)pi[pos]-min_pix_tmp)*(double)(i+1) ;
                  sum_yz += ((double)pi[pos]-min_pix_tmp)*(double)(j+1) ;
                  sum_z += (double)pi[pos]-min_pix_tmp ;
                  sum_x  += (double)(i+1) ;
                  sum_y  += (double)(j+1) ;
                  sum_pix += 1.0 ;
              }
              pos++;
          }
      }
      if (sum_z  < 0) sum_z  = 0; /* Can only become negative due to rounding */
      if (sum_xz < 0) sum_xz = 0; /* Can only become negative due to rounding */
      if (sum_yz < 0) sum_yz = 0; /* Can only become negative due to rounding */
      break ;
  }

#elif CPL_OPERATION == CPL_IMAGE_STATS_MINMAX

  case CPL_TYPE_T:
  {
      CPL_TYPE * pi;

      pi = (CPL_TYPE*)image->pixels ;

      min_pos = max_pos = firstgoodpos;
      min_pix = (double)pi[firstgoodpos] ;
      max_pix = (double)pi[firstgoodpos] ;
      for (j=lly-1 ; j<ury ; j++) {
          pos = (llx-1)+j*image->nx ;

          for (i=llx-1 ; i<urx ; i++) {
              if (badmap == NULL || !badmap[pos]) {
                 if (pi[pos] < min_pix) min_pix = (double)pi[min_pos = pos] ;
                 if (pi[pos] > max_pix) max_pix = (double)pi[max_pos = pos] ;
              }
              pos++;
          }
      }
      break ;
  }

#elif CPL_OPERATION == CPL_IMAGE_STATS_MEAN

  case CPL_TYPE_T:
  {
      CPL_TYPE * pi;

      pi = (CPL_TYPE*)image->pixels ;

      for (j=lly-1 ; j<ury ; j++) {
          pos = (llx-1)+j*image->nx ;

          for (i=llx-1 ; i<urx ; i++) {
              if (badmap == NULL || !badmap[pos]) {
                  pix_sum += (double)pi[pos] ;
                  abs_sum += (double)fabs((double)pi[pos]) ;
                  sqr_sum += ((double)pi[pos]) * ((double)pi[pos]) ;
              }
              pos++;
          }
      }

      break ;
  }

#elif CPL_OPERATION == CPL_IMAGE_STATS_MEANMINMAX

  case CPL_TYPE_T:
  {
      CPL_TYPE * pi;

      pi = (CPL_TYPE*)image->pixels ;
 
      min_pos = max_pos = firstgoodpos;
      min_pix = (double)pi[firstgoodpos] ;
      max_pix = (double)pi[firstgoodpos] ;
      for (j=lly-1 ; j<ury ; j++) {
          pos = (llx-1)+j*image->nx ;

          for (i=llx-1 ; i<urx ; i++) {
              if (badmap == NULL || !badmap[pos]) {
                  if (pi[pos] < min_pix) min_pix = (double)pi[min_pos = pos] ;
                  if (pi[pos] > max_pix) max_pix = (double)pi[max_pos = pos] ;
                  pix_sum += (double)pi[pos] ;
                  abs_sum += (double)fabs((double)pi[pos]) ;
                  sqr_sum += ((double)pi[pos]) * ((double)pi[pos]) ;
              }
              pos++;
          }
      }

      break ;
  }

#endif

#undef CPL_TYPE
#undef CPL_TYPE_T