/** @file  pmAstrometryUtils.c
 *
 *  @brief utility functions for transform and distort functions
 *
 *  @ingroup Astrometry
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2007-01-01 21:04:10 $
 *
 *  Copyright 2006 Institute for Astronomy, University of Hawaii
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <pslib.h>
#include "pmFPA.h"
#include "pmAstrometryUtils.h"

// given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y
// for which L,M = 0,0. tol is the allowed error on x,y.
psPlane *psPlaneTransformGetCenter (psPlaneTransform *trans, double tol)
{

    // crpix1,2 = X,Y(crval1,2)
    // start with linear solution for Xo,Yo:
    double R  = (trans->x->coeff[1][0]*trans->y->coeff[0][1] - trans->x->coeff[0][1]*trans->y->coeff[1][0]);
    double Xo = (trans->y->coeff[0][0]*trans->x->coeff[0][1] - trans->x->coeff[0][0]*trans->y->coeff[0][1])/R;
    double Yo = (trans->x->coeff[0][0]*trans->y->coeff[1][0] - trans->y->coeff[0][0]*trans->x->coeff[1][0])/R;

    // iterate to actual solution: requires small non-linear terms
    if (trans->x->nX > 1) {
        psPolynomial2D *XdX = psPolynomial2D_dX(NULL, trans->x);
        psPolynomial2D *XdY = psPolynomial2D_dY(NULL, trans->x);

        psPolynomial2D *YdX = psPolynomial2D_dX(NULL, trans->y);
        psPolynomial2D *YdY = psPolynomial2D_dY(NULL, trans->y);

        psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32);
        psVector *Beta = psVectorAlloc (2, PS_DATA_F32);

        /* this loop uses the Newton-Raphson method to solve for Xo,Yo
        * it needs the high order terms to be small 
        * Xo,Yo are in pixels;
        */
        double dPos = tol + 1;
        for (int i = 0; (dPos > tol) && (i < 20); i++) {
            // NOTE: order for Alpha is: [y][x]
            Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo);
            Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
            Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
            Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);

            Beta->data.F32[0] = psPolynomial2DEval (trans->x, Xo, Yo);
            Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo);

            psMatrixGJSolveF32 (Alpha, Beta);

            Xo -= Beta->data.F32[0];
            Yo -= Beta->data.F32[1];
            dPos = hypot(Beta->data.F32[0], Beta->data.F32[1]);
        }
        psFree (Alpha);
        psFree (Beta);
        psFree (XdX);
        psFree (XdY);
        psFree (YdX);
        psFree (YdY);
    }
    psPlane *center = psPlaneAlloc ();
    center->x = Xo;
    center->y = Yo;

    return center;
}

// convert a transformation L(x,y) to L'(x-xo,y-yo)
psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, double Xo, double Yo)
{

    // validate fit order

    if (output == NULL) {
        output = psPlaneTransformAlloc(input->x->nX, input->x->nY);
    }

    /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
     * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and 
     * noting that the constant term in each is the coefficient in the case of L(x,y) and is the 
     * value of L'(-xo,-yo) in the second case.
     */

    psPolynomial2D *tmp;

    psPolynomial2D *xPx = psPolynomial2DCopy (NULL, input->x);
    psPolynomial2D *yPx = psPolynomial2DCopy (NULL, input->y);

    for (int i = 0; i <= input->x->nX; i++) {
        psPolynomial2D *xPy = psPolynomial2DCopy (NULL, xPx);
        psPolynomial2D *yPy = psPolynomial2DCopy (NULL, yPx);
        for (int j = 0; j <= input->x->nY; j++) {
            output->x->mask[i][j] = input->x->mask[i][j];
            output->y->mask[i][j] = input->y->mask[i][j];
            output->x->coeff[i][j] = (output->x->mask[i][j]) ? 0 : psPolynomial2DEval (xPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);
            output->y->coeff[i][j] = (output->y->mask[i][j]) ? 0 : psPolynomial2DEval (yPy, Xo, Yo) / tgamma(i+1) / tgamma(j+1);

            // take the next derivative wrt y, catch output (is NULL on last pass)
            tmp = psPolynomial2D_dY(NULL, xPy);
            psFree (xPy);
            xPy = tmp;
            tmp = psPolynomial2D_dY(NULL, yPy);
            psFree (yPy);
            yPy = tmp;
        }
        // take the next derivative wrt x, catch output (is NULL on last pass)
        tmp = psPolynomial2D_dX(NULL, xPx);
        psFree (xPx);
        xPx = tmp;
        tmp = psPolynomial2D_dX(NULL, yPx);
        psFree (yPx);
        yPx = tmp;
    }
    return output;
}

// construct a psPlaneTransform which is the identify transformation
psPlaneTransform *psPlaneTransformIdentity (int order)
{

    psPlaneTransform *transform;

    if (order < 1)
        psAbort ("psastro", "invalid order");
    if (order > 3)
        psAbort ("psastro", "invalid order");

    // all coeffs and masks initially set to 0
    transform = psPlaneTransformAlloc (order, order);

    for (int i = 0; i <= order; i++) {
        for (int j = 0; j <= order; j++) {
            if (i + j > order) {
                transform->x->mask [i][j] = 1;
                transform->y->mask [i][j] = 1;
            }
        }
    }
    transform->x->coeff[1][0] = 1;
    transform->y->coeff[0][1] = 1;
    transform->x->mask[1][0] = 0;
    transform->y->mask[0][1] = 0;

    return transform;
}

// check that the given psPlaneTransform is the identity * (Xs,Ys)
bool psPlaneTransformIsDiagonal (psPlaneTransform *transform)
{

    int order;
    bool status;

    // we currently only support up to 3rd order polynomials
    if (transform->x->nX < 1)
        return false;
    if (transform->x->nY < 1)
        return false;
    if (transform->y->nX < 1)
        return false;
    if (transform->y->nY < 1)
        return false;

    if (transform->x->nX != transform->x->nY)
        return false;
    if (transform->y->nX != transform->y->nY)
        return false;

    // these are not actually valid tests
    if (transform->x->nX > 3)
        return false;
    if (transform->x->nY > 3)
        return false;
    if (transform->y->nX > 3)
        return false;
    if (transform->y->nY > 3)
        return false;

    status = true;
    order = transform->x->nX;
    for (int i = 0; i <= order; i++) {
        for (int j = 0; j <= order; j++) {
            if (i + j > order) {
                // high-order cross terms must be masked (eg, x^3 y^2)
                status &= transform->x->mask[i][j];
            } else {
                status &= !transform->x->mask[i][j];
                if ((i == 1) && (i + j == 1)) {
                    // linear, diagonal terms must be non-zero
                    status &= (fabs(transform->x->coeff[i][j]) > FLT_EPSILON);
                } else {
                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
                    status &= (fabs(transform->x->coeff[i][j]) < FLT_EPSILON);
                }
            }
        }
    }

    order = transform->y->nX;
    for (int i = 0; i <= order; i++) {
        for (int j = 0; j <= order; j++) {
            if (i + j > order) {
                // high-order cross terms must be masked (eg, x^3 y^2)
                status &= transform->y->mask[i][j];
            } else {
                status &= !transform->y->mask[i][j];
                if ((j == 1) && (i + j == 1)) {
                    // linear, diagonal terms must be 1.0
                    status &= (fabs(transform->y->coeff[i][j]) > FLT_EPSILON);
                } else {
                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
                    status &= (fabs(transform->y->coeff[i][j]) < FLT_EPSILON);
                }
            }
        }
    }
    return status;
}

// construct a psPlaneDistort which is the identify transformation
psPlaneDistort *psPlaneDistortIdentity (int order)
{

    psPlaneDistort *distort;

    if (order < 1)
        psAbort ("psastro", "invalid order");
    if (order > 3)
        psAbort ("psastro", "invalid order");

    // all coeffs and masks initially set to 0
    distort = psPlaneDistortAlloc (order, order, 0, 0);

    for (int i = 0; i <= order; i++) {
        for (int j = 0; j <= order; j++) {
            if (i + j > order) {
                distort->x->mask [i][j][0][0] = 1;
                distort->y->mask [i][j][0][0] = 1;
            }
        }
    }
    distort->x->coeff[1][0][0][0] = 1;
    distort->y->coeff[0][1][0][0] = 1;

    return distort;
}

// check that the given psPlaneDistort is the identity * (Xs,Ys)
bool psPlaneDistortIsDiagonal (psPlaneDistort *distort)
{

    int order;
    bool status;

    // we currently only support up to 3rd order polynomials
    if (distort->x->nX < 1)
        return false;
    if (distort->x->nY < 1)
        return false;
    if (distort->y->nX < 1)
        return false;
    if (distort->y->nY < 1)
        return false;

    if (distort->x->nX > 3)
        return false;
    if (distort->x->nY > 3)
        return false;
    if (distort->y->nX > 3)
        return false;
    if (distort->y->nY > 3)
        return false;

    if (distort->x->nZ > 0)
        return false;
    if (distort->x->nT > 0)
        return false;
    if (distort->y->nZ > 0)
        return false;
    if (distort->y->nT > 0)
        return false;

    if (distort->x->nX != distort->x->nY)
        return false;
    if (distort->y->nX != distort->y->nY)
        return false;

    status = true;
    order = distort->x->nX;
    for (int i = 0; i <= order; i++) {
        for (int j = 0; j <= order; j++) {
            if (i + j > order) {
                // high-order cross terms must be masked (eg, x^3 y^2)
                status &= distort->x->mask[i][j][0][0];
            } else {
                status &= !distort->x->mask[i][j][0][0];
                if ((i == 1) && (i + j == 1)) {
                    // linear, diagonal terms must be 1.0
                    status &= (fabs(distort->x->coeff[i][j][0][0]) > FLT_EPSILON);
                } else {
                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
                    status &= (fabs(distort->x->coeff[i][j][0][0]) < FLT_EPSILON);
                }
            }
        }
    }

    order = distort->y->nX;
    for (int i = 0; i <= order; i++) {
        for (int j = 0; j <= order; j++) {
            if (i + j > order) {
                // high-order cross terms must be masked (eg, x^3 y^2)
                status &= distort->y->mask[i][j][0][0];
            } else {
                status &= !distort->y->mask[i][j][0][0];
                if ((j == 1) && (i + j == 1)) {
                    // linear, diagonal terms must be 1.0
                    status &= (fabs(distort->y->coeff[i][j][0][0]) > FLT_EPSILON);
                } else {
                    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
                    status &= (fabs(distort->y->coeff[i][j][0][0]) < FLT_EPSILON);
                }
            }
        }
    }
    return status;
}
