Commit 3f8d82ca authored by Martin Uecker's avatar Martin Uecker

Merge tag 'v0.5.00'

version 0.5.00

Changes:
	- new tools: moba, looklocker, tgv, extract, fftrot
	- moba: non-linear model-based recon for T1 mapping (Xiaoqing Wang)
	- looklocker: post-processing for moba
	- extract: extract ranges from multi-dim. arrays
	- tgv: a tgv denoiser
	- fftrot: perform rotations using FFTs
	- join: append to non-existing file
	- pics: low-mem flag (-U)
	- pics: efficient basis functions (Sebastian Rosenzweig)
	- pics: pattern for non-Cartesian trajectories
	- nlinv: enlive on the gpu
	- twixread: read multi-raid files
	- twixread: mode for radial
	- traj: custom angle to traj (Aurélien Trotier)
	- traj: small golden angles (Jasper Schoormans)
	- traj: option to rotate trajectory
	- extract: multiple dims (Christian Holme)
	- slice: multiple dims (Christian Holme)
	- generic: NOEXEC_FLAG flag for non-executable stacks
	- generic: automatic rebuilds for makefile changes
	- generic: make it possible to run individual tests
	- library: use long in gridding code to support larger grids
	- library: zsmax functions for cpu and gpu
	- library: function for setting random initial state in wavthresh
	- library: lsqr/lad return operator_p to select regularization
	- library: generalize cuda fft a bit for efficiency
	- library: fix for CUDA strided memcpy
	- many other bug fixes and improvements
parents 3bc0e34e 5efc0a95
......@@ -59,7 +59,6 @@ matrix:
- liblapacke-dev
- libblocksruntime-dev
script:
- make test
- make utest
- name: "Makefile OpenMP"
......
Copyright (c) 2013-2018. The Regents of the University of California.
Copyright (c) 2013-2018. BART Developer Team and Contributors.
Copyright (c) 2013-2019. BART Developer Team and Contributors.
Copyright (c) 2012. Intel Corporation. (src/lapacke/)
All rights reserved.
......
# Copyright 2013-2015. The Regents of the University of California.
# Copyright 2015-2018. Martin Uecker <martin.uecker@med.uni-goettingen.de>
# Copyright 2015-2019. Martin Uecker <martin.uecker@med.uni-goettingen.de>
# All rights reserved. Use of this source code is governed by
# a BSD-style license which can be found in the LICENSE file.
......@@ -9,6 +9,9 @@ MAKESTAGE ?= 1
# silent make
#MAKEFLAGS += --silent
# auto clean on makefile updates
AUTOCLEAN?=1
# clear out all implicit rules and variables
MAKEFLAGS += -R
......@@ -23,6 +26,7 @@ SLINK?=0
DEBUG?=0
FFTWTHREADS?=1
ISMRMRD?=0
NOEXEC_STACK?=0
LOG_BACKEND?=0
LOG_SIEMENS_BACKEND?=0
......@@ -165,6 +169,7 @@ MODULES_pics = -lgrecon -lsense -liter -llinops -lwavelet -llowrank -lnoncart
MODULES_sqpics = -lsense -liter -llinops -lwavelet -llowrank -lnoncart
MODULES_pocsense = -lsense -liter -llinops -lwavelet
MODULES_nlinv = -lnoir -liter -lnlops -llinops
MODULES_moba = -lmoba -lnoir -liter -lnlops -llinops -lwavelet
MODULES_bpsense = -lsense -lnoncart -liter -llinops -lwavelet
MODULES_itsense = -liter -llinops
MODULES_ecalib = -lcalib
......@@ -177,21 +182,24 @@ MODULES_ccapply = -lcalib
MODULES_estvar = -lcalib
MODULES_nufft = -lnoncart -liter -llinops
MODULES_rof = -liter -llinops
MODULES_tgv = -liter -llinops
MODULES_bench = -lwavelet -llinops
MODULES_phantom = -lsimu
MODULES_bart = -lbox -lgrecon -lsense -lnoir -liter -llinops -lwavelet -llowrank -lnoncart -lcalib -lsimu -lsake -ldfwavelet -lnlops
MODULES_bart = -lbox -lgrecon -lsense -lnoir -liter -llinops -lwavelet -llowrank -lnoncart -lcalib -lsimu -lsake -ldfwavelet -lnlops -lmoba
MODULES_sake = -lsake
MODULES_wave = -liter -lwavelet -llinops -lsense
MODULES_traj = -lnoncart
MODULES_wave = -liter -lwavelet -llinops -llowrank
MODULES_threshold = -llowrank -liter -ldfwavelet -llinops -lwavelet
MODULES_fakeksp = -lsense -llinops
MODULES_lrmatrix = -llowrank -liter -llinops
MODULES_estdims = -lnoncart -llinops
MODULES_ismrmrd = -lismrm
MODULES_wavelet = -llinops -lwavelet
MODULES_wshfl = -llinops -lwavelet -liter -llowrank -llinops
MODULES_wshfl = -llinops -lwavelet -liter -llowrank
MAKEFILES = $(root)/Makefiles/Makefile.*
MAKEFILES = $(wildcard $(root)/Makefiles/Makefile.*)
ALLMAKEFILES = $(root)/Makefile $(wildcard $(root)/Makefile.* $(root)/*.mk $(root)/rules/*.mk $(root)/Makefiles/Makefile.*)
-include Makefile.$(NNAME)
-include Makefile.local
......@@ -233,6 +241,9 @@ CPPFLAGS += -g
CFLAGS += -g
endif
ifeq ($(NOEXEC_STACK),1)
CPPFLAGS += -DNOEXEC_STACK
endif
ifeq ($(PARALLEL),1)
MAKEFLAGS += -j
......@@ -243,6 +254,12 @@ ifeq ($(MAKESTAGE),1)
.PHONY: doc/commands.txt $(TARGETS)
default all clean allclean distclean doc/commands.txt doxygen test utest gputest $(TARGETS):
make MAKESTAGE=2 $(MAKECMDGOALS)
tests/test-%: force
make MAKESTAGE=2 $(MAKECMDGOALS)
force: ;
else
......@@ -506,20 +523,39 @@ $(UTARGETS): % : utests/utest.c utests/%.o $$(MODULES_%) $(MODULES)
clean:
rm -f `find $(srcdir) -name "*.o"`
rm -f utests/*.o
rm -f $(root)/utests/*.o
rm -f $(patsubst %, %, $(UTARGETS))
rm -f $(root)/lib/.*.lock
rm -f $(libdir)/.*.lock
allclean: clean
rm -f $(libdir)/*.a $(ALLDEPS)
rm -f $(patsubst %, %, $(TARGETS))
rm -f $(srcdir)/misc/version.inc
rm -rf doc/dx
rm -f doc/commands.txt
rm -rf $(root)/tests/tmp/*/
rm -rf $(root)/doc/dx
rm -f $(root)/doc/commands.txt
touch isclean
distclean: allclean
-include isclean
isclean: $(ALLMAKEFILES)
ifeq ($(AUTOCLEAN),1)
@echo "CONFIGURATION MODIFIED. RUNNING FULL REBUILD."
touch isclean
make allclean || rm isclean
else
ifneq ($(MAKECMDGOALS),allclean)
@echo "CONFIGURATION MODIFIED."
endif
endif
# automatic tests
# system tests
......
......@@ -72,7 +72,7 @@ The software tools in recon should run on any recent Linux distribution.
To install the required libraries on Debian and Ubuntu run:
$ sudo apt-get install gcc make libfftw3-dev liblapacke-dev libpng-dev libopenblas-dev
$ sudo apt-get install gcc make libfftw3-dev liblapacke-dev libpng-dev libopenblas-dev gfortran
(optional)
$ sudo apt-get install octave
......@@ -110,7 +110,7 @@ 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. If you use WSL, you need to
use the clang compiler.
either use the clang compiler or gcc with the NOEXEC_STACK option.
BART should also work with Cygwin:
......
......@@ -2,9 +2,10 @@
#
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
TNUM=fft fftmod fftshift noise bench threshold conv rss filter mandelbrot wavelet window var std fftrot
TRECO=pics pocsense sqpics itsense nlinv moba nufft rof tgv sake wave lrmatrix estdims estshift estdelay wavepsf wshfl
TCALIB=ecalib ecaltwo caldir walsh cc ccapply calmat svd estvar whiten
TMRI=homodyne poisson twixread fakeksp
TMRI=homodyne poisson twixread fakeksp looklocker
TSIM=phantom traj
TIO=toimg
(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
chemical shift imaging in muscular dystrophy with compressed sensing
and parallel imaging reconstruction.
Magn Reson Med 2014; 72:1610-1619.
Zhang T, Cheng JY, Potnick AG, Barth RA, Alley MT, Uecker M, Lustig M,
Pauly JM, Vasanawala SS.
Fast Pediatric 3D Free Breathing Abdominal Dynamic Contrast Enhanced MRI
with a High Spatiotemporal Resolution,
J Magn Reson Imaging 2015; 41:460-473.
Addy NO, Ingle RR, Wu HH, Hu BS, Nishimura DG.
High-resolution variable-density 3D cones coronary MRA.
Magn Reson Med 2015; 74:614-621.
Cheng JY, Zhang T, Ruangwattanapaisarn N, Alley MT, Uecker M, Pauly JM,
Lustig M, Vasanawala SS.
Free-Breathing Pediatric MRI with Nonrigid Motion Correction and Acceleration,
J Magn Reson Imaging 2015; 42:407-420.
Athalye V, Lustig M, Uecker M.
Parallel Magnetic Resonance Imaging as Approximation in a
Reproducing Kernel Hilbert Space,
Inverse Problems 2015; 31:045008.
Mann LW, Higgins DM, Peters CN, Cassidy S, Hodson KK, Coombs A,
Taylor R, Hollingsworth KG.
Accelerating MR Imaging Liver Steatosis Measurement Using Combined
Compressed Sensing and Parallel Imaging: A Quantitative Evaluation,
Radiology 2016; 278:245-256.
Cheng JY, Hanneman K, Zhang T, Alley MT, Lai P, Tamir JI, Uecker M,
Pauly JM, Lustig M, Vasanawala SS.
Comprehensive Motion-Compensated Highly-Accelerated 4D Flow MRI with
Ferumoxytol Enhancement for Pediatric Congenital Heart Disease.
J Magn Reson Imaging 2016; 43:1355-1368.
Tamir JI, Uecker M, Chen W, Lai P, Aleey MT, Vasanawala SS, Lustig M.
T2-Shuffling: Sharp, Multi-Contrast, Volumetric Fast Spin-Echo Imaging.
Magn Recon Med 2017; 77:180-195.
Uecker M, Lustig M.
Estimating Absolute-Phase Maps Using ESPIRiT and Virtual Conjugate Coils.
Magn Reson Med 2017; 77:1201-1207.
Cheng JY, Zhang T, Alley MT, Uecker M, Lustig M, Pauly JM, Vasanawala SS.
Comprehensive Multi-Dimensional MRI for the Simultaneous Assessment
of Cardiopulmonary Anatomy and Physiology.
Scientific Reports 2017; 7:5330.
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.
J Magn Reson Imaging 2017; 45:1700-1711.
Mazzoli V, Schoormans J, Froeling M, Sprengers AM, Coolen BF, Verdonschot N,
Strijkers GJ, Nederveen AJ.
Accelerated 4D self‐gated MRI of tibiofemoral kinematics.
NMR in Biomed 2017; 30:e3791.
Moghari MH, Uecker M, Roujol S, Sabbagh M, Geva T, Powell AJ.
Accelerated Whole-heart Magnetic Resonance Angiography Using
a Variable-Density Poisson-Disc Undersampling Pattern and
Compressed Sensing Reconstruction.
Magn Reson Med 2018; 79:761-769.
Peper ES, Strijkers GJ, Gazzola K, Potters WV, Motaal AG, Luirink IK, Hutten BA,
Wiegman A, van Ooij P, van den Born B-JH, Nederveen AJ, Coolen BF.
Regional assessment of carotid artery pulse wave velocity using compressed
sensing accelerated high temporal resolution 2D CINE PC MRI.
J Cardiovasc Magn Reson 2018; 20:1-12.
Mazzoli V, Gottwald LM, Peper ES, Froeling M, Coolen BF, Verdonschot N,
Sprengers AM, van Ooij P, Strijkers GJ, Nederveen AJ.
Accelerated 4D phase contrast MRI in skeletal muscles contraction.
Magn Reson Med 2018; 80:1799-1811.
Rosenzweig S, Holme HCM, Wilke RN, Voit D, Frahm J, Uecker M.
Simultaneous Multi-Slice Reconstruction Using Regularized Nonlinear Inversion:
SMS-NLINV.
Magn Reson Med 2018; 79:2057-2066.
Sanders J, Song H, Frank S, Bathala T, Venkatesan A, Anscher M, Tang C, Bruno T,
Wei W, Ma J.
Parallel imaging compressed sensing for accelerated imaging and improved SNR in
MRI-based prostate brachytherapy post-implant dosimetry.
Brachytherapy 2018; 17:816-824.
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.
Magn Reson Med 2019; 81:1566-1579.
de Jonge CS, Coolen BF, Peper ES, Motaal AG, Nio CY, Somers I, Strijkers GJ,
Stoker J, Nederveen AJ.
Evaluation of compressed sensing MRI for accelerated bowel motility imaging.
European Radiology Experimental 2019; 3:1-12.
Walheim J, Dillinger H, Kozerke.
Multipoint 5D flow cardiovascular magnetic resonance - accelerated cardiac-
and respiratory-motion resolved mapping of mean and turbulent velocities.
J Cardiovasc Magn Reson 2019; 21:42.
Wang X, Kohler F, Unterberg-Buchwald C, Lotz J, Frahm J, Uecker M.
Model-based myocardial T1 mapping with sparsity constraints using single-shot
inversion-recovery radial FLASH.
J Cardiovasc Magn Reson 2019; in press.
This diff is collapsed.
......@@ -12,23 +12,24 @@ function [varargout] = bart(cmd, varargin)
end
bart_path = getenv('TOOLBOX_PATH');
isWSL = false;
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
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
end
end
% clear the LD_LIBRARY_PATH environment variable (to work around
......@@ -60,20 +61,20 @@ function [varargout] = bart(cmd, varargin)
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
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
......
mobasrcs := $(wildcard $(srcdir)/moba/*.c)
mobaobjs := $(mobasrcs:.c=.o)
.INTERMEDIATE: $(mobaobjs)
lib/libmoba.a: libmoba.a($(mobaobjs))
UTARGETS += test_moba
MODULES_test_moba += -lmoba -lnlops -llinops
......@@ -18,5 +18,5 @@ lib/libnum.a: libnum.a($(numobjs))
UTARGETS += test_multind test_flpmath test_splines test_linalg test_polynom test_window
UTARGETS += test_blas test_mdfft
UTARGETS += test_blas test_mdfft test_ops test_ops_p
......@@ -159,7 +159,7 @@ calib_slice()
bart slice 2 $1 grasp_hybrid grasp1-$1
# extract first $CALIB spokes
bart extract 1 $(($SKIP + 0)) $(($SKIP + $CALIB - 1)) grasp1-$1 grasp2-$1
bart extract 1 $(($SKIP + 0)) $(($SKIP + $CALIB)) grasp1-$1 grasp2-$1
# reshape dimensions
bart reshape $(bart bitmask 0 1 2 3) 1 $READ $CALIB $COILS grasp2-$1 grasp3-$1
......@@ -174,7 +174,7 @@ recon_slice()
bart slice 2 $1 sens sens-$1
# extract spokes and split-off time dim
bart extract 1 $(($SKIP + 0)) $(($SKIP + $SPOKES * $PHASES - 1)) grasp1-$1 grasp2-$1
bart extract 1 $(($SKIP + 0)) $(($SKIP + $SPOKES * $PHASES)) grasp1-$1 grasp2-$1
bart reshape $(bart bitmask 1 2) $SPOKES $PHASES grasp2-$1 grasp1-$1
# move time dimensions to dim 10 and reshape
......
......@@ -18,7 +18,7 @@
#include "num/flpmath.h"
#include "num/rand.h"
#include "num/init.h"
#include "num/ops.h"
#include "num/ops_p.h"
#include "num/mdfft.h"
#include "num/fft.h"
......
......@@ -150,17 +150,13 @@ void crop_sens(const long dims[DIMS], complex float* ptr, bool soft, float crth,
*/
static float sure_crop(float var, const long evec_dims[5], complex float* evec_data, complex float* eptr, const long calreg_dims[4], const complex float* calreg)
{
// Must be in ascending order.
float start = 0.7;
float delta = 0.01;
long num_cvals = (long) (((1 - delta - start)/delta) + 1);
long num_maps = evec_dims[4];
// Construct low-resolution image
long im_dims[5];
md_select_dims(5, 15, im_dims, evec_dims);
complex float* im = md_alloc(5, im_dims, CFL_SIZE);
complex float* im = md_alloc_sameplace(5, im_dims, CFL_SIZE, calreg);
md_clear(5, im_dims, im, CFL_SIZE);
md_resize_center(5, im_dims, im, calreg_dims, calreg, CFL_SIZE);
ifftuc(5, im_dims, FFT_FLAGS, im, im);
......@@ -170,31 +166,31 @@ static float sure_crop(float var, const long evec_dims[5], complex float* evec_d
cropdims[4] = num_maps;
// Eigenvectors (M)
complex float* M = md_alloc(5, evec_dims, CFL_SIZE);
complex float* M = md_alloc_sameplace(5, evec_dims, CFL_SIZE, calreg);
md_copy(5, evec_dims, M, evec_data, CFL_SIZE);
// Temporary eigenvector holder to hold low resolution maps
complex float* LM = md_alloc(5, evec_dims, CFL_SIZE);
complex float* LM = md_alloc_sameplace(5, evec_dims, CFL_SIZE, calreg);
// Temporary holder for projection calreg
complex float* TC = md_alloc(5, calreg_dims, CFL_SIZE);
complex float* TC = md_alloc_sameplace(5, calreg_dims, CFL_SIZE, calreg);
// Temporary holder to hold low resolution calib maps
complex float* CM = md_alloc(5, cropdims, CFL_SIZE);
complex float* CM = md_alloc_sameplace(5, cropdims, CFL_SIZE, calreg);
// Eigenvalues (W)
long W_dims[5];
md_select_dims(5, 23, W_dims, evec_dims);
complex float* W = md_alloc(5, W_dims, CFL_SIZE);
complex float* W = md_alloc_sameplace(5, W_dims, CFL_SIZE, calreg);
md_copy(5, W_dims, W, eptr, CFL_SIZE);
// Place holder for the inner product result
complex float* ip = md_alloc(5, W_dims, CFL_SIZE);
complex float* ip = md_alloc_sameplace(5, W_dims, CFL_SIZE, calreg);
// Place holder for the projection result
complex float* proj = md_alloc(5, im_dims, CFL_SIZE);
complex float* proj = md_alloc_sameplace(5, im_dims, CFL_SIZE, calreg);
// Place holder for divergence term
long div_dims[5] = MD_INIT_ARRAY(5, 1);
complex float* div = md_alloc(5, div_dims, CFL_SIZE);
complex float* div = md_alloc_sameplace(5, div_dims, CFL_SIZE, calreg);
// Calculating strides.
long str1_ip[5];
......@@ -211,7 +207,7 @@ static float sure_crop(float var, const long evec_dims[5], complex float* evec_d
md_calc_strides(5, str1_proj, W_dims, CFL_SIZE);
md_calc_strides(5, str2_proj, evec_dims, CFL_SIZE);
md_calc_strides(5, stro_proj, im_dims, CFL_SIZE);
md_calc_strides(5, stro_proj, im_dims, CFL_SIZE);
long str1_div[5];
long str2_div[5];
......@@ -225,7 +221,6 @@ static float sure_crop(float var, const long evec_dims[5], complex float* evec_d
long tdims_proj[5];
for (unsigned int i = 0; i < 5; i++) {
assert((im_dims[i] == evec_dims[i]) || (1 == im_dims[i]) || (1 == evec_dims[i]));
assert(( W_dims[i] == evec_dims[i]) || (1 == W_dims[i]) || (1 == evec_dims[i]));
......@@ -234,56 +229,63 @@ static float sure_crop(float var, const long evec_dims[5], complex float* evec_d
}
// Starting parameter sweep with SURE.
float minMSE = 0;
float estMSE = 0;
float optVal = 0;
float c = 0;
for (int idx = 0; idx < num_cvals; idx++) {
estMSE = 0;
*div = 0;
md_clear(5, W_dims, ip, CFL_SIZE);
md_clear(5, im_dims, proj, CFL_SIZE);
md_clear(5, evec_dims, LM, CFL_SIZE);
md_clear(5, calreg_dims, TC, CFL_SIZE);
c = start + idx * delta;
// Cropping
crop_weight(evec_dims, M, crop_thresh_function, c, W);
float mse = -1;
float old_mse = 0;
// Projection (stored in proj)
md_zfmacc2(5, tdims_ip, stro_ip, ip, str1_ip, im, str2_ip, M);
md_zfmac2 (5, tdims_proj, stro_proj, proj, str1_proj, ip, str2_proj, M);
// Construct low resolution projection image.
fftuc(5, im_dims, FFT_FLAGS, proj, proj);
md_resize_center(5, calreg_dims, TC, im_dims, proj, CFL_SIZE);
md_resize_center(5, im_dims, proj, calreg_dims, TC, CFL_SIZE);
ifftuc(5, im_dims, FFT_FLAGS, proj, proj);
for (int jdx = 0; jdx < md_calc_size(5, im_dims); jdx++)
estMSE += powf(cabsf(im[jdx] - proj[jdx]), 2);
// Construct low-resolution maps
fftuc(5, evec_dims, FFT_FLAGS, LM, M);
md_resize_center(5, cropdims, CM, evec_dims, LM, CFL_SIZE);
md_resize_center(5, evec_dims, LM, cropdims, CM, CFL_SIZE);
ifftuc(5, evec_dims, FFT_FLAGS, LM, LM);
// Calculating SURE divergence using low resolution maps
md_zfmacc2(5, evec_dims, stro_div, div, str1_div, LM, str2_div, LM);
estMSE += 2 * var * creal(*div);
if ((0 == idx) || (estMSE < minMSE)) {
float s = -0.1;
float c = 0.99;
long ctr1 = 0;
long ctr2 = 0;
optVal = c;
minMSE = estMSE;
debug_printf(DP_INFO, "---------------------------------------------\n");
debug_printf(DP_INFO, "| CTR1 | CTR2 | Crop | Est. MSE |\n");
debug_printf(DP_INFO, "---------------------------------------------\n");
while (fabs(s) > 1E-4) {
ctr1 ++;
while (c < 0.999 && c > 0.001 && (ctr2 <= 1 || mse < old_mse)) {
ctr2 ++;
md_clear(5, W_dims, ip, CFL_SIZE);
md_clear(5, im_dims, proj, CFL_SIZE);
md_clear(5, div_dims, div, CFL_SIZE);
md_clear(5, evec_dims, M, CFL_SIZE);
md_clear(5, evec_dims, LM, CFL_SIZE);
md_clear(5, calreg_dims, TC, CFL_SIZE);
md_copy(5, evec_dims, M, evec_data, CFL_SIZE);
old_mse = mse;
mse = 0;
crop_weight(evec_dims, M, crop_thresh_function, c, W); // Cropping.
md_zfmacc2(5, tdims_ip, stro_ip, ip, str1_ip, im, str2_ip, M); // Projection.
md_zfmac2 (5, tdims_proj, stro_proj, proj, str1_proj, ip, str2_proj, M);
fftuc(5, im_dims, FFT_FLAGS, proj, proj); // Low res proj img.
md_resize_center(5, calreg_dims, TC, im_dims, proj, CFL_SIZE);
md_resize_center(5, im_dims, proj, calreg_dims, TC, CFL_SIZE);
ifftuc(5, im_dims, FFT_FLAGS, proj, proj);
for (int jdx = 0; jdx < md_calc_size(5, im_dims); jdx++)
mse += powf((float) (cabsf(im[jdx] - proj[jdx])), 2);
fftuc(5, evec_dims, FFT_FLAGS, LM, M); // low-res maps .
md_resize_center(5, cropdims, CM, evec_dims, LM, CFL_SIZE);
md_resize_center(5, evec_dims, LM, cropdims, CM, CFL_SIZE);
ifftuc(5, evec_dims, FFT_FLAGS, LM, LM);
md_zfmacc2(5, evec_dims, stro_div, div, str1_div, LM, str2_div, LM); // Calc SURE div using low res maps.
mse += (float) (2 * var * crealf(*div));
if (ctr2 == 1)
debug_printf(DP_INFO, "| %4ld | %4ld | %0.4f | %0.12e |\n", ctr1, ctr2, c, mse);
else
debug_printf(DP_INFO, "| | %4ld | %0.4f | %0.12e |\n", ctr2, c, mse);
c = c + s;
}
c -= s;
ctr2 = 0;
s = -s/2;
c += s;
}
c = c + s;
debug_printf(DP_INFO, "---------------------------------------------\n");
md_free(im);
md_free(TC);
......@@ -295,14 +297,9 @@ static float sure_crop(float var, const long evec_dims[5], complex float* evec_d
md_free(proj);
md_free(div);
// Smudge factor is to soften by little bit to improve robustness. This is to account for
// the sweeped thresholds possibly having a too large a step size between them and for any
// other inconsistencies.
float smudge = 0.99;
debug_printf(DP_DEBUG1, "Calculated c: %.3f\n", optVal);
debug_printf(DP_DEBUG1, "Smudge: %.3f\n", smudge);
debug_printf(DP_DEBUG1, "Calculated c: %.4f\n", c);
return smudge * optVal;
return c;
}
......@@ -556,7 +553,7 @@ void calib2(const struct ecalib_conf* conf, const long out_dims[DIMS], complex f
md_zsmul(DIMS, out_dims, out_data, out_data, sqrtf((float)channels));
}
float c = (conf->crop > 0) ? conf->crop : sure_crop(conf->var, out_dims, out_data, eptr, calreg_dims, data);
float c = (conf->crop >= 0.) ? conf->crop : sure_crop(conf->var, out_dims, out_data, eptr, calreg_dims, data);
debug_printf(DP_DEBUG1, "Crop maps... (c = %.2f)\n", c);
......
......@@ -200,7 +200,7 @@ void eigenmapscu(const long dims[5], _Complex float* optr, _Complex float* eptr,
size_t sharedMem = memPerPoint * pointsPerBlock;
eigenmapscu_kern<<<blocks, threads, sharedMem>>>(imgcov2_device_filled, imgcov2_device, optr_device, eptr_device, iter, x, y, z, N, M);
cudaThreadSynchronize();
cudaDeviceSynchronize();
cudaError_t cu_error = cudaGetLastError();
......
This diff is collapsed.
......@@ -97,6 +97,10 @@ int main_cc(int argc, char* argv[])
cal_data = extract_calib(caldims, calsize, in_dims, in_data, false);
}
if (0. == md_znorm(DIMS, caldims, cal_data))
debug_printf(DP_WARN, "Empty calibration region.\n");
if (ECC == cc_type)
debug_printf(DP_WARN, "Warning: ECC depends on a parameter choice rule for optimal results which is not implemented.\n");
......
......@@ -91,7 +91,7 @@ inline _hdev_ double abs(double2 z1) {
} while(0)
#define cuda_sync() do{ \
cuda (ThreadSynchronize()); \
cuda (DeviceSynchronize()); \
cuda (GetLastError()); \
} while(0)
......@@ -106,7 +106,7 @@ inline _hdev_ double abs(double2 z1) {
} while(0)
#define cuda_sync() do{ \
cuda (ThreadSynchronize()); \
cuda (DeviceSynchronize()); \
cuda (GetLastError()); \
} while(0)
......
......@@ -22,6 +22,7 @@
#include "num/multind.h"
#include "num/flpmath.h"
#include "num/ops_p.h"
#include "num/ops.h"
#include "misc/misc.h"
......
<
......@@ -160,7 +160,7 @@ int main_ecalib(int argc, char* argv[])
(conf.usegpu ? num_init_gpu : num_init)();
if ((conf.var < 0) && (conf.weighting || (conf.crop < 0)))
if ((conf.var < 0.) && (conf.weighting || (conf.crop < 0.)))
conf.var = estvar_calreg(conf.kdims, cal_dims, cal_data);