/*****************************************************************************
*
* 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: avtThresholdFilter.C
// ************************************************************************* //

#include <avtThresholdFilter.h>

#include <vtkCellData.h>
#include <vtkDataArray.h>
#include <vtkDataObject.h>
#include <vtkDataSet.h>
#include <vtkIdList.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkThreshold.h>
#include <vtkUnstructuredGrid.h>

#include <vtkVisItUtility.h>

#include <avtCallback.h>
#include <avtDataAttributes.h>
#include <avtIdentifierSelection.h>
#include <avtIntervalTree.h>
#include <avtMetaData.h>
#include <avtNamedSelection.h>
#include <avtStructuredMeshChunker.h>
#include <avtDataRangeSelection.h>

#include <DebugStream.h>
#include <ImproperUseException.h>
#include <NoDefaultVariableException.h>


// ****************************************************************************
//  Method: avtThreshold constructor
//
//  Programmer: Hank Childs
//  Creation:   March 19, 2005
//     
//  Modifications:
//
//    Hank Childs, Sun Mar 27 11:36:49 PST 2005
//    Don't initialize data members that have now been pushed into the base
//    class.
//
// ****************************************************************************

avtThresholdFilter::avtThresholdFilter()
{
    activeVarName = std::string("<unused>");
}


// ****************************************************************************
//  Method:  avtThresholdFilter::Create
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Tue Oct 23 16:38:18 PST 2001
//
// ****************************************************************************

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


// ****************************************************************************
//  Method:      avtThresholdFilter::SetAtts
//
//  Purpose: Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Tue Oct 23 16:38:18 PST 2001
//
//  Modifications:
//
//    Hank Childs, Thu Oct 25 09:17:54 PDT 2001
//    Register the active variable.
//
//    Mark Blair, Tue Mar  7 13:25:00 PST 2006
//    No longer recognizes the "default" variable.
//
//    Mark Blair, Tue Aug  8 17:47:00 PDT 2006
//    Now accommodates an empty list of threshold variables; does pass-through.
//
//    Mark Blair, Thu Sep 28 12:07:05 PDT 2006
//    Checks attributes for consistency.
//
//    Mark Blair, Wed Oct  4 17:45:48 PDT 2006
//    Makes better use of the "active variable".
//
//    Mark Blair, Tue Mar 13 19:51:29 PDT 2007
//    Now forces attribute consistency if attributes are inconsistent.
//
//    Mark Blair, Tue Apr 17 16:24:42 PDT 2007
//    Rewritten to support new Threshold GUI; no more "shown variable".
//
// ****************************************************************************

void
avtThresholdFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const ThresholdOpAttributes*)a;
    
    atts.SupplyMissingDefaultsIfAppropriate();
    
    if (!atts.AttributesAreConsistent()) atts.ForceAttributeConsistency();
    
    activeVarName = std::string("<unused>");
    
    stringVector curVarNames = atts.GetListedVarNames();

    if (curVarNames.size() > 0)
    {
        std::string defaultVariable = atts.GetDefaultVarName();
        std::string firstVarInList = curVarNames[0];
        
        if (firstVarInList == std::string("default"))
            firstVarInList = defaultVariable;

        if (firstVarInList != defaultVariable)
        {
            activeVarName = firstVarInList;
            SetActiveVariable(activeVarName.c_str());
        }
        else if (atts.GetDefaultVarIsScalar())
        {
            activeVarName = firstVarInList;
            SetActiveVariable(activeVarName.c_str());
        }
    }
}


// ****************************************************************************
//  Method: avtThresholdFilter::Equivalent
//
//  Purpose: Returns true if creating a new avtThresholdFilter with the given
//           parameters would result in an equivalent avtThresholdFilter.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Tue Oct 23 16:38:18 PST 2001
//
// ****************************************************************************

bool
avtThresholdFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(ThresholdOpAttributes*)a);
}


// ****************************************************************************
//  Method: avtThresholdFilter::ProcessOneChunk
//
//  Purpose: Sends the specified input and output through the Threshold filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: childs -- generated by xml2info
//  Creation:   Tue Oct 23 16:38:18 PST 2001
//
//  Modifications:
//
//    Hank Childs, Thu Oct 25 09:17:54 PDT 2001
//    Removed references to active variables.
//
//    Hank Childs, Tue Oct 30 09:29:13 PST 2001
//    Account for cases when we are not switching active variables.
//
//    Kathleen Bonnell, Fri Feb  8 11:03:49 PST 2002
//    vtkScalars has been deprecated in VTK 4.0, use vtkDataArray instead.
//
//    Hank Childs, Tue Sep 10 18:29:14 PDT 2002
//    Make memory management be self-contained in this routine.
//
//    Hank Childs, Sat Jun 21 09:25:34 PDT 2003
//    If we have poly data input, then we should have poly data output.
//
//    Hank Childs, Fri May  7 08:33:17 PDT 2004
//    If the variable is not a scalar, then issue a warning.
//
//    Hank Childs, Thu Jul 29 09:43:58 PDT 2004
//    Do not convert output to poly-data, since the base class will now
//    take care of this (when appropriate).
//
//    Hank Childs, Sat Mar 19 10:34:08 PST 2005
//    Use a structured mesh chunker if appropriate.
//
//    Hank Childs, Sun Mar 27 11:36:49 PST 2005
//    Renamed to ProcessOneChunk.  Modified to assume structured chunking.
//
//    Hank Childs, Sun Apr  3 13:06:01 PDT 2005
//    Fix bug with structured chunking.
//
//    Hank Childs, Tue Sep 13 09:03:12 PDT 2005
//    Add support for "PointsOnly".
//
//    Mark Blair, Tue Mar  7 13:25:00 PST 2006
//    Rewrote to support multi-variable thresholding.
//
//    Kathleen Bonnell, Tue May 16 15:49:24 PDT 2006
//    SetAttributeModeTo* deprecated, now use SetInputArrayToProcess.
//    Pass FieldData, as it is no longer done internally by vtkThreshold.
//
//    Mark Blair, Wed May 31 18:17:00 PDT 2006
//    SetInputScalars deprecated, now use SetInputArrayToProcess.
//
//    Mark Blair, Tue Aug  8 17:47:00 PDT 2006
//    Now accommodates an empty list of threshold variables; does pass-through.
//
//    Mark Blair, Thu Sep 28 12:07:05 PDT 2006
//    Checks attributes for consistency.
//
//    Mark Blair, Tue Mar 13 19:51:29 PDT 2007
//    Now forces attribute consistency if attributes are inconsistent.
//
//    Mark Blair, Tue Apr 17 16:24:42 PDT 2007
//    Rewritten to support new Threshold GUI; no more "shown variable".
//
//    Markus Glatter, Fri Aug 10 10:41:07 EDT 2007
//    Added avtDataRangeSelection entries to support contract-based filtering.
//
//    Hank Childs, Thu Sep 11 11:10:37 PDT 2008
//    Add fix for issue debugged by Gunther Weber ... reference counting can
//    go "off-by-one" when we have data selections in play.  This is because
//    the Threshold filter inherits from the structured mesh chunker, which
//    assumes the return value has an extra reference.
//
//    Kathleen Biagas, Mon Jan 28 11:09:02 PST 2013
//    Call Update on the fitler not the data object.
//
//    Eric Brugger, Wed Aug 20 17:06:55 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
//    Kevin Griffin, Thu Mar 23 08:25:21 PDT 2017
//    Add capability to threshold on multiple ranges (Feature #2646).
//
//    Kevin Griffin, Tue Feb 13 09:31:21 PDT 2018
//    Added check for incoming datasets with unknown mesh types
//    so they can be processed as point meshes. Fixes Bug #2505.
//
// ****************************************************************************

avtDataRepresentation *
avtThresholdFilter::ProcessOneChunk(avtDataRepresentation *in_dr, bool fromChunker)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    atts.SupplyMissingDefaultsIfAppropriate();

    if (!atts.AttributesAreConsistent()) atts.ForceAttributeConsistency();
    
    if (atts.GetListedVarNames().size() == 0)
    {
        in_ds->Register(NULL);
        return in_dr;
    }
    
    avtMeshType inputMeshType = GetInput()->GetInfo().GetAttributes().GetMeshType();

    if (atts.GetOutputMeshType() == ThresholdOpAttributes::PointMesh || inputMeshType == AVT_UNKNOWN_MESH)
    {
        vtkDataSet *out_ds = ThresholdToPointMesh(in_ds);

        avtDataRepresentation *out_dr = new avtDataRepresentation(out_ds,
                                                                  in_dr->GetDomain(),
                                                                  in_dr->GetLabel());
        out_ds->Delete();
        return out_dr;
    }

    if (fromChunker)
    {
        //
        // If in_ds is from the chunker, then the zones in in_ds are all
        // ones we identified that we wanted.  So just return them.
        //
        in_ds->Register(NULL);
        return in_dr;
    }
    
    vtkDataSet *curOutDataSet = in_ds;

    const stringVector curVariables    = atts.GetListedVarNames();
    const intVector    curZonePortions = atts.GetZonePortions();
    const doubleVector curLowerBounds  = atts.GetLowerBounds();
    const doubleVector curUpperBounds  = atts.GetUpperBounds();
    const stringVector curBoundsRange = atts.GetBoundsRange();
    
    const char *curVarName;
    char errMsg[1024];
    
    for (size_t curVarNum = 0; curVarNum < curVariables.size(); curVarNum++)
    {
        vtkThreshold *threshold = vtkThreshold::New();

        std::map<std::string,int>::iterator iterFind;
        bool bypassThreshold = false;
        iterFind = selIDs.find(curVariables[curVarNum]);
        if (iterFind != selIDs.end())
        {
            int selID = iterFind->second;
            if (GetInput()->GetInfo().GetAttributes().GetSelectionApplied(selID))
            {
                debug1 << "Bypassing Threshold operator because the database "
                       << "plugin claims to have applied selection on "
                       << curVariables[curVarNum] << endl;
                bypassThreshold = true;
            }
        }

        if (bypassThreshold == false)
        {
            curVarName = curVariables[curVarNum].c_str();
            
            if(atts.GetBoundsInputType() == ThresholdOpAttributes::Default || IsSimpleRange(curBoundsRange[curVarNum]))
            {
                threshold->SetInputData(curOutDataSet);
                // We registered curOutDataSet so it wouldn't be deleted.  But now that
                // we have fed it back into the threshold filter, we are done with it.
                // So decrement its reference count.
                if (curOutDataSet != in_ds)
                    curOutDataSet->Delete();
                
                threshold->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS,
                                                  vtkDataSetAttributes::SCALARS);
                
                if (curZonePortions[curVarNum] == (int)ThresholdOpAttributes::PartOfZone)
                {
                    threshold->AllScalarsOff();
                }
                else if (curZonePortions[curVarNum] == (int)ThresholdOpAttributes::EntireZone)
                {
                    threshold->AllScalarsOn();
                }
                else
                {
                    debug1 << "Invalid zone inclusion option encountered "
                    << "in Threshold operator attributes." << endl;
                    threshold->AllScalarsOff();
                }
                
                threshold->ThresholdBetween(curLowerBounds[curVarNum], curUpperBounds[curVarNum]);
                
                if (curOutDataSet->GetPointData()->GetArray(curVarName) != NULL)
                {
                    threshold->SetInputArrayToProcess(
                                                      0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, curVarName);
                }
                else if (curOutDataSet->GetCellData()->GetArray(curVarName) != NULL)
                {
                    threshold->SetInputArrayToProcess(
                                                      0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_CELLS, curVarName);
                }
                else
                {
                    threshold->Delete();
                    
                    sprintf (errMsg, "Data for variable \"%s\" is not currently available.", curVarName);
                    debug1 << errMsg << endl;
                    EXCEPTION1(VisItException, errMsg);
                }
                
                threshold->Update();
                curOutDataSet = threshold->GetOutput();
            }
            else
            {
                curOutDataSet = ThresholdOnRanges(curOutDataSet,
                                                  threshold,
                                                  curVarName,
                                                  curBoundsRange[curVarNum],
                                                  curZonePortions[curVarNum]);
            }
        }

        if (curOutDataSet->GetNumberOfCells() <= 0)
        {
            curOutDataSet = NULL;
            threshold->Delete();

            break;
        }

        if (bypassThreshold == false)
        {
            // We want this to stay in scope even after we delete the threshold
            // filter.  Register it.  If this goes on to another iteration, we'll
            // delete it after we connect it with the threshold filter again.
            curOutDataSet->Register(NULL);
        }

        threshold->Delete();
    }

    if (curOutDataSet == in_ds)
    {
        // The curOutDataSet equals the in_ds, meaning
        // that data selections allowed us to bypass any thresholding.
        // The structured mesh chunker will automatically decrement the
        // reference count.  So we need to add one.
        curOutDataSet->Register(NULL);
    }
    else if (curOutDataSet != NULL)
    {
        curOutDataSet->GetFieldData()->PassData(in_ds->GetFieldData());
    }
    
    avtDataRepresentation *out_dr = new avtDataRepresentation(curOutDataSet,
        in_dr->GetDomain(), in_dr->GetLabel());

    if (curOutDataSet != NULL)
        curOutDataSet->Delete();

    return out_dr;
}

// ****************************************************************************
//  Method: avtThresholdFilter::IsSimpleRange
//
//  Purpose: Determine if a range string consists of only one range or is
//           empty.
//
//  Arguments:
//      range           The range string to examine.
//
//  Returns: true if there is only one or no range in range string, otherwise
//           false
//
//  Programmer: Kevin Griffin
//  Creation:   Thu Mar 23 08:25:21 PDT 2017
//
//  Modifications:
// ****************************************************************************

bool
avtThresholdFilter::IsSimpleRange(const std::string range)
{
    if(range.empty())
    {
        return true;
    }
        
    
    std::size_t pos = range.find(",");
    
    return (pos == std::string::npos);
}

// ****************************************************************************
//  Method: avtThresholdFilter::CheckForMinMax
//
//  Purpose: Check the string for min or max.
//
//  Arguments:
//      str     The low or high range string to check
//
//  Returns: If the string contains 'min' or 'max' it will replace the string
//           with -1e+37 or 1e+37, respectively. Otherwise the original string
//           is returned.
//
//  Programmer: Kevin Griffin
//  Creation:   Thu Mar 23 08:25:21 PDT 2017
//
//  Modifications:
// ****************************************************************************
std::string
avtThresholdFilter::CheckForMinMax(const std::string str)
{
    size_t minPos = str.find("min");
    size_t maxPos = str.find("max");
    
    if(minPos != std::string::npos)
    {
        return std::string("-1e+37");
    }
    
    if(maxPos != std::string::npos)
    {
        return std::string("1e+37");
    }
    
    return str;
}

// ****************************************************************************
//  Method: avtThresholdFilter::GetRangeList
//
//  Purpose: Parse the range string and create a list of ranges.
//
//  Arguments:
//      rangeStr           The range string to parse.
//
//  Returns: A vector of ranges
//
//  Programmer: Kevin Griffin
//  Creation:   Thu Mar 23 08:25:21 PDT 2017
//
//  Modifications:
// ****************************************************************************

stringVector
avtThresholdFilter::GetRangeList(const std::string rangeStr)
{
    stringVector rangeTokens;
    
    if(rangeStr.empty())
    {
        return rangeTokens;
    }
    
    // Create list of ranges
    std::size_t pos = rangeStr.find(",");
    std::size_t startPos = 0;
    
    while(pos != std::string::npos)
    {
        rangeTokens.push_back(rangeStr.substr(startPos, pos-startPos));
        startPos = pos+1;
        pos = rangeStr.find(",", startPos);
    }
    
    rangeTokens.push_back(rangeStr.substr(startPos));
    
    return rangeTokens;
}

// ****************************************************************************
//  Method: avtThresholdFilter::ThresholdOnRanges
//
//  Purpose: Threshold the dataset on a custom range
//
//  Arguments:
//      in_ds           The input VTK dataset.
//      threshold       The VTK threshold operator.
//      varName         The threshold variable name.
//      rangeStr        The range string consisting of a comma separated list
//                      of min-max and/or minmax (e.g. "1-10, 13, 20-6")
//      curZonePortion  Contains the value of the "Show zone if" option
//
//  Returns: The dataset after the threshold range has been applied
//
//  Programmer: Kevin Griffin
//  Creation:   Thu Mar 23 08:25:21 PDT 2017
//
//  Modifications:
//    Kevin Griffin, Fri Mar 24 11:36:58 PDT 2017
//    Changed the range symbol from '-' to ':' so it wouldn't conflict with
//    negative numbers. Also allow the use of 'min' and 'max' in ranges.
//
// ****************************************************************************

vtkDataSet *
avtThresholdFilter::ThresholdOnRanges(vtkDataSet *in_ds,
                                      vtkThreshold *threshold,
                                      const char *varName,
                                      const std::string rangeStr,
                                      const int curZonePortion)
{
    // Create list of ranges
    stringVector rangeTokens = GetRangeList(rangeStr);
    
    // Setup Threshold Data Array
    const char *keeperName = "_avt_threshold_keeper";
    const vtkDataObject::FieldAssociations FIELD_ASSOC = (in_ds->GetPointData()->GetArray(varName) != NULL) ? vtkDataObject::FIELD_ASSOCIATION_POINTS
                                                                                                            : vtkDataObject::FIELD_ASSOCIATION_CELLS;
    
    vtkDataArray *curVarArray = (FIELD_ASSOC == vtkDataObject::FIELD_ASSOCIATION_POINTS) ? in_ds->GetPointData()->GetArray(varName)
                                                                                         : in_ds->GetCellData()->GetArray(varName);
    
    vtkIntArray *keeper = vtkIntArray::New();
    keeper->SetName(keeperName);
    
    int keeperSize = curVarArray->GetSize();
    keeper->SetNumberOfTuples(keeperSize);
    
    // Initialize keeper
    for(int i=0; i<keeperSize; i++)
    {
        keeper->SetTuple1(i, 0);
    }
    
    // Determine what data we want to keep based on the ranges
    for(int i=0; i<rangeTokens.size(); i++)
    {
        if(rangeTokens[i].empty())
        {
            continue;
        }
        
        double threshLow = 0;
        double threshHigh = 0;
        
        size_t pos = rangeTokens[i].find(":");
        if(pos != std::string::npos)
        {
            std::string lowStr = rangeTokens[i].substr(0, pos).c_str();
            std::string highStr = rangeTokens[i].substr(pos+1).c_str();
        
            threshLow = atof(CheckForMinMax(lowStr).c_str());
            threshHigh = atof(CheckForMinMax(highStr).c_str());
        }
        else // single number
        {
            threshLow = atof(CheckForMinMax(rangeTokens[i]).c_str());
            threshHigh = threshLow;
        }
        
        debug5 << "threshLow: " << threshLow << " threshHigh: " << threshHigh << endl;
        
        for(int j=0; j<curVarArray->GetSize(); j++)
        {
            double value = curVarArray->GetTuple1(j);
            
            if((value >= threshLow) && (value <= threshHigh))
            {
                keeper->SetTuple1(j, 1);
            }
        }
    }
    
    // Add threshold array to the dataset
    if(FIELD_ASSOC == vtkDataObject::FIELD_ASSOCIATION_POINTS)
    {
        in_ds->GetPointData()->AddArray(keeper);
    }
    else
    {
        in_ds->GetCellData()->AddArray(keeper);
    }
    
    // Ready to threshold
    threshold->SetInputData(in_ds);
    threshold->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS, vtkDataSetAttributes::SCALARS);
    
    if (curZonePortion == (int)ThresholdOpAttributes::PartOfZone)
    {
        threshold->AllScalarsOff();
    }
    else if (curZonePortion == (int)ThresholdOpAttributes::EntireZone)
    {
        threshold->AllScalarsOn();
    }
    else
    {
        debug1 << "Invalid zone inclusion option encountered "
        << "in Threshold operator attributes." << endl;
        threshold->AllScalarsOff();
    }
    
    threshold->ThresholdBetween(1, 1);
    threshold->SetInputArrayToProcess(0, 0, 0, FIELD_ASSOC, keeperName);
    threshold->Update();
    
    in_ds = threshold->GetOutput();
    
    if(FIELD_ASSOC == vtkDataObject::FIELD_ASSOCIATION_POINTS)
    {
        in_ds->GetPointData()->RemoveArray(keeperName);
    }
    else
    {
        in_ds->GetCellData()->RemoveArray(keeperName);
    }
    
    keeper->Delete();
    
    return in_ds;
}

// ****************************************************************************
//  Method: avtThresholdFilter::ThresholdToPointMesh
//
//  Purpose:
//      Does a threshold on a nodal quantity and then creates a point mesh.
//
//  Programmer: Hank Childs
//  Creation:   September 13, 2005
//
//  Modifications:
//
//    Mark Blair, Tue Mar  7 13:25:00 PST 2006
//    Rewrote to support multi-variable thresholding.
//
//    Hank Childs, Sat Jan 27 12:53:20 PST 2007
//    Only add points that are actually incident to cells.
//
//    Hank Childs, Sun Jan 31 11:33:26 PST 2010
//    No longer have check for default variable.  It figures itself out when
//    iterating over all variables.  Also, the logic I'm removing with this
//    change was faulty ... if you had a Pseudocolor plot of an expression and
//    you were thresholding by that expression, then the default var is the
//    mesh and that was causing the exception to fire.
//
//    Brad Whitlock, Wed Mar 21 12:00:03 PDT 2012
//    Switch to GetTuple1 since this operation is not common. Change how
//    points are allocated.
//
//    Kevin Griffin, Thu Mar 23 08:25:21 PDT 2017
//    Added support for variables with multiple threshold ranges
//    (Feature #2646).
//
//    Kevin Griffin, Fri Mar 24 11:36:58 PDT 2017
//    Changed the range symbol from '-' to ':' so it wouldn't conflict with
//    negative numbers. Also allow the use of 'min' and 'max' in ranges.
//
// ****************************************************************************

vtkDataSet *
avtThresholdFilter::ThresholdToPointMesh(vtkDataSet *in_ds)
{
    
    const stringVector curVariables = atts.GetListedVarNames();
    int curVarCount = (int)curVariables.size();
    int curVarNum;
    vtkPointData *inPointData = in_ds->GetPointData();
    
    std::vector<vtkDataArray *> valueArrays;

    for (curVarNum = 0; curVarNum < curVarCount; curVarNum++)
    {
        vtkDataArray *dataArray = inPointData->GetArray(curVariables[curVarNum].c_str());
        
        if (dataArray == NULL)
        {
            EXCEPTION1(VisItException,
                "All threshold variables must be nodal quantities "
                "when point mesh output is requested.");
        }

        valueArrays.push_back(dataArray);
    }

    const intVector    curZonePortions = atts.GetZonePortions();
    const doubleVector curLowerBounds  = atts.GetLowerBounds();
    const doubleVector curUpperBounds  = atts.GetUpperBounds();
    const stringVector curBoundsRange = atts.GetBoundsRange();
    
    int inPointCount = in_ds->GetNumberOfPoints();
    int plotPointCount = 0;
    int inPointID;
    double doubleValue;

    for (inPointID = 0; inPointID < inPointCount; inPointID++)
    {
        for (curVarNum = 0; curVarNum < curVarCount; curVarNum++)
        {
            doubleValue = valueArrays[curVarNum]->GetTuple1(inPointID);
            
            if((atts.GetBoundsInputType() == ThresholdOpAttributes::Default) || IsSimpleRange(curBoundsRange[curVarNum]))
            {
                if (doubleValue < curLowerBounds[curVarNum]) break;
                if (doubleValue > curUpperBounds[curVarNum]) break;
            }
            else
            {
                bool breakLoop = true;
                stringVector rangeTokens = GetRangeList(curBoundsRange[curVarNum]);
                
                for(int i=0; i<rangeTokens.size(); i++)
                {
                    if(rangeTokens[i].empty())
                    {
                        continue;
                    }
                    
                    double threshLow = 0;
                    double threshHigh = 0;
                    
                    size_t pos = rangeTokens[i].find(":");
                    if(pos != std::string::npos)
                    {
                        threshLow = atof(CheckForMinMax(rangeTokens[i].substr(0, pos)).c_str());
                        threshHigh = atof(CheckForMinMax(rangeTokens[i].substr(pos+1)).c_str());
                    }
                    else // single number
                    {
                        threshLow = atof(CheckForMinMax(rangeTokens[i]).c_str());
                        threshHigh = threshLow;
                    }
                    
                    if((doubleValue >= threshLow) && (doubleValue <= threshHigh))
                    {
                        breakLoop = false;
                        break;
                    }
                }
                
                if(breakLoop)
                {
                    break;
                }
            }
        }

        if (curVarNum >= curVarCount) plotPointCount++;
    }

    if (plotPointCount == 0) return NULL;

    vtkUnstructuredGrid *outputMesh = vtkUnstructuredGrid::New();
    vtkPoints *outMeshPoints = vtkVisItUtility::NewPoints(in_ds);
    vtkPointData *outPointData = outputMesh->GetPointData();
    int outPointID = 0;
    vtkIdType vertexIDs[1];
    double pointXYZ[3];

    outMeshPoints->SetNumberOfPoints(plotPointCount);
    outputMesh->SetPoints(outMeshPoints);
    outputMesh->Allocate(plotPointCount*(curVarCount+1));
    outPointData->CopyAllocate(inPointData, plotPointCount*curVarCount);
    outMeshPoints->Delete();

    bool needToSeeIfPointIsIncidentToCell = true;
    if (in_ds->GetDataObjectType() == VTK_RECTILINEAR_GRID)
        needToSeeIfPointIsIncidentToCell = false;
    if (in_ds->GetDataObjectType() == VTK_STRUCTURED_GRID)
        needToSeeIfPointIsIncidentToCell = false;
    vtkIdList *idList = vtkIdList::New();
    for (inPointID = 0; inPointID < inPointCount; inPointID++)
    {
        for (curVarNum = 0; curVarNum < curVarCount; curVarNum++)
        {
            doubleValue = valueArrays[curVarNum]->GetTuple1(inPointID);
            
            if((atts.GetBoundsInputType() == ThresholdOpAttributes::Default) || IsSimpleRange(curBoundsRange[curVarNum]))
            {
                if (doubleValue < curLowerBounds[curVarNum]) break;
                if (doubleValue > curUpperBounds[curVarNum]) break;
            }
            else
            {
                bool breakLoop = true;
                stringVector rangeTokens = GetRangeList(curBoundsRange[curVarNum]);
                
                for(int i=0; i<rangeTokens.size(); i++)
                {
                    if(rangeTokens[i].empty())
                    {
                        continue;
                    }
                    
                    double threshLow = 0;
                    double threshHigh = 0;
                    
                    size_t pos = rangeTokens[i].find(":");
                    if(pos != std::string::npos)
                    {
                        threshLow = atof(CheckForMinMax(rangeTokens[i].substr(0, pos)).c_str());
                        threshHigh = atof(CheckForMinMax(rangeTokens[i].substr(pos+1)).c_str());
                    }
                    else // single number
                    {
                        threshLow = atof(CheckForMinMax(rangeTokens[i]).c_str());
                        threshHigh = threshLow;
                    }
                    
                    if((doubleValue >= threshLow) && (doubleValue <= threshHigh))
                    {
                        breakLoop = false;
                        break;
                    }
                }
                
                if(breakLoop)
                {
                    break;
                }
            }
        }

        bool shouldAdd = true;
        if (curVarNum < curVarCount)
            shouldAdd = false;

        if (shouldAdd && needToSeeIfPointIsIncidentToCell)
        {
            in_ds->GetPointCells(inPointID, idList);
            if (idList->GetNumberOfIds() <= 0)
                shouldAdd = false;
        }

        if (shouldAdd)
        {
            outPointData->CopyData(inPointData, inPointID, outPointID);
            in_ds->GetPoint(inPointID, pointXYZ);
            outMeshPoints->SetPoint(outPointID, pointXYZ);
            vertexIDs[0] = outPointID++;
            outputMesh->InsertNextCell(VTK_VERTEX, 1, vertexIDs);
        }
    }

    idList->Delete();
    return outputMesh;
}


// ****************************************************************************
//  Function: UpdateNeighborCells
//
//  Purpose: Updates cells that are neighbors to a point with a value.
//
//  Programmer: Hank Childs
//  Creation:   March 19, 2005
//
// ****************************************************************************

static void UpdateNeighborCells(int pt, const int *pt_dims,
           avtStructuredMeshChunker::ZoneDesignation d,
           std::vector<avtStructuredMeshChunker::ZoneDesignation> &designation)
{
    int I = pt % pt_dims[0];
    int J = (pt / pt_dims[0]) % pt_dims[1];
    int K = pt / (pt_dims[0]*pt_dims[1]);

    for (int i = 0 ; i < 8 ; i++)
    {
        int sI = I - (i & 1 ? 1 : 0);
        int sJ = J - (i & 2 ? 1 : 0);
        int sK = K - (i & 4 ? 1 : 0);
        if (sI < 0)
            continue;
        if (sI >= (pt_dims[0]-1))
            continue;
        if (sJ < 0)
            continue;
        if (sJ >= (pt_dims[1]-1))
            continue;
        if (sK < 0)
            continue;
        if (sK >= (pt_dims[2]-1))
            continue;

        int cell = sK*(pt_dims[0]-1)*(pt_dims[1]-1) + sJ*(pt_dims[0]-1) + sI;
        designation[cell] = d;
    }
}


// ****************************************************************************
//  Method: avtThresholdFilter::GetAssignments
//
//  Purpose: Gets the assignments for each zone.
//
//  Programmer: Hank Childs
//  Creation:   March 27, 2005
//
//  Modifications:
//    Mark Blair, Tue Mar  7 13:25:00 PST 2006
//    Reworked to support multi-variable thresholding.
//
//    Brad Whitlock, Wed Mar 21 11:52:40 PDT 2012
//    Support more than float.
//
// ****************************************************************************

void
avtThresholdFilter::GetAssignments(vtkDataSet *in_ds, const int *dims,
    std::vector<avtStructuredMeshChunker::ZoneDesignation> &d)
{
    const stringVector curVariables    = atts.GetListedVarNames();
    const intVector    curZonePortions = atts.GetZonePortions();
    const doubleVector curLowerBounds  = atts.GetLowerBounds();
    const doubleVector curUpperBounds  = atts.GetUpperBounds();

    vtkIdType zoneCount = in_ds->GetNumberOfCells();
    vtkIdType pointCount = in_ds->GetNumberOfPoints();
    size_t curVarNum;
    vtkIdType zoneNum, pointNum;

    std::vector<avtStructuredMeshChunker::ZoneDesignation> curVarZDs(zoneCount);

    for (zoneNum = 0; zoneNum < zoneCount; zoneNum++)
        d[zoneNum] = avtStructuredMeshChunker::RETAIN;

    for (curVarNum = 0; curVarNum < curVariables.size(); curVarNum++)
    {
        vtkDataArray *dataArray = NULL;
        const char *curVarName = curVariables[curVarNum].c_str();

        bool varIsPointData;
        if ((dataArray = in_ds->GetPointData()->GetArray(curVarName)) != NULL)
        {
            varIsPointData = true;
        }
        else if ((dataArray = in_ds->GetCellData()->GetArray(curVarName)) != NULL)
        {
            varIsPointData = false;
        }
        else
        {
            char errMsg[1024];
            sprintf (errMsg,
                "Data for variable \"%s\" is not currently available.",
                curVarName);
            EXCEPTION1(VisItException, errMsg);
        }

        double lowerBound = curLowerBounds[curVarNum];
        double upperBound = curUpperBounds[curVarNum];
        double varValue;
        const float *varValues = (const float *)dataArray->GetVoidPointer(0);

        if (varIsPointData)
        {
            if (curZonePortions[curVarNum] == ThresholdOpAttributes::PartOfZone)
            {
                for (zoneNum = 0; zoneNum < zoneCount; zoneNum++)
                    curVarZDs[zoneNum] = avtStructuredMeshChunker::DISCARD;

                if(dataArray->GetDataType() == VTK_FLOAT)
                {
                    for (pointNum = 0; pointNum < pointCount; pointNum++)
                    {
                        varValue = (double)varValues[pointNum];

                        if ((varValue >= lowerBound) && (varValue <= upperBound))
                        {
                            UpdateNeighborCells(pointNum, dims,
                            avtStructuredMeshChunker::RETAIN, curVarZDs);
                        }
                    }
                }
                else
                {
                    for (pointNum = 0; pointNum < pointCount; pointNum++)
                    {
                        varValue = dataArray->GetTuple1(pointNum);

                        if ((varValue >= lowerBound) && (varValue <= upperBound))
                        {
                            UpdateNeighborCells(pointNum, dims,
                            avtStructuredMeshChunker::RETAIN, curVarZDs);
                        }
                    }
                }
            }
            else
            {
                for (zoneNum = 0; zoneNum < zoneCount; zoneNum++)
                    curVarZDs[zoneNum] = avtStructuredMeshChunker::RETAIN;

                if(dataArray->GetDataType() == VTK_FLOAT)
                {
                    for (pointNum = 0; pointNum < pointCount; pointNum++)
                    {
                        varValue = (double)varValues[pointNum];

                        if ((varValue < lowerBound) || (varValue > upperBound))
                        {
                            UpdateNeighborCells(pointNum, dims,
                            avtStructuredMeshChunker::DISCARD, curVarZDs);
                        }
                    }
                }
                else
                {
                    for (pointNum = 0; pointNum < pointCount; pointNum++)
                    {
                        varValue = dataArray->GetTuple1(pointNum);

                        if ((varValue < lowerBound) || (varValue > upperBound))
                        {
                            UpdateNeighborCells(pointNum, dims,
                            avtStructuredMeshChunker::DISCARD, curVarZDs);
                        }
                    }
                }
            }
        }
        else
        {
            if(dataArray->GetDataType() == VTK_FLOAT)
            {
                for (zoneNum = 0; zoneNum < zoneCount; zoneNum++)
                {
                    varValue = (double)varValues[zoneNum];

                    if ((varValue >= lowerBound) && (varValue <= upperBound))
                        curVarZDs[zoneNum] = avtStructuredMeshChunker::RETAIN;
                    else
                        curVarZDs[zoneNum] = avtStructuredMeshChunker::DISCARD;
                }
            }
            else
            {
                for (zoneNum = 0; zoneNum < zoneCount; zoneNum++)
                {
                    varValue = dataArray->GetTuple1(zoneNum);

                    if ((varValue >= lowerBound) && (varValue <= upperBound))
                        curVarZDs[zoneNum] = avtStructuredMeshChunker::RETAIN;
                    else
                        curVarZDs[zoneNum] = avtStructuredMeshChunker::DISCARD;
                }
            }
        }

        for (zoneNum = 0; zoneNum < zoneCount; zoneNum++)
        {
            if (curVarZDs[zoneNum] == avtStructuredMeshChunker::DISCARD)
                d[zoneNum] = avtStructuredMeshChunker::DISCARD;
        }
    }
}


// ****************************************************************************
//  Method: avtThresholdFilter::UpdateDataObjectInfo
//
//  Purpose: Indicates the zones no longer correspond to the original problem.
//
//  Programmer: Hank Childs
//  Creation:   October 23, 2001
//
//  Modifications:
//
//    Hank Childs, Tue Sep 13 09:07:05 PDT 2005
//    Add support for "PointsOnly".
//
//    Mark Blair, Tue Mar  7 13:25:00 PST 2006
//    Rewrote to support multi-variable thresholding.
//
//    Kathleen Bonnell, Mon May  1 08:50:46 PDT 2006 
//    Set OrigElementsRequiredForPick. 
//
// ****************************************************************************

void
avtThresholdFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().InvalidateZones();

    if (atts.GetOutputMeshType() == ThresholdOpAttributes::PointMesh)
    {
        GetOutput()->GetInfo().GetAttributes().SetTopologicalDimension(0);
    }

    GetOutput()->GetInfo().GetAttributes().SetOrigElementsRequiredForPick(true);
}


// *****************************************************************************
//  Method: avtThresholdFilter::PreExecute
//
//  Purpose: Determine if there is a "default" variable to work with.
//
//  Programmer: Hank Childs
//  Creation:   August 28, 2002
//
//  Modifications:
//
//    Hank Childs, Thu Aug 29 08:12:22 PDT 2002
//    If the variable was not the 'default' variable, the test was sometimes
//    wrong.
//
//    Hank Childs, Sat Mar 19 10:29:37 PST 2005
//    Initialize chunkedStructuredMesh.
//
//    Hank Childs, Sun Mar 27 11:49:20 PST 2005
//    Moved chunkedStructuredMesh to base class.
//
//    Mark Blair, Tue Mar  7 13:25:00 PST 2006
//    No longer recognizes the "default" variable.
//
//    Mark Blair, Tue Aug  8 17:47:00 PDT 2006
//    Now accommodates an empty list of threshold variables; does pass-through.
//
//    Mark Blair, Tue Sep 18 17:06:28 PDT 2007
//    Removes any variable from the threshold list that's not a scalar or not in
//    the input.  Also determines whether or not default variable is a scalar.
//
//    Hank Childs, Thu Oct  9 09:24:26 PDT 2008
//    Issue a better error message for non-scalars.
//
//    Hank Childs, Thu Jan 21 19:50:58 PST 2010
//    Don't issue an error message if the default vector is a vector, but it
//    isn't actually used.
//
// *****************************************************************************

void
avtThresholdFilter::PreExecute(void)
{
    avtPluginStructuredChunkDataTreeIterator::PreExecute();

    int inputVarCount = GetInput()->GetInfo().GetAttributes().GetNumberOfVariables();
    int inputVarNum;
    std::string inputVarName;
    stringVector inputVarNames;
    
    for (inputVarNum = 0; inputVarNum < inputVarCount; inputVarNum++)
    {
        inputVarName =
            GetInput()->GetInfo().GetAttributes().GetVariableName(inputVarNum);
        
        if (GetInput()->GetInfo().GetAttributes().GetVariableType(inputVarName.c_str())
            == AVT_SCALAR_VAR)
            inputVarNames.push_back(inputVarName);
    }
    
    inputVarCount = (int)inputVarNames.size();
    bool allScalars = true;
    stringVector curListedVarNames = atts.GetListedVarNames();
    for (size_t i = 0 ; i < curListedVarNames.size() ; i++)
    {
         bool foundMatch = false;
         for (size_t j = 0 ; j < inputVarNames.size() ; j++)
             if (inputVarNames[j] == curListedVarNames[i])
                 foundMatch = true;
         if (! foundMatch)
             allScalars = false;
    }
    
    if (! allScalars)
    {
        atts.SetDefaultVarIsScalar(false);
        static bool issuedWarning = false;
        if (!issuedWarning)
        {
            avtCallback::IssueWarning("The threshold operator was asked to "
               "threshold by a non-scalar.  Only scalars are supported "
               "for thresholding.  You can overcome this by using "
               "expressions.  For example, if you would like to threshold "
               "using a vector, then you can define an expression for "
               "vector magnitude and threshold by that.  This warning "
               "will only be issued once per VisIt session.");
            issuedWarning = true;
        }
    }

    intVector    curZonePortions   = atts.GetZonePortions();
    doubleVector curLowerBounds    = atts.GetLowerBounds();
    doubleVector curUpperBounds    = atts.GetUpperBounds();

    std::string listVarName;
    size_t listVarNum;
    bool changedTheList;
    bool atLeast1Change = false;
    
    do
    {
        changedTheList = false;
    
        for (listVarNum = 0; listVarNum < curListedVarNames.size(); listVarNum++)
        {
            listVarName = curListedVarNames[listVarNum];
    
            for (inputVarNum = 0; inputVarNum < inputVarCount; inputVarNum++)
            {
                if (inputVarNames[inputVarNum] == listVarName) break;
            }
            
            if (inputVarNum >= inputVarCount)
            {
                curListedVarNames.erase(curListedVarNames.begin() + listVarNum);
                curZonePortions.erase  (curZonePortions.begin()   + listVarNum);
                curLowerBounds.erase   (curLowerBounds.begin()    + listVarNum);
                curUpperBounds.erase   (curUpperBounds.begin()    + listVarNum);
                
                changedTheList = true; atLeast1Change = true;
                
                break;
            }
        }
    }
    while (changedTheList);
    
    if (atLeast1Change)
    {
        atts.SetListedVarNames(curListedVarNames);
        atts.SetZonePortions(curZonePortions);
        atts.SetLowerBounds(curLowerBounds);
        atts.SetUpperBounds(curUpperBounds);
    }
}


// ****************************************************************************
//  Method: avtThresholdFilter::ModifyContract
//
//  Purpose: Restrict the data processed by looking at the data extents.
//
//  Programmer: Hank Childs
//  Creation:   May 14, 2003
//
//  Modifications:
//
//    Mark C. Miller, Mon Oct 18 13:02:37 PDT 2004
//    Added code to pass variable name in call to GetDataExtents
//
//    Hank Childs, Sat Mar 19 10:34:08 PST 2005
//    Initialize data members for structured mesh chunking.
//
//    Hank Childs, Sun Mar 27 11:49:20 PST 2005
//    Moved data members for structured mesh chunking to base class.
//
//    Mark Blair, Tue Mar  7 13:25:00 PST 2006
//    Rewrote to support multi-variable thresholding.
//
//    Kathleen Bonnell, Mon May  1 08:50:46 PDT 2006 
//    Turn on Node & Zone numbers when appropriate. 
//
//    Mark Blair, Tue Aug  8 17:47:00 PDT 2006
//    Now accommodates an empty list of threshold variables; does pass-through.
//
//    Kathleen Bonnell, Mon Aug 14 16:40:30 PDT 2006
//    API change for avtIntervalTree.
//
//    Mark Blair, Thu Sep 28 12:07:05 PDT 2006
//    Checks attributes for consistency.
//
//    Mark Blair, Tue Oct  3 13:19:11 PDT 2006
//    Replace "default" with true variable name if scalar, otherwise delete it
//    from list of threshold variables entirely.
//
//    Mark Blair, Wed Oct  4 17:45:48 PDT 2006
//    Makes better use of the "active variable".
//
//    Mark Blair, Tue Mar 13 19:51:29 PDT 2007
//    Now forces attribute consistency if attributes are inconsistent.
//
//    Mark Blair, Tue Apr 17 16:24:42 PDT 2007
//    Rewritten to support new Threshold GUI; no more "shown variable".
//
//    Mark Blair, Tue Sep 18 17:06:28 PDT 2007
//    API change: New method name "SwitchDefaultVariableNameToTrueName".
//
//    Sean Ahern (Markus Glatter), Tue Sep 25 10:21:40 EDT 2007
//    Added avtDataRangeSelections for each variable.
//
// ****************************************************************************

avtContract_p
avtThresholdFilter::ModifyContract(avtContract_p in_spec)
{
    atts.SupplyMissingDefaultsIfAppropriate();
    
    if (!atts.AttributesAreConsistent()) atts.ForceAttributeConsistency();

    const char *pipelineVar = in_spec->GetDataRequest()->GetVariable();
    const char *activeVar = activeVarName.c_str();
    
    if (activeVarName == std::string("<unused>"))
        activeVar = pipelineVar;
        
    atts.SetDefaultVarName(std::string(pipelineVar));

    atts.SwitchDefaultVariableNameToTrueName();
    
    if (atts.GetListedVarNames().size() == 0) return in_spec;

    avtContract_p outSpec = new avtContract(in_spec);

    const char *curListedVar;
    const std::vector<CharStrRef> curSecondaryVars =
        outSpec->GetDataRequest()->GetSecondaryVariables();
    const stringVector curListedVars = atts.GetListedVarNames();
    size_t listedVarNum;
    size_t secVarNum;

    for (listedVarNum = 0; listedVarNum < curListedVars.size(); listedVarNum++)
    {
        curListedVar = curListedVars[listedVarNum].c_str();
        
        if ((strcmp(curListedVar, pipelineVar) != 0) &&
            (strcmp(curListedVar, activeVar  ) != 0))
        {
            for (secVarNum = 0; secVarNum < curSecondaryVars.size(); secVarNum++)
            {
                if (strcmp(*curSecondaryVars[secVarNum],curListedVar) == 0)
                {
                    break;
                }
            }

            if (secVarNum >= curSecondaryVars.size())
            {
                outSpec->GetDataRequest()->AddSecondaryVariable(curListedVar);
            }
        }
    }

    bool atLeastOneTree = false;
    avtIntervalTree *it;
    size_t varDomNum;
    size_t curDomNum;
    int    curDomain;
    intVector varDomains;
    intVector curDomains;
    intVector outDomains;

    doubleVector curLowerBounds = atts.GetLowerBounds();
    doubleVector curUpperBounds = atts.GetUpperBounds();
    double lowerBound, upperBound;

    for (listedVarNum = 0; listedVarNum < curListedVars.size(); listedVarNum++)
    {
        curListedVar = curListedVars[listedVarNum].c_str();
        
        if ((it = GetMetaData()->GetDataExtents(curListedVar)) != NULL)
        {
            lowerBound = curLowerBounds[listedVarNum];
            upperBound = curUpperBounds[listedVarNum];

            if (atLeastOneTree)
            {
                it->GetElementsListFromRange(&lowerBound, &upperBound, varDomains);

                for (curDomNum = 0; curDomNum < curDomains.size(); curDomNum++)
                {
                    if ((curDomain = curDomains[curDomNum]) != -1)
                    {
                        for (varDomNum=0; varDomNum<varDomains.size(); varDomNum++)
                        {
                            if (varDomains[varDomNum] == curDomain) break;
                        }

                        if (varDomNum >= varDomains.size())
                        {
                            curDomains[curDomNum] = -1;
                        }
                    }
                }
            }
            else
            {
                it->GetElementsListFromRange(&lowerBound, &upperBound, curDomains);
                atLeastOneTree = true;
            }
        }
    }

    if (atLeastOneTree)
    {
        for (curDomNum = 0; curDomNum < curDomains.size(); curDomNum++)
        {
            if (curDomains[curDomNum] != -1)
            {
                outDomains.push_back(curDomains[curDomNum]);
            }
        }

        outSpec->GetDataRequest()->GetRestriction()->RestrictDomains(outDomains);
    }

    if (outSpec->GetDataRequest()->MayRequireZones() ||
        outSpec->GetDataRequest()->MayRequireNodes())
    {
        // Turn on both Nodes and Zones, to prevent another re-execution if
        // user switches between zone and node pick.
        outSpec->GetDataRequest()->TurnZoneNumbersOn();
        outSpec->GetDataRequest()->TurnNodeNumbersOn();
    }

    // Add avtDataRangeSelection entries to data specs
    // (Enables contract-based filtering if supported by db)
    for (listedVarNum = 0; listedVarNum < curListedVars.size(); listedVarNum++)
    {
        lowerBound = curLowerBounds[listedVarNum];
        upperBound = curUpperBounds[listedVarNum];

        std::string curVar(curListedVars[listedVarNum]);

        selIDs[curVar] = outSpec->GetDataRequest()->AddDataSelection(
            new avtDataRangeSelection(curVar, lowerBound, upperBound)
            );

        debug1 << "Added variable " << curVar
            << " as Data Selection with range ("
            << lowerBound << ", " << upperBound << ")" << endl;
    }

    return outSpec;
}


// ****************************************************************************
//  Method: avtThresholdFilter::CreateNamedSelection
//
//  Purpose:
//      Creates a named selection.  This will only do something useful if
//      the upstream database can return an ID list when given a selection.
//
//  Programmer: Hank Childs
//  Creation:   February 23, 2009
//
//  Modifications:
//  
//    Gunther H. Weber, Mon Apr 27 20:41:18 PDT 2009
//    Fix a crash due to deleting a reference we don't own. (Analogous to
//    a fix Hank did on Mon Apr 6 17:13:58 PDT 2009 in 
//    avtParallelCoordinatesFilter::CreateDBAcceleratedNamedSelection 
//
//    Brad Whitlock, Fri Oct 28 11:00:26 PDT 2011
//    Change named selection API.
//
//    Brad Whitlock, Thu Mar 15 14:19:04 PDT 2012
//    Set the id variable for the floating point named selection.
//
// ****************************************************************************

avtNamedSelection *
avtThresholdFilter::CreateNamedSelection(avtContract_p c, const std::string &s)
{
    if (! GetInput()->GetInfo().GetValidity().GetZonesPreserved())
    {
        // Zones have been removed upstream, so the direct-to-database query
        // will be invalid.  Give up.
        return NULL;
    }

    std::vector<avtDataSelection *> drs;

    const stringVector curListedVars = atts.GetListedVarNames();
    doubleVector curLowerBounds = atts.GetLowerBounds();
    doubleVector curUpperBounds = atts.GetUpperBounds();
    double lowerBound, upperBound;

    for (size_t listedVarNum = 0; listedVarNum < curListedVars.size(); listedVarNum++)
    {
        lowerBound = curLowerBounds[listedVarNum];
        upperBound = curUpperBounds[listedVarNum];

        std::string var = curListedVars[listedVarNum];
        avtDataRangeSelection *sel = new avtDataRangeSelection(var, lowerBound, 
                                                               upperBound);
        drs.push_back(sel);
    }

    avtIdentifierSelection *ids = GetMetaData()->GetIdentifiers(drs);
    avtNamedSelection *rv = NULL;
    if (ids != NULL)
    {
        avtFloatingPointIdNamedSelection *fpns = new avtFloatingPointIdNamedSelection(s);
        fpns->SetIdentifiers(ids->GetIdentifiers());
        fpns->SetIdVariable(ids->GetIdVariable());
        rv = fpns;
    }

    // Don't delete ids, since it is being cached at the DB level and we don't
    // own this reference.
    // delete ids;

    for (size_t i = 0 ; i < drs.size() ; i++)
        delete drs[i];

    return rv;
}


