#include "startReconstructions.h"
#include "Data/TemperatureData.h"
#include "Function.h"
#include "MatlabWriter.h"
#include "config/config.h"
#include "log/log.h"
#include "common/DICOMExporter.h"
#include "common/getMeasurementMetaData.h"
#include "common/getGeometryInfo.h"
#include "common/estimatePulseLength.h"
#include "log/notify.h"
#include "transmissionReconstruction/startTransmissionReconstruction.h"

#include <string>
#include <iostream>
#include "Parser/Parser.h"
#include "Parser/Data/MetaData.h"
#include "Parser/ShotList/ShotList.h"

#include "Matrix.h"
#include "MatlabReader.h"
#include "Function1D.h"
#include "Function2D.h"
#include "src/common/ceMatchedFilterHandling.h"
#include "src/reflectionReconstruction/preprocessData/estimateOffset.h"
#include "src/reflectionReconstruction/preprocessData/precalcImageParameters.h"
#include "src/reflectionReconstruction/preprocessData/preprocessTransmissionReconstructionForReflection.h"
#include "src/reflectionReconstruction/startReflectionReconstruction.h"


using namespace Recon;
using namespace Aurora;

int Recon::startReconstructions(const std::string& aDataPath, const std::string& aDataRefPath, const std::string& aOutputPath)
{
    MatlabWriter writer(aOutputPath);
    Parser dataParser(aDataPath);
    Parser refParser(aDataRefPath);
    Recon::DICOMExporter exporter(dataParser.getPatientData());
    exporter.setExportBasePath(aOutputPath);
    if(!dataParser.getShotList()->isValid())
    {
        RECON_INFO("startReconstructions with {0}.",aDataPath);

        RECON_INFO("data path is invalid.");
        return -9;
    }

    if(!refParser.getShotList()->isValid())
    {
        RECON_INFO("empty water path is invalid.");
        return -10;
    }
    //init total used receiver/emitter/motorPos
    Matrix motorPosTotal, slList, snList, rlList, rnList;
    if(transParams::runTransmissionReco && reflectParams::runReflectionReco)
    {
        motorPosTotal = auroraUnion(reflectParams::motorPos, transParams::motorPos);
        slList = auroraUnion(reflectParams::senderTasList, transParams::senderTasList);
        snList = auroraUnion(reflectParams::senderElementList, transParams::senderElementList);
        rlList = auroraUnion(reflectParams::receiverTasList, transParams::receiverTasList);
        rnList = auroraUnion(reflectParams::receiverElementList, transParams::receiverElementList);
    }
    else if(reflectParams::runReflectionReco)
    {
        motorPosTotal = reflectParams::motorPos;
        slList = reflectParams::senderTasList;
        snList = reflectParams::senderElementList;
        rlList = reflectParams::receiverTasList;
        rnList = reflectParams::receiverElementList;
    }
    else if(transParams::runTransmissionReco)
    {
        motorPosTotal = transParams::motorPos;
        slList = transParams::senderTasList;
        snList = transParams::senderElementList;
        rlList = transParams::receiverTasList;
        rnList = transParams::receiverElementList;
    }
    else
    {
        return -12;
    }

    //getMeasurementMetaData
    float maxNumTAS = max(auroraUnion(slList, rlList)).getData()[0];
    MeasurementInfo expInfo = loadMeasurementInfos(&dataParser);
    TempInfo temp = getTemperatureInfo(&dataParser, maxNumTAS);
    CEInfo ce = getCEInfo(&dataParser, expInfo);
    TransFormInfo transformationInfo = getTransformationMatrix(&dataParser, motorPosTotal);
    Matrix transformationMatrices = transformationInfo.rotationMatrix;
    Matrix motorPosAvailable = transformationInfo.motorPos;

    Matrix motorPosAvailableRef;
    MeasurementInfo expInfoRef;
    TempInfo tempRef;
    CEInfo ceRef;
    Matrix transformationMatricesRef;
    //Recon::notifyProgress(1);
    if(transParams::runTransmissionReco)
    {
        expInfoRef = loadMeasurementInfos(&refParser);
        tempRef = getTemperatureInfo(&refParser, maxNumTAS);
        ceRef = getCEInfo(&refParser, expInfoRef);
        transformationInfo = getTransformationMatrix(&refParser, motorPosTotal);
        transformationMatricesRef = transformationInfo.rotationMatrix;
        motorPosAvailableRef = transformationInfo.motorPos;
        if(transformationMatricesRef.isNull())
        {
            Matrix motorPos1 = Matrix::fromRawData(new float[1] {1}, 1);
            transformationMatricesRef = getTransformationMatrix(&refParser, motorPos1).rotationMatrix;
        }
        else
        {
            Matrix mpRef_inter = intersect(motorPosAvailableRef, transParams::motorPos);
            //extractMatchingReferenceMotorPosition
            for(int i=0; i<motorPosAvailable.getDimSize(0); ++i)
            {
                if (sum(mpRef_inter == motorPosAvailable.getData()[i],FunctionDirection::All).getData()[0] == 0)
                {
                    motorPosAvailableRef.getData()[i] = max(mpRef_inter).getData()[0];
                    int a = 0;
                }
            }
        }
    }
    else
    {
        transformationMatricesRef = Matrix();
        motorPosAvailableRef = Matrix();
    }
    //Recon::notifyProgress(2);
    if(!ce.ce.isNull() && !ceRef.ce.isNull())
    {
        Matrix isEqual = (ce.ce == ceRef.ce);
        if (sum(isEqual, FunctionDirection::All).getData()[0] == isEqual.getDataSize())
        {
            std::cout<<"CEs are not equal."<<std::endl;
        }
    }
    else
    {
        std::cout<<"CEs is null matrix."<<std::endl;
    }

    GeometryInfo geom = getGeometryInfo(motorPosAvailable, transformationMatrices, rlList, rnList, slList, snList);
    PreComputes preComputes;
    if((reflectParams::matchedFilterCeAScan && reflectParams::runReflectionReco) || transParams::runTransmissionReco)
    {
        if(ce.measuredCEUsed)
        {
            preComputes.matchedFilter = createMatchedFilter(ce.ce, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
        }
        else
        {
            preComputes.matchedFilter = createMatchedFilter(ce.ceRef, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
        }
    }
    //Recon::notifyProgress(3);
    if(expInfo.sampleRate != reflectParams::aScanReconstructionFrequency)
    {
        reflectParams::expectedAScanDataLength = ceil(expInfo.numberSamples * ((float)reflectParams::aScanReconstructionFrequency / expInfo.sampleRate));
    }
    //Recon::notifyProgress(4);
    TransmissionReconstructionResult transmissionResult;
    bool sosAvailable = false;
    bool attAvailable = false;
    if(transParams::runTransmissionReco)
    {
        if(ceRef.measuredCEUsed)
        {
            preComputes.matchedFilterRef = createMatchedFilter(ceRef.ce, ceRef.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
        }
        else
        {
            preComputes.matchedFilterRef = createMatchedFilter(ceRef.ceRef, ceRef.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
        }
        preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
        preComputes.measuredCEUsed = ce.measuredCEUsed;
        preComputes.measuredCE_TASIndices = ce.tasIndices;
        preComputes.measuredCE_receiverIndices = ce.receiverIndices;
        preComputes.offset = estimateOffset(expInfo, ce, preComputes.matchedFilter);

        expInfo.matchedFilter = preComputes.matchedFilter;
        expInfoRef.matchedFilter = preComputes.matchedFilterRef;

        CeEstimatePulseLength ceEstimatePulseLength;
        CeEstimatePulseLength ceEstimatePulseLengthRef;
        if(ce.measuredCEUsed)
        {
            ceEstimatePulseLength.ce = mean(ce.ce, FunctionDirection::Row);
            ceEstimatePulseLengthRef.ce = mean(ce.ce, FunctionDirection::Row);
            ceEstimatePulseLength.ce_sf = reflectParams::aScanReconstructionFrequency;
            ceEstimatePulseLengthRef.ce_sf = reflectParams::aScanReconstructionFrequency;
        }
        else
        {
            ceEstimatePulseLength.ce = ce.ceRef;
            ceEstimatePulseLengthRef.ce = ceRef.ceRef;
            ceEstimatePulseLength.ce_sf = reflectParams::aScanReconstructionFrequency;
            ceEstimatePulseLengthRef.ce_sf = reflectParams::aScanReconstructionFrequency;
        }

        transParams::pulseLengthSamples = estimatePulseLength(ceEstimatePulseLength, reflectParams::aScanReconstructionFrequency);
        transParams::pulseLengthRefSamples = estimatePulseLength(ceEstimatePulseLengthRef, reflectParams::aScanReconstructionFrequency);
        Matrix iMp;
        Matrix mp_inter = intersect(motorPosAvailable, transParams::motorPos, iMp);
        float* mpRef_interData = Aurora::malloc(iMp.getDataSize());
        for(int i=0; i<iMp.getDataSize(); ++i)
        {
            mpRef_interData[i] = motorPosAvailableRef[iMp[i] - 1];
        }
        Matrix mpRef_inter = Matrix::New(mpRef_interData, iMp.getDataSize());
        Matrix slList_inter = intersect(slList, transParams::senderTasList);
        Matrix snList_inter = intersect(snList, transParams::senderElementList);
        Matrix rlList_inter = intersect(rlList, transParams::receiverTasList);
        Matrix rnList_inter = intersect(rnList, transParams::receiverElementList);
        transParams::aScanReconstructionFrequency = reflectParams::aScanReconstructionFrequency;
        transParams::gpuSelectionList = reconParams::gpuSelectionList;

        GeometryInfo geomRef = getGeometryInfo(motorPosAvailableRef, transformationMatricesRef, rlList, rnList, slList, snList);
        RECON_INFO("Start transmissionRecostruction.");
        //Recon::notifyProgress(5);
        transmissionResult = startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, temp, tempRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
        attAvailable = true;
        sosAvailable = true;
        exporter.exportDICOM(transmissionResult.recoSOS, DICOMExporter::SOS);
        exporter.exportDICOM(transmissionResult.recoATT, DICOMExporter::ATT);
        //  writer.setMatrix(transmissionResult.recoSOS, "sos");
        //  writer.setMatrix(transmissionResult.recoATT, "att");
    }

    if(reflectParams::runReflectionReco)
    {
        Matrix recoSOS = transmissionResult.recoSOS;
        Matrix recoATT = transmissionResult.recoATT;
        precalcImageParameters(geom);   
        //Recon::notifyProgress(21);   
        //检测可使用内存是否足够,输出警报用,todo
        //checkEnvAndMemory(reflectParams.imageInfos.IMAGE_XYZ);
        auto preProcessData = preprocessTransmissionReconstructionForReflection(recoSOS, recoATT, transmissionResult.ddmis, geom, temp);
        //Recon::notifyProgress(22);
        Matrix mp_inter = intersect(motorPosAvailable, reflectParams::motorPos);
        Matrix slList_inter = intersect(slList, reflectParams::senderTasList);
        Matrix snList_inter = intersect(snList, reflectParams::senderElementList);
        Matrix rlList_inter = intersect(rlList, reflectParams::receiverTasList);
        Matrix rnList_inter = intersect(rnList, reflectParams::receiverElementList);
        preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
        preComputes.measuredCEUsed = ce.measuredCEUsed;
        preComputes.measuredCE_TASIndices = ce.tasIndices;
        preComputes.measuredCE_receiverIndices = ce.receiverIndices;
        preComputes.offset = estimateOffset(expInfo, ce, preComputes.matchedFilter);

        reflectParams::gpuSelectionList = reconParams::gpuSelectionList;
        //Recon::notifyProgress(25);
        RECON_INFO("Start reflectionRecostruction.");
        Matrix env = startReflectionReconstruction(&dataParser, preProcessData.saftMode, mp_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, geom, preProcessData.transRecos, expInfo, preComputes);
         writer.setMatrix(env, "reflect");
        exporter.exportDICOM(env, Recon::DICOMExporter::REFL);
        // Recon::notifyProgress(99);
    }
        //  writer.write();
    return 0;

}