/* $Id: cpl_image_basic_body.h,v 1.63 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 defined CPL_CLASS && CPL_CLASS == CPL_CLASS_DOUBLE
#define CPL_TYPE double
#define CPL_TYPE_T CPL_TYPE_DOUBLE
#define CPL_MATH_ABS fabs
#define CPL_MATH_TYPE double

#elif defined CPL_CLASS && CPL_CLASS == CPL_CLASS_FLOAT
#define CPL_TYPE float
#define CPL_TYPE_T CPL_TYPE_FLOAT
#define CPL_MATH_ABS fabs
#define CPL_MATH_TYPE double

#elif defined CPL_CLASS && CPL_CLASS == CPL_CLASS_INT
#define CPL_TYPE int
#define CPL_TYPE_T CPL_TYPE_INT
#define CPL_MATH_ABS abs
#define CPL_MATH_TYPE int

#else
#undef CPL_TYPE
#undef CPL_TYPE_T
#undef CPL_MATH_ABS
#undef CPL_MATH_TYPE
#undef CPL_TYPE_ADD
#endif

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

#if CPL_OPERATION == CPL_IMAGE_BASIC_DECLARE

/*-----------------------------------------------------------------------------
                            Private Function prototypes
 -----------------------------------------------------------------------------*/

static cpl_image *
CPL_TYPE_ADD(cpl_image_collapse_window_create)(const cpl_image *,
                                               int, int, int, int, int);

#elif CPL_OPERATION == CPL_IMAGE_BASIC_DEFINE

/*----------------------------------------------------------------------------*/
/**
  @brief    Collapse an image region along its rows or columns.
  @param    self        Image of type CPL_TYPE to collapse.
  @param    llx         lower left x coord.
  @param    lly         lower left y coord
  @param    urx         upper right x coord
  @param    ury         upper right y coord
  @param    direction   Collapsing direction.
  @return   a newly allocated image or NULL in error case
  @note This static function will assert() on invalid image input
  @see      cpl_image_collapse_window_create()
  
  llx, lly, urx, ury are the image region coordinates in FITS convention.
  Those specified bounds are included in the collapsed region.

  The returned image must be deallocated using cpl_image_delete().

  Possible #_cpl_error_code_ set in this function:
  - CPL_ERROR_ILLEGAL_INPUT if the specified window is not valid
 */
/*----------------------------------------------------------------------------*/
static cpl_image *
CPL_TYPE_ADD(cpl_image_collapse_window_create)(const cpl_image * self,
                                               int               llx,
                                               int               lly,
                                               int               urx,
                                               int               ury,
                                               int               direction)
{

    cpl_image        * other;
    const CPL_TYPE   * pi   = (const CPL_TYPE*)self->pixels;
    const cpl_binary * bpmi = self->bpm ? cpl_mask_get_data(self->bpm) : NULL;
    CPL_TYPE         * po;
    cpl_binary       * bpmo = NULL;
    const int          n1x  = 1 + urx - llx;
    const int          n1y  = 1 + ury - lly;
    int                i, j;


    cpl_ensure(direction == 0 || direction == 1, CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(llx >= 1,        CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(lly >= 1,        CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(urx <= self->nx, CPL_ERROR_ILLEGAL_INPUT, NULL);
    cpl_ensure(ury <= self->ny, CPL_ERROR_ILLEGAL_INPUT, NULL);

    assert( self->type == CPL_TYPE_T );

    /* Let pi and bpmi point to first pixel to collapse */
    pi += self->nx * (lly-1) + (llx - 1);
    if (bpmi != NULL) bpmi += self->nx * (lly-1) + (llx - 1);

    if (direction == 0) {
        const double r1y = (double)n1y;
        int * nok;

        other = cpl_image_new(n1x, 1, CPL_TYPE_T);
        po = (CPL_TYPE*)other->pixels;

        /* To obtain a stride-1 access of pi[], some temporary storage is needed
           - this may be a disadvantage for very small images */
        nok = (int*)cpl_calloc(n1x, sizeof(int));

        for (j=0; j < n1y; j++) {
            for (i=0; i < n1x; i++) {
                if (bpmi == NULL || !bpmi[i]) {
                    po[i] += pi[i];
                    nok[i]++; /* Count good pixels */
                }
            }
            pi += self->nx;
            if (bpmi != NULL) bpmi += self->nx;
        }
        if (bpmi != NULL) {
            for (i=0; i < n1x; i++) {
                if (nok[i] == 0) {
                    /* assert(po[i] == 0.0); */
                    if (bpmo == NULL)
                        bpmo = cpl_mask_get_data(cpl_image_get_bpm(other));
                    bpmo[i] = CPL_BINARY_1;
                } else if (nok[i] < n1y) {
                    po[i] *= r1y / (double)nok[i];
                }
            }
        }
        cpl_free(nok);
    } else if (direction == 1) {
        const double r1x = (double)n1x;
        double sum;
        int    nok;

        other = cpl_image_new(1, n1y, CPL_TYPE_T);
        po = (CPL_TYPE*)other->pixels;

        for (j=0; j < n1y; j++) {
            sum = 0.0;
            nok = 0;
            for (i=0; i < n1x; i++) {
                if (bpmi == NULL || !bpmi[i]) {
                    sum += pi[i];
                    nok++;
                }
            }
            if (nok == 0) {
                /* assert(po[j] == 0.0); */
                if (bpmo == NULL)
                    bpmo = cpl_mask_get_data(cpl_image_get_bpm(other));
                bpmo[j] = CPL_BINARY_1;
            } else {
                po[j] = nok == n1x ? sum : sum * r1x / (double)nok;
            }
            pi += self->nx;
            if (bpmi != NULL) bpmi += self->nx;
        }
    } else {
        assert(0);
    }

#if CPL_CLASS != CPL_CLASS_INT
    /* FIXME: Counts also bad pixels */
    cpl_tools_add_flops(n1x * n1y);
#endif

    return other;
}


#elif CPL_OPERATION == CPL_IMAGE_BASIC_ASSIGN

    cpl_image       *   image_out ;
    register float  *   pf1,
                    *   pf2,
                    *   outpf ;
    register int    *   pi1,
                    *   pi2,
                    *   outpi ;
    register double *   pd1,
                    *   pd2,
                    *   outpd ;
    int                 i ;

    cpl_ensure(image1 && image2, CPL_ERROR_NULL_INPUT, NULL);
    /* Input data images shall have the same sizes */
    cpl_ensure(image1->nx == image2->nx && image1->ny == image2->ny,
               CPL_ERROR_ILLEGAL_INPUT, NULL);

    /* Switch on the first passed image type */
    switch (image1->type) {
        case CPL_TYPE_INT:
            image_out = cpl_image_new(image1->nx, image1->ny, CPL_TYPE_INT);
            pi1 = (int*)image1->pixels ;
            outpi = (int*)image_out->pixels ;
            /* Switch on the second passed image type */
            switch (image2->type) {
                case CPL_TYPE_INT:
                    pi2 = (int*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpi++, *pi1++, *pi2++);
                    break ;
                case CPL_TYPE_FLOAT:
                    pf2 = (float*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpi++, *pi1++, *pf2++);
                    break ;
                case CPL_TYPE_DOUBLE:
                    pd2 = (double*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpi++, *pi1++, *pd2++);
                    break ;
                default:
                    cpl_image_delete(image_out) ;
                    cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ;
            }
            break ;
        case CPL_TYPE_FLOAT:
            image_out = cpl_image_new(image1->nx, image1->ny, CPL_TYPE_FLOAT);
            pf1 = (float*)image1->pixels ;
            outpf = (float*)image_out->pixels ;
            /* Switch on the second passed image type */
            switch (image2->type) {
                case CPL_TYPE_INT:
                    pi2 = (int*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpf++, *pf1++, *pi2++);
                    break ;
                case CPL_TYPE_FLOAT:
                    pf2 = (float*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpf++, *pf1++, *pf2++);
                    break ;
                case CPL_TYPE_DOUBLE:
                    pd2 = (double*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpf++, *pf1++, *pd2++);
                    break ;
                default:
                    cpl_image_delete(image_out) ;
                    cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ;
            }
            cpl_tools_add_flops( image_out->nx * image_out->ny );
            break ;
        case CPL_TYPE_DOUBLE:
            image_out = cpl_image_new(image1->nx, image1->ny, CPL_TYPE_DOUBLE);
            pd1 = (double*)image1->pixels ;
            outpd = (double*)image_out->pixels ;
            /* Switch on the second passed image type */
            switch (image2->type) { 
                case CPL_TYPE_INT:
                    pi2 = (int*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpd++, *pd1++, *pi2++);
                    break ;
                case CPL_TYPE_FLOAT:
                    pf2 = (float*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpd++, *pd1++, *pf2++);
                    break ;
                case CPL_TYPE_DOUBLE:
                    pd2 = (double*)image2->pixels ;
                    for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                        CPL_OPERATOR(*outpd++, *pd1++, *pd2++);
                    break ;
                default:
                    cpl_image_delete(image_out) ;
                    cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ;
            }
            cpl_tools_add_flops( image_out->nx * image_out->ny );
            break ;
        default:
          cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ;
    }

    /* Handle bad pixels map */
    if (image1->bpm == NULL && image2->bpm == NULL) {
        image_out->bpm = NULL ;
    } else if (image1->bpm == NULL) {
        image_out->bpm = cpl_mask_duplicate(image2->bpm) ;
    } else if (image2->bpm == NULL) {
        image_out->bpm = cpl_mask_duplicate(image1->bpm) ;
    } else {
        image_out->bpm = cpl_mask_duplicate(image1->bpm) ;
        cpl_mask_or(image_out->bpm, image2->bpm) ;
    }

    return image_out ;

#elif  CPL_OPERATION == CPL_IMAGE_BASIC_ASSIGN_LOCAL

    register float  *   pf1,
                    *   pf2 ;
    register int    *   pi1,
                    *   pi2 ;
    register double *   pd1,
                    *   pd2 ;
    register int        i ;

    cpl_ensure_code(im1, CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(im2, CPL_ERROR_NULL_INPUT);
    /* Input data images shall have the same sizes */
    cpl_ensure_code(im1->nx == im2->nx, CPL_ERROR_ILLEGAL_INPUT);
    cpl_ensure_code(im1->ny == im2->ny, CPL_ERROR_ILLEGAL_INPUT);

    assert( im1->pixels );
    assert( im2->pixels );

    /* Switch on the first passed image type */
    switch (im1->type) {
        case CPL_TYPE_INT:
            pi1 = (int*)im1->pixels ;
            /* Switch on the second passed image type */
            switch (im2->type) {
                case CPL_TYPE_INT:
                    pi2 = (int*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im1->ny) ; i++) 
                        CPL_OPERATOR(pi1[i], pi2[i]);
                    break ;
                case CPL_TYPE_FLOAT:
                    pf2 = (float*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im1->ny) ; i++) 
                        CPL_OPERATOR(pi1[i], pf2[i]);
                    break ;
                case CPL_TYPE_DOUBLE:
                    pd2 = (double*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im1->ny) ; i++) 
                        CPL_OPERATOR(pi1[i], pd2[i]);
                    break ;
                default:
                    cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH);
            }
            break ;
        case CPL_TYPE_FLOAT:
            pf1 = (float*)im1->pixels ;
            /* Switch on the second passed image type */
            switch (im2->type) {
                case CPL_TYPE_INT:
                    pi2 = (int*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im1->ny) ; i++) 
                        CPL_OPERATOR(pf1[i], pi2[i]);
                    break ;
                case CPL_TYPE_FLOAT:
                    pf2 = (float*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im1->ny) ; i++)
                        CPL_OPERATOR(pf1[i], pf2[i]);
                    break ;
                case CPL_TYPE_DOUBLE:
                    pd2 = (double*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im1->ny) ; i++) 
                        CPL_OPERATOR(pf1[i], pd2[i]);
                    break ;
                default:
                    cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH);
            }
            cpl_tools_add_flops( im1->nx * im1->ny );
            break ;
        case CPL_TYPE_DOUBLE:
            pd1 = (double*)im1->pixels ;
            /* Switch on the second passed image type */
            switch (im2->type) {
                case CPL_TYPE_INT:
                    pi2 = (int*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im1->ny) ; i++) 
                        CPL_OPERATOR(pd1[i], pi2[i]);
                    break ;
                case CPL_TYPE_FLOAT:
                    pf2 = (float*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im2->ny) ; i++) 
                        CPL_OPERATOR(pd1[i], pf2[i]);
                    break ;
                case CPL_TYPE_DOUBLE:
                    pd2 = (double*)im2->pixels ;
                    for (i=0 ; i<(im1->nx * im2->ny) ; i++)
                        CPL_OPERATOR(pd1[i], pd2[i]);
                    break ;
                default:
                    cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH);
            }
            cpl_tools_add_flops( im1->nx * im1->ny );
            break ;
        default:
            cpl_ensure_code(0, CPL_ERROR_TYPE_MISMATCH) ;
    }

    /* Handle bad pixels map */
    if (im2->bpm != NULL) {
        if (im1->bpm == NULL) {
            im1->bpm = cpl_mask_duplicate(im2->bpm) ;
        } else {
            cpl_mask_or(im1->bpm, im2->bpm) ;
        }
    }

    return CPL_ERROR_NONE;

#elif CPL_OPERATION == CPL_IMAGE_BASIC_OP_SCALAR
    
    case CPL_TYPE_T:
    {
        CPL_TYPE * pio = (CPL_TYPE*)image->pixels ;
        for (i=0 ; i < (image->nx * image->ny) ; i++)
            CPL_OPERATOR(pio[i], scalar);
        if (CPL_TYPE_T != CPL_TYPE_INT)
            cpl_tools_add_flops( image->nx * image->ny );
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_SQRT
    
    case CPL_TYPE_T:
    {
        CPL_TYPE * pio = (CPL_TYPE*)image->pixels ;
        for (i=0 ; i < (image->nx * image->ny) ; i++)
            pio[i] = sqrt(pio[i]);
        if (CPL_TYPE_T != CPL_TYPE_INT)
            cpl_tools_add_flops( image->nx * image->ny );
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_THRESHOLD

    case CPL_TYPE_T:
    {
        CPL_TYPE * pi;
        pi = (CPL_TYPE*)image_in->pixels ;
        for (i=0 ; i<(image_in->nx * image_in->ny) ; i++) {
            if (pi[i] > hi_cut) pi[i] = (CPL_TYPE)assign_hi_cut ;
            else if (pi[i] < lo_cut) pi[i] = (CPL_TYPE)assign_lo_cut ;
        }
        break ;
    } 

#elif CPL_OPERATION == CPL_IMAGE_BASIC_ABS

      case CPL_TYPE_T:
      {
          CPL_TYPE * pio = (CPL_TYPE*)image->pixels;
          for (i = 0; i < image->nx * image->ny ; i++)
              pio[i] = (CPL_TYPE)CPL_MATH_ABS((CPL_MATH_TYPE)pio[i]) ;
          if (CPL_TYPE_T != CPL_TYPE_INT)
              cpl_tools_add_flops( image->nx * image->ny );
          break ;
      }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_AVERAGE

    
    case CPL_TYPE_T:
    {
        CPL_TYPE * pi1;
        CPL_TYPE * po;
        image_out = cpl_image_new(image_1->nx, image_1->ny, CPL_TYPE_T);
        pi1 = (CPL_TYPE*)image_1->pixels ;
        po = (CPL_TYPE*)image_out->pixels ;
        /* Switch on second passed image type */
        switch (image_2->type) {
            case CPL_TYPE_INT:
                pii2 = (int*)image_2->pixels ;
                for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                    po[i] = (CPL_TYPE)(0.5 * (pi1[i] + pii2[i])) ;
                break ;
            case CPL_TYPE_FLOAT:
                pfi2 = (float*)image_2->pixels ;
                for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                    po[i] = (CPL_TYPE)(0.5 * (pi1[i] + pfi2[i])) ;
                break ;
            case CPL_TYPE_DOUBLE:
                pdi2 = (double*)image_2->pixels ;
                for (i=0 ; i<(image_out->nx * image_out->ny) ; i++)
                    po[i] = (CPL_TYPE)(0.5 * (pi1[i] + pdi2[i])) ;
                break ;
            default:
                cpl_image_delete(image_out) ;
                cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL) ;
        }
        if (CPL_TYPE_T != CPL_TYPE_INT)
            cpl_tools_add_flops( 2 * image_out->nx * image_out->ny );
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_EXTRACT

    case CPL_TYPE_T:
    {
        const int outlx = urx - llx + 1;
        const int outly = ury - lly + 1;

        /* Output pixel buffer */
        void    * po    = cpl_malloc(outlx * outly * sizeof(CPL_TYPE));

        const cpl_error_code error
            = cpl_tools_copy_window(po, in->pixels, sizeof(CPL_TYPE),
                                    in->nx, in->ny, llx, lly, urx, ury);

        if (error) {
            cpl_free(po);
        } else {
            self = CPL_TYPE_ADD(cpl_image_wrap)(outlx, outly, (CPL_TYPE*)po);
        }

        break;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_EXTRACTROW

    case CPL_TYPE_T:
    {
        CPL_TYPE * pi;
        pi = (CPL_TYPE*)image_in->pixels ; 
        for (i=0 ; i<image_in->nx ; i++) {
            out_data[i] = (double)pi[i+(pos-1)*image_in->nx] ;
        }
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_EXTRACTCOL

    case CPL_TYPE_T:
    {
        CPL_TYPE * pi;
        pi = (CPL_TYPE*)image_in->pixels ; 
        for (i=0 ; i<image_in->ny ; i++) {
            out_data[i] = (double)pi[pos-1+i*image_in->nx] ;
        }
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_COLLAPSE

    case CPL_TYPE_T:
    {
        const CPL_TYPE   * pi   = (CPL_TYPE*)inimage->pixels;
        const cpl_binary * bpmi = inimage->bpm == NULL ? NULL
            : cpl_mask_get_data(inimage->bpm);
        CPL_TYPE         * po;
        cpl_binary       * bpmo = NULL;
        double             sum;
        cpl_boolean        ok;
        int                i, j;


        if (direction == 0) {
            image1d = cpl_image_new(inimage->nx, 1, CPL_TYPE_T);
            po = (CPL_TYPE*)image1d->pixels;

            for (i=0; i < inimage->nx; i++) {
                sum = 0.0;
                ok = CPL_FALSE;
                for (j=0; j < inimage->ny; j++) {
                    if (bpmi == NULL || !bpmi[inimage->nx * j + i]) {
                        sum += pi[inimage->nx * j + i];
                        ok = CPL_TRUE;
                    }
                }
                po[i] = sum;
                if (!ok) {
                    if (bpmo == NULL)
                        bpmo = cpl_mask_get_data(cpl_image_get_bpm(image1d));
                    bpmo[i] = CPL_BINARY_1;
                }
            }
        } else if (direction == 1) {
            image1d = cpl_image_new(1, inimage->ny, CPL_TYPE_T);
            po = (CPL_TYPE*)image1d->pixels;
            for (j=0; j < inimage->ny; j++) {
                sum = 0.0;
                ok = CPL_FALSE;
                for (i=0; i < inimage->nx; i++) {
                    if (bpmi == NULL || !bpmi[inimage->nx * j + i]) {
                        sum += pi[inimage->nx * j + i];
                        ok = CPL_TRUE;
                    }
                }
                po[j] = sum;
                if (!ok) {
                    if (bpmo == NULL)
                        bpmo = cpl_mask_get_data(cpl_image_get_bpm(image1d));
                    bpmo[j] = CPL_BINARY_1;
                }
            }
        } else {
            cpl_ensure(0, CPL_ERROR_TYPE_MISMATCH, NULL);
        }
        cpl_tools_add_flops( image1d->nx * image1d->ny);

        break ;
    }
#elif CPL_OPERATION == CPL_IMAGE_BASIC_COLLAPSE_MEDIAN

    case CPL_TYPE_T:
    {
        CPL_TYPE * pi;
        CPL_TYPE * po;
        pi = (CPL_TYPE*)in->pixels ;
        if (direction == 1) {
            /* Collapsing the image in the x direction */
            collapsed = cpl_image_new(1, in->ny, CPL_TYPE_T);
            po = (CPL_TYPE*)collapsed->pixels ;
            width = in->nx - discard_lo - discard_hi ;
            line = cpl_malloc(width * sizeof(double)) ;
            for (j=0 ; j<in->ny ; j++) {
                for (i=0 ; i<width ; i++) {
                    line[i] = (double)pi[(i+discard_lo)+j*in->nx] ;
                }
                po[j] = (CPL_TYPE)cpl_tools_get_median_double(line, width);
            }
            cpl_free(line);
        } else if (direction==0) {
            /* Collapsing the image in the y direction */
            collapsed = cpl_image_new(in->nx, 1, CPL_TYPE_T);
            po = (CPL_TYPE*)collapsed->pixels ;
            width = in->ny - discard_lo - discard_hi ;
            line = cpl_malloc(width * sizeof(double)) ;
            for (i=0 ; i<in->nx ; i++) {
                for (j=0 ; j<width ; j++) {
                    line[j] = (double)pi[i+(j+discard_lo)*in->nx] ;
                }
                po[i] = (CPL_TYPE)cpl_tools_get_median_double(line, width);
            }
            cpl_free(line);
        } else {
            collapsed = NULL ;
        }
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_ROTATE_INT_LOCAL

    case CPL_TYPE_T:
    {
        /* rot is 0, 1, 2 or 3. */
        switch(rot) {
            case 1: {
                CPL_TYPE * pi = (CPL_TYPE *)self->pixels;
                int        i,j;

                if (self->nx == self->ny) {
                    /* If self->nx is even,
                       then there is a multiple of 4 pixels to move.
                       If self->nx is odd,
                       then there is also a multiple of 4 pixels to move,
                       since the center pixel does not move. */
                    /* The first four pixels to move are the corner ones,
                       followed by their neighbors
                       - the last four pixels to move are the center ones. */
                    for (j = 0; j < self->ny/2; j++) {
                        for (i = j; i < self->nx-1-j; i++) {
                            const CPL_TYPE tmp = pi[i + j * self->nx];

                            pi[i + j * self->nx]
                                = pi[(self->ny-1-j) + i * self->nx];

                            pi[(self->ny-1-j) + i * self->nx]
                                = pi[(self->nx-1-i) + (self->ny-1-j) * self->nx];

                            pi[(self->nx-1-i) + (self->ny-1-j) * self->nx]
                                = pi[j + (self->nx-1-i) * self->nx];

                            pi[j + (self->nx-1-i) * self->nx] = tmp;
                        }
                    }
                } else if (self->nx == 1) {
                    /* No pixels need to move - just swap nx and ny */
                    self->nx = self->ny;
                    self->ny = 1;
                } else {
                    /* Duplicate the input image :-( */
                    cpl_image      * tmp = cpl_image_duplicate(self);
                    const CPL_TYPE * pt = (const CPL_TYPE *)tmp->pixels;

                    self->nx = tmp->ny;
                    self->ny = tmp->nx;
                    pi += ((self->ny)-1)* (self->nx);
                    for (j=0; j<(self->nx); j++) {
                        for (i=0; i<(self->ny); i++) {
                            *pi = *pt++;
                            pi -= (self->nx);
                        }
                        pi += (self->nx)*(self->ny)+1;
                    }
                    cpl_image_delete(tmp);
                }
                break;
            }
            case 2: {
                CPL_TYPE * pi = (CPL_TYPE *)self->pixels;
                int        i;
                int        j = self->nx*self->ny-1;

                for (i = 0; i < j; i++, j--) {
                    const CPL_TYPE tmp = pi[i];
                    pi[i] = pi[j];
                    pi[j] = tmp;
                }
                break;
            }
            case 3: {
                CPL_TYPE * pi = (CPL_TYPE *)self->pixels;
                int        i,j;
                if (self->nx == self->ny) {
                    /* See case 1. */
                    for (j = 0; j < self->ny/2; j++) {
                        for (i = j; i < self->nx-1-j; i++) {
                            const CPL_TYPE tmp = pi[i + j * self->nx];

                            pi[i + j * self->nx]
                                = pi[j + (self->nx-1-i) * self->nx];

                            pi[j + (self->nx-1-i) * self->nx]
                                = pi[(self->nx-1-i) + (self->ny-1-j) * self->nx];

                            pi[(self->nx-1-i) + (self->ny-1-j) * self->nx]
                                = pi[(self->ny-1-j) + i * self->nx];

                            pi[(self->ny-1-j) + i * self->nx] = tmp;
                        }
                    }
                } else if (self->ny == 1) {
                    /* No pixels need to move - just swap nx and ny */
                    self->ny = self->nx;
                    self->nx = 1;
                } else {
                    /* Duplicate the input image :-( */
                    cpl_image      * tmp = cpl_image_duplicate(self);
                    const CPL_TYPE * pt = (const CPL_TYPE *)tmp->pixels;

                    self->nx = tmp->ny;
                    self->ny = tmp->nx;
                    pi += (self->nx)-1;
                    for (j=0; j<(self->nx); j++) {
                        for (i=0; i<(self->ny); i++) {
                            *pi = *pt++;
                            pi += (self->nx);
                        }
                        pi -= (self->nx)*(self->ny)+1;
                    }
                    cpl_image_delete(tmp);
                }
                break;
            }
            default:
                break;
        }
        break;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_SHIFT_INT_LOCAL

    case CPL_TYPE_T:
    {
        CPL_TYPE * pi;
        CPL_TYPE * pt;
        pt = (CPL_TYPE *)tmp_im->pixels ;
        pi = (CPL_TYPE *)im->pixels ;
        for (j=0 ; j<ny ; j++) {
            for (i=0 ; i<nx ; i++) {
                if ((i-x_shift>=0) && (i-x_shift<nx) && 
                        (j-y_shift>=0) && (j-y_shift<ny)) {
                    pi[i+j*nx] = pt[(i-x_shift)+(j-y_shift)*nx] ;
                } else {
                    pi[i+j*nx] = (CPL_TYPE)0.00 ;
                }
            }
        }
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_COPY

    case CPL_TYPE_T:
    {
    
        CPL_TYPE * p1 ;
        CPL_TYPE * p2 ;
        p1 = (CPL_TYPE *)im1->pixels ;
        p2 = (CPL_TYPE *)im2->pixels ;
        for (j=ypos-1 ; j<ury ; j++) {
            for (i=xpos-1 ; i<urx ; i++) {
                im1_pos = i+j*nx1 ;
                im2_pos = (i-xpos+1) + (j-ypos+1)*nx2 ;
                p1[im1_pos] = p2[im2_pos] ;
            }
        }
        break ;
}

#elif CPL_OPERATION == CPL_IMAGE_BASIC_FLIP_LOCAL

    case CPL_TYPE_T:
    {
        CPL_TYPE * pi;
        CPL_TYPE * pt;
        pt = (CPL_TYPE *)tmp_im->pixels ;
        pi = (CPL_TYPE *)im->pixels ;
        switch(angle) {
            case 0:
                pt += (ny-1)*nx ;
                for (j=0 ; j<ny ; j++) {
                    for (i=0 ; i<nx ; i++) *pi++ = *pt++ ;
                    pt -= 2*nx ;
                }
                break ;
            case 2:
                pt += (nx-1) ;
                for (j=0 ; j<ny ; j++) {
                    for (i=0 ; i<nx ; i++) *pi++ = *pt-- ;
                    pt += 2*nx ;
                }
                break ;
            case 1:
                im->nx = ny ;
                im->ny = nx ;
                for (j=0 ; j<nx ; j++) {
                    for (i=0 ; i<ny ; i++) {
                        *pi++ = *pt ;
                        pt += nx ;
                    }
                    pt -= (nx*ny-1) ;
                }
                break ;
            case 3:
                im->nx = ny ;
                im->ny = nx ;
                pt += (nx*ny-1) ;
                for (j=0 ; j<nx ; j++) {
                    for (i=0 ; i<ny ; i++) {
                        *pi++ = *pt ;
                        pt -= nx ;
                    }
                    pt += (nx*ny-1) ;
                }
                break ;
            default:
                cpl_image_delete(tmp_im) ;
                cpl_ensure_code(0, CPL_ERROR_ILLEGAL_INPUT) ;
                break ;
        }
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_BASIC_MOVE_PIXELS

    case CPL_TYPE_T:
    {
        /* Duplicate the input image */
        cpl_image * tmp_im = cpl_image_duplicate(im) ;
        CPL_TYPE  * po ;
        CPL_TYPE  * pi ;

        cpl_ensure_code(tmp_im, cpl_error_get_code());

        /* Get pointer to the data */
        pi = (CPL_TYPE *)tmp_im->pixels ;
        po = (CPL_TYPE *)im->pixels ;
                        
        /* Move the pixels */
        for (j=0 ; j<nb_cut ; j++) {
            for (i=0 ; i<nb_cut ; i++) {
                tile_x = (new_pos[i+j*nb_cut]-1) % nb_cut ;
                tile_y = (int)((new_pos[i+j*nb_cut]-1) / nb_cut) ;
                for (l=0 ; l<tile_sz_y ; l++) {
                    for (k=0 ; k<tile_sz_x ; k++) {
                        opos=(k+i*tile_sz_x) + im->nx*(l+j*tile_sz_y);
                        npos=(k+tile_x*tile_sz_x) + 
                            im->nx*(l+tile_y*tile_sz_y);
                        po[npos] = pi[opos] ;
                    }
                }
            }
        }
        cpl_image_delete(tmp_im);
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_CHK_POSITIVE
    
    case CPL_TYPE_T:
    {
        CPL_TYPE * pio = (CPL_TYPE*)image->pixels ;
        for (i=0 ; ok && i < (image->nx * image->ny) ; i++)
            if (pio[i] <= 0) ok = 0;
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_CHK_INTEGER

    case CPL_TYPE_T:
    {
        CPL_TYPE * pio = (CPL_TYPE*)image->pixels ;
        for (i=0 ; ok && i < (image->nx * image->ny) ; i++)
            if (pio[i] != ceil(pio[i])) ok = 0;
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_CHK_NON_ZERO

    case CPL_TYPE_T:
    {
        CPL_TYPE * pio = (CPL_TYPE*)image->pixels ;
        for (i=0 ; ok && i < (image->nx * image->ny) ; i++)
            if (pio[i] == 0) ok = 0;
        break ;
    }

#elif CPL_OPERATION == CPL_IMAGE_CHK_NON_NEGATIVE

    case CPL_TYPE_T:
    {
        CPL_TYPE * pio = (CPL_TYPE*)image->pixels ;
        for (i=0 ; ok && i < (image->nx * image->ny) ; i++)
            if (pio[i] < 0) ok = 0;
        break ;
    }

#else
#error "Undefined CPL Operation"
#endif

#undef CPL_TYPE
#undef CPL_TYPE_T
#undef CPL_MATH_ABS
#undef CPL_MATH_TYPE

#undef CPL_TYPE_ADD