/* $Id: cpl_geom_img_body.h,v 1.2 2007/07/06 13:57:49 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
#define CPL_IMAGE_GET_DATA cpl_image_get_data_double
#define CPL_IMAGE_GET_DATA_CONST cpl_image_get_data_double_const

#elif CPL_CLASS == CPL_CLASS_FLOAT
#define CPL_TYPE float
#define CPL_TYPE_T CPL_TYPE_FLOAT
#define CPL_IMAGE_GET_DATA cpl_image_get_data_float
#define CPL_IMAGE_GET_DATA_CONST cpl_image_get_data_float_const

#else
#undef CPL_TYPE
#undef CPL_TYPE_T
#undef CPL_IMAGE_GET_DATA
#endif

#if CPL_OPERATION == CPL_GEOM_IMA_OFFSET_SAA

    case CPL_TYPE_T:
    {
        const CPL_TYPE *   pi ;
        CPL_TYPE       *   po ;

        /* Create output image */
        final = cpl_image_new(nx, ny, CPL_TYPE_T) ;
        po = CPL_IMAGE_GET_DATA(final) ;
        
        /* Triple loop */
        for (j=0 ; j<ny ; j++) {
            cpl_msg_progress(cpl_func, j, ny, "Shift and add") ;
            for (i=0 ; i<nx ; i++) {
                ncontrib = 0 ;
                for (p=0 ; p<nima ; p++) {
                    x = (double)i - offs_x[p] + start_x ;
                    y = (double)j - offs_y[p] + start_y ;

                    px = (int)x ;
                    py = (int)y ;
        
                    /* If original pixel is present in the current plane */
                    if (px>1 && px<(sizex-2) && py>1 && py<(sizey-2)) {
                        /* Interpolate from input image */
                        pos = px + py * sizex ;
                        /* Avoid overhead of function call here */
                        pi = CPL_IMAGE_GET_DATA_CONST(cpl_imagelist_get_const(ilist, p));
                        for (k=0 ; k<16 ; k++) {
                            neighbors[k] = (double)pi[pos+leaps[k]] ;
                        }
                        /* Which tabulated value index shall be used? */
                        tabx = (int)(0.5+(x-(double)px)*(double)(tabsperpix)) ;
                        taby = (int)(0.5+(y-(double)py)*(double)(tabsperpix)) ;
        
                        /* Compute resampling coefficients */
                        /* rsc[0..3] in x, rsc[4..7] in y  */

                        rsc[0] = interp_kernel[tabsperpix + tabx] ;
                        rsc[1] = interp_kernel[tabx] ;
                        rsc[2] = interp_kernel[tabsperpix - tabx] ;
                        rsc[3] = interp_kernel[2*tabsperpix - tabx] ;
                        rsc[4] = interp_kernel[tabsperpix + taby] ;
                        rsc[5] = interp_kernel[taby] ;
                        rsc[6] = interp_kernel[tabsperpix - taby] ;
                        rsc[7] = interp_kernel[2*tabsperpix - taby] ;

                        sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
                                (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;

                        /* Compute interpolated pixel now   */
                        interp= rsc[4] * (  rsc[0]*neighbors[0] +
                                            rsc[1]*neighbors[1] +
                                            rsc[2]*neighbors[2] +
                                            rsc[3]*neighbors[3] ) +
                                rsc[5] * (  rsc[0]*neighbors[4] +
                                            rsc[1]*neighbors[5] +
                                            rsc[2]*neighbors[6] +
                                            rsc[3]*neighbors[7] ) +
                                rsc[6] * (  rsc[0]*neighbors[8] +
                                            rsc[1]*neighbors[9] +
                                            rsc[2]*neighbors[10] +
                                            rsc[3]*neighbors[11] ) +
                                rsc[7] * (  rsc[0]*neighbors[12] +
                                            rsc[1]*neighbors[13] +
                                            rsc[2]*neighbors[14] +
                                            rsc[3]*neighbors[15] ) ;
                        interp /= sumrs ;
                        acc[ncontrib] = interp ;
                        ncontrib ++ ;
                    }
                }
                finpix = 0 ;
                if ((ncontrib>0) && (ncontrib>rtot)) {
                    cpl_tools_sort_double(acc, ncontrib) ;
                    for (p=rmin ; p<(ncontrib-rmax) ; p++) 
                        finpix += (double)acc[p];
                    finpix /= (double)(ncontrib-rtot) ;
                }
                po[i+j*nx] = (CPL_TYPE)finpix ;
                pcontrib[i+j*nx] = (int)(ncontrib-rtot) ;
            }
        }
        break ;
    }

#elif CPL_OPERATION == CPL_GEOM_IMA_OFFSET_XCORR

    case CPL_TYPE_T:
    {
        const CPL_TYPE   *   pi1 ;
        const CPL_TYPE   *   pi2 ;
        const CPL_TYPE   *   reg1 ;
        const CPL_TYPE   *   reg2 ;

        /* Move into the buffers, up to the requested searching place   */
        pi1 = CPL_IMAGE_GET_DATA_CONST(im1) ;
        pi2 = CPL_IMAGE_GET_DATA_CONST(im2) ;
        pi1 += im1_posx + im1_posy * nx1 ;
        pi2 += im2_posx + im2_posy * nx2 ;
        
        inc1 = nx1 - (2*m_hx+1) ;
        inc2 = nx2 - (2*m_hx+1) ;

        /* Compute the cross-correlation factors and keep the minimum one */
        for (l=-s_hy ; l<=s_hy ; l++) {
            for (k=-s_hx ; k<=s_hx ; k++) {
                sum = 0 ;
                reg1 = pi1 + k - m_hx + (l-m_hy) * nx1 ;
                reg2 = pi2 - m_hx - m_hy * nx2 ;
                /* One cross-correlation measure */
                for (j=-m_hy ; j<=m_hy ; j++) {
                    for (i=-m_hx ; i<=m_hx ; i++) {
                        value = (double)(*reg1++)-(double)(*reg2++) ;
                        value *= value ;
                        sum += value ;
                    }
                    reg1 += inc1 ;
                    reg2 += inc2 ;
                }
                correl[s_hx+k+(2*s_hx+1)*(s_hy+l)] = sum * inv_surface ;
                /* Keep it if it is the minimum */
                if (sum < sum_min) {
                    l_min = l ;
                    k_min = k ;
                    sum_min = sum ;
                }
            }
        }
        break ;
    }

#endif

#undef CPL_TYPE
#undef CPL_TYPE_T
#undef CPL_IMAGE_GET_DATA
#undef CPL_IMAGE_GET_DATA_CONST