Merging in Oscillate and Gridding steps (#14)

* Committing baseline

* Adding low rank penalty  and executable

* Adding LRobj object

* Adding baseline oscillateRecon.cpp code prior to gridding merge

* Fixing CMakeLists.txt for real this time.

* Checking in code from oscillate run

* Turing on GPU

* Adding some optimizations to pcSENSE

* Update README.md

* Update and rename CHANGELOG to CHANGELOG.md 

Adding gridding and Oscillate info.

Co-authored-by: Curtis <clj@localhost.localdomain>

Author: acerjanic

.devcontainer/devcontainer.json

@@ -0,0 +1,44 @@
+// For format details, see https://aka.ms/vscode-remote/devcontainer.json or the definition README at
+// https://github.com/microsoft/vscode-dev-containers/tree/master/containers/docker-existing-dockerfile
+{
+	"name": "PGDevEnv-PGI",
+
+	// Sets the run context to one level up instead of the .devcontainer folder.
+	"context": "..",
+
+	// Update the 'dockerFile' property if you aren't using the standard 'Dockerfile' filename.
+	"dockerFile": "../docker/powergrid-dev/Dockerfile",
+
+	// The optional 'runArgs' property can be used to specify additional runtime arguments.
+	"runArgs": [
+		"--gpus all"
+		// Uncomment the next line to use Docker from inside the container. See https://aka.ms/vscode-remote/samples/docker-in-docker for details.
+		// "-v","/var/run/docker.sock:/var/run/docker.sock",
+
+		// Uncomment the next line if you will be using a ptrace-based debugger like C++, Go, and Rust.
+		// "--cap-add=SYS_PTRACE", "--security-opt", "seccomp=unconfined"
+
+		// You may want to add a non-root user to your Dockerfile. On Linux, this will prevent
+		// new files getting created as root. See https://aka.ms/vscode-remote/containers/non-root-user
+		// for the needed Dockerfile updates and then uncomment the next line.
+		// "-u", "vscode"
+	],
+
+	// Use 'settings' to set *default* container specific settings.json values on container create. 
+	// You can edit these settings after create using File > Preferences > Settings > Remote.
+	"settings": { 
+		// This will ignore your local shell user setting for Linux since shells like zsh are typically 
+		// not in base container images. You can also update this to an specific shell to ensure VS Code 
+		// uses the right one for terminals and tasks. For example, /bin/bash (or /bin/ash for Alpine).
+		"terminal.integrated.shell.linux": null
+	},
+
+	// Uncomment the next line if you want to publish any ports.
+	// "appPort": [],
+
+	// Uncomment the next line to run commands after the container is created - for example installing git.
+	// "postCreateCommand": "apt-get update && apt-get install -y git",
+
+	// Add the IDs of extensions you want installed when the container is created in the array below.
+	"extensions": []
+}

.gitignore

@@ -31,6 +31,8 @@ build/*
 *.exe
 *.out
 *.app
+data/*
+
 
 # Archive files
 *.zip

ACKNOWLEDGEMENTS

@@ -0,0 +1,19 @@
+# Change Log
+All notable changes to this project will be documented in this file.
+
+## [Unreleased]
+
+## [v1.1.0] 2020-03-18
+## Updating with direct gridding reconstruction and low rank reconstructions
+* Implements calculation of density compensation functions in 2D/3D using method of Pipe et al.
+* Implements direct gridding reconstruction of coil images and sum-of-squares (SoS) reconstruction for fully sampled data
+* Implements phase corrected low rank reconstruction for Oscillate MRE
+    G McIlvain, AM Cerjanic, AG Christodoulou, MDJ McGarry, CL Johnson, “OSCILLATE: A Low-Rank Approach for Accelerated Magnetic Resonance Elastography,” to be presented at the 28th Annual Meeting of the International Society for Magnetic Resonance in Medicine, Sydney, Australia, April 18-23, 2020.
+
+## [v1.0.1] 2020-03-18
+## Updating with tagged docker hub builds
+* Implements pcSENSE reconstructions with both DFT and NUFFT based field correction
+* Implements 2D/3D field corrected non-Cartesian reconstructions
+
+## [v0.1.0] 2016-05-08
+### Initial release

CHANGELOG

@@ -98,7 +98,8 @@ set(PG_FILES
 	${PROJECT_SOURCE_DIR}/PowerGrid/pcSenseTimeSeg.cpp
 	${PROJECT_SOURCE_DIR}/PowerGrid/reconSolve.cpp
 	${PROJECT_SOURCE_DIR}/PowerGrid/acqTracking.cpp
-	${PROJECT_SOURCE_DIR}/PowerGrid/acqTracking.h
+	#${PROJECT_SOURCE_DIR}/PowerGrid/acqTracking.h
+	${PROJECT_SOURCE_DIR}/PowerGrid/reconSupport/directRecon.cpp
 )
 
 add_library(PowerGrid SHARED ${PG_FILES})
@@ -122,7 +123,7 @@ else()
 	add_definitions(-DBOOST_NO_ALIGNMENT)
 		if(OPENACC_GPU)
 			add_definitions(-DOPENACC_GPU)
-			set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Mnollvm -Mcuda -Minfo=accel -fast -acc -ta=tesla:cc35,cc60,cuda${CUDAVER},flushz -pgf90libs")
+			set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mp -tp=px -Mnollvm -Mcuda -Minfo=accel -fast -acc -ta=tesla:cc35,cc60,cuda${CUDAVER},flushz -pgf90libs")
 		elseif(OPENACC_MP)
 			add_definitions(-DOPENACC_MP)
 			set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Mnollvm -O1 -Minfo=accel -acc  -ta=multicore -pgf90libs")
@@ -133,9 +134,9 @@ else()
     endif()
 endif()
 
-#add_definitions(-DARMA_NO_DEBUG)
+add_definitions(-DARMA_NO_DEBUG)
 add_definitions(-DARMA_DONT_USE_HDF5)
-add_definitions(-DARMA_DONT_USE_OPENMP)
+#add_definitions(-DARMA_DONT_USE_OPENMP)
 if(OPENACC_GPU)
 	add_definitions(-DUSE_NVTX)
 endif()
@@ -158,11 +159,14 @@ set(POWERGRID_PCSENSE_TIMESEG
 add_executable(PowerGridPcSenseTimeSeg ${POWERGRID_PCSENSE_TIMESEG})
 target_link_libraries(PowerGridPcSenseTimeSeg ${LIBS} PowerGrid)
 
+add_subdirectory(PowerGrid/lowRank)
+add_subdirectory(PowerGrid/reconSupport)
+
 if (MPISupport)
     add_subdirectory(PowerGrid/MPI)
 endif (MPISupport)
 
 # Manage installing targets to sane locations in /opt
 install(TARGETS PowerGridPcSense DESTINATION bin)
 install(TARGETS PowerGridIsmrmrd DESTINATION bin)
-install(TARGETS PowerGrid DESTINATION lib)
+install(TARGETS PowerGrid DESTINATION lib)

CHANGELOG.md

@@ -209,158 +209,133 @@ inline Col<complex<T1>> Gnufft<T1>::operator/(const Col<complex<T1>> &d) const {
   Col<CxT1> temp(reinterpret_cast<CxT1 *>(pGridData), n1, false, true);
   return temp; // Return a vector of type T1
 }
-/*
+
 template <typename T1>
-Col<complex<T1>> Gnufft<T1>::trimmedForwardOp(
-    const Col<complex<T1>> &d,
-    const Col<complex<T1>> &tempInterp) const // Don't change these arguments
-{
-  // Let's trim the operations to avoid data overhead and transfers
-  // Basically if we know that the data points are zero, they have no impact
-  // on the transform
+inline Col<complex<T1>> Gnufft<T1>::forwardSpatialInterp(const Col<complex<T1>> &d) const {
   uword dataLength = this->n2;
-  uvec dataMaskTrimmed = find(abs(tempInterp) > 0);
-
-  Col<T1> kxTrimmed = kx.elem(dataMaskTrimmed);
-  Col<T1> kyTrimmed = ky.elem(dataMaskTrimmed);
 
-  Col<T1> kzTrimmed = kz.elem(dataMaskTrimmed);
-  uword dataLengthTrimmed = kxTrimmed.n_rows;
-  // std::cout << "Length of DataMaskTrimmed = " << dataLengthTrimmed <<
-  // std::endl;
-  // cout << "Separating real and imaginary data." << endl;
-
-  Col<T1> realData = real(d).eval();
-  Col<T1> imagData = imag(d).eval();
-
-  // Now we grab the data out of armadillo with the memptr() function
-  // This returns a pointer of the type of the elements of the
-  // array/vector/matrix/cube (3d matrix)
-  // Armadillo uses column major like MATLAB and Fortran, but different from
-  // 2D C++ arrays which are row major.
-
-  T1 *realDataPtr = realData.memptr();
-  T1 *imagDataPtr = imagData.memptr();
-
-  // cout << "Allocating memory for transformed data." << endl;
-  Col<T1> realXformedDataTrimmed(dataLengthTrimmed);
-  Col<T1> imagXformedDataTrimmed(dataLengthTrimmed);
+  const T1 *dataPtr = reinterpret_cast<const T1 *>(d.memptr());
 
-  // realXformedData.zeros();
-  // imagXformedData.zeros();
-  // cout << "Grabbing pointers for new memory." << endl;
+  parameters<T1> params;
+  params.sync = 0;
+  params.binsize = 128;
+  params.useLUT = 1;
+  params.kernelWidth = kernelWidth;
+  params.gridOS = gridOS;
+  params.imageSize[0] = Nx; // gridSize is gridOS times larger than imageSize.
+  params.imageSize[1] = Ny;
+  params.imageSize[2] = Nz;
+  params.gridSize[0] = std::ceil(gridOS * (T1)Nx);
+  params.gridSize[1] = std::ceil(gridOS * (T1)Ny);
+  if (params.gridSize[0] % 2) // 3D case, gridOS is adjusted on the z dimension:
+    params.gridSize[0] += 1; // That why we need to make sure here that the xy
+  if (params.gridSize[1] % 2) // dimensions have even sizes.
+    params.gridSize[1] += 1;
+  params.gridSize[2] = (Nz == 1) ? Nz : (std::ceil(gridOS * (T1)Nz)); // 2D or 3D
+  params.numSamples = dataLength;
+
+  unsigned int n = params.numSamples;
+  
+  memcpy(pGridData_os, dataPtr, sizeof(T1) * 2 * gridNumElems);
+  #pragma acc update device(pGridData_os[0:2 * gridNumElems])
 
-  T1 *realXformedDataPtr = realXformedDataTrimmed.memptr();
-  T1 *imagXformedDataPtr = imagXformedDataTrimmed.memptr();
+  //#pragma acc parallel loop present(pSamples [0:2 * n])
+    for (int ii = 0; ii < 2 * n; ii++) {
+        pSamples[ii] = (T1)0.0;
+    }
+    #pragma acc update device(pSamples [0:2 * n])
+
+    if (Nz == 1) {
+      gridding_forward_2D(n, params, kx,
+                         ky, beta, pSamples, LUT,
+                          sizeLUT, pGridData_os);
+    } else {
+      gridding_forward_3D(n, params, kx,
+                      ky, kz, beta, pSamples, LUT,
+                      sizeLUT, pGridData_os);
+    }
 
-  // Process data here, like calling a brute force transform, dft...
-  // I assume you create the pointers to the arrays where the transformed data
-  // will be stored
-  // realXformedDataPtr and imagXformedDataPtr and they are of type float*
+  #pragma acc update host(pSamples [0:2 * n])
 
-  // cout << "About to call the forward gridding routine." << endl;
-  cufftHandle *nPlan = const_cast<cufftHandle *>(&plan);
-  computeFd_CPU_Grid<T1>(dataLengthTrimmed, kxTrimmed,
-                         kyTrimmed, kzTrimmed, realDataPtr,
-                         imagDataPtr, Nx, Ny, Nz, gridOS, realXformedDataPtr,
-                         imagXformedDataPtr, kernelWidth, beta, LUT, sizeLUT,
-                         stream, nPlan, pGridData, pGridData_d, pGridData_os,
-                         pGridData_os_d, pSamples);
+  Col<CxT1> temp(reinterpret_cast<CxT1 *>(pSamples), n, false, true);
+  return temp; // Return a vector of type T1
 
-  // To return data, we need to put our data back into Armadillo objects
-  // We are telling the object how long it is because it will copy the data
-  // back into managed memory
-  // realXformedData(realXformedDataPtr, dataLength);
-  // imagXformedData(imagXformedDataPtr, dataLength);
-
-  // We can free the realDataXformPtr and imagDataXformPtr at this point and
-  // Armadillo will manage armadillo object memory as things change size or go
-  // out of scope and need to be destroyed
-  Col<complex<T1>> XformedDataTrimmed(dataLengthTrimmed);
-  Col<complex<T1>> XformedData(dataLength);
-  XformedDataTrimmed.set_real(realXformedDataTrimmed);
-  XformedDataTrimmed.set_imag(imagXformedDataTrimmed);
-  XformedData.elem(dataMaskTrimmed) = XformedDataTrimmed;
-
-  return conv_to<Col<complex<T1>>>::from(
-      XformedData); // Return a vector of type T1
 }
 
-// Adjoint transform operation
 template <typename T1>
-Col<complex<T1>>
-Gnufft<T1>::trimmedAdjointOp(const Col<complex<T1>> &d,
-                             const Col<complex<T1>> &tempInterp) const {
-
-  // uword dataLength = n2;
-  // Let's trim the operations to avoid data overhead and transfers
-  // Basically if we know that the data points are zero, they have no impact
-  // on the transform
+inline Col<complex<T1>> Gnufft<T1>::adjointSpatialInterp(const Col<complex<T1>> &d) const {
 
   uword dataLength = this->n2;
-  uvec dataMaskTrimmed = find(abs(tempInterp) > 0);
-  uword dataLengthTrimmed = dataMaskTrimmed.n_rows;
-  // std::cout << "Length of DataMaskTrimmed = " << dataLengthTrimmed <<
-  // std::endl;
 
-  Col<complex<T1>> dTrimmed = d.elem(dataMaskTrimmed);
-  Col<T1> kxTrimmed = kx.elem(dataMaskTrimmed);
-  Col<T1> kyTrimmed = ky.elem(dataMaskTrimmed);
-  Col<T1> kzTrimmed = kz.elem(dataMaskTrimmed);
+  const T1 *dataPtr = reinterpret_cast<const T1 *>(d.memptr());
 
-  Col<T1> realData = real(dTrimmed).eval();
-  Col<T1> imagData = imag(dTrimmed).eval();
+  parameters<T1> params;
+  params.sync = 0;
+  params.binsize = 128;
+  params.useLUT = 1;
+  params.kernelWidth = kernelWidth;
+  params.gridOS = gridOS;
+  params.imageSize[0] = Nx; // gridSize is gridOS times larger than imageSize.
+  params.imageSize[1] = Ny;
+  params.imageSize[2] = Nz;
+  params.gridSize[0] = std::ceil(gridOS * (T1)Nx);
+  params.gridSize[1] = std::ceil(gridOS * (T1)Ny);
+  if (params.gridSize[0] % 2) // 3D case, gridOS is adjusted on the z dimension:
+    params.gridSize[0] += 1; // That why we need to make sure here that the xy
+  if (params.gridSize[1] % 2) // dimensions have even sizes.
+    params.gridSize[1] += 1;
+  params.gridSize[2] = (Nz == 1) ? Nz : (std::ceil(gridOS * (T1)Nz)); // 2D or 3D
+  params.numSamples = dataLength;
+
+  unsigned int n = params.numSamples;
+
+  ReconstructionSample<T1>* samples; // Input Data
+  // allocate samples
+  samples = (ReconstructionSample<T1>*)malloc(
+  params.numSamples * sizeof(ReconstructionSample<T1>));
+
+  if (samples == NULL) {
+    printf("ERROR: Unable to allocate memory for input data\n");
+    exit(1);
+  }
 
-  T1 *realDataPtr = realData.memptr();
-  T1 *imagDataPtr = imagData.memptr();
+  //
+  for (int i = 0; i < params.numSamples; i++) {
 
-  Col<T1> realXformedData(n1);
-  Col<T1> imagXformedData(n1);
+    samples[i].kX = kx[i];
+    samples[i].kY = ky[i];
+    samples[i].kZ = kz[i];
 
-  // realXformedData.zeros();
-  // imagXformedData.zeros();
+    samples[i].real = dataPtr[2 * i];
+    samples[i].imag = dataPtr[2 * i + 1];
 
-  T1 *realXformedDataPtr = realXformedData.memptr();
-  T1 *imagXformedDataPtr = imagXformedData.memptr();
-  // Process data here, like calling a brute force transform, dft...
-  // I assume you create the pointers to the arrays where the transformed data
-  // will be stored
-  // realXformedDataPtr and imagXformedDataPtr and they are of type float*
+    samples[i].sdc = (T1)1.0;
+    // samples[i].t = t[i];
+  }
 
-  // T2 gridOS = 2.0;
-  cufftHandle *nPlan = const_cast<cufftHandle *>(&plan);
-  computeFH_CPU_Grid<T1>(dataLengthTrimmed, kxTrimmed,
-                         kyTrimmed, kzTrimmed, realDataPtr,
-                         imagDataPtr, Nx, Ny, Nz, gridOS, realXformedDataPtr,
-                         imagXformedDataPtr, kernelWidth, beta, LUT, sizeLUT,
-                         stream, nPlan, pGridData, pGridData_d, pGridData_os,
-                         pGridData_os_d);
-
-  //iftCpu<T2>(realXformedDataPtr,imagXformedDataPtr,
-  //           realDataPtr, imagDataPtr, kx.memptr(),
-  //           ky.memptr(), kz.memptr(),
-  //           ix.memptr(), iy.memptr(), iz.memptr(),
-  //           FM.memptr(), t.memptr(),
-  //           this->n2, this->n1
-  //           );
-
-  // realXformedData(realXformedDataPtr, dataLength);
-  // imagXformedData(imagXformedDataPtr, dataLength);
-
-  // We can free the realDataXformPtr and imagDataXformPtr at this point and
-  // Armadillo will manage armadillo object memory as things change size or go
-  // out of scope and need to be destroyed
-
-  Col<complex<T1>> XformedData(n1);
-  XformedData.set_real(realXformedData);
-  XformedData.set_imag(imagXformedData);
-  // XformedData.elem(dataMaskTrimmed) = XformedDataTrimmed;
-  // savemat("/shared/mrfil-data/data/PowerGridTest/64_64_16_4coils/ggrid.mat","img",XformedData);
-
-  return conv_to<Col<complex<T1>>>::from(
-      XformedData); // Return a vector of type T1
+  #pragma acc enter data copyin(samples [0:n])
+
+  //#pragma acc parallel loop 
+    for (int i = 0; i < gridNumElems; i++) {
+        pGridData_os[2 * i] = (T1)0.0;
+        pGridData_os[2 * i + 1] = (T1)0.0;
+    }
+  
+  #pragma acc update device(pGridData_os[0:2*gridNumElems])
+
+  if (Nz == 1) {
+    gridding_adjoint_2D(n, params, beta, samples,
+                        LUT, sizeLUT, pGridData_os);
+  } else {
+    gridding_adjoint_3D(n, params, beta, samples,
+                        LUT, sizeLUT, pGridData_os);
+  }
+
+  #pragma acc update host(pGridData_os[0:2*gridNumElems])
+  Col<CxT1> temp(reinterpret_cast<CxT1 *>(pGridData_os), gridNumElems, false, true);
+  return temp; // Return a vector of type T1
 }
-*/
+
 // Explicit Instantiation
 template class Gnufft<float>;
 template class Gnufft<double>;