/*
* Copyright (c) 2000 - 2017, Lawrence Livermore National Security, LLC
* Produced at the Lawrence Livermore National Laboratory
* LLNL-CODE-442911
* All rights reserved.
*
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

#include <vtkSphericalTransform.h>
#include <vtkGeoProjection.h>
#include <vtkGeoTransform.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkStructuredGrid.h>
#include <vtkRectilinearGrid.h>
#include <vtkVisItUtility.h>
#include <DebugStream.h>

// ************************************************************************* //
//  File: avtCartographicProjectionFilter.C
// ************************************************************************* //

#include <avtCartographicProjectionFilter.h>


// ****************************************************************************
//  Method: avtCartographicProjectionFilter constructor
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

avtCartographicProjectionFilter::avtCartographicProjectionFilter()
{
  //this->ProjType = CartographicProjectionAttributes::Hammer;
}


// ****************************************************************************
//  Method: avtCartographicProjectionFilter destructor
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
//  Modifications:
//
// ****************************************************************************

avtCartographicProjectionFilter::~avtCartographicProjectionFilter()
{
}


// ****************************************************************************
//  Method:  avtCartographicProjectionFilter::Create
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

avtFilter *
avtCartographicProjectionFilter::Create()
{
    return new avtCartographicProjectionFilter();
}


// ****************************************************************************
//  Method:      avtCartographicProjectionFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

void
avtCartographicProjectionFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const CartographicProjectionAttributes*)a;
}


// ****************************************************************************
//  Method: avtCartographicProjectionFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtCartographicProjectionFilter with the given
//      parameters would result in an equivalent avtCartographicProjectionFilter.
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
// ****************************************************************************

bool
avtCartographicProjectionFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(CartographicProjectionAttributes*)a);
}


// ****************************************************************************
//  Method: avtCartographicProjectionFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the CartographicProjection filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: jfavre -- generated by xml2avt
//  Creation:   Fri Aug 31 10:12:22 PDT 2012
//
//  Modifications:
//    Kathleen Biagas, Thu Oct 18 10:09:56 PDT 2012
//    Preserve coordinate data type.
//
//    Jean Favre, Mon Jan  7 17:02:19 CET 2013
//    Added a Central Meridian
//    Added support for PolyData (typically, they are shapefiles for continent lines)
//
//    Eric Brugger, Wed Jul 23 08:30:25 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtCartographicProjectionFilter::ExecuteData(avtDataRepresentation *in_dr)
{
/* The projection mappings in VTK, based on the libproj4 library, expect data with coordinates
*  in the range [-180, 180] for longitude, and [-90, 90] for latitude.
*  Further, mapped grids cannot sit on Rectilinear Grids. So, we start by converting rectilinear grids
*  to structured grids (i.e. fully curvilinear)
*  A new grid is created based on the input, coordinates are projected, then remapped to the classic range
*  of [-180, 180] for longitude, and [-90, 90] for latitude
*/
  //
  // Get the VTK data set.
  //
  vtkDataSet *in_ds = in_dr->GetDataVTK();

  int do_type = in_ds->GetDataObjectType();
  int dims[3], numPts = in_ds->GetNumberOfPoints();
  vtkPointSet *ds;
  vtkCellArray *ca;
  double in_bounds[6], out_bounds[6], tol = 10., p_pt[3], c_pt[3];

  vtkPoints *inPts = vtkVisItUtility::NewPoints(in_ds);
  inPts->Allocate(numPts);
  inPts->SetNumberOfPoints(numPts);

  for(int i=0; i < numPts; i++)
    {
    in_ds->GetPoint(i, c_pt);
// the next two lines used for our Geo-physics grid
    //c_pt[0] = (c_pt[2] - M_PI)*179.999/M_PI;
    //c_pt[1] = (c_pt[1] - M_PI_2)*89.999/M_PI_2;
// this is used for our Geo-physics grid
    c_pt[2] = 0.0;
    inPts->SetPoint(i, c_pt);
    }

  inPts->GetBounds(in_bounds);
  debug4 << "Input Bounds[6] = " << in_bounds[0] << " "  << in_bounds[1] << " " << in_bounds[2] << " " << in_bounds[3] << " " << in_bounds[4] << " " << in_bounds[5] <<  endl;

  vtkGeoProjection *proj = vtkGeoProjection::New();
  proj->SetCentralMeridian(atts.GetCentralMeridian());
  debug4 << "Central Meridian = " << proj->GetCentralMeridian() <<   endl;
  vtkGeoTransform *geoXform = vtkGeoTransform::New();

// Here we could have as many of the 100+ projections available.
// I'm not sure people want to see a pull-down list with 100+ items
  
  proj->SetName(atts.ProjectionID_ToString(atts.GetProjectionID()).c_str());
  debug4 << "setting name of projection " << atts.ProjectionID_ToString(atts.GetProjectionID()) <<   endl;
  geoXform->SetDestinationProjection(proj);

  vtkPoints *newPoints = vtkPoints::New(inPts->GetDataType());
  geoXform->TransformPoints(inPts, newPoints);

  newPoints->GetBounds(out_bounds);
  debug4 << "Output Bounds[6] = " << out_bounds[0] << " "  << out_bounds[1] << " " << out_bounds[2] << " " << out_bounds[3] << " " << out_bounds[4] << " " << out_bounds[5] <<  endl;

// rescale coordinates to the original set of longitude and latitude ranges.
// this is motivated by the fact that the continents shapefile provided by Brad for the SC'12 tutorial
// does not range fully to the North Pole. It was wrong to stretch its projection to +90 degrees.
  double alpha;
  for(int i=0; i < numPts; i++)
    {
    newPoints->GetPoint(i, c_pt);
    alpha = (c_pt[0] - out_bounds[0]) / (out_bounds[1]-out_bounds[0]);
    c_pt[0] = alpha * in_bounds[1] + (1.0 - alpha) * in_bounds[0];
    alpha = (c_pt[1] - out_bounds[2]) / (out_bounds[3]-out_bounds[2]);
    c_pt[1] = alpha * in_bounds[3] + (1.0 - alpha) * in_bounds[2];
    // only touch the X and Y coordinates. Ignore Z
    newPoints->SetPoint(i, c_pt);
    }
  // mark the coordinates as Modified to force a call to vtkPoints::ComputeBounds()
  // otherwise a ResetView() has trouble resetting.
  newPoints->Modified();
  newPoints->GetBounds(out_bounds);
  debug4 << "Output Bounds[6] = " << out_bounds[0] << " "  << out_bounds[1] << " " << out_bounds[2] << " " << out_bounds[3] << " " << out_bounds[4] << " " << out_bounds[5] <<  endl;

  vtkCellArray *ca_n;
  int changeOfSigns = 0, k;

  switch(do_type) {
    case VTK_STRUCTURED_GRID:
      ds = vtkStructuredGrid::New();
      static_cast<vtkStructuredGrid *>(in_ds)->GetDimensions(dims);
// make sure 3rd dimension is set in case we were given 2D datasets (e.g pressure.cdf)
      dims[2] = 1;
      static_cast<vtkStructuredGrid *>(ds)->SetDimensions(dims);
      ds->ShallowCopy(in_ds);
    break;
    case VTK_RECTILINEAR_GRID:
      ds = vtkStructuredGrid::New();
      static_cast<vtkRectilinearGrid *>(in_ds)->GetDimensions(dims);
// make sure 3rd dimension is set in case we were given 2D datasets (e.g pressure.cdf)
      dims[2] = 1;
      static_cast<vtkStructuredGrid *>(ds)->SetDimensions(dims);
      ds->ShallowCopy(in_ds);
    break;
    case VTK_POLY_DATA:
// some special treatment is done here for polylines which - when projected -
// "fall on the other side of the Earth".
// Detect an line segment within the polyline which has a very long length and split.
      ds = vtkPolyData::New();
      ca_n = vtkCellArray::New();
      static_cast<vtkPolyData *>(ds)->SetPolys(ca_n);
      ca_n->Delete();

      ca = static_cast<vtkPolyData *>(in_ds)->GetPolys();
      vtkIdType npts, *pts;

      ca->InitTraversal();
// for each polygon, change for big changes in coordinates and split lines
      for(int i =0; i < ca->GetNumberOfCells (); i++)
         {
         changeOfSigns = 0;
         ca->GetNextCell(npts, pts);
         k = npts-1;
// start from end and split if necessary
         for(int j =npts-1; j >0; j--)
           {
           newPoints->GetPoint(pts[j-1], p_pt); // previous pt
           newPoints->GetPoint(pts[j], c_pt);    // current pt
//           if(((p_pt[0] > 0) && (c_pt[0] < 0) || (p_pt[0] < 0) && (c_pt[0] > 0)))
// compute a 2d distance 
           if(sqrt((p_pt[0] - c_pt[0])*(p_pt[0] - c_pt[0]) + (p_pt[1] - c_pt[1])*(p_pt[1] - c_pt[1])) > tol )
             {
             changeOfSigns++;
             ca_n->InsertNextCell(k-j+1, &pts[j]); k = j-1;
             }
//           cerr << "pts["<< pts[j] << "] = " << x[0] << " "  << x[1] << " " << x[2] << endl;
           }
        if(changeOfSigns == 0)
          {
          ca_n->InsertNextCell(npts, pts);
          } // full polygon (original)
        else
          {
          ca_n->InsertNextCell(k+1, pts); // what is left-over after all the splits
          }
         }
    break;
    default:
      debug4 << "not supported for this grid type"  <<endl;
      return in_dr;
  }

  ds->SetPoints(newPoints);

  if(do_type != VTK_POLY_DATA)
  {
    debug4 << "dims = " << dims[0] << " x "  << dims[1] << " x " << dims[2] <<  endl;
  }

  if(do_type == VTK_RECTILINEAR_GRID)
    inPts->Delete();

  newPoints->Delete();
  proj->Delete();
  geoXform->Delete();

  avtDataRepresentation *out_dr = new avtDataRepresentation(ds,
      in_dr->GetDomain(), in_dr->GetLabel());

  ds->Delete();

  return out_dr;
}

// ****************************************************************************
// Method: avtCartographicProjectionFilter::UpdateDataObjectInfo
//
// Purpose:
//   Update the data object information.
//
// Programmer: Jean M. Favre
// Creation:   Mon Apr  7 16:11:17 PDT 2014
//
// Modifications:
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtCartographicProjectionFilter::UpdateDataObjectInfo(void)
{
  avtDataAttributes &dataAtts = GetOutput()->GetInfo().GetAttributes();
  dataAtts.SetXUnits("");
  dataAtts.SetXLabel("Longitude");
  dataAtts.SetYUnits("");
  dataAtts.SetYLabel("Latitude");
  dataAtts.SetSpatialDimension(2);

  dataAtts.AddFilterMetaData("CartographicProjection");
}
