/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2018-     Equinor ASA
//  Copyright (C) 2018-     Ceetron Solutions AS
//
//  Adapted from work by Paul D. Bourke named "conrec"
//
//  http://paulbourke.net/papers/conrec/.
//
//  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 "cafContourLines.h"

#include <algorithm>
#include <cmath>
#include <list>

const int caf::ContourLines::s_castab[3][3][3] = { { { 0, 0, 8 }, { 0, 2, 5 }, { 7, 6, 9 } },
                                                   { { 0, 3, 4 }, { 1, 3, 1 }, { 4, 3, 0 } },
                                                   { { 9, 6, 7 }, { 5, 2, 0 }, { 8, 0, 0 } } };

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void caf::ContourLines::create( const std::vector<double>&            dataXY,
                                const std::vector<double>&            xCoords,
                                const std::vector<double>&            yCoords,
                                const std::vector<double>&            contourLevels,
                                std::vector<std::vector<cvf::Vec2d>>* polygons )
{
    CVF_ASSERT( !contourLevels.empty() );
    int                 nContourLevels = static_cast<int>( contourLevels.size() );
    std::vector<int>    sh( 5, 0 );
    std::vector<double> h( 5, 0.0 ), xh( 5, 0.0 ), yh( 5, 0.0 );

    int nx = static_cast<int>( xCoords.size() );
    int ny = static_cast<int>( yCoords.size() );

    CVF_ASSERT( static_cast<int>( dataXY.size() ) == nx * ny );

    polygons->resize( nContourLevels );

    int im[4] = { 0, 1, 1, 0 }, jm[4] = { 0, 0, 1, 1 };

    for ( int j = ( ny - 2 ); j >= 0; j-- )
    {
        for ( int i = 0; i < nx - 1; i++ )
        {
            double temp1, temp2;
            temp1       = std::min( saneValue( gridIndex1d( i, j, nx ), dataXY, contourLevels ),
                              saneValue( gridIndex1d( i, j + 1, nx ), dataXY, contourLevels ) );
            temp2       = std::min( saneValue( gridIndex1d( i + 1, j, nx ), dataXY, contourLevels ),
                              saneValue( gridIndex1d( i + 1, j + 1, nx ), dataXY, contourLevels ) );
            double dmin = std::min( temp1, temp2 );
            temp1       = std::max( saneValue( gridIndex1d( i, j, nx ), dataXY, contourLevels ),
                              saneValue( gridIndex1d( i, j + 1, nx ), dataXY, contourLevels ) );
            temp2       = std::max( saneValue( gridIndex1d( i + 1, j, nx ), dataXY, contourLevels ),
                              saneValue( gridIndex1d( i + 1, j + 1, nx ), dataXY, contourLevels ) );
            double dmax = std::max( temp1, temp2 );
            // Using dmax <= contourLevels[0] as a deviation from Bourke because it empirically
            // Reduces gridding artifacts in our code.
            if ( dmax <= contourLevels[0] || dmin > contourLevels[nContourLevels - 1] ) continue;

            for ( int k = 0; k < nContourLevels; k++ )
            {
                if ( contourLevels[k] < dmin || contourLevels[k] > dmax ) continue;
                for ( int m = 4; m >= 0; m-- )
                {
                    if ( m > 0 )
                    {
                        double value = saneValue( gridIndex1d( i + im[m - 1], j + jm[m - 1], nx ), dataXY, contourLevels );
                        if ( value == invalidValue( contourLevels ) )
                        {
                            h[m] = invalidValue( contourLevels );
                        }
                        else
                        {
                            h[m] = value - contourLevels[k];
                        }
                        xh[m] = xCoords[i + im[m - 1]];
                        yh[m] = yCoords[j + jm[m - 1]];
                    }
                    else
                    {
                        h[0]  = 0.25 * ( h[1] + h[2] + h[3] + h[4] );
                        xh[0] = 0.5 * ( xCoords[i] + xCoords[i + 1] );
                        yh[0] = 0.5 * ( yCoords[j] + yCoords[j + 1] );
                    }
                    if ( h[m] > 0.0 )
                        sh[m] = 1;
                    else if ( h[m] < 0.0 )
                        sh[m] = -1;
                    else
                        sh[m] = 0;
                }

                /*
                   Note: at this stage the relative heights of the corners and the
                   centre are in the h array, and the corresponding coordinates are
                   in the xh and yh arrays. The centre of the box is indexed by 0
                   and the 4 corners by 1 to 4 as shown below.
                   Each triangle is then indexed by the parameter m, and the 3
                   vertices of each triangle are indexed by parameters m1,m2,and m3.
                   It is assumed that the centre of the box is always vertex 2
                   though this isimportant only when all 3 vertices lie exactly on
                   the same contour level, in which case only the side of the box
                   is drawn.
                      vertex 4 +-------------------+ vertex 3
                               | \               / |
                               |   \    m-3    /   |
                               |     \       /     |
                               |       \   /       |
                               |  m=2    X   m=2   |       the centre is vertex 0
                               |       /   \       |
                               |     /       \     |
                               |   /    m=1    \   |
                               | /               \ |
                      vertex 1 +-------------------+ vertex 2
                */
                /* Scan each triangle in the box */
                for ( int m = 1; m <= 4; m++ )
                {
                    int m1 = m;
                    int m2 = 0;
                    int m3 = ( m != 4 ) ? m + 1 : 1;

                    double x1 = 0.0, x2 = 0.0, y1 = 0.0, y2 = 0.0;
                    int    case_value = s_castab[sh[m1] + 1][sh[m2] + 1][sh[m3] + 1];
                    if ( case_value == 0 ) continue;

                    switch ( case_value )
                    {
                        case 1: /* Line between vertices 1 and 2 */
                            x1 = xh[m1];
                            y1 = yh[m1];
                            x2 = xh[m2];
                            y2 = yh[m2];
                            break;
                        case 2: /* Line between vertices 2 and 3 */
                            x1 = xh[m2];
                            y1 = yh[m2];
                            x2 = xh[m3];
                            y2 = yh[m3];
                            break;
                        case 3: /* Line between vertices 3 and 1 */
                            x1 = xh[m3];
                            y1 = yh[m3];
                            x2 = xh[m1];
                            y2 = yh[m1];
                            break;
                        case 4: /* Line between vertex 1 and side 2-3 */
                            x1 = xh[m1];
                            y1 = yh[m1];
                            x2 = xsect( m2, m3, h, xh, yh );
                            y2 = ysect( m2, m3, h, xh, yh );
                            break;
                        case 5: /* Line between vertex 2 and side 3-1 */
                            x1 = xh[m2];
                            y1 = yh[m2];
                            x2 = xsect( m3, m1, h, xh, yh );
                            y2 = ysect( m3, m1, h, xh, yh );
                            break;
                        case 6: /* Line between vertex 3 and side 1-2 */
                            x1 = xh[m3];
                            y1 = yh[m3];
                            x2 = xsect( m1, m2, h, xh, yh );
                            y2 = ysect( m1, m2, h, xh, yh );
                            break;
                        case 7: /* Line between sides 1-2 and 2-3 */
                            x1 = xsect( m1, m2, h, xh, yh );
                            y1 = ysect( m1, m2, h, xh, yh );
                            x2 = xsect( m2, m3, h, xh, yh );
                            y2 = ysect( m2, m3, h, xh, yh );
                            break;
                        case 8: /* Line between sides 2-3 and 3-1 */
                            x1 = xsect( m2, m3, h, xh, yh );
                            y1 = ysect( m2, m3, h, xh, yh );
                            x2 = xsect( m3, m1, h, xh, yh );
                            y2 = ysect( m3, m1, h, xh, yh );
                            break;
                        case 9: /* Line between sides 3-1 and 1-2 */
                            x1 = xsect( m3, m1, h, xh, yh );
                            y1 = ysect( m3, m1, h, xh, yh );
                            x2 = xsect( m1, m2, h, xh, yh );
                            y2 = ysect( m1, m2, h, xh, yh );
                            break;
                        default:
                            break;
                    }

                    /* Finally draw the line */
                    polygons->at( k ).push_back( cvf::Vec2d( x1, y1 ) );
                    polygons->at( k ).push_back( cvf::Vec2d( x2, y2 ) );
                } /* m */
            } /* k - contour */
        } /* i */
    } /* j */
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<caf::ContourLines::ListOfLineSegments> caf::ContourLines::create( const std::vector<double>& dataXY,
                                                                              const std::vector<double>& xPositions,
                                                                              const std::vector<double>& yPositions,
                                                                              const std::vector<double>& contourLevels )
{
    std::vector<std::vector<cvf::Vec2d>> contourLineSegments;
    caf::ContourLines::create( dataXY, xPositions, yPositions, contourLevels, &contourLineSegments );

    std::vector<ListOfLineSegments> listOfSegmentsPerLevel( contourLevels.size() );

    for ( size_t i = 0; i < contourLevels.size(); ++i )
    {
        size_t nPoints   = contourLineSegments[i].size();
        size_t nSegments = nPoints / 2;
        if ( nSegments >= 3u ) // Need at least three segments for a closed polygon
        {
            ListOfLineSegments unorderedSegments;
            for ( size_t j = 0; j < contourLineSegments[i].size(); j += 2 )
            {
                unorderedSegments.push_back( std::make_pair( cvf::Vec3d( contourLineSegments[i][j] ),
                                                             cvf::Vec3d( contourLineSegments[i][j + 1] ) ) );
            }
            listOfSegmentsPerLevel[i] = unorderedSegments;
        }
    }

    return listOfSegmentsPerLevel;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double caf::ContourLines::contourRange( const std::vector<double>& contourLevels )
{
    CVF_ASSERT( !contourLevels.empty() );
    return std::max( 1.0e-6, contourLevels.back() - contourLevels.front() );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double caf::ContourLines::invalidValue( const std::vector<double>& contourLevels )
{
    return contourLevels.front() - 1000.0 * contourRange( contourLevels );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double caf::ContourLines::saneValue( int index, const std::vector<double>& dataXY, const std::vector<double>& contourLevels )
{
    CVF_ASSERT( index >= 0 && index < static_cast<int>( dataXY.size() ) );

    // Place all invalid values below the bottom contour level.
    if ( dataXY[index] == -std::numeric_limits<double>::infinity() ||
         dataXY[index] == std::numeric_limits<double>::infinity() )
    {
        return invalidValue( contourLevels );
    }
    return dataXY[index];
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double caf::ContourLines::xsect( int                        p1,
                                 int                        p2,
                                 const std::vector<double>& h,
                                 const std::vector<double>& xh,
                                 const std::vector<double>& yh )
{
    return ( h[p2] * xh[p1] - h[p1] * xh[p2] ) / ( h[p2] - h[p1] );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double caf::ContourLines::ysect( int                        p1,
                                 int                        p2,
                                 const std::vector<double>& h,
                                 const std::vector<double>& xh,
                                 const std::vector<double>& yh )
{
    return ( h[p2] * yh[p1] - h[p1] * yh[p2] ) / ( h[p2] - h[p1] );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int caf::ContourLines::gridIndex1d( int i, int j, int nx )
{
    return j * nx + i;
}