/*****************************************************************************
 * $CAMITK_LICENCE_BEGIN$
 *
 * CamiTK - Computer Assisted Medical Intervention ToolKit
 * (c) 2001-2016 Univ. Grenoble Alpes, CNRS, TIMC-IMAG UMR 5525 (GMCAO)
 *
 * Visit http://camitk.imag.fr for more information
 *
 * This file is part of CamiTK.
 *
 * CamiTK is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License version 3
 * only, as published by the Free Software Foundation.
 *
 * CamiTK 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 Lesser General Public License version 3 for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * version 3 along with CamiTK.  If not, see <http://www.gnu.org/licenses/>.
 *
 * $CAMITK_LICENCE_END$
 ****************************************************************************/

#include "CellProperties.h"
#include "Cell.h"
#include "PhysicalModel.h"

#include <math.h>

// pmlschema includes
#include <Cell.hxx>

//----------------------- Constructors -----------------------
Cell::Cell (PhysicalModel *p, const StructureProperties::GeometricType t)
    : StructuralComponent (p) {
    this->Structure::properties = new CellProperties (p, t);
    parentSC = NULL;
}

Cell::Cell (PhysicalModel *p, physicalModel::Cell xmlCell, StructuralComponent* sc) : StructuralComponent (p) {
    // Add the basic cell structure properties
    StructureProperties::GeometricType cellType = StructureProperties::INVALID;
    switch (xmlCell.cellProperties().type()) {
    case physicalModel::GeometricType::ATOM:
        cellType = StructureProperties::ATOM;
        break;
    case physicalModel::GeometricType::LINE:
        cellType = StructureProperties::LINE;
        break;
    case physicalModel::GeometricType::TRIANGLE:
        cellType = StructureProperties::TRIANGLE;
        break;
    case physicalModel::GeometricType::QUAD:
        cellType = StructureProperties::QUAD;
        break;
    case physicalModel::GeometricType::TETRAHEDRON:
        cellType = StructureProperties::TETRAHEDRON;
        break;
    case physicalModel::GeometricType::WEDGE:
        cellType = StructureProperties::WEDGE;
        break;
    case physicalModel::GeometricType::PYRAMID:
        cellType = StructureProperties::PYRAMID;
        break;
    case physicalModel::GeometricType::HEXAHEDRON:
        cellType = StructureProperties::HEXAHEDRON;
        break;
    case physicalModel::GeometricType::POLY_LINE:
        cellType = StructureProperties::POLY_LINE;
        break;
    case physicalModel::GeometricType::POLY_VERTEX:
        cellType = StructureProperties::POLY_VERTEX;
        break;
    default: // INVALID
        break;
    }
    this->Structure::properties = new CellProperties (p, cellType, xmlCell.cellProperties());

    // its direct parent component
    parentSC = sc;

    // number of structures (optional)
    if (xmlCell.nrOfStructures().present())
        plannedNumberOfStructures (xmlCell.nrOfStructures().get().value());

    // color (optional)
    if (xmlCell.color().present())
        setColor (xmlCell.color().get().r(),
                  xmlCell.color().get().g(),
                  xmlCell.color().get().b(),
                  xmlCell.color().get().a());

    // parse the cell's atom references
    physicalModel::Cell::atomRef_sequence xmlAtoms = xmlCell.atomRef();
    for (physicalModel::Cell::atomRef_iterator atomRefIt =  xmlAtoms.begin();
            atomRefIt != xmlAtoms.end();
            atomRefIt++) {
        physicalModel::AtomRef xmlAtomRef = (*atomRefIt);
        Atom* a = p->getAtom (xmlAtomRef.index());
        if (a)
            addStructure (a);
        else
            std::cerr << "Cell::Cell: Cannot find Atom for reference: " << xmlAtomRef.index() << std::endl;
    }
}

Cell::Cell (PhysicalModel *p, const StructureProperties::GeometricType t, const unsigned int ind)
    : StructuralComponent (p) {
    this->Structure::properties = new CellProperties (p, t, ind);
}

//----------------------- Destructor -----------------------
Cell::~Cell() {
    // delete the structural component properties
    deleteProperties();

    // if (Structure::properties)
    delete (CellProperties*) Structure::properties;

    if (StructuralComponent::atomList)
        delete StructuralComponent::atomList;

    StructuralComponent::atomList = NULL;

    // delete all children
    deleteAllStructures();

    // tell all parents that I am going away to the paradise of pointers
    removeFromParents();
}

//----------------------- setIndex -----------------------
bool Cell::setIndex (const unsigned int index) {
    // set the property
    Structure::setIndex (index);
    // tell the physical model about the change (and return true if insertion was ok)
    return getPhysicalModel()->addGlobalIndexCellPair (std::GlobalIndexStructurePair (index, this));
}

// --------------- makePrintData ---------------
bool Cell::makePrintData (const StructuralComponent * sc) {
    // if the sc (=the object that is calling this method) is the first in the
    // mySC list to be a exclusive component, it means the cell has not been written yet. If it is not
    // the first exclusive component in the mySC list, the cell could simply writes its cell ref
    bool isExclusive = false;
    // search the first exclusive
    unsigned int i = 0;

    while (!isExclusive && i < getNumberOfStructuralComponents()) {
        isExclusive = getStructuralComponent (i)->isExclusive();

        if (!isExclusive)
            i++;
    }

    // post condition here:
    // - either isExclusive && i is the index of the first exclusive
    // - or !isExclusive && there are no exclusive sc containing this cell

    // print the cell data (and not the cellRef) only if :
    // - the first exclusive component is sc
    // - there are no exclusive components for this cell and sc is the first SC
    if (isExclusive)
        return (sc == getStructuralComponent (i));
    else
        return (sc == getStructuralComponent (0));
}

// --------------- xmlPrint ---------------
void Cell::xmlPrint (std::ostream &o, const StructuralComponent *sc) {

    // do what's our duty
    if (makePrintData (sc)) {
        // print out all the information...
        o << "<cell>" << std::endl;

        // ...the properties...
        ( (CellProperties *) Structure::properties)->xmlPrint (o);

        // ... and the color (if it is not the default one) ...

        if ( ( (StructuralComponentProperties *) StructuralComponent::properties)->getColor() != StructuralComponentProperties::DEFAULT) {
            o << "<color r=\"" << ( (StructuralComponentProperties *) StructuralComponent::properties)->getRed();
            o << "\" g=\"" << ( (StructuralComponentProperties *) StructuralComponent::properties)->getGreen();
            o << "\" b=\"" << ( (StructuralComponentProperties *) StructuralComponent::properties)->getBlue();
            o << "\" a=\"" << ( (StructuralComponentProperties *) StructuralComponent::properties)->getAlpha() << "\"/>" << std::endl;
        }

        // optimize the memory allocation for reading
        o << "<nrOfStructures value=\"" << structures.size() << "\"/>" << std::endl;

        // ... and finally the ref to the atoms
        for (unsigned int i = 0; i < structures.size(); i++) {
            o << "<atomRef index=\"" << structures[i]->getIndex() << "\"/>" << std::endl;
        }

        o << "</cell>" << std::endl;
    } else {
        // print only out the cellRef
        o << "<cellRef index=\"" << getIndex() << "\"/>" << std::endl;
    }
}

// ------------------ deleteAllStructures ---------------------
void Cell::deleteAllStructures() {
    // here the list elements (i.e. atoms) are NOT deleted
    // but the list is cleared and cell removed from sc list of atoms
    for (unsigned int i = 0; i < structures.size(); i++) {
        Atom* a= (Atom*) structures[i];
        a->removeStructuralComponent (this);
    }
    structures.clear();
}

// ------------------ getProperties ---------------------
CellProperties * Cell::getProperties() {
    return ( (CellProperties *) Structure::properties);
}

// ------------------ getTriangleFaces ---------------------
StructuralComponent * Cell::getTriangleFaces() {

    StructuralComponent * facets = new StructuralComponent (getPhysicalModel());
    Cell * face;

    switch (getProperties()->getType()) {
    case StructureProperties::HEXAHEDRON :
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (2));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (1));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (7));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (7));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (5));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (5));
        face->addStructure (getStructure (4));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (7));
        face->addStructure (getStructure (6));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (6));
        face->addStructure (getStructure (2));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (6));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (6));
        face->addStructure (getStructure (5));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (5));
        face->addStructure (getStructure (6));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (6));
        face->addStructure (getStructure (7));
        facets->addStructure (face);
        return facets;
    case StructureProperties::TETRAHEDRON :
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (2));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (1));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        return facets;
    case StructureProperties::WEDGE :
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (5));
        face->addStructure (getStructure (4));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (1));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (4));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (5));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (1));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (5));
        facets->addStructure (face);
        return facets;
    case StructureProperties::QUAD :
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (2));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::TRIANGLE);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        return facets;
    case StructureProperties::TRIANGLE :
        facets->addStructure (this);
        return facets;
    default :
        return NULL;
    }
}

// ------------------ getQuadFaces ---------------------
StructuralComponent * Cell::getQuadFaces() {

    StructuralComponent * facets = new StructuralComponent (getPhysicalModel());
    Cell * face;

    switch (getProperties()->getType()) {
    case StructureProperties::HEXAHEDRON :
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (1));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (7));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (5));
        face->addStructure (getStructure (4));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (7));
        face->addStructure (getStructure (6));
        face->addStructure (getStructure (2));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (6));
        face->addStructure (getStructure (5));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (5));
        face->addStructure (getStructure (6));
        face->addStructure (getStructure (7));
        facets->addStructure (face);
        return facets;
    case StructureProperties::WEDGE :
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (5));
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (1));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (1));
        face->addStructure (getStructure (4));
        face->addStructure (getStructure (3));
        facets->addStructure (face);
        face = new Cell (getPhysicalModel(), StructureProperties::QUAD);
        face->addStructure (getStructure (2));
        face->addStructure (getStructure (0));
        face->addStructure (getStructure (3));
        face->addStructure (getStructure (5));
        facets->addStructure (face);
        return facets;
    case StructureProperties::QUAD :
        facets->addStructure (this);
        return facets;
    default :
        return NULL;
    }

}

// ------------------ normal ---------------------
double * Cell::normal() {
    if (getProperties()->getType() == StructureProperties::QUAD || getProperties()->getType() == StructureProperties::TRIANGLE) {
        double * N = new double[3];
        double v1[3], v2[3];
        double posi[3], posip1[3];
        ( (Atom*) getStructure (0))->getPosition (posip1);
        ( (Atom*) getStructure (2))->getPosition (posi);
        v1[0] = posi[0] - posip1[0];
        v1[1] = posi[1] - posip1[1];
        v1[2] = posi[2] - posip1[2];
        ( (Atom*) getStructure (1))->getPosition (posi);
        v2[0] = posi[0] - posip1[0];
        v2[1] = posi[1] - posip1[1];
        v2[2] = posi[2] - posip1[2];

        N[0] = v1[1] * v2[2] - v1[2] * v2[1];
        N[1] = v1[2] * v2[0] - v1[0] * v2[2];
        N[2] = v1[0] * v2[1] - v1[1] * v2[0];

        double norm = sqrt (N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);

        N[0] /= norm;
        N[1] /= norm;
        N[2] /= norm;
        return N;
    }

    return NULL;

}

// ------------------ surface ---------------------
double Cell::surface() {
    if (getProperties()->getType() == StructureProperties::TRIANGLE) {
        double A[3] = {0.0, 0.0, 0.0 };
        unsigned int nbElem;

        nbElem = Cell::getNumberOfStructures();

        double posi[3], posip1[3];
        ( (Atom*) getStructure (0))->getPosition (posip1);

        for (unsigned int i = 0; i < nbElem; i++) {
            posi[0] = posip1[0];
            posi[1] = posip1[1];
            posi[2] = posip1[2];
            ( (Atom*) getStructure ( (i + 1) % nbElem))->getPosition (posip1);
            //Cross(posi, posip1, inter);
            A[0] += posi[1] * posip1[2] - posi[2] * posip1[1];
            A[1] += posi[2] * posip1[0] - posi[0] * posip1[2];
            A[2] += posi[0] * posip1[1] - posi[1] * posip1[0];
        }

        // here A is in fact twice the theoritical area vector
        A[0] /= 2.0;
        A[1] /= 2.0;
        A[2] /= 2.0;

        // face normal :
        double * N = normal();

        double surface = N[0] * A[0] + N[1] * A[1] + N[2] * A[2];

        return surface > 0 ? surface : -surface;
    } else {
        StructuralComponent * facets = getTriangleFaces();

        if (!facets)
            return 0.0;

        double surface = 0.0;

        for (unsigned int i = 0 ; i < facets->getNumberOfCells() ; i++) {
            surface += facets->getCell (i)->surface();
        }
        delete facets;
        return surface;
    }
}


// ------------------ volume ---------------------
double Cell::volume() {
    StructuralComponent * facets = getTriangleFaces();

    if (!facets || getProperties()->getType() == StructureProperties::QUAD /* || getProperties()->getType() == StructureProperties::TRIANGLE*/)
        return 0.0;

    double vol = 0.0;

    Cell * face;

    // volume can be computed with an closed triangular mesh with all normals pointing outwards
    if(getProperties()->getType() == StructureProperties::TRIANGLE) {

        //TODO: check if the mesh is closed with normals pointing outwards ?
        for (unsigned int i = 0; i < facets->getNumberOfCells(); i++) {
            face = facets->getCell (i);

            double p1[3], p2[3], p3[3];
            ( (Atom*) face->getStructure (0))->getPosition (p1);
            ( (Atom*) face->getStructure (1))->getPosition (p2);
            ( (Atom*) face->getStructure (2))->getPosition (p3);
            vol+= face->signedVolumeTriangle(p1,p2,p3);
        }
        // do not check if the volume is negative because it CAN be negative
        // the volume of an enclosed triangular mesh is the sum of positive and negative
        // volumes of teatrahedra from origin to each triangular face.


    } else {

        for (unsigned int i = 0; i < facets->getNumberOfCells(); i++) {
            face = facets->getCell (i);
            double pos[3];
            ( (Atom*) face->getStructure (0))->getPosition (pos);
            double * N = face->normal();
            vol += face->surface() * (pos[0] * N[0] + pos[1] * N[1] + pos[2] * N[2]);
        }
        vol /= 3.0;
        if (vol < 0)
            vol = -vol;
    }

    delete facets;
    return vol;

}


//-----------signedVolumeTriangle--------------------------------
double Cell::signedVolumeTriangle(double p1[3], double p2[3], double p3[3]) {

    double v321 = p3[0]*p2[1]*p1[2];
    double v231 = p2[0]*p3[1]*p1[2];
    double v312 = p3[0]*p1[1]*p2[2];
    double v132 = p1[0]*p3[1]*p2[2];
    double v213 = p2[0]*p1[1]*p3[2];
    double v123 = p1[0]*p2[1]*p3[2];
    double vol = (1.0/6.0)*(-v321 + v231 + v312 - v132 - v213 + v123);
    return vol;
}