[med-svn] [bart] 01/14: added lrdecom to do multi-scale low rank decomposition
Martin Uecker
uecker-guest at moszumanska.debian.org
Sun Dec 6 23:10:35 UTC 2015
This is an automated email from the git hooks/post-receive script.
uecker-guest pushed a commit to annotated tag v0.2.09d
in repository bart.
commit 53b631b4359b85293784ae17e23ac115ff401a09
Author: Frank Ong <frankong at berkeley.edu>
Date: Sun Nov 29 14:16:19 2015 -0800
added lrdecom to do multi-scale low rank decomposition
---
Makefile | 3 +-
src/lrdecom.c | 300 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 302 insertions(+), 1 deletion(-)
diff --git a/Makefile b/Makefile
index 7c76031..600ccf7 100644
--- a/Makefile
+++ b/Makefile
@@ -107,7 +107,7 @@ ismrm.top ?= /usr/local/ismrmrd/
TBASE=show slice crop resize join transpose zeros ones flip circshift extract repmat bitmask reshape version
TFLP=scale conj fmac saxpy sdot spow cpyphs creal normalize cdf97 relnorm pattern nrmse
TNUM=fft fftmod fftshift noise bench threshold conv rss
-TRECO=pics pocsense rsense bpsense itsense nlinv nufft rof sake wave
+TRECO=pics pocsense rsense bpsense itsense nlinv nufft rof sake wave lrdecom
TCALIB=ecalib ecaltwo caldir walsh cc calmat svd estvar
TMRI=homodyne poisson twixread fakeksp
TSIM=phantom traj
@@ -147,6 +147,7 @@ MODULES_sake += -lsake
MODULES_wave += -liter -lwavelet2 -llinops -lsense
MODULES_threshold += -llowrank -llinops -lwavelet2 -liter -ldfwavelet
MODULES_fakeksp += -lsense -llinops
+MODULES_lrdecom = -llowrank -liter -llinops
-include Makefile.$(NNAME)
-include Makefile.local
diff --git a/src/lrdecom.c b/src/lrdecom.c
new file mode 100644
index 0000000..ab00589
--- /dev/null
+++ b/src/lrdecom.c
@@ -0,0 +1,300 @@
+/* Copyright 2013. The Regents of the University of California.
+ * All rights reserved. Use of this source code is governed by
+ * a educational/research license which can be found in the
+ * LICENSE file.
+ *
+ * Authors:
+ * 2014 Frank Ong <frankong at berkeley.edu>
+ */
+
+
+#define _GNU_SOURCE
+#include <stdlib.h>
+#include <assert.h>
+#include <stdbool.h>
+#include <complex.h>
+#include <stdio.h>
+#include <math.h>
+#include <unistd.h>
+
+#include "num/multind.h"
+#include "num/flpmath.h"
+#include "num/init.h"
+#include "num/ops.h"
+
+#include "linops/linop.h"
+
+#include "iter/iter.h"
+#include "iter/lsqr.h"
+#include "iter/thresh.h"
+
+#include "lowrank/lrthresh.h"
+#include "linops/sum.h"
+#include "iter/prox.h"
+#include "linops/someops.h"
+
+#include "misc/debug.h"
+#include "misc/mri.h"
+#include "misc/mmio.h"
+#include "misc/misc.h"
+
+struct s_data{
+ long size;
+};
+
+// x = (z1 + z2)/2
+
+static void sum_xupdate( const void* _data, float rho, complex float* dst, const complex float* src )
+{
+ UNUSED(rho);
+
+ const struct s_data* data = (const struct s_data*) _data;
+
+ for( int i = 0; i < data->size; i++)
+ dst[i] = src[i] / 2.;
+}
+
+static void sum_xupdate_free( const void* data )
+{
+ free( (void*) data);
+}
+
+
+
+static void usage(const char* name, FILE* fd)
+{
+ fprintf(fd, "Usage: %s [-options] <input> <output>\n", name);
+}
+
+static void help(void)
+{
+ printf( "\n"
+ "Perform multi-scale low rank decomposition.\n"
+ "-i\t\tmaximum iterations.\n"
+ "-m\t\tflags to specify which dimensions are reshaped to matrix columns.\n"
+ "-f\t\tflags to specify which dimensions to perform multi-scale partition.\n"
+ "-j scale\tblock size scaling from one scale to the next one.\n"
+ "-k block-size\tsmallest block size\n"
+ "-N\t\tadd noise scale to account for Gaussian noise.\n"
+ "-s\t\tperform low rank + sparse decomposition.\n"
+ "-l block-size\tperform locally low rank soft thresholding with specified block size.\n"
+ "-o <output2>\tsummed over all non-noise scales to create a denoised output.\n"
+ "\n");
+}
+
+
+int main_lrdecom(int argc, char* argv[])
+{
+ double start_time = timestamp();
+
+ bool use_gpu = false;
+
+ int maxiter = 100;
+ float rho = 0.25;
+ int blkskip = 2;
+ _Bool randshift = true;
+ unsigned long mflags = 1;
+ unsigned long flags = ~0;
+ const char* sum_str = NULL;
+ _Bool noise = false;
+
+ _Bool llr = false;
+ long llrblk = 8;
+ _Bool ls = false;
+ _Bool hogwild = false;
+ _Bool fast = true;
+ long initblk = 1;
+ int remove_mean = 0;
+
+ int c;
+ while (-1 != (c = getopt(argc, argv, "uvNi:p:m:j:k:o:hnl:sf:gHF"))) {
+ switch(c) {
+
+ case 'u':
+ remove_mean = 1;
+ break;
+
+ case 'v':
+ remove_mean = 2;
+ break;
+
+ case 'H':
+ hogwild = true;
+ break;
+
+ case 'k':
+ initblk = atoi(optarg);
+ break;
+
+ case 'o':
+ sum_str = strdup(optarg);
+ break;
+
+
+ case 'i':
+ maxiter = atoi(optarg);
+ break;
+
+ case 'j':
+ blkskip = atoi(optarg);
+ break;
+
+ case 'p':
+ rho = atof(optarg);
+ break;
+
+ case 'N':
+ noise = true;
+ break;
+
+ case 'n':
+ randshift = false;
+ break;
+
+ case 'l':
+ llr = true;
+ llrblk = atoi(optarg);
+ break;
+
+ case 's':
+ ls = true;
+ break;
+
+ case 'f':
+ flags = labs(atol(optarg));
+ break;
+
+ case 'm':
+ mflags = labs(atol(optarg));
+ break;
+
+ case 'g':
+ use_gpu = true;
+ break;
+
+ case 'h':
+ usage(argv[0], stdout);
+ help();
+ exit(0);
+
+ default:
+ usage(argv[0], stderr);
+ exit(1);
+ }
+ }
+
+ if (argc - optind != 2) {
+
+ usage(argv[0], stderr);
+ exit(1);
+ }
+
+
+ long idims[DIMS];
+ long odims[DIMS];
+
+ // Load input
+ complex float* idata = load_cfl(argv[optind + 0], DIMS, idims);
+
+ // Get levels and block dimensions
+ long blkdims[MAX_LEV][DIMS];
+ long levels;
+ if (llr)
+ levels = llr_blkdims(blkdims, flags, idims, llrblk);
+ else if (ls)
+ levels = ls_blkdims(blkdims, idims);
+ else
+ levels = multilr_blkdims(blkdims, flags, idims, blkskip, initblk);
+
+ if (noise)
+ add_lrnoiseblk( &levels, blkdims, idims );
+ debug_printf(DP_INFO, "Number of levels: %ld\n", levels);
+
+ // Get outdims
+ md_copy_dims(DIMS, odims, idims);
+ odims[LEVEL_DIM] = levels;
+ complex float* odata = create_cfl(argv[optind + 1], DIMS, odims);
+ md_clear( DIMS, odims, odata, sizeof(complex float) );
+
+
+ // Initialize algorithm
+ void* iconf;
+
+ struct iter_admm_conf mmconf;
+ memcpy(&mmconf, &iter_admm_defaults, sizeof(struct iter_admm_conf));
+ mmconf.maxiter = maxiter;
+ mmconf.rho = rho;
+ mmconf.hogwild = hogwild;
+ mmconf.fast = fast;
+
+ iconf = &mmconf;
+
+
+ // Initialize operators
+
+ const struct linop_s* sum_op = sum_create( odims, use_gpu );
+ const struct operator_p_s* sum_prox = prox_lineq_create( sum_op, idata );
+ const struct operator_p_s* lr_prox = lrthresh_create(odims, randshift, mflags, (const long (*)[])blkdims, 1., noise, remove_mean, use_gpu);
+
+ assert(use_gpu == false);
+
+ (use_gpu ? num_init_gpu : num_init)();
+
+ if (use_gpu)
+ debug_printf(DP_INFO, "GPU reconstruction\n");
+
+ // put into iter2 format
+ unsigned int num_funs = 2;
+ const struct linop_s* eye_op = linop_identity_create(DIMS, odims);
+ const struct linop_s* ops[2] = { eye_op, eye_op };
+ const struct operator_p_s* prox_ops[2] = { sum_prox, lr_prox };
+ long size = 2 * md_calc_size(DIMS, odims);
+ struct s_data s_data = { size / 2 };
+
+ const struct operator_p_s* sum_xupdate_op = operator_p_create( DIMS, odims, DIMS, odims, (void*) &s_data, sum_xupdate, sum_xupdate_free );
+
+
+ // do recon
+
+ iter2_admm( iconf,
+ NULL,
+ num_funs,
+ prox_ops,
+ ops,
+ sum_xupdate_op,
+ size, (float*) odata, NULL,
+ NULL, NULL, NULL );
+
+
+
+ // Sum
+ if (sum_str)
+ {
+ complex float* sdata = create_cfl(sum_str, DIMS, idims);
+ long istrs[DIMS];
+ long ostrs[DIMS];
+
+ md_calc_strides(DIMS, istrs, idims, sizeof(complex float));
+ md_calc_strides(DIMS, ostrs, odims, sizeof(complex float));
+
+ md_clear(DIMS, idims, sdata, sizeof(complex float));
+ odims[LEVEL_DIM]--;
+ md_zaxpy2(DIMS, odims, istrs, sdata, 1./sqrt(levels), ostrs, odata);
+ odims[LEVEL_DIM]++;
+ unmap_cfl(DIMS, idims, sdata);
+ }
+
+
+ // Clean up
+ unmap_cfl(DIMS, idims, idata);
+ unmap_cfl(DIMS, odims, odata);
+ linop_free( sum_op );
+ operator_p_free( sum_prox );
+ operator_p_free( lr_prox );
+
+
+ double end_time = timestamp();
+ debug_printf(DP_INFO, "Total Time: %f\n", end_time - start_time);
+ exit(0);
+}
+
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/bart.git
More information about the debian-med-commit
mailing list