From 7ffdb48acf04e7d980ed701f09f2997609426b88 Mon Sep 17 00:00:00 2001 From: Alex Cerjanic Date: Sun, 21 Jul 2019 18:07:56 -0500 Subject: [PATCH] Adding inferred time seg option (#12) --- PowerGrid/MPI/PowerGridPcSenseMPI_TS.cpp | 31 ++++++++------ PowerGrid/PowerGridIsmrmrd.cpp | 51 +++++++++++++++--------- PowerGrid/PowerGridPcSenseTimeSeg.cpp | 51 ++++++++++-------------- 3 files changed, 72 insertions(+), 61 deletions(-) diff --git a/PowerGrid/MPI/PowerGridPcSenseMPI_TS.cpp b/PowerGrid/MPI/PowerGridPcSenseMPI_TS.cpp index eddee5c..62f063a 100755 --- a/PowerGrid/MPI/PowerGridPcSenseMPI_TS.cpp +++ b/PowerGrid/MPI/PowerGridPcSenseMPI_TS.cpp @@ -76,18 +76,6 @@ int main(int argc, char** argv) return 1; } - if (TimeSegmentationInterp.compare("hanning") == 0) { - type = 1; - } else if (TimeSegmentationInterp.compare("minmax") == 0) { - type = 2; - } else if (TimeSegmentationInterp.compare("histo") == 0) { - type = 3; - } else { - std::cout << "Did not recognize temporal interpolator selection. " << std::endl - << "Acceptable values are hanning or minmax." << std::endl; - return 1; - } - } catch (boost::program_options::error& e) { std::cerr << "Error: " << e.what() << std::endl; std::cout << desc << std::endl; @@ -264,6 +252,25 @@ int main(int argc, char** argv) getCompleteISMRMRDAcqData(d, acqTrack, NSlice, NRep, NAvg, NEcho, NPhase, data, kx, ky, kz, tvec); + + // Deal with the number of time segments + if(!vm.count("TimeSegments")) { + switch(type) { + case 1: + L = ceil((arma::max(tvec) - arma::min(tvec))/2E-3); + break; + case 2: + L = ceil((arma::max(tvec) - arma::min(tvec))/3E-3); + break; + case 3: + L = ceil((arma::max(tvec) - arma::min(tvec))/3E-3); + break; + default: + L = 0; + } + std::cout << "Info: Setting L = " << L << " by default." << std::endl; + } + PMap = getISMRMRDCompletePhaseMap(d, NSlice, NSet, NRep, NAvg, NPhase, NEcho, NSeg, (uword)(Nx * Ny * Nz)); SMap = getISMRMRDCompleteSENSEMap >(d, sen, NSlice, (uword)(Nx * Ny * Nz)); diff --git a/PowerGrid/PowerGridIsmrmrd.cpp b/PowerGrid/PowerGridIsmrmrd.cpp index 21b7704..56f1f44 100755 --- a/PowerGrid/PowerGridIsmrmrd.cpp +++ b/PowerGrid/PowerGridIsmrmrd.cpp @@ -72,17 +72,7 @@ int main(int argc, char **argv) { std::cout << desc << std::endl; return 1; } - /* - if(precisionString.compare("double") ==0) { - typedef double PGPrecision; - } else if(precisionString.compare("float") == 0) { - typedef float PGPrecision; - } else { - typedef double PGPrecision; - std::cout << "Did not recognize precision option. Defaulting to - double precision." << std::endl; - } - */ + if (FourierTrans.compare("DFT") == 0) { FtType = 2; @@ -90,16 +80,20 @@ int main(int argc, char **argv) { } else if (FourierTrans.compare("NUFFT") == 0) { FtType = 1; - if (TimeSegmentationInterp.compare("hanning") == 0) { + if(!vm.count("TimeSegmentationInterp")) { type = 1; - } else if (TimeSegmentationInterp.compare("minmax") == 0) { - type = 2; - } else if (TimeSegmentationInterp.compare("histo") == 0) { - type = 3; } else { - std::cout << "Did not recognize temporal interpolator selection. " << std::endl - << "Acceptable values are hanning or minmax." << std::endl; - return 1; + if (TimeSegmentationInterp.compare("hanning") == 0) { + type = 1; + } else if (TimeSegmentationInterp.compare("minmax") == 0) { + type = 2; + } else if (TimeSegmentationInterp.compare("histo") == 0) { + type = 3; + } else { + std::cout << "Did not recognize temporal interpolator selection. " << std::endl + << "Acceptable values are hanning or minmax." << std::endl; + return 1; + } } } else if (FourierTrans.compare("DFTGrads") == 0) { FtType = 3; @@ -220,6 +214,25 @@ int main(int argc, char **argv) { getCompleteISMRMRDAcqData(d, acqTrack, NSlice, NRep, NAvg, NEcho, NPhase, data, kx, ky, kz, tvec); + // Deal with the number of time segments + if(!vm.count("TimeSegments")) { + switch(type) { + case 1: + L = ceil((arma::max(tvec) - arma::min(tvec))/2E-3); + break; + case 2: + L = ceil((arma::max(tvec) - arma::min(tvec))/3E-3); + break; + case 3: + L = ceil((arma::max(tvec) - arma::min(tvec))/3E-3); + break; + default: + L = 0; + } + std::cout << "Info: Setting L = " << L << " by default." << std::endl; + } + + std::cout << "Number of elements in kx = " << kx.n_rows << std::endl; std::cout << "Number of elements in ky = " << ky.n_rows << std::endl; std::cout << "Number of elements in kz = " << kz.n_rows << std::endl; diff --git a/PowerGrid/PowerGridPcSenseTimeSeg.cpp b/PowerGrid/PowerGridPcSenseTimeSeg.cpp index 691fe2a..0d030eb 100644 --- a/PowerGrid/PowerGridPcSenseTimeSeg.cpp +++ b/PowerGrid/PowerGridPcSenseTimeSeg.cpp @@ -74,28 +74,6 @@ int main(int argc, char **argv) std::cout << desc << std::endl; return 1; } - /* - if(precisionString.compare("double") ==0) { - typedef double PGPrecision; - } else if(precisionString.compare("float") == 0) { - typedef float PGPrecision; - } else { - typedef double PGPrecision; - std::cout << "Did not recognize precision option. Defaulting to - double precision." << std::endl; - } - */ - if (TimeSegmentationInterp.compare("hanning") == 0) { - type = 1; - } else if (TimeSegmentationInterp.compare("minmax") == 0) { - type = 2; - } else if (TimeSegmentationInterp.compare("histo") == 0) { - type = 3; - } else { - std::cout << "Did not recognize temporal interpolator selection. " << std::endl - << "Acceptable values are hanning, minmax, or histo." << std::endl; - return 1; - } } catch (boost::program_options::error& e) { @@ -114,9 +92,6 @@ int main(int argc, char **argv) acqTracking* acqTrack; processISMRMRDInput(rawDataFilePath, d, hdr, FM, sen, acqTrack); - //savemat("testFM.mat", "FM", FM); - //savemat("testSen.mat", "sen", sen); - std::cout << "Number of elements in SENSE Map = " << sen.n_rows << std::endl; std::cout << "Number of elements in Field Map = " << FM.n_rows << std::endl; @@ -198,6 +173,26 @@ int main(int argc, char **argv) getCompleteISMRMRDAcqData(d, acqTrack, NSlice, NRep, NAvg, NEcho, NPhase, data, kx, ky, kz, tvec); + + // Deal with the number of time segments + if(!vm.count("TimeSegments")) { + switch(type) { + case 1: + L = ceil((arma::max(tvec) - arma::min(tvec))/2E-3); + break; + case 2: + L = ceil((arma::max(tvec) - arma::min(tvec))/3E-3); + break; + case 3: + L = ceil((arma::max(tvec) - arma::min(tvec))/3E-3); + break; + default: + L = 0; + } + std::cout << "Info: Setting L = " << L << " by default." << std::endl; + } + + PMap = getISMRMRDCompletePhaseMap(d, NSlice, NSet, NRep, NAvg, NPhase, NEcho, NSeg, (uword) (Nx*Ny*Nz)); SMap = getISMRMRDCompleteSENSEMap>(d, sen, NSlice, (uword) (Nx*Ny*Nz)); @@ -212,18 +207,14 @@ int main(int argc, char **argv) std::cout << "Number of columns in data = " << data.n_cols << std::endl; - //Gnufft A(kx.n_rows, (float) 2.0, Nx, Ny, Nz, kx, ky, kz, ix, - // iy, iz); - //Gdft A(kx.n_rows, Nx*Ny*Nz,kx,ky,kz,ix,iy,iz,FM,tvec); pcSenseTimeSeg S_DWI(kx, ky, kz, Nx, Ny, Nz, nc, tvec, L, type, SMap, FMap, 0-PMap); - //pcSENSE> Sg(A, sen, kx.n_rows, Nx*Ny*Nz, nc); QuadPenalty R(Nx, Ny, Nz, beta, dims2penalize); ImageTemp = reconSolve, QuadPenalty>(data, S_DWI, R, kx, ky, kz, Nx, Ny, Nz, tvec, NIter); - //writeISMRMRDImageData(d, ImageTemp, Nx, Ny, Nz); + writeNiftiMagPhsImage(filename,ImageTemp,Nx,Ny,Nz); } }