/*****************************************************************************
*
* 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.
*
*****************************************************************************/

// ************************************************************************* //
//  File: avtSimilarityTransformFilter.C
// ************************************************************************* //

#include <avtSimilarityTransformFilter.h>

#include <avtExtents.h>
#include <vtkDataArray.h>
#include <vtkPointData.h>
#include <vtkRectilinearGrid.h>

#include <BadVectorException.h>


// ****************************************************************************
//  Method: avtSimilarityTransformFilter constructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Mon Sep 24 13:29:16 PST 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Initialize invM.
//
// ****************************************************************************

avtSimilarityTransformFilter::avtSimilarityTransformFilter()
{
    M = NULL;
    invM = NULL;
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter destructor
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Mon Sep 24 13:29:16 PST 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Delete invM.
//
// ****************************************************************************

avtSimilarityTransformFilter::~avtSimilarityTransformFilter()
{
    if (M)
        M->Delete();
    M = NULL;
    if (invM)
        invM->Delete();
    invM = NULL;
}


// ****************************************************************************
//  Method:      avtSimilarityTransformFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: meredith -- generated by xml2info
//  Creation:   Mon Sep 24 13:29:16 PST 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Delete invM.
//
//    Kathleen Bonnell, Wed May 21 11:38:23 PDT 2003   
//    Check for bad axis of rotation. 
//
//    Dave Pugmire, Fri May 14 08:04:43 EDT 2010
//    Flag for vector transformations.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const SimilarityTransformAttributes*)a;
    if (M)
        M->Delete();
    M = NULL;
    if (invM)
        invM->Delete();
    invM = NULL;

    SetVectorTransform(atts.GetTransformVectors());

    const double *axis = atts.GetRotateAxis();
    if (axis[0] == 0. && axis[1] == 0. && axis[2] == 0.)
    {
        EXCEPTION1(BadVectorException, "Rotation Axis");
        return;
    }
}


// ****************************************************************************
//  Method:  Jeremy Meredith, Mon Sep 24 14:24:14 PDT 2001
//
//  Purpose:
//    Setup the vtk matrix from the transform attributs
//
//  Programmer:  Jeremy Meredith
//  Creation:    September 24, 2001
//
//  Modifications:
//    Kathleen Bonnell, Thu Apr 10 11:07:48 PDT 2003
//    Compute the inverse matrix. 
//
// ****************************************************************************

void
avtSimilarityTransformFilter::SetupMatrix()
{
    if (M)
        return;

    M = vtkMatrix4x4::New();
    M->Identity();
    invM = vtkMatrix4x4::New();
    invM->Identity();

    //
    //  A place to store the inverse transforms for each type.
    //
    vtkMatrix4x4 *IR = vtkMatrix4x4::New(); // inverse Rotate
    IR->Identity();
    vtkMatrix4x4 *IS = vtkMatrix4x4::New(); // inverse Scale
    IS->Identity();
    vtkMatrix4x4 *IT = vtkMatrix4x4::New(); // inverse Transform
    IT->Identity();

    //
    // Do the rotation
    //
    if (atts.GetDoRotate())
    {
        double oX = atts.GetRotateOrigin()[0];
        double oY = atts.GetRotateOrigin()[1];
        double oZ = atts.GetRotateOrigin()[2];

        double X = atts.GetRotateAxis()[0];
        double Y = atts.GetRotateAxis()[1];
        double Z = atts.GetRotateAxis()[2];
        double p = atts.GetRotateAmount();

        // Convert to radians
        if (atts.GetRotateType() == atts.Deg)
            p *= 3.1415926535898/180.;

        // Normalize the axis vector
        double len = sqrt(X*X + Y*Y + Z*Z);
        if (len)
        {
            X /= len;
            Y /= len;
            Z /= len;
        }

        // Compose the quaternion
        double sin_p = sin(p / 2.);
        double cos_p = cos(p / 2.);

        double q[4] = { X*sin_p, Y*sin_p, Z*sin_p, cos_p };

        vtkMatrix4x4 *T1  = vtkMatrix4x4::New();
        vtkMatrix4x4 *R   = vtkMatrix4x4::New();
        vtkMatrix4x4 *T2  = vtkMatrix4x4::New();
        vtkMatrix4x4 *tmp = vtkMatrix4x4::New();

        // Pre-translate
        T1->Identity();
        (*T1)[0][3] = -oX;
        (*T1)[1][3] = -oY;
        (*T1)[2][3] = -oZ;
        
        // Rotate
        R->Identity();
        (*R)[0][0] = 1.0 - 2.0 * (q[1] * q[1] + q[2] * q[2]);
        (*R)[0][1] = 2.0 * (q[0] * q[1] - q[2] * q[3]);
        (*R)[0][2] = 2.0 * (q[2] * q[0] + q[1] * q[3]);

        (*R)[1][0] = 2.0 * (q[0] * q[1] + q[2] * q[3]);
        (*R)[1][1] = 1.0 - 2.0 * (q[2] * q[2] + q[0] * q[0]);
        (*R)[1][2] = 2.0 * (q[1] * q[2] - q[0] * q[3]);

        (*R)[2][0] = 2.0 * (q[2] * q[0] - q[1] * q[3]);
        (*R)[2][1] = 2.0 * (q[1] * q[2] + q[0] * q[3]);
        (*R)[2][2] = 1.0 - 2.0 * (q[1] * q[1] + q[0] * q[0]);

        // Post-translate
        T2->Identity();
        (*T2)[0][3] = oX;
        (*T2)[1][3] = oY;
        (*T2)[2][3] = oZ;

        // apply it
        vtkMatrix4x4::Multiply4x4(T1,  M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(R,   M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  M,  tmp);   M->DeepCopy(tmp);

        // Create the inverse rotation.
        vtkMatrix4x4 *R_t= vtkMatrix4x4::New();
        vtkMatrix4x4::Transpose(R, R_t);

        vtkMatrix4x4::Multiply4x4(T1,  IR,  tmp);   IR->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(R_t, IR,  tmp);   IR->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  IR,  tmp);   IR->DeepCopy(tmp);

        R_t->Delete();
        T1 ->Delete();
        R  ->Delete();
        T2 ->Delete();
        tmp->Delete();
    }
    //
    // Do the scale
    //
    if (atts.GetDoScale())
    {
        double oX = atts.GetScaleOrigin()[0];
        double oY = atts.GetScaleOrigin()[1];
        double oZ = atts.GetScaleOrigin()[2];

        double X = atts.GetScaleX();
        double Y = atts.GetScaleY();
        double Z = atts.GetScaleZ();

        vtkMatrix4x4 *T1  = vtkMatrix4x4::New();
        vtkMatrix4x4 *S   = vtkMatrix4x4::New();
        vtkMatrix4x4 *T2  = vtkMatrix4x4::New();
        vtkMatrix4x4 *tmp = vtkMatrix4x4::New();

        // Pre-translate
        T1->Identity();
        (*T1)[0][3] = -oX;
        (*T1)[1][3] = -oY;
        (*T1)[2][3] = -oZ;

        // Scale
        S->Identity();
        (*S)[0][0] = X;
        (*S)[1][1] = Y;
        (*S)[2][2] = Z;

        // Post-translate
        T2->Identity();
        (*T2)[0][3] = oX;
        (*T2)[1][3] = oY;
        (*T2)[2][3] = oZ;

        // apply it
        vtkMatrix4x4::Multiply4x4(T1,  M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(S,   M,  tmp);   M->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  M,  tmp);   M->DeepCopy(tmp);

        // create the inverse scale
        vtkMatrix4x4 *S_i= vtkMatrix4x4::New();
        S_i->Identity();
        (*S_i)[0][0] = (X != 0. ? 1./X : 0.);
        (*S_i)[1][1] = (Y != 0. ? 1./Y : 0.);
        (*S_i)[2][2] = (Z != 0. ? 1./Z : 0.);
        vtkMatrix4x4::Multiply4x4(T1,  IS,  tmp);   IS->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(S_i, IS,  tmp);   IS->DeepCopy(tmp);
        vtkMatrix4x4::Multiply4x4(T2,  IS,  tmp);   IS->DeepCopy(tmp);

        S_i->Delete(); 
        T1 ->Delete();
        S  ->Delete();
        T2 ->Delete();
        tmp->Delete();
    }
    //
    // Do the translation
    //
    if (atts.GetDoTranslate())
    {
        vtkMatrix4x4 *T = vtkMatrix4x4::New();
        vtkMatrix4x4 *tmp = vtkMatrix4x4::New();

        // Translate
        (*T)[0][3] = atts.GetTranslateX();
        (*T)[1][3] = atts.GetTranslateY();
        (*T)[2][3] = atts.GetTranslateZ();

        // apply it
        vtkMatrix4x4::Multiply4x4(T,  M,   tmp);   M->DeepCopy(tmp);

        // Inverse translation
        (*IT)[0][3] = -atts.GetTranslateX();
        (*IT)[1][3] = -atts.GetTranslateY();
        (*IT)[2][3] = -atts.GetTranslateZ();

        T  ->Delete();
        tmp->Delete();
    }

    //
    // Compose the separate inverse matrices into one.
    //  IT * IS * IR
    // 
    vtkMatrix4x4 *tmp = vtkMatrix4x4::New();
    vtkMatrix4x4::Multiply4x4(IT, invM, tmp); invM->DeepCopy(tmp);
    vtkMatrix4x4::Multiply4x4(IS, invM, tmp); invM->DeepCopy(tmp);
    vtkMatrix4x4::Multiply4x4(IR, invM, tmp); invM->DeepCopy(tmp);

    tmp->Delete();

    IR  ->Delete();
    IS  ->Delete();
    IT  ->Delete();
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::PerformRestriciton
//
//  Purpose:
//    Turn on Zone numbers flag if needed, so that original cell array
//    will be propagated throught the pipeline.
//
//  Programmer: Kathleen Bonnell
//  Creation:   November 28, 2001
//
//  Modifications:
//    Kathleen Bonnell, Wed Jun 19 12:28:10 PDT 2002
//    Don't turn off Zone numbers if they have been turned on elsewhere in
//    the pipeline.
//
//    Kathleen Bonnell, Wed Jun 19 13:42:37 PDT 2002
//    Completely removed the code turning off zone numbers.  Why set a flag
//    to false if it is already false?  False is the default setting.
//
//    Kathleen Bonnell, Wed Jun  2 09:08:35 PDT 2004 
//    Added code to turn on node numbers when required. 
//
// ****************************************************************************

avtContract_p
avtSimilarityTransformFilter::ModifyContract(avtContract_p spec)
{
    avtContract_p rv = new avtContract(spec);
    if (rv->GetDataRequest()->MayRequireZones())
    {
        rv->GetDataRequest()->TurnZoneNumbersOn();
    }
    if (rv->GetDataRequest()->MayRequireNodes())
    {
        rv->GetDataRequest()->TurnNodeNumbersOn();
    }
    return rv;
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::UpdateDataObjectInfo
//
//  Purpose:
//      If a 2D plot is being revolved into 3D space, then indicate that it is
//      a 3D plot.
//
//  Programmer: Hank Childs
//  Creation:   March 7, 2003
//
//  Modifications:
//  
//    Hank Childs, Thu Jun 17 10:44:20 PDT 2004
//    Set the topological dimension to be higher so that we know we need
//    to apply the facelist filter.  Also transform the extents.
//    
//    Hank Childs, Thu Jul 28 09:07:31 PDT 2005
//    Fix memory leak.
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::UpdateDataObjectInfo(void)
{
    //
    // The base class does some good work about setting extents.  Use that.
    //
    avtTransform::UpdateDataObjectInfo();

    if (atts.GetDoRotate())
    {
        double X = atts.GetRotateAxis()[0];
        double Y = atts.GetRotateAxis()[1];
        if (X != 0. || Y != 0.)
        {
            avtDataAttributes &inAtts = GetInput()->GetInfo().GetAttributes();
            avtDataAttributes &outAtts= GetOutput()->GetInfo().GetAttributes();
            if (inAtts.GetSpatialDimension() < 3)
            {
                //
                // Tell the new plot that it has spatial dimension 3.  Also
                // lie and say that its topological dimension is also 3.  This
                // will cause the facelist filter to execute.
                //
                outAtts.SetSpatialDimension(3);
                if (inAtts.GetTopologicalDimension() == 2)
                    outAtts.SetTopologicalDimension(3);

                double extents[6];
                vtkMatrix4x4 *t = GetTransform();

                //
                // Transform the extents into 3D.
                //
                if (inAtts.GetOriginalSpatialExtents()->HasExtents())
                {
                    avtExtents *newOriginalSpatial = new avtExtents(3);
                    inAtts.GetOriginalSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newOriginalSpatial->Set(extents);
                    *(outAtts.GetOriginalSpatialExtents()) = *newOriginalSpatial;
                    outAtts.GetOriginalSpatialExtents()->Transform(t);
                    delete newOriginalSpatial;
                }

                if (inAtts.GetThisProcsOriginalSpatialExtents()->HasExtents())
                {
                    avtExtents *newThisProcsOriginalSpatial = new avtExtents(3);
                    inAtts.GetThisProcsOriginalSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newThisProcsOriginalSpatial->Set(extents);
                    *(outAtts.GetThisProcsOriginalSpatialExtents()) = 
                                                     *newThisProcsOriginalSpatial;
                    outAtts.GetThisProcsOriginalSpatialExtents()->Transform(t);
                    delete newThisProcsOriginalSpatial;
                }

                if (inAtts.GetDesiredSpatialExtents()->HasExtents())
                {
                    avtExtents *newDesiredSpatial = new avtExtents(3);
                    inAtts.GetDesiredSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newDesiredSpatial->Set(extents);
                    *(outAtts.GetDesiredSpatialExtents()) = *newDesiredSpatial;
                    outAtts.GetDesiredSpatialExtents()->Transform(t);
                    delete newDesiredSpatial;
                }

                if (inAtts.GetActualSpatialExtents()->HasExtents())
                {
                    avtExtents *newActualSpatial = new avtExtents(3);
                    inAtts.GetActualSpatialExtents()->CopyTo(extents);
                    extents[4] = extents[5] = 0.;
                    newActualSpatial->Set(extents);
                    *(outAtts.GetActualSpatialExtents()) = *newActualSpatial;
                    outAtts.GetActualSpatialExtents()->Transform(t);
                    delete newActualSpatial;
                }
            }
        }
    }
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::PostExecute
//
//  Purpose:
//      This is called to set the inverse transformation matrix in the output.  
//
//  Programmer: Kathleen Bonnell 
//  Creation:   April 10, 2003 
//
//  Modifications:
//    Kathleen Bonnell, Wed Jun  2 09:08:35 PDT 2004
//    Set both InvTransform and Transform in the output.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::PostExecute()
{
    GetOutput()->GetInfo().GetAttributes().SetInvTransform((*invM)[0]);
    GetOutput()->GetInfo().GetAttributes().SetTransform((*M)[0]);
}


// ****************************************************************************
//  Method: avtSimilarityTransformFilter::TransformData
//
//  Purpose:
//    Performs Translation and Scaling on the point-data scalars of a 
//    1D Rectilinear grid representing a curve.
//
//  Programmer: Kathleen Bonnell 
//  Creation:   March 28, 2008 
//
//  Modifications:
//    Kathleen Biagas, Wed Aug 22 16:00:03 MST 2012
//    Preserve data type of the scalars.
//
// ****************************************************************************

void
avtSimilarityTransformFilter::TransformData(vtkRectilinearGrid *rgrid)
{
    double yScaleVal = atts.GetScaleY();
    double yTransVal = atts.GetTranslateY();
    bool doScale = atts.GetDoScale() && yScaleVal != 1.;
    bool doTrans = atts.GetDoTranslate() && yTransVal != 0.;

    if (doScale || doTrans)
    {
        vtkDataArray *scalars = rgrid->GetPointData()->GetScalars();
        vtkDataArray *transScalars = scalars->NewInstance();
        transScalars->SetNumberOfTuples(scalars->GetNumberOfTuples());
        transScalars->SetName(scalars->GetName());
        transScalars->DeepCopy(scalars);
        for (vtkIdType i = 0; i < scalars->GetNumberOfTuples(); i++)
        {
            double val = transScalars->GetTuple1(i);
            if (doScale)
                val *= yScaleVal;
            if (doTrans)
                val += yTransVal;
            transScalars->SetTuple1(i, val);
        }
        rgrid->GetPointData()->SetScalars(transScalars);
        transScalars->Delete();
    }
}
