/*****************************************************************************
*
* Copyright (c) 2000 - 2018, 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.
*
*****************************************************************************/

// ************************************************************************* //
//  File: avtLineoutFilter.C
// ************************************************************************* //

#include <avtLineoutFilter.h>

#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCellIntersections.h>
#include <vtkDataSet.h>
#include <vtkDataSetRemoveGhostCells.h>
#include <vtkGenericCell.h>
#include <vtkIdList.h>
#include <vtkInformation.h>
#include <vtkLineoutFilter.h>
#include <vtkMath.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkStreamingDemandDrivenPipeline.h>
#include <vtkUnsignedIntArray.h>
#include <vtkVisItUtility.h>

#include <avtDatasetExaminer.h>
#include <avtExtents.h>
#include <avtMetaData.h>
#include <avtVector.h>
#include <avtIntervalTree.h>
#include <ImproperUseException.h>
#include <InvalidDimensionsException.h>
#include <UnexpectedValueException.h>
#include <DebugStream.h>

#include <maptypes.h>
#include <set>
#include <vector>

struct IntersectionPoint
{
    double x[3];
};

struct CellInfo
{
    int origCell;
    int origDomain;
    intVector currCell;
    std::vector<IntersectionPoint> isect;
} ;


// ****************************************************************************
//  Method: avtLineoutFilter constructor
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Thu Apr 25 16:01:28 PST 2002
//
//  Modifications:
//    Kathleen Bonnell, Wed Oct 20 17:20:38 PDT 2004
//    Initialize useOriginalCells.
//
//    Kathleen Bonnell, Mon Aug 14 18:12:09 PDT 2006
//    Initialize ndims.
//
//    Hank Childs, Thu Jan 24 09:44:45 PST 2008
//    Initialized new data members.
//
//    Hank Childs, Fri Jan 25 09:59:29 PST 2008
//    Remove ignoreGlobal, which was unused.
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
// ****************************************************************************

avtLineoutFilter::avtLineoutFilter()
{
    OverrideOriginalSpatialExtents();
    useOriginalCells = false;
    ndims = 2;

    point1[0] = 0.;
    point1[1] = 0.;
    point1[2] = 0.;
    point2[0] = 1.;
    point2[1] = 1.;
    point2[2] = 0.;
    samplingOn   = false;
    numberOfSamplePoints = 50;
}


// ****************************************************************************
//  Method: avtLineoutFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the Lineout filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Thu Apr 25 16:01:28 PST 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Jul 12 17:28:31 PDT 2002
//    No longer send YScale to vtkLineoutFilter, it is not needed.
//
//    Hank Childs, Tue Sep 10 16:46:57 PDT 2002
//    Re-work memory management.
//
//    Kathleen Bonnell, Tue Dec 23 10:18:06 PST 2003
//    Set vtkLineoutFilter's UpdateGhostLevel, so that ghost levels can be
//    ignored.  Ensure output has points.
//
//    Kathleen Bonnell, Thu Jul 29 09:55:49 PDT 2004
//    Moved code to Sampling method, added call to NoSampling.
//
//    Hank Childs, Thu Jan 24 09:44:45 PST 2008
//    Make use of new data members.
//
//    Hank Childs, Fri Jan 25 09:59:29 PST 2008
//    Fix stupid mistake (missing negation of samplingOn when stripping out
//    the LineoutAttributes).
//
//    Eric Brugger, Mon Jul 21 14:09:11 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Kathleen Biagas, Wed Nov 12 08:38:55 PST 2014
//    Account for NULL out_ds, don't create out_dr in that case.
//
// ****************************************************************************

avtDataRepresentation *
avtLineoutFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set and domain number.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();
    int domain = in_dr->GetDomain();

    vtkDataSet *out_ds = NULL;
    if (!samplingOn)
        out_ds = NoSampling(in_ds, domain);
    else
        out_ds = Sampling(in_ds, domain);

    avtDataRepresentation *out_dr = NULL;
    if (out_ds != NULL)
    {
        out_dr = new avtDataRepresentation(out_ds, in_dr->GetDomain(),
                                           in_dr->GetLabel());
        out_ds->Delete();
    }

    return out_dr;
}


// ****************************************************************************
//  Method: avtLineoutFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Allows the filter to change its output's data object information, which
//      is a description of the data object.
//
//  Programmer: Kathleen Bonnell
//  Creation:   April 25, 202
//
//  Modifications:
//    Kathleen Bonnell, Tue Mar 23 08:48:33 PST 2004
//    Set X and Y axis labels.
//
//    Brad Whitlock, Thu Jul 22 17:20:05 PST 2004
//    Set the Y units.
//
//    Kathleen Bonnell, Thu Jan  6 10:34:57 PST 2005
//    Remove TRY-CATCH block in favor of testing for ValidActiveVariable.
//
// ****************************************************************************

void
avtLineoutFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetAttributes().SetTopologicalDimension(1);
    GetOutput()->GetInfo().GetAttributes().SetXLabel("Distance");
    GetOutput()->GetInfo().GetAttributes().SetYLabel("Value");
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetValidity().InvalidateSpatialMetaData();

    if (GetInput()->GetInfo().GetAttributes().ValidActiveVariable())
    {
        std::string units(GetInput()->GetInfo().GetAttributes().GetVariableUnits());
        if(units != "")
            GetOutput()->GetInfo().GetAttributes().SetYUnits(units);
    }
}


// ****************************************************************************
//  Method: avtLineoutFilter::VerifyInput
//
//  Purpose:
//      Verifies that the input is 2D data, throws an exception if not.
//
//  Programmer: Kathleen Bonnell
//  Creation:   April 26, 2002
//
//  Modifications:
//    Kathleen Bonnell, Mon Dec 23 11:50:54 PST 2002
//    Modified to test topological dimenions, so that lineouts of point-var
//    or lines data will be disallowed.  Lineouts of 3d now allowed.
//
//    Kathleen Bonnell, Mon Aug 14 18:12:09 PDT 2006
//    Grab spatial dimension.
//
// ****************************************************************************

void
avtLineoutFilter::VerifyInput(void)
{
    if  (GetInput()->GetInfo().GetAttributes().GetTopologicalDimension() < 2)
    {
        EXCEPTION2(InvalidDimensionsException, "Lineout", "2D or 3D");
    }
    ndims = GetInput()->GetInfo().GetAttributes().GetSpatialDimension();
}




// ****************************************************************************
//  Method: avtLineoutFilter::ModifyContract
//
//  Purpose:
//      Calculates the restriction on the meta-data and the line endpoints.
//
//  Arguments:
//      spec    The current pipeline specification.
//
//  Returns:    The new specification.
//
//  Programmer: Kathleen Bonnell
//  Creation:   December 19, 2003
//
//  Modifications:
//    Kathleen Bonnell, Wed Oct 20 17:20:38 PDT 2004
//    Set useOriginalCells.
//
//   Kathleen Bonnell, Mon Aug 14 16:40:30 PDT 2006
//   API change for avtIntervalTree.
//
//   Hank Childs, Thu Jan 24 09:44:45 PST 2008
//   Make use of new data members.
//
//   Hank Childs, Thu Dec 30 22:42:23 PST 2010
//   Add support for Lineout when it creates the output variable.
//
//   Brad Whitlock, Wed Jul 27 13:20:53 PDT 2011
//   Fix trailing slash problem.
//
//   Kathleen Biagas, Mon Feb  9 15:34:07 PST 2015
//   Do not use interval tree if spatial meta data was invalidated.
//
// ****************************************************************************


avtContract_p
avtLineoutFilter::ModifyContract(avtContract_p in_contract)
{
    avtContract_p rv;
    const char *opLineout = "operators/Lineout/";
    if (strncmp(pipelineVariable, opLineout, strlen(opLineout)) == 0)
    {
        const char *var = pipelineVariable + strlen(opLineout);
        avtDataRequest_p dr = new avtDataRequest(in_contract->GetDataRequest(), var);
        rv = new avtContract(in_contract, dr);
    }
    else
    {
        rv = new avtContract(in_contract);
    }

    useOriginalCells = false;
    if (!GetInput()->GetInfo().GetValidity().GetZonesPreserved())
    {
        rv->GetDataRequest()->TurnZoneNumbersOn();
        useOriginalCells = true;
        return rv;
    }

    if (!GetInput()->GetInfo().GetValidity().GetSpatialMetaDataPreserved())
    {
        return rv;
    }

    //
    // Get the interval tree.
    //
    avtIntervalTree *it = GetMetaData()->GetSpatialExtents();
    if (it == NULL)
    {
        return rv;
    }

    double rayDir[3] = {point2[0]-point1[0], point2[1]-point1[1],
                        point2[2]-point1[2]};

    intVector domains;
    it->GetElementsList(point1, rayDir, domains);
    rv->GetDataRequest()->GetRestriction()->RestrictDomains(domains);

    return rv;
}


// ****************************************************************************
//  Method: avtLineoutFilter::PostExecute
//
//  Purpose:
//      Cleans up after the execution.  This manages extents.
//
//  Programmer: Kathleen Bonnell
//  Creation:   January 14, 2004
//
//  Modifications:
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
// ****************************************************************************

void
avtLineoutFilter::PostExecute(void)
{
    avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();
    outAtts.GetOriginalSpatialExtents()->Clear();
    outAtts.GetDesiredSpatialExtents()->Clear();

    double bounds[6];
    avtDataset_p ds = GetTypedOutput();
    avtDatasetExaminer::GetSpatialExtents(ds, bounds);
    outAtts.GetThisProcsOriginalSpatialExtents()->Set(bounds);
}


// ****************************************************************************
//  Method: avtLineoutFilter::CreateRGrid
//
//  Purpose:
//    Creates a vtkPolyData (vertices) from  the passed information.
//    The x-coordinate is determined using the distance of each point (pts)
//    from the origin of the line (pt1).  The y-coordinate is determined by
//    the scalar value at each intersected cell (cells).  If the scalars
//    are point-centered, then they are averaged for each cell.
//
//  Arguments:
//    ds        The input dataset.
//    pt1       The origin of the line (used to calculate distances).
//    pt2       The endpoint of the line (used to calculate distances).
//    pts       Cell center points.
//    cells     A list of intersected cells.
//
//  Returns:    The output poly data.
//
//  Programmer: Kathleen Bonnell
//  Creation:   July 27, 2004
//
//  Modifications:
//    Kathleen Bonnell, Wed Oct 20 17:20:38 PDT 2004
//    Added pt2 arg.  Pts being passed are now cell centers, so find
//    closest point along the line defined by pt1 and pt2, and use it to
//    calculate distance.
//
//    Kathleen Bonnell, Mon Jul 31 10:15:00 PDT 2006
//    Create RectilinearGrid instead of PolyData for curve representation.
//
//    Kathleen Bonnell, Mon Sep 11 16:47:08 PDT 2006
//    Removed calculation of ClosestPointOnLine, no longer using cell centers.
//
//    Kathleen Bonnell, Thu Mar  6 09:07:33 PST 2008
//    Removed unused variable.
//
//    Brad Whitlock, Mon Jan  7 14:15:45 PST 2013
//    Use VTK_FLOAT as the minimum precision for the rgrid and then adjust
//    upward to VTK_DOUBLE if the input coordinates or scalar use that type.
//    This fixes a bug that could cause the rgrid to be created using int
//    coordinates if the field was int. That becomes a problem with small
//    coordinates!
//
//    Kathleen Biagas, Mon Feb 26 09:24:23 PST 2018
//    Ensure output has a name attached to the PD array.
//
// ****************************************************************************

vtkRectilinearGrid *
avtLineoutFilter::CreateRGrid(vtkDataSet *ds, double *pt1, double *pt2,
                              vtkPoints *pts, vtkIdList *cells)
{
    bool pointData = true;
    vtkDataArray *scalars = ds->GetPointData()->GetScalars();

    if (scalars == NULL)
    {
        pointData = false;
        scalars = ds->GetCellData()->GetScalars();
        if (scalars == NULL)
        {
            return NULL;
        }
    }

    int npts = pts->GetNumberOfPoints();
    int ptsType = pts->GetDataType();
    int sType = scalars->GetDataType();
    int dType = VTK_FLOAT;
    if(ptsType == VTK_DOUBLE || sType == VTK_DOUBLE)
        dType = VTK_DOUBLE;
    vtkRectilinearGrid *rgrid = vtkVisItUtility::Create1DRGrid(0, dType);
    vtkDataArray *outXC = rgrid->GetXCoordinates();
    vtkDataArray *outVal = outXC->NewInstance();
    outVal->SetName(scalars->GetName());
    rgrid->GetPointData()->SetScalars(outVal);
    outVal->Delete();

    vtkIdList *ptIds = vtkIdList::New();
    double currentPoint[3];
    double sum = 0.;
    int i, j;
    double newX = 0., oldX = -1., newVal = 0.;
    bool requiresSort = false;
    for (i = 0; i < npts; i++)
    {
        pts->GetPoint(i, currentPoint);
        newX = sqrt(vtkMath::Distance2BetweenPoints(pt1, currentPoint));
        if (newX < oldX)
        {
            requiresSort = true;
        }
        oldX = newX;
        if (pointData)
        {
            sum = 0;
            ds->GetCellPoints(cells->GetId(i), ptIds);
            int numCellPts = ptIds->GetNumberOfIds();
            for (j = 0; j < numCellPts; j++)
                sum += scalars->GetTuple1(ptIds->GetId(j));
            if (numCellPts > 0)
               sum /= (double) numCellPts;

            newVal = sum;
        }
        else
        {
            newVal = scalars->GetTuple1(cells->GetId(i));
        }
        outXC->InsertNextTuple1(newX);
        outVal->InsertNextTuple1(newVal);
    }
    ptIds->Delete();
    rgrid->SetDimensions(outXC->GetNumberOfTuples(), 1 , 1);

    if (requiresSort)
    {
        vtkDataArray *sortedXC = outXC->NewInstance();
        vtkDataArray *sortedVal = outVal->NewInstance();
        sortedVal->SetName(outVal->GetName());
        DoubleIntMap sortedIds;
        double x;
        for (i = 0; i < outXC->GetNumberOfTuples(); i++)
        {
            x = outXC->GetTuple1(i);
            sortedIds.insert(DoubleIntMap::value_type(x, i));
        }
        DoubleIntMap::iterator it;
        for (it = sortedIds.begin(); it != sortedIds.end(); it++)
        {
            sortedXC->InsertNextTuple1(outXC->GetTuple1((*it).second));
            sortedVal->InsertNextTuple1(outVal->GetTuple1((*it).second));
        }
        rgrid->SetXCoordinates(sortedXC);
        rgrid->GetPointData()->SetScalars(sortedVal);
        sortedXC->Delete();
        sortedVal->Delete();
        rgrid->SetDimensions(sortedXC->GetNumberOfTuples(), 1, 1);
    }

    return rgrid;
}


// ****************************************************************************
//  Method: avtLineoutFilter::NoSampling
//
//  Purpose:
//    Peforms a lineout by intersecting cells.
//
//  Arguments:
//    in_ds      The input dataset.
//    domain     The domain number.
//
//  Returns:       The output dataset.
//
//  Programmer: Kathleen Bonnell
//  Creation:   July 27, 2004
//
//  Modifications:
//    Kathleen Bonnell, Wed Oct 20 17:20:38 PDT 2004
//    Use different CreatePolys method if we must use Original Cell Numbers.
//    Tell the locator to ignore ghost zones.
//
//    Hank Childs, Fri Mar 11 17:04:37 PST 2005
//    Fix memory leak.
//
//    Kathleen Bonnell, Thu May 19 11:34:05 PDT 2005
//    Added logic to determine if the cell-centers or the intersection points
//    should be used in creating the curve.
//
//    Kathleen Bonnell, Mon Jul 31 10:15:00 PDT 2006
//    Create RectilinearGrid instead of PolyData for the curves.
//
//    Kathleen Bonnell, Mon Aug 14 18:12:09 PDT 2006
//    Use avtIntervalTree and vtkCellIntersections(new) instead of
//    vtkVisItCellLocator, for better accuracy.
//
//    Kathleen Bonnell, Mon Aug 21 13:34:18 PDT 2006
//    Tell IntervalTree to NOT do communication (via arg to constructor).
//
//    Kathleen Bonnell, Mon Sep 11 16:47:08 PDT 2006
//    Removed calculation of cell centers.
//
//    Hank Childs, Thu Jan 24 09:44:45 PST 2008
//    Make use of new data members.
//
//    Kathleen Bonnell, Thu Jun 11 08:29:51 PDT 2009
//    Calculate the min side-length of a cell, to be used as a tolerance when
//    determining if PointsEqual.
//
//    Brad Whitlock, Tue Mar 26 10:11:39 PDT 2013
//    Make sure we have original cells. If we don't, back up to another lineout
//    method instead of failing.
//
//    Eric Brugger, Mon Jul 21 14:09:11 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

double
GetMinRange(double *b, int ndims)
{
    double r1 = 1e38;
    double r2 = 1e38;
    double r3 = 1e38;
    r1 = b[1] - b[0];
    r2 = b[3] - b[2];
    if (ndims == 3)
        r3 = b[5] - b[4];
    return (r1<r2 ? (r1 <r3 ? r1 : r3) : (r2 < r3 ? r2 : r3));
}

vtkDataSet *
avtLineoutFilter::NoSampling(vtkDataSet *in_ds, int domain)
{
    bool rgrid = in_ds->GetDataObjectType() == VTK_RECTILINEAR_GRID;
    vtkDataSet *rv = NULL;

    //
    // Use interval tree to find intersected cells
    //

    int ncells = in_ds->GetNumberOfCells();
    avtIntervalTree itree(ncells, ndims, false);
    double bounds[6];
    double tol = 1e38;
    for (int i = 0; i < ncells; i++)
    {
        in_ds->GetCellBounds(i, bounds);
        itree.AddElement(i, bounds);
        double min = GetMinRange(bounds, ndims);
        if (min < tol)
            tol = min;
    }
    tol = tol/2.;
    itree.Calculate();

    intVector isectedCells;
    doubleVector isectedPts;
    intVector lineCells;
    doubleVector linePts;
    itree.GetElementsListFromLine(point1, point2, lineCells, linePts, &tol);

    if (lineCells.size() == 0)
    {
        debug5 << "avtIntervalTree returned NO intersected cells for domain "
               << domain << "." << endl;
        return NULL;
    }

    vtkUnsignedCharArray *ghosts =
       (vtkUnsignedCharArray *)in_ds->GetCellData()->GetArray("avtGhostZones");

    vtkCellIntersections *cellIntersections = NULL;
    vtkGenericCell *cell = NULL;
    double isect[3], p1[3], t;
    if (!rgrid)
    {
        cellIntersections = vtkCellIntersections::New();
        cellIntersections->SetTestCoPlanar(true);
        cell = vtkGenericCell::New();
    }

    for (size_t i = 0; i < lineCells.size(); i++)
    {
        //
        // Want to skip ghost zones!!
        //
        if (ghosts && ghosts->GetComponent(lineCells[i], 0))
        {
            continue;
        }

        bool doMore = true;
        if (!rgrid)
        {
            //
            // do a bit more checking if not a rectlinearGrid, as the interval
            // tree only calculated intersections with a cell's bounding box.
            //

            in_ds->GetCell(lineCells[i], cell);
            double isect2[3];
            // Check both directions of the line, because we don't want
            // to intersect only a node.
            if (cellIntersections->CellIntersectWithLine(cell, point1, point2,
                                        t, isect) &&
                cellIntersections->CellIntersectWithLine(cell, point2, point1,
                                        t, isect2))
            {
                if (vtkVisItUtility::PointsEqual(isect, isect2, &tol))
                {
                    // Discard single-point intersections that are equivalent
                    // to an endpoint.
                    doMore = false;
                    if (vtkVisItUtility::CellContainsPoint(cell, point1))
                    {
                        if (!vtkVisItUtility::PointsEqual(point1, isect, &tol))
                        {
                            p1[0] = point1[0];
                            p1[1] = point1[1];
                            p1[2] = point1[2];
                            doMore = true;
                        }
                    }
                    else if (vtkVisItUtility::CellContainsPoint(cell, point2))
                    {
                        if (!vtkVisItUtility::PointsEqual(point2, isect, &tol))
                        {
                            p1[0] = point2[0];
                            p1[1] = point2[1];
                            p1[2] = point2[2];
                            doMore = true;
                        }
                    }
                } // points equal
                else
                {
                    p1[0] = isect[0];
                    p1[1] = isect[1];
                    p1[2] = isect[2];
                    doMore = true;
                }
            } // isects both directions
            else
            {
                doMore = false;
            }
        } // rgrid
        else
        {
            p1[0] = linePts[i*3];
            p1[1] = linePts[i*3+1];
            p1[2] = linePts[i*3+2];
        }
        if (doMore)
        {
            bool dupFound = false;
            for (size_t j = 0; j < isectedCells.size() && !dupFound; j++)
            {
                double p2[3];
                p2[0] = isectedPts[j*3];
                p2[1] = isectedPts[j*3+1];
                p2[2] = isectedPts[j*3+2];
                if (vtkVisItUtility::PointsEqual(p1, p2, &tol))
                {
                    dupFound = true;
                    if (lineCells[i] < isectedCells[j])
                        isectedCells[j] = lineCells[i];
                }
            }
            if (!dupFound)
            {
                isectedCells.push_back(lineCells[i]);
                isectedPts.push_back(p1[0]);
                isectedPts.push_back(p1[1]);
                isectedPts.push_back(p1[2]);
            }
        }
    }

    vtkIdList *cells = vtkIdList::New();
    vtkPoints *pts = vtkVisItUtility::NewPoints(in_ds);
    double cpt[3];
    for (size_t i = 0; i < isectedCells.size(); i++)
    {
        cpt[0] = isectedPts[i*3];
        cpt[1] = isectedPts[i*3+1];
        cpt[2] = isectedPts[i*3+2];
        pts->InsertNextPoint(cpt);
        cells->InsertNextId(isectedCells[i]);
    }

    // If we want original cells, be sure that we have them since it's
    // possible that filters before lineout might not be capable of preserving
    // the original cells.
    if (useOriginalCells &&
        in_ds->GetCellData()->GetArray("avtOriginalCellNumbers") != NULL)
    {
        rv = CreateRGridFromOrigCells(in_ds, point1, point2, pts, cells);
    }
    else
    {
        rv = CreateRGrid(in_ds, point1, point2, pts, cells);
    }

    if (rv->GetNumberOfCells() == 0 ||
        rv->GetNumberOfPoints() == 0)
    {
        debug5 << "avtIntervalTree returned empty DS for domain "
                   << domain << "." << endl;
        rv->Delete();
        rv = NULL;
    }

    cells->Delete();
    pts->Delete();

    if (!rgrid)
    {
        cellIntersections->Delete();
        cell->Delete();
    }

    return rv;
}


// ****************************************************************************
//  Method: avtLineoutFilter::Sampling
//
//  Purpose:
//    Executes the lineout using a sampling method.
//
//  Arguments:
//    in_ds      The input dataset.
//    domain     The domain number.
//
//  Returns:       The output dataset.
//
//  Programmer: kbonnell -- generated by xml2info
//  Creation:   Thu Apr 25 16:01:28 PST 2002
//
//  Modifications:
//    Kathleen Bonnell, Fri Jul 12 17:28:31 PDT 2002
//    No longer send YScale to vtkLineoutFilter, it is not needed.
//
//    Hank Childs, Tue Sep 10 16:46:57 PDT 2002
//    Re-work memory management.
//
//    Kathleen Bonnell, Tue Dec 23 10:18:06 PST 2003
//    Set vtkLineoutFilter's UpdateGhostLevel, so that ghost levels can be
//    ignored.  Ensure output has points.
//
//    Kathleen Bonnell, Tue Jul 27 10:18:14 PDT 2004
//    Moved from 'ExecueData' method.
//
//    Hank Childs, Wed Sep  8 19:57:21 PDT 2004
//    Remove ghost zones before doing a lineout.  This is because the
//    vtkLineoutFilter is a bit touchy about ghost zone values.  If you have
//    values > 1, then it can do interpolations to nodal data that can
//    mistakenly identify real zones as ghost.
//
//    Hank Childs, Thu Jan  4 09:51:34 PST 2007
//    Manually force an update to the ghost filter.  If we don't do this,
//    then the lineout filter gets bad data.  I believe this is because
//    the data set remove ghost cells filter doesn't know what it's real
//    output is until it updates.
//
//    Kathleen Biagas, Fri Jan 25 16:04:46 PST 2013
//    Call Update on the filter, not the data object.
//
//    Eric Brugger, Mon Jul 21 14:09:11 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Kathleen Biagas, Mon Aug 15 14:09:55 PDT 2016
//    VTK-8, API for updating GhostLevel changed.
//
// ****************************************************************************

vtkDataSet *
avtLineoutFilter::Sampling(vtkDataSet *in_ds, int domain)
{
    vtkDataSetRemoveGhostCells *ghosts = vtkDataSetRemoveGhostCells::New();
    ghosts->SetInputData(in_ds);

    vtkLineoutFilter *filter = vtkLineoutFilter::New();

    filter->SetInputConnection(ghosts->GetOutputPort());
    filter->SetPoint1(point1);
    filter->SetPoint2(point2);
    filter->SetNumberOfSamplePoints(numberOfSamplePoints);
    filter->GetInformation()->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), 0);
    filter->Update();
    vtkPolyData *outPolys = filter->GetOutput();
    outPolys->Register(NULL);

    vtkDataSet *rv = outPolys;
    if (outPolys->GetNumberOfCells() == 0 ||
        outPolys->GetNumberOfPoints() == 0)
    {
        debug5 << "vtkLineoutFilter returned empty DS for domain "
               << domain << "." << endl;
        rv->Delete();
        rv = NULL;
    }

    filter->Delete();
    ghosts->Delete();

    return rv;
}

// ****************************************************************************
//  Method: avtLineoutFilter::CreatePolysFromOrigCells
//
//  Purpose:
//    Creates a vtkPolyData (vertices) from  the passed information.
//    The x-coordinate is determined using the distance of each point (pts)
//    from the origin of the line (pt1).  The y-coordinate is determined by
//    the scalar value at each intersected cell (cells).  If the scalars
//    are point-centered, then they are averaged for each cell.
//
//    This method uses the avtOriginalCellNumbers array to combine values
//    that come from cells that share the same original cell. (e.g. a Slice,
//    which will triangulate quads -- we only want 1 intersection point
//    for the quad, not 1 for each triangle.
//
//  Arguments:
//    ds        The input dataset.
//    pt1       The origin of the line (used to calculate distances).
//    pts       Intersection points along the line.
//    cells     A list of intersected cells.
//
//  Returns:    The output poly data.
//
//  Programmer: Kathleen Bonnell
//  Creation:   July 27, 2004
//
//  Modifications:
//    Brad Whitlock, Wed Nov 3 10:16:32 PDT 2004
//    Fixed on win32.
//
//    Kathleen Bonnell, Mon Jul 31 10:15:00 PDT 2006
//    Create RectilinearGrid instead of PolyData for the curves.
//
//    Kathleen Bonnell, Mon Aug 14 18:12:09 PDT 2006
//    Remove unnecessary if(origCells) statements.
//
//    Kathleen Bonnell, Mon Sep 11 16:47:08 PDT 2006
//    Removed calculation of ClosestPointOnLine, no longer using cell centers.
//
//    Kathleen Bonnell, Thu Mar  6 09:07:33 PST 2008
//    Removed unused variable.
//
//    Jeremy Meredith, Wed Sep  3 10:28:07 EDT 2008
//    Renamed Point struct to IntersectionPoint to avoid namespace conflict
//    with a file in common/state.
//
//    Brad Whitlock, Mon Jan  7 14:15:45 PST 2013
//    Use VTK_FLOAT as the minimum precision for the rgrid and then adjust
//    upward to VTK_DOUBLE if the input coordinates or scalar use that type.
//    This fixes a bug that could cause the rgrid to be created using int
//    coordinates if the field was int. That becomes a problem with small
//    coordinates!
//
// ****************************************************************************

vtkRectilinearGrid *
avtLineoutFilter::CreateRGridFromOrigCells(vtkDataSet *ds, double *pt1,
    double *pt2, vtkPoints *pts, vtkIdList *cells)
{
    bool pointData = true;
    vtkDataArray *scalars = ds->GetPointData()->GetScalars();
    if (scalars == NULL)
    {
        pointData = false;
        scalars = ds->GetCellData()->GetScalars();
        if (scalars == NULL)
        {
            return NULL;
        }
    }

    size_t npts = pts->GetNumberOfPoints();
    vtkUnsignedIntArray *origCells = vtkUnsignedIntArray::SafeDownCast(
        ds->GetCellData()->GetArray("avtOriginalCellNumbers"));

    if (!origCells)
    {
       EXCEPTION2(UnexpectedValueException, "avtOriginalCellNumbers",
                  "nada.  Internal error.");
    }
    int currentCell;
    int origCell;
    int origDomain;
    double center[3];
    std::vector<CellInfo> cellInfoList;
    bool dup = false;
    for (size_t i = 0; i < npts; i++)
    {
        currentCell = cells->GetId(i);
        pts->GetPoint(i, center);

        origCell = (int)origCells->GetComponent(currentCell, 1);
        origDomain = (int)origCells->GetComponent(currentCell, 0);
        dup = false;
        for (size_t j = 0; j < cellInfoList.size() && !dup; j++)
        {
            if ((origCell == cellInfoList[j].origCell) &&
                (origDomain == cellInfoList[j].origDomain))
            {
                dup = true;
                cellInfoList[j].currCell.push_back(currentCell);
                IntersectionPoint p;
                p.x[0]  = center[0];
                p.x[1]  = center[1];
                p.x[2]  = center[2];
                cellInfoList[j].isect.push_back(p);
            }
        }
        if (!dup)
        {
            CellInfo a;
            a.origCell = origCell;
            a.origDomain = origDomain;
            a.currCell.push_back(currentCell);
            IntersectionPoint p;
            p.x[0]  = center[0];
            p.x[1]  = center[1];
            p.x[2]  = center[2];
            a.isect.push_back(p);
            cellInfoList.push_back(a);
        }
    }

    double newX = 0., oldX = -1, newVal = 0.;

    int ptsType = pts->GetDataType();
    int sType = scalars->GetDataType();
    int dType = VTK_FLOAT;
    if(ptsType == VTK_DOUBLE || sType == VTK_DOUBLE)
        dType = VTK_DOUBLE;
    vtkRectilinearGrid *rgrid = vtkVisItUtility::Create1DRGrid(0, dType);

    vtkDataArray *outXC = rgrid->GetXCoordinates();
    vtkDataArray *outVal = outXC->NewInstance();
    rgrid->GetPointData()->SetScalars(outVal);
    outVal->Delete();
    vtkIdList *ptIds = vtkIdList::New();
    double sum = 0.;
    bool requiresSort = false;
    CellInfo cellInfo;

    for (size_t i = 0; i < cellInfoList.size(); i++)
    {
        cellInfo = cellInfoList[i];
        size_t nDups = cellInfo.currCell.size();

        if (nDups == 1)
        {
        newX = sqrt(vtkMath::Distance2BetweenPoints(pt1, cellInfo.isect[0].x));
        }
        else
        {
            sum = 0.;
            for (size_t j = 0; j < nDups; j++)
            {
              sum += vtkMath::Distance2BetweenPoints(pt1, cellInfo.isect[j].x);
            }
            newX = sqrt(sum/(double)nDups);
        }
        if (newX < oldX)
        {
            requiresSort = true;
        }
        oldX = newX;
        if (pointData)
        {
            if (nDups == 1)
            {
                sum = 0;
                ds->GetCellPoints(cellInfoList[i].currCell[0], ptIds);
                vtkIdType numCellPts = ptIds->GetNumberOfIds();
                for (vtkIdType j = 0; j < numCellPts; j++)
                    sum += scalars->GetTuple1(ptIds->GetId(j));
                if (numCellPts > 0)
                   sum /= (double) numCellPts;
                newVal = sum;
            }
            else
            {
                sum = 0;
                std::set<int> uniquePts;
                for (size_t j = 0; j < nDups; j++)
                {
                    ds->GetCellPoints(cellInfoList[i].currCell[j], ptIds);
                    vtkIdType numCellPts = ptIds->GetNumberOfIds();
                    for (vtkIdType k = 0; k < numCellPts; k++)
                    {
                        if (uniquePts.count(ptIds->GetId(k)) == 0)
                        {
                            uniquePts.insert(ptIds->GetId(k));
                            sum += scalars->GetTuple1(ptIds->GetId(j));
                        }
                    }
                }
                if (uniquePts.size() > 0)
                   sum /= (double) uniquePts.size();
                newVal = sum;
            }
        }
        else
        {
            if (nDups == 1)
            {
                newVal = scalars->GetTuple1(cellInfoList[i].currCell[0]);
            }
            else
            {
                newVal = 0.;
                for (size_t j = 0; j < nDups; j++)
                    newVal += scalars->GetTuple1(cellInfoList[i].currCell[j]);
                newVal /= (double) cellInfoList[i].currCell.size();
            }
        }
        outXC->InsertNextTuple1(newX);
        outVal->InsertNextTuple1(newVal);
    }
    ptIds->Delete();
    rgrid->SetDimensions(outXC->GetNumberOfTuples(), 1, 1);

    if (requiresSort)
    {
        vtkDataArray *sortedXC = outXC->NewInstance();
        vtkDataArray *sortedVal = outVal->NewInstance();
        DoubleIntMap sortedIds;
        double x;
        for (int i = 0; i < outXC->GetNumberOfTuples(); i++)
        {
            x = outXC->GetTuple1(i);
            sortedIds.insert(DoubleIntMap::value_type(x, i));
        }
        DoubleIntMap::iterator it;
        for (it = sortedIds.begin(); it != sortedIds.end(); it++)
        {
            sortedXC->InsertNextTuple1(outXC->GetTuple1((*it).second));
            sortedVal->InsertNextTuple1(outVal->GetTuple1((*it).second));
        }
        rgrid->SetXCoordinates(sortedXC);
        rgrid->GetPointData()->SetScalars(sortedVal);
        rgrid->SetDimensions(sortedXC->GetNumberOfTuples(), 1, 1);
        sortedXC->Delete();
        sortedVal->Delete();
    }

    return rgrid;
}

