/*
   For more information, please see: http://software.sci.utah.edu

   The MIT License

   Copyright (c) 2009 Scientific Computing and Imaging Institute,
   University of Utah.

   License for the specific language governing rights and limitations under
   Permission is hereby granted, free of charge, to any person obtaining a
   copy of this software and associated documentation files (the "Software"),
   to deal in the Software without restriction, including without limitation
   the rights to use, copy, modify, merge, publish, distribute, sublicense,
   and/or sell copies of the Software, and to permit persons to whom the
   Software is furnished to do so, subject to the following conditions:

   The above copyright notice and this permission notice shall be included
   in all copies or substantial portions of the Software.

   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
   OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
   DEALINGS IN THE SOFTWARE.
*/


/*
 *  FieldlineAnalyzerLib.cc:
 *
 *  Written by:
 *   Allen Sanderson
 *   SCI Institute
 *   University of Utah
 *   September 2006
 *
 */

#include <FieldlineAnalyzerLib.h>

#include <algorithm>
#include <iostream>
#include <map>

#include <float.h>

#ifdef STRAIGHTLINE_SKELETON
#include "dir_graph.h"
#endif

#define SIGN(x) ((x) < 0.0 ? (int) -1 : (int) 1)


#ifdef COMMENT_OUT
/////////// Begin Guoning Code


/*
   This routine computes the mean and standard deviation
        for a sequence of winding pairs.
        Added by Guoning 10/05/2010
*/
void
standard_deviation_analysis( vector< int > t_windings,
                             vector< int > p_windings,
                             int truncate_size,
                             double &mean_saftyFactor,
                             double &max_saftyFactor,
                             double &min_saftyFactor,
                             vector< double > &deviations,
                             double &std_deviation)
{
  vector < double > vals;
  unsigned int i;
        
  /*cerr << " numbers of values in t_windings and p_windings are (" <<
    t_windings.size() << ", " << p_windings.size() << ")" << std::endl;
    std::cerr << " The input pairs are " << std::endl;

    for (i=0; i<t_windings.size(); i++)
    {
    std::cerr << "(" << t_windings[i] << ", "
    << p_windings[i] << ")" << std::endl;
    }
  */

  // we first compute the estimated slope (i.e. safety factor) for each pair
  for (i=0; i<t_windings.size(); i++)
  {
    if (p_windings[i] == 0)
      vals.push_back(0.);
    else
      vals.push_back((double)t_windings[i]/(double)p_windings[i]);
  }

  // we need to sort the obtaining safety factors in 'vals' and remove the front few ones
  
  //cerr << "start sorting the vals ..." << std::endl;
  double *temp_vals_array = new double[vals.size()];
  for (i=0; i<vals.size(); i++)
    temp_vals_array[i]=vals[i];
  quickSort <double> (temp_vals_array, 0, vals.size()-1);
  //cerr << "finish sorting. " << std::endl;
  for (i=0; i<vals.size(); i++)
    vals[i]=temp_vals_array[i];
        
  delete [] temp_vals_array;
        
  // Here we may remove the max and min or multiple max and min_dist
  //int truncate_size = 2;
  for (i=0; i<truncate_size; i++)
  {
    if (vals.size()>2)
      vals.erase(vals.begin());
  }
  
  for (i=0; i<truncate_size; i++)
  {
    if (vals.size()>2)
      vals.erase(vals.begin()+vals.size()-1);
  }
  
  // search for the max and min slope for the later envelope construction
  min_saftyFactor = 10000.;
  max_saftyFactor = 0.;
  for (i=0; i<vals.size(); i++)
  {
    if (vals[i]>max_saftyFactor) max_saftyFactor = vals[i];
    if (vals[i]<min_saftyFactor) min_saftyFactor = vals[i];
  }

  // second, we compute the mean
  mean_saftyFactor = 0.;
  for (i=0; i<vals.size(); i++)
    mean_saftyFactor += vals[i];

  mean_saftyFactor /= vals.size();
  //cerr << "current mean safetyFactor is " << mean_saftyFactor << std::endl;

  // third, we compute the deviation of each estimated slope to the mean
  deviations.clear();
  for (i=0; i<vals.size(); i++)
    deviations.push_back(vals[i]-mean_saftyFactor);
  //cerr << "current number of values in vals is " << vals.size() << std::endl;
                
  ////////////////////////////////////////////////////////////////////
  // additional step, we remove the maximal positive and negative deviation

/*      double max_positive, max_negative;
        max_positive = -10000.;
        max_negative = 10000.;
        for (i=0; i<vals.size(); i++)
        {
                if (max_positive<vals[i]) max_positive=vals[i];

                if (max_negative>vals[i]) max_negative=vals[i];
        }
        
        for (i=0; i<vals.size(); i++)
        {
                if (vals[i]>=max_positive) vals.erase(vals.begin()+i);

                if (vals[i]<=max_negative) vals.erase(vals.begin()+i);
        }

        min_saftyFactor = 10000.;
        max_saftyFactor = 0.;
        for (i=0; i<vals.size(); i++)
        {
                if (vals[i]>max_saftyFactor) max_saftyFactor = vals[i];
                if (vals[i]<min_saftyFactor) min_saftyFactor = vals[i];
        }

        // second, we compute the mean
        mean_saftyFactor = 0.;
        for (i=0; i<vals.size(); i++)
                mean_saftyFactor += vals[i];

        mean_saftyFactor /= vals.size();
        //cerr << "current mean safetyFactor is " << mean_saftyFactor << std::endl;
        
        // third, we compute the deviation of each estimated slope to the mean
        deviations.clear();
        for (i=0; i<vals.size(); i++)
                deviations.push_back(vals[i]-mean_saftyFactor);
        
        //cerr << "current number of values in vals is " << vals.size() << std::endl;
*/

  ////////////////////////////////////////////////////////////////////
  // compute the square of each deviation
  for (i=0; i<deviations.size(); i++)
    deviations[i] *= deviations[i];

  // fourth, we compute the standard deviation by summing all the
  // above deviations and taking the square root of the sum
  std_deviation = 0.;
  for (i=0; i<deviations.size(); i++)
    std_deviation += deviations[i];
  
  std_deviation = sqrt(std_deviation);
}


/*
  This routine is used to removed redundant poloidal entries
*/
void
windingPairs_cleanUp(vector <int> &toroidal, vector <int> &poloidal)
{
  unsigned int i;
  
  // we first remove all the zero poloidal winding
  for (i=0; i<poloidal.size(); i++)
  {
    if (poloidal[0]==0)
    {
      poloidal.erase(poloidal.begin());
      toroidal.erase(toroidal.begin());
    }
    else
      break;
  }
  
  int cur_poloidal = poloidal[0];
  
  for (i=1; i<toroidal.size(); i++)
  {
    if (poloidal[i]==cur_poloidal)
    {
      poloidal.erase(poloidal.begin()+i);
      toroidal.erase(toroidal.begin()+i);
      i--;
    }
    
    else
      cur_poloidal = poloidal[i];
  }
}

/*
  For better estimation of safety factor, we remove the first half of
  the pairs.  Note that this routine should be called before calling
  the windingPairs_cleanUp to remove the redundant pairs (with the
  same poloidal windings, i.e. flat).
*/
void
remove_first_half_pairs(vector <int> &toroidal, vector <int> &poloidal)
{
  vector <int> half_tor;
  vector <int> half_pol;
  int half_pos = toroidal.size()/2;
  unsigned int i;

  // Take the second half
  for (i=half_pos; i<toroidal.size(); i++)
  {
    half_tor.push_back (toroidal[i]);
    half_pol.push_back (poloidal[i]);
  }

  // Store the second half to the original arrays
  toroidal.clear();
  poloidal.clear();
  for (i=0; i<half_tor.size(); i++)
  {
    toroidal.push_back (half_tor[i]);
    poloidal.push_back (half_pol[i]);
  }

  half_tor.clear();
  half_pol.clear();
}

#endif

// A simple linear regression using linear lease square fitting
template< class TYPE >
void least_square_fit( std::vector< TYPE > p_windings,
                       double &bestfit_safetyFactor)
{
  unsigned int start = (unsigned int)p_windings.size() *.75;
  unsigned int stop  = (unsigned int)p_windings.size();

  unsigned int cc = 0;

  double t_mean = 0;
  for (unsigned int i=start; i<stop; i++, ++cc)
    t_mean += i;

  t_mean /= (double) cc;

  double t_square_sum = 0;
  for (unsigned int i=start; i<stop; i++)
    t_square_sum += (double) (i * i);

  double SS_xx = t_square_sum-cc*t_mean*t_mean;

  double p_mean = 0;
  for (unsigned int i=start; i<stop; i++)
    p_mean += p_windings[i];

  p_mean /= (double) cc;

  double p_square_sum = 0;
  for (unsigned int i=start; i<stop; i++)
    p_square_sum += p_windings[i]*p_windings[i];

  double SS_xy = 0;

  for (unsigned int i=start; i<stop; i++)
    SS_xy += (i-t_mean)*(p_windings[i]-p_mean);

  bestfit_safetyFactor = SS_xx/SS_xy;
}


/////////// End Guoning Code







Point FieldlineLib::interpert( Point lastPt, Point currPt, double t ) {

  return Point( Vector( lastPt ) + Vector( currPt - lastPt ) * t );
}


// Is point 0 (v0) clockwise or counter clockwise of point 1 (v1)
int FieldlineLib::ccwXY( Vector v0, Vector v1 ) {
    
  if( v0.x * v1.y - v0.y * v1.x > FLT_MIN ) return  1;    // CCW
  if( v0.y * v1.x - v0.x * v1.y > FLT_MIN ) return -1;    // CW
  if( v0.x * v1.x < 0.0 || v0.y * v1.y < 0.0 ) return -1; // CW
    
  if( v0.x*v0.x+v0.y*v0.y >=
      v1.x*v1.x+v1.y*v1.y ) return 0;                     //  ON LINE
    
  return 1;                                               //  CCW
}


// Is point 0 (v0) clockwise or counter clockwise of point 1 (v1)
int FieldlineLib::ccwXZ( Vector v0, Vector v1 ) {
    
  if( v0.x * v1.z - v0.z * v1.x > FLT_MIN ) return  1;    // CCW
  if( v0.z * v1.x - v0.x * v1.z > FLT_MIN ) return -1;    // CW
  if( v0.x * v1.x < 0.0 || v0.z * v1.z < 0.0 ) return -1; // CW
    
  if( v0.x*v0.x+v0.z*v0.z >=
      v1.x*v1.x+v1.z*v1.z ) return 0;                     //  ON LINE
    
  return 1;                                               //  CCW
}


int FieldlineLib::intersect( Point l0_p0, Point l0_p1,
                             Point l1_p0, Point l1_p1 )
{
  if( (l0_p0.x < l1_p0.x && l0_p0.x < l1_p1.x && 
       l0_p1.x < l1_p0.x && l0_p1.x < l1_p1.x) ||

      (l0_p0.x > l1_p0.x && l0_p0.x > l1_p1.x && 
       l0_p1.x > l1_p0.x && l0_p1.x > l1_p1.x) ||

      (l0_p0.z < l1_p0.z && l0_p0.z < l1_p1.z && 
       l0_p1.z < l1_p0.z && l0_p1.z < l1_p1.z) ||

      (l0_p0.z > l1_p0.z && l0_p0.z > l1_p1.z && 
       l0_p1.z > l1_p0.z && l0_p1.z > l1_p1.z) )
  //  Lines do not intersect no bounding box intersection.
  return 0;

  //  See if the lines intersect.    
  if( ( ccwXZ( Vector(l0_p1-l0_p0), Vector(l1_p0-l0_p0)) *
        ccwXZ( Vector(l0_p1-l0_p0), Vector(l1_p1-l0_p0) ) <= 0 ) &&
      ( ccwXZ( Vector(l1_p1-l1_p0), Vector(l0_p0-l1_p0)) *
        ccwXZ( Vector(l1_p1-l1_p0), Vector(l0_p1-l1_p0) ) <= 0 ) ) {
        
    //  See if there is not a shared point.
    if( l0_p0 != l1_p0 && l0_p0 != l1_p1 &&
        l0_p1 != l1_p0 && l0_p1 != l1_p1 )
      return 1;
        
    //  See if there is a shared point.
    else if( l0_p0 == l1_p0 || l0_p0 == l1_p1 ||
             l0_p1 == l1_p0 || l0_p1 == l1_p1 )
      return 2;
        
    //  There must be a point that is on the other line.
    else
      return 3;
  }
    
  //  Lines do not intersect.
  return 0;
}


void
FieldlineLib::convexHull( std::vector< std::pair< Point, unsigned int > > &hullPts,
                          unsigned int &m, // starting index
                          unsigned int npts, // number of points
                          int dir )
{
  // Three or less points so which is a hull
  if( npts - m < 3 ) {
    m = npts;
    return;
  }

  unsigned int min = m;

  // Find the point with the minimum z value.
  for( unsigned int i=m; i<npts; i++ ) {
    if( hullPts[ min ].first.z > hullPts[i].first.z )
      min = i;
  }

  // Store the minimum point so that we know when we are done.
  hullPts[npts] = hullPts[min];

  do {

    // Make the point with the smallest z value first.
    if( min != m ) {
      std::pair< Point, unsigned int > tmpPt = hullPts[m];
      hullPts[m] = hullPts[min];
      hullPts[min] = tmpPt;
    }

    m++;
    min = m;

    // Find the next point that is the farthest to the right of all others.
    for( unsigned int i=min+1; i<npts+1; i++ ) {

      Vector v0 = (Vector) hullPts[m-1].first - (Vector) hullPts[min].first;
      Vector v1 = (Vector) hullPts[i  ].first - (Vector) hullPts[min].first;

      int CCW = ccwXZ( v0, v1 );

      // The next point that is the farthest to the right of all
      // others will be in a clockwise direction (i.e. the convex
      // hull is in the counter clockwise direction.
      if( CCW == dir )
        min = i;

      // In this case the points are co-linear so take the point
      // that is closest to the starting point.
      else if( hullPts[i].second != hullPts[m-1].second && CCW == 0 ) {
        v1 = (Vector) hullPts[m-1].first - (Vector) hullPts[i].first;

        if( v0.length2() > v1.length2() )
          min = i;
      }
    }

    // Stop when the first point is found again.
  } while( min != npts );

//   for( unsigned int i=0; i<m; ++i ) {
//     if( verboseFlag )
//       std::cerr << "convexHull " << hullPts[i].second << "   "
//                 << hullPts[i].first << std::endl;
//   }

//   if( verboseFlag )
//     std::cerr << std::endl;
}


bool FieldlineLib::hullCheck( std::vector< Point > &points, int &direction )
{
  // If one, two, or three points the ordering makes no difference and
  // it is convex.
  if( points.size() <= 2 ) {
      return true;
  }
  else if( points.size() == 3 ) {
    return ccwXZ( points[0] - points[1], points[2] - points[1] );
  }

  std::vector< std::pair < Point, unsigned int > > pts;
  std::vector< std::pair < Point, unsigned int > > hullPts;

  // Store the points and their original order in a temporary vector.
  for( unsigned int i=0; i<points.size(); i++ )
    pts.push_back( std::pair< Point, unsigned int >( points[i], i ) );

  unsigned int npts = chainHull_2D( pts, hullPts, 0 );

  unsigned int cc = 0;

  for( unsigned int i=0; i<npts-1; i++ )
  {
    // The hull points are returned in a clockwise order. So if the
    // point ordering is increasing then the point ordering is also in
    // clockwise.
    if( hullPts[i].second < hullPts[i+1].second+1)
      ++cc;
  }
   
  // The hull is found in a clockwise direction if the point count is
  // increasing.
  direction = ((cc > npts/2) ? 1 : -1);

  return (npts == points.size());
}


// Determine if a number is prime - return 1 otherwise the GCD.
unsigned int FieldlineLib::isPrime( unsigned int a )
{
  for( unsigned int i=a/2; i>1; --i )
  {
    if( a % i == 0 )
      return i;
  }

  return 1;
}


// Find the great comon denominator.
unsigned int FieldlineLib::GCD( unsigned int a, unsigned int b )
{
  if( a <= 0 || b <= 0 )
    return 0;

  if( a < b )
    { unsigned int tmp = a; a = b; b = tmp; }

  if( a % b )
    return GCD(b, a % b);
  else
    return b;
}


// Find the great comon denominator in a list of numbers.
unsigned int FieldlineLib::GCD( std::vector< unsigned int > values,
                                unsigned int &freq,
                                unsigned int minGCD )
{
  if( values.size() == 1 )
  {
    freq = 1;
    return values[0];
  }

  // Find the greatest common denominator between each value in the list.
  std::map< int, int > GCDCount;  
  std::map< int, int>::iterator ic;
  
  for( unsigned int i=0; i<values.size(); ++i )
  {
    for( unsigned int j=i+1; j<values.size(); ++j )
    {
      unsigned int gcd = GCD( values[i], values[j] );
      
      // Find the GCD excluding those smaller than the min.
      if( gcd >= minGCD )
      {
        ic = GCDCount.find( gcd );
        
        if( ic == GCDCount.end() )
          GCDCount.insert( std::pair<int, int>( gcd, 1) );
        else (*ic).second++;
      }
    }
  }
    
  // Find the most frequent greatest common denominator
  unsigned int gcd = 1;
  freq = 0;

  ic = GCDCount.begin();

  while( ic != GCDCount.end() )
  {
    // Two GCD have the same count so take the larger GCD.
    if( freq == (unsigned int)(*ic).second && gcd < (unsigned int)(*ic).first )
    {
      gcd = (*ic).first;
    }

    // GCD with a larger count
    else if( freq < (unsigned int)(*ic).second )
    {
      gcd = (*ic).first;
      
      freq = (*ic).second;
    }

    ++ic;
  }
  
  return gcd;
}


// Find the resonance via the great common denominator in a list of
// samples.
unsigned int FieldlineLib::
ResonanceCheck( std::vector< std::pair< unsigned int, double > > &stats,
                unsigned int baseResonance,
                unsigned int max_samples )
{
  unsigned int freq;
  std::vector< unsigned int > values;

  unsigned int max_groups = (unsigned int)stats.size() / max_samples;
  unsigned int mult = 0;
  double minPercent = 0.6;

  std::map< int, int > GCDCount;  
  std::map<int, int>::iterator ic;

  // For a secondary resonance to exists the group times the
  // baseResonance should equal the resonance.
  for( unsigned int group=max_groups; group>1; --group )
  {
    unsigned int num_entries = (unsigned int)stats.size() / group;

    values.resize( num_entries );

    for( unsigned int i=0; i<num_entries; ++i )
    {
      values[i] = stats[i].first;
    }

    // The group is not the GCD so find the GCD.
     unsigned int gcd = GCD( values, freq );

//     std::cerr << group  << "  " << gcd << std::endl;

    // Store the GCD and counts just in case while in the chaotic
    // regime we are not able to get the proper GCD,
    // i.e. gcd == group * baseResonance fails.
    ic = GCDCount.find( gcd );
        
    if( ic == GCDCount.end() )
      GCDCount.insert( std::pair<int, int>( gcd, 1) );
    else (*ic).second++;

    if( gcd == group * baseResonance )
    {
      unsigned int cc = 0;

      // Make sure the majority of the integer mutliples can be
      // divided by the resonance. Sometimes odd gcd values can sneak
      // in as well as odd resonances when near a separtrice.
      for( unsigned int i=0; i<num_entries; ++i )
      {
        if( values[i] % gcd == 0 )
          ++cc;
      }

      if( minPercent < (double) cc / (double) num_entries)
      {
        // If no previous resonance match take it and quit.
        if( mult == 0 )
        {
          minPercent = (double) cc / (double) num_entries;
          
          mult = group;
        }

        // Previous resonance it may be a secondary resonance but now
        // we found the primary.
        else if( isPrime(group) > 1 )
        {
          minPercent = (double) cc / (double) num_entries;
          
          mult = group;   
        }
        else
          continue;

        if( verboseFlag )
          std::cerr << "GCD resonance " << group << "  "
                    << minPercent << std::endl;
        break;
      }
    }
    else
    {
      // Assume the GCD is the group so as not to miss a possible
      // secondary resonance.
      unsigned int gcd = group * baseResonance;
      unsigned int cc = 0;

      // Make sure the majority of the integer mutliples can be
      // divided by the resonance. Sometimes odd gcd values can sneak
      // in as well as odd resonances when near a separtrice.
      for( unsigned int i=0; i<num_entries; ++i )
      {
        if( values[i] % gcd == 0 )
          ++cc;
      }

      if( minPercent < (double) cc / (double) num_entries)
      {
        minPercent = (double) cc / (double) num_entries;

        mult = group;

        if( verboseFlag )
          std::cerr << "Possible group resonance " << group << "  "
                    << minPercent << std::endl;
      }
      else
      {
        if( verboseFlag )
          std::cerr << "Failed group resonance " << group << "  "
                    << (double) cc / (double) num_entries
                    << std::endl;
      }
    }
  }

  // Probably chaotic so just use the most frequent GCD. 
  if( mult == 0 )
  {
    // Find the most frequent greatest common denominator
    freq = 0;

    ic = GCDCount.begin();

    while( ic != GCDCount.end() )
    {
      // Two GCD have the same count so take the larger GCD.
      if( freq == (unsigned int)(*ic).second && mult < (unsigned int)(*ic).first )
      {
        mult = (*ic).first;
      }
      
      // GCD with a larger count
      else if( freq < (unsigned int)(*ic).second )
      {
        mult = (*ic).first;
        
        freq = (*ic).second;
      }
      
      ++ic;
    }
  }
  
  // The multiplier is not greater than 1 then there is no second
  // order resonance present.
  if( mult > 1 )
  {
    unsigned int entries = (unsigned int)stats.size() / mult;

    // Divide the stats up into groups and keep the first two groups.
    if( mult > 2 )
    {
      if( verboseFlag )
        std::cerr << "Splitting into " << mult << " groups " << std::endl;
      
      unsigned int num_entries = 2 * (unsigned int)stats.size() / mult;
      
      if( num_entries && num_entries < stats.size() )
      {
        if( verboseFlag )
          std::cerr << "Erasing entries " << num_entries << " to "
               << stats.size() << std::endl;
        
        stats.resize( num_entries );
      }
    }
    else
    {
      if( verboseFlag )
        std::cerr << "No splitting needed for " << mult
                  << " groups " << std::endl;
    }

    // If present, thresholding should leave only higher order
    // resonances.
    thresholdStats( stats, true, 2 );
        
//  values.resize( stats.size() );
    values.resize( entries );
        
//  for( unsigned int i=0; i<stats.size(); ++i )
    for( unsigned int i=0; i<entries; ++i )
      values[i] = stats[i].first;
    
    return mult * baseResonance;
  }
  else
    return 1;
}

void FieldlineLib::
thresholdStats( std::vector< std::pair< unsigned int, double > >& stats,
                bool erase,
                unsigned int checkType )
{
  Otsu otsu;
  double maxVet, cutoff;

  otsu.getOtsuThreshold2( stats, cutoff, maxVet );

//  unsigned int cutoffIndex = stats.size();
  unsigned int cutoffIndex;

  if( erase ) cutoffIndex = (unsigned int)stats.size() / 2;
  else        cutoffIndex = (unsigned int)stats.size();

  for( unsigned int i=0; i<(unsigned int)stats.size(); ++i )
  {
//    if( stats[i].second <= cutoff ) 
    if( i < cutoffIndex ) 
    {
      if( verboseFlag )
        std::cerr << "period  " << stats[i].first << "  "
             <<  (checkType == 2 ? "distance  " :  "variance  " )
             << stats[i].second << "  "
             << std::endl;
    }
    else
    {
      cutoffIndex = i;
      break;
    }

    if( stats[i].second <= cutoff && stats[i+1].second > cutoff )
    {
      if( verboseFlag )
        std::cerr << "vet max  " << maxVet << "  "
                  << "Threshold " << cutoff << std::endl;
      }
  }

  for( unsigned int i=cutoffIndex; i<stats.size(); ++i )
  {
    if( verboseFlag )
    {
      // Print the cutoff entries
      std::cerr << "period  " << stats[i].first << "  "
           <<  (checkType == 2 ? "distance  " :  "variance  " )
           << stats[i].second
           << (erase ? "  **" : "  *") << std::endl;
    }

    if( stats[i].second <= cutoff && stats[i+1].second > cutoff )
    {
      if( verboseFlag )
        std::cerr << "vet max  " << maxVet << "  "
                  << "Threshold " << cutoff << std::endl;
    }
  }

  if( erase && cutoffIndex > 1 )
    stats.resize(cutoffIndex);
}


Point FieldlineLib::circle(Point &pt1, Point &pt2, Point &pt3)
{
  if (!IsPerpendicular(pt1, pt2, pt3) )
    return CalcCircle(pt1, pt2, pt3);   
  else if (!IsPerpendicular(pt1, pt3, pt2) )
    return CalcCircle(pt1, pt3, pt2);   
  else if (!IsPerpendicular(pt2, pt1, pt3) )
    return CalcCircle(pt2, pt1, pt3);   
  else if (!IsPerpendicular(pt2, pt3, pt1) )
    return CalcCircle(pt2, pt3, pt1);   
  else if (!IsPerpendicular(pt3, pt2, pt1) )
    return CalcCircle(pt3, pt2, pt1);   
  else if (!IsPerpendicular(pt3, pt1, pt2) )
    return CalcCircle(pt3, pt1, pt2);   
  else
    return Point(-1,-1,-1);
}


// Check the given point are perpendicular to x or y axis 
bool FieldlineLib::IsPerpendicular(Point &pt1, Point &pt2, Point &pt3)
{
  double d21z = pt2.z - pt1.z;
  double d21x = pt2.x - pt1.x;
  double d32z = pt3.z - pt2.z;
  double d32x = pt3.x - pt2.x;
        
  // checking whether the line of the two pts are vertical
  if (fabs(d21x) <= FLT_MIN && fabs(d32z) <= FLT_MIN)
    return false;
    
  else if (fabs(d21z) < FLT_MIN ||
           fabs(d32z) < FLT_MIN ||
           fabs(d21x) < FLT_MIN ||
           fabs(d32x)<= FLT_MIN )
    return true;

  else
    return false;
}

Point FieldlineLib::CalcCircle(Point &pt1, Point &pt2, Point &pt3)
{
  double d21z = pt2.z - pt1.z;
  double d21x = pt2.x - pt1.x;
  double d32z = pt3.z - pt2.z;
  double d32x = pt3.x - pt2.x;
        
  if (fabs(d21x) < FLT_MIN && fabs(d32z ) < FLT_MIN ) {

    return Point( ( 0.5*(pt2.x + pt3.x) ),
                  ( pt1.y ),
                  ( 0.5*(pt1.z + pt2.z) ) );
  }
        
  // IsPerpendicular() assure that xDelta(s) are not zero
  double aSlope = d21z / d21x;
  double bSlope = d32z / d32x;

  // checking whether the given points are colinear.    
  if (fabs(aSlope-bSlope) > FLT_MIN) {
    
    // calc center
    double x = (aSlope*bSlope*(pt1.z - pt3.z) +
                bSlope*(pt1.x + pt2.x) -
                aSlope*(pt2.x + pt3.x) ) / (2.0* (bSlope-aSlope) );

    return Point( x,
                  pt1.y,
                  -(x - (pt1.x+pt2.x)/2.0) / aSlope + (pt1.z+pt2.z)/2.0 );
  }

  return Point(0,0,0);
}


bool
FieldlineLib::IntersectCheck( std::vector< Point >& points,
                              unsigned int nbins,
                              unsigned int offset )
{
  // If the offset and point ordering have the same directions then
  // the next group is the offset. Otherwise if they are in the
  // opposite direction then nbins-offset is the next group.
  //int offsetDir = Dot( v0, v1 ) > 0.0 ? 1 : -1;
  for( unsigned int i=0, j=nbins; i<nbins && j<points.size(); ++i, ++j )
  {
    Point l0_p0 = points[i];
    Point l0_p1 = points[j];

    // The neighbor groups
    unsigned int skipNext     = (i + offset + nbins) % nbins;
    unsigned int skipPrevious = (i - offset + nbins) % nbins;

    for( unsigned int k=0, l=nbins; k<nbins && l<points.size(); ++k, ++l )
    {
      // Do not check the segment against itself or its immediate
      // neighbors.
      if( k == i || k == skipNext || k == skipPrevious )
        continue;

      Point l1_p0 = points[k];
      Point l1_p1 = points[l];

      if( intersect( l0_p0, l0_p1, l1_p0, l1_p1 ) == 1)
      {
//      std::cerr << nbins << "  "<< offset << "    "
//                << skipNext << "  " << skipPrevious << "     "
//                << i << "  " << j << "  " << k << "  " << l << std::endl;

        return false;
      }
    }
  }

  return true;
}


unsigned int FieldlineLib::Blankinship( unsigned int toroidalWinding,
                                        unsigned int poloidalWinding,
                                        unsigned int skip    /* = 1 */ )
{
  unsigned int offset = GCD( toroidalWinding, poloidalWinding );

  if( toroidalWinding == poloidalWinding )
    skip = 0;

  else if( toroidalWinding > 1 && poloidalWinding != 0 ) {
    //  To find the skip find the mutual primes via the
    //  Blankinship Algorithm.
    for( ; skip<toroidalWinding; skip++ )
      if( (skip * poloidalWinding) % toroidalWinding == offset )
        break;
    
    if( skip == toroidalWinding )
      skip = 0;
    
  } else {
    skip = 0;
  }

  return skip;
}


template< class TYPE >
void FieldlineLib::safetyFactorStats( std::vector< TYPE > &poloidalWindingCounts,
                                      double &averageSafetyFactor,
                                      double &stdDev )
{
  unsigned int start = poloidalWindingCounts.size() * .75;
  unsigned int stop  = poloidalWindingCounts.size();

  unsigned int cc = 0;

  averageSafetyFactor = 0;

  for( unsigned int i=start; i<stop; ++i, ++cc )
    averageSafetyFactor += (double) i / (double) poloidalWindingCounts[i];

  averageSafetyFactor /= (double) cc;

  stdDev = 0;

  for( unsigned int i=start; i<stop; ++i )
  {
    double diff =
      averageSafetyFactor - (double) i / (double) poloidalWindingCounts[i];

    stdDev += diff * diff;
  }

  stdDev /= (double) cc;

  stdDev = sqrt(stdDev);
}


// Smallest value is first.
bool compareWindingPairStats( const WindingPairStat s0,
                              const WindingPairStat s1 )
{
  return (s0 > s1);
}


void FieldlineLib::
SortWindingPairStats( std::vector< WindingPairStat > &windingPairStats,
                      bool reverseOrder )
{
  // Now sort the results.
  std::sort( windingPairStats.begin(), windingPairStats.end(),
             compareWindingPairStats );

  if( reverseOrder )
    reverse( windingPairStats.begin(), windingPairStats.end() );
}


void FieldlineLib::
RankWindingPairStats( std::vector< WindingPairStat > &windingPairStats,
                      bool LT )
{
  if( windingPairStats.size() < 2 )
    return;

  // Now rank the results using a numerical ordering while accounting
  // for ties.
  unsigned int rank = 0;

  windingPairStats[0].ranking = rank;

  for( unsigned int i=1; i<windingPairStats.size(); ++i )
  {
    if( ( LT && windingPairStats[i].stat < windingPairStats[i-1].stat) ||
        (!LT && windingPairStats[i].stat > windingPairStats[i-1].stat) )
      ++rank;

    windingPairStats[i].ranking = rank;
  }
}


void FieldlineLib::
poloidalWindingCheck( std::vector< unsigned int > &poloidalWindingCounts,
                      std::vector< WindingPairStat > &windingPairStats )
{
  windingPairStats.clear();

  unsigned int nsets = (unsigned int)poloidalWindingCounts.size();

  // The premise is that for a given toroidal winding the poloidal
  // winding should be consistent between each Nth punction point,
  // where N is the toroidal winding. For instance, if the toroidal
  // winding is 5 and the poloidal winding is 2. Then the pattern
  // could be:

  // 0 1 1 1 2 - 2 3 3 3 4 - 4 5 5 5 6

  // In this case the different between every 5th value (the toroidal
  // winding) should be 2 (the poloidal winding).

  unsigned int maxToroidalWinding = (unsigned int)poloidalWindingCounts.size() / 2;

  for( unsigned int toroidalWinding=1;
       toroidalWinding <= maxToroidalWinding;
       ++toroidalWinding )
  {

    std::map< int, int > differenceCount;
    std::map< int, int >::iterator ic;

    // Find all the differences and count each one.
    for( unsigned int i=0; i<nsets-toroidalWinding; ++i)
    {
      // Get the poloidal winding between two counts.
      unsigned int poloidalWinding =
        poloidalWindingCounts[i+toroidalWinding] - poloidalWindingCounts[i];

      // Find this difference in the list.
      ic = differenceCount.find( poloidalWinding );

      // Not found, new difference.
      if( ic == differenceCount.end() )
        differenceCount.insert( std::pair<int, int>( poloidalWinding, 1) );
      // Found this difference, increment the count.
      else (*ic).second++;
    }

    // Find the difference that occurs most often.
    unsigned int nMatches = 0;
    unsigned int poloidalWinding = 0;
    
    ic = differenceCount.begin();
    
    while( ic != differenceCount.end() )
    {
      if( nMatches < (unsigned int)(*ic).second )
      {
        poloidalWinding = (*ic).first;
        nMatches = (*ic).second;
      }

      ++ic;
    }

    double consistency = (double) nMatches / (double) (nsets-toroidalWinding);

    unsigned int t = toroidalWinding;
    unsigned int p = poloidalWinding;

    // If the toroidalWinding and poloidalWinding have a common
    // denominator find the greatest denominator and remove it.
//    if( t != p )
    if( 0 )
    {
      for( unsigned int d=p; d>1; --d) {
        if( t % d == 0 && p % d == 0 ) {
          t /= d;
          p /= d;
          
          d = p;
        }
      }
    }

    bool lowOrder = false;

    // Keep the low order toroidalWinding / poloidalWinding
    if( t != toroidalWinding && p != poloidalWinding )
    {
      for( unsigned int i=0; i<windingPairStats.size(); ++i )
      {
        if( windingPairStats[i].toroidal == t &&
            windingPairStats[i].poloidal == p )
        {
          lowOrder = true;

          // If the consistency happens to be higher for the higher
          // order keep it instead. Typically the lower order math is
          // better.
          if( windingPairStats[i].stat < consistency )
            windingPairStats[i].stat = consistency;
          
          break;
        }
      }
    }

    // Did not find a lower order match so record the set.
    if( ! lowOrder && toroidalWinding > 0 && poloidalWinding > 0 )
    {
      WindingPairStat windingPairStat;
      windingPairStat.toroidal = t;
      windingPairStat.poloidal = p;
      windingPairStat.stat = consistency;
      windingPairStat.ranking = 0;

      windingPairStats.push_back( windingPairStat );
    }
  }

  // Now sort the results.
  SortWindingPairStats( windingPairStats );
  RankWindingPairStats( windingPairStats );
}


double FieldlineLib::
calculateSumOfSquares( std::vector< Point >& points,
                       unsigned int period,
                       unsigned int checkType )
{
  // Find the sum of squares for each of the periods. If the period is
  // correct the sum of squares should be a small compared to a wrong
  // value.

  double sum = 0;
  double tSamples = 0;

  if( checkType == 4 )
  {
    // Get the sum of squares for each bin.
    for( unsigned int j=0; j<points.size(); ++j )
    {
      double len = (points[j]-points[(j+period)%points.size()]).length();
      sum += len;
      ++tSamples;
    }
  }
  else
  {
    for( unsigned int i=0; i<period; ++i )
    {
      Vector centroid(0,0,0);
      double length = 0;
      double nSamples = 0;
    
      if( checkType != 2 && checkType != 4 )
      {
        // Find the centroid of the points based on the period.
        for( unsigned int j=i; j<points.size(); j+=period )
        {
          // Centroid or Z difference.
          if( checkType == 0 || checkType == 1 )
          {
            centroid += points[j];
            ++nSamples;
          }

          // Find the average length of the distance between points
          // based on the period.
          else if( checkType == 3 && j>period )
          {
            length += (points[j]-points[j-period]).length();
            ++nSamples;
          }
        }
      
        if( checkType == 0 || checkType == 1 )
          centroid /= (double) nSamples;
      
        else if( checkType == 3 )
          length   /= (double) nSamples;
      }

      if( checkType == 2 || nSamples > 1 )
      {
        // Get the sum of squares for each bin.
        for( unsigned int j=i; j<points.size(); j+=period )
        {
          // Centroid difference
          if( checkType == 0 )
          {
            Vector diff = points[j] - centroid;
            sum += diff.length2();
            ++tSamples;
          }
          // Z difference
          else if( checkType == 1 )
          {
            double diff = points[j].z - centroid.z;
            sum += (diff * diff);
            ++tSamples;
          }

          // Length sum
          else if( checkType == 2 && j>period )
          {
            double len = (points[j]-points[j-period]).length();
            sum += len;
            ++tSamples;
          }
          
          // Length difference
          else if( checkType == 3 && j>period )
          {
            double diff = (points[j]-points[j-period]).length() - length;
            sum += (diff * diff);
            ++tSamples;
          }
        }
      }
    }
  }

  if( checkType == 2 || checkType == 4 )
    return sum;
  else
    return sum / tSamples;
}


// Largest value is first in the list.
int
compareSecondPair( const std::pair< unsigned int, double > s0,
                   const std::pair< unsigned int, double > s1 )
{
  return ( s0.second < s1.second );
}


void FieldlineLib::
periodicityStats( std::vector< Point >& points,
                  std::vector< std::pair< unsigned int, double > >& stats,
                  unsigned int max_period,
                  unsigned int checkType )
{
  stats.clear();

  // Find the base period variance.
  unsigned int best_period = (unsigned int)points.size();
  double best_var = 1.0e9;

  double base_var = calculateSumOfSquares( points, 1, checkType );

  if( verboseFlag )
    std::cerr << "Base variance  " << base_var << std::endl;
  
  // Find the period with the best variance.
  if( max_period == 0 )
    max_period = points.size() / 2.0;

  for( unsigned int i=1; i<=max_period; ++i ) 
  {
    double var = calculateSumOfSquares( points, i, checkType );
    
    stats.push_back( std::pair< unsigned int, double > (i, var ) );

    if( best_var > var ) 
    {
      best_var = var;
      best_period = i;
    }
  }

  if( verboseFlag )
    std::cerr << "Best period " << best_period << "  "
              <<  (checkType == 2 ? "distance  " :  "variance  " )
              << calculateSumOfSquares( points, best_period, checkType )
              << std::endl;

  if( stats.size() == 0 )
    stats.push_back( std::pair< unsigned int, double > (best_period, best_var ) );

  // Now sort the results.
  std::sort( stats.begin(), stats.end(), compareSecondPair );

//   for( unsigned int i=0; i<stats.size(); ++i )
//   {
//     std::cerr << stats[i].first << "  " << stats[i].second << "  "
//       << 0 << "  " << stats[i].second << "  "
//       << (i<stats.size()/2 ? "  " : "**")
//       << std::endl;
//   }
}


bool FieldlineLib::
rationalCheck( std::vector< Point >& points,
               unsigned int toroidalWinding,
               unsigned int &nnodes,
               float delta)
{
  // Look at the distance between the centroid of each toroidal group
  // and the points that are in it. If it is smaller then the distance
  // between the points that make up the fieldine is assumed to be on
  // a rational surface.

  // Note: this distance, delta, is defined indirectly based on the
  // integration step size.

  for( unsigned int i=0; i<toroidalWinding; i++ ) {

    // Get the local centroid for the toroidal group.
    Vector localCentroid(0,0,0);

    unsigned int npts = 0;

    for( unsigned int j=i; j<points.size(); j+=toroidalWinding, ++npts )
      localCentroid += (Vector) points[j];

    localCentroid /= (float) npts;

    for( unsigned int j=i; j<points.size(); j+=toroidalWinding ) {

      Vector vec = (Vector) points[j] - localCentroid;

      if( vec.length() > delta ) {
        return false;
      }
    }
  }

  nnodes = 1;

  return true;
}


bool FieldlineLib::
islandChecks( std::vector< Point >& points,
              unsigned int toroidalWinding,
              unsigned int &islands,
              unsigned int &nnodes,
              bool &complete )
{
  islands = 0;
  nnodes = 0;
  complete = 0;

  unsigned int npts;

  // Get the centroid for all of the points.
  Vector globalCentroid(0,0,0);

  for( unsigned int i=0; i<points.size(); i++ )
    globalCentroid += (Vector) points[i];

  globalCentroid /= (float) points.size();

  if( toroidalWinding == 1 )
  {
    // Check to see if it overlaps itself - it should or more points
    // need to be obtained.
    npts = (unsigned int)points.size() / toroidalWinding;

    for( unsigned int j=toroidalWinding, k=j+toroidalWinding;
         k<points.size();
         j+=toroidalWinding, k+=toroidalWinding )
    {     
      // See if the test point is between the first secton.
      Vector v0 = (Vector) points[0] - (Vector) points[k];
      Vector v1 = (Vector) points[1] - (Vector) points[k];
        
      if( Dot( v0, v1 ) < 0.0 )
      {
        npts = k / toroidalWinding;
        complete = true;
        break;
      }

      // See if the first point is between the test section.
      v0 = (Vector) points[0] - (Vector) points[j];
      v1 = (Vector) points[0] - (Vector) points[k];
      
      if( Dot( v0, v1 ) < 0.0)
      {
        npts = k / toroidalWinding;
        complete = true;
        break;
      }
    }
  }
  else
  {
    npts = toroidalWinding;
  }

 std::vector< Point > tmp_points;
  int direction;

  tmp_points.resize( npts );

  for( unsigned int j=toroidalWinding; j<npts; ++j )
    tmp_points[j] = points[j];

  bool convex = hullCheck( tmp_points, direction );

  // A convex hull for a toroidalWinding of 1 is defined as not
  // being an island.
  if( toroidalWinding == 1 && convex )
  {
//  std::cerr << "toroidalWinding of 1 and convex hull" << std::endl;
    return false;
  }

//  std::cerr << "Convex  " << convex << std::endl;

  // Check for islands. Islands will exists if there is a change in
  // direction of the connected points relative to a base point. If
  // the hull is convex the base point may the centroid of all of the
  // points or based upon a point that is perpendicular to the
  // principal axis of the group of points.
 std::vector< unsigned int > nodes(toroidalWinding);

  for( unsigned int i=0; i<toroidalWinding; i++ )
  {
    Vector baseCentroid;

    // If the hull is convex the global centroid can be used because
    // all of the islands will surround it in a radial manner.
    if( convex )
    {
      baseCentroid = globalCentroid;
    }
    // Otherwise use an offset from the local centroid of each group.
    else // !concave
    {
      // Get the local centroid for the group.
      Vector localCentroid(0,0,0);

      unsigned int npts = 0;

      for( unsigned int j=i; j<points.size(); j+=toroidalWinding, npts++ )
        localCentroid += (Vector) points[j];

      localCentroid /= (float) npts;

//      std::cerr << "localCentroid  " << localCentroid << std::endl;

      // Get the principal axes of the points.
      float Ixx = 0.0;
      float Ixz = 0.0;
      float Izz = 0.0;

      float maxDist = 0;

      // Get the moments of intertial for each point. It assumed that
      // everything is in the Y plane as such the moments of intertial
      // along the Y axis are zero.
      for( unsigned int j=i; j<points.size(); j+=toroidalWinding )
      {

        Vector vec = (Vector) points[j] - localCentroid;

        if( maxDist < vec.length() )
          maxDist = vec.length();

        Ixx += vec.z*vec.z;
        Ixz -= vec.x*vec.z;
        Izz += vec.x*vec.x;
      }

      // Short cut to the principal axes because the Y moments of
      // intertial are zero.
      float alpha = atan( 2.0 * Ixz / (Ixx - Izz) ) / 2.0;

      float Pxx = Ixx + Ixz * sin(alpha       )/cos(alpha       );
      float Pzz = Izz + Ixz * cos(alpha+M_PI/2)/sin(alpha+M_PI/2);

//     std::cerr << "Principal axes of intertial "
//       << Pxx << "  " << Pzz << "  " << Pxx/Pzz << "  " << Pzz/Pxx << std::endl;

      // For a toroidalWinding of 1 if the moment of interia of one axis is
      // similar to the other axis then use the centroid directly.
      if( toroidalWinding == 1 && Pxx / Pzz < 5.0 && Pzz / Pxx < 5.0 )
      {
//      std::cerr << "Using local centroid " << std::endl;
        baseCentroid = localCentroid;
      }
      else
      {
        // Find the approximate center if points were projected onto a
        // circle.
        unsigned int npts = 0;
        Vector center(0,0,0);
 
        for( unsigned int j=i+toroidalWinding;
             j<points.size()-toroidalWinding;
             j+=toroidalWinding ) {
          
          unsigned int j_1 = j - toroidalWinding;
          unsigned int j1  = j + toroidalWinding;


          center += (Vector) circle( points[j_1],
                                     points[j  ],
                                     points[j1 ] );

          ++npts;
        }

        center /= (float) npts;

//      std::cerr << "center  " << center << std::endl;

        // Use the principal axes to get an offset from the local
        // centroid which gives a point outside the island.

        // The direction along the axis is determined by where the
        // center is located.

//      std::cerr << "localCentroid  " << localCentroid << std::endl;

        if( Pxx > Pzz )
        {
          if( Dot( center - localCentroid, 
                   Vector( cos(alpha), 0, sin(alpha) ) ) )
            baseCentroid = localCentroid +
              Vector(  cos(alpha), 0, sin(alpha) ) * maxDist;
          else
            baseCentroid = localCentroid -
              Vector(  cos(alpha), 0, sin(alpha) ) * maxDist;
        }
        else
        {
          if( Dot( center - localCentroid, 
                   Vector( -sin(alpha), 0, cos(alpha) ) ) )
            baseCentroid = localCentroid +
              Vector( -sin(alpha), 0, cos(alpha) ) * maxDist;
          else
            baseCentroid = localCentroid -
              Vector( -sin(alpha), 0, cos(alpha) ) * maxDist;
        }
      }
    }

//    std::cerr << "baseCentroid  " << baseCentroid << std::endl;

    unsigned int turns = 0;
    unsigned int firstTurn = 0;
    unsigned int   midTurn = 0;
    unsigned int  lastTurn = 0;

    // Get the direction based on the first two points.
    Vector v0 = (Vector) points[i                ] - baseCentroid;
    Vector v1 = (Vector) points[i+toroidalWinding] - baseCentroid;

    bool localCCW, lastCCW = (ccwXZ( v0, v1 ) == 1);

    // Get the direction based on the remaining points.
    for( unsigned int j=i+2*toroidalWinding;
         j<points.size();
         j+=toroidalWinding )
    {
      v0 = v1;
      v1 = (Vector) points[j] - baseCentroid;

      // Points are parallel - this result should not occur as the
      // base centroid point should be on the concave side of the
      // island.
      if (ccwXZ( v0, v1 ) == 0 )
        localCCW = lastCCW;
      else
        localCCW = (ccwXZ( v0, v1 ) == 1);

      // A change in direction indicates that an island is present.
      if( localCCW != lastCCW )
      {
        lastCCW = localCCW;

        // Record the number of turns. Ideally three turns will be
        // found which indicates that the island is complete. Although
        // the island could be complete with only two turns.
        ++turns;

        if( turns == 1 )
          ++islands;

        if( turns == 1 )     firstTurn = j - toroidalWinding;
        else if( turns == 2 )  midTurn = j - toroidalWinding;
        else if( turns == 3 ) lastTurn = j - toroidalWinding;

        if( turns == 3 )
          break;
      }
    }

    unsigned int overlap = 0;

    // Approximate number of points in the islands. 
    nodes[i] = (unsigned int)points.size() / toroidalWinding;

    // Check to see if the island overlaps itself. Only use this if
    // there are two or more turns. Otherwise the it can appear that
    // it overlaps itself earlier.
    if( turns >= 2 )
    {
      // Okay initial approximatation. It really depending on where
      // the turns are relative to each other/
      if( turns == 3 )
        nodes[i] = (lastTurn - firstTurn) / toroidalWinding + 1;

      // Look for the overlap between the ends of the island which are
      // flatter.
      unsigned int base, start;      
      unsigned int quarter_dist =
        ((midTurn - firstTurn) / 4 / toroidalWinding) * toroidalWinding;

      // If the start is in the middle section use it.
      if( quarter_dist <= firstTurn && firstTurn <= quarter_dist * 3 )
      {
        base = i;
        start = midTurn;
      }

      // If the start is before the middle section use the first
      // middle section.
      else if( firstTurn >= quarter_dist * 2 )
      {
        base = firstTurn - quarter_dist * 2;
        start = midTurn;
      }

      // If the start is after the middle section use the next middle
      // secction.
      else
      {
        base = firstTurn + quarter_dist;

        if( lastTurn )
          start = lastTurn;
        else
          start = midTurn;
      }

      if( turns && verboseFlag )
        std::cerr << "Island " << i << " - "
                  << " nodes " << nodes[i]
                  << " turns " << turns
                  << " first " << firstTurn
                  << " mid " << midTurn
                  << " last " << lastTurn
                  << " base " << base
                  << " start " << start
                  << std::endl;

      for( unsigned int j=start, k=j+toroidalWinding;
           k<points.size();
           j+=toroidalWinding, k+=toroidalWinding ) {

        // See if the test point is between the first secton.
        v0 = ((Vector) points[base                ] -
              (Vector) points[k                   ]);
        v1 = ((Vector) points[base+toroidalWinding] -
              (Vector) points[k                   ]);
        
        if( Dot( v0, v1 ) < 0.0 ) {
          overlap = k;
          nodes[i] = (overlap-base) / toroidalWinding;
          turns = 3;
          complete = true;
          break;
        }

        // See if the first point is between the test section.
        v0 = (Vector) points[base] - (Vector) points[j];
        v1 = (Vector) points[base] - (Vector) points[k];
        
        if( Dot( v0, v1 ) < 0.0 ) {
          overlap = k;
          nodes[i] = (overlap-base) / toroidalWinding;
          turns = 3;
          complete = true;
          break;
        }
      }

      // Do a distance search based on the overlaping sections having
      // two points in close proximity while at the same time the
      // overlaping sections have approximately the same length.
      if( 1 || !complete )
      {
        unsigned int stop;

        // Because the curve may overlap the start multiple times
        if( lastTurn )
        {
          stop = lastTurn +
            ((midTurn - firstTurn) / 4 / toroidalWinding) * toroidalWinding;
        }
        else
        {
          stop = (unsigned int)points.size()-toroidalWinding;
        }

        // Get the maximum distance from the first point so that the
        // search is not conducted near neighbor points.
        double max_gap = 0;
        
        for( unsigned int j=i+toroidalWinding;
             j<points.size()-toroidalWinding;
             j+=toroidalWinding )
        {
          double tmp_gap = (points[i] - points[j]).length();

          if( max_gap < tmp_gap )
            max_gap = tmp_gap;
        }

        max_gap *= 0.5;

        double min_gap = 1.0e9;  // Distance between points.
        unsigned int min_gap_index = 0;

        double min_len = 1.0e9;  // Difference in segments.
        unsigned int min_len_index = 0;

        double min_dist = 1.0e9; // Minimum distance.
        unsigned int min_dist_index = 0;

        double min_dist2 = 1.0e9;// Second minimum distance.
        unsigned int min_dist2_index = 0;

        bool past_max_gap = false;

        // Base length of the first segment.
        double len = (points[i] - points[i+toroidalWinding]).length();

        for( unsigned int j=i+toroidalWinding;
             j<points.size()-toroidalWinding && j<stop;
             j+=toroidalWinding )
        {
          double tmp_gap = (points[i] - points[j]).length();

          if( tmp_gap > max_gap )
            past_max_gap = true;

          double tmp_len = (points[j] - points[j+toroidalWinding]).length();
          double dist = sqrt(tmp_gap*tmp_gap + (tmp_len-len)*(tmp_len-len));

          if( past_max_gap )
          {
            if( min_gap > tmp_gap )
            {
              min_gap = tmp_gap;
              min_gap_index = j;
            }

            if( min_len > tmp_len )
            {
             min_len = tmp_len;
             min_len_index = j;
            }
          
            if( min_dist > dist )
            {
//          std::cerr << j << "  "
//               << tmp_gap << "  "
//               << fabs(tmp_len - len) << "  " 
//               << dist << "  "
//               << std::endl;

              min_dist2 = min_dist;
              min_dist2_index = min_dist_index;

              min_dist = dist;
              min_dist_index = j;
            }

            else if( min_dist2 > dist )
            {
//          std::cerr << j << "  "
//               << tmp_gap << "  "
//               << fabs(tmp_len - len) << "  " 
//               << dist << "  "
//               << std::endl;

              min_dist2 = dist;
              min_dist2_index = j;
            }
          }
        }

        if( !complete )
        {
          overlap = min_dist_index;
          nodes[i] = (min_dist_index-i)/toroidalWinding;
          turns = 3;
          complete = true;
        }

        std::cerr << "gap " << min_gap_index << "  "
                  << "len " << min_len_index << "  "
                  << "dist " << min_dist_index << "  "
                  << "dist2 " << min_dist2_index << "  "
                  << "nodes " << (min_dist_index-i)/toroidalWinding
                  << std::endl;
      }
    }
    
    if( turns && verboseFlag )
      std::cerr << "Island " << i << " - "
                << " nodes " << nodes[i]
                << " turns " << turns 
                << " first " << firstTurn
                << " mid " << midTurn
                << " last " << lastTurn
                << " overlap " << overlap
                << " complete " << complete
                << std::endl;

    nnodes += nodes[i];

    // Started on a surface but jumpped to an island - i.e. a
    // separtrice
    if( !complete && lastTurn &&
        ((int) (firstTurn / toroidalWinding) - (int) (nodes[i]/2)) > 2 )
    {
      --islands;

      if( turns && verboseFlag )
        std::cerr << "  !!!!!!!!!!!!!!! Separatrice !!!!!!!!!!!!!!!"
                  << std::endl;
    }
  }

  // Get the average number of nodes.
  nnodes /= (float) toroidalWinding;

  if( islands ) {

    unsigned int cc = 0;
    unsigned int new_nnodes = 0;

    for( unsigned int i=0; i<toroidalWinding; i++ )
    {
      if( nodes[i] < nnodes - 1 || nnodes + 1 < nodes[i] )
      {
        cc++;
      }
      else
        new_nnodes += nodes[i];
    }

    if( cc == 1 )
    {
      nnodes = new_nnodes / (float) (toroidalWinding-cc);
    }

    for( unsigned int i=0; i<toroidalWinding; i++ )
    {
      if( nodes[i] < nnodes - 1 || nnodes + 1 < nodes[i] )
      {
        if( verboseFlag )
        {        
          std::cerr << "Appears to be islands but not self consistent, "
                    << "average number of nodes  " << nnodes << "   (  ";
          
          for( unsigned int i=0; i<toroidalWinding; i++ )
            std::cerr << nodes[i] << "  ";

          std::cerr << ")" << std::endl;
        }

        break;
      }
    }  
  }

  return (bool) islands;
}


void
FieldlineLib::getPunctures( std::vector< Point > &ptList,
                            Vector planeN,
                            std::vector< Point > &puncturePts )
{
  unsigned int startIndex = 0;

  // Set up the plane equation.
  Vector planePt(0,0,0);
  double plane[4];
            
  plane[0] = planeN.x;
  plane[1] = planeN.y;
  plane[2] = planeN.z;
  plane[3] = planePt.dot(planeN);
  
  
  // So to get the winding groups consistant start examining
  // the fieldline in the same place for each plane.
  Point lastPt, currPt = ptList[startIndex];
  double lastDist, currDist = planeN.dot( currPt ) - plane[3];
            
  for( unsigned int i=startIndex+1; i<ptList.size(); ++i )
  {
    lastPt = currPt;
    currPt = Vector(ptList[i]);
                
    lastDist = currDist;
    currDist = Dot( planeN, currPt ) - plane[3];
    
    // First look at only points that intersect the plane.
    if( SIGN(lastDist) != SIGN(currDist) ) 
    {
      Vector dir(currPt-lastPt);
                    
      double dot = Dot(planeN, dir);
                    
      // If the segment is in the same direction as the plane then
      // find where it intersects the plane.
      if( dot > 0.0 )
      {
        Vector w = lastPt - planePt;
        
        double t = -Dot(planeN, w ) / dot;
        
        Point point = Point(lastPt + dir * t);
        
        puncturePts.push_back( point );
      }
    }
  }
} 


void
FieldlineLib::getFieldlineBaseValues( std::vector< Point > &ptList,
                                      std::vector< double > &timeList,
                                      std::vector< Point > &poloidal_puncture_pts,
                                      std::vector< Point > &ridgeline_points,
                                      std::vector< double > &rotationalSums,
                                      std::vector< unsigned int > &poloidalWindingCounts,
                                      float &delta,
                                      unsigned int OLineToroidalWinding )
{
  double toroidalBase = OLineToroidalWinding ? OLineToroidalWinding : 1;

  Point axisPt( 0.1, 0.1, 0.2 );

  unsigned int totalPoloidalPuncturePts = 0;

  poloidal_puncture_pts.clear();
  ridgeline_points.clear();

  rotationalSums.clear();
  poloidalWindingCounts.clear();

  delta = 0.0;

  Point tmpPt, currPt = ptList[0], nextPt = ptList[1];
  double /*tmpTime,*/ currTime = timeList[0], nextTime = timeList[1];
  Vector planePt( 0, 0, 0 );

  if( verboseFlag )
  {
    std::cerr << "-----------------------------------------------------------------"
         << std::endl
         << "Analyzing  " << ptList[0] << "  " << timeList[0] << "  "
         << "with  " << ptList.size() << " fieldline points"
         << std::endl;
  }

  // Set up the Y plane equation as the base poloidal puncture
  // analysis takes place in the X-Z plane.
  Vector planeNY( 0, 1, 0 );
  double planeY[4];

  planeY[0] = planeNY.x;
  planeY[1] = planeNY.y;
  planeY[2] = planeNY.z;
  planeY[3] = Dot( planePt, planeNY );

  double currDistY, nextDistY = Dot(planeNY, nextPt) - planeY[3];

  // Rotational sum
  double phiBase = atan2(currPt.y, currPt.x);
  double currPhi;
  double nextPhi = -1; ///TODO: check on uninitialized var

  double phiTotal = phiBase;

  double dPhi, dTheta;
  double rotationalSum = 0;
  double ridgelineRotationalSum = 0;

  double nextR, currR, deltaR = 0;
  double nextZ, currZ, deltaZ = 0;

  double prevTheta=0, currTheta=0, nextTheta=0;
        
  bool zeroAngle = false;
  
  double maxZ = currPt.z;


  if( OLineToroidalWinding )
  {
    // Axis point for the current point.
    axisPt = getAxisPt( currPt, phiBase, toroidalBase );

    // Current theta relative to the axis.
    tmpPt = currPt - axisPt;
    currZ = tmpPt[2]; tmpPt[2] = 0;
    currR = tmpPt.length();
    
    currTheta = atan2( currZ, currR );

    // Current and next phi angle.
    currPhi = atan2(currPt.y, currPt.x);
    nextPhi = atan2(nextPt.y, nextPt.x);

    dPhi = nextPhi - currPhi;

    if (dPhi > M_PI)
      dPhi -= (2.0*M_PI);
    
    else if (dPhi < -M_PI)
      dPhi += (2.0*M_PI);

    phiBase += dPhi;
    phiTotal += dPhi;
    
    if( phiBase > toroidalBase * 2.0 * M_PI )
      phiBase -= toroidalBase * 2.0 * M_PI;
    
    else if( phiBase < -toroidalBase * 2.0 * M_PI )
      phiBase += toroidalBase * 2.0 * M_PI;

    // Axis point for the next point.
    axisPt = getAxisPt( nextPt, phiBase, toroidalBase );
    
    // Next theta relative to the axis.
    tmpPt = nextPt - axisPt;
    nextZ = tmpPt[2]; tmpPt[2] = 0;
    nextR = tmpPt.length();
    
    nextTheta = atan2( nextZ, nextR );

    // Rotational amount.
    dTheta = nextTheta - currTheta;
      
    if (dTheta > M_PI)
      dTheta -= (2.0*M_PI);
    
    else if (dTheta < -M_PI)
      dTheta += (2.0*M_PI);
    
    rotationalSum += dTheta;
    ridgelineRotationalSum += dTheta;
  }
  else
  {
    // Current R Z values relative to the origin
    tmpPt = currPt;
    currZ = tmpPt[2]; tmpPt[2] = 0;
    currR = tmpPt.length();

    // Next R Z values relative to the origin
    tmpPt = nextPt;
    nextZ = tmpPt[2]; tmpPt[2] = 0;
    nextR = tmpPt.length();

    // Difference R Z values
    deltaR = nextR - currR;
    deltaZ = nextZ - currZ;

    // Current theta value.
    currTheta = atan2(deltaZ, deltaR);

    if( deltaR == 0 && deltaZ == 0 )
      zeroAngle = true;
  }


  // Now collect the points.
  unsigned int npts = 0;

  for( unsigned int i=2; i<ptList.size(); ++i)
  {    
    currPt = nextPt;
    currTime = nextTime;

    nextPt = ptList[i];
    nextTime = timeList[i];

    if( maxZ < currPt.z )
    {
      maxZ = currPt.z;
    }

    // Save the distance between points to use for finding periodic
    // fieldlines (i.e. rational surfaces and re-connection).
    Vector s = (Vector) nextPt - (Vector) currPt;

    double ds = s.length();

    delta += ds;
    ++npts;

    // Poloidal plane distances.
    currDistY = nextDistY;
    nextDistY = Dot( planeNY, nextPt ) - planeY[3];

    // First look at only points that intersect the poloidal plane.
    if( SIGN(currDistY) != SIGN(nextDistY) ) 
    {
      Vector dir(nextPt-currPt);
      
      double dot = Dot(planeNY, dir);
      
      // If the segment is in the same direction as the poloidal plane
      // then find where it intersects the plane.
      if( dot > 0.0 )
      {
        Vector w = (Vector) currPt - planePt;
        
        double t = -Dot(planeNY, w ) / dot;
        
        Point point = Point(currPt + dir * t);

        bool savePt = false;

        // Double Poincare puncture test.
        if( puncturePeriod )
        {
          double time = currTime + (nextTime-currTime) * t;

          // Get the number of periods traversed.
          double nPeriods = time / puncturePeriod;

          // Get the factional part - should be close to zero for an
          // even period.
          double intPart, fracPart = modf(nPeriods, &intPart);

          if( fracPart < puncturePeriodTolerance ||
              1.0-puncturePeriodTolerance < fracPart )
          {
            savePt = true;
          }
        }
        else
          savePt = true;

        
        if( savePt &&
            totalPoloidalPuncturePts % (unsigned int) toroidalBase == 0 )
        {
          poloidal_puncture_pts.push_back( point );
          
          poloidalWindingCounts.push_back( fabs(rotationalSum) /
                                           (2.0 * M_PI) );

//        std::cerr << poloidalWindingCounts.size() << "  "
//                  << rotationalSum << "  "
//                  << fabs(rotationalSum) / (2.0 * M_PI)
//                  << std::endl;

          rotationalSums.push_back( fabs(rotationalSum) );
        }
            
        ++totalPoloidalPuncturePts;
      }
    }

    // Calculate the absolute phi angle which will happen when
    // crossing the -π to π boundary or visa-versa.
    if( OLineToroidalWinding )
    {
      currPhi = nextPhi;
      nextPhi = atan2(nextPt.y, nextPt.x);

      dPhi = nextPhi - currPhi;

      if (dPhi > M_PI)
        dPhi -= (2.0*M_PI);
      
      else if (dPhi < -M_PI)
        dPhi += (2.0*M_PI);

      phiBase += dPhi;
      phiTotal += dPhi;

      if( phiBase > toroidalBase * 2.0 * M_PI )
        phiBase -= toroidalBase * 2.0 * M_PI;

      else if( phiBase < -toroidalBase * 2.0 * M_PI )
        phiBase += toroidalBase * 2.0 * M_PI;

      axisPt = getAxisPt( nextPt, phiBase, toroidalBase );

      currTheta = nextTheta;

      tmpPt = nextPt - axisPt;
      nextZ = tmpPt[2]; tmpPt[2] = 0;
      nextR = tmpPt.length();

      nextTheta = atan2( nextZ, nextR );

      dTheta = nextTheta - currTheta;
      
      if (dTheta > M_PI)
        dTheta -= (2.0*M_PI);
      
      else if (dTheta < -M_PI)
        dTheta += (2.0*M_PI);

      dTheta = fabs( dTheta );
    }
    else
    {
      // Values for the rotational transform summation.
      currR = nextR;
      currZ = nextZ;
      
      tmpPt = nextPt;
      nextZ = tmpPt[2]; tmpPt[2] = 0;
      nextR = tmpPt.length();
      
      // Calculate the rotational transform summation.
      deltaR = nextR - currR;
      deltaZ = nextZ - currZ;
      
      // Check for a zero angle coming in on the right.
      if( !zeroAngle )
      {
        prevTheta = currTheta;
        
        if( deltaR == 0 && deltaZ == 0 )
          zeroAngle = true;
      }
      
      // In zero mode so look for a good angle coming in on the right.
      else if( zeroAngle && (deltaR != 0 || deltaZ != 0) )
      {
        zeroAngle = false;
      }
    
      if( !zeroAngle )
      {
        currTheta = atan2(deltaZ, deltaR);
        
        dTheta = currTheta - prevTheta;
      
        if (dTheta > M_PI)
          dTheta -= (2.0*M_PI);
        
        else if (dTheta < -M_PI)
          dTheta += (2.0*M_PI);
      }
      else
        dTheta = 0;
    }

    rotationalSum += dTheta;
    ridgelineRotationalSum += dTheta;

    if( fabs(ridgelineRotationalSum) >= 2.0 * M_PI )
    {
      ridgelineRotationalSum -=
        2.0 * M_PI * SIGN(ridgelineRotationalSum);

      ridgeline_points.push_back( Point( (float) ridgeline_points.size(),
                                         0,
                                         maxZ) );
      maxZ = 0;
    }
  }

  if( verboseFlag )
    std::cerr << "phiTotal " << phiTotal << "  "
              << "number of rotations " << phiTotal/(2.0*M_PI)
              << std::endl;

  // At this point all of the poloidal and toroidal puncture points
  // have been found.
  if( verboseFlag )
    std::cerr << poloidal_puncture_pts.size() << " poloidal puncture pts, "
              << ridgeline_points.size() << " ridgeline pts " << std::endl;

  // Get the average distance between puncture points for finding
  // periodic fieldlines (i.e. rational surfaces and re-connection).
  delta  /= (float) npts;
}


Point FieldlineLib::getAxisPt( Point pt, double phiTest, double toroidalBase )
{
  if( OLineAxisPts.empty() )
  {
    OLineAxisIndex = 1;

    unsigned int index = 1;

    Point currPt,  nextPt;

    double phiBase, phiTotal = 0; ///TODO check for fix of uninitialized variable
    double dPhi;

    double currPhi = atan2(currPt.y, currPt.x);
    double nextPhi = atan2(nextPt.y, nextPt.x);

    FILE *fp = fopen( OLineAxisFileName.c_str(), "r" );

    if( fp == NULL )
    {
      std::cerr << "Trying to perform O-line analysis "
                << "but the O-line axis file can not be openned." << std::endl;

      return Point(0,0,0);
    }

    while( !feof(fp) )
    {
      currPt = nextPt;
      currPhi = nextPhi;

      if( fscanf( fp, "%lf %lf %lf", &nextPt.x, &nextPt.y, &nextPt.z ) != 3 )
      {
        if( feof(fp) == 0 )
        {
          std::cerr << "Trying to perform O-line analysis "
                    << "but the O-line axis file can not be read." << std::endl;
          
          OLineAxisPts.clear();
          OLineAxisPhiAngles.clear();
          
          return Point(0,0,0);
        }
        else if( feof(fp) )
          break;
      }

      OLineAxisPts.push_back( nextPt );

      if( OLineAxisPhiAngles.empty() )
      {
        phiBase = atan2(nextPt.y, nextPt.x);
        phiTotal = 0; //phiBase;
      }
      else
      {
        nextPhi = atan2(nextPt.y, nextPt.x);

        dPhi = nextPhi - currPhi;
        
        if (dPhi > M_PI)
          dPhi -= (2.0*M_PI);
        
        else if (dPhi < -M_PI)
          dPhi += (2.0*M_PI);
        
        phiBase += dPhi;
        phiTotal += dPhi;
      }

      if( phiTotal/(2.0*M_PI) / toroidalBase < 1.00 )
        ++index;

      OLineAxisPhiAngles.push_back( phiBase );
    }
     
    fclose (fp);

    if( verboseFlag )
      std::cerr << "total points " << OLineAxisPhiAngles.size() << "  "
                << "phiTotal " << phiTotal << "  "
                << "number of rotations " << phiTotal/(2.0*M_PI)
                << std::endl;

    if( fabs( phiTotal/(2.0*M_PI) / toroidalBase ) < 1.0 )
    {
      std::cerr << "Too few points in the axis "
                << "expected " << toroidalBase << " rotations, "
                << "got " << phiTotal/(2.0*M_PI) <<  " rotations."
                << std::endl;

      if( index > OLineAxisPhiAngles.size() ) 
        index = (unsigned int)OLineAxisPhiAngles.size();

      if( index < OLineAxisPhiAngles.size() ) 
        OLineAxisPhiAngles.erase( OLineAxisPhiAngles.begin()+index,
                                  OLineAxisPhiAngles.end() );
    }

    else if( index < OLineAxisPhiAngles.size() )
    {
      std::cerr << "Too too many points in the axis "
                << "expected " << toroidalBase << " rotations, "
                << "got " << phiTotal/(2.0*M_PI) <<  " rotations. "
                << "Truncating"
                << std::endl;

      OLineAxisPhiAngles.erase( OLineAxisPhiAngles.begin()+index,
                                OLineAxisPhiAngles.end() );
    }
  }

  Point axisPt( 0, 0, 0 );

  for( unsigned int i=1; i<OLineAxisPts.size(); ++i )
  {
    if( (phiTest >= 0 &&
         OLineAxisPhiAngles[OLineAxisIndex-1] <= phiTest &&
         phiTest < OLineAxisPhiAngles[OLineAxisIndex]) ||
        (phiTest <= 0 &&
         OLineAxisPhiAngles[OLineAxisIndex-1] >= phiTest &&
         phiTest > OLineAxisPhiAngles[OLineAxisIndex]) )
    {
      double t =
        (phiTest - OLineAxisPhiAngles[OLineAxisIndex-1]) /
        (OLineAxisPhiAngles[OLineAxisIndex]-OLineAxisPhiAngles[OLineAxisIndex-1]);
      
      axisPt = OLineAxisPts[OLineAxisIndex-1] +
        (OLineAxisPts[OLineAxisIndex]-OLineAxisPts[OLineAxisIndex-1]) * t;

      return axisPt;
    }
    
    ++OLineAxisIndex;
    
    if( OLineAxisIndex == OLineAxisPts.size() )
      OLineAxisIndex = 1;
  }

  unsigned int i = (unsigned int)OLineAxisPhiAngles.size()-1;
  
  if( phiTest < OLineAxisPhiAngles[0] && OLineAxisPhiAngles[i] > 0 )
    phiTest += 2.0 * M_PI;
  else if( phiTest > OLineAxisPhiAngles[0] && OLineAxisPhiAngles[i] < 0 )
    phiTest -= 2.0 * M_PI;

  for( unsigned int i=1; i<OLineAxisPts.size(); ++i )
  {
    if( (phiTest >= 0 &&
         OLineAxisPhiAngles[OLineAxisIndex-1] <= phiTest &&
         phiTest < OLineAxisPhiAngles[OLineAxisIndex]) ||
        (phiTest <= 0 &&
         OLineAxisPhiAngles[OLineAxisIndex-1] >= phiTest &&
         phiTest > OLineAxisPhiAngles[OLineAxisIndex]) )
    {
      double t =
        (phiTest - OLineAxisPhiAngles[OLineAxisIndex-1]) /
        (OLineAxisPhiAngles[OLineAxisIndex]-OLineAxisPhiAngles[OLineAxisIndex-1]);
      
      axisPt = OLineAxisPts[OLineAxisIndex-1] +
        (OLineAxisPts[OLineAxisIndex]-OLineAxisPts[OLineAxisIndex-1]) * t;

      return axisPt;
    }
    
    ++OLineAxisIndex;
    
    if( OLineAxisIndex == OLineAxisPts.size() )
      OLineAxisIndex = 1;
  }

  i = (unsigned int)OLineAxisPhiAngles.size()-1;
  
  std::cerr << "Can not find axis point for point " << pt
            << " at phi " << phiTest << "    "
            << OLineAxisPhiAngles[0] << "  " << OLineAxisPhiAngles[1] << "  "
            << OLineAxisPhiAngles[i-1] << "  " << OLineAxisPhiAngles[i] << "  "
            << OLineAxisIndex
            << std::endl;
  
  return axisPt;
}


void FieldlineLib::
GetBaseWindingPairs( std::vector< unsigned int > &poloidalWindingCounts,
                     std::vector< Point > &poloidal_puncture_pts,
                     std::vector< WindingPairStat > &baseWindingPairStats,
                     double &windingPairConfidence,
                     unsigned int &toroidalWindingMax,
                     unsigned int &poloidalWindingMax,
                     unsigned int &drawableBaseWindingPairIndex )
{
  // Check the consistency of the poloidal winding counts. 
  poloidalWindingCheck( poloidalWindingCounts, baseWindingPairStats );

  // Report the winding number pairs.
  std::vector< WindingPairStat >::iterator iter = baseWindingPairStats.begin();
  
  // Get the first set that passes the intersection test and passes
  // the user setable match limit. Default is 0.90 (90%)
  for( unsigned int i=0; i<baseWindingPairStats.size(); ++i, ++iter )
  {
    if( baseWindingPairStats[i].stat < windingPairConfidence )
      break;

    double local_safetyFactor =
      (double) baseWindingPairStats[i].toroidal /
      (double) baseWindingPairStats[i].poloidal;

    if( toroidalWindingMax < baseWindingPairStats[i].toroidal )
      toroidalWindingMax = baseWindingPairStats[i].toroidal;

    if( poloidalWindingMax < baseWindingPairStats[i].poloidal )
      poloidalWindingMax = baseWindingPairStats[i].poloidal;

    unsigned int windingGroupOffset =
      Blankinship( baseWindingPairStats[i].toroidal,
                   baseWindingPairStats[i].poloidal );

    if( IntersectCheck( poloidal_puncture_pts,
                        baseWindingPairStats[i].toroidal,
                        windingGroupOffset ) )
    {
      // Record the index of the first winding pair that does not self
      // intersect.
      if( drawableBaseWindingPairIndex == (unsigned int)-1 ) ///TODO: -1 cannot happen with unsigned int so casting it -1 of unsigned int which should be max value (CHECK & FIX!)
        drawableBaseWindingPairIndex = i;
        
      if( verboseFlag )
        std::cerr << "**Drawable winding pair "
                  << baseWindingPairStats[i].toroidal << ","
                  << baseWindingPairStats[i].poloidal << "  ("
                  << local_safetyFactor << " - "
                  << "consistency " << 100.0 * baseWindingPairStats[i].stat
                  << "%)" << "  "
                  << windingGroupOffset << "  "
                  << std::endl;
    }
    else
    {
      // Debug info
      if( verboseFlag )
        std::cerr << "Undrawable winding pair "
                  << baseWindingPairStats[i].toroidal << ","
                  << baseWindingPairStats[i].poloidal << "  ("
                  << local_safetyFactor << " - "
                  << "consistency " << 100.0 * baseWindingPairStats[i].stat
                  << "%)" << "  "
                  << windingGroupOffset << "  "
                  << std::endl;
    }
  }

  // Remove the winding number sets that are below the limit.
  if( iter != baseWindingPairStats.end() )
    baseWindingPairStats.erase( iter, baseWindingPairStats.end() );
}


void FieldlineLib::
GetPeriodWindingPairs( std::vector< WindingPairStat > &baseWindingPairStats,
                       std::vector< WindingPairStat > &periodWindingPairStats,
                       std::vector< std::pair< unsigned int,
                                               double > > &toroidalStats,
                       std::vector< std::pair< unsigned int,
                                               double > > &poloidalStats )
{
  bool pairFound;

  for( unsigned int i=0; i<baseWindingPairStats.size(); ++i )
  {
    pairFound = false;

    // Look for a toroidal winding
    for( unsigned int j=0; j<toroidalStats.size(); ++j )
    {
      if( toroidalStats[j].first == baseWindingPairStats[i].toroidal )
      {
        // Look for a poloidal winding
        for( unsigned int k=0; k<poloidalStats.size(); ++k )
        {
          if( poloidalStats[k].first == baseWindingPairStats[i].poloidal &&
           
              // Make sure the ratio of both periods is the same. This
              // ratio is important for island chains.
              toroidalStats[j].first / baseWindingPairStats[i].toroidal ==
              poloidalStats[k].first / baseWindingPairStats[i].poloidal )
          {
            pairFound = true;

            WindingPairStat windingPairStat = baseWindingPairStats[i];
            
            windingPairStat.stat = sqrt((double) (j*j+k*k));
            windingPairStat.ranking = 0;

            periodWindingPairStats.push_back( windingPairStat );
            break;
          }     
        }
      }

      if( pairFound )
        break;
    }
  }

  // Now sort the results based on the index Euclidian distance.
  SortWindingPairStats( periodWindingPairStats, true );
  RankWindingPairStats( periodWindingPairStats, false );

  for( unsigned int i=0; i<periodWindingPairStats.size(); ++i )
  {
    if( verboseFlag && i<10 )
      std::cerr << "Period based winding pair:  " 
                << periodWindingPairStats[i].toroidal << ","
                << periodWindingPairStats[i].poloidal << "  "
                << "Distance  " << periodWindingPairStats[i].stat << "  "
                << "Rank  " << periodWindingPairStats[i].ranking << "  "
                << std::endl;
  }
}


double
FieldlineLib::fieldlinePeriod( std::vector< Point > &ptList,
                               std::vector< double > &timeList,
                               double maxPeriod )
{
  int nTimeSteps = (double) (maxPeriod / (timeList[1] - timeList[0]));

  double bestPeriod = 0;
  double maxDotProd = 0;

  std::vector< double > dotProds(nTimeSteps+1);

  dotProds[0] = 0;

  for( unsigned int i=1; i<=nTimeSteps; i++ )
  {
    dotProds[i] = 0;

    avtVector v0 = ptList[1] - ptList[0];

    for( unsigned int j=1, k=2; k<ptList.size(); ++j, ++k )
    {
      avtVector v1 = ptList[k] - ptList[j];

      dotProds[i] += Dot(v0, v1);

      v0 = v1;
    }

    if( maxDotProd < dotProds[i] )
    {
      bestPeriod = i*(timeList[1] - timeList[0]);

      maxDotProd = dotProds[i];

      std::cerr << i << "  " << i*(timeList[1] - timeList[0]) << "  "
                << dotProds[i] << "  " 
                << std::endl;
    }
  }

  return bestPeriod;
}


// int
// compareSecondWindingPair( const std::pair< WindingPair, unsigned int > s0,
//                           const std::pair< WindingPair, unsigned int > s1 )
// {
// //  return ( s0.second < s1.second );
//   return ( s0.first.toroidal < s1.first.toroidal );
// }


void
FieldlineLib::fieldlineProperties( std::vector< Point > &ptList,
                                   std::vector< double > &timeList,
                                   FieldlineProperties &fi,
                                   unsigned int overrideToroidalWinding,
                                   unsigned int overridePoloidalWinding,
                                   unsigned int maxToroidalWinding,
                                   double windingPairConfidence,
                                   double rationalSurfaceFactor,
                                   bool detectIslandCenters,
                                   unsigned int OLineToroidalWinding,
                                   std::string OLineAxisFilename )
{
  if( ptList.empty() )
  {
    fi.type = FieldlineProperties::UNKNOWN_TYPE;

    fi.analysisState = FieldlineProperties::TERMINATED;

    fi.safetyFactor = 0;
    fi.toroidalWinding = 0;
    fi.poloidalWinding = 0;
    fi.poloidalWindingP = 0;
    fi.toroidalResonance = 0;
    fi.poloidalResonance = 0;
    fi.windingPairs.clear();
    fi.topWindingPairs.clear();
    fi.windingGroupOffset = 0;
    fi.islands = 0;
    fi.islandGroups = 0;
    fi.nnodes  = 0;

    fi.nPuncturesNeeded  = 0;

    fi.seedPoints.clear();

    return;
  }

  std::vector< Point > poloidal_puncture_pts;
  std::vector< Point > poloidal_puncture_pts2;
  std::vector< Point > ridgeline_points;

  std::vector< double > rotationalSums;
  std::vector< unsigned int > poloidalWindingCounts;

  // Save the distance between fieldline points to use for finding
  // periodic fieldlines (i.e. rational surfaces and re-connection).
  float rationalSurfaceTolerance = 0.0;

  // Some info for when doing O-Line axis analysis.
  double toroidalBase = OLineToroidalWinding ? OLineToroidalWinding : 1;
  OLineAxisFileName = OLineAxisFilename;

  getFieldlineBaseValues( ptList,
                          timeList,
                          poloidal_puncture_pts,
                          ridgeline_points,
                          rotationalSums,
                          poloidalWindingCounts,
                          rationalSurfaceTolerance,
                          OLineToroidalWinding );

  // Check to see if puncture points were actually added.
  if( fi.analysisState == FieldlineProperties::ADDING_POINTS &&
      fi.numPunctures == poloidal_puncture_pts.size() )
  {
    fi.analysisState = FieldlineProperties::TERMINATED;

    return;
  }

  fi.numPunctures = (unsigned int)poloidal_puncture_pts.size();

  if( ptList.empty() ||
      poloidal_puncture_pts.empty() ||
      ridgeline_points.empty() )
  {
    fi.type = FieldlineProperties::UNKNOWN_TYPE;

    fi.analysisState = FieldlineProperties::TERMINATED;

    fi.safetyFactor = 0;
    fi.toroidalWinding = 0;
    fi.poloidalWinding = 0;
    fi.poloidalWindingP = 0;
    fi.toroidalResonance = 0;
    fi.poloidalResonance = 0;
    fi.windingPairs.clear();
    fi.topWindingPairs.clear();
    fi.windingGroupOffset = 0;
    fi.islands = 0;
    fi.islandGroups = 0;
    fi.nnodes  = 0;

    fi.nPuncturesNeeded  = 0;

    fi.seedPoints.clear();

    return;
  }

  rationalSurfaceTolerance *= rationalSurfaceFactor;

  int helicity = ccwXZ( ptList[1] - ptList[0], ptList[1] - ptList[2] );

  // Get the safety factor.
//   for( unsigned int i=0; i<poloidal_puncture_pts.size(); ++i )
//     std::cerr << i << "  "
//            << (2.0 * M_PI * (i+1)) / fabs(rotationalSums[i]) << std::endl;

  // Last pair estimation for the safety factor.
  double LRS_SafetyFactor = (2.0 * M_PI * poloidal_puncture_pts.size()) /
    fabs(rotationalSums[poloidal_puncture_pts.size()-1]);

  // Average estimation for the safety factor.
//   double averageRotationalSum, stdDev;

//   safetyFactorStats( rotationalSums, averageRotationalSum, stdDev );
//   averageRotationalSum *= 2.0 * M_PI;

  // Least squares estimation for the safety factor.
  double LSRS_SafetyFactor;

  least_square_fit( rotationalSums, LSRS_SafetyFactor );
  LSRS_SafetyFactor *= 2.0 * M_PI;

  double safetyFactor = LSRS_SafetyFactor;

  if( verboseFlag )
  {
    std::cerr << "Limit Rotational Sum Safety Factor    "
              << LRS_SafetyFactor << std::endl
//            << "Average Rotational Sum Safety Factor         "
//            << averageRotationalSum << std::endl
              << "Least Square Rotational Sum Safety Factor    "
              << LSRS_SafetyFactor << std::endl

              << "Using safety factor " << safetyFactor << std::endl;
  }

  // Start the analysis.
  FieldlineProperties::FieldlineType type =
    FieldlineProperties::UNKNOWN_TYPE;
  FieldlineProperties::AnalysisState analysisState =
    FieldlineProperties::UNKNOWN_ANALYSIS;
  FieldlineProperties::SearchState searchState = fi.searchState;
//    FieldlineProperties::UNKNOWN_SEARCH;

  unsigned int toroidalWinding = 0, poloidalWinding = 0;
  unsigned int toroidalWindingP = 0, poloidalWindingP = 0;
  unsigned int toroidalResonance = 1, poloidalResonance = 1;
  unsigned int toroidalPeriod = 0, poloidalPeriod = 0;

  unsigned int windingGroupOffset = 0;
  unsigned int nnodes = 0;
  unsigned int islands = 0, islandGroups = 0;
  unsigned int nPuncturesNeeded = 0;

  unsigned int drawableBaseWindingPairIndex = -1;

  unsigned int toroidalWindingMax = 0, poloidalWindingMax = 0;

  std::vector< WindingPairStat > baseWindingPairStats,
    //approximateWindingPairs,
    periodWindingPairStats, mergedWindingPairStats;

  double searchDelta = fi.searchDelta;
  avtVector searchNormal = fi.searchNormal;
  Point lastSeedPoint = fi.lastSeedPoint;

  std::map< WindingPair, unsigned int > topWindingPairs = fi.topWindingPairs;

  std::vector< Point > islandSeedPts;

  GetBaseWindingPairs( poloidalWindingCounts, poloidal_puncture_pts,
                       baseWindingPairStats,
                       windingPairConfidence,
                       toroidalWindingMax, poloidalWindingMax,
                       drawableBaseWindingPairIndex );


  // No drawable winding pair or match consistency is less than the
  // user set value. Run more expensive tests to identify the
  // fieldline.
  if( drawableBaseWindingPairIndex == (unsigned int)-1 ) ///TODO: unsigned int cannot be -1 so switch to max value (CHECK & FIX!)
  {
    if( verboseFlag )
      std::cerr << "Poor consistency - probably chaotic" << std::endl;

    if( overrideToroidalWinding && overridePoloidalWinding )
    {
      windingGroupOffset = Blankinship( overrideToroidalWinding,
                                        overridePoloidalWinding );

      fi.analysisState = FieldlineProperties::TERMINATED;
      fi.type = FieldlineProperties::CHAOTIC;
    }
    else
    {
      windingGroupOffset = 0;

      fi.analysisState = FieldlineProperties::TERMINATED;
      fi.type = FieldlineProperties::CHAOTIC;
    }
    
    fi.safetyFactor = safetyFactor;
    fi.toroidalWinding = overrideToroidalWinding;
    fi.poloidalWinding = overridePoloidalWinding;
    fi.toroidalWindingP = 0;
    fi.poloidalWindingP = 0;
    fi.toroidalPeriod = 0;
    fi.poloidalPeriod = 0;
    fi.toroidalResonance = 1;
    fi.poloidalResonance = 1;
    fi.windingPairs.clear();
    fi.windingGroupOffset = windingGroupOffset;
    fi.islands = 0;
    fi.islandGroups = 0;
    fi.nnodes  = 0;

    fi.nPuncturesNeeded  = 0;

    fi.seedPoints.clear();

    return;
  }

  std::vector< std::pair< unsigned int, double > > toroidalStats;
  std::vector< std::pair< unsigned int, double > > toroidalStats2;
  std::vector< std::pair< unsigned int, double > > poloidalStats;

  // Find the best toroidal periodicity. For a flux surface the period
  // will be the toroidal winding number. For an island chain the
  // period will be the toroidal winding number times the number of
  // nodes.
  if( verboseFlag )
    std::cerr << "Toroidal Winding via "
              << poloidal_puncture_pts.size() << "  "
              << "poloidal punctures, "
              << "max period " << toroidalWindingMax
              << std::endl;

  periodicityStats( poloidal_puncture_pts, toroidalStats,
                    toroidalWindingMax, 2 );

  // Find the best poloidal periodicity. For a flux surface the period
  // will be the poloidal winding number. For an island chain the
  // period will be the poloidal winding number times the number of
  // nodes.
  if( verboseFlag )
    std::cerr << "Poloidal Winding via "
              << ridgeline_points.size() << "  "
              << "ridgeline points, "
              << "max period " << poloidalWindingMax
              << std::endl;
  
  periodicityStats( ridgeline_points, poloidalStats, poloidalWindingMax, 1 );
  

  // Form a second winding number list that is ranked based on the
  // euclidian distance of each of the period lists.
  GetPeriodWindingPairs( baseWindingPairStats, periodWindingPairStats,
                         toroidalStats, poloidalStats );


  // Merge the base and period winding pairs together based on the
  // Euclidian distance via the index.
  for( unsigned int i=0; i<baseWindingPairStats.size(); ++i )
  {
    int ii = baseWindingPairStats[i].ranking;

    // Search for the same sibling pair in the best rational
    // approximation winding pair list.
    int jj = -1;

    for( unsigned int j=0; j<periodWindingPairStats.size(); ++j )
    {
      if( baseWindingPairStats[i].toroidal == periodWindingPairStats[j].toroidal &&
          baseWindingPairStats[i].poloidal == periodWindingPairStats[j].poloidal )
      {
        jj = periodWindingPairStats[j].ranking;
        break;
      }
    }

    // Found a matching pair so compute the index Euclidian distance.
    if( jj != -1 )
    {
      WindingPairStat windingPairStat = baseWindingPairStats[i];

      windingPairStat.stat = sqrt((double)(ii*ii+jj*jj));
      windingPairStat.ranking = 0;
            
      mergedWindingPairStats.push_back( windingPairStat );
    }
  }        

  // Now sort the results based on the Euclidian distance via the index.
  SortWindingPairStats( mergedWindingPairStats, true );
  RankWindingPairStats( mergedWindingPairStats, false );


  // Perform resonance checks on the toriodal and poloidal stats.
  if( mergedWindingPairStats.size() > 1 )
  {
    // The base values are the first order resonances if present.
    toroidalResonance = ResonanceCheck( toroidalStats, 1, 3 );
    poloidalResonance = ResonanceCheck( poloidalStats, 1, 3 );

    // It is possible when near the chaotic regime that the resonance
    // will not fall out from the toroidal and poloidal stats. As
    // such, try to get it from the winding pairs which while limited
    // should be more consistant.
    if( (float) mergedWindingPairStats[0].toroidal /
        (float) mergedWindingPairStats[0].poloidal !=
        (float) toroidalResonance / (float) poloidalResonance )
    {
      if( verboseFlag )
        std::cerr << "Windings found "
                  << mergedWindingPairStats[0].toroidal << "  "
                  << mergedWindingPairStats[0].poloidal << "  "
                  << "Resonances found "
                  << toroidalResonance << "  "
                  << poloidalResonance << std::endl;

      // Can not remember this case ... and why???????????? The first
      // test assures that the toroidal resonance is "good". While the
      // second test assures that the first pair is not the resonance.
      if ( mergedWindingPairStats[0].toroidal % toroidalResonance == 0 &&
           mergedWindingPairStats[0].poloidal != poloidalResonance )
      {
        // Get GCD from winding pairs ...
        unsigned int freq;
        std::vector< unsigned int > values;
        unsigned int minResonance;

        unsigned int lastDrawable = 0;

        if( toroidalResonance == 1 && poloidalResonance == 1 )
          minResonance = 2;
        else
          minResonance = 2;

        values.resize( mergedWindingPairStats.size() );

        // Check the toroidal windings ...
        for( unsigned int i=0; i<mergedWindingPairStats.size(); ++i )
        {
          values[i] = mergedWindingPairStats[i].toroidal;

          windingGroupOffset = Blankinship( mergedWindingPairStats[i].toroidal,
                                            mergedWindingPairStats[i].poloidal );

          bool drawable = IntersectCheck( poloidal_puncture_pts,
                                          mergedWindingPairStats[i].toroidal,
                                          windingGroupOffset );

          if( drawable )
          {
            std::cerr << "Drawable " << i << "  winding pair  "
                      << mergedWindingPairStats[i].toroidal << "  "
                      << mergedWindingPairStats[i].poloidal << "  "
                      << std::endl;

            lastDrawable = i;
          }
        }
        
        if( lastDrawable == 0 )
          lastDrawable = (unsigned int)mergedWindingPairStats.size();


        values.resize( lastDrawable );

        toroidalResonance = GCD( values, freq, minResonance );

        // Check the poloidal windings ...
        for( unsigned int i=0; i<mergedWindingPairStats.size(); ++i )
          values[i] = mergedWindingPairStats[i].poloidal;

        values.resize( lastDrawable );

        poloidalResonance = GCD( values, freq, minResonance );

        if( verboseFlag )
          std::cerr << lastDrawable << "  " << mergedWindingPairStats.size()
                    << "  Winding pair resonance "
                    << toroidalResonance << "  "
                    << poloidalResonance << std::endl;

        // Still no match so give up.
        if( (float) mergedWindingPairStats[0].toroidal /
            (float) mergedWindingPairStats[0].poloidal !=
            (float) toroidalResonance / (float) poloidalResonance )
        {
          if( verboseFlag )
            std::cerr << "Resonances don't match first winding pair - chaotic??" << std::endl;

          toroidalResonance = 1;
          poloidalResonance = 1;
        }
      }
      else
      {
        toroidalResonance = 1;
        poloidalResonance = 1;
      }
    }
  }
  else
  {
    toroidalResonance = 1;
    poloidalResonance = 1;
  }

  // For a 1,1 case a 2,2 resonance can be found as the winding pairs
  // can easily be divided evenly between ood, even values. As such,
  // throwout the case where a 2,2 resonce is found for a 1,1 surface.
  if( toroidalResonance == poloidalResonance &&
      toroidalResonance == 2 && poloidalResonance == 2 )
  {
    toroidalResonance = 1;
    poloidalResonance = 1;
  }

  // Get the resonance GCD which gives the indication that there are
  // secondary islands.
  unsigned int resonanceGCD = GCD( toroidalResonance, poloidalResonance );

  if( verboseFlag )
    std::cerr << "Winding Pair "
              << mergedWindingPairStats[0].toroidal << ","
              << mergedWindingPairStats[0].poloidal << "  "
              << "GCD = " << GCD( mergedWindingPairStats[0].toroidal,
                                  mergedWindingPairStats[0].poloidal ) << "   "
              << "Toroial, Poloidal resonances = "
              << toroidalResonance << "," << poloidalResonance << "  "
              << "GCD = " << resonanceGCD << "   "
              << std::endl;


  int drawableRank  = -1;  // Rank of the first drawable widing pair
  int drawableIndex = -1;  // Index of the first drawable widing pair
  std::vector< unsigned int > drawableIndexs;

  std::vector< WindingPair > windingPairs;

  // Loop through all the merged winding pairs and check for the cusp
  // condition. As such the rotational sum will be off by 1 or 2. We
  // can figure out the correct poloidal winding using the inverse
  // Blankenship via the toroidal windong the winding group offset.
  for( unsigned int i=0; i<mergedWindingPairStats.size(); ++i )
  {
    WindingPair windingPair( mergedWindingPairStats[i].toroidal,
                             mergedWindingPairStats[i].poloidal );

    windingPairs.push_back( windingPair );

    if( i < 10 )
    {
      std::map< WindingPair, unsigned int >::iterator ic;
      
      ic = topWindingPairs.find( windingPair );
      
      // Not found, new windingPair.
      if( ic == topWindingPairs.end() )
        topWindingPairs.insert( std::pair< WindingPair,
                                unsigned int >( windingPair, 1) );
      
      // Found this windingPair, increment the count.
      else
        (*ic).second++;
    }

    windingGroupOffset = Blankinship( mergedWindingPairStats[i].toroidal,
                                      mergedWindingPairStats[i].poloidal );

    bool drawable = IntersectCheck( poloidal_puncture_pts,
                                    mergedWindingPairStats[i].toroidal,
                                    windingGroupOffset );

    // If islands check the based period.
    unsigned int gcd = toroidalResonance == 1 ?
      1 : GCD( mergedWindingPairStats[i].toroidal,
               mergedWindingPairStats[i].poloidal );

    // Does not appear to work with islands about the O-point.
    if( toroidalBase == 1.0 && drawable &&
        mergedWindingPairStats[i].toroidal/gcd >= 5 &&

        // This eliminates islands that are not classified as islands
        // because of chaos.
        GCD( mergedWindingPairStats[i].toroidal,
             mergedWindingPairStats[i].poloidal ) !=
        mergedWindingPairStats[i].poloidal )
    {
      // Load in the the first N puncture points where N is the toroidal
      // winding.
      poloidal_puncture_pts2 = poloidal_puncture_pts;

      poloidal_puncture_pts2.resize(mergedWindingPairStats[i].toroidal/gcd);

      bool tmpVerboseFlag = verboseFlag;

      verboseFlag = false;

      // Find the closest neighbor offset index by looking for the
      // minimum distance between points. The minimum distance will
      // occur when the two points are adjacent.
      periodicityStats( poloidal_puncture_pts2, toroidalStats2,
                        mergedWindingPairStats[i].toroidal/gcd/2, 4 );

      verboseFlag = tmpVerboseFlag;

      // The winding group offsest should match the closest neighbor
      // offset but due to cusps it may not.
      if( windingGroupOffset != toroidalStats2[0].first && 
          windingGroupOffset != (mergedWindingPairStats[i].toroidal/gcd -
                                 toroidalStats2[0].first) )
      {      
        if( verboseFlag && (drawable || i < 10) )
          std::cerr << "Examining "
                    << mergedWindingPairStats[i].toroidal << ","
                    << mergedWindingPairStats[i].poloidal << "  "
                    << "Offset via "
                    << poloidal_puncture_pts2.size() << "  "
                    << "poloidal points, "
                    << "max period " << mergedWindingPairStats[i].toroidal/gcd/2
                    << std::endl;

        // Redo the stats for verbose.
        periodicityStats( poloidal_puncture_pts2, toroidalStats2,
                          mergedWindingPairStats[i].toroidal/gcd/2, 4 );
        
        if( verboseFlag && (drawable || i < 10) )
          std::cerr << "Rotational sum error, "
                    << "expected windingGroupOffset " << toroidalStats2[0].first
                    << " or " << mergedWindingPairStats[i].toroidal/gcd - toroidalStats2[0].first
                    << " got " << windingGroupOffset << std::endl
                    << "Incorrect poloidal winding is " << mergedWindingPairStats[i].poloidal/gcd
                    << "  ";

        // Low order Blankinship offset which is an indication of what
        // value to check.

        unsigned int wMax;
        unsigned int windingGroupOffset;

        // At the present time we believe this error only happens for
        // cases where q < 1. As such, we can not directly use the
        // Blankinship value as it is a low order value. Such that we
        // believe that the correct poloidal winding for these cases
        // will be mergedWindingPairs[i].toroidal + wMax. Where wMax
        // is the lowest of the two possible values.
 

        // The best windingGroupOffset can be either either side
        // i.e. the offset or toroidalWinding - offset. As such take
        // the one that results in the lowest order Balnkinship Value.

        if( Blankinship( mergedWindingPairStats[i].toroidal/gcd,
                         toroidalStats2[0].first ) <
            Blankinship( mergedWindingPairStats[i].toroidal/gcd,
                         mergedWindingPairStats[i].toroidal/gcd -
                         toroidalStats2[0].first ) )
        {
          windingGroupOffset = toroidalStats2[0].first;
        }
        else
        {
          windingGroupOffset = (mergedWindingPairStats[i].toroidal/gcd -
                                toroidalStats2[0].first );
        }

        // Lowest order Balnkinship value.
        wMax = Blankinship( mergedWindingPairStats[i].toroidal/gcd,
                            windingGroupOffset );

        if( verboseFlag )
          std::cerr << "windingGroupOffset " << windingGroupOffset << "  "
                    << "wMax " << wMax << std::endl;

        wMax *= gcd;

        // Search through the possible choices to find an offset.
        // NOTE: ARS believes that the correct poloidal winding for
        // these cases will be mergedWindingPairs[i].toroidal + wMax
        if( mergedWindingPairStats[i].poloidal >
            mergedWindingPairStats[i].toroidal )
          mergedWindingPairStats[i].poloidal =
            mergedWindingPairStats[i].toroidal + wMax;
        else if( abs((int)(mergedWindingPairStats[i].poloidal-wMax)) <
                 abs((int)(mergedWindingPairStats[i].poloidal-(mergedWindingPairStats[i].toroidal-wMax)) ) )
          mergedWindingPairStats[i].poloidal = wMax;
        else
          mergedWindingPairStats[i].poloidal =
            mergedWindingPairStats[i].toroidal - wMax;

        if( verboseFlag && (drawable || i < 10) )
          std::cerr << "Correct poloidal winding is "
                    << mergedWindingPairStats[i].poloidal
//                  << ( w == wMax ? " as predicted" : "" )
                    << std::endl;

//        unsigned int w;


//         if( w > wMax )
//         {
//           if( verboseFlag && (drawable || i < 10) )
//           {
//             std::cerr << "Can not correct poloidal winding." << std::endl;
          
//             for( w=1; w<=wMax; ++w )
//             {
//               if( verboseFlag && (drawable || i < 10) )
//                 std::cerr << "Checking poloidal value " << w;
              
//               std::cerr << "  Failed  "
//                         << Blankinship( mergedWindingPairStats[i].toroidal,
//                                         mergedWindingPairStats[i].toroidal+w )
//                         << std::endl;
//             }
//           }
//         }
      }
    }

    // Now have valid winding pairs.
    if( verboseFlag && (drawable || i < 10) )
      std::cerr << "Final " << i << "  "
                << (drawable ? "Drawable " : "Rejected ") 
                << "winding pair:  " 
                << mergedWindingPairStats[i].toroidal << ","
                << mergedWindingPairStats[i].poloidal << "  "
                << "Distance  " << mergedWindingPairStats[i].stat << "  "
                << "Rank  " << mergedWindingPairStats[i].ranking << "  ";

    if( (drawableRank == -1 ||
         drawableRank == mergedWindingPairStats[i].ranking) &&

        // Ignore the user requested if too large
        (maxToroidalWinding == 0 ||
         mergedWindingPairStats[i].toroidal <= maxToroidalWinding) &&

        // Keep only those that are drawable.
        drawable )
    {
      drawableIndexs.push_back( i );

      // If a tie take the higher order winding.
      if( toroidalWinding < mergedWindingPairStats[i].toroidal )
      {
        toroidalWinding = mergedWindingPairStats[i].toroidal;

        drawableRank = mergedWindingPairStats[i].ranking;
        drawableIndex = i;
      }
    }

    if( verboseFlag && (drawable || i<10) )
      std::cerr << std::endl;
  }

  // Sort the winding pairs to look for a bimodal set.
  // std::vector< std::pair< WindingPair, unsigned int > >
  //   topWindingPairList( topWindingPairs.begin(), topWindingPairs.end() );
        
  // std::sort( topWindingPairList.begin(), topWindingPairList.end(),
  //            compareSecondWindingPair );
    
//   for( unsigned int j=0; j<topWindingPairList.size(); ++j )
//     std::cerr << "LINE " << __LINE__ << "  "
//            << topWindingPairList[j].first.toroidal << "  "
//            << topWindingPairList[j].first.poloidal << "  "
//            << topWindingPairList[j].second << "  "
//            << std::endl;


  // No winding pair found is drawable. 
  if( drawableIndex == -1 )
  {
    // Look for an overriding winding.
    if( overrideToroidalWinding && overridePoloidalWinding )
    {
      toroidalWinding = overrideToroidalWinding;
      poloidalWinding = overridePoloidalWinding;

      windingGroupOffset = Blankinship( toroidalWinding, poloidalWinding );

      nnodes = (unsigned int)poloidal_puncture_pts.size() / 2;

      // Check to see if the fieldline is periodic. I.e. on a rational
      // surface.  If within "delta" of the distance the fieldline is
      // probably on a rational surface.

      // NOTE: Define the rational based on the lowest order of teh
      // surface. Thus divide by the GCD.
      if( rationalCheck( poloidal_puncture_pts,
                         toroidalWinding / GCD( toroidalWinding,
                                                poloidalWinding ),
                         nnodes, rationalSurfaceTolerance ) ) 
      {
        type = FieldlineProperties::RATIONAL;
        islands = 0;
        
        analysisState = FieldlineProperties::OVERRIDE;
        
        fi.rationalSurfaceTolerance = rationalSurfaceTolerance;

        if( verboseFlag )
          std::cerr << "Appears to be a rational surface "
                    << rationalSurfaceTolerance << std::endl;
      }
      else
      {
        fi.analysisState = FieldlineProperties::OVERRIDE;
        fi.type = FieldlineProperties::CHAOTIC;
      }

      fi.safetyFactor = safetyFactor;
      fi.toroidalWinding = toroidalWinding;
      fi.poloidalWinding = poloidalWinding;
      fi.toroidalWindingP = 0;
      fi.poloidalWindingP = 0;
      fi.toroidalPeriod = 0;
      fi.poloidalPeriod = 0;
      fi.toroidalResonance = 1;
      fi.poloidalResonance = 1;
      fi.windingPairs = windingPairs;
      fi.topWindingPairs = topWindingPairs;
      fi.windingGroupOffset = windingGroupOffset;
      fi.islands = 0;
      fi.islandGroups = 0;
      fi.nnodes = poloidal_puncture_pts.size() / 2;

      fi.nPuncturesNeeded  = 0;

      fi.seedPoints.clear();

      return;
    }
    else
    {
      if( poloidal_puncture_pts.size() < fi.maxPunctures )
      {
        if( verboseFlag )
          std::cerr << "Garbage matches adding more points" << std::endl;
        
        fi.nPuncturesNeeded = fi.maxPunctures * 1.25;
        fi.analysisState = FieldlineProperties::ADDING_POINTS;
      }
      else
      {
        fi.nPuncturesNeeded = 0;

        fi.analysisState = FieldlineProperties::TERMINATED;
      }

      return;
    }
  }

  // At this point there is a valid winding pair. It is not known if
  // it is an island or flux surface.

  toroidalPeriod   = mergedWindingPairStats[drawableIndex].toroidal;
  poloidalPeriod   = mergedWindingPairStats[drawableIndex].poloidal;

  toroidalWinding  = mergedWindingPairStats[drawableIndex].toroidal;
  poloidalWinding  = mergedWindingPairStats[drawableIndex].poloidal;

  toroidalWindingP = mergedWindingPairStats[drawableIndex].toroidal;
  poloidalWindingP = mergedWindingPairStats[drawableIndex].poloidal;

  float local_safetyFactor = (float) toroidalWinding / (float) poloidalWinding;

  // The windingGCD only has meaning when islands are present and is
  // related to the number of points in each island.

  // nnodes = windingGCD / resonanceGCD;

  unsigned int windingGCD = GCD( mergedWindingPairStats[drawableIndex].toroidal,
                                 mergedWindingPairStats[drawableIndex].poloidal );


  // Check for primary islands and secondary islands.  NOTE: Even with
  // islands the poloidal resonance can be one as such only check
  // the toroidal resonance
  if( (type == FieldlineProperties::UNKNOWN_TYPE ||
       type == FieldlineProperties::ISLAND_PRIMARY_CHAIN ||
       type == FieldlineProperties::ISLAND_SECONDARY_CHAIN ) &&

        // Acceptance case for an island
      ( (toroidalResonance > 1 /* && poloidalResonance >= 1 */) || // Always true.

        // Also acceptance a 1,1 as an island
        (toroidalWinding == poloidalWinding &&
         toroidalResonance == poloidalResonance) ) )
  {
    // Set the windings to reflect the resonances which is the number
    // of islands.
    if( toroidalWinding == poloidalWinding )
    {
      toroidalWinding  = 1;
      poloidalWinding  = 1;

      toroidalWindingP = toroidalWinding;
      poloidalWindingP = poloidalWinding;
    }
    else
    {
      toroidalWinding  = toroidalResonance;
      poloidalWinding  = poloidalResonance;

      toroidalWindingP = toroidalResonance;
      poloidalWindingP = poloidalResonance;
    }

    // The number of islands is always the toroidal resonance
    // regardless if a primary or secondary island.
    islands      = toroidalResonance;
    islandGroups = toroidalResonance / resonanceGCD;

    // If the resonance GCD is 1 then only one island per group thus a
    // primary island chain.
    if( resonanceGCD == 1 )
    {
      type = FieldlineProperties::ISLAND_PRIMARY_CHAIN;

      if( verboseFlag )
        std::cerr << "Primary resonances = "
                  << toroidalResonance << "," << poloidalResonance << " with "
                  << islands << " islands "
                  << std::endl;
    }

    // The resonance GCD (aka second order resonance) is the number of
    // secondary islands around a primary island.
    else // if( resonanceGCD > 1 )
    {
      type = FieldlineProperties::ISLAND_SECONDARY_CHAIN;

      if( verboseFlag )
        std::cerr << "Secondary resonances = "
                  << toroidalResonance << "," << poloidalResonance << " with "
                  << islands << " islands "
                  << "(" << islands / resonanceGCD << " groups with "
                  << resonanceGCD << " islands each)"
                  << std::endl;
    }

    // When the drawable winding pair is the resonance then potentially
    // in a chaotic area as no other better solutions are drawable.
    if( toroidalResonance == mergedWindingPairStats[drawableIndex].toroidal &&
        poloidalResonance == mergedWindingPairStats[drawableIndex].poloidal )
    {

      if( toroidalResonance == poloidalResonance )
        nnodes = toroidalResonance;

      // The best guestimate of number of nodes will be the GCD of the
      // best winding pair. If there is a tie it does not matter as it
      // is a guestimate.

      // Note: when the island is intering the chaotic zone the nnodes
      // will not be stable between two tracings of the boundary.
      else
        nnodes = GCD( mergedWindingPairStats[0].toroidal,
                      mergedWindingPairStats[0].poloidal ) / resonanceGCD;

      // Less than the maximum number punctures allowed so add more
      // puncture points.
      if( poloidal_puncture_pts.size() < fi.maxPunctures )
      {
        // Get the number of nodes per island.
        unsigned int nodes = (unsigned int)poloidal_puncture_pts.size() / islands / 2;
        
        // Add two more puncture points per island within and island.
        if( type == FieldlineProperties::ISLAND_SECONDARY_CHAIN )
          nPuncturesNeeded = (nodes + 2) * islands * 2;
        // Add five more puncture points per island.
        else //if( type == FieldlineProperties::ISLAND_PRIMARY_CHAIN )
          nPuncturesNeeded = (nodes + 5) * islands * 2;

        analysisState = FieldlineProperties::ADDING_POINTS;

        if( verboseFlag )
          std::cerr << std::endl
                    << "Island: Not enough puncture points yet; "
                    << "need " << nPuncturesNeeded << " "
                    << "have " << poloidal_puncture_pts.size() << " "
                    << "asking for " << nPuncturesNeeded << " puncture points"
                    << std::endl;
      }

      // Reached the user set maximum number of punctures allowed.
      else //if( poloidal_puncture_pts.size() >= fi.maxPunctures )
      {
        analysisState = FieldlineProperties::TERMINATED;

        if( verboseFlag )
          std::cerr << "Potentially within the chaotic regime." << std::endl;
      }
    }

    // When the correct number of points is found the winding pair
    // will be drawable and the winding GCD divided by the resonance GCD
    // will be the number of points in the cross section of each
    // island within an island.

    // Note: for a simple island chain the resonance GCD will be 1.
    else // Possible solution.
    {
      nnodes = windingGCD / resonanceGCD;

      unsigned int nnodesPlus1 = nnodes + 1;

      // For a 1:1 island the nnodes will get stuck in a local minimum
      // of 1,2, or 3. At least that is our observation. As such, add
      // points to see if the analysis can get out of the local
      // minimum.
      if( type == FieldlineProperties::ISLAND_PRIMARY_CHAIN &&
          toroidalWinding == poloidalWinding &&
          nnodes <= 5 )
      {
        nnodes = (unsigned int)poloidal_puncture_pts.size() / toroidalWinding / 2;

        nPuncturesNeeded = poloidal_puncture_pts.size() * 1.05;

        if( nPuncturesNeeded - poloidal_puncture_pts.size() < 5 )
          nPuncturesNeeded = (unsigned int)poloidal_puncture_pts.size() + 5;

        analysisState = FieldlineProperties::ADDING_POINTS;

        if( verboseFlag )
          std::cerr << "Island: Local minimum 0, not enough puncture points; "
                    << "need " << nPuncturesNeeded << " "
                    << "have " << poloidal_puncture_pts.size() << " "
                    << "asking for " << nPuncturesNeeded << " puncture points"
                    << std::endl;
      }

      // Try to get at least four points per island which doubled so
      // the analysis has enough points to work on.
      else if( nnodes < 4 )
      {
        if( 4 * islands * 2 > poloidal_puncture_pts.size() )
        {
          nPuncturesNeeded = 4 * islands * 2;
          analysisState = FieldlineProperties::ADDING_POINTS;

          if( verboseFlag )
            std::cerr << "Adding puncture points for four points per island; "
                      << "have " << poloidal_puncture_pts.size() << " "
                      << "asking for " << nPuncturesNeeded << " puncture points"
                      << std::endl;

        }
        else
        {
          nPuncturesNeeded = 0;

          analysisState = FieldlineProperties::TERMINATED;
        }
      }

      // Get enough points so that the full toroidal and poloidal
      // periods can be analyzed.
      else if( poloidal_puncture_pts.size() < 2*(toroidalResonance*nnodesPlus1) ||
               ridgeline_points.size()      < 2*(poloidalResonance*nnodesPlus1) )
      {
        // For the toroidal period allow for one more possible period
        // to be exaimed.
        if( nPuncturesNeeded < 2.0 * (toroidalResonance*nnodesPlus1) )
        {
          nPuncturesNeeded = 2.0 * (toroidalResonance*nnodesPlus1);
          analysisState = FieldlineProperties::ADDING_POINTS;
        
          if( verboseFlag )
            std::cerr << "Analysis: Not enough puncture points; "
                      << "need " << nPuncturesNeeded << " "
                      << "have " << poloidal_puncture_pts.size() << " "
                      << "asking for " << nPuncturesNeeded << " puncture points"
                      << std::endl;
        }
        
        // For the poloidal period allow for one more possible period
        // to be exaimed which requires two more additional toroidal
        // punctures.
        if( nPuncturesNeeded <
            2.0 * (poloidalResonance*(nnodes+2)) * local_safetyFactor + 0.5)
        {
          nPuncturesNeeded =
            2.0 * (poloidalResonance*(nnodes+2)) * local_safetyFactor + 0.5;
          analysisState = FieldlineProperties::ADDING_POINTS;

          if( verboseFlag )
            std::cerr << "Analysis: Not enough ridgeline points; "
                      << "need " << 2*(poloidalResonance*nnodes+2) << " "
                      << "have " << ridgeline_points.size() << " "
                      << "asking for " << nPuncturesNeeded << " puncture points"
                      << std::endl;
        }
      }

      // Check to see if the fieldline is periodic. I.e. on a rational
      // surface.  If within "delta" of the distance the fieldline is
      // probably on a rational surface.

      // NOTE: Define the rational based on the lowest order of teh
      // surface. Thus divide by the GCD.
      if( rationalCheck( poloidal_puncture_pts,
                         toroidalWinding / GCD(toroidalWinding, poloidalWinding ),
                         nnodes, rationalSurfaceTolerance ) ) 
      {
        type = FieldlineProperties::O_POINT;
        analysisState = FieldlineProperties::COMPLETED;

        if( verboseFlag )
          std::cerr << "Appears to be an O point "
                    << rationalSurfaceTolerance << std::endl;

        fi.rationalSurfaceTolerance = rationalSurfaceTolerance;

        if( detectIslandCenters )
        {
          islandSeedPts.clear();

          for( unsigned int i=0; i<islands; ++i )
          {
            islandSeedPts.push_back( poloidal_puncture_pts[i] );
          }

          analysisState = FieldlineProperties::ADD_BOUNDARY_POINT;
        }
      }

      if( type != FieldlineProperties::O_POINT &&
          nPuncturesNeeded == 0 )
        analysisState = FieldlineProperties::COMPLETED;
    }

    windingGroupOffset = Blankinship( toroidalWinding, poloidalWinding );

    // Surface or island so check for the secondary axis case. To date
    // we have seen secondary axis in 1,1 surfaces and islands and in
    // higher order islands. But not in higher order surfaces.
    if( analysisState == FieldlineProperties::COMPLETED &&
        //        toroidalWinding == poloidalWinding &&
        (type == FieldlineProperties::ISLAND_PRIMARY_CHAIN ||
         type == FieldlineProperties::ISLAND_SECONDARY_CHAIN) )
    {
      // QUESTION - is it possible to have secondary islands within
      // the 1,1 that have an secondary axis. YES
      if( toroidalWinding == poloidalWinding )
      {
        if( toroidalWinding > 1 ) 
          toroidalWindingP = toroidalWinding;
        else
          toroidalWindingP = nnodes;
      }
      else
      {
        toroidalWindingP = nnodes;
      }

      if( verboseFlag )
        std::cerr << "Checking for a secondary axis "
                  << "using a base toroidal winding of "
                  << toroidalWindingP << std::endl; 
      
      std::vector< std::pair< Point, unsigned int > > hullPts;
      
      hullPts.resize( toroidalWindingP+1 );

      // Get the points from the first island to be used for creating
      // the convex hull.
      for(unsigned int i=0, j=0; i<toroidalWindingP*islands; i+=islands, ++j )
      {
        hullPts[j] =
          std::pair< Point, unsigned int >( poloidal_puncture_pts[i], j );
      }
      
      // input - starting index, output - number of points on the hull.
      unsigned int nHullPts = 0;
      
      convexHull( hullPts, nHullPts, toroidalWindingP, helicity );

      if( nHullPts != toroidalWindingP )
      {
        if( verboseFlag )
          std::cerr << "The surface does not have a convex hull, "
                    << toroidalWindingP-nHullPts << " point(s) are missing."
                    << std::endl; 
        
        convexHull( hullPts, nHullPts, toroidalWindingP, -helicity );

        if( nHullPts != toroidalWindingP )
          if( verboseFlag )
            std::cerr << "The second surface does not have a convex hull, "
                      << toroidalWindingP-nHullPts << " point(s) are missing."
                      << std::endl; 
      }
      
      std::map< unsigned int, unsigned int > offsets;
      std::map< unsigned int, unsigned int >::iterator ic;
      
      unsigned int offset = 0; ///TODO: check on fix for uninitialized warning
      
      // Find all the differences and count each one.
      for(unsigned int i=0; i<nHullPts; ++i )
      {
        // Index of the neighboring point.
        unsigned int i1 = (i+1) % toroidalWindingP;

        // Index offset from one puncture point to it's neighbor.
        offset = (hullPts[i1].second - hullPts[i].second + toroidalWindingP) %
          toroidalWindingP;
        
        // Find this offset in the list.
        ic = offsets.find( offset );
        
        // Not found, new offset.
        if( ic == offsets.end() )
          offsets.insert( std::pair< unsigned int, unsigned int >(offset, 1) );
        // Found this offset, increment the count.
        else
          (*ic).second++;
      }

      if( offsets.size() != 1 )
        if( verboseFlag )
          std::cerr << "Multiple offsets  ";

      // Find the offset that occurs most often.
      unsigned int nMatches = 0;
      
      ic = offsets.begin();
      
      while( ic != offsets.end() )
      {
        if( offsets.size() > 1 )
          if( verboseFlag && (*ic).second > 1 )
            std::cerr << (*ic).first << " (" << (*ic).second << ")  ";

        // Keep the offset with most occurances.
        if( nMatches < (*ic).second )
        {
          offset = (*ic).first;
          nMatches = (*ic).second;
        }

        ++ic;
      }

      if( offsets.size() > 1 )
        if( verboseFlag )
          std::cerr << std::endl;

      // Secondary rotation around the nonaxisymmetric island.
      if( offset != 1 && offset != toroidalWindingP-1 )
      {
        toroidalWindingP = toroidalWindingP;
        poloidalWindingP = Blankinship( toroidalWindingP, offset ) + toroidalWindingP;

        windingGroupOffset = offset;

//        nnodes = 1;

        if(type == FieldlineProperties::ISLAND_PRIMARY_CHAIN)
          type = FieldlineProperties::ISLAND_PRIMARY_SECONDARY_AXIS;
        
        else if(type == FieldlineProperties::ISLAND_SECONDARY_CHAIN)
          type = FieldlineProperties::ISLAND_SECONDARY_SECONDARY_AXIS;
        
        if( verboseFlag )
          std::cerr << "Secondary poloidal rotation  "
                    << toroidalWindingP << "," << poloidalWindingP << " ("
                    << ((float) toroidalWindingP / (float) poloidalWindingP) << ")  "
                    << "with offset " << offset << std::endl;
      }
      else
      {
        toroidalWindingP = toroidalWinding;
        poloidalWindingP = poloidalWinding;
      }

      if( detectIslandCenters )
      {
        std::vector< Point > axis;

        if( searchState == FieldlineProperties::NO_SEARCH ||
            searchState == FieldlineProperties::ISLAND_O_POINT )
        {
          if( type == FieldlineProperties::ISLAND_PRIMARY_CHAIN ||
              type == FieldlineProperties::ISLAND_SECONDARY_CHAIN )
            findIslandCenters( poloidal_puncture_pts,
                               islands,
                               toroidalWinding,
                               nnodes,
                               (unsigned int)poloidal_puncture_pts.size(),
                               islandSeedPts,
                               axis );

          else if( type == FieldlineProperties::ISLAND_PRIMARY_SECONDARY_AXIS ||
                   type == FieldlineProperties::ISLAND_SECONDARY_SECONDARY_AXIS )
            findIslandCenters( poloidal_puncture_pts,
                               islands,
                               islands * windingGroupOffset, //offset
                               toroidalWindingP,             //nnodes
                               islands * toroidalWindingP,   //modulo
                               islandSeedPts,
                               axis );
          
          if( islandSeedPts.empty() )
            analysisState = FieldlineProperties::COMPLETED;
          else
          {
            analysisState = FieldlineProperties::ADD_O_POINTS;

            searchNormal = axis[0];
            searchNormal.normalize();
            searchDelta = 2.0 * axis[1].length();
            
            lastSeedPoint = islandSeedPts[0];
            
            std::cerr << "LINE " << __LINE__ << "  "
                      << islandSeedPts[0] << "  "
                      << axis[0] << "  "
                      << axis[0].length() << "  "
                      << axis[1] << "  "
                      << axis[1].length() << "  "
                      << lastSeedPoint << "  "
                      << searchNormal << "  "
                      << searchDelta << std::endl;
          }
        }
        else if( searchState == FieldlineProperties::ISLAND_BOUNDARY_SEARCH )
        {
          analysisState = FieldlineProperties::ADD_BOUNDARY_POINT;
        }
        else
          analysisState = FieldlineProperties::COMPLETED;
      }
    }
    else
    {
      toroidalWindingP = toroidalWinding;
      poloidalWindingP = poloidalWinding;
    }
  }


  // Check to see if the fieldline is periodic. I.e. on a rational
  // surface.  If within "delta" of the distance the fieldline is
  // probably on a rational surface.
  else if( rationalCheck( poloidal_puncture_pts,
                          toroidalWinding / GCD(toroidalWinding, poloidalWinding ),
                          nnodes, rationalSurfaceTolerance ) ) 
  {
    if( GCD(toroidalWinding, poloidalWinding ) == 1 )
    {
      type = FieldlineProperties::RATIONAL;
      analysisState = FieldlineProperties::COMPLETED;

      islands = 0;
      islandGroups = 0;

      fi.rationalSurfaceTolerance = rationalSurfaceTolerance;

      if( verboseFlag )
        std::cerr << "Appears to be a rational surface "
                  << rationalSurfaceTolerance << std::endl;
    }
    else // if( GCD(toroidalWinding, poloidalWinding ) > 1 )
    {
      type = FieldlineProperties::O_POINT;
      analysisState = FieldlineProperties::COMPLETED;

      islands = toroidalWinding;
      islandGroups = 1;

      fi.rationalSurfaceTolerance = rationalSurfaceTolerance;

      if( verboseFlag )
        std::cerr << "Appears to be an O point "
                  << rationalSurfaceTolerance << std::endl;

      if( detectIslandCenters )
      {
        islandSeedPts.clear();
        
        for( unsigned int i=0; i<islands; ++i )
        {
          islandSeedPts.push_back( poloidal_puncture_pts[i] );
        }
        
        analysisState = FieldlineProperties::ADD_BOUNDARY_POINT;
      }
    }

    windingGroupOffset = Blankinship( toroidalWinding, poloidalWinding );

    if( windingGroupOffset == 0 &&
        toroidalWinding < poloidalWinding &&
        GCD(toroidalWinding, poloidalWinding) )
      windingGroupOffset = 1;
  }

  // At this point assume the surface is irrational and thus a flux
  // surface.
  else
  {
    type = FieldlineProperties::FLUX_SURFACE;
    islands = 0;
    islandGroups = 0;

    windingGroupOffset = Blankinship( toroidalWinding, poloidalWinding );

    // Get the direction based on the first two points in a group.
    Vector v0 = (Vector) poloidal_puncture_pts[toroidalWinding] -
      (Vector) poloidal_puncture_pts[0];

    // Get the direction based on the first points from adjacent groups.
    Vector v1 = (Vector) poloidal_puncture_pts[windingGroupOffset] -
      (Vector) poloidal_puncture_pts[0];
      
    // If the windingGroupOffset and point ordering are opposite in
    // directions then the previous group is the
    // windingGroupOffset. Otherwise is they have the same direction
    // then toroidalWinding-windingGroupOffset is the previous group.
    int offsetDir = (Dot( v0, v1 ) > 0.0 ? 1 : -1);
      
    std::vector< unsigned int > nodes(toroidalWinding);

    nnodes = 0;

    // Find the first point from another group that overlaps the first
    // group. This only works if there is an overlap between groups.
    for( unsigned int i=0; i<toroidalWinding; ++i ) 
    {
      nodes[i] = (unsigned int)poloidal_puncture_pts.size() / toroidalWinding;

      // The next group
      unsigned int j =
        (i+offsetDir*windingGroupOffset + toroidalWinding) % toroidalWinding;
        
      Vector firstPoint = (Vector) poloidal_puncture_pts[j];
      Vector nextPoint  = (Vector) poloidal_puncture_pts[j+toroidalWinding];
        
      Vector v0 = nextPoint - firstPoint;

      for( unsigned int k=i; k<poloidal_puncture_pts.size(); k+=toroidalWinding ) 
      {
        Vector  testPoint = (Vector) poloidal_puncture_pts[k];
        
        Vector v1 = testPoint - firstPoint;
        
        if( (poloidalWinding > 2 || v1.length() / v0.length() < 5.0)
            && Dot( v0, v1 ) > 0.0 )
        {
          analysisState = FieldlineProperties::COMPLETED;
          
          nodes[i] = k / toroidalWinding;

          break;
        }
      }

      if( nodes[i] < 1 )
        nodes[i] = (unsigned int)poloidal_puncture_pts.size() / toroidalWinding;

      nnodes += nodes[i];
    }

    // Get the average number of nodes.
    nnodes /= (float) toroidalWinding;

    unsigned int cc = 0;
    unsigned int new_nnodes = 0;
      
    for( unsigned int i=0; i<toroidalWinding; i++ )
    {
      if( nodes[i] < nnodes - 1 || nnodes + 1 < nodes[i] )
      {
        cc++;
      }
      else
        new_nnodes += nodes[i];
    }
    
    if( cc == 1 )
    {
      nnodes = new_nnodes / (float) (toroidalWinding-cc);
    }
      
    for( unsigned int i=0; i<toroidalWinding; i++ )
    {
      if( nodes[i] < nnodes - 1 || nnodes + 1 < nodes[i] )
      {
        if( verboseFlag )
        {
          std::cerr << "Appears to be a flux surface but not self consistent, "
                    << "average number of nodes  " << nnodes << "   (  ";
          
          for( unsigned int i=0; i<toroidalWinding; i++ )
            std::cerr << nodes[i] << "  ";

          std::cerr << ")" << std::endl;
        }

        break;
      }
    }

    if( analysisState != FieldlineProperties::COMPLETED )
    {
      unsigned int additionalPts = 0;

      // Get the direction based on the first two points in a group.
      Vector v0 = (Vector) poloidal_puncture_pts[toroidalWinding] -
        (Vector) poloidal_puncture_pts[0];

      // Get the direction based on the first points from adjacent groups.
      Vector v1 = (Vector) poloidal_puncture_pts[windingGroupOffset] -
        (Vector) poloidal_puncture_pts[0];

      // If the windingGroupOffset and point ordering are opposite in
      // directions then the previous group is the
      // windingGroupOffset. Otherwise is they have the same direction
      // then toroidalWinding-windingGroupOffset is the previous
      // group.
      int offsetDir = (Dot( v0, v1 ) > 0.0 ? 1 : -1);

      for( unsigned int i=0; i<toroidalWinding; ++i ) 
      {
        // The next group
        unsigned int j =
          (i+offsetDir*windingGroupOffset + toroidalWinding) % toroidalWinding;

        Vector firstPoint = (Vector) poloidal_puncture_pts[j];
          
        Vector  lastPoint =
          (Vector) poloidal_puncture_pts[i+(nnodes-1)*toroidalWinding];
        Vector  prevPoint =
          (Vector) poloidal_puncture_pts[i+(nnodes-2)*toroidalWinding];
          
        unsigned int needPts = ( (firstPoint-lastPoint).length() /
                                 (prevPoint-lastPoint).length() + 0.5 );
          
        if( additionalPts < needPts )
          additionalPts = needPts;
      }

      if( additionalPts )
      {
        if( nPuncturesNeeded == 0 )
        {
          nPuncturesNeeded = (nnodes+additionalPts) * toroidalWinding + 1;
          analysisState = FieldlineProperties::ADDING_POINTS;
        
          if( verboseFlag )
            std::cerr << "Too few puncture points, at least "
                      << (nnodes+additionalPts) * toroidalWinding
                      << " are needed to complete the boundary."
                      << std::endl;
        }
      }
      else
        analysisState = FieldlineProperties::COMPLETED;
    }
  }

  // The user has set the toroidal winding get the poloidal winding
  // based on the data found.
  if( overrideToroidalWinding )
  {
    toroidalWinding = overrideToroidalWinding;

    if( overridePoloidalWinding )
    {
      poloidalWinding  = overridePoloidalWinding;
      poloidalWindingP = overridePoloidalWinding;

      double local_safetyFactor =
        (double) toroidalWinding / (double) poloidalWinding;

      if( verboseFlag )
        std::cerr << "overriding" << std::endl
                  << "**using**   winding pair "
                  << toroidalWinding << ","
                  << poloidalWinding << "  ("
                  << local_safetyFactor << " - "
                  << fabs(safetyFactor - local_safetyFactor) << ")  "
                  << std::endl;
    }
    else
    {
      unsigned int nsets = (unsigned int)poloidalWindingCounts.size();

      std::map< unsigned int, unsigned int > differenceCount;
      std::map< unsigned int, unsigned int >::iterator ic;
      
      // Find all the differences and count each one.
      for( unsigned int i=0; i<nsets-toroidalWinding; ++i)
      {
        // Get the poloidal winding between two counts.
        poloidalWinding = poloidalWindingP =
          poloidalWindingCounts[i+toroidalWinding] - poloidalWindingCounts[i];
        
        // Find this difference in the list.
        ic = differenceCount.find( poloidalWinding );
        
        // Not found, new difference.
        if( ic == differenceCount.end() )
          differenceCount.insert( std::pair<int, int>( poloidalWinding, 1) );
        // Found this difference, increment the count.
        else (*ic).second++;
      }
      
      // Find the difference that occurs most often.
      unsigned int nMatches = 0;
      
      ic = differenceCount.begin();
      
      while( ic != differenceCount.end() )
      {
        if( nMatches < (*ic).second )
        {
          poloidalWinding  = (*ic).first;
          poloidalWindingP = (*ic).first;
          nMatches = (*ic).second;
        }
        
        ++ic;
      }

      double local_safetyFactor =
        (double) toroidalWinding / (double) poloidalWinding;

      double consistency =
        (double) nMatches / (double) (nsets-toroidalWinding);

      if( verboseFlag )
        std::cerr << "overriding" << std::endl
                  << "**using**   winding pair "
                  << toroidalWinding << ","
                  << poloidalWinding << "  ("
                  << local_safetyFactor << " - "
                  << fabs(safetyFactor - local_safetyFactor) << ")  "
                  << "consistency "
                  << 100.0 * consistency
                  << " %)" << std::endl;
    }

    if( poloidalWinding == poloidalWindingP )
      windingGroupOffset = Blankinship( toroidalWinding, poloidalWinding );

    nnodes = (unsigned int)poloidal_puncture_pts.size() / toroidalWinding;

    analysisState = FieldlineProperties::OVERRIDE;
  }

  // A catch all in case the somehow the number of puncture points
  // exceeds the maximum allowed.
  if( poloidal_puncture_pts.size() > fi.maxPunctures )
  {
    if( verboseFlag )
      std::cerr << "FORCE TERMINATING  " << nPuncturesNeeded << std::endl;

    nPuncturesNeeded = 0;

    analysisState = FieldlineProperties::TERMINATED;
  }

  // Record the analysis.
  fi.analysisState = analysisState;
  fi.searchState = searchState;

  fi.type = type;

  fi.safetyFactor = safetyFactor;
  fi.toroidalWinding = toroidalWinding * toroidalBase;
  fi.poloidalWinding = poloidalWinding * toroidalBase;
  fi.toroidalPeriod = toroidalPeriod;
  fi.poloidalPeriod = poloidalPeriod;
  fi.toroidalWindingP = toroidalWindingP;
  fi.poloidalWindingP = poloidalWindingP;
  fi.toroidalResonance = toroidalResonance;
  fi.poloidalResonance = poloidalResonance;
  fi.windingPairs = windingPairs;
  fi.topWindingPairs = topWindingPairs;
  fi.windingGroupOffset = windingGroupOffset;
  fi.islands = islands;
  fi.islandGroups = islandGroups;
  fi.nnodes = nnodes;
  fi.nPuncturesNeeded = nPuncturesNeeded;

  fi.searchDelta   = searchDelta;
  fi.searchNormal  = searchNormal;
  fi.lastSeedPoint = lastSeedPoint;

  fi.seedPoints.clear();

  // For an island O Point search save the geometric centers.
  if( !(islandSeedPts.empty()) )
  {
    fi.seedPoints.resize( islandSeedPts.size() );

    for( unsigned int i=0; i<islandSeedPts.size(); ++i )
    {
      fi.seedPoints[i] = islandSeedPts[i];
    }
  }

//   FieldlineProperties tmp;

//   if( analysisState == FieldlineProperties::COMPLETED )
//     fieldlineProperties2( ptList, tmp );
}


void
FieldlineLib::fieldlineProperties2( std::vector< Point > &ptList,
                                    double rationalSurfaceFactor,
                                    FieldlineProperties &fi )
{
  /*
 std::vector< Point > poloidal_puncture_pts;
 std::vector< Point > ridgeline_points;

 std::vector< double > rotationalSums;
 std::vector< unsigned int > poloidalWindingCounts;

 float delta = 0.0;

  getPunctures( ptList, Vector(0,1,0), poloidal_puncture_pts );

  if( ptList.empty() && poloidal_puncture_pts.empty() )
  {
    fi.type = FieldlineProperties::UNKNOWN_TYPE;
    fi.analysisState = FieldlineProperties::UNKNOWN_ANALYSIS;

    fi.safetyFactor = 0;
    fi.toroidalWinding = 0;
    fi.poloidalWinding = 0;
    fi.poloidalWindingP = 0;
    fi.toroidalResonance = 0;
    fi.poloidalResonance = 0;
    fi.windingPairs.clear();
    fi.windingGroupOffset = 0;
    fi.islands = 0;
    fi.islandGroups = 0;
    fi.nnodes  = 0;

    fi.nPuncturesNeeded  = 0;

    fi.seedPoints.clear();

    return;
  }

  // Start the analysis.
  FieldlineProperties::FieldlineType type = FieldlineProperties::UNKNOWN_TYPE;
  FieldlineProperties::AnalysisState analysisState = FieldlineProperties::UNKNOWN_ANALYSIS;

  unsigned int toroidalWinding = 0;
  unsigned int toroidalResonance = 1;

  unsigned int nnodes = 0;
  unsigned int islands = 0, islandGroups = 0;
  unsigned int nPuncturesNeeded = 0;

  unsigned int toroidalWindingMax = poloidal_puncture_pts.size() / 2;

  // Find the best toroidal periodicity. For a flux surface the period
  // will be the toroidal winding number. For an island chain the
  // period will be the toroidal winding number times the number of
  // nodes.
  if( verboseFlag )
    std::cerr << "Toroidal Winding via "
         << poloidal_puncture_pts.size() << "  "
         << "poloidal punctures, "
         << "max period " << toroidalWindingMax
         << std::endl;

 std::vector< std::pair< unsigned int, double > > toroidalStats;

  periodicityStats( poloidal_puncture_pts, toroidalStats,
                    toroidalWindingMax, 2 );


  int drawableRank  = -1;
  int drawableIndex = -1;
  std::vector< unsigned int > drawableIndexs;

  for( unsigned int i=0; i<toroidalStats.size(); ++i )
  {
    bool drawable = false;

    for( unsigned int j=1; j<toroidalStats[i].first; ++j )
    {
      if( IntersectCheck( poloidal_puncture_pts, toroidalStats[i].first, j ) )
      {
        drawable = true;
        break;
      }
    }

    if( verboseFlag && (drawable || i < 10) )
      std::cerr << "Final "
           << (drawable ? "Drawable " : "Rejected ") 
           << "toroidal winding :  " 
           << toroidalStats[i].first << std::endl;

    if( (drawableRank == -1 ||
         drawableRank == toroidalStats[i].second) &&

        // Keep only those that are drawable.
        drawable )
    {
      drawableIndexs.push_back( i );

      // If a tie take the higher order winding.
      if( toroidalWinding < toroidalStats[i].first )
      {
        toroidalWinding = toroidalStats[i].first;

        drawableRank = toroidalStats[i].first;
        drawableIndex = i;
      }
    }
  }

  if( drawableIndex == -1 )
  {
    toroidalWinding = 0;

    if( verboseFlag )
      std::cerr << "Garbage matches adding more points" << std::endl;

    fi.nPuncturesNeeded = poloidal_puncture_pts.size() * 1.25;

    return;
  }

  if( toroidalStats.size() == 1 )
  {
    toroidalResonance = 1;
  }
  else
  {
    // The base values are the primary resonances if present.
    toroidalResonance = ResonanceCheck( toroidalStats, 1, 3 );
  }

  if( verboseFlag )
    std::cerr << "Toroidal Winding " << toroidalWinding << "  "
         << "Toroial Resonances = " << toroidalResonance << "  "
         << std::endl;

  // Check for islands and islands around islands.
  // NOTE: Even with islands the poloidalFirstResonance can be one as such
  // only check the toroidalFirstResonance
  if( (type == FieldlineProperties::UNKNOWN_TYPE ||
       type == FieldlineProperties::ISLAND_PRIMARY_CHAIN ||
       type == FieldlineProperties::ISLAND_SECONDARY_CHAIN ) &&

      toroidalResonance > 1 ) // && poloidalResonance >= 1 ) // Always true.
  {
    // Set the windings to reflect the resonances which is the number
    // of islands.
    toroidalWinding = toroidalResonance;

    // The number of islands is always the toroidal resonance.
    islands = toroidalResonance;

    // If the resonance GCD is 1 then only one island per group thus a
    // simple island chain.
    if( isPrime( toroidalResonance ) == 1 )
    {
      type = FieldlineProperties::ISLAND_PRIMARY_CHAIN;

      if( verboseFlag )
        std::cerr << "Primary resonances = "
             << toroidalResonance << " with "
             << islands << " islands "
             << std::endl;
    }
    // The resonance GCD (aka secondary resonance) is the number of
    // smaller islands around an island.
    else
    {
      type = FieldlineProperties::ISLAND_SECONDARY_CHAIN;

      if( verboseFlag )
        std::cerr << "Secondary resonances = "
             << toroidalResonance << " with "
             << islands << " islands "
             << std::endl;
    }

    // When the drawable winding pair is the resonance then potentially
    // in a chaotic area as no other better solutions are drawable.
    if( toroidalResonance == toroidalStats[drawableIndex].first )
    {
      // The best guestimate of number of nodes will be the GCD of the
      // best winding pair. If there is a tie it does not matter as it
      // is a guestimate.

      // Note: when the island is intering the chaotic zone the nnodes
      // will not be stable between two tracings of the boundary.
      nnodes = toroidalStats[0].first / toroidalResonance;

      // Less than the maximum number punctures allowed so add more
      // puncture points.
      if( poloidal_puncture_pts.size() < fi.maxPunctures )
      {
        analysisState = FieldlineProperties::ADDING_POINTS;

        // Get the number of nodes per island.
        unsigned int nodes = poloidal_puncture_pts.size() / islands / 2;
        
        // Add two more puncture points per island for a secondary
        // chain.
        if( type == FieldlineProperties::ISLAND_SECONDARY_CHAIN )
          nPuncturesNeeded = (nodes + 2) * islands * 2;
        // Add five more puncture points per island for a primary
        // chain.
        else //if( type == FieldlineProperties::ISLAND_PRIMARY_CHAIN )
          nPuncturesNeeded = (nodes + 5) * islands * 2;

        if( verboseFlag )
          std::cerr << std::endl
               << "Island: Not enough puncture points yet; "
               << "need " << nPuncturesNeeded << " "
               << "have " << poloidal_puncture_pts.size() << " "
               << "asking for " << nPuncturesNeeded << " puncture points"
               << std::endl;
      }

      // Reached the user set maximum number of punctures allowed.
      else //if( poloidal_puncture_pts.size() >= fi.maxPunctures )
      {
        analysisState = FieldlineProperties::TERMINATED;

        if( verboseFlag )
          std::cerr << "Potentially within the chaotic regime." << std::endl;
      }
    }

    // When the correct number of points is found the winding pair
    // will be drawable and the winding GCD divided by the resonance GCD
    // will be the number of points in the cross section of each
    // island within an island.

    // Note: for a simple island chain the resonance GCD will 1.
    else // Possible solution.
    {
      nnodes = toroidalStats[0].first / toroidalResonance;

      unsigned int nnodesPlus1 = nnodes + 1;

      // Try to get at least four points per island.
      if( nnodes < 4 )
      {
        analysisState = FieldlineProperties::ADDING_POINTS;

        nPuncturesNeeded = 4 * islands * 2;

        if( nPuncturesNeeded <= poloidal_puncture_pts.size() )
        {
          nPuncturesNeeded = 0;
          analysisState = FieldlineProperties::TERMINATED;
        }

        else
        {
          if( verboseFlag )
            std::cerr << "Adding puncture points for four points per island; "
                 << "have " << poloidal_puncture_pts.size() << " "
                 << "asking for " << nPuncturesNeeded << " puncture points"
                 << std::endl;
        }
      }

      // For a 1:1 island the nnodes will get stuck in a local minimum
      // of 1,2, or 3. At least that is our observation. As such, add
      // points to see if the analysis can get out of the local
      // minimum.
      else if( type == FieldlineProperties::ISLAND_PRIMARY_CHAIN &&
               toroidalWinding == 1 && nnodes <= 5 )
      {
        nnodes = poloidal_puncture_pts.size() / toroidalWinding / 2;

        analysisState = FieldlineProperties::ADDING_POINTS;

        nPuncturesNeeded = poloidal_puncture_pts.size() * 1.05;

        if( nPuncturesNeeded - poloidal_puncture_pts.size() < 5 )
          nPuncturesNeeded = poloidal_puncture_pts.size() + 5;

        if( verboseFlag )
          std::cerr << "Island:: Local minimum 1, not enough puncture points; "
               << "need " << nPuncturesNeeded << " "
               << "have " << poloidal_puncture_pts.size() << " "
               << "asking for " << nPuncturesNeeded << " puncture points"
               << std::endl;
      }

      // Get enough points so that the full toroidal and poloidal
      // periods can be analyzed.
      else if( poloidal_puncture_pts.size() <
               2*(toroidalResonance*nnodesPlus1) )
      {
        analysisState = FieldlineProperties::ADDING_POINTS;

        // For the toroidal period allow for one more possible period
        // to be exaimed.
        if( nPuncturesNeeded < 2.0 * (toroidalResonance*nnodesPlus1) )
        {
          nPuncturesNeeded = 2.0 * (toroidalResonance*nnodesPlus1);
        
          if( verboseFlag )
            std::cerr << "Analysis:: Not enough puncture points; "
                 << "need " << nPuncturesNeeded << " "
                 << "have " << poloidal_puncture_pts.size() << " "
                 << "asking for " << nPuncturesNeeded << " puncture points"
                 << std::endl;
        }
      }

      else if( nPuncturesNeeded == 0 )
      {
        analysisState = FieldlineProperties::COMPLETED;
      }

      // Check to see if the fieldline is periodic. I.e. on a rational
      // surface.  If within "delta" of the distance the fieldline is
      // probably on a rational surface.
      if( rationalCheck( poloidal_puncture_pts,
                         toroidalWinding,
                         nnodes, delta*rationalSurfaceFactor ) ) 
      {
        type = FieldlineProperties::O_POINT;
        analysisState = FieldlineProperties::COMPLETED;
        
        if( verboseFlag )
          std::cerr << "Appears to be an O point "
                    << delta*rationalSurfaceFactor << std::endl;
      }
    }
  }

  // Check to see if the fieldline is periodic. I.e. on a rational
  // surface.  If within "delta" of the distance the fieldline is
  // probably on a rational surface.
  else if( rationalCheck( poloidal_puncture_pts,
                          toroidalWinding,
                          nnodes, delta*rationalSurfaceFactor ) ) 
  {
    type = FieldlineProperties::RATIONAL;
    islands = 0;
    islandGroups = 0;

    analysisState = FieldlineProperties::COMPLETED;

    if( verboseFlag )
      std::cerr << "Appears to be a rational surface "
                << delta*rationalSurfaceFactor << std::endl;
  }

  // At this point assume the surface is irrational.
  else
  {
    type = FieldlineProperties::FLUX_SURFACE;
    islands = 0;
    islandGroups = 0;

    // Get the average number of nodes.
    nnodes = poloidal_puncture_pts.size() / toroidalWinding;
  }

  // Record the analysis.
  fi.analysisState = analysisState;

  fi.type = type;

  fi.toroidalWinding = toroidalWinding;
  fi.toroidalResonance = toroidalResonance;
  fi.islands = islands;
  fi.nnodes  = nnodes;

  if( poloidal_puncture_pts.size() > fi.maxPunctures &&
      (analysisState != FieldlineProperties::COMPLETED ||
       analysisState != FieldlineProperties::TERMINATED) )
  {
    if( verboseFlag )
      std::cerr << "1 TERMINATING  " << nPuncturesNeeded << std::endl;

    analysisState = FieldlineProperties::TERMINATED;
    fi.nPuncturesNeeded = 0;
  }
  else
    fi.nPuncturesNeeded = nPuncturesNeeded + (nPuncturesNeeded ? 1 : 0);

  fi.seedPoints.clear();
  */
}


unsigned int
FieldlineLib::
islandProperties( std::vector< Point > &points,
                  Vector &baseCentroid,
                  unsigned int &startIndex,
                  unsigned int &middleIndex,
                  unsigned int &stopIndex,
                  unsigned int &nodes )
{
  // Get the local centroid for the group.
  Vector localCentroid(0,0,0);

  for( unsigned int i=0; i<points.size(); i++ )
    localCentroid += (Vector) points[i];

  localCentroid /= (float) points.size();

  // Get the principal axes of the points.
  float Ixx = 0.0;
  float Ixz = 0.0;
  float Izz = 0.0;

  double maxDist = 0;

  // Get the moments of intertial for each point. It assumed that
  // everything is in the Y plane as such there the moments of
  // intertial along the Y axis are zero.
  for( unsigned int i=0; i<points.size(); i++ ) {

    Vector vec = (Vector) points[i] - localCentroid;
    
    if( maxDist < vec.length() )
      maxDist = vec.length();
    
    Ixx += vec.z*vec.z;
    Ixz -= vec.x*vec.z;
    Izz += vec.x*vec.x;
  }
  
  // Short cut to the principal axes because the Y moments of
  // intertial are zero.
  float alpha = atan( 2.0 * Ixz / (Ixx - Izz) ) / 2.0;

  //       std::cerr << "PRINCIPAL AXES " << alpha * 180.0 / M_PI << "    "
  //       << Ixx + Ixz * sin(alpha       )/cos(alpha       ) << "    "
  //       << Izz + Ixz * cos(alpha+M_PI/2)/sin(alpha+M_PI/2) << std::endl;

  // Use the principal axes to get an offset from the local
  // centroid which gives a point outside the island.
  baseCentroid = localCentroid;
  
  if( Ixx + Ixz * sin(alpha       )/cos(alpha       ) >
      Izz + Ixz * cos(alpha+M_PI/2)/sin(alpha+M_PI/2) )
    baseCentroid -= Vector(  cos(alpha), 0, sin(alpha) ) * maxDist;
  else
    baseCentroid -= Vector( -sin(alpha), 0, cos(alpha) ) * maxDist;

  // Determine if islands exists. If an island exists there will be
  // both clockwise and counterclockwise sections when compared to the
  // global centroid.
  unsigned int turns = 0;

  startIndex = middleIndex = stopIndex = 0;

  Vector v0 = (Vector) points[0] - baseCentroid;
  Vector v1 = (Vector) points[1] - baseCentroid;

  bool lastCCW = (ccwXZ(v0, v1) == 1);
  v0 = v1;

  for( unsigned int j=2; j<points.size(); j++ ) {

    v1 = (Vector) points[j] - baseCentroid;

    bool CCW = (ccwXZ(v0, v1) == 1);
    v0 = v1;

    if( CCW != lastCCW ) {
      turns++;

      if( turns == 1 )      startIndex  = j - 1;
      else if( turns == 2 ) middleIndex = j - 1;
      else if( turns == 3 ) stopIndex   = j - 1;

      if( turns == 3 )
        break;

      lastCCW = CCW;
    }
  }

  if( turns == 0 ) {
    stopIndex   = (unsigned int)points.size() - 1;
    middleIndex = 0;
    startIndex  = 0;

    nodes = (unsigned int)points.size();
  } else if( turns == 1 ) {
    stopIndex   = (unsigned int)points.size() - 1;
    middleIndex = startIndex;
    startIndex  = 0;

    nodes = (unsigned int)points.size();
  } else {

    if( turns == 3 ) {

      // Check for a negative epsilon.
      double length0 = 0;
      double length1 = 0;

      for( unsigned j=0, k=nodes, l=nodes+1; l<points.size(); j++, k++, l++ ) {
        length0 += ((Vector) points[j] - (Vector) points[k]).length();
        length1 += ((Vector) points[j] - (Vector) points[l]).length();
      }

      if( length0 < length1 )
        stopIndex--;

      if( 2*startIndex == middleIndex + 1 ) {
        // First point is actually the start point.
//      std::cerr << "First point is actually the start point.\n";

        stopIndex   = middleIndex;
        middleIndex = startIndex;
        startIndex  = 0;
      }
    } else if( turns == 2 ) {

      if( 2*startIndex == middleIndex + 1 ) {
        // First point is actually the start point.
//      std::cerr << "First point is actually the start point.\n";

        stopIndex   = middleIndex;
        middleIndex = startIndex;
        startIndex  = 0;
      
        turns = 3;

      } else if( points.size() < 2 * (middleIndex - startIndex) - 1 ) {
        // No possible over lap.
//      std::cerr <<  "islandProperties - No possible over lap.\n";

        stopIndex = startIndex + (unsigned int)points.size() - 1;
      }
    }

    nodes = stopIndex - startIndex + 1;

    for( unsigned int j=middleIndex, k=middleIndex+1;
         j<points.size() && k<points.size();
         ++j, ++k ) {

      // See if the test point is between the first section.
      if( Dot( (Vector) points[0] - (Vector) points[j],
               (Vector) points[1] - (Vector) points[j] )
          < 0.0 ) {
        stopIndex = startIndex + j;
        nodes = j;
        turns = 3;
//      std::cerr << "islandProperties - A point overlaps the first section at " << j-1 << std::endl;
        break;
      }

      // See if the first point is between the test section.
      if( Dot( (Vector) points[j] - (Vector) points[0],
               (Vector) points[k] - (Vector) points[0] )
          < 0.0 ) {
        stopIndex = startIndex + j;
        nodes = j;
        turns = 3;
//      std::cerr <<  "islandProperties - First point overlaps another section after " << j-1 << std::endl;
        break;
      }
    }
      
    // No overlap found
    if( turns == 2 ) {
      stopIndex = startIndex + (unsigned int)points.size() - 1;
      nodes = (unsigned int)points.size();
//      std::cerr << "islandProperties - No overlap found\n";
    }
  }

  return turns;
}


unsigned int
FieldlineLib::
surfaceOverlapCheck( std::vector<std::vector< Point > > &bins,
                     unsigned int windingGroupOffset )
{
  unsigned int toroidalWinding = bins.size();

  unsigned int nnodes = (unsigned int) bins[0].size();
  
  // First make sure none of the groups overlap themselves.
  for( unsigned int i=0; i<toroidalWinding; ++i )
  {
    if( nnodes > bins[i].size() )
      nnodes = (unsigned int) bins[i].size();

    for( unsigned int j=2; j<nnodes; j++ )
    {
      Vector v0 = (Vector) bins[i][0] - (Vector) bins[i][j];
      Vector v1 = (Vector) bins[i][1] - (Vector) bins[i][j];
      
      if( Dot( v0, v1 ) < 0.0 )
      {
        nnodes = j;
        break;
      }

      Vector midPt = ((Vector) bins[i][j] + (Vector) bins[i][j-1]) / 2.0;

      v0 = (Vector) bins[i][0] - midPt;
      v1 = (Vector) bins[i][1] - midPt;
      
      if( Dot( v0, v1 ) < 0.0 )
      {
        nnodes = j;
        break;
      }

    }
  }

  if( toroidalWinding == 1 || nnodes == 1 )
    return nnodes;

  // If the offset and point ordering are opposite in directions then
  // the previous group is the windingGroupOffset. Otherwise if they
  // have the same direction then -windingGroupOffset is the previous
  // group.
  Vector intra = (Vector) bins[                 0][1] - (Vector) bins[0][0];
  Vector inter = (Vector) bins[windingGroupOffset][0] - (Vector) bins[0][0];

  int offset = (Dot( intra, inter ) < 0.0) ?
    windingGroupOffset : toroidalWinding-windingGroupOffset;

  // Second make sure none of the groups overlap each other.
  for( unsigned int i=0; i<toroidalWinding; ++i )
  {
    // The previous group
    unsigned int j = (i + offset) % toroidalWinding;

    // Check for a point in the previous group being between the first
    // two points in the current group.
    for( unsigned int k=0; k<nnodes; ++k )
    {
      Vector v0 = (Vector) bins[i][0] - (Vector) bins[j][k];
      Vector v1 = (Vector) bins[i][1] - (Vector) bins[j][k];
      
      if( Dot( v0, v1 ) < 0.0 ) {
        nnodes = k;
        // std::cerr << "adjacent overlap1 " << i << "  " << j << "  " << k
        //        << std::endl;
      }
    }

    // Check for a point in the current group being between two points
    // in the previous group.
    for( unsigned int k=1; k<nnodes; ++k )
    {
      Vector v0 = (Vector) bins[j][k  ] - (Vector) bins[i][0];
      Vector v1 = (Vector) bins[j][k-1] - (Vector) bins[i][0];
      
      if( Dot( v0, v1 ) < 0.0 ) {
        nnodes = k;
        // std::cerr << "adjacent overlap2 " << i << "  " << j << "  " << k
        //        << std::endl;
        break;
      }
    }
  }

  return nnodes;
}


unsigned int
FieldlineLib::
surfaceGroupCheck( std::vector<std::vector< Point > > &bins,
                   unsigned int i,
                   unsigned int j,
                   unsigned int nnodes )
{
  unsigned int nodes = nnodes;

  while( nodes < bins[i].size() ) {
    // Check to see if the first overlapping point is really a
    // fill-in point. This happens because the spacing between
    // toroidalWinding groups varries between groups.
    Vector v0 = (Vector) bins[j][0      ] - (Vector) bins[i][nodes];
    Vector v1 = (Vector) bins[i][nodes-1] - (Vector) bins[i][nodes];
    
    if( Dot( v0, v1 ) < 0.0 )
      nodes++;
    else
      break;
  }

  return nodes;
}


unsigned int
FieldlineLib::
removeOverlap( std::vector< std::vector< std::vector < Point > > > &bins,
               unsigned int windingGroupOffset )
{
  unsigned int nnodes = 0;
  unsigned int nSections = bins.size();

  if( verboseFlag )
    std::cerr << "removeOverlap - "
              << "number of sections " << nSections
              << std::endl;

  
  for( unsigned int s=0; s<nSections; ++s )
  {
    unsigned int toroidalWinding = bins[s].size();

    // Get the centroid for this section.
    unsigned int cc = 0;
    Vector globalCentroid = Vector(0,0,0);
    
    for( unsigned int i=0; i<toroidalWinding; ++i )
    {
      for( unsigned int j=0; j<bins[s][i].size(); ++j )
      {
        globalCentroid += (Vector) bins[s][i][j];
        ++cc;
      }
    }
  
    globalCentroid /= (double) cc;
    
    // This gives the minimal number of nodes for each winding group.
    nnodes = surfaceOverlapCheck( bins[s], windingGroupOffset );
    
    if( nnodes == 0 )
    {
      if( verboseFlag )
        std::cerr << "removeOverlap - "
                  << "Surface properties returned ZERO NODES for surface "
                  << std::endl;

      nnodes = (unsigned int) bins[s][0].size();

      for( unsigned int i=1; i<toroidalWinding; i++ )
        {
        if( nnodes > bins[s][i].size())
          nnodes = (unsigned int) bins[s][i].size();
        }
    }
    
    // If the offset and point ordering are opposite in directions
    // then the next group is the -windingGroupOffset. Otherwise if
    // they have the same direction then windingGroupOffset is the
    // next group.
    Vector intra = ((Vector) bins[s][0][1] -
                    (Vector) bins[s][0][0]);
    Vector inter = ((Vector) bins[s][windingGroupOffset][0] -
                    (Vector) bins[s][0][0]);
    
    int offset = (Dot( intra, inter ) < 0.0) ?
      toroidalWinding-windingGroupOffset : windingGroupOffset;
    
    for( unsigned int i=0; i<toroidalWinding; ++i )
    {
      // The next group
      unsigned int j = (i + offset) % toroidalWinding;
      
      // Add back in any nodes that may not overlap.
      unsigned int nodes = surfaceGroupCheck( bins[s], i, j, nnodes );
      
      // No more than one point should be added.
      if( nodes > nnodes+1 )
      {
        if( verboseFlag )
          std::cerr << "removeOverlap - Surface " << i << "  "
                    << "nnodes mismatch " << nnodes << "  "
                    << nodes << std::endl;
      }
      
      // Erase all of the overlapping points.
      bins[s][i].erase( bins[s][i].begin()+nodes, bins[s][i].end() );
    }
  }
  
  return nnodes;
}


unsigned int
FieldlineLib::
smoothCurve( std::vector< std::vector < Point > > &bins,
             unsigned int &nnodes,
             unsigned int toroidalWinding,
             unsigned int offset,
             unsigned int island )
{
  Vector globalCentroid = Vector(0,0,0);;

  for( unsigned int i=0; i<toroidalWinding; i++ )
    for( unsigned int j=0; j<nnodes; j++ )
      globalCentroid += (Vector) bins[i][j];
  
  globalCentroid /= (toroidalWinding*nnodes);

  unsigned int add = 2;

  if( island ) {

    for( unsigned int i=0; i<toroidalWinding; i++ ) {
//      for( unsigned int s=0; s<add; s++ )
      {
       std::vector< std::pair< Point, unsigned int > > newPts;

        newPts.resize( add*nnodes );

        for( unsigned int j=0; j<add*nnodes; j++ )
          newPts[j] = std::pair< Point, unsigned int > (Point(0,0,0), 0 );
        
        for( unsigned int j=1; j<nnodes-1; j++ ) {

          unsigned int j_1 = (j-1+nnodes) % nnodes;
          unsigned int j1  = (j+1+nnodes) % nnodes;

          Vector v0 = (Vector) bins[i][j1] - (Vector) bins[i][j  ];
          Vector v1 = (Vector) bins[i][j ] - (Vector) bins[i][j_1];

          if( verboseFlag )
            std::cerr << i << "  smooth  " << j_1 << "  "
                      << j << " "  << j1 << "  "
                      << ( v0.length() > v1.length() ?
                           v0.length() / v1.length() :
                           v1.length() / v0.length() ) << std::endl;

          if( Dot( v0, v1 ) > 0 &&
              ( v0.length() > v1.length() ?
                v0.length() / v1.length() :
                v1.length() / v0.length() ) < 10.0 ) {

            Vector center = (Vector) circle( bins[i][j_1],
                                             bins[i][j  ],
                                             bins[i][j1 ] );

            double rad = ((Vector) bins[i][j] - center).length();


            for( unsigned int s=0; s<add; s++ ) {
              Vector midPt = (Vector) bins[i][j_1] +
                (double) (add-s) / (double) (add+1) *
                ((Vector) bins[i][j] - (Vector) bins[i][j_1] );
                

              Vector midVec = midPt - center;

              midVec.normalize();

              newPts[add*j+s].first += center + midVec * rad;
              newPts[add*j+s].second += 1;

              midPt = (Vector) bins[i][j] +
                (double) (add-s) / (double) (add+1) *
                ((Vector) bins[i][j1] - (Vector) bins[i][j] );

              midVec = midPt - center;

              midVec.normalize();

              newPts[add*j1+s].first += center + midVec * rad;
              newPts[add*j1+s].second += 1;
            }
          }
        }

        for( unsigned int j=nnodes-1; j>0; j-- ) {

          for( unsigned int s=0; s<add; s++ ) {

            unsigned int k = add * j + s;

            if( newPts[k].second > 0 ) {
              
              newPts[k].first /= newPts[k].second;
              
//            std::cerr << i << " insert " << j << "  " << newPts[k] << std::endl;
              
              bins[i].insert( bins[i].begin()+j, newPts[k].first );
            }
          }
        }

        for( unsigned int s=0; s<add; s++ ) {

          unsigned int k = add - 1 - s;

          if( newPts[k].second > 0 ) {
              
            newPts[k].first /= newPts[k].second;
              
//            std::cerr << i << " insert " << 0.0<< "  " << newPts[k] << std::endl;
              
            bins[i].push_back( newPts[k].first );
          }
        }
      }
    }

  } else {

    for( unsigned int i=0; i<toroidalWinding; i++ ) {

      if( bins[i].size() < 2 )
        continue;

      // Index of the next toroidalWinding group
      unsigned int j = (i+offset)%toroidalWinding;

      // Insert the first point from the next toroidalWinding so the curve
      // is contiguous.
      bins[i].push_back( bins[j][0] );

      //for( unsigned int s=0; s<add; s++ )
      {
        unsigned int nodes = (unsigned int)bins[i].size();

       std::vector<std::pair< Point, unsigned int > > newPts(add*nodes);

        for( unsigned int j=0; j<add*nodes; j++ )
          newPts[j] = std::pair< Point, unsigned int > (Point(0,0,0), 0 );
        
        for( unsigned int j=1; j<nodes-1; j++ ) {

          unsigned int j_1 = j - 1;
          unsigned int j1  = j + 1;

          Vector v0 = (Vector) bins[i][j1] - (Vector) bins[i][j  ];
          Vector v1 = (Vector) bins[i][j ] - (Vector) bins[i][j_1];

          //      std::cerr << i << " smooth " << j_1 << " "  << j << " "  << j1 << "  "
          //           << ( v0.length() > v1.length() ?
          //                v0.length() / v1.length() :
          //                v1.length() / v0.length() ) << std::endl;

          if( Dot( v0, v1 ) > 0 &&
              ( v0.length() > v1.length() ?
                v0.length() / v1.length() :
                v1.length() / v0.length() ) < 10.0 ) {

            Vector center = (Vector) circle( bins[i][j_1],
                                             bins[i][j  ],
                                             bins[i][j1 ] );

            double rad = ((Vector) bins[i][j] - center).length();


            for( unsigned int s=0; s<add; s++ ) {
              Vector midPt = (Vector) bins[i][j_1] +
                (double) (add-s) / (double) (add+1) *
                ((Vector) bins[i][j] - (Vector) bins[i][j_1] );
                

              Vector midVec = midPt - center;

              midVec.normalize();

              newPts[add*j+s].first += center + midVec * rad;
              newPts[add*j+s].second += 1;

              midPt = (Vector) bins[i][j] +
                (double) (add-s) / (double) (add+1) *
                ((Vector) bins[i][j1] - (Vector) bins[i][j] );

              midVec = midPt - center;

              midVec.normalize();

              newPts[add*j1+s].first += center + midVec * rad;
              newPts[add*j1+s].second += 1;
            }
          }
        }

        for( int j=nodes-1; j>=0; j-- ) {

          for( unsigned int s=0; s<add; s++ ) {

            unsigned int k = add * j + s;

            if( newPts[k].second > 0 ) {
              
              newPts[k].first /= newPts[k].second;
              
              if( verboseFlag )
                std::cerr << i << "  insert  " << j << "  "
                          << newPts[k].first << std::endl;
              
              bins[i].insert( bins[i].begin()+j, newPts[k].first );
            }
          }
        }
      }

      // Remove the last point so it is possilble to see the groups.
      bins[i].erase( bins[i].end() );
    }
  }

  return toroidalWinding*(add+1)*nnodes;
}


unsigned int
FieldlineLib::
mergeOverlap( std::vector< std::vector < Point > > &bins,
              unsigned int &nnodes,
              unsigned int toroidalWinding,
              unsigned int windingGroupOffset,
              unsigned int island )
{
  Vector globalCentroid = Vector(0,0,0);;

  for( unsigned int i=0; i<toroidalWinding; i++ )
    for( unsigned int j=0; j<nnodes; j++ )
      globalCentroid += (Vector) bins[i][j];
  
  globalCentroid /= (toroidalWinding*nnodes);
    
  if( island ) {

    std::vector< std::vector< Point > > tmp_bins;

    tmp_bins.resize( toroidalWinding );

    for( unsigned int i=0; i<toroidalWinding; i++ ) {
      
      unsigned int startIndex;
      unsigned int middleIndex;
      unsigned int stopIndex;
      unsigned int nodes;
        
     std::vector< Point > points;
      points.resize( bins[i].size() );

      for( unsigned int j=0; j<bins[i].size(); j++ )
        points[j] = bins[i][j];

      unsigned int turns =
        islandProperties( points, globalCentroid,
                          startIndex, middleIndex, stopIndex, nodes );

      // Merge only if there are overlapping points.
      if( turns == 3) {

        if( nnodes == bins[i].size() )
          continue;

        // Store the overlapping points.
        for( unsigned int j=nnodes; j<bins[i].size(); j++ )
          tmp_bins[i].push_back( bins[i][j] );

        if( verboseFlag )
          std::cerr << i << " stored extra points "
                    << tmp_bins[i].size() << std::endl;

        // Erase all of the overlapping points.
        bins[i].erase( bins[i].begin()+nnodes, bins[i].end() );
          
        // Insert the first point so the curve is contiguous.
        bins[i].insert( bins[i].begin()+nnodes, bins[i][0] );

        unsigned int index_prediction = 1;
        unsigned int prediction_true  = 0;
        unsigned int prediction_false = 0;

        unsigned int modulo = (unsigned int)bins[i].size() - 1;
 
        // Insert the remaining points.
        for( unsigned int j=0; j<tmp_bins[i].size(); j++ ) {

          Vector v0 = (Vector) bins[i][0] -
            (Vector) tmp_bins[i][j];

          double angle = 0;
          double length = 99999;
          unsigned int angleIndex = 0;

          for( unsigned int k=1; k<bins[i].size(); k++ ) {

            Vector v1 = (Vector) bins[i][k] -
              (Vector) tmp_bins[i][j];

            double ang = acos( Dot(v0, v1) / (v0.length() * v1.length()) );

            if( angle < ang ) {
              angle = ang;
              angleIndex = k;
            }

            if( length < v1.length() ) {
              length = v1.length();
            }

            // Go on.
            v0 = v1;
          }

          // Insert it between the other two.
          if( angle > M_PI / 3.0 )
            bins[i].insert( bins[i].begin()+angleIndex, tmp_bins[i][j] );

          if( verboseFlag )
            std::cerr << i << "  " << modulo << "  " << j + nnodes
                      << "  Prediction " << index_prediction
                      << " actual " << angleIndex << "  "
                      << (angleIndex == index_prediction) << std::endl;

          // Check to see if the prediction and the actual index are
          // the same.
          if( angleIndex == index_prediction )
            prediction_true++;
          else // if( angleIndex != index_prediction )
            prediction_false++;

          // Predict where the next insertion will take place.
          if( (j+1) % modulo == 0 )
            index_prediction = 1 + (j+1) / modulo;
          else
            index_prediction = angleIndex + (j+1) / modulo + 2;
        }

        if( verboseFlag )
          std::cerr << "ToroidalWinding " << i << " inserted "
                    << prediction_true+prediction_false << " nodes "
                    << " True " << prediction_true
                    << " False " << prediction_false << std::endl;

        // If more of the predictions are incorrect than correct
        // insert based on the predictions.
        if( 0 && prediction_true < prediction_false ) {

          if( verboseFlag )
            std::cerr << "ToroidalWinding " << i << " bad predicted insertion ";

          unsigned int cc = 0;

          for( unsigned int j=0; j<tmp_bins[i].size(); j++ ) {

           std::vector< Point >::iterator inList =
              find( bins[i].begin(), bins[i].end(), tmp_bins[i][j] );
              
            if( inList != bins[i].end() ) {
              bins[i].erase( inList );

              cc++;
            }
          }

          if( verboseFlag )
            std::cerr << "removed " << cc << " points" << std::endl;

          unsigned int index = 1;
            
          for( unsigned int j=0; j<tmp_bins[i].size(); j++ ) {
            
            // Insert it between the other two.
            bins[i].insert( bins[i].begin()+index, tmp_bins[i][j] );

            if( verboseFlag )
              std::cerr << i << "  " << modulo << "  " << j + nnodes
                        << " actual " << index << std::endl;

            // Predict where the next insertion will take place.
            if( (j+1) % modulo == 0 )
              index = 1 + (j+1) / modulo;
            else
              index += (j+1) / modulo + 2;
          }
        }

        unsigned int start0  = 0;
        unsigned int end0    = 0;
        unsigned int start1  = 0;
        unsigned int end1    = 0;

        if( prediction_true > prediction_false ) {
          // See if any of the segments cross.
          for( unsigned int j=0; 0 && j<bins[i].size()-1; j++ ) {
              
            Point l0_p0 = bins[i][j];
            Point l0_p1 = bins[i][j+1];
              
            for( unsigned int k=j+2; k<bins[i].size()-1; k++ ) {
                
              Point l1_p0 = bins[i][k];
              Point l1_p1 = bins[i][k+1];
                
              if( intersect( l0_p0, l0_p1, l1_p0, l1_p1 ) == 1 ) {
                if( start0 == 0 ) {
                  start0  = j + 1;
                  end1    = k + 1;
                } else {
                  end0   = j + 1;
                  start1 = k + 1;

                  if( verboseFlag )
                    std::cerr << " merge self intersection " 
                              << start0 << "  " << end0 << "  "
                              << start1 << "  " << end1 << std::endl;

                  if( 0 ) {
                    std::vector < Point > tmp_bins[2];

                    for( unsigned int l=start0; l<end0; l++ )
                      tmp_bins[0].push_back( bins[i][l] );

                    for( unsigned int l=start1; l<end1; l++ )
                      tmp_bins[1].push_back( bins[i][l] );

                    bins[i].erase( bins[i].begin()+start1,
                                   bins[i].begin()+end1 );

                    bins[i].erase( bins[i].begin()+start0,
                                   bins[i].begin()+end0 );

                    for( unsigned int l=0; l<tmp_bins[1].size(); l++ )
                      bins[i].insert( bins[i].begin()+start0,
                                      tmp_bins[1][l] );

                    for( unsigned int l=0; l<tmp_bins[0].size(); l++ )
                      bins[i].insert( bins[i].begin() + start1 -
                                      tmp_bins[0].size() +
                                      tmp_bins[1].size(),
                                      tmp_bins[0][l] );
                  }

                  start0 = 0;
                }
              }
            }
          }
        }
      }
    }
  } else {

    std::vector< std::vector < Point > > tmp_bins(toroidalWinding);

   // This gives the minimal number of nodes for each group.
    nnodes = surfaceOverlapCheck( bins, windingGroupOffset );

    if( nnodes == 0 )
    {
      for( unsigned int i=1; i<toroidalWinding; i++ ) {
        if( nnodes > bins[i].size())
          nnodes = (unsigned int)bins[i].size();
      }
    }

    // If the offset and point ordering are opposite in directions
    // then the next group is the -windingGroupOffset. Otherwise if
    // they have the same direction then windingGroupOffset is the
    // next group.
    Vector intra = (Vector) bins[                 0][1] - (Vector) bins[0][0];
    Vector inter = (Vector) bins[windingGroupOffset][0] - (Vector) bins[0][0];

    int offset = (Dot( intra, inter ) < 0.0) ?
      toroidalWinding-windingGroupOffset : windingGroupOffset;

    for( unsigned int i=0; i<toroidalWinding; i++ )
    {
      // The next group
      unsigned int j = (i + offset) % toroidalWinding;

      // Add back in any nodes that may not overlap.
      unsigned int nodes = surfaceGroupCheck( bins, i, j, nnodes );

      // No more than one point should added.
      if( nodes > nnodes+1 )
        if( verboseFlag )
          std::cerr << "Surface fill " << i
                    << " nnodes mismatch " << nodes << std::endl;

      // Store the overlapping points.
      for( unsigned int j=nodes; j<bins[i].size(); j++ )
        tmp_bins[i].push_back( bins[i][j] );

      if( verboseFlag )
        std::cerr << i << " stored extra points "
                  << tmp_bins[i].size() << std::endl;

      // Erase all of the overlapping points.
      bins[i].erase( bins[i].begin()+nodes, bins[i].end() );

      // Insert the first point from the next toroidalWinding so the curve
      // is contiguous.
      bins[i].push_back( bins[(i+offset)%toroidalWinding][0] );
    }


    for( unsigned int i=0; i<toroidalWinding; i++ )
    {
      unsigned int toroidalWinding_prediction = (i+offset)%toroidalWinding;
      unsigned int index_prediction = 1;
      unsigned int prediction_true  = 0;
      unsigned int prediction_false = 0;

      for( unsigned int i0=0; i0<tmp_bins[i].size(); i0++ ) {

        double angle = 0;
        unsigned int index_wd = 0;
        unsigned int index_pt = 0;

        for( unsigned int j=0; j<toroidalWinding; j++ ) {

          Vector v0 = (Vector) bins[j][0] -
            (Vector) tmp_bins[i][i0];

          for( unsigned int j0=1; j0<bins[j].size(); j0++ ) {
            Vector v1 = (Vector) bins[j][j0] -
              (Vector) tmp_bins[i][i0];
        
            double ang = acos( Dot(v0, v1) / (v0.length() * v1.length()) );

            if( angle < ang ) {

              angle = ang;
              index_wd = j;
              index_pt = j0;
            }

            // Go on.
            v0 = v1;
          }
        }

        // Insert it between the other two.
        bins[index_wd].insert( bins[index_wd].begin()+index_pt,
                               tmp_bins[i][i0] );

        if( verboseFlag )
          std::cerr << "ToroidalWinding prediction "
                    << toroidalWinding_prediction
                    << " actual " << index_wd
                    << "  Index prediction  " << index_prediction
                    << " actual " << index_pt << "  "
                    << (index_wd == toroidalWinding_prediction &&
                        index_pt == index_prediction) << std::endl;

        // Check to see if the prediction of where the point was inserted
        // is correct;
        if( index_wd == toroidalWinding_prediction &&
            index_pt == index_prediction )
          prediction_true++;
        else 
          prediction_false++;

        // Prediction of where the next insertion will take place.
        index_prediction = index_pt + 2;
        toroidalWinding_prediction = index_wd;
      }

      if( verboseFlag )
        std::cerr << "ToroidalWinding " << i << " inserted "
             << prediction_true+prediction_false << " nodes "
             << " True " << prediction_true
             << " False " << prediction_false << std::endl;
    }

    // Remove the last point so it is possilble to see the groups.
    for( unsigned int i=0; i<toroidalWinding; i++ )
      bins[i].erase( bins[i].end() );
  }

  // Update the approximate node count.
  nnodes = 9999;

  for( unsigned int i=0; i<toroidalWinding; i++ )
    if( nnodes > bins[i].size() )
      nnodes = (unsigned int)bins[i].size();

  return nnodes;
}






void PCA( Skeleton::PointVector verts,
          std::vector< Skeleton::Point > &ev )
{
  // Estimate the center of the point set
  Skeleton::Point pos_sum(0,0);

  unsigned int num_points = 0;

  for (Skeleton::PointVector::iterator SL=verts.begin(); SL!= verts.end(); ++SL)
  {
    pos_sum.x += (*SL).x;
    pos_sum.y += (*SL).y;
    ++num_points;
  }

  pos_sum.x /= num_points;
  pos_sum.y /= num_points;

  // Estimate the pca of the set of vertices
  float mat[2][2] = {{0.0, 0.0}, {0.0, 0.0}};

  for( Skeleton::PointVector::iterator SL=verts.begin(); SL != verts.end(); ++SL )
  {
    float off[2]; 
    
    off[0] = (*SL).x - pos_sum.x;
    off[1] = (*SL).y - pos_sum.y;
    
    mat[0][0] += off[0] * off[0];
    mat[0][1] += off[0] * off[1];
    mat[1][1] += off[1] * off[1];
  }

  mat[1][0] = mat[0][1];

  for( unsigned int i=0; i<2; ++i )
    for( unsigned int j=0; j<2; ++j )
      mat[i][j] /= num_points;


  double la = mat[0][0];
  double lb = mat[0][1];
  double lc = mat[1][0];
  double ld = mat[1][1];

  double A = 1.0;
  double B = -(la + ld);
  double C = (la * ld - lb * lc);

  double delta = B*B - 4.0*A*C;

  double evalues[2] = {0.0, 0.0} ; ///TODO: check on fix for uninitialized variables
  Skeleton::Point evec[2];

  if( delta >= 0.0 )
  {
    evalues[0] = (-B - sqrt(delta)) / 2.0;
    evalues[1] = (-B + sqrt(delta)) / 2.0;
    
    // For real eigen values, calculate the eigen vectors
    evec[0].x = mat[1][1] - evalues[0];
    evec[0].y = -mat[1][0];

    evec[1].x = -(mat[1][1] - evalues[1]);
    evec[1].y = mat[1][0];

    evec[0].normalize();
    evec[1].normalize();
  }
  else
  {
    evec[0].x = evec[0].y = 0.;
    evec[1].x = evec[1].y = 0.;
  }

  if( evalues[0] < evalues[1] )
  {
    Skeleton::Point axis;

    axis.x = evec[0].x;
    axis.y = evec[0].y;
    ev.push_back(axis);

    axis.x = evec[1].x;
    axis.y = evec[1].y;
    ev.push_back(axis);
  }
  else // if( evalues[0] >= evalues[1] )
  {
    Skeleton::Point axis;

    axis.x = evec[1].x;
    axis.y = evec[1].y;
    ev.push_back(axis);

    axis.x = evec[0].x;
    axis.y = evec[0].y;
    ev.push_back(axis);
  }
}


// ****************************************************************************
//  Method: FieldlineLib::findIslandCenters
//
//  Purpose: Finds the geometric centers of an island chain.
//
//  Arguments:
//
//  Returns:      void
//
//  Programmer: Allen Sanderson
//  Creation:   Wed Feb 25 09:52:11 EST 2009
//
//  Modifications:
//
// ****************************************************************************

void
FieldlineLib::findIslandCenters( std::vector< Point > &puncturePts,
                                 unsigned int islands,
                                 unsigned int offset,
                                 unsigned int nnodes,
                                 unsigned int moduloValue,
                                 std::vector< Point > &centers,
                                 std::vector< Vector > &axis )
{
#ifdef STRAIGHTLINE_SKELETON

  // Loop through each toroidial group
  for( unsigned int i=0; i<islands; ++i ) 
  {
    if( verboseFlag )
      std::cerr << "Island " << i << "  ";

    // temp storage incase the order needs to be reversed.
   std::vector< Point > tmp_points;

    bool selfIntersect = false;

    // Loop through each point in toroidial group
    for( unsigned int j=i, jc=0;
         jc<nnodes;
         j+=offset, ++jc )
    {
//       if( verboseFlag )
//         std::cerr << "hull points  " << jc << "  " << j%moduloValue << "  " 
//                   <<  puncturePts[j%moduloValue]
//                   << std::endl;
     
      tmp_points.push_back( puncturePts[j%moduloValue] );
    }

    // Points for the convex hull check
    for( unsigned int j=0; j<tmp_points.size()-1; ++j )
    {
      // Check for self intersections
      for( unsigned int k=j+2; k<tmp_points.size()-1; ++k )
      {
        int result;

        if((result = intersect( tmp_points[j], tmp_points[j+1],
                                tmp_points[k], tmp_points[k+1] )))
        {
          if( verboseFlag )
            std::cerr << " hull self intersects  "
                      << j << "  " << k << std::endl;
          
          // Self intersection found so skip this island.
          selfIntersect = true;
          break;
        }
      }

      if( selfIntersect )
        break;
    }
 
    // Self intersection found so skip this island as the skeleton
    // will not work.
    if( selfIntersect )
      continue;

    // Get the convex hull and the direction.
    int direction = 0;
    hullCheck( tmp_points, direction );
      
    // Store the points a 2D vector.
    Skeleton::PointVector pointVec;

    for( unsigned int j=0; j<tmp_points.size()-1; ++j )
      pointVec.push_back( Skeleton::Point( tmp_points[j].x,
                                           tmp_points[j].z ) );
      
    // If the points are clockwise reverse them as the skeleton needs
    // the points to be in a counter clockwise direction.
    if( direction == 1 )
      reverse( pointVec.begin(), pointVec.end() );

    if( verboseFlag )
      std::cerr << " Skeleton check ";
      
    Skeleton::Skeleton s (Skeleton::makeSkeleton (pointVec));


    if( 1 )
    {
      // Convert the skeleton to a graph
      DIR_Graph *skl_graph = new DIR_Graph(0, 0);
      skl_graph->convert_SKL_to_graph (s);
      
      // Compute the approximate geometric center 
      Skeleton::Point center = skl_graph->get_centroid();

      centers.push_back( Point( center.x, 0, center.y ) );
    }
    else
    {
      // Compute the approximate geometric center 
      Point center = findSkeletonCenter( s, (unsigned int)pointVec.size() );

      centers.push_back( center );
    }



//    Skeleton::PointVector pointVec;   // load the boundary point of the hull

    pointVec.clear();

    for( unsigned int j=0; j<tmp_points.size()-1; ++j )
      pointVec.push_back( Skeleton::Point( tmp_points[j].x,
                                           tmp_points[j].z ) );
 
    std::vector< Skeleton::Point > axes;

    PCA( pointVec, axes );  // the first axes is the shortest axis that we want to march along
    
    Point minBoundaryPoint, maxBoundaryPoint;
    double minDistance = 1.0e9, maxDistance = -1.0e9;

    for( unsigned int j=0; j<tmp_points.size()-1; ++j )
    {
      double tmp = avtVector( centers[i] - tmp_points[j] ).length();
      
      if( minDistance > tmp )
      {
        minDistance = tmp;
        minBoundaryPoint = tmp_points[j];
      }
      else if( maxDistance < tmp )
      {
        maxDistance = tmp;
        maxBoundaryPoint = tmp_points[j];
      }
    }

    axis.push_back( Vector( axes[0].x, 0, axes[0].y ) *
                    (minBoundaryPoint-centers[i]).length() );

    axis.push_back( Vector( axes[1].x, 0, axes[1].y ) *
                    (maxBoundaryPoint-centers[i]).length() );
  }
#endif
}


// ****************************************************************************
//  Method: FieldlineLib::findSkeletonCenter
//
//  Purpose: Finds the geometric center of skeleton
//
//  Arguments:
//
//  Returns:    Point
//
//  Programmer: Allen Sanderson
//  Creation:   Wed Feb 25 09:52:11 EST 2009
//
//  Modifications:
//
// ****************************************************************************
Point
FieldlineLib::findSkeletonCenter( Skeleton::Skeleton &s,
                                  unsigned int nHullPts )
{
  Point center;

    double cordLength = 0;
    int nEndpoints = 0;

    std::list<Skeleton::SkeletonLine>::iterator SL;

    std::map< int, int > indexCount;
    std::map< int, int >::iterator ic;

    int startIndex = -1, endIndex = -1; ///TODO: check for fix of uninitialized variable

    // It is possible to get more than 2 end points. As such, loop
    // through leaving at least half the segmetns while naively
    // pruning the shortest end point segments in hopes that a valid
    // skeleton remains.
    unsigned int cc = 0;

    while( cc < s.size() / 2 )
    {
      SL = s.begin();
 
      cordLength = 0;
      nEndpoints = 0;

      indexCount.clear();
     
      // Remove all of the points on the boundary while counting the
      // number of time each point is used.
      while ( SL != s.end() )
      {
        // Delete all of the hull points only happens on the first
        // iteration.
        if( (unsigned int)(*SL).lower.vertex->ID < nHullPts ||
            (unsigned int)(*SL).higher.vertex->ID < nHullPts ) 
        {
          // Boundary point so delete it.
          std::list<Skeleton::SkeletonLine>::iterator deleteSL = SL;
          ++SL;
          s.erase( deleteSL );
        }
        else
        {
          // Count the number of times each end point is
          // used. Hopefully each is used twice with the exception of
          // the skeleton end points.
        
          // Lower end point
          int index = (*SL).lower.vertex->ID;
        
          ic = indexCount.find( index );
        
          if( ic == indexCount.end() )
            indexCount.insert( std::pair<int, int>( index,1) );
          else (*ic).second++;
        
          // Higher end point
          index = (*SL).higher.vertex->ID;
        
          ic = indexCount.find( index );
        
          if( ic == indexCount.end() )
            indexCount.insert( std::pair<int, int>( index,1) );
          else (*ic).second++;
        
          ++SL;
        }
      }                   
        
      // Loop through a second time and find the end segment that has
      // the smallest length while getting the cord length of the
      // remaining interior segments
      std::list<Skeleton::SkeletonLine>::iterator minSL;
            
      double minCordLength = 1.0e12;
    
      SL = s.begin();

      while ( SL != s.end() )
      {
        double tmp = sqrt( ((*SL).higher.vertex->point.x-
                            (*SL).lower.vertex->point.x) *
                           ((*SL).higher.vertex->point.x-
                            (*SL).lower.vertex->point.x) +
                           
                           ((*SL).higher.vertex->point.y -
                            (*SL).lower.vertex->point.y) * 
                           ((*SL).higher.vertex->point.y -
                            (*SL).lower.vertex->point.y) );

        ic = indexCount.find( (*SL).lower.vertex->ID );
        
        if( verboseFlag )
          std::cerr << (*SL).lower.vertex->ID << " ("
                    << (*ic).second << ")  ";
        
        ic = indexCount.find( (*SL).higher.vertex->ID );

        if( verboseFlag )
          std::cerr << (*SL).higher.vertex->ID << " ("
                    << (*ic).second << ")   "
                    <<(*SL).lower.vertex->point << "  "
                    << (*SL).higher.vertex->point << "  "
                    << tmp
                    << std::endl;

        // Running cord length.
        cordLength += tmp;
        
        // Lower end point
        int index = (*SL).lower.vertex->ID;
        
        ic = indexCount.find( index );
        
        if( (*ic).second == 1 )
        {
          if( minCordLength > tmp )
          {
            minCordLength = tmp;
            minSL = SL;
          }
        }

        // Higher end point
        index = (*SL).higher.vertex->ID;
        
        ic = indexCount.find( index );
        
        if( (*ic).second == 1 )
        {
          if( minCordLength > tmp )
          {
            minCordLength = tmp;
            minSL = SL;
          }
        }
        
        ++SL;
      }


      // Loop through the counts and determine if there are two end
      // points which describe the skeleton.
      ic = indexCount.begin();
      
      while( ic != indexCount.end() )
      {
        // Find the end points that are used but once. 
        // There should only be two but if not we can try trimming.
        if( (*ic).second == 1 )
        {
          if( nEndpoints == 1 )
            startIndex = (*ic).first;
          else if( nEndpoints == 0 )
            endIndex = (*ic).first;
          
          ++nEndpoints;
        }

        ++ic;
      }

      if( nEndpoints < 2 )
      {
        if( verboseFlag )
          std::cerr << "Not enough end points " << nEndpoints << std::endl;
        continue;
      }
      else if( nEndpoints > 2 )
      {
        if( verboseFlag )
          std::cerr << "Too many end points " << nEndpoints << "  "
                    << "removing the smallest segment "
                    << (*minSL).lower.vertex->ID << "  "
                    << (*minSL).higher.vertex->ID << "  "
                    << std::endl;
        
        s.erase( minSL );
      }
      else //if( nEndpoints == 2 )
      {
        break;
      }

      ++cc;
    }

    // Have two end points now check for and remove double segments.
    if( nEndpoints == 2 )
    {
      ic = indexCount.begin();
    
      while( ic != indexCount.end() )
      {
        if( (*ic).second > 2 )
        {
           if( verboseFlag )
             std::cerr << "double segment ??? " << (*ic).first << std::endl;
        
          SL = s.begin();
        
          while( SL != s.end() &&
               (*SL).lower.vertex->ID != (*ic).first &&
               (*SL).higher.vertex->ID != (*ic).first )
          SL++;
        
          // Have an index so process
          if( SL != s.end() )
          {
            // Remove double segments;
            if( (*SL).lower.vertex->ID == (*ic).first &&
                (*indexCount.find( (*SL).higher.vertex->ID )).second > 2 )
            {
              if( verboseFlag )
                std::cerr << "removing double segment ??? "
                          << (*ic).first << "  "
                          << (*indexCount.find( (*SL).higher.vertex->ID )).first
                          << std::endl;

              (*ic).second--;
              (*indexCount.find( (*SL).higher.vertex->ID )).second--;
              
              s.erase( SL );
            }
          
            // Remove double segments;
            if( (*SL).higher.vertex->ID == (*ic).first &&
                (*indexCount.find( (*SL).lower.vertex->ID )).second > 2 )
            {
              if( verboseFlag )
                std::cerr << "removing double segment ??? "
                          << (*ic).first << "  "
                          << (*indexCount.find( (*SL).lower.vertex->ID )).first
                          << std::endl;
              
              (*ic).second--;
              (*indexCount.find( (*SL).lower.vertex->ID )).second--;
              
              s.erase( SL );
            }
          }                    
        }
      
        ++ic;
      }
    
      if( verboseFlag )
        std::cerr << " New skeleton "
                  << "start index " << startIndex
                  << "  end index " << endIndex
                  << std::endl;
      
      double currentCord = 0;

      int nextIndex = startIndex;


      // Loop through all of the segments which are described
      // as set of paired points.
      while ( s.begin() != s.end() )
      {
        // Index of the first point.
        SL = s.begin();

        while( SL != s.end() &&
               (*SL).lower.vertex->ID != nextIndex &&
               (*SL).higher.vertex->ID != nextIndex )
          SL++;
                    
        // Have an index so process
        if( SL != s.end() )
        {
//      std::cerr << "Island " << j
//                << " Line segment " << nlines++
//                << " index " << nextIndex;

          int lastIndex = nextIndex;

          // Index of the leading point.
          if( (*SL).lower.vertex->ID == nextIndex )
            nextIndex = (*SL).higher.vertex->ID;
          else // if( (*SL).higher.vertex->ID == startIndex )
            nextIndex = (*SL).lower.vertex->ID;
            
//      std::cerr << " index " << nextIndex
//                << std::endl;

          double localCord = sqrt( ((*SL).higher.vertex->point.x-
                                    (*SL).lower.vertex->point.x) *
                                   ((*SL).higher.vertex->point.x-
                                    (*SL).lower.vertex->point.x) +
                                     
                                   ((*SL).higher.vertex->point.y -
                                    (*SL).lower.vertex->point.y) * 
                                   ((*SL).higher.vertex->point.y -
                                    (*SL).lower.vertex->point.y) );

          // Check to see if the segment spans the mid cord.
          if( currentCord < cordLength/2 &&
              cordLength/2.0 <= currentCord + localCord )
          {
            double t = (cordLength/2.0 - currentCord) / localCord;

            // Already updated the nextIndex so use the lastIndex
            // for figuring out the closest to the cordlength center.
            if( (*SL).lower.vertex->ID == lastIndex )
            {
              center.x = (*SL).lower.vertex->point.x + t *
                ((*SL).higher.vertex->point.x-(*SL).lower.vertex->point.x);
              center.y = 0;
              center.z = (*SL).lower.vertex->point.y + t *
                ((*SL).higher.vertex->point.y-(*SL).lower.vertex->point.y);
            }

            else //if( (*SL).higher.vertex->ID == lastIndex )
            {
              center.x = (*SL).higher.vertex->point.x + t *
                ((*SL).lower.vertex->point.x-(*SL).higher.vertex->point.x);
              center.y = 0;
              center.z = (*SL).higher.vertex->point.y + t *
                ((*SL).lower.vertex->point.y-(*SL).higher.vertex->point.y);
            }
          }
        
          // Update the running cord length so that the mid cord is
          // not calculated again.
          currentCord += localCord;
        
          // Remove this point from the list so that it is
          // not touched again.
          std::list<Skeleton::SkeletonLine>::iterator deleteSL = SL;
          ++SL;
          s.erase( deleteSL );
        }
        else if( nextIndex != endIndex )
        {
//      if( verboseFlag )
//        std::cerr << "Did not find end index  "
//                  << nextIndex << "  " <<  endIndex
//                  << std::endl;

          break;
        }
      }
    }
    else
    {
        if( verboseFlag )
          std::cerr << "Dispite trimming too many end points." << std::endl;
    }

    if( verboseFlag )
      std::cerr << "O Point " << center << std::endl;

    return center;
}


//===================================================================
// Adapted from Tolga Birdal

// Compute the q values in the equation
double Otsu::Px( unsigned int init, unsigned int end,
                std::vector< unsigned int > &histo)
{
  int sum = 0;

  for (unsigned int i=init; i<end; ++i)
    sum += histo[i];

  return (double) sum;
}

// Compute the mean values in the equation (mu)
double Otsu::Mx( unsigned int init, unsigned int end,
                std::vector< unsigned int > &histo)
{
  int sum = 0;

  for (unsigned int i=init; i<end; ++i)
    sum += i * histo[i];
      
  return (double) sum;
}

// Find the maximum element in a vector
unsigned int Otsu::findMaxVet( std::vector< double > &vet, double &maxVet)
{
  maxVet = 0;
  unsigned int index = 0;
      
  for (unsigned int i=1; i<vet.size()-1; ++i)
  {
    if( maxVet < vet[i] )
    {
      maxVet = vet[i];
      index = i;
    }
  }

  return index;
}
    
// Compute the histogram
void Otsu::getHistogram( std::vector< std::pair< unsigned int, double > >& stats,
                         std::vector< unsigned int > &histo )
{
  unsigned int nbins = (unsigned int)stats.size() / 2;

  if( nbins < 10 )
    nbins = (unsigned int)stats.size();

  double min = 1.0e9, max = -1.0e9;
    
  for( unsigned int i=0; i<stats.size(); ++i )
  {
    if( min > stats[i].second )
      min = stats[i].second;
    
    if( max < stats[i].second )
      max = stats[i].second;
  }

  histo.resize( nbins );
    
  for( unsigned int i=0; i<nbins; ++i )
    histo[i] = 0;
    
  for( unsigned int i=0; i<stats.size(); ++i )
  {
    unsigned int index =
      ( ((double) nbins-1.0) * (stats[i].second-min) / (max-min) );
      
    histo[index]++;
  }
  
//   WindingPair zeropair[2];

//   zeropair[0] = WindingPair(0,0);
//   zeropair[1] = WindingPair(0,0);

//   unsigned int zerostart=nbins, zerostop=nbins;

//   for( unsigned int i=0; i<nbins; ++i )
//   {
//     if( histo[i] == 0 && zerostart == nbins )
//     {
//       zerostart = i;
//       zerostop  = i;
//     }
//     else if( histo[i] != 0 && zerostart != nbins )
//     {
//       zerostop = i-1;
      
//       if( zerostop-zerostart > zeropair[0].second-zeropair[0].first )
//       {
//         zeropair[1] = zeropair[0];
//         zeropair[0] = WindingPair(zerostart, zerostop);
//       }
//       else if( zerostop-zerostart > zeropair[1].second-zeropair[1].first )
//       {
//         zeropair[1] = WindingPair(zerostart, zerostop);
//       }
      
//       zerostart = nbins;
//     }

//     int diff;
    
//     if( i )
//       diff = abs((int)histo[i] - (int)histo[i-1]);
//     else
//       diff = 0;
    
//     std::cerr << i << "  " << histo[i] << "  " << diff << "  " << histo[i] << std::endl;
//   }

//   std::cerr << std::endl
//        << "zero pairs  "
//        << zeropair[0].first << "-" << zeropair[0].second << "  "
//        << zeropair[0].second-zeropair[0].first+1 << "     "
//        << zeropair[1].first << "-" << zeropair[1].second << "  "
//        << zeropair[1].second-zeropair[1].first+1 << std::endl;
}

// find otsu threshold
void Otsu::
getOtsuThreshold2( std::vector< std::pair< unsigned int, double > >& stats,
                   double &threshold, double &maxVet )
{
  if( stats.size() == 1 )
  {
    threshold = stats[0].second;
    maxVet = 0;
    
    return;
  }
  
  else if( stats.size() == 2 )
  {
    threshold = (stats[0].second+stats[1].second)/2.0;
    maxVet = 0;
    
    return;
  }

 std::vector< unsigned int > histo;
  getHistogram( stats, histo );
  
  unsigned int nbins = (unsigned int)histo.size();
  
  unsigned int index = 0; ///TODO: check on uninitialized var warning
  maxVet = 0;
  
  // Loop through all possible t values and maximize between class variance.
  for( unsigned int i=1; i<nbins-1; ++i )
  {
    double p0 = Px(0,     i, histo);
    double p1 = Px(i, nbins, histo);
    
    double p01 = p0 * p1;
    
    if (p01 == 0) p01 = 1.0;
    
    double m0 = Mx(0,     i, histo);
    double m1 = Mx(i, nbins, histo);

    double diff = m0 * p1 - m1 * p0;
      
    double vet = diff * diff / p01;

    if( maxVet < vet )
    {
      maxVet = vet;
      index = i;
    }
  }

  // Calculate the min and max from the values to get the threshold.
  double min = 1.0e9, max = -1.0e9;
  
  for( unsigned int i=0; i<stats.size(); ++i )
  {
    if( min > stats[i].second )
      min = stats[i].second;
    
    if( max < stats[i].second )
      max = stats[i].second;
  }

  threshold = min + (double) index / ((double) nbins - 1.0) * (max-min);
}


// Find the Otsu threshold for three division
void Otsu::
getOtsuThreshold3( std::vector< std::pair< unsigned int, double > >& stats,
                   double &threshold0, double &threshold1,
                   double &maxVet )
{
 std::vector< unsigned int > histo;
  getHistogram( stats, histo );

  unsigned int nbins = (unsigned int)histo.size();

  unsigned int index0 = 0; ///TODO: check on fix for uninitialized warning
  unsigned int index1 = 0; ///TODO: check on fix for uninitialized warning

  maxVet = 0;

  // Loop through all possible t values and maximize between class variance
  for ( unsigned int i=1; i<nbins-2; ++i )
  {
    double p0 = Px(0, i, histo);
    double m0 = Mx(0, i, histo);

    for ( unsigned int j=i+1; j<nbins-1; ++j )
    {
      double p1 = Px(i,     j, histo);
      double p2 = Px(j, nbins, histo);

      double p01 = p0 * p1;
      double p12 = p1 * p2;

      if (p01 == 0) p01 = 1;
      if (p12 == 0) p12 = 1;

      double m1 = Mx(i,     j, histo);
      double m2 = Mx(j, nbins, histo);
        
      double diff0 = (m0 * p1 - m1 * p0);
      double vet0 = diff0 * diff0 / p01;

      double diff1 = (m1 * p2 - m2 * p1);
      double vet1 = diff1 * diff1 / p12;

      if( maxVet < vet0 + vet1 )
      {
        maxVet = vet0 + vet1;
        index0 = i;
        index1 = j;
      }
    }
  }

  // Calculate the min and max from the values to get the thresholds.
  double min = 1.0e9, max = -1.0e9;
    
  for( unsigned int i=0; i<stats.size(); ++i )
  {
    if( min > stats[i].second )
      min = stats[i].second;
      
    if( max < stats[i].second )
      max = stats[i].second;
  }

  threshold0 = min + (double) index0 / ((double) nbins-1.0) * (max-min);
  threshold1 = min + (double) index1 / ((double) nbins-1.0) * (max-min);
}



//===================================================================

//===================================================================
// Copyright 2001, softSurfer (www.softsurfer.com)
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.

//===================================================================
// Isleft(): tests if a point is Left|On|Right of an infinite line.
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and P1
//            =0 for P2 on the line
//            <0 for P2 right of the line
//    See: the January 2001 Algorithm on Area of Triangles
inline float
isLeft( Point P0, Point P1, Point P2 )
{
    return (P1.x - P0.x)*(P2.z - P0.z) - (P2.x - P0.x)*(P1.z - P0.z);
}

//===================================================================
// ptcmp(): 
//     Input:  two points to compare
//
//     Return: 1 if point 0 is "greater than" point 1
//            -1 if point 0 is "less than" point 1
//             0 if point 0 is "equal to" point 1
inline int ptcmp( const void* v0, const void* v1 )
{
  Point *p0 = (Point *) v0;
  Point *p1 = (Point *) v1;

  if( p0->x > p1->x ) return 1;
  else if( p0->x < p1->x ) return -1;
  else if( p0->y > p1->y ) return 1;
  else if( p0->y < p1->y ) return -1;
  else return 0;
}

//===================================================================
// chainHull_2D(): Andrew's monotone chain 2D convex hull algorithm
//     Input:  pts = an array of 2D points 
//
//     Output: hullPts = an array of the convex hull vertices
//                       in a clockwise order
//     Return: the number of points in hullPts
int chainHull_2D( std::vector< std::pair< Point, unsigned int > > &pts,
                  std::vector< std::pair< Point, unsigned int > > &hullPts,
                  int direction  )
{
    //  Presorted by increasing x- and y-coordinates
    qsort( &(pts[0]), pts.size(),
           sizeof( std::pair< Point, unsigned int > ), ptcmp );

    int n = (unsigned int)pts.size();

    hullPts.resize( n+1 );

    // the output array hullPts[] will be used as the stack
    int    bot=0, top=(-1);  // indices for bottom and top of the stack
    int    i;                // array scan index

    // Get the indices of points with min x-coord and min|max y-coord
    int minmin = 0, minmax;
    float xmin = pts[0].first.x;
    for (i=1; i<n; i++)
        if (pts[i].first.x != xmin) break;
    minmax = i-1;
    if (minmax == n-1) {       // degenerate case: all x-coords == xmin
        hullPts[++top] = pts[minmin];
        if (pts[minmax].first.z != pts[minmin].first.z) // a nontrivial segment
            hullPts[++top] = pts[minmax];
        hullPts[++top] = pts[minmin];           // add polygon endpoint
        return top+1;
    }

    // Get the indices of points with max x-coord and min|max y-coord
    int maxmin, maxmax = n-1;
    float xmax = pts[n-1].first.x;
    for (i=n-2; i>=0; i--)
        if (pts[i].first.x != xmax) break;
    maxmin = i+1;

    // Compute the lower hull on the stack H
    hullPts[++top] = pts[minmin];      // push minmin point onto stack
    i = minmax;
    while (++i <= maxmin)
    {
        // the lower line joins pts[minmin] with pts[maxmin]
        if (isLeft( pts[minmin].first,
                    pts[maxmin].first, pts[i].first) >= 0 && i < maxmin)
            continue;          // ignore pts[i] above or on the lower line

        while (top > 0)        // there are at least 2 points on the stack
        {
            // test if pts[i] is left of the line at the stack top
            if (isLeft( hullPts[top-1].first,
                        hullPts[top].first, pts[i].first) > 0)
                break;         // pts[i] is a new hull vertex
            else
                top--;         // pop top point off stack
        }
        hullPts[++top] = pts[i];       // push pts[i] onto stack
    }

    // Next, compute the upper hull on the stack H above the bottom hull
    if (maxmax != maxmin)      // if distinct xmax points
        hullPts[++top] = pts[maxmax];  // push maxmax point onto stack
    bot = top;                 // the bottom point of the upper hull stack
    i = maxmin;
    while (--i >= minmax)
    {
        // the upper line joins pts[maxmax] with pts[minmax]
        if (isLeft( pts[maxmax].first,
                    pts[minmax].first, pts[i].first) >= 0 && i > minmax)
            continue;          // ignore pts[i] below or on the upper line

        while (top > bot)    // at least 2 points on the upper stack
        {
            // test if pts[i] is left of the line at the stack top
            if (isLeft( hullPts[top-1].first,
                        hullPts[top].first, pts[i].first) > 0)
                break;         // pts[i] is a new hull vertex
            else
                top--;         // pop top point off stack
        }
        hullPts[++top] = pts[i];       // push pts[i] onto stack
    }

    if (minmax != minmin)
      hullPts[++top] = pts[minmin];  // push joining endpoint onto stack

    hullPts.resize( top+1 );

    return (unsigned int)hullPts.size();
}

//===================================================================
