.travis.yml

@@ -53,7 +53,7 @@ matrix:
         - sudo apt-get update -qq
         - sudo apt-get install -qq liblapacke-dev liblapack-dev
 
-    - env: OMP=0
+    - env: OMP=0 CFLAGS=-fblocks LDFLAGS=-lBlocksRuntime
       os: linux
       compiler: clang
       sudo: required
@@ -67,7 +67,7 @@ matrix:
         - make utest
       before_install:
         - sudo apt-get update -qq
-        - sudo apt-get install -qq liblapacke-dev
+        - sudo apt-get install -qq liblapacke-dev libblocksruntime-dev
 
     - env: OMP=1
       os: linux
@@ -99,7 +99,7 @@ matrix:
         - sudo apt-get update -qq
         - sudo apt-get install -qq liblapack-dev
 
-    - env: CUDA=1 CUDA_BASE=/usr/local/cuda-7.5/
+    - env: CUDA=1 CUDA_BASE=/usr/local/cuda-7.5/ CFLAGS=-fblocks LDFLAGS=-lBlocksRuntime
       os: linux
       compiler: clang-3.5
       sudo: required
@@ -114,7 +114,7 @@ matrix:
         - sudo apt-get update -qq
         - sudo apt-get install -qq clang-3.5
         - sudo apt-get install -qq --no-install-recommends cuda-drivers cuda-core-7.5 cuda-cudart-dev-7.5 cuda-cufft-dev-7.5 cuda-cublas-dev-7.5
-        - sudo apt-get install -qq liblapacke-dev
+        - sudo apt-get install -qq liblapacke-dev libblocksruntime-dev
 
     - env: BUILD_NAME="CMake"
       os: linux
@@ -156,7 +156,6 @@ matrix:
 
     - env: BUILD_NAME="CMake + CUDA" CUDA_BASE=/usr/local/cuda-7.5/
       os: linux
-      compiler: clang-3.5
       sudo: required
       dist: trusty
       addons:
@@ -196,16 +195,6 @@ matrix:
         - CC=gcc CXX=g++ cmake -DLINALG_VENDOR=LAPACKE -DLAPACKE_DIR=/usr/lib -DUSE_MATLAB=OFF -DUSE_CUDA=ON -DCUDA_TOOLKIT_ROOT_DIR=${CUDA_BASE} ..
         - make
 
-    - env: MACPORTS=0
-      os: osx
-      compiler: gcc-4.8
-      script:
-        - make test
-        # make utest
-      before_install:
-        - brew update
-        - brew install fftw gcc48 homebrew/science/openblas
-
 script:
   - make bart
   - make all

CMakeLists.txt

@@ -11,7 +11,7 @@ enable_language(C)
 macro(use_c99)
   if (CMAKE_VERSION VERSION_LESS "3.1")
     if (CMAKE_C_COMPILER_ID STREQUAL "GNU")
-      set (CMAKE_C_FLAGS "--std=gnu99 ${CMAKE_C_FLAGS}")
+      set (CMAKE_C_FLAGS "--std=gnu11 ${CMAKE_C_FLAGS}")
     endif ()
   else ()
     set (CMAKE_C_STANDARD 99)
@@ -264,20 +264,20 @@ install(TARGETS bartsupport
         ARCHIVE DESTINATION lib/static)
 
 
-## PASS #1:  Build stand-alone programs
-## Generate stand alone programs
-file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/StandAloneCode)
-foreach(curr_prog ${ALLPROGS})
-  CONFIG_BY_REPLACEMENT( ${CMAKE_CURRENT_LIST_DIR}/src/main.c ${CMAKE_CURRENT_BINARY_DIR}/StandAloneCode/${curr_prog}.c
-                         "/* Generated by cmake */"
-                         "main_real" "main_${curr_prog}")
-  bart_add_executable(${curr_prog} ${CMAKE_CURRENT_BINARY_DIR}/StandAloneCode/${curr_prog}.c ${CMAKE_CURRENT_LIST_DIR}/src/${curr_prog}.c)
-  target_link_libraries(${curr_prog} bartsupport )
-  install(TARGETS ${curr_prog}
-        RUNTIME DESTINATION bin
-        LIBRARY DESTINATION lib
-        ARCHIVE DESTINATION lib/static)
-endforeach()
+### PASS #1:  Build stand-alone programs
+### Generate stand alone programs
+#file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/StandAloneCode)
+#foreach(curr_prog ${ALLPROGS})
+#  CONFIG_BY_REPLACEMENT( ${CMAKE_CURRENT_LIST_DIR}/src/main.c ${CMAKE_CURRENT_BINARY_DIR}/StandAloneCode/${curr_prog}.c
+#                         "/* Generated by cmake */"
+#                        "main_real" "main_${curr_prog}")
+#  bart_add_executable(${curr_prog} ${CMAKE_CURRENT_BINARY_DIR}/StandAloneCode/${curr_prog}.c ${CMAKE_CURRENT_LIST_DIR}/src/${curr_prog}.c)
+#  target_link_libraries(${curr_prog} bartsupport )
+#  install(TARGETS ${curr_prog}
+#        RUNTIME DESTINATION bin
+#        LIBRARY DESTINATION lib
+#        ARCHIVE DESTINATION lib/static)
+#endforeach()
 
 ## More crafty file manipulation so that we maintain backward comaptibility
 ## with the Makefile
@@ -327,7 +327,7 @@ install(TARGETS bart
 
 #==============================================
 # TODO: Matlab code
-option(USE_MATLAB "Specify if the optional matlab programs should be built" ON)
+option(USE_MATLAB "Specify if the optional matlab programs should be built" OFF)
 if(USE_MATLAB)
   find_package(Matlab REQUIRED)
   if(MATLAB_FOUND)

Makefile

@@ -24,7 +24,7 @@ FFTWTHREADS?=1
 ISMRMRD?=0
 
 DESTDIR ?= /
-PREFIX ?= usr/
+PREFIX ?= usr/local/
 
 BUILDTYPE = Linux
 UNAME = $(shell uname -s)
@@ -91,7 +91,7 @@ CFLAGS ?= $(OPT) -Wmissing-prototypes
 CXXFLAGS ?= $(OPT)
 
 ifeq ($(BUILDTYPE), MacOSX)
-	CC ?= gcc-mp-4.7
+	CC ?= gcc-mp-6
 else
 	CC ?= gcc
 	# for symbols in backtraces

README

@@ -49,6 +49,9 @@ Updates and further information can be found here:
 GCC compiler, the FFTW library, and optionally CUDA.
 (see recon/Makefile to turn options on or off)
 
+It should be possible to use the clang compiler, but older
+version before version 4.0 are known to cause problems.
+
 The software can be used in combination with Matlab or octave.
 
 
@@ -81,11 +84,11 @@ To install the required libraries on Debian and Ubuntu run:
 ### 2.1.2. Mac OS X
 
 Xcode is required and it is recommended to install a newer version 
-of gcc (4.7 seems to work) from MacPorts (http://www.macports.org/).
+of gcc from MacPorts (http://www.macports.org/).
 
 
     $ sudo port install fftw-3-single
-    $ sudo port install gcc47
+    $ sudo port install gcc6
     $ sudo port install libpng
     $ sudo port install openblas
 

doc/references.txt

@@ -12,11 +12,15 @@ Software Toolbox and Programming Library for Compressed Sensing and
 Parallel Imaging,
 ISMRM Workshop on Data Sampling and Image Reconstruction, Sedona 2013.
 
+Tamir JI, Ong F, Cheng JY, Uecker M, Lustig M. Generalized Magnetic
+Resonance Image Reconstruction using The Berkeley Advanced
+Reconstruction Toolbox,
+ISMRM Workshop on Data Sampling and Image Reconstruction, Sedona 2016.
 
 
 
 - sensitivity-encoded parallel imaging -
-(commands: itsense, pocsense, bpsense, rsense, pics)
+(commands: itsense, pocsense, bpsense, pics)
 
 
 Ra JB and Rim CY.
@@ -58,6 +62,11 @@ trajectories.
 Annual Meeting of the ISMRM, Glasgow 2001,
 In: Proc Intl Soc Mag Reson Med 9; 767.
 
+Uecker M, Zhang S, Frahm J.
+Nonlinear Inverse Reconstruction for Real-time MRI of the
+Human Heart Using Undersampled Radial FLASH.
+Magn Reson Med 2010; 63:1456-1462.
+
 Ong F, Uecker M, Jiang W, Lustig M.
 Fast Non-Cartesian Reconstruction with Pruned Fast Fourier Transform.
 Annual Meeting ISMRM, Toronto 2015,
@@ -119,6 +128,11 @@ Calibrationless Parallel Imaging Reconstruction Based on Structured
 Low-Rank Matrix Completion.
 Magn Reson Med 2014; 72:959-970.
 
+Holme HCM, Rosenzweig S, Ong F, Wilke RN, Lustig M, Uecker M.
+ENLIVE: An Efficient Nonlinear Method for Calibrationless and
+Robust Parallel Imaging.
+arXiv:1706.09780 [physics.med-ph]
+
 
 
 
@@ -169,7 +183,7 @@ Magn Reson Med 2009; 61:145-152.
 
 
 - sparsity transforms, variational penalties, regularization -
-(commands: cdf97, rof, lrmatrix, pocsense, rsense, pics)
+(commands: cdf97, rof, lrmatrix, pocsense, pics)
 
 
 Rudin LI, Osher S, Fatemi E.
@@ -177,7 +191,7 @@ Nonlinear total variation based noise removal algorithms,
 Physica D: Nonlinear Phenomena 1992; 60:259-268.
 
 Figueiredo MAT and Nowak RD.
-An EM algorithm for wavelet-based image restoration
+An EM algorithm for wavelet-based image restoration.
 IEEE Trans Image Process 2003; 12:906-916.
 
 Ong F, Uecker M, Tariq U, Hsiao A, Alley MT, Vasanawala SS, Lustig M.
@@ -186,8 +200,7 @@ Magn Reson Med 2015; 73:828-842.
 
 Ong F, Lustig M.
 Beyond low rank + sparse: Multi-scale low rank matrix decomposition,
-preprint 2015; arXiv:1507.08751.
-
+IEEE J Sel Topics Signal Process 2016; 10:672-687.
 
 
 
@@ -283,10 +296,25 @@ Radiology 2016; 278:245-256.
 Cheng JY, Hanneman K, Zhang T, Alley MT, Lai P, Tamir JI, Uecker M,
 Pauly JM, Lustig M, Vasanawala SS.
 Comprehensive Motion-Compensated Highly-Accelerated 4D Flow MRI with
-Ferumoxytol Enhancement for Pediatric Congenital Heart Disease,
-J Magn Reson Imaging, Epub (2015)
+Ferumoxytol Enhancement for Pediatric Congenital Heart Disease.
+J Magn Reson Imaging 2016; 43:1355-1368.
 
 Tamir JI, Uecker M, Chen W, Lai P, Aleey MT, Vasanawala SS, Lustig M.
-T2-Shuffling: Sharp, Multi-Contrast, Volumetric Fast Spin-Echo Imaging
-Magn Recon Med; in press (2015).
+T2-Shuffling: Sharp, Multi-Contrast, Volumetric Fast Spin-Echo Imaging.
+Magn Recon Med 2017; 77:180-195.
+
+Uecker M, Lustig M.
+Estimating Absolute-Phase Maps Using ESPIRiT and Virtual Conjugate Coils.
+Magn Reson Med 2017; 77:1201-1207.
+
+Moghari MH, Uecker M, Roujol S, Sabbagh M, Geva T, Powell AJ.
+Accelerated Whole-heart Magnetic Resonance Angiography Using
+a Variable-Density Poisson-Disc Undersampling Pattern and
+Compressed Sensing Reconstruction.
+Magn Reson Med 2018; 79:761-769.
+
+Cheng JY, Zhang T, Alley MT, Uecker M, Lustig M, Pauly JM, Vasanawala SS.
+Comprehensive Multi-Dimensional MRI for the Simultaneous Assessment
+of Cardiopulmonary Anatomy and Physiology.
+Scientific Reports 2017; 7:5330.
 

rules/misc.mk

@@ -19,4 +19,4 @@ DOTHIS := $(shell $(root)/rules/update-version.sh)
 $(srcdir)/misc/version.o: $(srcdir)/misc/version.inc
 
 
-UTARGETS += test_pattern
+UTARGETS += test_pattern test_types

rules/num.mk

@@ -18,5 +18,5 @@ lib/libnum.a: libnum.a($(numobjs))
 
 
 UTARGETS += test_multind test_flpmath test_splines test_linalg test_polynom test_window
-UTARGETS += test_blas
+UTARGETS += test_blas test_mdfft
 

scripts/espirit_econ.sh

+#!/bin/bash
+# Copyright 2018. Martin Uecker.
+# All rights reserved. Use of this source code is governed by
+# a BSD-style license which can be found in the LICENSE file.
+#
+# Authors:
+# 2018 Martin Uecker <martin.uecker@med.uni-goettingen.de>
+#
+# Memory-saving ESPIRiT
+#
+set -e
+
+LOGFILE=/dev/stdout
+
+title=$(cat <<- EOF
+	ESPIRiT-ECON
+EOF
+)
+
+helpstr=$(cat <<- EOF
+-l logfile
+-h help
+EOF
+)
+
+usage="Usage: $0 [-h] <kspace> <output>"
+
+echo "$title"
+echo
+
+while getopts "hl:" opt; do
+        case $opt in
+	h)
+		echo "$usage"
+		echo
+		echo "$helpstr"
+		exit 0
+	;;
+	l)
+		LOGFILE=$(readlink -f "$OPTARG")
+	;;
+        \?)
+        	echo "$usage" >&2
+		exit 1
+        ;;
+        esac
+done
+
+shift $((OPTIND - 1))
+
+
+if [ $# -lt 2 ] ; then
+
+        echo "$usage" >&2
+        exit 1
+fi
+
+
+export PATH=$TOOLBOX_PATH:$PATH
+
+input=$(readlink -f "$1")
+output=$(readlink -f "$2")
+
+if [ ! -e $input.cfl ] ; then
+	echo "Input file does not exist." >&2
+	echo "$usage" >&2
+	exit 1
+fi
+
+if [ ! -e $TOOLBOX_PATH/bart ] ; then
+        echo "\$TOOLBOX_PATH is not set correctly!" >&2
+	exit 1
+fi
+
+
+#WORKDIR=$(mktemp -d)
+# Mac: http://unix.stackexchange.com/questions/30091/fix-or-alternative-for-mktemp-in-os-x
+WORKDIR=`mktemp -d 2>/dev/null || mktemp -d -t 'mytmpdir'`
+trap 'rm -rf "$WORKDIR"' EXIT
+cd $WORKDIR
+
+# start group for redirection of output to the logfile
+{
+
+XX=$(bart show -d0 $input)
+YY=$(bart show -d1 $input)
+ZZ=$(bart show -d2 $input)
+
+DIM=2
+# To decouple along another dimension:
+# 1. change DIM
+# 2. replace ZZ below
+# 3. change the ecaltwo command
+
+bart ecalib -1 $input eon
+
+# zero-pad
+bart fft $(bart bitmask ${DIM}) eon eon_fft
+bart resize -c ${DIM} ${ZZ} eon_fft eon_fft2
+bart fft -i $(bart bitmask ${DIM}) eon_fft2 eon
+
+for i in `seq -w 0 $(($ZZ - 1))` ; do
+	
+	bart slice ${DIM} $i eon sl
+	bart ecaltwo ${XX} ${YY} 1 sl sens-$i.coo
+done
+
+# # join slices back together
+bart join ${DIM} sens-*.coo $output
+
+} > $LOGFILE
+
+exit 0

src/bart.c

@@ -60,11 +60,12 @@ int main_bart(int argc, char* argv[])
 			exit(1);
 		}
 
-		const char* tpath[3] = {
-
+		const char* tpath[] = {
+#ifdef TOOLBOX_PATH_OVERRIDE
 			getenv("TOOLBOX_PATH"),
-			"/usr/lib/bart/commands/",
+#endif
 			"/usr/local/lib/bart/commands/",
+			"/usr/lib/bart/commands/",
 		};
 
 		for (unsigned int i = 0; i < ARRAY_SIZE(tpath); i++) {

src/bench.c

 /* Copyright 2014. The Regents of the University of California.
- * Copyright 2015-2017. Martin Uecker.
+ * Copyright 2015-2018. Martin Uecker.
  * All rights reserved. Use of this source code is governed by
  * a BSD-style license which can be found in the LICENSE file.
  * 
  * Authors: 
- * 2014-2017 Martin Uecker <martin.uecker@med.uni-goettingen.de>
+ * 2014-2018 Martin Uecker <martin.uecker@med.uni-goettingen.de>
  * 2014 Jonathan Tamir <jtamir@eecs.berkeley.edu>
  */
 
@@ -19,6 +19,8 @@
 #include "num/rand.h"
 #include "num/init.h"
 #include "num/ops.h"
+#include "num/mdfft.h"
+#include "num/fft.h"
 
 #include "wavelet/wavthresh.h"
 
@@ -414,6 +416,94 @@ static double bench_wavelet(long scale)
 }
 
 
+static double bench_generic_mdfft(long dims[DIMS], unsigned long flags)
+{
+	complex float* x = md_alloc(DIMS, dims, CFL_SIZE);
+	complex float* y = md_alloc(DIMS, dims, CFL_SIZE);
+
+	md_gaussian_rand(DIMS, dims, x);
+
+	double tic = timestamp();
+
+	md_fft(DIMS, dims, flags, 0u, y, x);
+
+	double toc = timestamp();
+
+	md_free(x);
+	md_free(y);
+
+	return toc - tic;
+}
+
+static double bench_mdfft(long scale)
+{
+	long dims[DIMS] = { 1, 128 * scale, 128 * scale, 1, 1, 4, 1, 4 };
+	return bench_generic_mdfft(dims, 6ul);
+}
+
+
+
+static double bench_generic_fft(long dims[DIMS], unsigned long flags)
+{
+	complex float* x = md_alloc(DIMS, dims, CFL_SIZE);
+	complex float* y = md_alloc(DIMS, dims, CFL_SIZE);
+
+	md_gaussian_rand(DIMS, dims, x);
+
+	double tic = timestamp();
+
+	fft(DIMS, dims, flags, y, x);
+
+	double toc = timestamp();
+
+	md_free(x);
+	md_free(y);
+
+	return toc - tic;
+}
+
+
+
+static double bench_fft(long scale)
+{
+	long dims[DIMS] = { 1, 256 * scale, 256 * scale, 1, 1, 16, 1, 8 };
+	return bench_generic_fft(dims, 6ul);
+}
+
+
+
+
+static double bench_generic_fftmod(long dims[DIMS], unsigned long flags)
+{
+	complex float* x = md_alloc(DIMS, dims, CFL_SIZE);
+	complex float* y = md_alloc(DIMS, dims, CFL_SIZE);
+
+	md_gaussian_rand(DIMS, dims, x);
+
+	double tic = timestamp();
+
+	fftmod(DIMS, dims, flags, y, x);
+
+	double toc = timestamp();
+
+	md_free(x);
+	md_free(y);
+
+	return toc - tic;
+}
+
+
+
+static double bench_fftmod(long scale)
+{
+	long dims[DIMS] = { 1, 256 * scale, 256 * scale, 1, 1, 16, 1, 16 };
+	return bench_generic_fftmod(dims, 6ul);
+}
+
+
+
+
+
 enum bench_indices { REPETITION_IND, SCALE_IND, THREADS_IND, TESTS_IND, BENCH_DIMS };
 
 typedef double (*bench_fun)(long scale);
@@ -472,6 +562,9 @@ const struct benchmark_s {
 	{ bench_copy1,		"copy 1" },
 	{ bench_copy2,		"copy 2" },
 	{ bench_wavelet,	"wavelet soft thresh" },
+	{ bench_mdfft,		"(MD-)FFT" },
+	{ bench_fft,		"FFT" },
+	{ bench_fftmod,		"fftmod" },
 };
 
 

src/cabs.c

@@ -11,6 +11,7 @@
 #include <complex.h>
 #include <stdio.h>
 
+#include "num/init.h"
 #include "num/flpmath.h"
 
 #include "misc/mmio.h"
@@ -28,6 +29,8 @@ int main_cabs(int argc, char* argv[])
 {
 	mini_cmdline(&argc, argv, 2, usage_str, help_str);
 
+	num_init();
+
 	long dims[DIMS];
 
 	complex float* idata = load_cfl(argv[1], DIMS, dims);

src/calib/calib.c

@@ -514,6 +514,7 @@ void calib2(const struct ecalib_conf* conf, const long out_dims[DIMS], complex f
 
 	if (conf->rotphase) {
 
+		// rotate the the phase with respect to the first principle component
 		long scc_dims[DIMS] = MD_INIT_ARRAY(DIMS, 1);
 		scc_dims[COIL_DIM] = channels;
 		scc_dims[MAPS_DIM] = channels;
@@ -564,7 +565,6 @@ void calib2(const struct ecalib_conf* conf, const long out_dims[DIMS], complex f
 	debug_printf(DP_DEBUG1, "Fix phase...\n");
 
 
-	// rotate the the phase with respect to the first principle component
 	fixphase2(DIMS, out_dims, COIL_DIM, rot[0], out_data, out_data);
 
 	md_free(imgcov);

src/calib/estvar.c

@@ -108,24 +108,24 @@ static char* file_name(const long kernel_dims[3], const long calreg_dims[4]) {
  *  L              - Number of elements in E.
  *  E              - Load simulated noise singular values to.
  */
-static int load_noise_sv(const long kernel_dims[3], const long calreg_dims[4], long L, float* E) {
-
+static int load_noise_sv(const long kernel_dims[3], const long calreg_dims[4], long L, float* E)
+{
     char* name = file_name(kernel_dims, calreg_dims);
     FILE* fp = fopen(name, "rb");
 
     if (!fp) {
-        free(name);
+
+        xfree(name);
         return 0;
     }
 
     int c = fread(E, sizeof(float), L, fp);
     assert(c == L);
     
-    free(name);
+    xfree(name);
     fclose(fp);
 
     return 1;
-
 }
 
 /**
@@ -139,13 +139,14 @@ static int load_noise_sv(const long kernel_dims[3], const long calreg_dims[4], l
  *  L              - Number of elements in E.
  *  E              - Load simulated noise singular values to.
  */
-static void save_noise_sv(const long kernel_dims[3], const long calreg_dims[4], long L, float* E) {
-
+static void save_noise_sv(const long kernel_dims[3], const long calreg_dims[4], long L, float* E)
+{
     char* name = file_name(kernel_dims, calreg_dims);
     FILE* fp = fopen(name, "wb");
 
     if (!fp) {
-        free(name);
+
+        xfree(name);
         return;
     }
 
@@ -153,7 +154,6 @@ static void save_noise_sv(const long kernel_dims[3], const long calreg_dims[4],
 
     free(name);
     fclose(fp);
-
 }
 
 /**
@@ -288,6 +288,7 @@ extern float estvar_calreg(const long kernel_dims[3], const long calreg_dims[4],
     PTR_ALLOC(complex float[N][N], vec);
     covariance_function(kernel_dims, N, *vec, calreg_dims, calreg);
     lapack_eig(N, tmpE, *vec);
+    PTR_FREE(vec);
 
     for (int idx = 0; idx < L; idx ++)
         S[idx] = sqrtf(tmpE[N-idx-1]);
@@ -296,13 +297,12 @@ extern float estvar_calreg(const long kernel_dims[3], const long calreg_dims[4],
 
 }
 
-extern float estvar_kspace(long N, const long kernel_dims[3], const long calib_size[3], const long kspace_dims[N], const complex float* kspace) {
-
+extern float estvar_kspace(long N, const long kernel_dims[3], const long calib_size[3], const long kspace_dims[N], const complex float* kspace)
+{
     long calreg_dims[N];
     complex float* calreg = NULL;
     calreg = extract_calib(calreg_dims, calib_size, kspace_dims, kspace, false);
     float variance = estvar_calreg(kernel_dims, calreg_dims, calreg);
     md_free(calreg);
     return variance;
-
 }

src/dfwavelet/prox_dfwavelet.c

@@ -274,7 +274,7 @@ struct prox_4pt_dfwavelet_data* prepare_prox_4pt_dfwavelet_data(const long im_di
 	md_calc_strides(DIMS, strs, data->tim_dims, CFL_SIZE);
 
 	data->w_op = linop_wavelet_create(DIMS, FFT_FLAGS, data->tim_dims, strs, min_size, false);
-        data->wthresh_op = prox_unithresh_create(DIMS, data->w_op, lambda, MD_BIT(data->flow_dim), use_gpu);
+        data->wthresh_op = prox_unithresh_create(DIMS, data->w_op, lambda, MD_BIT(data->flow_dim));
 
         return PTR_PASS(data);
 }

src/ecalib.c

@@ -65,6 +65,7 @@ int main_ecalib(int argc, char* argv[])
 		OPT_SET('W', &conf.weighting, "soft-weighting of the singular vectors."),
 		OPT_SET('I', &conf.intensity, "intensity correction"),
 		OPT_SET('1', &one, "perform only first part of the calibration"),
+		OPT_CLEAR('P', &conf.rotphase, "Do not rotate the phase with respect to the first principal component"),
 		OPT_CLEAR('O', &conf.orthiter, "()"),
 		OPT_FLOAT('b', &conf.perturb, "", "()"),
 		OPT_SET('V', &print_svals, "()"),

src/ecaltwo.c

@@ -47,7 +47,7 @@ int main_ecaltwo(int argc, char* argv[])
 
 		OPT_FLOAT('c', &conf.crop, "crop_value", "Crop the sensitivities if the eigenvalue is smaller than {crop_value}."),
 		OPT_LONG('m', &maps, "maps", "Number of maps to compute."),
-		OPT_SET('S', &conf.softcrop, "()"),
+		OPT_SET('S', &conf.softcrop, "Create maps with smooth transitions (Soft-SENSE)."),
 		OPT_CLEAR('O', &conf.orthiter, "()"),
 		OPT_SET('g', &conf.usegpu, "()"),
 	};

src/grecon/optreg.c

@@ -160,6 +160,13 @@ bool opt_reg(void* ptr, char c, const char* optarg)
 			assert(2 == ret);
 			regs[r].xflags = 0u;
 		}
+		else if (strcmp(rt, "S") == 0) {
+
+			regs[r].xform = POS;
+			regs[r].lambda = 0u;
+			regs[r].xflags = 0u;
+			regs[r].jflags = 0u;
+		}
 		else if (strcmp(rt, "Q") == 0) {
 
 			regs[r].xform = L2IMG;
@@ -306,10 +313,9 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 		case NIHTWAV:
 		{
 			debug_printf(DP_INFO, "NIHT with wavelets regularization: k = %d%% of total elements in each wavelet transform\n", regs[nr].k);
-			if (use_gpu){
-				debug_printf(DP_WARN, "GPU operation is not currently implemented for NIHT.\nContinuing with CPU.\n");
-				use_gpu = false; // not implemented, TODO: implement NIHT with gpu
-			}
+
+			if (use_gpu)
+				error("GPU operation is not currently implemented for NIHT.\n");
 
 			long img_strs[N];
 			md_calc_strides(N, img_strs, img_dims, CFL_SIZE);
@@ -344,24 +350,24 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 				debug_printf(DP_DEBUG3,"%d ", wav_dims[i]);
 			debug_printf(DP_DEBUG3, "]\n");
 			
-			prox_ops[nr] = prox_niht_thresh_create(N, wav_dims, K, regs[nr].jflags, use_gpu );			
+			prox_ops[nr] = prox_niht_thresh_create(N, wav_dims, K, regs[nr].jflags);
 			break;
 		}
 
 		case NIHTIM:
 		{
 			debug_printf(DP_INFO, "NIHT regularization in the image domain: k = %d%% of total elements in image vector\n", regs[nr].k);
-			if (use_gpu){
-				debug_printf(DP_WARN, "GPU operation is not currently implemented for NIHT.\nContinuing with CPU.\n");
-				use_gpu = false; // not implemented, TODO: implement NIHT with gpu
-			}
+
+			if (use_gpu)
+				error("GPU operation is not currently implemented for NIHT.\n");
 
 			long thresh_dims[N];
 			md_select_dims(N, regs[nr].xflags, thresh_dims, img_dims);		
 			unsigned int K = (md_calc_size(N, thresh_dims) / 100) * regs[nr].k;
 			debug_printf(DP_INFO, "k = %d%%, actual K = %d\n", regs[nr].k, K);
 
-			prox_ops[nr] = prox_niht_thresh_create(N, img_dims, K, regs[nr].jflags, use_gpu );
+			trafos[nr] = linop_identity_create(DIMS, img_dims);
+			prox_ops[nr] = prox_niht_thresh_create(N, img_dims, K, regs[nr].jflags);
 			debug_printf(DP_INFO, "NIHTIM initialization complete\n");
 			break;
 		}
@@ -372,7 +378,7 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 			trafos[nr] = linop_grad_create(DIMS, img_dims, regs[nr].xflags);
 			prox_ops[nr] = prox_thresh_create(DIMS + 1,
 					linop_codomain(trafos[nr])->dims,
-					regs[nr].lambda, regs[nr].jflags | MD_BIT(DIMS), use_gpu);
+					regs[nr].lambda, regs[nr].jflags | MD_BIT(DIMS));
 			break;
 
 		case LAPLACE:
@@ -394,13 +400,17 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 			trafos[nr] = linop_conv_create(DIMS, regs[nr].xflags, CONV_TRUNCATED, CONV_SYMMETRIC, img_dims, img_dims, krn_dims, krn);
 			prox_ops[nr] = prox_thresh_create(DIMS,
 					linop_codomain(trafos[nr])->dims,
-					regs[nr].lambda, regs[nr].jflags, use_gpu);
+					regs[nr].lambda, regs[nr].jflags);
 			break;
 
 		case LLR:
 
 			debug_printf(DP_INFO, "lowrank regularization: %f\n", regs[nr].lambda);
 
+			if (use_gpu)
+				error("GPU operation is not currently implemented for lowrank regularization.\n");
+
+
 			// add locally lowrank penalty
 			levels = llr_blkdims(blkdims, regs[nr].jflags, img_dims, llr_blk);
 
@@ -418,7 +428,7 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 			int remove_mean = 0;
 
 			trafos[nr] = linop_identity_create(DIMS, img_dims);
-			prox_ops[nr] = lrthresh_create(img_dims, randshift, regs[nr].xflags, (const long (*)[DIMS])blkdims, regs[nr].lambda, false, remove_mean, use_gpu);
+			prox_ops[nr] = lrthresh_create(img_dims, randshift, regs[nr].xflags, (const long (*)[DIMS])blkdims, regs[nr].lambda, false, remove_mean);
 			break;
 
 		case MLR:
@@ -445,7 +455,7 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 			linop_free(decom_op);
 			linop_free(tmp_op);
 #else
-			debug_printf(DP_WARN, "multi-scale lowrank regularization not yet supported: %f\n", regs[nr].lambda);
+			error("multi-scale lowrank regularization not supported.\n");
 #endif
 
 			break;
@@ -454,7 +464,7 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 			debug_printf(DP_INFO, "l1 regularization of imaginary part: %f\n", regs[nr].lambda);
 
 			trafos[nr] = linop_rdiag_create(DIMS, img_dims, 0, &(complex float){ 1.i });
-			prox_ops[nr] = prox_thresh_create(DIMS, img_dims, regs[nr].lambda, regs[nr].jflags, use_gpu);
+			prox_ops[nr] = prox_thresh_create(DIMS, img_dims, regs[nr].lambda, regs[nr].jflags);
 			break;
 
 		case IMAGL2:
@@ -468,9 +478,17 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 			debug_printf(DP_INFO, "l1 regularization: %f\n", regs[nr].lambda);
 
 			trafos[nr] = linop_identity_create(DIMS, img_dims);
-			prox_ops[nr] = prox_thresh_create(DIMS, img_dims, regs[nr].lambda, regs[nr].jflags, use_gpu);
+			prox_ops[nr] = prox_thresh_create(DIMS, img_dims, regs[nr].lambda, regs[nr].jflags);
+			break;
+
+		case POS:
+			debug_printf(DP_INFO, "non-negative constraint\n");
+
+			trafos[nr] = linop_identity_create(DIMS, img_dims);
+			prox_ops[nr] = prox_nonneg_create(DIMS, img_dims);
 			break;
 
+
 		case L2IMG:
 			debug_printf(DP_INFO, "l2 regularization: %f\n", regs[nr].lambda);
 
@@ -482,7 +500,7 @@ void opt_reg_configure(unsigned int N, const long img_dims[N], struct opt_reg_s*
 			debug_printf(DP_INFO, "l1 regularization of Fourier transform: %f\n", regs[nr].lambda);
 
 			trafos[nr] = linop_fft_create(DIMS, img_dims, regs[nr].xflags);
-			prox_ops[nr] = prox_thresh_create(DIMS, img_dims, regs[nr].lambda, regs[nr].jflags, use_gpu);
+			prox_ops[nr] = prox_thresh_create(DIMS, img_dims, regs[nr].lambda, regs[nr].jflags);
 			break;
 		}
 
@@ -498,6 +516,5 @@ void opt_reg_free(struct opt_reg_s* ropts, const struct operator_p_s* prox_ops[N
 
 		operator_p_free(prox_ops[nr]);
 		linop_free(trafos[nr]);
-
 	}
 }

src/grecon/optreg.h

@@ -18,7 +18,7 @@ enum algo_t { CG, IST, FISTA, ADMM, NIHT, PRIDU };
 
 struct reg_s {
 
-	enum { L1WAV, NIHTWAV, NIHTIM, TV, LLR, MLR, IMAGL1, IMAGL2, L1IMG, L2IMG, FTL1, LAPLACE } xform;
+	enum { L1WAV, NIHTWAV, NIHTIM, TV, LLR, MLR, IMAGL1, IMAGL2, L1IMG, L2IMG, FTL1, LAPLACE, POS } xform;
 
 	unsigned int xflags;
 	unsigned int jflags;

src/homodyne.c

-/* Copyright 2014-2015. The Regents of the University of California.
+/* Copyright 2014-2017. The Regents of the University of California.
  * Copyright 2015-2016. Martin Uecker.
  * All rights reserved. Use of this source code is governed by
  * a BSD-style license which can be found in the LICENSE file.
  *
  * Authors:
  * 2013-2016 Martin Uecker <martin.uecker@med.uni-goettingen.de>
- * 2015 Jonathan Tamir <jtamir@eecs.berkeley.edu>
+ * 2015, 2017 Jonathan Tamir <jtamir@eecs.berkeley.edu>
  */
 
 #include <stdbool.h>
@@ -20,6 +20,7 @@
 #include "misc/mmio.h"
 #include "misc/misc.h"
 #include "misc/opts.h"
+#include "misc/nested.h"
 
 
 
@@ -34,13 +35,10 @@ static const char help_str[] = "Perform homodyne reconstruction along dimension
 struct wdata {
 
 	float frac;
-	float alpha;
 	int pfdim;
 	long wdims[DIMS];
 	long wstrs[DIMS];
 	complex float* weights;
-
-	bool clear;
 };
 
 
@@ -74,15 +72,10 @@ static float homodyne_filter(long N, float frac, float alpha, bool clear, long p
 	return ret;
 }
 
-static void comp_weights(void* _data, const long pos[])
-{
-	struct wdata* data = _data;
-	data->weights[md_calc_offset(DIMS, data->wstrs, pos) / CFL_SIZE] 
-		= homodyne_filter(data->wdims[data->pfdim], data->frac, data->alpha, data->clear, pos[data->pfdim]);
-}
+
 
 static complex float* estimate_phase(struct wdata wdata, unsigned int flags,
-		unsigned int N, const long dims[N], const complex float* idata)
+		unsigned int N, const long dims[N], const complex float* idata, bool center_fft)
 {
 
 	long cdims[N];
@@ -97,7 +90,7 @@ static complex float* estimate_phase(struct wdata wdata, unsigned int flags,
 	md_resize_center(N, dims, phase, cdims, center, CFL_SIZE);
 	md_free(center);
 
-	ifftuc(N, dims, flags, phase, phase);
+	(center_fft ? ifftuc : ifftu)(N, dims, flags, phase, phase);
 	md_zphsr(N, dims, phase, phase);
 
 	return phase;
@@ -105,10 +98,10 @@ static complex float* estimate_phase(struct wdata wdata, unsigned int flags,
 
 static void homodyne(struct wdata wdata, unsigned int flags, unsigned int N, const long dims[N],
 		const long strs[N], complex float* data, const complex float* idata,
-		const long pstrs[N], const complex float* phase)
+		const long pstrs[N], const complex float* phase, bool center_fft)
 {
 	md_zmul2(N, dims, strs, data, strs, idata, wdata.wstrs, wdata.weights);
-	ifftuc(N, dims, flags, data, data);
+	(center_fft ? ifftuc : ifftu)(N, dims, flags, data, data);
 
 	md_zmulc2(N, dims, strs, data, strs, data, pstrs, phase);
 	md_zreal(N, dims, data, data);
@@ -120,6 +113,7 @@ int main_homodyne(int argc, char* argv[])
 {
 	bool clear = false;
 	bool image = false;
+	bool center_fft = true;
 	const char* phase_ref = NULL;
 
 	float alpha = 0.;
@@ -132,6 +126,7 @@ int main_homodyne(int argc, char* argv[])
 		OPT_SET('I', &image, "Input is in image domain"),
 		OPT_SET('C', &clear, "Clear unacquired portion of kspace"),
 		OPT_STRING('P', &phase_ref, "phase_ref>", "Use <phase_ref> as phase reference"),
+		OPT_CLEAR('n', &center_fft, "use uncentered ffts"),
 	};
 
 	cmdline(&argc, argv, 4, 4, usage_str, help_str, ARRAY_SIZE(opts), opts);
@@ -150,7 +145,7 @@ int main_homodyne(int argc, char* argv[])
 
 	if (image) {
 		complex float* ksp_in = md_alloc(N, dims, CFL_SIZE);
-		fftuc(N, dims, FFT_FLAGS, ksp_in, idata);
+		(center_fft ? fftuc : fftu)(N, dims, FFT_FLAGS, ksp_in, idata);
 		md_copy(N, dims, idata, ksp_in, CFL_SIZE);
 		md_free(ksp_in);
 	}
@@ -165,10 +160,14 @@ int main_homodyne(int argc, char* argv[])
 	md_select_dims(N, MD_BIT(pfdim), wdata.wdims, dims);
 	md_calc_strides(N, wdata.wstrs, wdata.wdims, CFL_SIZE);
 	wdata.weights = md_alloc(N, wdata.wdims, CFL_SIZE);
-	wdata.alpha = alpha;
-	wdata.clear = clear;
 
-	md_loop(N, wdata.wdims, &wdata, comp_weights);
+	NESTED(void, comp_weights, (const long pos[]))
+	{
+		wdata.weights[md_calc_offset(DIMS, wdata.wstrs, pos) / CFL_SIZE]
+			= homodyne_filter(wdata.wdims[pfdim], frac, alpha, clear, pos[pfdim]);
+	};
+
+	md_loop(N, wdata.wdims, comp_weights);
 
 	long pstrs[N];
 	long pdims[N];
@@ -176,7 +175,7 @@ int main_homodyne(int argc, char* argv[])
 
 	if (NULL == phase_ref) {
 
-		phase = estimate_phase(wdata, FFT_FLAGS, N, dims, idata);
+		phase = estimate_phase(wdata, FFT_FLAGS, N, dims, idata, center_fft);
 		md_copy_dims(N, pdims, dims);
 	}
 	else
@@ -184,7 +183,7 @@ int main_homodyne(int argc, char* argv[])
 
 	md_calc_strides(N, pstrs, pdims, CFL_SIZE);
 
-	homodyne(wdata, FFT_FLAGS, N, dims, strs, data, idata, pstrs, phase);
+	homodyne(wdata, FFT_FLAGS, N, dims, strs, data, idata, pstrs, phase, center_fft);
 
 	md_free(wdata.weights);
 

src/iter/admm.c

-/* Copyright 2014-2016. The Regents of the University of California.
+/* Copyright 2014-2018. The Regents of the University of California.
  * Copyright 2016-2017. Martin Uecker.
  * All rights reserved. Use of this source code is governed by
  * a BSD-style license which can be found in the LICENSE file.
  *
  * Authors:
  * 2014-2017 Martin Uecker <martin.uecker@med.uni-goettingen.de>
- * 2014-2016 Jonathan Tamir <jtamir@eecs.berkeley.edu>
+ * 2014-2016, 2017 Jon Tamir <jtamir@eecs.berkeley.edu>
  * 
  *
  *
@@ -20,7 +20,6 @@
  * via finite element approximation
  * Computers & Mathematics with Applications, 2:17-40 (1976)
  *
- *
  * Afonso MA, Bioucas-Dias JM, Figueiredo M. An Augmented Lagrangian Approach to
  * the Constrained Optimization Formulation of Imaging Inverse Problems,
  * IEEE Trans Image Process, 20:681-695 (2011)
@@ -29,6 +28,9 @@
  * Statistical Learning via the Alternating Direction Method of Multipliers,
  * Foundations and Trends in Machine Learning, 3:1-122 (2011)
  *
+ * Wohlberg B. ADMM Penalty Parameter Selection by Residual Balancing,
+ * arXiv:1704.06209 (2017)
+ *
  */
 
 #include <math.h>
@@ -197,6 +199,7 @@ void admm(const struct admm_plan_s* plan,
 		GH_usum = vops->allocate(N);
 
 	float rho = plan->rho;
+	float tau = plan->tau;
 
 	struct admm_normaleq_data ndata = {
 
@@ -388,6 +391,9 @@ void admm(const struct admm_plan_s* plan,
 		float s_norm = 0.;
 		float r_norm = 0.;
 
+		double s_scaling = 1.;
+		double r_scaling = 1.;
+
 		if (plan->dynamic_rho || !plan->fast) {
 
 			s_norm = rho * vops->norm(N, s);
@@ -406,34 +412,38 @@ void admm(const struct admm_plan_s* plan,
 			for (unsigned int j = 0; j < num_funs; j++)
 				n2 += pow(vops->norm(z_dims[j], z[j]), 2.);
 
-			double n = MAX(MAX(n1, n2), n3);
+			r_scaling = sqrt(MAX(MAX(n1, n2), n3));
+			s_scaling = rho * vops->norm(N, GH_usum);
 
 			long M = sum_long_array(num_funs, z_dims);
 
-			float eps_pri = plan->ABSTOL * sqrt(M) + plan->RELTOL * sqrt(n);
-			float eps_dual = plan->ABSTOL * sqrt(N) + plan->RELTOL * rho * vops->norm(N, GH_usum);
+			float eps_pri = plan->ABSTOL * sqrt(M) + plan->RELTOL * r_scaling;
+			float eps_dual = plan->ABSTOL * sqrt(N) + plan->RELTOL * s_scaling;
 
 
 			struct admm_history_s history;
 
 			history.s_norm = s_norm;
 			history.r_norm = r_norm;
+			history.s_scaling = s_scaling;
+			history.r_scaling = r_scaling;
 			history.eps_pri = eps_pri;
 			history.eps_dual = eps_dual;
 			history.rho = rho;
+			history.tau = tau;
 			history.numiter = i;
 			history.nr_invokes = ndata.nr_invokes;
 
 			iter_history(monitor, CAST_UP(&history));
 
 			if (0 == i)
-				debug_printf(DP_DEBUG2, "%3s\t%3s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n",
-					"iter", "cgiter", "rho", "r norm", "eps pri",
+				debug_printf(DP_DEBUG2, "%3s\t%3s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n",
+					"iter", "cgiter", "rho", "tau", "r norm", "eps pri",
 					"s norm", "eps dual", "obj", "relMSE");
 
 
-			debug_printf(DP_DEBUG2, "%3d\t%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.4f\n",
-				history.numiter, history.nr_invokes, history.rho,
+			debug_printf(DP_DEBUG2, "%3d\t%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.4f\n",
+				history.numiter, history.nr_invokes, history.rho, history.tau,
 				history.r_norm, history.eps_pri, history.s_norm, history.eps_dual,
 				(NULL == monitor) ? -1. : monitor->obj,
 				(NULL == monitor) ? -1. : monitor->err);
@@ -456,13 +466,33 @@ void admm(const struct admm_plan_s* plan,
 
 		assert(!(plan->dynamic_rho && plan->hogwild));
 
+		if (plan->dynamic_tau) {
+
+			double t = sqrt(r_norm / s_norm);
+			if (plan->tau_max > t && 1 <= t)
+				tau = t;
+			else if (1 > t && (1/plan->tau_max) < t)
+				tau = 1. / t;
+			else
+				tau = plan->tau_max;
+		}
+
 		if (plan->dynamic_rho) {
 
-			if (r_norm > plan->mu * s_norm)
-				sc = plan->tau;
+			double r = r_norm;
+			double s = s_norm;
+
+			if (plan->relative_norm) {
+
+				r /= r_scaling;
+				s /= s_scaling;
+			}
+
+			if (r > plan->mu * s)
+				sc = tau;
 			else
-			if (s_norm > plan->mu * r_norm)
-				sc = 1. / plan->tau;
+			if (s > plan->mu * r)
+				sc = 1. / tau;
 		}
 
 		if (plan->hogwild) {