From e9aaa0953d2113c241d7d0960393bba95fb9068b Mon Sep 17 00:00:00 2001 From: kradchen Date: Mon, 19 Jun 2023 09:33:38 +0800 Subject: [PATCH] Improve preprocessAscanBlock --- .../preprocessAScanBlockForReflection.cpp | 483 ++++++++++-------- 1 file changed, 272 insertions(+), 211 deletions(-) diff --git a/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp index ae69a25..1626c4b 100644 --- a/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp +++ b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -74,6 +75,12 @@ namespace Recon { } + void fakeFree(void*){} + Aurora::Matrix getPartMatrixColRef(const Aurora::Matrix & aMatrix,int col){ + std::shared_ptr ptr{aMatrix.getData()+aMatrix.getDimSize(0)*col*aMatrix.getValueType(),fakeFree}; + return Aurora::Matrix(ptr,{aMatrix.getDimSize(0),1,1},aMatrix.getValueType()); + } + Aurora::Matrix generateGaussWindow(size_t winLength) { auto n = Aurora::linspace(0,5,winLength); @@ -106,7 +113,7 @@ namespace Recon { std::vector performSignalProcessing( - const Aurora::Matrix &blockedAScans, + Aurora::Matrix &blockedAScans, const Aurora::Matrix &blockedSL, const Aurora::Matrix &blockedRL, const Aurora::Matrix &blockedSenderPosition, @@ -123,244 +130,298 @@ namespace Recon { Aurora::nantoval(blockedGain_t, 0); valid = valid*blockedGain_t; } - - auto _blockedAScans = blockedAScans/blockedGain; - if (reflectParams::removeDCOffset == 1) - { - _blockedAScans = _blockedAScans-Aurora::mean(_blockedAScans); - } int winLength = 100; - if (reflectParams::removeDCOffset == 1) + auto ascanMapValue = Aurora::zeros(1,blockedAScans.getDimSize(1),1); + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { - auto meanAS = mean(_blockedAScans.block(0,winLength-1,_blockedAScans.getDimSize(0) -1 -winLength)); - _blockedAScans = _blockedAScans-meanAS.getScalar(); - } - _blockedAScans.setBlockValue(0,0,winLength-1,0); - _blockedAScans.setBlockValue(0,_blockedAScans.getDimSize(0)-winLength-1,_blockedAScans.getDimSize(0)-1,0); - auto maxVal = max(abs(_blockedAScans)); - auto stdVal = Aurora::std(_blockedAScans); - auto meanVal = Aurora::mean(abs(_blockedAScans)); - auto ascanMapValue = meanVal*stdVal; - auto temp = (stdVal == 0) * (maxVal != 0); - Aurora::compareSet(valid,temp,1,0,Aurora::EQ); - if (reflectParams::findDefects==1) - { - valid = valid*(maxVal/meanVal >= reflectParams::epsilon); - } - // debugOutput.maxVal = maxVal; - // debugOutput.stdVal = stdVal; - // debugOutput.meanVal = meanVal; - // debugOutput.snr = maxVal ./ meanVal; - // debugOutput.snrPass = (maxVal./meanVal) < flags.dataSelection.epsilon; - if (reflectParams::suppressSameHead == 1) - { - for (size_t i = 0; i < blockedSL.getDataSize(); i++) + auto _blockedAScans = getPartMatrixColRef(blockedAScans,i); + _blockedAScans = _blockedAScans/blockedGain[i]; + if (reflectParams::removeDCOffset == 1) { - if (blockedSL[i] == blockedRL[i]) - { - auto begin_ptr = _blockedAScans.getData()+i*_blockedAScans.getDimSize(0); - int count = reflectParams::suppressSameHeadLength* sizeof(double); - std::memset(begin_ptr,0,count); - } + _blockedAScans = _blockedAScans-Aurora::mean(_blockedAScans).getScalar(); } - + + if (reflectParams::removeDCOffset == 1) + { + auto meanAS = mean(_blockedAScans.block(0,winLength-1,_blockedAScans.getDimSize(0) -1 -winLength)); + _blockedAScans = _blockedAScans-meanAS.getScalar(); + } + _blockedAScans.setBlockValue(0,0,winLength-1,0); + _blockedAScans.setBlockValue(0,_blockedAScans.getDimSize(0)-winLength-1,_blockedAScans.getDimSize(0)-1,0); + auto maxVal = max(abs(_blockedAScans)); + auto stdVal = Aurora::std(_blockedAScans); + auto meanVal = Aurora::mean(abs(_blockedAScans)); + ascanMapValue[i] = meanVal.getScalar()*stdVal.getScalar(); + auto temp = (stdVal == 0) * (maxVal != 0); + valid[i] = temp[0]==1?0:valid[i]; + // Aurora::compareSet(valid,temp,1,0,Aurora::EQ); + if (reflectParams::findDefects==1) + { + valid[i] = (valid[i]*(maxVal/meanVal >= reflectParams::epsilon)).getScalar(); + } + // debugOutput.maxVal = maxVal; + // debugOutput.stdVal = stdVal; + // debugOutput.meanVal = meanVal; + // debugOutput.snr = maxVal ./ meanVal; + // debugOutput.snrPass = (maxVal./meanVal) < flags.dataSelection.epsilon; + if (reflectParams::suppressSameHead == 1) + { + if (blockedSL[i] == blockedRL[i]) + { + auto begin_ptr = _blockedAScans.getData(); + int count = reflectParams::suppressSameHeadLength* sizeof(double); + std::memset(begin_ptr,0,count); + } + } + blockedAScans(Aurora::$,i)= _blockedAScans; } - auto fx = fft(_blockedAScans); - Aurora::Matrix fh; + auto fx_all = fft(blockedAScans); + if (reflectParams::useCorrelation && reflectParams::matchedFilterCeAScan) { - if (!preComputes.measuredCEUsed) + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { - fh = preComputes.matchedFilter.block(1, 0, nAScans-1); + Aurora::Matrix fx = getPartMatrixColRef(fx_all,i); + Aurora::Matrix fh; + if (!preComputes.measuredCEUsed) + { + fh = preComputes.matchedFilter.block(1, 0, nAScans-1); + } + else{ + fh = preComputes.matchedFilter(Aurora::$,blockedChannels[i]-1).toMatrix(); + } + double* value1 = Aurora::malloc(fx.getDataSize()); + vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1); + double* value2 = Aurora::malloc(fx.getDataSize()); + vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1); + double* realData = Aurora::malloc(fx.getDataSize()); + vdAdd(fx.getDataSize(), value1, value2, realData); + Aurora::Matrix real = Aurora::Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1)); + + vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1); + vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1); + double* imagData = Aurora::malloc(fx.getDataSize()); + vdSub(fx.getDataSize(), value1, value2, imagData); + Aurora::Matrix image = Aurora::Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1)); + + double* complexData = Aurora::malloc(real.getDataSize(), true); + cblas_dcopy(real.getDataSize(), real.getData(), 1 , complexData ,2); + cblas_dcopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2); + Aurora::Matrix complex = Aurora::Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex); + Aurora::free(value1); + Aurora::free(value2); + fx_all(Aurora::$,i) = complex; } - else{ - fh = preComputes.matchedFilter(Aurora::$,blockedChannels[0]-1).toMatrix(); - } - double* value1 = Aurora::malloc(fx.getDataSize()); - vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1); - double* value2 = Aurora::malloc(fx.getDataSize()); - vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1); - double* realData = Aurora::malloc(fx.getDataSize()); - vdAdd(fx.getDataSize(), value1, value2, realData); - Aurora::Matrix real = Aurora::Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1)); + blockedAScans = Aurora::real(ifft(fx_all)); - vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1); - vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1); - double* imagData = Aurora::malloc(fx.getDataSize()); - vdSub(fx.getDataSize(), value1, value2, imagData); - Aurora::Matrix image = Aurora::Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1)); - - double* complexData = Aurora::malloc(real.getDataSize(), true); - cblas_dcopy(real.getDataSize(), real.getData(), 1 , complexData ,2); - cblas_dcopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2); - Aurora::Matrix complex = Aurora::Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex); - Aurora::free(value1); - Aurora::free(value2); - _blockedAScans = Aurora::real(ifft(complex)); - - _blockedAScans.setBlockValue(0, 0, winLength-1, 0); - _blockedAScans.setBlockValue(0, _blockedAScans.getDimSize(0)-winLength, _blockedAScans.getDimSize(0)-1, 0); - _blockedAScans = fft(_blockedAScans); + blockedAScans.setBlockValue(0, 0, winLength-1, 0); + blockedAScans.setBlockValue(0, blockedAScans.getDimSize(0)-winLength, blockedAScans.getDimSize(0)-1, 0); + blockedAScans = fft(blockedAScans); } else{ - _blockedAScans = fx; + blockedAScans = fx_all; } auto exponent = offsetSignalFourier(preComputes, nSamples); // exponent.forceReshape(1, exponent.getDataSize(), 1); - _blockedAScans =_blockedAScans*exponent; - _blockedAScans = real(ifft(_blockedAScans)); + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) + { + Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); + blockedAScans(Aurora::$,i) =_blockedAScans*exponent; + } + blockedAScans = real(ifft(blockedAScans)); if (reflectParams::removeTransmissionSignal == 1) { - auto distanceEPosRPos = - sqrt( - ((blockedSenderPosition(0, Aurora::$).toMatrix() - blockedReceiverPosition(0, Aurora::$).toMatrix())^2) + - ((blockedSenderPosition(1, Aurora::$).toMatrix() - blockedReceiverPosition(1, Aurora::$).toMatrix())^2) + - ((blockedSenderPosition(2, Aurora::$).toMatrix() -blockedReceiverPosition(2, Aurora::$).toMatrix())^2)); - auto windowLength = reflectParams::windowLength; - auto start = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[1])*reflectParams::aScanReconstructionFrequency); - auto ende = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[0])*reflectParams::aScanReconstructionFrequency); - auto area = ceil(reflectParams::expectedPulseLength*info.Wavelength+1); - start = start- round(area*0.1); - ende = ende + round(area*0.9); - Aurora::compareSet(start, 1, 2, Aurora::NG); - Aurora::compareSet(start,nSamples,nSamples,Aurora::GT); - Aurora::compareSet(ende,nSamples,nSamples,Aurora::GT); - - for (size_t idx = 0; idx < nAScans; idx++) + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { - auto partData = _blockedAScans(Aurora::$,idx).toMatrix().block(0, start[idx]-1, ende[idx]-1); - - auto l_partData = partData.getDataSize(); - auto energyMove = _createDiffVector(partData,partData[0],0); - energyMove = _blockOp(energyMove,0); - energyMove = (energyMove^2)*0.5; - auto energyPot = partData^2; - - auto energySum = energyMove + energyPot; - - auto mean_energySum = Aurora::zeros(l_partData + windowLength,1); - for (size_t i = 0; i < windowLength; i++) + Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); + auto distanceEPosRPos = + sqrt( + ((blockedSenderPosition(0, i).toMatrix() - blockedReceiverPosition(0, i).toMatrix())^2) + + ((blockedSenderPosition(1, i).toMatrix() - blockedReceiverPosition(1, i).toMatrix())^2) + + ((blockedSenderPosition(2, i).toMatrix() -blockedReceiverPosition(2, i).toMatrix())^2)); + auto windowLength = reflectParams::windowLength; + auto start = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[1])*reflectParams::aScanReconstructionFrequency); + auto ende = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[0])*reflectParams::aScanReconstructionFrequency); + auto area = ceil(reflectParams::expectedPulseLength*info.Wavelength+1); + start = start- round(area*0.1); + ende = ende + round(area*0.9); + Aurora::compareSet(start, 1, 2, Aurora::NG); + Aurora::compareSet(start,nSamples,nSamples,Aurora::GT); + Aurora::compareSet(ende,nSamples,nSamples,Aurora::GT); + + // for (size_t idx = 0; idx < 1; idx++) { - double* begin = mean_energySum.getData() + i; - int length = l_partData; - vdAddI(length, energySum.getData(), 1, begin, 1, begin, 1); - } - mean_energySum = mean_energySum.block(0,windowLength/2-1,mean_energySum.getDataSize()-1-round((windowLength-0.001)/2)); - - auto g = _createDiffVector(mean_energySum,0,0); - g = _blockOp(g,0); - auto f = _createDiffVector(g,g[0],0); - f = _blockOp(f,0); - Aurora::compareSet(f, 0, 0, Aurora::LT); - Aurora::compareSet(f,g,0,0,Aurora::NG); - Aurora::compareSet(f,0,0,Aurora::NG); - auto mean_energySum_f = f; - long scheitelRow=0,scheitelCol=0; - max(g,Aurora::Column,scheitelRow,scheitelCol); - long indexRow=0, indexCol=0; - //不知道对不对？？？ - auto scheitel = scheitelRow+scheitelCol; - max(mean_energySum_f.block(0,0,scheitel),Aurora::Column,indexRow,indexCol); - size_t index = indexRow+indexCol; - double median_mean_energySum = Aurora::median(mean_energySum).getScalar(); - bool flag = false; - if (mean_energySum[index] > median_mean_energySum) - { - for (size_t i = index; i > 0 ; i--) + auto partData = _blockedAScans.block(0, start[0]-1, ende[0]-1); + + auto l_partData = partData.getDataSize(); + auto energyMove = _createDiffVector(partData,partData[0],0); + energyMove = _blockOp(energyMove,0); + energyMove = (energyMove^2)*0.5; + auto energyPot = partData^2; + + auto energySum = energyMove + energyPot; + + auto mean_energySum = Aurora::zeros(l_partData + windowLength,1); + for (size_t k = 0; k < windowLength; k++) { - if(mean_energySum[i] median_mean_energySum) + { + for (size_t k = index; k > 0 ; k--) + { + if(mean_energySum[k] hIndex; + for (size_t k = 0; k < hIndex_s.getDataSize(); k++) + { + if (hIndex_s[k]>0)hIndex.push_back(k); + } + if (hIndex.empty()){ + hIndex.push_back(mean_energySum.getDataSize()-scheitel-1); + } + int index2 = hIndex[0]+scheitel-1+10; + //???存在没有index？？？matlab中有这一分支 + auto sig_begin = index + start[0] ; + //在这里+1 调整会matlab的index便于后续计算的正确性 + auto sig_end = index2 + start[0] +1; + int winLength2 = 100; + winLength2 = sig_end+winLength2+200-sig_begin+1; + auto win = generateGaussWindow(winLength2); + double* begin = _blockedAScans.getData() + (size_t)sig_begin-1; + int length = winLength2; + auto winReverse = Aurora::zeros(1,win.getDataSize()); + std::reverse_copy(win.getData(), win.getData()+win.getDataSize(), winReverse.getData()); + vdMulI(length, begin, 1, winReverse.getData(), 1, begin, 1); + winLength2 = sig_end; + win = generateGaussWindow(winLength2); + begin = _blockedAScans.getData() ; + length = sig_end; + vdMulI(length, begin, 1, win.getData(), 1, begin, 1); + // _blockedAScans = fft(_blockedAScans); + // blockedAScans(Aurora::$,i) =_blockedAScans; } - auto hIndex_s =mean_energySum.block(0, scheitelCol, mean_energySum.getDataSize()-1)<(0.1*mean_energySum[scheitel]); - std::vector hIndex; - for (size_t i = 0; i < hIndex_s.getDataSize(); i++) - { - if (hIndex_s[i]>0)hIndex.push_back(i); - } - if (hIndex.empty()){ - hIndex.push_back(mean_energySum.getDataSize()-scheitel-1); - } - int index2 = hIndex[0]+scheitel-1+10; - //???存在没有index？？？matlab中有这一分支 - auto sig_begin = index + start[idx] ; - //在这里+1 调整会matlab的index便于后续计算的正确性 - auto sig_end = index2 + start[idx] +1; - - winLength = sig_end+winLength+200-sig_begin+1; - auto win = generateGaussWindow(winLength); - double* begin = _blockedAScans.getData() + (size_t)sig_begin-1; - int length = winLength; - auto winReverse = Aurora::zeros(1,win.getDataSize()); - std::reverse_copy(win.getData(), win.getData()+win.getDataSize(), winReverse.getData()); - vdMulI(length, begin, 1, winReverse.getData(), 1, begin, 1); - winLength = sig_end; - win = generateGaussWindow(winLength); - begin = _blockedAScans.getData() ; - length = sig_end; - vdMulI(length, begin, 1, win.getData(), 1, begin, 1); } } - _blockedAScans = fft(_blockedAScans); + // auto __blockedAScans = fft(blockedAScans(Aurora::$,0).toMatrix()); + blockedAScans = fft(blockedAScans); if(reflectParams::useOptPulse==1){ auto n = nSamples; auto nHalf = round((double)n/2); - _blockedAScans(0,Aurora::$) = std::complex{0,0}; - _blockedAScans.setBlockComplexValue(0,nHalf,_blockedAScans.getDimSize(0)-1, std::complex{0,0}); - _blockedAScans.setBlock(0,1,nHalf-1,_blockedAScans.block(0,1,nHalf-1)*2); - - - _blockedAScans = ifft(_blockedAScans, n); - - _blockedAScans = abs(_blockedAScans.block(0,0,nSamples-1)); - auto help = _blockedAScans.deepCopy(); - help(0,Aurora::$) = 0; - help.setBlock(0, 1, help.getDimSize(0)-1, _blockOp(_blockedAScans, 1)); + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { - auto help_bbegin = help.block(0, 0, help.getDimSize(0)-2); - auto help_bend = help.block(0, 1, help.getDimSize(0)-1); - help_bend = (help_bend>0)*(help_bbegin<0); - auto tempBlock = Aurora::zeros(_blockedAScans.getDimSize(0),_blockedAScans.getDimSize(1),_blockedAScans.getDimSize(2)); - tempBlock.setBlock(0, 0, tempBlock.getDimSize(0)-2, help_bend); - _blockedAScans = _blockedAScans*tempBlock; + Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); + + _blockedAScans(0,Aurora::$) = std::complex{0,0}; + _blockedAScans.setBlockComplexValue(0,nHalf,_blockedAScans.getDimSize(0)-1, std::complex{0,0}); + _blockedAScans.setBlock(0,1,nHalf-1,_blockedAScans.block(0,1,nHalf-1)*2); } - help = _blockedAScans.deepCopy(); - Aurora::compareSet(help,0,0,Aurora::LT); + + blockedAScans = ifft(blockedAScans, n); + + blockedAScans = abs(blockedAScans); + + auto help_all = blockedAScans.deepCopy(); + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) + { + Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); + Aurora::Matrix help = _blockedAScans.deepCopy(); + help(0,Aurora::$) = 0; + help.setBlock(0, 1, help.getDimSize(0)-1, _blockOp(_blockedAScans, 1)); + { + auto help_bbegin = help.block(0, 0, help.getDimSize(0)-2); + auto help_bend = help.block(0, 1, help.getDimSize(0)-1); + help_bend = (help_bend>0)*(help_bbegin<0); + auto tempBlock = Aurora::zeros(_blockedAScans.getDimSize(0),_blockedAScans.getDimSize(1),_blockedAScans.getDimSize(2)); + tempBlock.setBlock(0, 0, tempBlock.getDimSize(0)-2, help_bend); + blockedAScans(Aurora::$,i) = _blockedAScans*tempBlock; + } + help = _blockedAScans.deepCopy(); + Aurora::compareSet(help,0,0,Aurora::LT); - if(reflectParams::limitNumPulsesTo>0) - { - - if(size(help,1)>reflectParams::limitNumPulsesTo) + if(reflectParams::limitNumPulsesTo>0) { - auto help2 = Aurora::sort(help); - auto temp = repmat(help2(help2.getDimSize(0)-reflectParams::limitNumPulsesTo,Aurora::$).toMatrix(),help.getDimSize(0),1); - Aurora::compareSet(help,temp,0,Aurora::LT); - } - } + if(size(help,1)>reflectParams::limitNumPulsesTo) + { + auto help2 = Aurora::sort(help); + auto temp = repmat(help2(help2.getDimSize(0)-reflectParams::limitNumPulsesTo,Aurora::$).toMatrix(),help.getDimSize(0),1); + Aurora::compareSet(help,temp,0,Aurora::LT); + } + + } - if(reflectParams::normalizePeaks) - { - Aurora::compareSet(help,0,1,Aurora::GT); + if(reflectParams::normalizePeaks) + { + Aurora::compareSet(help,0,1,Aurora::GT); + } + help_all(Aurora::$,i) = help; } - - - - help= real(ifft(fft(help)*(preComputes.sincPeak_ft.block(1,0,nAScans-1)))); - Aurora::compareSet(_blockedAScans,0,help,Aurora::NL); + help_all = fft(help_all); + + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < help_all.getDimSize(1); i++) + { + Aurora::Matrix help = getPartMatrixColRef(help_all,i); + help_all(Aurora::$,i) = help*(preComputes.sincPeak_ft); + } + + help_all = real(ifft(help_all)); + + #pragma omp parallel for num_threads(32) + for (size_t i = 0; i < help_all.getDimSize(1); i++) + { + Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); + Aurora::Matrix help = getPartMatrixColRef(help_all,i); + Aurora::compareSet(_blockedAScans,0,help,Aurora::NL); + } } else{ - _blockedAScans = real(ifft(_blockedAScans)); + blockedAScans = real(ifft(blockedAScans)); } - result.push_back(_blockedAScans); + result.push_back(blockedAScans); result.push_back(valid); - ascanMapValue[0] = (float)ascanMapValue[0]; + for (size_t i = 0; i < ascanMapValue.getDataSize(); i++) + { + ascanMapValue[i] = (float)ascanMapValue[i]; + } + result.push_back(ascanMapValue); return result; } @@ -389,21 +450,21 @@ namespace Recon { std::cerr<<"error USCT II only!!"<