#include "lvox3_stepcomputepad.h"

//In/Out
#include "ct_result/model/inModel/ct_inresultmodelgroup.h"
#include "ct_result/model/outModel/ct_outresultmodelgroup.h"
#include "ct_result/ct_resultgroup.h"
#include "ct_model/tools/ct_modelsearchhelper.h"
#include "ct_model/inModel/tools/ct_instdmodelpossibility.h"
#include "ct_result/model/inModel/ct_inresultmodelgrouptocopy.h"
#include "ct_result/model/outModel/tools/ct_outresultmodelgrouptocopypossibilities.h"

// Tools
#include "ct_view/ct_stepconfigurabledialog.h"
#include "ct_turn/inTurn/tools/ct_inturnmanager.h"
#include "tools/lvox3_utils.h"

// Drawables
#include "tools/lvox3_gridtype.h"
#include "qdebug.h"

// Alias for indexing in models
#define DEF_resultIn_grids "r"
#define DEF_groupIn_grids "gr"
#define DEF_itemIn_deltaHits "deltaHits"
#define DEF_itemIn_deltaTheoritical "deltaTheoritical"
#define DEF_itemIn_deltaBefore "deltaBefore"

// Alias for indexing out models
#define DEF_SearchInResult "r"
#define DEF_itemOut_PADCF "PAD_CF"
#define DEF_itemOut_PADMLE "PAD_MLE"
#define DEF_itemOut_PADBiasCorMLE "PAD_BiasCorMLE"
#define DEF_itemOut_PADBiasCorBL "PAD_BiasCorBL"
#define DEF_itemOut_PADBL "PAD_BL"
#define DEF_itemOut_PADUnePathBiasCorrBL "PAD_UnePathBiasCorBL"

#define DEF_itemOut_statDeltaSum "DeltaSum"
#define DEF_itemOut_statDeltaEffSum "DeltaEffSum"
#define DEF_itemOut_statDeltaEffSquareSum "DeltaEffSquareSum"
#define DEF_itemOut_statFreeSum "FreeSum"
#define DEF_itemOut_statFreeEffSum "FreeEffSum"
#define DEF_itemOut_stat1_zleDelta_effSum "1_zleDelta_effSum"


#include <math.h>
//# define M_PI           3.14159265358979323846  /* pi ?*/

// Constructor : initialization of parameters
LVOX3_StepComputePAD::LVOX3_StepComputePAD(CT_StepInitializeData &dataInit) : CT_AbstractStep(dataInit)
{
    _CF = false;
    _MLE = false;
    _BCMLE = false;
    _BL = false;
    _BCBL = false;
    _BC2BL = false;
    _EraseInfinity = false;
    _PADlimit = 100; // 100 m2/voxel
    _NThreshold = 5;
    _isRelativeThreshold = false;
    _grilleSup = false;
}

// Step description (tooltip of contextual menu)
QString LVOX3_StepComputePAD::getStepDescription() const
{
    return tr("3b - Estimation du Plant Area Density");
}

// Step copy method
CT_VirtualAbstractStep* LVOX3_StepComputePAD::createNewInstance(CT_StepInitializeData &dataInit)
{
    return new LVOX3_StepComputePAD(dataInit);
}

//////////////////// PROTECTED METHODS //////////////////

// Creation and affiliation of IN models
void LVOX3_StepComputePAD::createInResultModelListProtected()
{
    CT_InResultModelGroupToCopy *resultModel = createNewInResultModelForCopy(
            DEF_SearchInResult, tr("Grids"), "", true);

    resultModel->setZeroOrMoreRootGroup();
    resultModel->addGroupModel("", DEF_groupIn_grids, CT_AbstractItemGroup::staticGetType(), tr("Group"));

    resultModel->addItemModel(DEF_groupIn_grids, DEF_itemIn_deltaHits, LVOX3_Grid3D<int>::staticGetType(), tr("Stats hits"));
    resultModel->addItemModel(DEF_groupIn_grids, DEF_itemIn_deltaTheoritical, LVOX3_Grid3D<int>::staticGetType(), tr("Stats th."));
    resultModel->addItemModel(DEF_groupIn_grids, DEF_itemIn_deltaBefore, LVOX3_Grid3D<int>::staticGetType(), tr("Stats before"));

}

// Creation and affiliation of OUT models
void LVOX3_StepComputePAD::createOutResultModelListProtected()
{
    CT_OutResultModelGroupToCopyPossibilities *resultModel = createNewOutResultModelToCopy(DEF_SearchInResult);

    if(resultModel != NULL){
        if(_CF)
            resultModel->addItemModel(DEF_groupIn_grids, _PADCF_ModelName, new lvox::Grid3Df(), tr("PAD-ContactFreq"));
        if(_MLE)
            resultModel->addItemModel(DEF_groupIn_grids, _PADMLE_ModelName, new lvox::Grid3Df(), tr("PAD-MLE"));
        if(_BCMLE)
            resultModel->addItemModel(DEF_groupIn_grids, _PADBCMLE_ModelName, new lvox::Grid3Df(), tr("PAD-BiasCorMLE"));
        if(_BL)
            resultModel->addItemModel(DEF_groupIn_grids, _PADBL_ModelName, new lvox::Grid3Df(), tr("PAD-BeerLamb"));
        if(_BCBL)
            resultModel->addItemModel(DEF_groupIn_grids, _PADBCBL_ModelName, new lvox::Grid3Df(), tr("PAD-BiasCorBeerLamb"));
        if(_BC2BL)
            resultModel->addItemModel(DEF_groupIn_grids, _PADBC2BL_ModelName, new lvox::Grid3Df(), tr("PAD-UnequalPathBeerLamb"));
        if(_grilleSup)
        {
            resultModel->addItemModel(DEF_groupIn_grids, _FreeSum_ModelName, new lvox::Grid3Df(), tr("FreePath Sum"));
            resultModel->addItemModel(DEF_groupIn_grids, _FreeEffSum_ModelName, new lvox::Grid3Df(), tr("Effective FreePath Sum"));
            resultModel->addItemModel(DEF_groupIn_grids, _DeltaSum_ModelName, new lvox::Grid3Df(), tr("Delta Sum"));
            resultModel->addItemModel(DEF_groupIn_grids, _DeltaEffSum_ModelName, new lvox::Grid3Df(), tr("Effective Delta Sum"));
            resultModel->addItemModel(DEF_groupIn_grids, _DeltaEffSquareSum_ModelName, new lvox::Grid3Df(), tr("Effective Delta Square Sum"));
            resultModel->addItemModel(DEF_groupIn_grids, _1ZleDeltaEffSum_ModelName, new lvox::Grid3Df(), tr("1_zleDelta Effective Sum"));
        }
    }
}

// Semi-automatic creation of step parameters DialogBox
void LVOX3_StepComputePAD::createPostConfigurationDialog()
{
    CT_StepConfigurableDialog *configDialog = newStandardPostConfigurationDialog();

    //configDialog->addBool(tr("    -> Ecraser les valeurs infinies "),"", "", _EraseInfinity);
    //configDialog->addDouble(tr("   -> Effacer les valeurs > à"), "", 0, 100, 2, _PADlimit);

    configDialog->addTitle(tr("Attenuation coefficient estimators :"),"");
    configDialog->addBool(tr("Contact Frequency"),"","",_CF);
    configDialog->addBool(tr("MLE"),"","",_MLE);
    configDialog->addBool(tr("Bias-corrected MLE"),"","",_BCMLE);
    configDialog->addBool(tr("Beer-Lambert"),"","",_BL);
    configDialog->addBool(tr("Bias-corrected Beer-Lambert (equal path)"),"","",_BCBL);
    configDialog->addBool(tr("Bias-corrected Beer-Lambert (unequal path)"),"","",_BC2BL);

    configDialog->addText(tr("Configuration options :"),"","","");
    configDialog->addInt(tr("Compute estimators for N larger than"),"",0,9999,_NThreshold,"");
    configDialog->addBool(tr("Is relative threshold, use Nb/Nt instead of N"),"","",_isRelativeThreshold);
    configDialog->addBool(tr("Create Grids Corresponding to the different summations"),"","",_grilleSup);

    //configDialog->addStringChoice(tr("Niveau de confiance"),"",_confidenceLevelCollection.keys(), _confidenceLevel);
}


void LVOX3_StepComputePAD::compute()
{

    float errorVal=-1; // a attenuation coefficient of -0.5 lead to a PAD of -1 (to track error)

    QList<CT_ResultGroup*> outResultList = getOutResultList();
    CT_ResultGroup* resultOut_grids = outResultList.at(0);

    CT_ResultGroupIterator itGrp(resultOut_grids, this, DEF_groupIn_grids);
    while (itGrp.hasNext() && !isStopped())
    {
        CT_AbstractItemGroup *group = (CT_AbstractItemGroup*) itGrp.next();

        const lvox::Grid3Di* itemIn_deltaH = (const lvox::Grid3Di*) group->firstItemByINModelName(this, DEF_itemIn_deltaHits);
        const lvox::Grid3Di* itemIn_deltaT = (const lvox::Grid3Di*) group->firstItemByINModelName(this, DEF_itemIn_deltaTheoritical);
        const lvox::Grid3Di* itemIn_deltaB = (const lvox::Grid3Di*) group->firstItemByINModelName(this, DEF_itemIn_deltaBefore);

        if (itemIn_deltaH != NULL && itemIn_deltaT != NULL && itemIn_deltaB != NULL)
        {
            float xres = itemIn_deltaH->xresolution();
            float yres = itemIn_deltaH->yresolution();
            float zres = itemIn_deltaH->zresolution();
            float NAd = itemIn_deltaH->NA();

            _lambda1 = itemIn_deltaT->getLambda1();

            qDebug()<<"x array size theo before hit"<< itemIn_deltaT->xArraySize() << itemIn_deltaB->xArraySize() << itemIn_deltaH->xArraySize();
            qDebug()<<"y array size theo before hit"<< itemIn_deltaT->yArraySize() << itemIn_deltaB->yArraySize() << itemIn_deltaH->yArraySize();
            qDebug()<<"z array size theo before hit"<< itemIn_deltaT->zArraySize() << itemIn_deltaB->zArraySize() << itemIn_deltaH->zArraySize();
            qDebug()<<"min z theo before hit"<< itemIn_deltaT->minZ() << itemIn_deltaB->minZ() << itemIn_deltaH->minZ();

            size_t xdim = itemIn_deltaH->xdim();
            size_t ydim = itemIn_deltaH->ydim();
            size_t zdim = itemIn_deltaH->zdim();
            lvox::Grid3Df* itemOut_PADCF = new lvox::Grid3Df(_PADCF_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_PADMLEBias = new lvox::Grid3Df(_PADBCMLE_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_PADMLESimple = new lvox::Grid3Df(_PADMLE_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_PADBLEqual = new lvox::Grid3Df(_PADBCBL_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_PADBLBasic = new lvox::Grid3Df(_PADBL_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_PADBLUnequal = new lvox::Grid3Df(_PADBC2BL_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_deltaSum = new lvox::Grid3Df(_DeltaSum_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_deltaEffSum = new lvox::Grid3Df(_DeltaEffSum_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_zSum = new lvox::Grid3Df(_FreeSum_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_zEffSum = new lvox::Grid3Df(_FreeEffSum_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_deltaSquareSum = new lvox::Grid3Df(_DeltaEffSquareSum_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);
            lvox::Grid3Df* itemOut_1_zleDelta_effSum = new lvox::Grid3Df(_1ZleDeltaEffSum_ModelName.completeName(), resultOut_grids,
                                                                 itemIn_deltaH->minX(), itemIn_deltaH->minY(), itemIn_deltaH->minZ(),
                                                                 xdim, ydim, zdim,
                                                                 xres,yres,zres, NAd, NAd);

            if(_BCMLE)
                group->addItemDrawable(itemOut_PADMLEBias);
            if(_CF)
                group->addItemDrawable(itemOut_PADCF);
            if(_MLE)
                group->addItemDrawable(itemOut_PADMLESimple);
            if(_BCBL)
                group->addItemDrawable(itemOut_PADBLEqual);
            if(_BL)
                group->addItemDrawable(itemOut_PADBLBasic);
            if(_BC2BL)
                group->addItemDrawable(itemOut_PADBLUnequal);
            if(_grilleSup)
            {
                group->addItemDrawable(itemOut_zSum);
                group->addItemDrawable(itemOut_zEffSum);
                group->addItemDrawable(itemOut_deltaSum);
                group->addItemDrawable(itemOut_deltaEffSum);
                group->addItemDrawable(itemOut_deltaSquareSum);
                group->addItemDrawable(itemOut_1_zleDelta_effSum);
            }

            /**@TODO: VTN
            * if Ni = 0 && N > 0 then pad should be 0 instead of errorVal
            **/
            for(int i = 0; i < itemIn_deltaH->nCells(); i++)
            {
                int Nt = itemIn_deltaT->valueAtIndex(i);
                int Ni = itemIn_deltaH->valueAtIndex(i);
                int Nb = itemIn_deltaB->valueAtIndex(i);
                int N = Nt - Nb;

                //now switch to relative Nthreshold is %
                // (Nt - N)/Nt < threshold/100.0
                /*
                double passRate = static_cast<double>(N)/Nt;

                if(Ni == 0 || passRate <= _NThreshold/100.0)
                {
                    // passRate <= _NThreshold/100.0
                    double noDataOrNoMaterialValue = errorVal;
                    // Ni == 0
                    if (passRate > _NThreshold / 100.0){
                        noDataOrNoMaterialValue = 0;
                    }*/
                //shoud be called by the same function isNotOcclusion(VoxelData) as in merged grids
                if(Ni == 0 || isOccluded(N, Nt))
                {
                    double noDataOrNoMaterialValue = errorVal;
                    if (!isOccluded(N, Nt)){
                        //qDebug()<<"occluded Ni=" << Ni << " N="<< N << " Nt=" << Nt;
                        noDataOrNoMaterialValue = 0;
                    }
                    if(_BCMLE)
                        itemOut_PADMLEBias->setValueAtIndex(i, noDataOrNoMaterialValue);
                    if(_CF)
                        itemOut_PADCF->setValueAtIndex(i, noDataOrNoMaterialValue);
                    if(_MLE)
                        itemOut_PADMLESimple->setValueAtIndex(i, noDataOrNoMaterialValue);
                    if(_BCBL)
                        itemOut_PADBLEqual->setValueAtIndex(i, noDataOrNoMaterialValue);
                    if(_BL)
                        itemOut_PADBLBasic->setValueAtIndex(i, noDataOrNoMaterialValue);
                    if(_BC2BL)
                        itemOut_PADBLUnequal->setValueAtIndex(i, noDataOrNoMaterialValue);
                    if(_grilleSup)
                    {
                        itemOut_zSum->setValueAtIndex(i, noDataOrNoMaterialValue);
                        itemOut_zEffSum->setValueAtIndex(i, noDataOrNoMaterialValue);
                        itemOut_deltaSum->setValueAtIndex(i, noDataOrNoMaterialValue);
                        itemOut_deltaEffSum->setValueAtIndex(i, noDataOrNoMaterialValue);
                        itemOut_deltaSquareSum->setValueAtIndex(i, noDataOrNoMaterialValue);
                        itemOut_1_zleDelta_effSum->setValueAtIndex(i, noDataOrNoMaterialValue);
                    }

                }
                else
                {
                    //qDebug() <<itemIn_deltaT->valueAtIndex(i)<< itemIn_deltaT->getSommaDelta(i);
                    //delta: path formed by the intersection of the ray and the voxel faces
                    //z    : path formed by the intersection of the ray and objects surface z <= delta
                    double deltaSum = itemIn_deltaT->getSommaDelta(i) - itemIn_deltaB->getSommaDelta(i);/*itemIn_deltaH->getSommaDelta(i) +*/
                    double deltaEffSum = itemIn_deltaT->getSommaEffectiveDelta(i) - itemIn_deltaB->getSommaEffectiveDelta(i);/*itemIn_deltaH->getSommaEffectiveDelta(i) +*/
                    double zSum = itemIn_deltaH->getSommaFree(i) + itemIn_deltaT->getSommaDelta(i) - itemIn_deltaB->getSommaDelta(i) - itemIn_deltaH->getSommaDelta(i);
                    double zEffSum = itemIn_deltaH->getSommaEffectiveFree(i) + itemIn_deltaT->getSommaEffectiveDelta(i) - itemIn_deltaB->getSommaEffectiveDelta(i) - itemIn_deltaH->getSommaEffectiveDelta(i);
                    double deltaSquareSum = itemIn_deltaT->getSommaDeltaSquare(i)-itemIn_deltaB->getSommaDeltaSquare(i);/*itemIn_deltaH->getSommaDeltaSquare(i)+*/
                    //qDebug() <<itemIn_deltaT->getSommaDeltaSquare(i)<<itemIn_deltaB->getSommaDeltaSquare(i);                    
                    double I = static_cast<double>(Ni)/ N;
                    //qDebug() << I << itemIn_deltaH->valueAtIndex(i) << N;
                    double deltaMean = deltaSum/N;
                    double freeMean = zSum/N;
                    /*
                    if (zSum < 0){
                        qDebug() <<"###"<<i<< " Hit free: " << itemIn_deltaH->getSommaFree(i) << " Theory delta: " << itemIn_deltaT->getSommaDelta(i) << " Before delta: "<< -itemIn_deltaB->getSommaDelta(i) << " Hits delta: "<< -itemIn_deltaH->getSommaDelta(i);
                    }*/
                    double effectiveFreePath = zEffSum/N;
                    double effectiveDeltaPath = deltaEffSum/N;
                    double effectiveFreePathHits = itemIn_deltaH->getSommaEffectiveHits(i)/N;

                    double deltaVarEff = deltaSquareSum/N - effectiveDeltaPath*effectiveDeltaPath;
                    double alpha = deltaVarEff / effectiveDeltaPath;
                    //qDebug() << "alpha="<<_alpha<<deltaSquareSum/N<<_effectiveDeltaPath*_effectiveDeltaPath;
                    //qDebug() << "alpha="<<_alpha<<deltaSquareSum/N<<_deltaMean;
                    //qDebug() << I << deltaMean << freeMean << effectiveFreePath << effectiveDeltaPath << effectiveFreePathHits;


                    //Calculs des Eq.5.8.25.27.32
                    // FP Modified this section to make sure inconsistent statistics lead to PAD of -1 (with G=0.5)
                    double lambdaCF = I/deltaMean;
                    if(deltaMean<=0)
                        lambdaCF = errorVal*0.5;
                    //Eq. 5
                    double lambdaHat = -((log(1- I))/deltaMean);
                    if(isinf(lambdaHat)||deltaMean<=0)
                        lambdaHat = errorVal*0.5;
                    //Eq. 8
                    double lambdaTilda = I / freeMean;
                    if(freeMean <= 0)
                        lambdaTilda = errorVal*0.5;

                    //Eq. 25 : Beer Lambert
                    double capLambdaHat = errorVal*0.5;
                    if(effectiveDeltaPath > 0) {
                        if(I < 1)   {
                            capLambdaHat = (-log(1 - I) - I/(2*N*(1 - I)))/effectiveDeltaPath;
                        } else if(I >= 1) {
                            capLambdaHat = log(2*N + 2)/effectiveDeltaPath;
                        }
                        /*
                        else { //case I>1...
                            //VTN: in case of I > 1
                            capLambdaHat = errorVal*0.5;
                        }*/
                    }
                    // Eq. 27
                    double capLambdaHat2 = capLambdaHat;
                    if (capLambdaHat > errorVal*0.5) {
                        if (1-2*alpha*capLambdaHat >= 0.05)   { // case of correction
                            capLambdaHat2 = 1/alpha*(1 - sqrt(1 - 2*alpha*capLambdaHat));
                        }
                        /*not need
                        else {
                            capLambdaHat2 = capLambdaHat;
                        }*/
                    }

                    //qDebug() <<Ni<< Nb<< Nt<< N<< I<< "effectiveFreePath: "<< effectiveFreePath;
                    //qDebug() << "_effectiveFreePath_Hits: "<< _effectiveFreePath_Hits;
                    //qDebug() << "_deltaMean: " << _deltaMean;
                    //qDebug()<< "Ni:" <<Ni<<"Nb:" << Nb<<"Nt:" << Nt<<"N:" << N<<"RDI:" << I<< "Beer Lambert:"<< capLambdaHat2 << "effectiveDeltaPath:" << effectiveDeltaPath << "deltaEffSum" << deltaEffSum;
                    //qDebug() << "Theory Eff delta" << itemIn_deltaT->getSommaEffectiveDelta(i) << "Before Eff delta"<< itemIn_deltaB->getSommaEffectiveDelta(i);
                    /**
                     * Van-Tho Nguyen
                     * Refactor to avoid divison by zero or very small double number
                     **/
                    double capLambdaTilda = 0;
                    if (effectiveFreePath <= std::numeric_limits<double>::epsilon()){
                        capLambdaTilda = errorVal*0.5;
                        Eigen::Vector3d cCenter;
                        itemIn_deltaH->getCellCenterCoordinates(i, cCenter);
                        qDebug() <<"#vid:"<<i<< "v coords:"<< cCenter.x()<< cCenter.y()<< cCenter.z()
                                 << "Z:"<< effectiveFreePath
                                 << "Hit eff free:" << itemIn_deltaH->getSommaEffectiveFree(i)
                                 << "Theory eff delta:" << itemIn_deltaT->getSommaEffectiveDelta(i)
                                 << "Before eff delta:"<< -itemIn_deltaB->getSommaEffectiveDelta(i)
                                 << "Hits effdelta: "<< -itemIn_deltaH->getSommaEffectiveDelta(i)
                                 << "Nt ni nb:"<< itemIn_deltaT->valueAtIndex(i) << itemIn_deltaH->valueAtIndex(i)<<itemIn_deltaB->valueAtIndex(i)
                                 << "nt - ni -nb:" << itemIn_deltaT->valueAtIndex(i) - itemIn_deltaH->valueAtIndex(i) - itemIn_deltaB->valueAtIndex(i);
                    }else{

                        capLambdaTilda = I / effectiveFreePath - effectiveFreePathHits / ( N*effectiveFreePath * effectiveFreePath );
                        if (capLambdaTilda < 0){
                            capLambdaTilda = 0; // force to be equal to zero if bias correction lead to negative value (this is not an error here)
                        }
                    }

                    //qDebug() << capLambdaTilda << lambdaTilda << capLambdaHat << lambdaHat << capLambdaHat2;

                    // Attenuation coefficients are converted to PAD assuming a G function equal to 0.5
                    if(_CF)itemOut_PADCF->setValueAtIndex(i,lambdaCF * (1/0.5));
                    if(_BCMLE)itemOut_PADMLEBias->setValueAtIndex(i,capLambdaTilda * (1/0.5));
                    if(_MLE)itemOut_PADMLESimple->setValueAtIndex(i,lambdaTilda* (1/0.5));
                    if(_BCBL)itemOut_PADBLEqual->setValueAtIndex(i,capLambdaHat* (1/0.5));
                    if(_BL)itemOut_PADBLBasic->setValueAtIndex(i,lambdaHat* (1/0.5));
                    if(_BC2BL)itemOut_PADBLUnequal->setValueAtIndex(i,capLambdaHat2* (1/0.5));
                    if(_grilleSup)
                    {
                        itemOut_zSum->setValueAtIndex(i,zSum);
                        itemOut_zEffSum->setValueAtIndex(i,zEffSum);
                        itemOut_deltaSum->setValueAtIndex(i,deltaSum);
                        itemOut_deltaEffSum->setValueAtIndex(i,deltaEffSum);
                        itemOut_deltaSquareSum->setValueAtIndex(i,deltaSquareSum);
                        itemOut_1_zleDelta_effSum->setValueAtIndex(i,itemIn_deltaH->getSommaEffectiveHits(i));
                    }
                }
            }

            //Visualisation
            itemOut_PADCF->computeMinMax();
            itemOut_PADMLEBias->computeMinMax();
            itemOut_PADMLESimple->computeMinMax();
            itemOut_PADBLEqual->computeMinMax();
            itemOut_PADBLBasic->computeMinMax();
            itemOut_PADBLUnequal->computeMinMax();
            itemOut_zSum->computeMinMax();
            itemOut_zEffSum->computeMinMax();
            itemOut_deltaSum->computeMinMax();
            itemOut_deltaEffSum->computeMinMax();
            itemOut_deltaSquareSum->computeMinMax();
            itemOut_1_zleDelta_effSum->computeMinMax();
        }
    }
}

void LVOX3_StepComputePAD::computeStatistique(int i)
{

}

//Préléminaires pour futures vérifications (Présentement non utilisé)
void LVOX3_StepComputePAD::computeVariance()
{
    //_estimatedAttenuationCoefficient = (_I/_effectiveFreePath) - (_effectiveFreePath_Hits/(_N*pow(_effectiveFreePath,2)));
    //_estimatedVariance = (_I/_effectiveFreePath_Hits)*pow(1-(_effectiveFreePath_Hits/(_N*_I*_effectiveFreePath)),2); // Equation Non-fini
}

//Préléminaires pour futures vérifications (Présentement non utilisé)
void LVOX3_StepComputePAD::computeIntervals()
{
    /*
    double confidanceIntervalTop = 0;
    double confidanceIntervalBot = 0;

    bool ok = false;
    double gaussianDistribution = 1 - (_confidenceLevel.toDouble(&ok))/2;

    if(_estimatedAttenuationCoefficient*_deltaMean <= 0.5)
    {
        double IC = (_I + pow(gaussianDistribution,2)/(2*_N)) / (1+pow(gaussianDistribution,2)/(_N));
        double numShotsC = _N + pow(gaussianDistribution,2);
        double attenuationCoefficientC = (IC/_effectiveFreePath) - (_effectiveFreePath_Hits/(numShotsC*pow(_effectiveFreePath,2)));
        double estimatedVarianceC = (IC/_effectiveFreePath_Hits)*pow(1-(_effectiveFreePath_Hits/(numShotsC*IC*_effectiveFreePath)),2); // Equation Non-fini

        confidanceIntervalBot = attenuationCoefficientC -(gaussianDistribution * sqrt(estimatedVarianceC));
        confidanceIntervalTop = attenuationCoefficientC +(gaussianDistribution * sqrt(estimatedVarianceC));
    }
    else
    {
        confidanceIntervalBot = _estimatedAttenuationCoefficient -(gaussianDistribution * sqrt(_estimatedVariance));
        confidanceIntervalTop = _estimatedAttenuationCoefficient +(gaussianDistribution * sqrt(_estimatedVariance));
    }
    */
}

//Préléminaires pour futures vérifications (Présentement non utilisé)
void LVOX3_StepComputePAD::computeBias()
{
    //double L1 = _lambda1 * deltaMean;

    //Table 2

    //Bias <1%
    bool Bias1 = true;
    //1
    /*
    if(L1 <= 0.01 && _numShots >=100){}else{Bias1 = false;}
    */
    //2
    /*
    if(L1 <= 0.01 && _N >= 3){}else{Bias1 = false;}
    if(L1 <= 0.1 && _N >= 5){}else{Bias1 = false;}
    if(L1 <= 0.2 && _N >= 15){}else{Bias1 = false;}
    if(L1 <= 0.3 && _N >= 30){}else{Bias1 = false;}
    /*
    //3
    /*
    if(L1 <= 2 && L1<=0.01 && _numShots >= 100){}else{Bias1 = false;}
    */
    //4
    /*
    if(L1 <= 0.5 && L1 <=0.2 && _numShots >= 7){}else{Bias1 = false;}
    if(L1 <= 1 && L1 <= 0.2 && _numShots >= 10){}else{Bias1 = false;}
    if(L1 <= 1.5 && L1 <=0.2 && _numShots >= 15){}else{Bias1 = false;}
    if(L1 <= 2 && L1 <= 0.05 && _numShots >= 40){}else{Bias1 = false;}
    if(L1 <= 2.5 && L1 <= 0.005 && _numShots >= 75){}else{Bias1 = false;}
    if(L1 <= 3 && L1 <= 0.001 && _numShots >= 75){}else{Bias1 = false;}
    */

    //Bias <5%
    bool Bias5 = true;
    //1
    /*
    if(L1 <= 0.01 && _numShots >=20){}else{Bias5 = false;}
    if(L1 <= 0.05 && _numShots >=30){}else{Bias5 = false;}
    */
    //2
    /*
    if(L1 <= 0.05 && _N >= 3){}else{Bias5 = false;}
    if(L1 <= 0.1 && _N >= 5){}else{Bias5 = false;}
    if(L1 <= 0.3 && _N >= 10){}else{Bias5 = false;}
    */
    //3
    /*
    if(L1<=0.05 && _numShots >= 40){}else{Bias5 = false;}
    */
    //4
    /*
    if(L1 <= 1 && L1 <= 0.1 && _numShots >= 5){}else{Bias5 = false;}
    if(L1 <= 1.5 && L1 <=0.2 && _numShots >= 10){}else{Bias5 = false;}
    if(L1 <= 2 && L1 <= 0.2 && _numShots >= 15){}else{Bias5 = false;}
    if(L1 <= 2.5 && L1 <= 0.2 && _numShots >= 40){}else{Bias5 = false;}
    if(L1 <= 3 && L1 <= 0.1 && _numShots >= 75){}else{Bias5 = false;}
    */

    //Bias <10%
    bool Bias10 = true;
    //1
    /*
    if(L1 <= 0.01 && _numShots >=10){}else{Bias10 = false;}
    if(L1 <= 0.1 && _numShots >=20){}else{Bias10 = false;}
    */
    //2
    /*
    if(L1 <= 0.1 && _N >= 3){}else{Bias10 = false;}
    if(L1 <= 0.2 && _N >= 5){}else{Bias10 = false;}
    if(L1 <= 0.3 && _N >= 7){}else{Bias10 = false;}
    if(L1 <= 0.5 && _N >= 10){}else{Bias10 = false;}
    */
    //3
    /*
    if(L1 <= 2 && L1<=0.01 && _numShots >= 10){}else{Bias10 = false;}
    if(L1 <= 2.5 && L1<=0.1 && _numShots >= 20){}else{Bias10 = false;}
    if(L1 <= 3 && L1<=0.1 && _numShots >= 30){}else{Bias10 = false;}
    */
    //4
    /*
    if(L1 <= 1.5 && L1 <=0.2 && _numShots >= 5){}else{Bias10 = false;}
    if(L1 <= 2 && L1 <= 0.3 && _numShots >= 10){}else{Bias10 = false;}
    if(L1 <= 2.5 && L1 <=0.3 && _numShots >= 20){}else{Bias10 = false;}
    if(L1 <= 3 && L1 <= 0.3 && _numShots >= 40){}else{Bias10 = false;}
    */

}

bool LVOX3_StepComputePAD::isOccluded(int enteredShots, int thereticalShots){
    if(_isRelativeThreshold){
        return enteredShots*100.0/thereticalShots <= _NThreshold;
    }else{
        return enteredShots <= _NThreshold;
    }
}