diff --git a/Makefile b/Makefile
index 2bb17d82a..c7f5d5344 100644
--- a/Makefile
+++ b/Makefile
@@ -95,11 +95,11 @@ ismrm.top := /usr/local/ismrmrd/
 # Main build targets
 
 TBASE=slice crop resize join transpose zeros ones flip circshift extract repmat bitmask
-TFLP=scale conj fmac saxpy sdot conv
-TNUM=fft fftmod noise bench threshold creal
+TFLP=scale conj fmac saxpy sdot spow cpyphs creal
+TNUM=fft fftmod noise bench threshold conv
 TRECO=sense pocsense rsense bpsense itsense nlinv nufft rof nusense
-TCALIB=ecalib caldir walsh
-TMRI=rss homodyne pattern scc poisson twixread
+TCALIB=ecalib caldir walsh cc
+TMRI=rss homodyne pattern poisson twixread
 TSIM=phantom traj
 TARGETS = $(TBASE) $(TFLP) $(TNUM) $(TRECO) $(TCALIB) $(TMRI) $(TSIM) $(TIO)
 
@@ -120,6 +120,7 @@ MODULES_itsense = -liter -llinops
 MODULES_ecalib = -lcalib
 MODULES_caldir = -lcalib
 MODULES_walsh = -lcalib
+MODULES_cc = -lcalib
 MODULES_nufft = -lnoncart -liter -llinops
 MODULES_rof = -liter -llinops
 MODULES_bench = -lwavelet2 -llinops
diff --git a/README b/README
index f5e0f2eba..6a56bab14 100644
--- a/README
+++ b/README
@@ -278,6 +278,11 @@ saxpy scale <input1> <input2> <output>
 	Multiply input1 with scale factor and add input2.
 
 
+spow exponent <input> <output>
+
+	Raise array to the power of {exponent}. The exponent can be a complex number.
+
+
 fmac [-A] [-C] [-s bitmask] <input1> <input2> <output>
 
 	Multiply and accumulate.
@@ -292,6 +297,11 @@ conj <input> <output>
 	Compute complex conjugate.
 
 
+cpyphs <input> <output
+
+	Copy phase from <input> to <output>.
+
+
 fft [-u] [-i] bitmask <input> <output>
 
 	Performs a fast Fourier transform (FFT) along dimensions
@@ -380,14 +390,14 @@ ecaltwo [-c crop] [-m maps] x y z <input> <sensitivities> [<ev_maps>]
 	-m maps		Number of maps to compute.
 
 	
-scc [-v] [-A] [-r cal_size] [-P num_coeffs] <kspace> <coeff>|<proj_kspace>
+cc [-A] [-r cal_size] [-P num_coeffs] <kspace> <coeff>|<proj_kspace>
 
 	Performs simple coil compression.
 
 	-P N	perform compression to N virtual channels
 	-r S	size of calibration region
 	-A	use all data to compute coefficients
-	-v	verbose
+	-S|G|E	type: SVD, Geometric, ESPIRiT
 	-h	help
 
 
@@ -419,7 +429,7 @@ nusense [-l1/-l2] [-r lambda] <traj> <kspace> <sensitivities> <output>
 	-r lambda	regularization parameter
 
 
-bpsense [-g] [-r l2lambda] [-c] [-e eps] [-u rho] <kspace> <sensitivities> <output>
+bpsense [-r l2lambda] [-c] [-e eps] [-u rho] <kspace> <sensitivities> <output>
 
 	Perform basis pursuit denoising for SENSE/ESPIRiT reconstruction:
 	min_x ||T x||_1 + lambda/2 ||x||_2^2 subject to: ||y - Ax||_2 <= eps
@@ -429,7 +439,7 @@ bpsense [-g] [-r l2lambda] [-c] [-e eps] [-u rho] <kspace> <sensitivities> <outp
 	-u rho	ADMM penalty parameter
 	-c	real-value constraint
 	-t	use TV norm
-	-F	truth image
+	-T	truth image
 
 
 nlinv [-i iterations] <kspace> <output> [<sens>]
@@ -479,7 +489,7 @@ rof lambda bitmask <input> <output>
 
 poisson <outfile>
 
-Computes Poisson-disc sampling pattern.
+	Computes Poisson-disc sampling pattern.
 
 	-X	size dimension 0 (readout)
 	-Y	size dimension 1 (phase)
@@ -526,6 +536,8 @@ twixread [...] [-a A] <dat file> <output>
 	-s S	number of slices
 	-c C	number of channels
 	-a A	total number of ADCs
+	-L	use linectr offset
+	-P	use partctr offset
 	-h	help
 
 
diff --git a/src/bpsense.c b/src/bpsense.c
index db4b40091..8b31ef223 100755
--- a/src/bpsense.c
+++ b/src/bpsense.c
@@ -43,7 +43,7 @@
 
 static void usage(const char* name, FILE* fd)
 {
-	fprintf(fd, "Usage: %s [-g] [-r l2lambda] [-c] [-e eps] [-u rho] <kspace> <sensitivities> <output>\n", name);
+	fprintf(fd, "Usage: %s [-r l2lambda] [-c] [-e eps] [-u rho] <kspace> <sensitivities> <output>\n", name);
 }
 
 static void help(void)
@@ -57,7 +57,7 @@ static void help(void)
 		"-u rho\tADMM penalty parameter\n"
 		"-c\treal-value constraint\n"
 		"-t\tuse TV norm\n"
-		"-F\ttruth image\n");
+		"-T\ttruth image\n");
 }
 
 
@@ -75,18 +75,20 @@ int main(int argc, char* argv[])
 
 	bool usegpu = false;
 	const char* psf = NULL;
-	char image_truth_fname[100];
-	_Bool im_truth = false;
-	_Bool use_tvnorm = false;
+	const char* image_truth_fname = NULL;
+	bool im_truth = false;
+	bool use_tvnorm = false;
 
 	double start_time = timestamp();
 
 	int c;
-	while (-1 != (c = getopt(argc, argv, "F:r:e:i:u:p:tcgh"))) {
+	while (-1 != (c = getopt(argc, argv, "T:r:e:i:u:p:tcgh"))) {
 		switch(c) {
-		case 'F':
+
+		case 'T':
 			im_truth = true;
-			sprintf(image_truth_fname, "%s", optarg);
+			image_truth_fname = strdup(optarg);
+			assert(NULL != image_truth_fname);
 			break;
 
 		case 'r':
@@ -292,6 +294,12 @@ int main(int argc, char* argv[])
 	unmap_cfl(N, ksp_dims, kspace_data);
 	unmap_cfl(N, img_dims, image);
 
+	if (im_truth) {
+
+		free((void*)image_truth_fname);
+		unmap_cfl(DIMS, img_dims, image_truth);
+	}
+
 	double end_time = timestamp();
 	debug_printf(DP_INFO, "Total Time: %f\n", end_time - start_time);
 
diff --git a/src/caldir.c b/src/caldir.c
index 9a6defe8f..0357892de 100644
--- a/src/caldir.c
+++ b/src/caldir.c
@@ -36,9 +36,9 @@ int main_caldir(int argc, char* argv[])
 {
 	mini_cmdline(argc, argv, 3, usage_str, help_str);
 
-	long dims[KSPACE_DIMS];
+	long dims[DIMS];
 
-	complex float* in_data = load_cfl(argv[2], KSPACE_DIMS, dims);
+	complex float* in_data = load_cfl(argv[2], DIMS, dims);
 
 	int calsize_ro = atoi(argv[1]);
 	long calsize[3] = { calsize_ro, calsize_ro, calsize_ro };
@@ -46,10 +46,10 @@ int main_caldir(int argc, char* argv[])
 	assert((dims[0] == 1) || (calsize_ro < dims[0]));
 	assert(1 == dims[4]);
 	
-	complex float* out_data = create_cfl(argv[3], KSPACE_DIMS, dims);
+	complex float* out_data = create_cfl(argv[3], DIMS, dims);
 
 
-	long caldims[KSPACE_DIMS];
+	long caldims[DIMS];
 	complex float* cal_data = extract_calib(caldims, calsize, dims, in_data, false);
 
 	printf("Calibration region %ldx%ldx%ld\n", caldims[0], caldims[1], caldims[2]);
@@ -60,8 +60,8 @@ int main_caldir(int argc, char* argv[])
 
 	md_free(cal_data);
 
-	unmap_cfl(KSPACE_DIMS, dims, (void*)out_data);
-	unmap_cfl(KSPACE_DIMS, dims, (void*)in_data);
+	unmap_cfl(DIMS, dims, (void*)out_data);
+	unmap_cfl(DIMS, dims, (void*)in_data);
 
 	exit(0);
 }
diff --git a/src/calib/calib.c b/src/calib/calib.c
index 5fc2a3a7b..121b811a0 100644
--- a/src/calib/calib.c
+++ b/src/calib/calib.c
@@ -64,6 +64,10 @@ static void eigen_herm3(int M, int N, float val[M], complex float matrix[N][N])
 
 
 
+/*
+ * rotate phase jointly along dim so that the 0-th slice along dim has phase = 0
+ *
+ */
 void fixphase(unsigned int N, const long dims[N], unsigned int dim, complex float* out, const complex float* in)
 {
 	assert(dim < N);
diff --git a/src/calib/cc.c b/src/calib/cc.c
new file mode 100644
index 000000000..c921baa4f
--- /dev/null
+++ b/src/calib/cc.c
@@ -0,0 +1,204 @@
+/* Copyright 2014. The Regents of the University of California.
+ * All rights reserved. Use of this source code is governed by
+ * a BSD-style license which can be found in the LICENSE file.
+ *
+ * Authors:
+ * 2012-2014	Martin Uecker <uecker@eecs.berkeley.edu>
+ * 2013		Dara Bahri <dbahri123@gmail.com>
+ *
+ *
+ * Huang F, Vijayakumar S, Li Y, Hertel S, Duensing GR. A software channel
+ * compression technique for faster reconstruction with many channels.
+ * Magn Reson Imaging 2008; 26:133-141.
+ *
+ * Buehrer M, Pruessmann KP, Boesiger P, Kozerke S. Array compression for MRI
+ * with large coil arrays. Magn Reson Med 2007, 57: 1131–1139.
+ *
+ * Zhang T, Pauly JM, Vasanawala SS, Lustig M. Coil compression for
+ * accelerated imaging with cartesian sampling. Magn Reson Med 2013;
+ * 69:571-582.
+ *
+ * Bahri D, Uecker M, Lustig M. ESPIRiT-Based Coil Compression for
+ * Cartesian Sampling,  Annual Meeting ISMRM, Salt Lake City 2013,
+ * In: Proc Intl Soc Mag Reson Med 21:2657
+ */
+ 
+#include <complex.h>
+#include <assert.h>
+
+#include "num/multind.h"
+#include "num/flpmath.h"
+#include "num/fft.h"
+#include "num/lapack.h"
+#include "num/la.h"
+
+#include "misc/debug.h"
+#include "misc/mri.h"
+
+#include "calib/calib.h"
+
+#include "cc.h"
+
+
+void scc(const long out_dims[DIMS], complex float* out_data, const long caldims[DIMS], const complex float* cal_data)
+{
+	int channels = caldims[COIL_DIM];
+
+	assert(1 == md_calc_size(3, out_dims));
+	assert(out_dims[COIL_DIM] == channels);
+	assert(out_dims[MAPS_DIM] == channels);
+
+	complex float tmp[channels][channels];
+
+
+	size_t csize = md_calc_size(3, caldims);
+	gram_matrix(channels, tmp, csize, (const complex float (*)[csize])cal_data);
+
+	float vals[channels];
+	eigendecomp(channels, vals, tmp);
+
+	md_flip(DIMS, out_dims, MAPS_FLAG, out_data, tmp, CFL_SIZE);
+
+
+	debug_printf(DP_DEBUG1, "Coefficients:");
+
+	for (int i = 0; i < channels; i++)
+		debug_printf(DP_DEBUG1, " %.3f", vals[channels - 1 - i] / vals[channels - 1]);
+
+	debug_printf(DP_DEBUG1, "\n");
+}
+
+
+
+static void align1(int M, int N, complex float out[M][N], const complex float in1[M][N], const complex float in2[M][N])
+{
+	assert(M <= N);
+#if 1
+	complex float in1T[N][M];
+	complex float C[M][M];
+	complex float U[M][M];
+	complex float VH[M][M];
+	float S[M];
+	complex float P[M][M];
+
+	mat_adjoint(M, N, in1T, in1);	// A_{x-1}^H
+	mat_mul(M, N, M, C, in2, in1T);	// C = A_{x} A_{x-1}^H
+	// VH and U are switched here because SVD uses column-major arrays
+	svd(M, M, VH, U, S, C);		// U S V^H = C
+	mat_mul(M, M, M, C, U, VH);	// U V^H
+	mat_adjoint(M, M, P, C);	// P_x = V U^H
+	mat_mul(M, M, N, out, P, in2);	// A_{x} <- P_x A_{x}
+#else
+	mat_copy(M, N, out, in2);
+#endif
+}
+
+
+
+void align_ro(const long dims[DIMS], complex float* odata, const complex float* idata)
+{
+	int ro = dims[READ_DIM];
+
+	long tmp_dims[DIMS];
+	md_select_dims(DIMS, ~READ_FLAG, tmp_dims, dims);
+
+	complex float* tmp1 = md_alloc(DIMS, tmp_dims, CFL_SIZE);
+	complex float* tmp2 = md_alloc(DIMS, tmp_dims, CFL_SIZE);
+	complex float* tmp3 = md_alloc(DIMS, tmp_dims, CFL_SIZE);
+
+	md_copy_block(DIMS, (long[DIMS]){ [READ_DIM] = 0 }, tmp_dims, tmp1, dims, idata, CFL_SIZE);
+
+	for (int i = 0; i < ro - 1; i++) {
+
+		md_copy_block(DIMS, (long[DIMS]){ [READ_DIM] = i + 1 }, tmp_dims, tmp2, dims, idata, CFL_SIZE);
+
+		align1(tmp_dims[MAPS_DIM], tmp_dims[COIL_DIM],
+				MD_CAST_ARRAY2(      complex float, DIMS, tmp_dims, tmp3, COIL_DIM, MAPS_DIM),
+				MD_CAST_ARRAY2(const complex float, DIMS, tmp_dims, tmp1, COIL_DIM, MAPS_DIM),
+				MD_CAST_ARRAY2(const complex float, DIMS, tmp_dims, tmp2, COIL_DIM, MAPS_DIM));
+
+		md_copy(DIMS, tmp_dims, tmp1, tmp3, CFL_SIZE);
+
+		md_copy_block(DIMS, (long[DIMS]){ [READ_DIM] = i + 1 }, dims, odata, tmp_dims, tmp3, CFL_SIZE);
+	}
+
+	md_free(tmp1);
+	md_free(tmp2);
+	md_free(tmp3);
+}
+
+
+void gcc(const long out_dims[DIMS], complex float* out_data, const long caldims[DIMS], const complex float* cal_data)
+{
+	int ro = out_dims[READ_DIM];
+
+	// zero pad calibration region along readout and FFT
+
+	long tmp_dims[DIMS];
+	md_copy_dims(DIMS, tmp_dims, caldims);
+	tmp_dims[READ_DIM] = ro;
+	complex float* tmp = md_alloc(DIMS, tmp_dims, CFL_SIZE);
+
+	md_resizec(DIMS, tmp_dims, tmp, caldims, cal_data, CFL_SIZE);
+	ifftuc(DIMS, tmp_dims, READ_FLAG, tmp, tmp);
+
+	// apply scc at each readout location
+
+	long tmp2_dims[DIMS];
+	md_select_dims(DIMS, ~READ_FLAG, tmp2_dims, tmp_dims);
+	complex float* tmp2 = md_alloc(DIMS, tmp2_dims, CFL_SIZE);
+
+	long out2_dims[DIMS];
+	md_select_dims(DIMS, ~READ_FLAG, out2_dims, out_dims);
+	complex float* out2 = md_alloc(DIMS, out2_dims, CFL_SIZE);
+
+	for (int i = 0; i < ro; i++) {
+
+		long pos[DIMS] = { [READ_DIM] = i };
+		md_copy_block(DIMS, pos, tmp2_dims, tmp2, tmp_dims, tmp, CFL_SIZE);
+
+		scc(out2_dims, out2, tmp2_dims, tmp2);
+
+		md_copy_block(DIMS, pos, out_dims, out_data, out2_dims, out2, CFL_SIZE);
+	}
+
+	md_free(out2);
+	md_free(tmp2);
+}
+
+
+
+void ecc(const long out_dims[DIMS], complex float* out_data, const long caldims[DIMS], const complex float* cal_data)
+{
+	int channels = caldims[COIL_DIM];
+
+	assert(1 == out_dims[PHS1_DIM]);
+	assert(1 == out_dims[PHS2_DIM]);
+	assert(out_dims[COIL_DIM] == channels);
+	assert(out_dims[MAPS_DIM] == channels);
+
+	struct ecalib_conf conf;
+	memcpy(&conf, &ecalib_defaults, sizeof(struct ecalib_conf));
+
+	conf.threshold = 0.001;
+	conf.crop = 0.;
+	conf.kdims[0] = 6;
+	conf.kdims[1] = 1;
+	conf.kdims[2] = 1;
+//	conf.numsv = L;
+	conf.weighting = false;
+	conf.orthiter = false;
+	conf.perturb = 0.;
+
+	long map_dims[DIMS];
+	md_select_dims(DIMS, ~MAPS_FLAG, map_dims, out_dims);
+        complex float* emaps = md_alloc(DIMS, map_dims, CFL_SIZE);
+
+	int K = conf.kdims[0] * caldims[COIL_DIM];
+	float svals[K];
+	calib(&conf, out_dims, out_data, emaps, K, svals, caldims, cal_data); 
+
+	md_free(emaps);
+}
+
+
diff --git a/src/calib/cc.h b/src/calib/cc.h
new file mode 100644
index 000000000..b61ebca0c
--- /dev/null
+++ b/src/calib/cc.h
@@ -0,0 +1,8 @@
+
+#include "misc/mri.h"
+
+extern void scc(const long out_dims[DIMS], complex float* out_data, const long caldims[DIMS], const complex float* cal_data);
+extern void gcc(const long out_dims[DIMS], complex float* out_data, const long caldims[DIMS], const complex float* cal_data);
+extern void ecc(const long out_dims[DIMS], complex float* out_data, const long caldims[DIMS], const complex float* cal_data);
+extern void align_ro(const long dims[DIMS], complex float* odata, const complex float* idata);
+
diff --git a/src/calib/walsh.c b/src/calib/walsh.c
index a78e7de90..9bbf8d437 100644
--- a/src/calib/walsh.c
+++ b/src/calib/walsh.c
@@ -28,7 +28,7 @@
 
 	
 
-void walsh(const long bsize[3], const long dims[KSPACE_DIMS], complex float* sens, const long caldims[KSPACE_DIMS], const complex float* data)
+void walsh(const long bsize[3], const long dims[DIMS], complex float* sens, const long caldims[DIMS], const complex float* data)
 {
 	assert(1 == caldims[MAPS_DIM]);
 	assert(1 == dims[MAPS_DIM]);
@@ -37,8 +37,8 @@ void walsh(const long bsize[3], const long dims[KSPACE_DIMS], complex float* sen
 	int cosize = channels * (channels + 1) / 2;
 	assert(dims[COIL_DIM] == cosize);
 
-	long dims1[KSPACE_DIMS];
-	md_copy_dims(KSPACE_DIMS, dims1, dims);
+	long dims1[DIMS];
+	md_copy_dims(DIMS, dims1, dims);
 	dims1[COIL_DIM] = channels;
 	
 	long kdims[4];
@@ -46,17 +46,17 @@ void walsh(const long bsize[3], const long dims[KSPACE_DIMS], complex float* sen
 	kdims[1] = MIN(bsize[1], dims[1]);
 	kdims[2] = MIN(bsize[2], dims[2]);
 
-	md_resizec(KSPACE_DIMS, dims1, sens, caldims, data, sizeof(complex float));
-	ifftc(KSPACE_DIMS, dims1, FFT_FLAGS, sens, sens);
+	md_resizec(DIMS, dims1, sens, caldims, data, CFL_SIZE);
+	ifftc(DIMS, dims1, FFT_FLAGS, sens, sens);
 
-	long odims[KSPACE_DIMS];
-	md_copy_dims(KSPACE_DIMS, odims, dims1);
+	long odims[DIMS];
+	md_copy_dims(DIMS, odims, dims1);
 
 	for (int i = 0; i < 3; i++)
 		odims[i] = dims[i] + kdims[i] - 1;
 
-	complex float* tmp = md_alloc(KSPACE_DIMS, odims, sizeof(complex float));
-	md_resizec(KSPACE_DIMS, odims, tmp, dims1, sens, sizeof(complex float));
+	complex float* tmp = md_alloc(DIMS, odims, CFL_SIZE);
+	md_resizec(DIMS, odims, tmp, dims1, sens, CFL_SIZE);
 
 	// FIXME: we should have the option to compute this from a periodic 
 	// extension 
diff --git a/src/calib/walsh.h b/src/calib/walsh.h
index bb8415240..5cb43caed 100644
--- a/src/calib/walsh.h
+++ b/src/calib/walsh.h
@@ -5,5 +5,5 @@
  
 #include "misc/mri.h"
 
-extern void walsh(const long bsize[3], const long dims[KSPACE_DIMS], _Complex float* sens, const long caldims[KSPACE_DIMS], const _Complex float* data);
+extern void walsh(const long bsize[3], const long dims[DIMS], _Complex float* sens, const long caldims[DIMS], const _Complex float* data);
 
diff --git a/src/scc.c b/src/cc.c
similarity index 53%
rename from src/scc.c
rename to src/cc.c
index 9ec05e58a..2c339d685 100644
--- a/src/scc.c
+++ b/src/cc.c
@@ -1,19 +1,10 @@
-/* Copyright 2013. The Regents of the University of California.
- * All rights reserved. Use of this source code is governed by 
+/* Copyright 2013-2014. The Regents of the University of California.
+ * All rights reserved. Use of this source code is governed by
  * a BSD-style license which can be found in the LICENSE file.
  *
- * Authors: 
- * 2012-2013	Martin Uecker <uecker@eecs.berkeley.edu>
+ * Authors:
+ * 2012-2014	Martin Uecker <uecker@eecs.berkeley.edu>
  * 2013		Jonathan Tamir <jtamir@eecs.berkeley.edu>
- * 
- *
- * Huang F, Vijayakumar S, Li Y, Hertel S, Duensing GR. A software channel
- * compression technique for faster reconstruction with many channels.
- * Magn Reson Imaging 2008; 26:133-141.
- *
- * Buehrer M, Pruessmann KP, Boesiger P, Kozerke S. Array compression for MRI
- * with large coil arrays. Magn Reson Med 2007, 57: 1131–1139. 
- * 
  */
 
 #define _GNU_SOURCE
@@ -24,22 +15,23 @@
 #include <stdio.h>
 #include <unistd.h>
 
-#include "misc/io.h"
 #include "misc/mmio.h"
 #include "misc/mri.h"
 #include "misc/misc.h"
+#include "misc/debug.h"
 
 #include "num/multind.h"
 #include "num/flpmath.h"
-#include "num/lapack.h"
-#include "num/la.h"
+#include "num/fft.h"
+
+#include "calib/cc.h"
 
 
 
 
 static void usage(const char* name, FILE* fd)
 {
-	fprintf(fd, 	"Usage: %s [-v] [-A] [-r cal_size] [-P num_coeffs]"
+	fprintf(fd, 	"Usage: %s [-A] [-r cal_size] [-P num_coeffs]"
 			" <kspace> <coeff>|<proj_kspace>\n", name);
 }
 
@@ -51,23 +43,23 @@ static void help(void)
 		"-P N\tperform compression to N virtual channels\n"
 		"-r S\tsize of calibration region\n"
 		"-A\tuse all data to compute coefficients\n"
-		"-v\tverbose\n"
+		"-S|G|E\ttype: SVD, Geometric, ESPIRiT\n"
 		"-h\thelp\n");
 }
 
 
 
 
-int main_scc(int argc, char* argv[])
+int main_cc(int argc, char* argv[])
 {
-	long calsize[3] = { 24, 24, 24 }; 
+	long calsize[3] = { 24, 24, 24 };
 	bool proj = false;
-	bool verbose = false;
 	long P = 0;
 	bool all = false;
+	enum cc_type { SCC, GCC, ECC } cc_type = SCC;
 
 	int c;
-	while (-1 != (c = getopt(argc, argv, "r:R:p:P:vAh"))) {
+	while (-1 != (c = getopt(argc, argv, "r:R:SGEP:Ah"))) {
 
 		switch (c) {
 
@@ -90,8 +82,16 @@ int main_scc(int argc, char* argv[])
 			P = atoi(optarg);
 			break;
 
-		case 'v':
-			verbose = true;
+		case 'S':
+			cc_type = SCC;
+			break;
+
+		case 'G':
+			cc_type = GCC;
+			break;
+
+		case 'E':
+			cc_type = ECC;
 			break;
 
 		case 'h':
@@ -111,7 +111,7 @@ int main_scc(int argc, char* argv[])
 		usage(argv[0], stderr);
 		exit(1);
 	}
-		
+
 
 	long in_dims[DIMS];
 
@@ -126,13 +126,9 @@ int main_scc(int argc, char* argv[])
 	long out_dims[DIMS] = MD_INIT_ARRAY(DIMS, 1);
 	out_dims[COIL_DIM] = channels;
 	out_dims[MAPS_DIM] = channels;
+	out_dims[READ_DIM] = (SCC == cc_type) ? 1 : in_dims[READ_DIM];
 
-	complex float* out_data = NULL;
-	
-	if (proj)
-		out_data = md_alloc(DIMS, out_dims, CFL_SIZE);
-	else
-		out_data = create_cfl(argv[optind + 1], DIMS, out_dims);
+	complex float* out_data = (proj ? anon_cfl : create_cfl)(argv[optind + 1], DIMS, out_dims);
 
 
 	long caldims[DIMS];
@@ -148,34 +144,23 @@ int main_scc(int argc, char* argv[])
 		cal_data = extract_calib(caldims, calsize, in_dims, in_data, false);
 	}
 
-	
-
-	complex float* tmp = xmalloc(channels * channels * CFL_SIZE);
-	size_t csize = md_calc_size(3, caldims);
-	gram_matrix(channels, (complex float (*)[channels])tmp, csize, (const complex float (*)[csize])cal_data);
+	if (ECC == cc_type)
+		debug_printf(DP_WARN, "Warning: ECC depends on a parameter choice rule for optimal results which is not implemented.\n");
 
-	float vals[channels];
-	eigendecomp(channels, vals, (complex float (*)[])tmp);
-	md_flip(DIMS, out_dims, MAPS_FLAG, out_data, tmp, CFL_SIZE);
 
-	free(tmp);
-
-	if (verbose) {
-
-		printf("Coefficients:");
-
-		for (int i = 0; i < channels; i++)
-			printf(" %.3f", vals[channels - 1 - i] / vals[channels - 1]);
-
-		printf("\n");
+	switch (cc_type) {
+	case SCC: scc(out_dims, out_data, caldims, cal_data); break;
+	case GCC: gcc(out_dims, out_data, caldims, cal_data); break;
+	case ECC: ecc(out_dims, out_data, caldims, cal_data); break;
 	}
 
+	if (!all)
+		md_free(cal_data);
 
 
 	if (proj) {
-	
-		if (verbose)
-			printf("Compressing to %ld virtual coils...\n", P);
+
+		debug_printf(DP_DEBUG1, "Compressing to %ld virtual coils...\n", P);
 
 		long trans_dims[DIMS];
 		md_copy_dims(DIMS, trans_dims, in_dims);
@@ -187,23 +172,50 @@ int main_scc(int argc, char* argv[])
 		md_select_dims(DIMS, ~COIL_FLAG, fake_trans_dims, in_dims);
 		fake_trans_dims[MAPS_DIM] = P;
 
-		md_zmatmulc(DIMS, fake_trans_dims, trans_data, out_dims, out_data, in_dims, in_data);
+		long out2_dims[DIMS];
+		md_copy_dims(DIMS, out2_dims, out_dims);
+		out2_dims[MAPS_DIM] = P;
 
-		unmap_cfl(DIMS, trans_dims, trans_data);
-	}
+		if (SCC != cc_type) {
 
-	printf("Done.\n");
+			complex float* in2_data = anon_cfl(NULL, DIMS, in_dims);
 
-	if (!all)
-		md_free(cal_data);
+			ifftuc(DIMS, in_dims, READ_FLAG, in2_data, in_data);
+
+			unmap_cfl(DIMS, in_dims, in_data);
+			in_data = in2_data;
+
+
+			complex float* out2 = anon_cfl(NULL, DIMS, out2_dims);
+			align_ro(out2_dims, out2, out_data);
+
+			unmap_cfl(DIMS, out_dims, out_data);
+			out_data = out2;
+		}
+
+		md_zmatmulc(DIMS, fake_trans_dims, trans_data, out2_dims, out_data, in_dims, in_data);
+
+		if (SCC != cc_type) {
+
+			fftuc(DIMS, trans_dims, READ_FLAG, trans_data, trans_data);
 
-	unmap_cfl(DIMS, in_dims, in_data);
+			unmap_cfl(DIMS, out2_dims, out_data);
 
-	if (proj)
-		md_free(out_data);
-	else
+		} else {
+
+			unmap_cfl(DIMS, out_dims, out_data);
+		}
+
+		unmap_cfl(DIMS, trans_dims, trans_data);
+		unmap_cfl(DIMS, in_dims, in_data);
+
+	} else {
+
+		unmap_cfl(DIMS, in_dims, in_data);
 		unmap_cfl(DIMS, out_dims, out_data);
+	}
 
+	printf("Done.\n");
 	exit(0);
 }
 
diff --git a/src/cpyphs.c b/src/cpyphs.c
new file mode 100644
index 000000000..fa57ce7a2
--- /dev/null
+++ b/src/cpyphs.c
@@ -0,0 +1,48 @@
+/* Copyright 2014. The Regents of the University of California.
+ * 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 Martin Uecker <uecker@eecs.berkeley.edu>
+ */
+
+#include <stdlib.h>
+#include <assert.h>
+#include <stdio.h>
+#include <complex.h>
+#include <assert.h>
+
+#include "num/multind.h"
+#include "num/flpmath.h"
+
+#include "misc/mmio.h"
+#include "misc/misc.h"
+
+
+
+#ifndef DIMS
+#define DIMS 16
+#endif
+
+
+static const char* usage_str = "<input> <output";
+static const char* help_str = "Copy phase from <input> to <output>.\n";
+
+
+int main_cpyphs(int argc, char* argv[argc])
+{
+	mini_cmdline(argc, argv, 2, usage_str, help_str);
+
+	long dims[DIMS];
+	
+	complex float* in_data = load_cfl(argv[1], DIMS, dims);
+	complex float* out_data = create_cfl(argv[2], DIMS, dims);
+
+	md_zphsr(DIMS, dims, out_data, in_data);
+
+	unmap_cfl(DIMS, dims, in_data);
+	unmap_cfl(DIMS, dims, out_data);
+	exit(0);
+}
+
+
diff --git a/src/ecalib.c b/src/ecalib.c
index f4855a5b7..454d3744d 100644
--- a/src/ecalib.c
+++ b/src/ecalib.c
@@ -26,6 +26,10 @@
 
 #include "calib/calib.h"
 
+#ifndef CFL_SIZE
+#define CFL_SIZE sizeof(complex float)
+#endif
+
 
 
 static void usage(const char* name, FILE* fd)
@@ -190,8 +194,8 @@ int main_ecalib(int argc, char* argv[])
 		for (int i = 3; i < N; i++)
 			cal_dims[i] = ksp_dims[i];
 
-		cal_data = md_alloc(5, cal_dims, sizeof(complex float));
-		md_resizec(5, cal_dims, cal_data, ksp_dims, in_data, sizeof(complex float));
+		cal_data = md_alloc(5, cal_dims, CFL_SIZE);
+		md_resizec(5, cal_dims, cal_data, ksp_dims, in_data, CFL_SIZE);
 	 }
 
 
@@ -233,13 +237,13 @@ int main_ecalib(int argc, char* argv[])
 		long cov_dims[4];
 
 		calone_dims(&conf, cov_dims, channels);
-		complex float* imgcov = md_alloc(4, cov_dims, sizeof(complex float));
+		complex float* imgcov = md_alloc(4, cov_dims, CFL_SIZE);
 
 
 		calone(&conf, cov_dims, imgcov, K, svals, cal_dims, cal_data);
 
 		complex float* out = create_cfl(argv[optind + 1], 4, cov_dims);
-		md_copy(4, cov_dims, out, imgcov, sizeof(complex float));
+		md_copy(4, cov_dims, out, imgcov, CFL_SIZE);
 		unmap_cfl(4, cov_dims, out);
 
 //		caltwo(crthr, out_dims, out_data, emaps, cov_dims, imgcov, NULL, NULL);
@@ -272,25 +276,18 @@ int main_ecalib(int argc, char* argv[])
 		map_dims[COIL_DIM] = 1;
 		map_dims[MAPS_DIM] = maps;
 
+		const char* emaps_file = NULL;
 
-		complex float* out_data = create_cfl(argv[optind + 1], N, out_dims);
-		complex float* emaps;
-
-		if (argc - optind == 3) 
-			emaps = create_cfl(argv[optind + 2], N, map_dims);
-		else
-			emaps = md_alloc(N, map_dims, sizeof(complex float));
-
+		if (3 == argc - optind)
+			emaps_file = argv[optind + 2];
 
+		complex float* out_data = create_cfl(argv[optind + 1], N, out_dims);
+		complex float* emaps = (emaps_file ? create_cfl : anon_cfl)(emaps_file, N, map_dims);
 
 		calib(&conf, out_dims, out_data, emaps, K, svals, cal_dims, cal_data); 
 
 		unmap_cfl(N, out_dims, out_data);
-
-		if (argc - optind == 3) 
-			unmap_cfl(N, map_dims, emaps);
-		else
-			md_free(emaps);
+		unmap_cfl(N, map_dims, emaps);
 	}
 
 
diff --git a/src/ecaltwo.c b/src/ecaltwo.c
index dd0f188d8..cb0ff869c 100644
--- a/src/ecaltwo.c
+++ b/src/ecaltwo.c
@@ -21,6 +21,10 @@
 #include "misc/mmio.h"
 #include "misc/mri.h"
 
+#ifndef CFL_SIZE
+#define CFL_SIZE sizeof(complex float)
+#endif
+
 
 
 static void usage(const char* name, FILE* fp)
@@ -89,9 +93,9 @@ int main_ecaltwo(int argc, char* argv[])
 		exit(1);
 	}
 
-	long in_dims[KSPACE_DIMS];
+	long in_dims[DIMS];
 
-	complex float* in_data = load_cfl(argv[optind + 3], KSPACE_DIMS, in_dims);
+	complex float* in_data = load_cfl(argv[optind + 3], DIMS, in_dims);
 
 	int channels = 0;
 
@@ -103,8 +107,8 @@ int main_ecaltwo(int argc, char* argv[])
 	assert(maps <= channels);
 
 
-	long out_dims[KSPACE_DIMS] = { [0 ... KSPACE_DIMS - 1] = 1 };
-	long map_dims[KSPACE_DIMS] = { [0 ... KSPACE_DIMS - 1] = 1 };
+	long out_dims[DIMS] = { [0 ... DIMS - 1] = 1 };
+	long map_dims[DIMS] = { [0 ... DIMS - 1] = 1 };
 	
 	out_dims[0] = atoi(argv[optind + 0]);
 	out_dims[1] = atoi(argv[optind + 1]);
@@ -124,23 +128,23 @@ int main_ecaltwo(int argc, char* argv[])
 	map_dims[4] = maps;
 
 
-	complex float* out_data = create_cfl(argv[optind + 4], KSPACE_DIMS, out_dims);
+	complex float* out_data = create_cfl(argv[optind + 4], DIMS, out_dims);
 	complex float* emaps;
 
 	if (6 == argc - optind)
-		emaps = create_cfl(argv[optind + 5], KSPACE_DIMS, map_dims);
+		emaps = create_cfl(argv[optind + 5], DIMS, map_dims);
 	else
-		emaps = md_alloc(KSPACE_DIMS, map_dims, sizeof(complex float));
+		emaps = md_alloc(KSPACE_DIMS, map_dims, CFL_SIZE);
 
 	caltwo(&conf, out_dims, out_data, emaps, in_dims, in_data, NULL, NULL);
 
 	printf("Done.\n");
 
-	unmap_cfl(KSPACE_DIMS, in_dims, in_data);
-	unmap_cfl(KSPACE_DIMS, out_dims, out_data);
+	unmap_cfl(DIMS, in_dims, in_data);
+	unmap_cfl(DIMS, out_dims, out_data);
 
 	if (6 == argc - optind)
-		unmap_cfl(KSPACE_DIMS, map_dims, emaps);
+		unmap_cfl(DIMS, map_dims, emaps);
 	else
 		md_free(emaps);
 
diff --git a/src/fmac.c b/src/fmac.c
index fd7a014c8..24de9b8ef 100644
--- a/src/fmac.c
+++ b/src/fmac.c
@@ -104,16 +104,16 @@ int main_fmac(int argc, char* argv[])
 
 	if (clear) {
 
-		md_clear(N, dimso, out, sizeof(complex float));
+		md_clear(N, dimso, out, CFL_SIZE);
 	}
 
 	long str1[N];
 	long str2[N];
 	long stro[N];
 
-	md_calc_strides(N, str1, dims1, sizeof(complex float));
-	md_calc_strides(N, str2, dims2, sizeof(complex float));
-	md_calc_strides(N, stro, dimso, sizeof(complex float));
+	md_calc_strides(N, str1, dims1, CFL_SIZE);
+	md_calc_strides(N, str2, dims2, CFL_SIZE);
+	md_calc_strides(N, stro, dimso, CFL_SIZE);
 
 	(conj ? md_zfmacc2 : md_zfmac2)(N, dims, stro, out, str1, data1, str2, data2);
 
diff --git a/src/iter/admm.c b/src/iter/admm.c
index 2e4e7742a..b1c5a481b 100644
--- a/src/iter/admm.c
+++ b/src/iter/admm.c
@@ -95,7 +95,7 @@ void admm(struct admm_history_s* history, const struct admm_plan_s* plan,
 	  void* obj_eval_data,
 	  float (*obj_eval)(const void*, const float*))
 {
-	_Bool fast = plan->fast;
+	bool fast = plan->fast;
 	double ABSTOL = plan->ABSTOL;
 	double RELTOL = plan->RELTOL;
 	float tau = plan->tau;
@@ -209,7 +209,6 @@ void admm(struct admm_history_s* history, const struct admm_plan_s* plan,
 		// update x
 		vops->clear(N, rhs);
 		vops->sub(M, r, z, u);
-
 		for (unsigned int j = 0; j < num_funs; j++) {
 			pos = md_calc_offset(j, fake_strs, z_dims);
 			plan->ops[j].adjoint(plan->ops[j].data, s, r + pos);
diff --git a/src/iter/italgos.c b/src/iter/italgos.c
index d767cb62e..7db6da168 100644
--- a/src/iter/italgos.c
+++ b/src/iter/italgos.c
@@ -156,8 +156,8 @@ static float ist_continuation(struct iter_data* itrdata, const float delta)
  * @param x initial estimate
  * @param b observations
  */
-void ist(unsigned int maxiter, float epsilon, float tau, 
-		float continuation, _Bool hogwild, long N, void* data,
+void ist(unsigned int maxiter, float epsilon, float tau,
+		float continuation, bool hogwild, long N, void* data,
 		const struct vec_iter_s* vops,
 		void (*op)(void* data, float* dst, const float* src), 
 		void (*thresh)(void* data, float lambda, float* dst, const float* src),
@@ -264,7 +264,7 @@ void ist(unsigned int maxiter, float epsilon, float tau,
  * @param b observations
  */
 void fista(unsigned int maxiter, float epsilon, float tau, 
-	   float continuation, _Bool hogwild, 
+	   float continuation, bool hogwild,
 	   long N, void* data,
 	   const struct vec_iter_s* vops,
 	   void (*op)(void* data, float* dst, const float* src), 
diff --git a/src/iter/prox.c b/src/iter/prox.c
index ee53ed610..68159c68b 100755
--- a/src/iter/prox.c
+++ b/src/iter/prox.c
@@ -408,7 +408,7 @@ const struct operator_p_s* prox_lineq_create(const struct linop_s* op, const com
 	pdata->op = op;
 
 	pdata->adj = md_alloc_sameplace(N, dims, CFL_SIZE, y);
-	linop_adjoint(op, N, dims, pdata->adj, N, dims, y);
+	linop_adjoint(op, N, dims, pdata->adj, N, linop_codomain(op)->dims, y);
 
 	pdata->tmp = md_alloc_sameplace(N, dims, CFL_SIZE, y);
 
diff --git a/src/iter/thresh.c b/src/iter/thresh.c
index 3ad9ee410..6b3ee2501 100755
--- a/src/iter/thresh.c
+++ b/src/iter/thresh.c
@@ -153,7 +153,7 @@ const struct operator_p_s* prox_thresh_create(unsigned int D, const long dim[D],
  * @param flags bitmask for joint soft-thresholding
  * @param gpu true if using gpu, false if using cpu
  */
-extern const struct operator_p_s* prox_unithresh_create(unsigned int D, const struct linop_s* unitary_op, const float lambda, const unsigned long flags, _Bool gpu)
+extern const struct operator_p_s* prox_unithresh_create(unsigned int D, const struct linop_s* unitary_op, const float lambda, const unsigned long flags, bool gpu)
 {
 	struct thresh_s* data = xmalloc(sizeof(struct thresh_s));
 
diff --git a/src/linops/sampling.c b/src/linops/sampling.c
index 5d27383e3..6bb5c6137 100644
--- a/src/linops/sampling.c
+++ b/src/linops/sampling.c
@@ -10,6 +10,7 @@
 
 #include "misc/mri.h"
 #include "misc/misc.h"
+#include "misc/debug.h"
 
 #include "num/flpmath.h"
 #include "num/multind.h"
@@ -46,7 +47,7 @@ struct linop_s* sampling_create(const long dims[DIMS], const long pat_dims[DIMS]
 	struct sampling_data_s* data = xmalloc(sizeof(struct sampling_data_s));
 
 	md_select_dims(DIMS, ~MAPS_FLAG, data->dims, dims); // dimensions of kspace
-	md_calc_strides(DIMS, data->strs, dims, CFL_SIZE);
+	md_calc_strides(DIMS, data->strs, data->dims, CFL_SIZE);
 	md_calc_strides(DIMS, data->pat_strs, pat_dims, CFL_SIZE);
 
 	data->pattern = pattern;
diff --git a/src/linops/ufft.c b/src/linops/ufft.c
index e0d448d86..eb58e893e 100644
--- a/src/linops/ufft.c
+++ b/src/linops/ufft.c
@@ -36,7 +36,8 @@
  * data structure for holding the undersampled fft data.
  */
 struct ufft_data {
-	_Bool use_gpu;
+
+	bool use_gpu;
 	unsigned int flags;
 
 	long ksp_dims[DIMS];
@@ -50,7 +51,7 @@ struct ufft_data {
         complex float* pat;
 };
 
-static struct ufft_data* ufft_create_data(const long ksp_dims[DIMS], const long pat_dims[DIMS], const complex float* pat, unsigned int flags, _Bool use_gpu);
+static struct ufft_data* ufft_create_data(const long ksp_dims[DIMS], const long pat_dims[DIMS], const complex float* pat, unsigned int flags, bool use_gpu);
 static void ufft_free_data(const void* _data );
 static void ufft_apply(const void* _data, complex float* dst, const complex float* src);
 static void ufft_apply_adjoint(const void* _data, complex float* dst, const complex float* src);
@@ -60,7 +61,7 @@ static void ufft_apply_pinverse(const void* _data, float rho, complex float* dst
 /**
  * Create undersampled/weighted fft operator
  */
-const struct linop_s* ufft_create(const long ksp_dims[DIMS], const long pat_dims[DIMS], const complex float* pat, unsigned int flags, _Bool use_gpu)
+const struct linop_s* ufft_create(const long ksp_dims[DIMS], const long pat_dims[DIMS], const complex float* pat, unsigned int flags, bool use_gpu)
 {
 	struct ufft_data* data = ufft_create_data( ksp_dims, pat_dims, pat, flags, use_gpu );
 
@@ -70,7 +71,7 @@ const struct linop_s* ufft_create(const long ksp_dims[DIMS], const long pat_dims
 }
 
 
-static struct ufft_data* ufft_create_data(const long ksp_dims[DIMS], const long pat_dims[DIMS], const complex float* pat, unsigned int flags, _Bool use_gpu)
+static struct ufft_data* ufft_create_data(const long ksp_dims[DIMS], const long pat_dims[DIMS], const complex float* pat, unsigned int flags, bool use_gpu)
 {
 	struct ufft_data* data = xmalloc(sizeof(struct ufft_data));
 
diff --git a/src/misc/mmio.c b/src/misc/mmio.c
index 036a00e80..f9f9be6fc 100644
--- a/src/misc/mmio.c
+++ b/src/misc/mmio.c
@@ -20,6 +20,8 @@
 #include <sys/mman.h>
 
 #include "num/multind.h"
+
+#include "misc/misc.h"
 #include "misc/io.h"
 
 #include "mmio.h"
@@ -68,6 +70,8 @@ float* create_coo(const char* name, unsigned int D, const long dims[D])
 
 
 
+
+
 complex float* create_zcoo(const char* name, unsigned int D, const long dimensions[D])
 {
 	long dims[D + 1];
@@ -111,6 +115,7 @@ complex float* create_cfl(const char* name, unsigned int D, const long dimension
 
 
 
+
 float* load_coo(const char* name, unsigned int D, long dims[D])
 {
 	int fd;
@@ -217,6 +222,20 @@ complex float* shared_cfl(unsigned int D, const long dims[D], const char* name)
 }
 
 
+complex float* anon_cfl(const char* name, unsigned int D, const long dims[D])
+{
+	UNUSED(name);
+	void* addr;
+	long T = md_calc_size(D, dims) * sizeof(complex float);
+
+	if (MAP_FAILED == (addr = mmap(NULL, T, PROT_READ|PROT_WRITE, MAP_ANONYMOUS|MAP_PRIVATE, -1, 0)))
+		abort();
+
+	return (complex float*)addr;
+}
+
+
+
 #if 0
 void* private_raw(size_t* size, const char* name)
 {
diff --git a/src/misc/mmio.h b/src/misc/mmio.h
index 9b9a149c3..7e59b7876 100644
--- a/src/misc/mmio.h
+++ b/src/misc/mmio.h
@@ -20,6 +20,7 @@ extern _Complex float* shared_cfl(unsigned int D, const long dims[__VLA(D)], con
 extern _Complex float* private_cfl(unsigned int D, const long dims[__VLA(D)], const char* name);
 extern void unmap_cfl(unsigned int D, const long dims[__VLA(D)], const _Complex float* x);
 
+extern _Complex float* anon_cfl(const char* name, unsigned int D, const long dims[__VLA(D)]);
 extern _Complex float* create_cfl(const char* name, unsigned int D, const long dimensions[__VLA(D)]);
 extern _Complex float* load_cfl(const char* name, unsigned int D, long dimensions[__VLA(D)]);
 
diff --git a/src/misc/mri.h b/src/misc/mri.h
index b2a07c28d..670dc9204 100644
--- a/src/misc/mri.h
+++ b/src/misc/mri.h
@@ -36,7 +36,9 @@ enum mri_dims {
 	SLICE_DIM,
 };
 
+#ifdef BERKELEY_SVN
 #define KSPACE_DIMS	16u
+#endif
 
 #ifndef DIMS
 #define DIMS		16u
diff --git a/src/noncart/nufft.c b/src/noncart/nufft.c
index 970b658fe..cb72bf5b0 100644
--- a/src/noncart/nufft.c
+++ b/src/noncart/nufft.c
@@ -45,16 +45,16 @@
  *
  */
 struct nufft_data {
+
 	void* cgconf;
 
-	_Bool use_gpu;
+	bool use_gpu;
 
 	float scale;
 
 	int nshifts;
 	float* shifts;
 	complex float* linphase;
-	
 
 	const complex float* traj;
 	const complex float* pat;
@@ -87,8 +87,8 @@ struct nufft_data {
 	long img_size;
 	long coilim_size;
 
-	_Bool toeplitz;
-	_Bool precond;
+	bool toeplitz;
+	bool precond;
 
 
 	float rho;
@@ -96,7 +96,7 @@ struct nufft_data {
 
 
 // Forward: from image to kspace
-static struct nufft_data* nufft_create_data( const long ksp_dims[DIMS], const long coilim_dims[DIMS], const complex float* traj, const complex float* pat, _Bool toeplitz, _Bool precond, void* cgconf, _Bool use_gpu);
+static struct nufft_data* nufft_create_data( const long ksp_dims[DIMS], const long coilim_dims[DIMS], const complex float* traj, const complex float* pat, bool toeplitz, bool precond, void* cgconf, bool use_gpu);
 static void nufft_free_data( const void* data );
 static void nufft_apply(const void* _data, complex float* dst, const complex float* src);
 static void nufft_apply_adjoint(const void* _data, complex float* dst, const complex float* src);
@@ -126,7 +126,7 @@ static void fill_linphases( struct nufft_data* data );
  * @param use_gpu       -     use gpu boolean
  *
  */
-struct linop_s* nufft_create( const long ksp_dims[DIMS], const long coilim_dims[DIMS], const complex float* traj, const complex float* pat, _Bool toeplitz, _Bool precond, void* cgconf, _Bool use_gpu)
+struct linop_s* nufft_create( const long ksp_dims[DIMS], const long coilim_dims[DIMS], const complex float* traj, const complex float* pat, bool toeplitz, bool precond, void* cgconf, bool use_gpu)
 {
 	struct nufft_data* data = nufft_create_data( ksp_dims, coilim_dims, traj, pat, toeplitz, precond, cgconf, use_gpu);
 
@@ -138,7 +138,7 @@ struct linop_s* nufft_create( const long ksp_dims[DIMS], const long coilim_dims[
 
 
 
-static struct nufft_data* nufft_create_data( const long ksp_dims[DIMS], const long coilim_dims[DIMS], const complex float* traj, const complex float* pat, _Bool toeplitz, _Bool precond, void* cgconf, _Bool use_gpu)
+static struct nufft_data* nufft_create_data( const long ksp_dims[DIMS], const long coilim_dims[DIMS], const complex float* traj, const complex float* pat, bool toeplitz, bool precond, void* cgconf, bool use_gpu)
 {
 	
 	struct nufft_data* data = (struct nufft_data*) xmalloc( sizeof( struct nufft_data ) );
@@ -212,11 +212,13 @@ static struct nufft_data* nufft_create_data( const long ksp_dims[DIMS], const lo
 	data->toeplitz = toeplitz;
 	data->precond = precond && toeplitz;
 
+	data->psf = md_alloc( DIMS, linphase_dims, CFL_SIZE );
+
 	if (toeplitz)
 	{
-		data->psf = md_alloc( DIMS, linphase_dims, CFL_SIZE );
 		fill_psf( data );
 	}
+	dump_cfl("psf", DIMS, linphase_dims, data->psf );
 
 
 	return data;
diff --git a/src/nufft.c b/src/nufft.c
index 7ff2e9c43..e0980e7df 100644
--- a/src/nufft.c
+++ b/src/nufft.c
@@ -78,7 +78,7 @@ int main_nufft(int argc, char* argv[])
 	md_singleton_dims(DIMS, coilim_dims);
 
 	int maxiter = 50;
-	float lambda = 0.001;
+	float lambda = 0.00;
 
 	const char* pat_str = NULL;
 
diff --git a/src/num/flpmath.c b/src/num/flpmath.c
index 91dc7f61d..d45708ebe 100644
--- a/src/num/flpmath.c
+++ b/src/num/flpmath.c
@@ -1167,6 +1167,10 @@ static void md_zmatmul2_priv(unsigned int D, const long out_dims[D], const long
 	long max_dims[D];
 	md_mmdims(D, max_dims, in_dims, out_dims);
 
+	long max2_dims[D];
+	md_mmdims(D, max2_dims, mat_dims, out_dims);
+	assert(md_check_compat(D, 0, max_dims, max2_dims));
+
 	md_clear2(D, out_dims, out_strs, dst, CFL_SIZE);
 	(conj ? md_zfmacc2 : md_zfmac2)(D, max_dims, out_strs, dst, in_strs, src, mat_strs, mat);
 }
diff --git a/src/num/lapack.c b/src/num/lapack.c
index e373eae9d..db42adf8f 100644
--- a/src/num/lapack.c
+++ b/src/num/lapack.c
@@ -35,12 +35,19 @@
 #include "num/lapack.h"
 
 
+/* ATTENTION: blas and lapack use column-major matrices
+ * while native C uses row-major. All matrices are
+ * transposed to what one would expect.
+  **/
+
+
 #ifdef USE_ACML
 #if 1
+// FIXME: check indices
 extern void cheev(char jobz, char uplo, long N, complex float a[N][N], long lda, float w[N], long* info);
 extern void zheev(char jobz, char uplo, long N, complex double a[N][N], long lda, double w[N], long* info);
-extern void cgesdd(const char jobz, long M, long N, complex float A[M][N], long lda, float* S, complex float U[M][N], long ldu, complex float VT[M][N], long ldvt, const long* info);
-extern void zgesdd(const char jobz, long M, long N, complex double A[M][N], long lda, double* S, complex double U[M][N], long ldu, complex double VT[M][N], long ldvt, const long* info);
+extern void cgesdd(const char jobz, long M, long N, complex float A[M][N], long lda, float* S, complex float U[M][N], long ldu, complex float VH[M][N], long ldvt, const long* info);
+extern void zgesdd(const char jobz, long M, long N, complex double A[M][N], long lda, double* S, complex double U[M][N], long ldu, complex double VH[M][N], long ldvt, const long* info);
 extern void cgesvd(char jobu, char jobvt, long M, long N, complex float a[M][N], long lda, float* S, complex float u[M][N], long ldu, complex float vt[M][N], long ldvt, long *info);
 extern void cgemm(const char transa, const char transb, long M, long N,  long K, const complex float* alpha, const complex float A[M][K], const long lda, const complex float B[K][N], const long ldb, const complex float* beta, complex float C[M][N], const long ldc );
 #else
@@ -52,9 +59,9 @@ extern void cgemm(const char transa, const char transb, long M, long N,  long K,
 #else
 extern void cheev_(const char jobz[1], const char uplo[1], const long* N, complex float a[*N][*N], const long* lda, float w[*N], complex float* work, const long* lwork, float* rwork, long* info);
 extern void zheev_(const char jobz[1], const char uplo[1], const long* N, complex double a[*N][*N], const long* lda, double w[*N], complex double* work, const long* lwork, double* rwork, long* info);
-extern void cgesdd_(const char jobz[1], const long* M, const long* N, complex float A[*M][*N], const long lda[1], float* S, complex float U[*M][*N], const long* ldu, complex float VT[*M][*N], const long* ldvt, complex float* work, const long* lwork, float* rwork, const long* iwork, const long* info);
-extern void zgesdd_(const char jobz[1], const long* M, const long* N, complex double A[*M][*N], const long lda[1], double* S, complex double U[*M][*N], const long* ldu, complex double VT[*M][*N], const long* ldvt, complex double* work, const long* lwork, double* rwork, const long* iwork, const long* info);
-extern void cgesvd_(const char jobu[1], const char jobvt[1], const long* M, const long* N, complex float A[*M][*N], const long* lda, float* s, complex float U[*M][*N], long* ldu, complex float VT[*M][*N], long* ldvt, complex float* work, long* lwork, float* rwork, const long* iwork, long* info);
+extern void cgesdd_(const char jobz[1], const long* M, const long* N, complex float A[*M][*N], const long lda[1], float* S, complex float U[*M][*N], const long* ldu, complex float VH[*M][*N], const long* ldvt, complex float* work, const long* lwork, float* rwork, const long* iwork, const long* info);
+extern void zgesdd_(const char jobz[1], const long* M, const long* N, complex double A[*M][*N], const long lda[1], double* S, complex double U[*M][*N], const long* ldu, complex double VH[*M][*N], const long* ldvt, complex double* work, const long* lwork, double* rwork, const long* iwork, const long* info);
+extern void cgesvd_(const char jobu[1], const char jobvt[1], const long* M, const long* N, complex float A[*M][*N], const long* lda, float* s, complex float U[*M][*N], long* ldu, complex float VH[*M][*N], long* ldvt, complex float* work, long* lwork, float* rwork, const long* iwork, long* info);
 extern void cgemm_(const char transa[1], const char transb[1], const long* M, const long* N, const long* K, const complex float* alpha, const complex float A[*M][*K], const long* lda, const complex float B[*K][*N], const long* ldb, const complex float* beta, complex float C[*M][*N], const long* ldc );
 #endif
 
@@ -127,7 +134,7 @@ void eigendecomp(long N, float eigenval[N], complex float matrix[N][N])
 
 
 
-void svd(long M, long N, complex float U[M][M], complex float VT[N][N], float S[(N > M) ? M : N], complex float A[M][N])
+void svd(long M, long N, complex float U[M][M], complex float VH[N][N], float S[(N > M) ? M : N], complex float A[M][N])
 {
         long info = 0;
 	//assert(M >= N);
@@ -135,27 +142,27 @@ void svd(long M, long N, complex float U[M][M], complex float VT[N][N], float S[
 #ifdef USE_CUDA
 #ifdef USE_CULA
 	if (cuda_ondevice( A )) {
-		culaDeviceCgesvd( 'A', 'A', M, N, (culaDeviceFloatComplex *)A, M, (culaDeviceFloat *)S, (culaDeviceFloatComplex *)U, M, (culaDeviceFloatComplex *)VT, N );
+		culaDeviceCgesvd( 'A', 'A', M, N, (culaDeviceFloatComplex *)A, M, (culaDeviceFloat *)S, (culaDeviceFloatComplex *)U, M, (culaDeviceFloatComplex *)VH, N );
 	} else
 #endif 
 #endif 
 	{
 #ifdef USE_ACML
-	cgesdd('A', M, N, A, M, S, U, M, VT, N, &info);
+	cgesdd('A', M, N, A, M, S, U, M, VH, N, &info);
 #else
 	long lwork = -1;
 	complex float work1[1];
 	float* rwork = xmalloc(MIN(M, N) * MAX(5 * MIN(M, N) + 7, 2 * MAX(M, N) + 2 * MIN(M, N) + 1) * sizeof(float));
 	long* iwork = xmalloc(8 * MIN(M, N) * sizeof(long));
 
-	cgesdd_("A", &M, &N, A, &M, S, U, &M, VT, &N, work1, &lwork, rwork, iwork, &info);
+	cgesdd_("A", &M, &N, A, &M, S, U, &M, VH, &N, work1, &lwork, rwork, iwork, &info);
 
 	if (0 != info)
 		goto err;
 
 	lwork = (int)work1[0];
 	complex float* work = xmalloc(MAX(1, lwork) * sizeof(complex float));
-	cgesdd_("A", &M, &N, A, &M, S, U, &M, VT, &N, work, &lwork, rwork, iwork, &info);
+	cgesdd_("A", &M, &N, A, &M, S, U, &M, VH, &N, work, &lwork, rwork, iwork, &info);
 	free(rwork);
 	free(iwork);
 	free(work);
@@ -174,7 +181,7 @@ void svd(long M, long N, complex float U[M][M], complex float VT[N][N], float S[
 
 void svd_econ(long M, long N, 
 	      complex float U[M][(N > M) ? M : N], 
-	      complex float VT[(N > M) ? M : N][N], 
+	      complex float VH[(N > M) ? M : N][N],
 	      float S[(N > M) ? M : N], complex float A[M][N])
 {
 	long info = 0;
@@ -184,7 +191,7 @@ void svd_econ(long M, long N,
 #ifdef USE_CUDA
 #ifdef USE_CULA
 	if (cuda_ondevice( A )) {
-		culaDeviceCgesvd( 'S', 'S', M, N, (culaDeviceFloatComplex *)A, M, (culaDeviceFloat *)S, (culaDeviceFloatComplex *)U, M, (culaDeviceFloatComplex *)VT, minMN );
+		culaDeviceCgesvd( 'S', 'S', M, N, (culaDeviceFloatComplex *)A, M, (culaDeviceFloat *)S, (culaDeviceFloatComplex *)U, M, (culaDeviceFloatComplex *)VH, minMN );
 	} else
 #endif 
 #endif 
@@ -192,21 +199,21 @@ void svd_econ(long M, long N,
 	{
 
 #ifdef USE_ACML
-	cgesvd('S', 'S', M, N, A, M, S, U, M, VT, minMN, &info);
+	cgesvd('S', 'S', M, N, A, M, S, U, M, VH, minMN, &info);
 #else
 	long lwork = -1;
 	complex float work1[1];
 	float* rwork = xmalloc(5 * N * sizeof(float));
 	long* iwork = xmalloc(8 * minMN * sizeof(long));
 
-	cgesvd_("S", "S", &M, &N, A, &M, S, U, &M, VT, &minMN, work1, &lwork, rwork, iwork, &info);
+	cgesvd_("S", "S", &M, &N, A, &M, S, U, &M, VH, &minMN, work1, &lwork, rwork, iwork, &info);
 
 	if(0 != info)
 		goto err;
 
 	lwork = (int)work1[0];
 	complex float* work = xmalloc(lwork * sizeof(complex float));
-	cgesvd_("S", "S", &M, &N, A, &M, S, U, &M, VT, &minMN, work, &lwork, rwork, iwork, &info);
+	cgesvd_("S", "S", &M, &N, A, &M, S, U, &M, VH, &minMN, work, &lwork, rwork, iwork, &info);
 
 	free(work);
 	free(iwork);
@@ -224,27 +231,27 @@ void svd_econ(long M, long N,
 }
 
 
-void svd_double(long M, long N, complex double U[M][M], complex double VT[N][N], double S[(N > M) ? M : N], complex double A[M][N])
+void svd_double(long M, long N, complex double U[M][M], complex double VH[N][N], double S[(N > M) ? M : N], complex double A[M][N])
 {
         long info = 0;
 	//assert(M >= N);
 	
 #ifdef USE_ACML
-	zgesdd('A', M, N, A, M, S, U, M, VT, N, &info);
+	zgesdd('A', M, N, A, M, S, U, M, VH, N, &info);
 #else
 	long lwork = -1;
 	complex double work1[1];
 	double* rwork = xmalloc(MIN(M, N) * MAX(5 * MIN(M, N) + 7, 2 * MAX(M, N) + 2 * MIN(M, N) + 1) * sizeof(double));
 	long* iwork = xmalloc(8 * MIN(M, N) * sizeof(long));
 
-	zgesdd_("A", &M, &N, A, &M, S, U, &M, VT, &N, work1, &lwork, rwork, iwork, &info);
+	zgesdd_("A", &M, &N, A, &M, S, U, &M, VH, &N, work1, &lwork, rwork, iwork, &info);
 
 	if (0 != info)
 		goto err;
 
 	lwork = (int)work1[0];
 	complex double* work = xmalloc(MAX(1, lwork) * sizeof(complex double));
-	zgesdd_("A", &M, &N, A, &M, S, U, &M, VT, &N, work, &lwork, rwork, iwork, &info);
+	zgesdd_("A", &M, &N, A, &M, S, U, &M, VH, &N, work, &lwork, rwork, iwork, &info);
 	free(rwork);
 	free(iwork);
 	free(work);
diff --git a/src/num/lapack.h b/src/num/lapack.h
index 8ee689431..d7f63b632 100644
--- a/src/num/lapack.h
+++ b/src/num/lapack.h
@@ -6,13 +6,14 @@
 #include <complex.h>
 
 extern void eigendecomp(long N, float eigenval[N], complex float matrix[N][N]);
-extern void svd(long M, long N, complex float U[M][M], complex float VT[N][N], float S[(N > M) ? M : N], complex float A[M][N]);
+extern void svd(long M, long N, complex float U[M][M], complex float VH[N][N], float S[(N > M) ? M : N], complex float A[N][M]);
 extern void svd_econ(long M, long N, 
 		     complex float U[M][(N > M) ? M : N], 
-		     complex float VT[(N > M) ? M : N][N], 
-		     float S[(N > M) ? M : N], complex float A[M][N]);
+		     complex float VH[(N > M) ? M : N][N],
+		     float S[(N > M) ? M : N],
+		     complex float A[N][M]);
 
 extern void eigendecomp_double(long N, double eigenval[N], complex double matrix[N][N]);
-extern void svd_double(long M, long N, complex double U[M][M], complex double VT[N][N], double S[(N > M) ? M : N], complex double A[M][N]);
+extern void svd_double(long M, long N, complex double U[M][M], complex double VH[N][N], double S[(N > M) ? M : N], complex double A[N][M]);
 extern void matrix_multiply(long M, long N, long K, complex float C[M][N], const complex float A[M][K], const complex float B[K][N]);
 extern void cgemm_sameplace(const char transa, const char transb, long M, long N,  long K, const complex float* alpha, const complex float A[M][K], const long lda, const complex float B[K][N], const long ldb, const complex float* beta, complex float C[M][N], const long ldc );
diff --git a/src/num/multind.c b/src/num/multind.c
index f97488368..eb23793a8 100644
--- a/src/num/multind.c
+++ b/src/num/multind.c
@@ -212,7 +212,7 @@ unsigned int md_calc_blockdim(unsigned int D, const long dim[D], const long str[
  * @param odims output dimensions
  * @param idims input dimensions
  */
-extern void md_select_dims(unsigned int D, unsigned long flags, long odims[D], const long idims[D])
+void md_select_dims(unsigned int D, unsigned long flags, long odims[D], const long idims[D])
 {
 	md_copy_dims(D, odims, idims);
 	
@@ -228,7 +228,7 @@ extern void md_select_dims(unsigned int D, unsigned long flags, long odims[D], c
  *
  * odims[i] = idims[i]
  */
-extern void md_copy_dims(unsigned int D, long odims[D], const long idims[D])
+void md_copy_dims(unsigned int D, long odims[D], const long idims[D])
 {
 	memcpy(odims, idims, D  * sizeof(long));
 }
@@ -240,7 +240,7 @@ extern void md_copy_dims(unsigned int D, long odims[D], const long idims[D])
  *
  * ostrs[i] = istrs[i]
  */
-extern void md_copy_strides(unsigned int D, long ostrs[D], const long istrs[D])
+void md_copy_strides(unsigned int D, long ostrs[D], const long istrs[D])
 {
 	memcpy(ostrs, istrs, D  * sizeof(long));
 }
@@ -252,7 +252,7 @@ extern void md_copy_strides(unsigned int D, long ostrs[D], const long istrs[D])
  *
  * dims[i] = val
  */
-extern void md_set_dims(unsigned int D, long dims[D], long val)
+void md_set_dims(unsigned int D, long dims[D], long val)
 {
 	for (unsigned int i = 0; i < D; i++)
 		dims[i] = val;
@@ -263,7 +263,7 @@ extern void md_set_dims(unsigned int D, long dims[D], long val)
 /**
  * returns whether or not @param pos is a valid index of an array of dimension @param dims
  */
-extern bool md_is_index(unsigned int D, const long pos[D], const long dims[D])
+bool md_is_index(unsigned int D, const long pos[D], const long dims[D])
 {
 	if (D == 0)
 		return true;
@@ -273,12 +273,24 @@ extern bool md_is_index(unsigned int D, const long pos[D], const long dims[D])
 
 
 
+/**
+ * return whether some other dimensions are >1
+ */
+bool md_check_dimensions(unsigned int N, const long dims[N], unsigned int flags)
+{
+	long d[N];
+	md_select_dims(N, ~flags, d, dims);
+	return (1 != md_calc_size(N, d));
+}
+
+
+
 /**
  * Set all dimensions to one
  *
  * dims[i] = 1
  */
-extern void md_singleton_dims(unsigned int D, long dims[D])
+void md_singleton_dims(unsigned int D, long dims[D])
 {
 	for (unsigned int i = 0; i < D; i++)
 		dims[i] = 1;
@@ -291,7 +303,7 @@ extern void md_singleton_dims(unsigned int D, long dims[D])
  *
  * dims[i] = 1
  */
-extern void md_singleton_strides(unsigned int D, long strs[D])
+void md_singleton_strides(unsigned int D, long strs[D])
 {
 	for (unsigned int i = 0; i < D; i++)
 		strs[i] = 0;
@@ -328,7 +340,7 @@ bool md_check_compat(unsigned int D, unsigned long flags, const long dim1[D], co
  * @param idims1 input 1 dimensions
  * @param idims2 input 2 dimensions
  */
-extern void md_min_dims(unsigned int D, unsigned long flags, long odims[D], const long idims1[D], const long idims2[D])
+void md_min_dims(unsigned int D, unsigned long flags, long odims[D], const long idims1[D], const long idims2[D])
 {
 	for (unsigned int i = 0; i < D; i++)
 		if ((flags & (1 << i)))
@@ -670,6 +682,42 @@ void md_swap(unsigned int D, const long dim[D], void* optr, void* iptr, size_t s
 
 
 
+/**
+ * Move a block from an array to another array (with strides)
+ *
+ */
+void md_move_block2(unsigned int D, const long dim[D], const long opos[D], const long odim[D], const long ostr[D], void* optr, const long ipos[D], const long idim[D], const long istr[D], const void* iptr, size_t size)
+{
+	for (unsigned int i = 0; i < D; i++) {
+
+		assert(dim[i] <= odim[i]);
+		assert(dim[i] <= idim[i]);
+		assert((0 <= opos[i]) && (opos[i] <= odim[i] - dim[i]));
+		assert((0 <= ipos[i]) && (ipos[i] <= idim[i] - dim[i]));
+	}
+
+	long ioff = md_calc_offset(D, istr, ipos);
+	long ooff = md_calc_offset(D, ostr, opos);
+
+	md_copy2(D, dim, ostr, optr + ooff, istr, iptr + ioff, size);
+}
+
+
+/**
+ * Move a block from an array to another array (without strides)
+ *
+ */
+void md_move_block(unsigned int D, const long dim[D], const long opos[D], const long odim[D], void* optr, const long ipos[D], const long idim[D], const void* iptr, size_t size)
+{
+	long istr[D];
+	long ostr[D];
+	md_calc_strides(D, istr, idim, size);
+	md_calc_strides(D, ostr, odim, size);
+
+	md_move_block2(D, dim, opos, odim, ostr, optr, ipos, idim, istr, iptr, size);
+}
+
+
 /**
  * Copy a block from an array to another array (with strides)
  *
@@ -688,31 +736,20 @@ void md_copy_block2(unsigned int D, const long pos[D], const long odim[D], const
 
 	for (unsigned int i = 0; i < D; i++) {
 
+		assert((idim[i] != odim[i]) || (0 == pos[i]));
+
 		dim[i] = MIN(odim[i], idim[i]);
 		ipos[i] = 0;
 		opos[i] = 0;
 
-		if (idim[i] != dim[i]) {
-
-			assert((0 <= pos[i]) && (pos[i] <= idim[i] - dim[i]));
+		if (idim[i] != dim[i])
 			ipos[i] = pos[i];
 
-		} else
-		if (odim[i] != dim[i]) {
-
-			assert((0 <= pos[i]) && (pos[i] <= odim[i] - dim[i]));
+		if (odim[i] != dim[i])
 			opos[i] = pos[i];
-
-		} else {
-
-			assert(0 == pos[i]);
-		}
 	}
-		
-	long ioff = md_calc_offset(D, istr, ipos);
-	long ooff = md_calc_offset(D, ostr, opos);
 
-	md_copy2(D, dim, ostr, optr + ooff, istr, iptr + ioff, size);
+	md_move_block2(D, dim, opos, odim, ostr, optr, ipos, idim, istr, iptr, size);
 }
 
 
diff --git a/src/num/multind.h b/src/num/multind.h
index ed6d6b18b..5cde4123d 100644
--- a/src/num/multind.h
+++ b/src/num/multind.h
@@ -50,6 +50,9 @@ extern void md_copy2(unsigned int D, const long dim[__VLA(D)], const long ostr[_
 extern void md_copy(unsigned int D, const long dim[__VLA(D)],  void* optr, const void* iptr, size_t size);
 extern void md_copy_block2(unsigned int D, const long pos[__VLA(D)], const long odim[__VLA(D)], const long ostr[__VLA(D)], void* optr, const long idim[__VLA(D)], const long istr[__VLA(D)], const void* iptr, size_t size);
 extern void md_copy_block(unsigned int D, const long pos[__VLA(D)], const long odim[__VLA(D)], void* optr, const long idim[__VLA(D)], const void* iptr, size_t size);
+extern void md_move_block2(unsigned int D, const long dim[__VLA(D)], const long opos[__VLA(D)], const long odim[__VLA(D)], const long ostr[__VLA(D)], void* optr, const long ipos[__VLA(D)], const long idim[__VLA(D)], const long istr[__VLA(D)], const void* iptr, size_t size);
+extern void md_move_block(unsigned int D, const long dim[__VLA(D)], const long opos[__VLA(D)], const long odim[__VLA(D)], void* optr, const long ipos[__VLA(D)], const long idim[__VLA(D)], const void* iptr, size_t size);
+
 extern void md_resize(unsigned int D, const long odim[__VLA(D)], void* optr, const long idim[__VLA(D)], const void* iptr, size_t size);
 extern void md_resizec(unsigned int D, const long odim[__VLA(D)], void* optr, const long idim[__VLA(D)], const void* iptr, size_t size);
 extern void md_fill2(unsigned int D, const long dim[__VLA(D)], const long str[__VLA(D)], void* ptr, const void* iptr, size_t size);
@@ -102,6 +105,7 @@ extern void md_singleton_strides(unsigned int D, long strs[__VLA(D)]);
 extern void md_set_dims(unsigned int D, long dims[__VLA(D)], long val);
 extern void md_min_dims(unsigned int D, unsigned long flags, long odims[__VLA(D)], const long idims1[__VLA(D)], const long idims2[__VLA(D)]);
 extern _Bool md_is_index(unsigned int D, const long pos[__VLA(D)], const long dims[__VLA(D)]);
+extern _Bool md_check_dimensions(unsigned int N, const long dims[__VLA(N)], unsigned int flags);
 extern void md_permute_dims(unsigned int D, const unsigned int order[__VLA(D)], long odims[__VLA(D)], const long idims[__VLA(D)]);
 extern void md_transpose_dims(unsigned int D, unsigned int dim1, unsigned int dim2, long odims[__VLA(D)], const long idims[__VLA(D)]);
 
@@ -110,6 +114,16 @@ extern void md_transpose_dims(unsigned int D, unsigned int dim1, unsigned int di
 #define MD_DIMS(...) MD_MAKE_ARRAY(long, __VA_ARGS__)
 
 
+#define MD_CAST_ARRAY2(T, N, dims, x, a, b) \
+	(assert(((a) < (b)) && !md_check_dimensions((N), (dims), (1 << (a)) | (1 << (b)))), \
+					*(T (*)[(dims)[b]][(dims)[a]])(x))
+#define MD_CAST_ARRAY3(T, N, dims, x, a, b, c) \
+	(assert(((a) < (b)) && ((b) < (c)) && !md_check_dimensions((N), (dims), (1 << (a)) | (1 << (b) | (1 << (c)))), \
+					*(T (*)[(dims)[c]][(dims)[b]][(dims)[a]])(x))
+
+
+
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/nusense.c b/src/nusense.c
index 7e30407ed..853e9810a 100644
--- a/src/nusense.c
+++ b/src/nusense.c
@@ -251,10 +251,10 @@ int main_nusense(int argc, char* argv[])
 	if (l1wav) {
 
 		long minsize[DIMS] = { [0 ... DIMS - 1] = 1 };
-		minsize[0] = MIN(img_dims[0], 16);
-		minsize[1] = MIN(img_dims[1], 16);
-		minsize[2] = MIN(img_dims[2], 16);
-		thresh_op = prox_wavethresh_create(DIMS, img_dims, FFT_FLAGS, minsize, lambda, randshift, use_gpu);
+		minsize[0] = MIN(img_dims[0], 8);
+		minsize[1] = MIN(img_dims[1], 8);
+		minsize[2] = MIN(img_dims[2], 8);
+		thresh_op = prox_wavethresh_create(DIMS, img_dims, 3, minsize, lambda, randshift, use_gpu);
 	}
 
 	// apply scaling
diff --git a/src/resize.c b/src/resize.c
index 6db304189..03b0c235d 100644
--- a/src/resize.c
+++ b/src/resize.c
@@ -4,6 +4,7 @@
  * 
  * Author; 
  * 2012-2013 Martin Uecker <uecker@eecs.berkeley.edu>
+ * 2014 Jonathan Tamir <jtamir@eecs.berkeley.edu>
  */
 
 #define _GNU_SOURCE
@@ -26,13 +27,13 @@
 
 static void usage(const char* name, FILE* fp)
 {
-	fprintf(fp, "Usage: %s [-c] dimension size <input> <output>\n", name);
+	fprintf(fp, "Usage: %s [-c] dim1 size1 ... dimn sizen <input> <output>\n", name);
 }
 
 static void help(void)
 {
 	printf(	"\n"
-		"Resizes an array along dimension to size by truncating or zero-padding.\n"
+		"Resizes an array along dimensions to sizes by truncating or zero-padding.\n"
 		"\n"
 		"-c\tcenter\n");
 }
@@ -62,7 +63,7 @@ int main_resize(int argc, char* argv[])
 		}
 	}
 
-	if (argc - optind != 4) {
+	if (argc - optind < 4) {
 
 		usage(argv[0], stderr);
 		exit(1);
@@ -70,23 +71,27 @@ int main_resize(int argc, char* argv[])
 
 	unsigned int N = DIMS;
 
-	unsigned int dim = atoi(argv[optind + 0]);
-	unsigned int count = atoi(argv[optind + 1]);
-
-	assert(dim < N);
-	assert(count >= 1);
+	int count = argc - optind - 2;
+	assert((count > 0) && (count % 2 == 0));
 
 	long in_dims[N];
 	long out_dims[N];
 
-	void* in_data = load_cfl(argv[optind + 2], N, in_dims);
+	void* in_data = load_cfl(argv[argc - 2], N, in_dims);
+	md_copy_dims(N, out_dims, in_dims);
 	
-	for (unsigned int i = 0; i < N; i++)
-		out_dims[i] = in_dims[i];
+	for (int i = 0; i < count; i += 2) {
+
+		unsigned int dim = atoi(argv[optind + i]);
+		unsigned int size = atoi(argv[optind + i + 1]);
 
-	out_dims[dim] = count;
+		assert(dim < N);
+		assert(size >= 1);
+
+		out_dims[dim] = size;
+	}
 
-	void* out_data = create_cfl(argv[optind + 3], N, out_dims);
+	void* out_data = create_cfl(argv[argc - 1], N, out_dims);
 
 	(center ? md_resizec : md_resize)(N, out_dims, out_data, in_dims, in_data, sizeof(complex float));
 
diff --git a/src/sake.c b/src/sake.c
index 59117a9de..c005cbb4d 100644
--- a/src/sake.c
+++ b/src/sake.c
@@ -34,7 +34,7 @@
 #include "sake/sake.h"
 
 #undef DIMS // FIXME
-#define DIMS 5
+#define DIMS 16
 
 // 8 channels: alpha 0.2, 50 iter
 
diff --git a/src/sake/sake.c b/src/sake/sake.c
index 2675a8cd7..c31d1cb29 100644
--- a/src/sake/sake.c
+++ b/src/sake/sake.c
@@ -4,6 +4,7 @@
  *
  * Authors:
  * 2013 Martin Uecker <uecker@eecs.berkeley.edu>
+ * 2014 Jonathan Tamir <jtamir@eecs.berkeley.edu>
  *
  *
  * Peter J. Shin, Peder E.Z. Larson, Michael A. Ohliger, Michael Elad,
@@ -36,7 +37,7 @@
 #include "sake.h"
 
 #undef DIMS // FIXME
-#define DIMS 5
+#define DIMS 16
 
 #if 0
 static float thresh(float lambda, float x)
@@ -79,15 +80,16 @@ static void ravine(unsigned int N, const long dims[N], float* ftp, complex float
         ft = (1.f + sqrtf(1.f + 4.f * ft * ft)) / 2.f;
         *ftp = ft;
 
-	md_swap(N, dims, xa, xb, sizeof(complex float));
+	md_swap(N, dims, xa, xb, CFL_SIZE);
 	complex float val = (1.f - tfo) / ft - 1.f;
+
 	long dims1[N];
-	for (unsigned int i = 0; i < N; i++)
-		dims1[i] = 1.;
+	md_singleton_dims(N, dims1);
+
 	long strs1[N];
 	long strs[N];
-	md_calc_strides(N, strs1, dims1, sizeof(complex float));
-	md_calc_strides(N, strs, dims, sizeof(complex float));
+	md_calc_strides(N, strs1, dims1, CFL_SIZE);
+	md_calc_strides(N, strs, dims, CFL_SIZE);
 
 	md_zfmac2(N, dims, strs, xa, strs1, &val, strs, xa);
 	val *= -1.;
@@ -99,8 +101,9 @@ static void ravine(unsigned int N, const long dims[N], float* ftp, complex float
 
 
 
-static void lowrank(float alpha, const long dims[5], complex float* matrix)
+static void lowrank(float alpha, int D, const long dims[D], complex float* matrix)
 {
+#if 0
 	long x = dims[0];
 	long y = dims[1];
 	long z = dims[2];
@@ -113,19 +116,48 @@ static void lowrank(float alpha, const long dims[5], complex float* matrix)
 	int kz = MIN(6, z);
 	
 	long calreg_dims[4] = { x, y, z, channels };
-	long kernel_dims[4] = { kx, ky, kz, channels };
+	long kern_dims[4] = { kx, ky, kz, channels };
 
 	debug_printf(DP_INFO, "%ld %ld %ld %ld\n", x, y, z, channels);
 
-	long calmat_dims[2] = { (x - kx + 1) * (y - ky + 1) * (z - kz + 1), md_calc_size(4, kernel_dims) };
+	long calmat_dims[2] = { (x - kx + 1) * (y - ky + 1) * (z - kz + 1), md_calc_size(4, kern_dims) };
 
-	complex float* calmat = md_alloc(2, calmat_dims, sizeof(complex float));
+	complex float* calmat = md_alloc(2, calmat_dims, CFL_SIZE);
 //	complex float* calmat = create_cfl("calmat", 2, calmat_dims);
 
 	long str[4];
-	md_calc_strides(4, str, calreg_dims, sizeof(complex float));
+	md_calc_strides(4, str, calreg_dims, CFL_SIZE);
+
+	casorati_matrix(4, kern_dims, calmat_dims, calmat, calreg_dims, str, matrix);
+#else
+	assert(1 == dims[MAPS_DIM]);
+
+	debug_printf(DP_INFO, "mat_dims = \t");
+	debug_print_dims(DP_INFO, D, dims);
+
+	long kern_min[4] = { 6, 6, 6, dims[COIL_DIM] };
+	long kern_dims[D];
+
+	md_set_dims(D, kern_dims, 1);
+	md_min_dims(4, ~0u, kern_dims, kern_min, dims);
+
+	debug_printf(DP_INFO, "kern_dims = \t");
+	debug_print_dims(DP_INFO, D, kern_dims);
 
-	casorati_matrix(4, kernel_dims, calmat_dims, calmat, calreg_dims, str, matrix);
+	long calmat_dims[2];
+	casorati_dims(D, calmat_dims, kern_dims, dims);
+
+	debug_printf(DP_INFO, "calmat_dims = \t");
+	debug_print_dims(DP_INFO, 2, calmat_dims);
+
+	complex float* calmat = md_alloc(2, calmat_dims, CFL_SIZE);
+
+	long str[D];
+	md_calc_strides(D, str, dims, CFL_SIZE);
+
+	casorati_matrix(D, kern_dims, calmat_dims, calmat, dims, str, matrix);
+
+#endif
 
 	int N = calmat_dims[0];
 	int M = calmat_dims[1];
@@ -138,8 +170,8 @@ static void lowrank(float alpha, const long dims[5], complex float* matrix)
 		long dimsU[2] = { N, N };
 		long dimsV[2] = { M, M };
 
-		complex float* U = md_alloc(2, dimsU, sizeof(complex float));
-		complex float* VT = md_alloc(2, dimsV, sizeof(complex float));
+		complex float* U = md_alloc(2, dimsU, CFL_SIZE);
+		complex float* VT = md_alloc(2, dimsV, CFL_SIZE);
 		//	complex float* U = create_cfl("U", 2, dimsU);
 		//	complex float* VT = create_cfl("VT", 2, dimsV);
 		float* S = xmalloc(MIN(N, M) * sizeof(float));
@@ -154,10 +186,10 @@ static void lowrank(float alpha, const long dims[5], complex float* matrix)
 		// put it back together
 		long dimU2[2] = { N, MIN(N, M) };
 		long dimV2[2] = { MIN(N, M), M };
-		complex float* U2 = md_alloc(2, dimU2, sizeof(complex float));
-		complex float* V2 = md_alloc(2, dimV2, sizeof(complex float));
-		md_resize(2, dimU2, U2, dimsU, U, sizeof(complex float));
-		md_resize(2, dimV2, V2, dimsV, VT, sizeof(complex float));
+		complex float* U2 = md_alloc(2, dimU2, CFL_SIZE);
+		complex float* V2 = md_alloc(2, dimV2, CFL_SIZE);
+		md_resize(2, dimU2, U2, dimsU, U, CFL_SIZE);
+		md_resize(2, dimV2, V2, dimsV, VT, CFL_SIZE);
 
 		for (int i = 0; i < M; i++) {
 
@@ -178,11 +210,17 @@ static void lowrank(float alpha, const long dims[5], complex float* matrix)
 		free(S);
 	}
 
+#if 0
 	//md_clear(5, dims, matrix, sizeof(complex float));
-	casorati_matrixH(4, kernel_dims, calreg_dims, str, matrix, calmat_dims, calmat);
+	casorati_matrixH(4, kern_dims, calreg_dims, str, matrix, calmat_dims, calmat);
 	md_zsmul(5, dims, matrix, matrix, 1. / (double)(kx * ky * kz)); // FIXME: not right at the border
 
 	//unmap_cfl(2, calmat_dims, calmat);
+#else
+	md_clear(D, dims, matrix, CFL_SIZE);
+	casorati_matrixH(D, kern_dims, dims, str, matrix, calmat_dims, calmat);
+	md_zsmul(D, dims, matrix, matrix, 1. / (double)md_calc_size(3, kern_dims)); // FIXME: not right at the border
+#endif
 	md_free(calmat);
 }
 
@@ -191,32 +229,31 @@ static void lowrank(float alpha, const long dims[5], complex float* matrix)
 void lrmc(float alpha, int iter, float lambda, int N, const long dims[N], complex float* out, const complex float* in)
 {
 	long dims1[N];
-	memcpy(dims1, dims, N * sizeof(long));
-	dims1[3] = 1;
+	md_select_dims(N, ~COIL_FLAG, dims1, dims);
 
-	md_copy(N, dims, out, in, sizeof(complex float));
+	md_copy(N, dims, out, in, CFL_SIZE);
 
-	complex float* pattern = md_alloc(N, dims1, sizeof(complex float));
+	complex float* pattern = md_alloc(N, dims1, CFL_SIZE);
 
-	assert(5 == N);
+	//assert(5 == N);
 	estimate_pattern(N, dims, COIL_DIM, pattern, in);
 
-	complex float* comp = md_alloc(N, dims1, sizeof(complex float));
+	complex float* comp = md_alloc(N, dims1, CFL_SIZE);
 	md_zfill(N, dims1, comp, 1.);
 
-	lowrank(-1., dims1, comp);
+	lowrank(-1., N, dims1, comp);
 
 #ifdef RAVINE
-	complex float* o = md_alloc(N, dims, sizeof(complex float));
-	md_clear(N, dims, o, sizeof(complex float));
+	complex float* o = md_alloc(N, dims, CFL_SIZE);
+	md_clear(N, dims, o, CFL_SIZE);
 	float fl = 1.;
 #endif
 
 	long strs1[N];
-	md_calc_strides(N, strs1, dims1, sizeof(complex float));
+	md_calc_strides(N, strs1, dims1, CFL_SIZE);
 
 	long strs[N];
-	md_calc_strides(N, strs, dims, sizeof(complex float));
+	md_calc_strides(N, strs, dims, CFL_SIZE);
 	
 
 	for (int i = 0; i < iter; i++) {
@@ -228,7 +265,7 @@ void lrmc(float alpha, int iter, float lambda, int N, const long dims[N], comple
 		else
 			data_consistency(dims, out, pattern, in, out);
 
-		lowrank(alpha, dims, out);
+		lowrank(alpha, N, dims, out);
 		md_zdiv2(N, dims, strs, out, strs, out, strs1, comp);
 #ifdef RAVINE
 		ravine(N, dims, &fl, out, o);
diff --git a/src/sense.c b/src/sense.c
index 5da82987a..aebc87e47 100644
--- a/src/sense.c
+++ b/src/sense.c
@@ -81,7 +81,7 @@ int main_sense(int argc, char* argv[])
 	float restrict_fov = -1.;
 	const char* pat_file = NULL;
 	const char* traj_file = NULL;
-	char image_truth_fname[100];
+	const char* image_truth_file = NULL;
 	bool im_truth = false;
 	bool scale_im = false;
 
@@ -90,7 +90,7 @@ int main_sense(int argc, char* argv[])
 	float admm_rho = iter_admm_defaults.rho;
 
 	int c;
-	while (-1 != (c = getopt(argc, argv, "Fq:l:r:s:i:u:o:O:f:t:cTImghp:Sd:H"))) {
+	while (-1 != (c = getopt(argc, argv, "Fq:l:r:s:i:u:o:O:f:t:cT:Imghp:Sd:H"))) {
 		switch(c) {
 
 		case 'H':
@@ -107,7 +107,8 @@ int main_sense(int argc, char* argv[])
 
 		case 'T':
 			im_truth = true;
-			sprintf(image_truth_fname, "%s", optarg);
+			image_truth_file = strdup(optarg);
+			assert(NULL != image_truth_file);
 			break;
 
 		case 'd':
@@ -216,8 +217,8 @@ int main_sense(int argc, char* argv[])
 	complex float* kspace = load_cfl(argv[optind + 0], DIMS, ksp_dims);
 	complex float* maps = load_cfl(argv[optind + 1], DIMS, map_dims);
 
-	for (unsigned int i = 0; i < DIMS; i++)
-		max_dims[i] = MAX(ksp_dims[i], map_dims[i]);
+	md_copy_dims(DIMS, max_dims, ksp_dims);
+	max_dims[MAPS_DIM] = map_dims[MAPS_DIM];
 
 	md_select_dims(DIMS, ~COIL_FLAG, pat_dims, ksp_dims);
 	md_select_dims(DIMS, ~COIL_FLAG, img_dims, max_dims);
@@ -339,7 +340,7 @@ int main_sense(int argc, char* argv[])
 
 	if (im_truth) {
 
-		image_truth = load_cfl(image_truth_fname, DIMS, img_truth_dims);
+		image_truth = load_cfl(image_truth_file, DIMS, img_truth_dims);
 		//md_zsmul(DIMS, img_dims, image_truth, image_truth, 1. / scaling);
 	}
 
@@ -417,6 +418,13 @@ int main_sense(int argc, char* argv[])
 	unmap_cfl(DIMS, ksp_dims, kspace);
 	unmap_cfl(DIMS, img_dims, image);
 
+	if (im_truth) {
+
+		free((void*)image_truth_file);
+		unmap_cfl(DIMS, img_dims, image_truth);
+	}
+
+
 	double end_time = timestamp();
 	debug_printf(DP_INFO, "Total Time: %f\n", end_time - start_time);
 	exit(0);
diff --git a/src/simu/shepplogan.c b/src/simu/shepplogan.c
index 5a3baaa2a..032edbdfe 100644
--- a/src/simu/shepplogan.c
+++ b/src/simu/shepplogan.c
@@ -8,6 +8,7 @@
 #include <complex.h>
 #include <math.h>
 #include <assert.h>
+#include <stdbool.h>
 
 #ifdef USE_GSL
 #include <gsl/gsl_specfunc.h>
@@ -149,7 +150,7 @@ complex double krectangle(const double center[2], const double axis[2], double a
 }
 
 
-complex double phantom(unsigned int N, const struct ellipsis_s arr[N], const double pos[2], _Bool ksp)
+complex double phantom(unsigned int N, const struct ellipsis_s arr[N], const double pos[2], bool ksp)
 {
 	complex double res = 0.;
 
@@ -159,7 +160,7 @@ complex double phantom(unsigned int N, const struct ellipsis_s arr[N], const dou
 	return res;
 }
 
-complex double phantomX(unsigned int N, const struct ellipsis_s arr[N], const double pos[2], _Bool ksp)
+complex double phantomX(unsigned int N, const struct ellipsis_s arr[N], const double pos[2], bool ksp)
 {
 	complex double res = 0.;
 
diff --git a/src/spow.c b/src/spow.c
new file mode 100644
index 000000000..906cb6fdd
--- /dev/null
+++ b/src/spow.c
@@ -0,0 +1,53 @@
+/* Copyright 2014. The Regents of the University of California.
+ * 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 Martin Uecker <uecker@eecs.berkeley.edu>
+ */
+
+#include <stdlib.h>
+#include <assert.h>
+#include <stdbool.h>
+#include <complex.h>
+#include <stdio.h>
+
+#include "num/multind.h"
+#include "num/flpmath.h"
+
+#include "misc/mmio.h"
+#include "misc/misc.h"
+
+#ifndef DIMS
+#define DIMS 16
+#endif
+
+static const char* usage_str = "exponent <input> <output>";
+static const char* help_str = "Raise array to the power of {exponent}. The exponent can be a complex number.\n";
+
+
+int main_spow(int argc, char* argv[argc])
+{
+	mini_cmdline(argc, argv, 3, usage_str, help_str);
+
+	complex float expo;
+
+	if (0 != parse_cfl(&expo, argv[1])) {
+
+		fprintf(stderr, "ERROR: exponent %s is not a number.\n", argv[1]);
+		exit(1);
+	}
+
+	const int N = DIMS;
+	long dims[N];
+	complex float* idata = load_cfl(argv[2], N, dims);
+	complex float* odata = create_cfl(argv[3], N, dims);
+		
+	md_zspow(N, dims, odata, idata, expo);
+
+	unmap_cfl(N, dims, idata);
+	unmap_cfl(N, dims, odata);
+	exit(0);
+}
+
+
diff --git a/src/twixread.c b/src/twixread.c
index a54ebb821..f7d4ce4f8 100644
--- a/src/twixread.c
+++ b/src/twixread.c
@@ -104,7 +104,7 @@ struct mdh2 {	// second part of mdh
 
 
 
-static void siemens_adc_read(bool vd, int fd, const long dims[DIMS], long pos[DIMS], complex float* buf)
+static int siemens_adc_read(bool vd, int fd, bool linectr, bool partctr, const long dims[DIMS], long pos[DIMS], complex float* buf)
 {
 	char scan_hdr[vd ? 192 : 0];
 	xread(fd, scan_hdr, sizeof(scan_hdr));
@@ -120,9 +120,9 @@ static void siemens_adc_read(bool vd, int fd, const long dims[DIMS], long pos[DI
 		if (0 == pos[COIL_DIM]) {
 
 			// TODO: rethink this
-			pos[PHS1_DIM] = mdh.sLC[0]; // - mdh.linectr;
+			pos[PHS1_DIM] = mdh.sLC[0] + (linectr ? mdh.linectr : 0);
 			pos[SLICE_DIM] = mdh.sLC[2];
-			pos[PHS2_DIM] = mdh.sLC[3]; // - mdh.partctr;
+			pos[PHS2_DIM] = mdh.sLC[3] + (partctr ? mdh.partctr : 0);
 			pos[TE_DIM] = mdh.sLC[4];
 			pos[TIME_DIM] = mdh.sLC[6];
 			pos[TIME2_DIM] = mdh.sLC[7];
@@ -130,15 +130,17 @@ static void siemens_adc_read(bool vd, int fd, const long dims[DIMS], long pos[DI
 
 		debug_print_dims(DP_DEBUG1, DIMS, pos);
 
-		if (dims[READ_DIM] != mdh.samples)
-			error("wrong number of samples");
+		if (dims[READ_DIM] != mdh.samples) {
 
-		assert(md_is_index(DIMS, pos, dims));
+			debug_printf(DP_WARN, "Wrong number of samples.\n");
+			return -1;
+		}
 
 		xread(fd, buf + pos[COIL_DIM] * dims[READ_DIM], dims[READ_DIM] * CFL_SIZE);
 	}
 
 	pos[COIL_DIM] = 0;
+	return 0;
 }
 
 
@@ -161,6 +163,8 @@ static void help(void)
 		"-s S\tnumber of slices\n"
 		"-c C\tnumber of channels\n"
 		"-a A\ttotal number of ADCs\n"
+		"-L\tuse linectr offset\n"
+		"-P\tuse partctr offset\n"
 		"-h\thelp\n");
 }
 
@@ -171,10 +175,13 @@ int main(int argc, char* argv[argc])
 	int c;
 	long adcs = 0;
 
+	bool linectr = false;
+	bool partctr = false;
+
 	long dims[DIMS];
 	md_singleton_dims(DIMS, dims);
 
-	while (-1 != (c = getopt(argc, argv, "x:y:z:s:c:a:h"))) {
+	while (-1 != (c = getopt(argc, argv, "x:y:z:s:c:a:PLh"))) {
 		switch (c) {
 
 		case 'x':
@@ -201,6 +208,14 @@ int main(int argc, char* argv[argc])
 			dims[COIL_DIM] = atoi(optarg);
 			break;
 
+		case 'P':
+			partctr = true;
+			break;
+
+		case 'L':
+			linectr = true;
+			break;
+
 		case 'h':
 			usage(argv[0], stdout);
 			help();
@@ -231,6 +246,7 @@ int main(int argc, char* argv[argc])
 	bool vd = siemens_meas_setup(ifd, &hdr);
 
 	complex float* out = create_cfl(argv[optind + 1], DIMS, dims);
+	md_clear(DIMS, dims, out, CFL_SIZE);
 
 	long adc_dims[DIMS];
 	md_select_dims(DIMS, READ_FLAG|COIL_FLAG, adc_dims, dims);
@@ -240,13 +256,24 @@ int main(int argc, char* argv[argc])
 	while (adcs--) {
 
 		long pos[DIMS] = { [0 ... DIMS - 1] = 0 };
-		siemens_adc_read(vd, ifd, dims, pos, buf);
+
+		if (-1 == siemens_adc_read(vd, ifd, linectr, partctr, dims, pos, buf)) {
+
+			debug_printf(DP_WARN, "Stopping.\n");
+			break;
+		}
 
 		debug_print_dims(DP_DEBUG1, DIMS, pos);
+
+		if (!md_is_index(DIMS, pos, dims)) {
+
+			debug_printf(DP_WARN, "Index out of bounds.\n");
+			continue;
+		}
+
 		md_copy_block(DIMS, pos, dims, out, adc_dims, buf, CFL_SIZE); 
 	}
 
-
 	md_free(buf);
 	unmap_cfl(DIMS, dims, out);
 	exit(0);
diff --git a/src/walsh.c b/src/walsh.c
index bfa9ee339..752164a6d 100644
--- a/src/walsh.c
+++ b/src/walsh.c
@@ -88,30 +88,30 @@ int main_walsh(int argc, char* argv[])
 	}
 
 
-	long dims[KSPACE_DIMS];
+	long dims[DIMS];
 
-	complex float* in_data = load_cfl(argv[optind + 0], KSPACE_DIMS, dims);
+	complex float* in_data = load_cfl(argv[optind + 0], DIMS, dims);
 
 	assert((dims[0] == 1) || (calsize[0] < dims[0]));
 	assert((dims[1] == 1) || (calsize[1] < dims[1]));
 	assert((dims[2] == 1) || (calsize[2] < dims[2]));
-	assert(1 == dims[4]);
+	assert(1 == dims[MAPS_DIM]);
 
-	long caldims[KSPACE_DIMS];
+	long caldims[DIMS];
 	complex float* cal_data = extract_calib(caldims, calsize, dims, in_data, false);
-	unmap_cfl(KSPACE_DIMS, dims, in_data);
+	unmap_cfl(DIMS, dims, in_data);
 
 	printf("Calibration region %ldx%ldx%ld\n", caldims[0], caldims[1], caldims[2]);
 
-	dims[3] = dims[3] * (dims[3] + 1) / 2;
-	complex float* out_data = create_cfl(argv[optind + 1], KSPACE_DIMS, dims);
+	dims[COIL_DIM] = dims[COIL_DIM] * (dims[COIL_DIM] + 1) / 2;
+	complex float* out_data = create_cfl(argv[optind + 1], DIMS, dims);
 
 	walsh(bsize, dims, out_data, caldims, cal_data);
 
 	printf("Done.\n");
 
 	md_free(cal_data);
-	unmap_cfl(KSPACE_DIMS, dims, out_data);
+	unmap_cfl(DIMS, dims, out_data);
 	exit(0);
 }
 
diff --git a/src/wavelet2/wavelet.c b/src/wavelet2/wavelet.c
index b69f5e826..6503b0a50 100644
--- a/src/wavelet2/wavelet.c
+++ b/src/wavelet2/wavelet.c
@@ -683,8 +683,8 @@ void softthresh_cpu(struct wavelet_plan_s* plan, data_t* coeff, scalar_t thresh)
 	int numMax = plan->numCoeff_tr;
 	int i;
 #pragma omp parallel for
-	for(i = plan->numCoarse_tr; i < numMax; i++)
-		//for(i = 0; i < numMax; i++)
+//	for(i = plan->numCoarse_tr; i < numMax; i++)
+	for(i = 0; i < numMax; i++)
 	{
 		scalar_t norm = cabsf(coeff[i]);
 		scalar_t red = norm - thresh;
@@ -944,7 +944,7 @@ const float wavelet2_cdf44[4][10] = {
 
 void circshift(struct wavelet_plan_s* plan, data_t* data) 
 {
-	_Bool shift = false;
+	bool shift = false;
 	int N = plan->numdims_tr;
 
 	for (int i = 0; i < N; i++)
@@ -962,7 +962,7 @@ void circshift(struct wavelet_plan_s* plan, data_t* data)
 
 void circunshift(struct wavelet_plan_s* plan, data_t* data) 
 {
-	_Bool shift = false;
+	bool shift = false;
 	int N = plan->numdims_tr;
 	long rand_shifts[N];