From 4913ab7e6092d2c12e6bec21e454327a18bebeda Mon Sep 17 00:00:00 2001
From: sunwen <sunwen@realcan.com>
Date: Wed, 14 Jun 2023 12:19:54 +0800
Subject: [PATCH] Update startReonstruction.

---
 src/startReconstructions.cpp                  | 130 +++++++++++-------
 .../reconstruction/reconstruction.h           |   1 -
 .../startTransmissionReconstruction.cpp       |  11 +-
 .../startTransmissionReconstruction.h         |  10 +-
 4 files changed, 95 insertions(+), 57 deletions(-)

diff --git a/src/startReconstructions.cpp b/src/startReconstructions.cpp
index a984dbd..b22a720 100644
--- a/src/startReconstructions.cpp
+++ b/src/startReconstructions.cpp
@@ -13,10 +13,14 @@
 #include "Parser/Data/MetaData.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;
@@ -26,8 +30,8 @@ void Recon::startReconstructions()
 {
     // std::string dataPath = "/home/AScans_Data/volunteer_20230328/20230328T123237/";
     // std::string refPath = "/home/AScans_Data/volunteer_20230328/20230328T123237/";
-    std::string dataPath = "/home/AScans_Data/ADW_TAS_Issue/20230512T135108/";
-    std::string refPath = "/home/AScans_Data/ADW_TAS_Issue/20230512T141626/";
+    std::string dataPath = "/home/AScans_Data/ADW_TAS_Issue/20230418T145123/";
+    std::string refPath = "/home/AScans_Data/ADW_TAS_Issue/20230418T141000/";
     std::string outputPath;
     Parser dataParser(dataPath);
     Parser refParser(refPath);
@@ -150,6 +154,9 @@ void Recon::startReconstructions()
         reflectParams::expectedAScanDataLength = ceil(expInfo.numberSamples * ((double)reflectParams::aScanReconstructionFrequency / expInfo.sampleRate));
     }
 
+    TransmissionReconstructionResult transmissionResult;
+    bool sosAvailable = false;
+    bool attAvailable = false;
     if(transParams::runTransmissionReco)
     {
         if(ceRef.measuredCEUsed)
@@ -160,56 +167,79 @@ void Recon::startReconstructions()
         {
             preComputes.matchedFilterRef = createMatchedFilter(ceRef.ceRef, ceRef.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
         }
+
+        preComputes.timeInterval = (double)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);
+        double* 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);
+
+        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;
     }
 
-
-
-    preComputes.timeInterval = (double)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
+    if(reflectParams::runReflectionReco)
     {
-        ceEstimatePulseLength.ce = ce.ceRef;
-        ceEstimatePulseLengthRef.ce = ceRef.ceRef;
-        ceEstimatePulseLength.ce_sf = reflectParams::aScanReconstructionFrequency;
-        ceEstimatePulseLengthRef.ce_sf = reflectParams::aScanReconstructionFrequency;
+        Matrix recoSOS = transmissionResult.recoSOS;
+        Matrix recoATT = transmissionResult.recoATT;
+        precalcImageParameters(geom);      
+        //检测可使用内存是否足够,输出警报用,todo
+        //checkEnvAndMemory(reflectParams.imageInfos.IMAGE_XYZ);
+        auto preProcessData = preprocessTransmissionReconstructionForReflection(recoSOS, recoATT, transmissionResult.ddmis, geom, temp);
+        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 = (double)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;
+        startReflectionReconstruction(&dataParser, preProcessData.saftMode, mp_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, geom, preProcessData.transRecos, expInfo, preComputes);
     }
-
-    transParams::pulseLengthSamples = estimatePulseLength(ceEstimatePulseLength, reflectParams::aScanReconstructionFrequency);
-    transParams::pulseLengthRefSamples = estimatePulseLength(ceEstimatePulseLengthRef, reflectParams::aScanReconstructionFrequency);
-
-    Matrix iMp;
-    Matrix mp_inter = intersect(motorPosAvailable, transParams::motorPos, iMp);
-    double* 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);
-    
-    startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, temp, tempRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
 }
\ No newline at end of file
diff --git a/src/transmissionReconstruction/reconstruction/reconstruction.h b/src/transmissionReconstruction/reconstruction/reconstruction.h
index 54fef9e..b57b509 100644
--- a/src/transmissionReconstruction/reconstruction/reconstruction.h
+++ b/src/transmissionReconstruction/reconstruction/reconstruction.h
@@ -12,7 +12,6 @@ struct DiscretizePositionValues{
 struct ArtResult{
     Aurora::Matrix outSOS;
     Aurora::Matrix outATT;
-    Aurora::Matrix outAll;
     bool bentReco;
     Aurora::Matrix allHitMaps;
 };
diff --git a/src/transmissionReconstruction/startTransmissionReconstruction.cpp b/src/transmissionReconstruction/startTransmissionReconstruction.cpp
index 987ffe1..1bc9f48 100644
--- a/src/transmissionReconstruction/startTransmissionReconstruction.cpp
+++ b/src/transmissionReconstruction/startTransmissionReconstruction.cpp
@@ -16,7 +16,7 @@
 using namespace Aurora;
 using namespace Recon;
 
-void Recon::startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
+TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
                                      const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
                                      const TempInfo& aTemp, const TempInfo& aTempRef, Recon::GeometryInfo& aGeom,
                                      GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
@@ -40,8 +40,9 @@ void Recon::startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, con
     Matrix reveiverList = removeDataFromArrays(positionValues.receiverCoordList, valid);
     auto transmissionReon = reconstructArt(tofData, attData, positionValues.dims, senderList, reveiverList, positionValues.res, aTemp.expectedSOSWater[0]);
 
-    MatlabWriter w1("/home/sun/transmissionATT.mat");
-    MatlabWriter w2("/home/sun/transmissionSOS.mat");
-    w1.write(transmissionReon.outATT, "ATT");
-    w2.write(transmissionReon.outSOS, "SOS");
+    TransmissionReconstructionResult result;
+    result.recoATT = transmissionReon.outATT;
+    result.recoSOS = transmissionReon.outSOS;
+    result.ddmis = positionValues.ddims;
+    return result;
 }                                    
\ No newline at end of file
diff --git a/src/transmissionReconstruction/startTransmissionReconstruction.h b/src/transmissionReconstruction/startTransmissionReconstruction.h
index aed0db1..57c8c66 100644
--- a/src/transmissionReconstruction/startTransmissionReconstruction.h
+++ b/src/transmissionReconstruction/startTransmissionReconstruction.h
@@ -4,12 +4,20 @@
 #include "Matrix.h"
 #include "../common/getMeasurementMetaData.h"
 #include "../common/getGeometryInfo.h"
+#include "transmissionReconstruction/detection/getTransmissionData.h"
 
 class Parser;
 
 namespace Recon
 {
-    void startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
+    struct TransmissionReconstructionResult
+    {
+        Aurora::Matrix recoSOS;
+        Aurora::Matrix recoATT;
+        Aurora::Matrix ddmis;
+    };
+
+    TransmissionReconstructionResult startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
                                          const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList, 
                                          const TempInfo& aTemp, const TempInfo& aTempRef, Recon::GeometryInfo& aGeom,
                                          GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,