/** @file  pmTrend2D.c
 *
 *  @author EAM, IfA
 *
 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2008-01-15 02:47:51 $
 *  Copyright 2004 Institute for Astronomy, University of Hawaii
 *
 */

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

# include <strings.h>
# include <pslib.h>
# include "pmTrend2D.h"

static void pmTrend2DFree (pmTrend2D *trend) {

    if (trend == NULL) 
        return;

    psFree (trend->stats);
    psFree (trend->poly);
    psFree (trend->map);
    return;
}

pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats)
{
    PS_ASSERT_PTR_NON_NULL(stats, NULL);
    if (mode == PM_TREND_MAP) {
        PS_ASSERT_PTR_NON_NULL(image, NULL);
    }

    pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
    psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);

    trend->map = NULL;
    trend->poly = NULL;
    trend->stats = psMemIncrRefCounter (stats);
    trend->mode = mode;
	
    switch (mode) {
      case PM_TREND_POLY_ORD:
	trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
	// set masking somehow
 	for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
	    for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
		if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
		    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
		} else {
		    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
		}
	    }
	}
	break;

      case PM_TREND_POLY_CHEB:
	trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, nXtrend, nYtrend);
	break;

      case PM_TREND_MAP: {
	  // binning defines the map scale relationship
	  psImageBinning *binning = psImageBinningAlloc();
	  binning->nXruff = nXtrend;
	  binning->nYruff = nYtrend;
	  binning->nXfine = image->numCols;
	  binning->nYfine = image->numRows;

	  trend->map = psImageMapAlloc (image, binning, stats);
	  psFree (binning);
	  break;
      }
      // XXX: Put a more graceful error here.
      default:
	psAbort ("error");
    }
    return (trend);
}

bool psMemCheckTrend2D(psPtr ptr)
{
    PS_ASSERT_PTR(ptr, false);
    return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmTrend2DFree);
}

pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats)
{
    if (mode == PM_TREND_MAP) {
        PS_ASSERT_PTR_NON_NULL(binning, NULL);
        PS_ASSERT_PTR_NON_NULL(stats, NULL);
    }
    pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
    psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);

    trend->map = NULL;
    trend->poly = NULL;
    trend->stats = psMemIncrRefCounter (stats);
    trend->mode = mode;
	
    switch (mode) {
      case PM_TREND_POLY_ORD:
	trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, binning->nXruff, binning->nYruff);
	// set masking somehow
 	for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
	    for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
		if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
		    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
		} else {
		    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
		}
	    }
	}
	break;

      case PM_TREND_POLY_CHEB:
	trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, binning->nXruff, binning->nYruff);
	break;

      case PM_TREND_MAP: {
	  // binning defines the map scale relationship
	  trend->map = psImageMapNoImageAlloc (binning, stats);
	  break;
      }

      default:
	psAbort ("error");
    }
    return (trend);
}

pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats)
{
    PS_ASSERT_PTR_NON_NULL(stats, NULL);
    pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
    psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);

    trend->map = NULL;
    trend->poly = NULL;
    trend->stats = psMemIncrRefCounter (stats);
    trend->mode = mode;
	
    switch (mode) {
      case PM_TREND_POLY_ORD:
	trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
	// set masking somehow
 	for (int nx = 0; nx < trend->poly->nX + 1; nx++) {
	    for (int ny = 0; ny < trend->poly->nY + 1; ny++) {
		if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) {
		    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_SET;
		} else {
		    trend->poly->coeffMask[nx][ny] = PS_POLY_MASK_NONE;
		}
	    }
	}
	break;

      case PM_TREND_POLY_CHEB:
	trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, nXtrend, nYtrend);
	break;

      case PM_TREND_MAP: {
	  // binning defines the map scale relationship
	  psImageBinning *binning = psImageBinningAlloc();
	  binning->nXfine = nXfield;
	  binning->nYfine = nYfield;
	  binning->nXruff = nXtrend;
	  binning->nYruff = nYtrend;

	  trend->map = psImageMapAlloc (NULL, binning, stats);
	  psFree (binning);
	  break;
      }

      default:
        // XXX: Put a more graceful error here.
	psAbort ("error");
    }
    return (trend);
}

bool pmTrend2DFit (pmTrend2D *trend, psVector *mask, psMaskType maskVal, psVector *x,
                   psVector *y, psVector *f, psVector *df)
{
    PS_ASSERT_PTR_NON_NULL(trend, false);
    PS_ASSERT_VECTOR_NON_NULL(x, false);
    PS_ASSERT_VECTOR_NON_NULL(y, false);
    PS_ASSERT_VECTOR_NON_NULL(f, false);

    bool status;
    switch (trend->mode) {
      case PM_TREND_POLY_ORD:
      case PM_TREND_POLY_CHEB:
        status = psVectorClipFitPolynomial2D (trend->poly, trend->stats, mask, maskVal, f, df, x, y);
	// we can use the API here which adjusts the polynomial order based on the number
	// of points in the image, and potentially based on the fractional range of the
	// data?
	break;

      case PM_TREND_MAP:
	// XXX supply fraction from trend elements
	// XXX need to add the API which adjusts the scale
	status = psImageMapClipFit (trend->map, trend->stats, mask, maskVal, x, y, f, df);
	break;

      default:
	psAbort ("error");
    }
    return status;
}

double pmTrend2DEval (pmTrend2D *trend, float x, float y)
{
    if (!trend) return 0.0;

    double result;
    switch (trend->mode) {
      case PM_TREND_POLY_ORD:
      case PM_TREND_POLY_CHEB:
	result = psPolynomial2DEval (trend->poly, x, y);
	break;

      case PM_TREND_MAP:
	result = psImageMapEval (trend->map, x, y);
	break;

      default:
	psAbort ("error");
    }
    return result;
}

psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y)
{
    PS_ASSERT_PTR_NON_NULL(trend, NULL);
    PS_ASSERT_VECTOR_NON_NULL(x, false);
    PS_ASSERT_VECTOR_NON_NULL(y, false);
    psVector *result;

    switch (trend->mode) {
      case PM_TREND_POLY_ORD:
      case PM_TREND_POLY_CHEB:
	result = psPolynomial2DEvalVector (trend->poly, x, y);
	break;

      case PM_TREND_MAP:
	result = psImageMapEvalVector (trend->map, x, y);
	break;

      default:
	psAbort ("error");
    }
    return result;
}

psString pmTrend2DModeToString (pmTrend2DMode mode) {
    
    psString name;

    switch (mode) {
      case PM_TREND_NONE:
	name = psStringCopy ("NONE");
	break;
      case PM_TREND_POLY_ORD:
	name = psStringCopy ("POLY_ORD");
	break;
      case PM_TREND_POLY_CHEB:
	name = psStringCopy ("POLY_CHEB");
	break;
      case PM_TREND_MAP:
	name = psStringCopy ("MAP");
	break;
      default:
        psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode\n");
	psAbort ("invalid mode %d", mode);
    }
    return name;
}

pmTrend2DMode pmTrend2DModeFromString (psString name) {

    if (!name) return PM_TREND_NONE;

    if (!strcasecmp (name, "NONE")) {
	return PM_TREND_NONE;
    }
    if (!strcasecmp (name, "POLY_ORD")) {
	return PM_TREND_POLY_ORD;
    }
    if (!strcasecmp (name, "POLY_CHEB")) {
	return PM_TREND_POLY_CHEB;
    }
    if (!strcasecmp (name, "MAP")) {
	return PM_TREND_MAP;
    }
    psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode %s\n", name);
    return PM_TREND_NONE;
}
