ACKNOWLEDGEMENTS

 
-Contributors:
+We want to acknowledge the following persons for supporting
+BART by contributing source code, testing, feedback, data,
+bug reports, etc.
+
+
+(alphabetrical)
 
 Marcus T. Alley
+Michael Anderson
+Jakob Asslaender
 Dara Bahri
+Soumick Chatterjee
 Joseph Y. Cheng
+Sofia Dimoudi
 Siddharth Iyer
-Miki Lustig
-Mark Murphy
+Hans Johnson
+Christian Holme
+Gregory R. Lee
+Evan G. Levine
+Tim Loderhose
+Michael Lustig
+Lyu Mengye
+Damien Nguyen
 Frank Ong
-Jonathan Tamir
+Sebastian Rosenzweig
+Nick Scholand
+David Smith
+Jonathan I. Tamir
+Michelle Tamir (logo)
+Johannes Toger
 Martin Uecker
 Shreyas S. Vasanawala
 Sana Vaziri
-Pat Virtue
+Patrick Virtue
+Simon Yeung
 Tao Zhang
 
-Logo designed by: Michelle Tamir
-
-
-Different parts of this software have been written for
-projects funded by:
-
-American Heart Association Grant 12BGIA9660006
-NIH Grant P41RR09784 and Grant R01EB009690
-UC Discovery Grant 193037
-Sloan Research Fellowship
-GE Healthcare
-
-A personal donation by David Donoho
-
-
-

CMakeLists.txt

@@ -542,10 +542,12 @@ endif()
 
 # ------------------------------------------------------------------------------
 
-# Include integration testing
-if(NOT BART_MEMONLY_CFL)
-  add_subdirectory(tests)
-endif()
+##Testing broke in 2c8c1a27 when adding a test.
+
+## Include integration testing
+#if(NOT BART_MEMONLY_CFL)
+#  add_subdirectory(tests)
+#endif()
 
 # ==============================================================================
 

LICENSE

-Copyright (c) 2013-2017. The Regents of the University of California.
-Copyright (c) 2013-2017. BART Developer Team and Contributors.
+Copyright (c) 2013-2018. The Regents of the University of California.
+Copyright (c) 2013-2018. BART Developer Team and Contributors.
 Copyright (c) 2012. Intel Corporation. (src/lapacke/)
 All rights reserved.
 

Makefile

@@ -407,6 +407,7 @@ endif
 ifeq ($(SLINK),1)
 # work around fortran problems with static linking
 LDFLAGS += -static -Wl,--whole-archive -lpthread -Wl,--no-whole-archive -Wl,--allow-multiple-definition
+LIBS += -lmvec
 BLAS_L += -llapack -lblas -lgfortran -lquadmath
 endif
 
@@ -490,14 +491,14 @@ endif
 
 .SECONDEXPANSION:
 $(TARGETS): % : src/main.c $(srcdir)/%.o $$(MODULES_%) $(MODULES)
-	$(LINKER) $(LDFLAGS) $(CFLAGS) $(CPPFLAGS) -Dmain_real=main_$@ -o $@ $+ $(FFTW_L) $(CUDA_L) $(BLAS_L) $(PNG_L) $(ISMRM_L) -lm
+	$(LINKER) $(LDFLAGS) $(CFLAGS) $(CPPFLAGS) -Dmain_real=main_$@ -o $@ $+ $(FFTW_L) $(CUDA_L) $(BLAS_L) $(PNG_L) $(ISMRM_L) $(LIBS) -lm
 #	rm $(srcdir)/$@.o
 
 UTESTS=$(shell $(root)/utests/utests-collect.sh ./utests/$@.c)
 
 .SECONDEXPANSION:
 $(UTARGETS): % : utests/utest.c utests/%.o $$(MODULES_%) $(MODULES)
-	$(CC) $(LDFLAGS) $(CFLAGS) $(CPPFLAGS) -DUTESTS="$(UTESTS)" -o $@ $+ $(FFTW_L) $(CUDA_L) $(BLAS_L) -lm
+	$(CC) $(LDFLAGS) $(CFLAGS) $(CPPFLAGS) -DUTESTS="$(UTESTS)" -o $@ $+ $(FFTW_L) $(CUDA_L) $(BLAS_L) $(LIBS) -lm
 
 
 # linker script version - does not work on MacOS X

README

@@ -103,8 +103,14 @@ of gcc from MacPorts (http://www.macports.org/).
 
 ### 2.1.3. Windows
 
+Note that Windows is currently not supported. Nevertheless, many people
+use BART on Windows. If you have trouble to get it to work, you may
+want to try an older version (e.g. version 4.00).
+
+
 The recommended way to use BART on Windows is with the Windows Subsystem for
-Linux (WSL) which is available for Windows 10.
+Linux (WSL) which is available for Windows 10. If you use WSL, you need to
+use the clang compiler.
 
 
 BART should also work with Cygwin:

build_targets.mk

 # Main build targets
 #
-TBASE=show slice crop resize join transpose squeeze flatten zeros ones flip circshift extract repmat bitmask reshape version delta copy casorati vec poly
-TFLP=scale invert conj fmac saxpy sdot spow cpyphs creal carg normalize cdf97 pattern nrmse mip avg cabs zexpj
+TBASE=show slice crop resize join transpose squeeze flatten zeros ones flip circshift extract repmat bitmask reshape version delta copy casorati vec poly index
+TFLP=scale invert conj fmac saxpy sdot spow cpyphs creal carg normalize cdf97 pattern nrmse mip avg cabs zexp
 TNUM=fft fftmod fftshift noise bench threshold conv rss filter mandelbrot wavelet window var std
 TRECO=pics pocsense sqpics itsense nlinv nufft rof sake wave lrmatrix estdims estshift estdelay wavepsf wshfl
 TCALIB=ecalib ecaltwo caldir walsh cc ccapply calmat svd estvar whiten

doc/bitmasks.txt

+
+A bitmask is a binary number where the individual bits
+are used to indicate something.
+
+For example, a bitmask is often used to select a subset
+of dimensions, e.g. if the FFT should be applited to
+dimensions 3 and 5 the corresponding bits at position
+3 and 5 are set in a bitmask:
+
+876543210  (position of the bit)
+000101000b (bitmask as binary number)
+
+In decimal this binary number is 40, so the command
+
+$ bart fft 40 input output
+
+would transform the dimensions 3 and 5 (counting from
+zero).
+
+The 'bitmask' command computes the required bitmask
+from the set of dimensions:
+
+$ bart bitmask 3 5
+40
+
+On the command-line both commands can also be combined:
+
+$ bart fft 'bart bitmask 3 5' input output
+

doc/references.txt

@@ -19,6 +19,35 @@ ISMRM Workshop on Data Sampling and Image Reconstruction, Sedona 2016.
 
 
 
+
+- reproducible publications using BART -
+
+Uecker M, Lustig M.
+Estimating Absolute-Phase Maps Using ESPIRiT and Virtual Conjugate Coils.
+Magn Reson Med 2017; 77:1201-1207.
+https://github.com/mrirecon/vcc-espirit
+
+Rosenzweig S, Holme HCM, Wilke RN, Voit D, Frahm J, Uecker M.
+Simultaneous Multi-Slice Reconstruction Using Regularized Nonlinear Inversion:
+SMS-NLINV.
+Magnetic Resonance in Medicine 2018; 79:2057-2066.
+https://github.com/mrirecon/sms-nlinv
+
+Rosenzweig S, Holme HCM, Uecker M.
+Simple Auto-Calibrated Gradient Delay Estimation From Few Spokes Using Radial
+Intersections (RING).
+Magnetic Resonance in Medicine, DOI:10.1002/mrm.27506 (2018)
+https://github.com/mrirecon/ring
+
+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]
+https://github.com/mrirecon/enlive
+
+
+
+
 - sensitivity-encoded parallel imaging -
 (commands: itsense, pocsense, bpsense, pics)
 
@@ -154,7 +183,7 @@ 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,
+ESPIRiT-Based Coil Compression for Cartesian Sampling.
 Annual Meeting ISMRM, Salt Lake City 2013,
 In: Proc Intl Soc Mag Reson Med 21; 2657.
 
@@ -162,7 +191,7 @@ In: Proc Intl Soc Mag Reson Med 21; 2657.
 
 
 - compressed sensing MRI -
-(commands: pocsense, rsense, pics)
+(commands: pocsense, pics)
 
 
 Block KT, Uecker M, and Frahm J.
@@ -221,8 +250,26 @@ In: Proc Intl Soc Mag Reson Med 17; 379.
 
 
 
+- trajectory correction -
+(commands: estdelay)
+
+
+Block KT, Uecker M.
+Simple Method for Adaptive Gradient-Delay Compensation in Radial MRI.
+Annual Meeting ISMRM, Montreal 2011,
+In: Proc. Intl. Soc. Mag. Reson. Med 19: 2816.
+
+Rosenzweig S, Holme HCM, Uecker M.
+Simple Auto-Calibrated Gradient Delay Estimation From Few Spokes Using Radial
+Intersections (RING).
+Magnetic Resonance in Medicine, DOI:10.1002/mrm.27506 (2018)
+
+
+
+
+
 - acceleration with graphical processing units -
-(commands: pocsense, rsense, nufft, pics)
+(commands: pocsense, nufft, pics, nlinv)
 
 
 Uecker M, Zhang S, Frahm J.
@@ -258,8 +305,9 @@ IEEE Trans Med Imaging 2012; 31:626-636.
 
 
 
-- applications -
 
+- applications -
+(an incomplete list of papers using BART...)
 
 Hollingsworth KG, Higgins DM, McCallum M, Ward L, Coombs A, Straub V.
 Investigating the quantitative fidelity of prospectively undersampled
@@ -318,3 +366,24 @@ Comprehensive Multi-Dimensional MRI for the Simultaneous Assessment
 of Cardiopulmonary Anatomy and Physiology.
 Scientific Reports 2017; 7:5330.
 
+Bao S, Tamir JI, Young JL, Tariq U, Uecker M, Lai P, Chen W, Lustig M,
+Vasanawala SS.
+Fast comprehensive single-sequence four-dimensional pediatric knee MRI with
+T2 shuffling.
+Journal of Magnetic Resonance Imaging 2017; 45:1700-1711.
+
+Rosenzweig S, Holme HCM, Wilke RN, Voit D, Frahm J, Uecker M.
+Simultaneous Multi-Slice Reconstruction Using Regularized Nonlinear Inversion:
+SMS-NLINV.
+Magnetic Resonance in Medicine 2018; 79:2057-2066.
+
+Roeloffs V, Rosenzweig S, Holme HCM, Uecker M, Frahm J.
+Frequency-modulated SSFP with radial sampling and subspace reconstruction:
+A time-efficient alternative to phase-cycled bSSFP.
+Magnetic Resonance in Medicine, DOI:10.1002/mrm.27505 (2018)
+
+Rosenzweig S, Holme HCM, Uecker M.
+Simple Auto-Calibrated Gradient Delay Estimation From Few Spokes Using Radial
+Intersections (RING).
+Magnetic Resonance in Medicine, DOI:10.1002/mrm.27506 (2018)
+

matlab/WSLPathCorrection.m

+%Soumick Chatterjee <soumick.chatterjee@ovgu.de>
+
+function [outData] = WSLPathCorrection(inData)
+    outData=inData;
+    for i = 'a':'z' %Replace drive letters with /mnt/
+        outData=strrep(outData,[i,':'],['/mnt/',i]); %if drive letter is supplied in lowercase
+        outData=strrep(outData,[upper(i),':'],['/mnt/',i]); %if drive letter is supplied as uppercase
+    end
+    outData = strrep(outData, '\', '/'); %Change windows filesep to linux filesep
+end

matlab/bart.m

-function [varargout] = bart(cmd, varargin);
+function [varargout] = bart(cmd, varargin)
 % BART	Call BART command from Matlab.
 %   [A, B] = bart('command', X, Y) call command with inputs X Y and outputs A B
 %
 % 2014-2016 Martin Uecker <uecker@med.uni-goettingen.de>
+% Edited for WSL, by Soumick Chatterjee <soumick.chatterjee@ovgu.de>
 
 
 	if nargin==0 || all(cmd==0)
@@ -11,12 +12,20 @@ function [varargout] = bart(cmd, varargin);
 	end
 
 	bart_path = getenv('TOOLBOX_PATH');
+    isWSL = false;
 
 	if isempty(bart_path)
 		if exist('/usr/local/bin/bart', 'file')
 			bart_path = '/usr/local/bin';
 		elseif exist('/usr/bin/bart', 'file')
 			bart_path = '/usr/bin';
+        else    
+            % Try to execute bart inside wsl, if it works, then it returns
+            % status 0
+            [bartstatus, ~] = system('wsl bart version -V');
+            if bartstatus==0
+                bart_path = '/usr/bin';
+                isWSL = true;
             else
                 error('Environment variable TOOLBOX_PATH is not set.');
             end
@@ -35,7 +44,7 @@ function [varargout] = bart(cmd, varargin);
 
 	in = cell(1, nargin - 1);
 
-	for i=1:nargin - 1,
+	for i=1:nargin - 1
 		in{i} = strcat(name, 'in', num2str(i));
 		writecfl(in{i}, varargin{i});
 	end
@@ -44,24 +53,32 @@ function [varargout] = bart(cmd, varargin);
 
 	out = cell(1, nargout);
 
-	for i=1:nargout,
+	for i=1:nargout
 		out{i} = strcat(name, 'out', num2str(i));
 	end
 
 	out_str = sprintf(' %s', out{:});
 
 	if ispc
+        if isWSL
+            % For WSL and modify paths
+            cmdWSL = WSLPathCorrection(cmd);
+            in_strWSL = WSLPathCorrection(in_str);
+            out_strWSL =  WSLPathCorrection(out_str);		
+            ERR = system(['wsl bart ', cmdWSL, ' ', in_strWSL, ' ', out_strWSL]);
+        else
             % For cygwin use bash and modify paths
             ERR = system(['bash.exe --login -c ', ...
                 strrep(bart_path, filesep, '/'), ...
                         '"', '/bart ', strrep(cmd, filesep, '/'), ' ', ...
                 strrep(in_str, filesep, '/'), ...
                         ' ', strrep(out_str, filesep, '/'), '"']);
+        end
 	else
 		ERR = system([bart_path, '/bart ', cmd, ' ', in_str, ' ', out_str]);
 	end
 
-	for i=1:nargin - 1,
+	for i=1:nargin - 1
 		if (exist(strcat(in{i}, '.cfl'),'file'))
 			delete(strcat(in{i}, '.cfl'));
 		end
@@ -71,7 +88,7 @@ function [varargout] = bart(cmd, varargin);
 		end
 	end
 
-	for i=1:nargout,
+	for i=1:nargout
 		if ERR==0
 			varargout{i} = readcfl(out{i});
 		end

python/bart.py

@@ -4,31 +4,42 @@
 #
 # Authors: 
 # 2016 Siddharth Iyer <sid8795@gmail.com>
+# 2018 Soumick Chatterjee <soumick.chatterjee@ovgu.de> , WSL Support
 
 import subprocess as sp
 import tempfile as tmp
 import cfl
 import os
+from wslsupport import PathCorrection
 
 def bart(nargout, cmd, *args):
 
     if type(nargout) != int or nargout < 0:
-        print("Usage: bart(<nargout>, <command>, <arguements...>)");
+        print("Usage: bart(<nargout>, <command>, <arguements...>)")
         return None
 
-    bart_path = os.environ['TOOLBOX_PATH'] + '/bart ';
+    try:
+        bart_path = os.environ['TOOLBOX_PATH'] + '/bart '
+    except:
+        bart_path = None
+    isWSL = False
 
     if not bart_path:
         if os.path.isfile('/usr/local/bin/bart'):
             bart_path = '/usr/local/bin'
         elif os.path.isfile('/usr/bin/bart'):
             bart_path = '/usr/bin'
+        else:
+            bartstatus = os.system('wsl bart version -V')
+            if bartstatus==0:
+                bart_path = '/usr/bin'
+                isWSL = True
             else:
                 raise Exception('Environment variable TOOLBOX_PATH is not set.')
 
     name = tmp.NamedTemporaryFile().name
 
-    nargin = len(args);
+    nargin = len(args)
     infiles = [name + 'in' + str(idx) for idx in range(nargin)]
     in_str = ' '.join(infiles)
 
@@ -38,8 +49,19 @@ def bart(nargout, cmd, *args):
     outfiles = [name + 'out' + str(idx) for idx in range(nargout)]
     out_str = ' '.join(outfiles)
 
-    #TODO: Windows option.
-    ERR = os.system(bart_path + '/bart ' + cmd + ' ' + in_str + ' ' + out_str);
+    if os.name =='nt':
+        if isWSL:
+            #For WSL and modify paths
+            cmdWSL = PathCorrection(cmd)
+            in_strWSL = PathCorrection(in_str)
+            out_strWSL =  PathCorrection(out_str)	
+            ERR = os.system('wsl bart ' + cmdWSL + ' ' + in_strWSL + ' ' + out_strWSL)
+        else:
+            #For cygwin use bash and modify paths
+            ERR = os.system('bash.exe --login -c ' + bart_path + '"/bart ' + cmd.replace(os.path.sep, '/') + ' ' + in_str.replace(os.path.sep, '/') + ' ' + out_str.replace(os.path.sep, '/') + '"')
+            #TODO: Test with cygwin, this is just translation from matlab code
+    else:
+        ERR = os.system(bart_path + '/bart ' + cmd + ' ' + in_str + ' ' + out_str)
 
     for elm in infiles:
         if os.path.isfile(elm + '.cfl'):

python/wslsupport.py

+import string
+import os
+
+def PathCorrection(inData):
+	outData=inData
+	for i in string.ascii_lowercase: #Replace drive letters with /mnt/
+		outData=outData.replace(i+':','/mnt/'+i) #if drive letter is supplied in lowercase
+		outData=outData.replace(i.upper()+':','/mnt/'+i) #if drive letter is supplied as uppercase
+	outData=outData.replace(os.path.sep, '/') #Change windows filesep to linux filesep
+
+	return outData

rules/noncart.mk

@@ -12,4 +12,6 @@ noncartobjs := $(noncartsrcs:.c=.o)
 lib/libnoncart.a: libnoncart.a($(noncartobjs))
 
 
+UTARGETS += test_nufft
+MODULES_test_nufft += -lnoncart -llinops
 

rules/simu.mk

@@ -13,6 +13,7 @@ simuobjs := $(simusrcs:.c=.o)
 lib/libsimu.a: libsimu.a($(simuobjs))
 
 
-UTARGETS += test_biot_savart
+UTARGETS += test_ode_bloch test_biot_savart
+MODULES_test_ode_bloch += -lsimu
 MODULES_test_biot_savart += -lsimu
 

src/calib/cc.c

@@ -191,6 +191,8 @@ void gcc(const long out_dims[DIMS], complex float* out_data, const long caldims[
 		md_free(out2);
 		md_free(tmp2);
 	}
+
+	md_free(tmp);
 }
 
 

src/estdelay.c

@@ -425,6 +425,8 @@ int main_estdelay(int argc, char* argv[])
 	for (unsigned int i = 0; i < N; i++)
 		angles[i] = M_PI + atan2f(crealf(traj1[3 * i + 0]), crealf(traj1[3 * i + 1]));
 
+	md_free(traj1);
+
 
 	long full_dims[DIMS];
 	const complex float* full_in = load_cfl(argv[2], DIMS, full_dims);
@@ -475,6 +477,7 @@ int main_estdelay(int argc, char* argv[])
 		complex float* im_pad = md_alloc(DIMS, pad_dims, CFL_SIZE);
 
 		md_resize_center(DIMS, pad_dims, im_pad, dims, im, CFL_SIZE);
+		md_free(im);
 
 		// Sinc filter in k-space (= crop FOV in image space)
 		long crop_dims[DIMS];
@@ -489,8 +492,11 @@ int main_estdelay(int argc, char* argv[])
 
 		complex float* im_pad2 = md_alloc(DIMS, pad_dims, CFL_SIZE);
 		md_zmul(DIMS, pad_dims, im_pad2, im_pad, mask);
+		md_free(im_pad);
+		md_free(mask);
 
 		fftuc(DIMS, pad_dims, PHS1_FLAG, k_pad, im_pad2);
+		md_free(im_pad2);
 
 		//--- Consider only center region ---
 
@@ -503,6 +509,7 @@ int main_estdelay(int argc, char* argv[])
 		long pos[DIMS] = { 0 };
 		pos[PHS1_DIM] = pad_dims[PHS1_DIM]/2 - (c_region/2);
 		md_copy_block(DIMS, pos, kc_dims, kc, pad_dims, k_pad, CFL_SIZE);
+		md_free(k_pad);
 
 		//--- Calculate intersections ---
 
@@ -516,6 +523,7 @@ int main_estdelay(int argc, char* argv[])
 		check_intersections(Nint, N, S, angles, idx, c_region);
 
 		bart_printf("%f:%f:%f\n\n", creal(S[0][0]) / pad_factor, creal(S[1][0]) / pad_factor, creal(S[2][0]) / pad_factor);
+		md_free(kc);
 	}
 
 	unmap_cfl(DIMS, full_dims, full_in);

src/estdims.c

@@ -17,9 +17,12 @@
 #include "misc/mmio.h"
 #include "misc/io.h"
 #include "misc/misc.h"
+
 #include "noncart/nufft.h"
 
 
+#define FFT_DIMS (MD_BIT(0)|MD_BIT(1)|MD_BIT(2))
+
 static const char usage_str[] = "<traj>";
 static const char help_str[] = "Estimate image dimension from non-Cartesian trajectory.\n"
 			"Assume trajectory scaled to -DIM/2 to DIM/2 (ie dk=1/FOV=1)\n";
@@ -40,7 +43,7 @@ int main_estdims(int argc, char* argv[])
 
 	long im_dims[N];
 	
-	estimate_im_dims(N, im_dims, traj_dims, traj);
+	estimate_im_dims(N, FFT_DIMS, im_dims, traj_dims, traj);
 
 	bart_printf("%ld %ld %ld\n", im_dims[0], im_dims[1], im_dims[2]);
 	

src/grecon/italgo.c

@@ -12,6 +12,7 @@
 
 #include <stdbool.h>
 #include <assert.h>
+#include <math.h>
 
 #include "misc/debug.h"
 #include "misc/types.h"
@@ -155,8 +156,12 @@ struct iter italgo_config(enum algo_t algo, int nr_penalties, const struct reg_s
 			*pdconf = iter_chambolle_pock_defaults;
 
 			pdconf->maxiter = maxiter;
-			pdconf->sigma = 1. * scaling;
-			pdconf->tau = 1. / pdconf->sigma;
+			pdconf->sigma = sqrtf(step);
+			pdconf->tau = sqrtf(step);
+
+			pdconf->sigma *= scaling;
+			pdconf->tau /= scaling;
+
 			pdconf->theta = 1.;
 			pdconf->decay = (hogwild ? .95 : 1.);
 			pdconf->tol = 1.E-4;

src/index.c

+/* 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: 
+ * 2017 Martin Uecker <martin.uecker@med.uni-goettingen.de>
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <complex.h>
+#include <string.h>
+
+#include "num/multind.h"
+#include "num/flpmath.h"
+#include "num/init.h"
+
+#include "misc/mmio.h"
+#include "misc/io.h"
+#include "misc/misc.h"
+
+
+static const char usage_str[] = "dim size name";
+static const char help_str[] = "Create an array counting from 0 to {size-1} in dimensions {dim}.\n";
+
+
+
+int main_index(int argc, char* argv[])
+{
+	mini_cmdline(&argc, argv, 3, usage_str, help_str);
+
+	num_init();
+
+	int N = atoi(argv[1]);
+	int s = atoi(argv[2]);
+
+	assert(N >= 0);
+	assert(s >= 0);
+
+	long dims[N + 1];
+
+	for (int i = 0; i < N; i++)
+		dims[i] = 1;
+
+	dims[N] = s;
+
+	complex float* x = create_cfl(argv[3], N + 1, dims);
+
+	for (int i = 0; i < s; i++)
+		x[i] = i;
+
+	unmap_cfl(N + 1, dims, x);
+	exit(0);
+}
+
+

src/iter/italgos.c

 /* Copyright 2013-2017. The Regents of the University of California.
- * Copyright 2016-2017. Martin Uecker.
+ * Copyright 2016-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:
- * 2012-2017 Martin Uecker <martin.uecker@med.uni-goettingen.de>
+ * 2012-2018 Martin Uecker <martin.uecker@med.uni-goettingen.de>
  * 2013-2014 Frank Ong <frankong@berkeley.edu>
  * 2013-2014,2017 Jon Tamir <jtamir@eecs.berkeley.edu>
  *
@@ -420,6 +420,7 @@ float conjgrad(unsigned int maxiter, float l2lambda, float epsilon,
 		debug_printf(DP_DEBUG3, "CG: early out\n");
 		goto cleanup;
 	}
+
 	for (i = 0; i < maxiter; i++) {
 
 		iter_monitor(monitor, vops, x);
@@ -439,24 +440,23 @@ float conjgrad(unsigned int maxiter, float l2lambda, float epsilon,
 		vops->axpy(N, x, +alpha, p);
 		vops->axpy(N, r, -alpha, Ap);
 	
-		rsnew = (float)pow(vops->norm(N, r), 2.);
+		rsnew = pow(vops->norm(N, r), 2.);
+
 		float beta = rsnew / rsold;
 		
 		rsold = rsnew;
 
-		if (rsnew <= eps_squared) {
-			//debug_printf(DP_DEBUG3, "%d ", i);
+		if (rsnew <= eps_squared)
 			break;
-		}
 
 		vops->xpay(N, beta, p, r);	// p = beta * p + r
-
 	}
 
 cleanup:
 	vops->del(Ap);
 	vops->del(p);
 	vops->del(r);
+
 	debug_printf(DP_DEBUG2, "\t cg: %3d\n", i);
 
 	return sqrtf(rsnew);
@@ -507,7 +507,9 @@ void irgnm(unsigned int iter, float alpha, float alpha_min, float redu, long N,
 		iter_op_p_call(inv, alpha, h, p);
 
 		vops->axpy(N, x, 1., h);
+
 		alpha = (alpha - alpha_min) / redu + alpha_min;
+
 		if (NULL != callback.fun)
 			iter_op_call(callback, x, x);
 	}
@@ -607,6 +609,7 @@ double power(unsigned int maxiter,
 	for (unsigned int i = 0; i < maxiter; i++) {
 
 		iter_op_call(op, u, u);		// r = A x
+
 		s = vops->norm(N, u);
 		vops->smul(N, 1. / s, u, u);
 	}
@@ -663,7 +666,6 @@ void chambolle_pock(unsigned int maxiter, float epsilon, float tau, float sigma,
 	vops->clear(M, u_old);
 
 
-
 	for (unsigned int i = 0; i < maxiter; i++) {
 
 		float lambda = (float)pow(decay, i);
@@ -676,8 +678,11 @@ void chambolle_pock(unsigned int maxiter, float epsilon, float tau, float sigma,
 		 */
 
 		iter_op_call(op_forw, u_old, x_avg);
+
 		vops->axpy(M, u_old, 1. / sigma, u); // (u + sigma * A(x)) / sigma
+
 		iter_op_p_call(prox1, 1. / sigma, u_new, u_old);
+
 		vops->axpbz(M, u_new, -1. * sigma, u_new, sigma, u_old);
 		vops->copy(M, u_old, u);
 		vops->axpbz(M, u, lambda, u_new, 1. - lambda, u_old);
@@ -689,9 +694,13 @@ void chambolle_pock(unsigned int maxiter, float epsilon, float tau, float sigma,
 		 * x = lambda * x + (1 - lambda * x0)
 		 */
 		vops->copy(N, x_old, x);
+
 		iter_op_call(op_adj, x_new, u);
+
 		vops->axpy(N, x, -1. * tau, x_new);
+
 		iter_op_p_call(prox2, tau, x_new, x);
+
 		vops->axpbz(N, x, lambda, x_new, 1. - lambda, x_old);
 
 		/* update x_avg
@@ -712,19 +721,6 @@ void chambolle_pock(unsigned int maxiter, float epsilon, float tau, float sigma,
 
 		if (epsilon > (res1 + res2))
 			break;
-
-#if 0 // buggy
-		if (res1 < 100 * res2) {
-
-			sigma /= 2;
-			tau *= 2;
-		}
-		else if (res2 > 100 * res1) {
-
-			sigma *= 2;
-			tau /= 2;
-		}
-#endif
 	}
 
 	debug_printf(DP_DEBUG3, "\n");