/*****************************************************************************/

// Copyright 2006-2007 Adobe Systems Incorporated

// All Rights Reserved.

//

// NOTICE:  Adobe permits you to use, modify, and distribute this file in

// accordance with the terms of the Adobe license agreement accompanying it.

/*****************************************************************************/

/* $Id: //mondo/dng_sdk_1_4/dng_sdk/source/dng_spline.cpp#1 $ */ 

/* $DateTime: 2012/05/30 13:28:51 $ */

/* $Change: 832332 $ */

/* $Author: tknoll $ */

/*****************************************************************************/

#include "dng_spline.h"

#include "dng_assertions.h"

#include "dng_exceptions.h"

/******************************************************************************/

dng_spline_solver::dng_spline_solver ()

  : X ()

  , Y ()

  , S ()

  {

  }

/******************************************************************************/

dng_spline_solver::~dng_spline_solver ()

  {

  }

/******************************************************************************/

void dng_spline_solver::Reset ()

  {

  X.clear ();

  Y.clear ();

  S.clear ();

  }

/******************************************************************************/

void dng_spline_solver::Add (real64 x, real64 y)

  {

  X.push_back (x);

  Y.push_back (y);

  }

/******************************************************************************/

void dng_spline_solver::Solve ()

  {

  // This code computes the unique curve such that:

  //    It is C0, C1, and C2 continuous

  //    The second derivative is zero at the end points

  int32 count = (int32) X.size ();

  DNG_ASSERT (count >= 2, "Too few points");

  int32 start = 0;

  int32 end   = count;

  real64 A =  X [start+1] - X [start];

  real64 B = (Y [start+1] - Y [start]) / A;

  S.resize (count);

  S [start] = B;

  int32 j;

  // Slopes here are a weighted average of the slopes

  // to each of the adjcent control points.

  for (j = start + 2; j < end; ++j)

    {

    real64 C = X [j] - X [j-1];

    real64 D = (Y [j] - Y [j-1]) / C;

    S [j-1] = (B * C + D * A) / (A + C);

    A = C;

    B = D;

    }

  S [end-1] = 2.0 * B - S [end-2];

  S [start] = 2.0 * S [start] - S [start+1];

  if ((end - start) > 2)

    {

    dng_std_vector<real64> E;

    dng_std_vector<real64> F;

    dng_std_vector<real64> G;

    E.resize (count);

    F.resize (count);

    G.resize (count);

    F [start] = 0.5;

    E [end-1] = 0.5;

    G [start] = 0.75 * (S [start] + S [start+1]);

    G [end-1] = 0.75 * (S [end-2] + S [end-1]);

    for (j = start+1; j < end - 1; ++j)

      {

      A = (X [j+1] - X [j-1]) * 2.0;

      E [j] = (X [j+1] - X [j]) / A;

      F [j] = (X [j] - X [j-1]) / A;

      G [j] = 1.5 * S [j];

      }

    for (j = start+1; j < end; ++j)

      {

      A = 1.0 - F [j-1] * E [j];

      if (j != end-1) F [j] /= A;

      G [j] = (G [j] - G [j-1] * E [j]) / A;

      }

    for (j = end - 2; j >= start; --j)

      G [j] = G [j] - F [j] * G [j+1];

    for (j = start; j < end; ++j)

      S [j] = G [j];

    }

  }

/******************************************************************************/

bool dng_spline_solver::IsIdentity () const

  {

  int32 count = (int32) X.size ();

  if (count != 2)

    return false;

  if (X [0] != 0.0 || X [1] != 1.0 ||

    Y [0] != 0.0 || Y [1] != 1.0)

    return false;

  return true;

  }

/******************************************************************************/

real64 dng_spline_solver::Evaluate (real64 x) const

  {

  int32 count = (int32) X.size ();

  // Check for off each end of point list.

  if (x <= X [0])

    return Y [0];

  if (x >= X [count-1])

    return Y [count-1];

  // Binary search for the index.

  int32 lower = 1;

  int32 upper = count - 1;

  while (upper > lower)

    {

    int32 mid = (lower + upper) >> 1;

    if (x == X [mid])

      {

      return Y [mid];

      }

    if (x > X [mid])

      lower = mid + 1;

    else

      upper = mid;

    }

  DNG_ASSERT (upper == lower, "Binary search error in point list");

  int32 j = lower;

  // X [j - 1] < x <= X [j]

  // A is the distance between the X [j] and X [j - 1]

  // B and C describe the fractional distance to either side. B + C = 1.

  // We compute a cubic spline between the two points with slopes

  // S[j-1] and S[j] at either end. Specifically, we compute the 1-D Bezier

  // with control values:

  //

  //    Y[j-1], Y[j-1] + S[j-1]*A, Y[j]-S[j]*A, Y[j]

  return EvaluateSplineSegment (x,

                  X [j - 1],

                  Y [j - 1],

                  S [j - 1],

                  X [j    ],

                  Y [j    ],

                  S [j    ]);

  }

/*****************************************************************************/