/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2018-     Statoil ASA
//
//  ResInsight is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
//  WARRANTY; without even the implied warranty of MERCHANTABILITY or
//  FITNESS FOR A PARTICULAR PURPOSE.
//
//  See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "Riv3dWellLogCurveGeomertyGenerator.h"

#include "RigCurveDataTools.h"
#include "RigWellPath.h"

#include "cafDisplayCoordTransform.h"
#include "cvfPrimitiveSetIndexedUInt.h"

#include <cmath>

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::ref<cvf::DrawableGeo>
    Riv3dWellLogCurveGeometryGenerator::createCurveLine(const caf::DisplayCoordTransform* displayCoordTransform,
                                                        const Rim3dWellLogCurve*          rim3dWellLogCurve) const
{
    std::vector<cvf::Vec3f> vertices;
    std::vector<cvf::uint>  indices;

    createCurveVerticesAndIndices(rim3dWellLogCurve, displayCoordTransform, &vertices, &indices);

    if (vertices.empty() || indices.empty())
    {
        return nullptr;
    }

    cvf::ref<cvf::PrimitiveSetIndexedUInt> indexedUInt = new cvf::PrimitiveSetIndexedUInt(cvf::PrimitiveType::PT_LINES);
    cvf::ref<cvf::UIntArray>               indexArray  = new cvf::UIntArray(indices);

    cvf::ref<cvf::DrawableGeo> drawable = new cvf::DrawableGeo();

    indexedUInt->setIndices(indexArray.p());
    drawable->addPrimitiveSet(indexedUInt.p());

    cvf::ref<cvf::Vec3fArray> vertexArray = new cvf::Vec3fArray(vertices);
    drawable->setVertexArray(vertexArray.p());

    return drawable;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::ref<cvf::DrawableGeo> Riv3dWellLogCurveGeometryGenerator::createGrid(const caf::DisplayCoordTransform* displayCoordTransform,
                                                                          const Rim3dWellLogCurve::DrawPlane drawPlane,
                                                                          double gridIntervalSize) const
{
    CVF_ASSERT(gridIntervalSize > 0);

    if (m_wellPathGeometry.isNull()) return nullptr;

    std::vector<cvf::Vec3d> wellPathPoints = m_wellPathGeometry->m_wellPathPoints;
    if (wellPathPoints.empty()) return nullptr;
    
    const cvf::Vec3d globalDirection = (wellPathPoints.back() - wellPathPoints.front()).getNormalized();
    const cvf::Vec3d up(0, 0, 1);

    std::vector<cvf::Vec3d> pointNormals;
    std::vector<cvf::Vec3d> gridPoints;

    double firstMd = m_wellPathGeometry->m_measuredDepths.front();
    double lastMd  = m_wellPathGeometry->m_measuredDepths.back();
    if (m_wellPathGeometry->m_measuredDepths.empty()) return nullptr;

    double md = firstMd;
    while (md <= lastMd)
    {
        cvf::Vec3d point = m_wellPathGeometry->interpolatedPointAlongWellPath(md);
        gridPoints.push_back(point);
        md += gridIntervalSize;
    }

    pointNormals = calculatePointNormals(drawPlane, gridPoints);
    if (pointNormals.size() != gridPoints.size()) return nullptr;

    std::vector<cvf::Vec3f> vertices;
    vertices.reserve(gridPoints.size() * 2);

    std::vector<cvf::uint> indices;
    indices.reserve(gridPoints.size() * 2);

    cvf::uint counter = 0;

    for (size_t i = 0; i < pointNormals.size(); i++)
    {
        vertices.push_back(cvf::Vec3f(displayCoordTransform->transformToDisplayCoord(gridPoints[i])));
        vertices.push_back(cvf::Vec3f(displayCoordTransform->transformToDisplayCoord(gridPoints[i] + pointNormals[i] * 100)));

        indices.push_back(counter++);
        indices.push_back(counter++);
    }

    cvf::ref<cvf::PrimitiveSetIndexedUInt> indexedUInt = new cvf::PrimitiveSetIndexedUInt(cvf::PrimitiveType::PT_LINES);
    cvf::ref<cvf::UIntArray>               indexArray  = new cvf::UIntArray(indices);

    cvf::ref<cvf::DrawableGeo> drawable = new cvf::DrawableGeo();

    indexedUInt->setIndices(indexArray.p());
    drawable->addPrimitiveSet(indexedUInt.p());

    cvf::ref<cvf::Vec3fArray> vertexArray = new cvf::Vec3fArray(vertices);
    drawable->setVertexArray(vertexArray.p());

    return drawable;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void Riv3dWellLogCurveGeometryGenerator::createCurveVerticesAndIndices(const Rim3dWellLogCurve*          rim3dWellLogCurve,
                                                                       const caf::DisplayCoordTransform* displayCoordTransform,
                                                                       std::vector<cvf::Vec3f>*          vertices,
                                                                       std::vector<cvf::uint>*           indices) const
{
    if (m_wellPathGeometry.isNull()) return;

    std::vector<double> resultValues;
    std::vector<double> mds;
    rim3dWellLogCurve->curveValuesAndMds(&resultValues, &mds);

    if (resultValues.empty()) return;
    CVF_ASSERT(resultValues.size() == mds.size());

    cvf::Vec3d globalDirection =
        (m_wellPathGeometry->m_wellPathPoints.back() - m_wellPathGeometry->m_wellPathPoints.front()).getNormalized();

    std::vector<cvf::Vec3d> interpolatedWellPathPoints;
    interpolatedWellPathPoints.reserve(mds.size());

    for (double md : mds)
    {
        interpolatedWellPathPoints.push_back(m_wellPathGeometry->interpolatedPointAlongWellPath(md));
    }

    std::vector<cvf::Vec3d> pointNormals = calculatePointNormals(rim3dWellLogCurve->drawPlane(), interpolatedWellPathPoints);
    if (interpolatedWellPathPoints.size() != pointNormals.size()) return;

    double maxResult = -HUGE_VAL;
    double minResult = HUGE_VAL;

    for (double result : resultValues)
    {
        if (!RigCurveDataTools::isValidValue(result, false)) continue;

        maxResult = std::max(result, maxResult);
        minResult = std::min(result, minResult);
    }

    vertices->resize(interpolatedWellPathPoints.size());

    double range  = maxResult - minResult;
    double factor = 60.0 / range;
    double offset = 30.0;

    if (minResult < 0)
    {
        offset += cvf::Math::abs(minResult * factor);
    }

    for (size_t i = 0; i < pointNormals.size(); i++)
    {
        cvf::Vec3d result(0, 0, 0);

        if (RigCurveDataTools::isValidValue(resultValues[i], false))
        {
            result = resultValues[i] * factor * pointNormals[i];
        }

        (*vertices)[i] = cvf::Vec3f(
            displayCoordTransform->transformToDisplayCoord(interpolatedWellPathPoints[i] + pointNormals[i] * offset + result));
    }

    std::vector<std::pair<size_t, size_t>> valuesIntervals =
        RigCurveDataTools::calculateIntervalsOfValidValues(resultValues, false);

    for (const std::pair<size_t, size_t>& interval : valuesIntervals)
    {
        for (size_t i = interval.first; i < interval.second; i++)
        {
            indices->push_back(cvf::uint(i));
            indices->push_back(cvf::uint(i + 1));
        }
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<cvf::Vec3d> Riv3dWellLogCurveGeometryGenerator::calculatePointNormals(Rim3dWellLogCurve::DrawPlane   drawPlane,
                                                                                  const std::vector<cvf::Vec3d>& points) const
{
    std::vector<cvf::Vec3d> pointNormals;
    
    if (m_wellPathGeometry.isNull()) return pointNormals;
    pointNormals.reserve(points.size());

    const cvf::Vec3d globalDirection = (points.back() - points.front()).getNormalized();
    const cvf::Vec3d up(0, 0, 1);

    for (const cvf::Vec3d point : points)
    {
        cvf::Vec3d p1 = cvf::Vec3d::UNDEFINED;
        cvf::Vec3d p2 = cvf::Vec3d::UNDEFINED;
        m_wellPathGeometry->twoClosestPoints(point, &p1, &p2);
        if (p1.isUndefined() || p2.isUndefined()) continue;

        cvf::Vec3d vecAlongPath = (p2 - p1).getNormalized();

        double dotProduct = up * vecAlongPath;

        cvf::Vec3d Ex;

        if (cvf::Math::abs(dotProduct) > 0.7071)
        {
            Ex = globalDirection;
        }
        else
        {
            Ex = vecAlongPath;
        }

        cvf::Vec3d Ey = (up ^ Ex).getNormalized();
        cvf::Vec3d Ez = (Ex ^ Ey).getNormalized();

        cvf::Vec3d normal;

        switch (drawPlane)
        {
        case Rim3dWellLogCurve::HORIZONTAL_LEFT:
            normal = -Ey;
            break;
        case Rim3dWellLogCurve::HORIZONTAL_RIGHT:
            normal = Ey;
            break;
        case Rim3dWellLogCurve::VERTICAL_ABOVE:
            normal = Ez;
            break;
        case Rim3dWellLogCurve::VERTICAL_BELOW:
            normal = -Ez;
            break;
        default:
            break;
        }

        pointNormals.push_back(normal);
    }

    return pointNormals;
}