Merge branch 'comonicon'

README.md

@@ -12,16 +12,48 @@ This project is an improvement of the [Python source code](http://www.neuroimagi
 
 - Langkammer, C., Bredies, K., Poser, B. A., Barth, M., Reishofer, G., Fan, A. P., ... & Ropele, S. (2015). Fast quantitative susceptibility mapping using 3D EPI and total generalized variation. Neuroimage, 111, 622-630.
 
-## Setup
+## Command line usage
+
+### Command line Setup
+
+1. Install [Julia](https://julialang.org/downloads/) (v1.9+ recommended)
+2. Make sure `julia` can be executed from the command line
+3. Download the single file tgv_qsm.jl (<a id="raw-url" href="https://github.com/korbinian90/QuantitativeSusceptibilityMappingTGV.jl/blob/main/tgv_qsm.jl">Download tgv_qsm.jl</a>)
+
+### Run script
+
+```bash
+julia <folder>/tgv_qsm.jl --help
+```
+
+On the first usage, the script will download all dependencies.
+
+### Optional configuration
+
+Under Linux: Make the file executable with `chmod +x tgv_qsm.jl` and run via
+
+```bash
+<folder>/tgv_qsm.jl --help
+```
+
+## Run in Julia
+
+### Setup
 
 1. Install [Julia](https://julialang.org/downloads/) (v1.9+ recommended)
 2. Install this package  
     Open the julia REPL and type:
 
+    ```julia
+    julia> ] # enters package manager
+    pkg> add QuantitativeSusceptibilityMappingTGV MriResearchTools
+    ```
+
+    or
+
     ```julia
     import Pkg
-    Pkg.add(url="https://github.com/korbinian90/QuantitativeSusceptibilityMappingTGV.jl")
-    Pkg.add("MriResearchTools") # for nifti handling
+    Pkg.add(["QuantitativeSusceptibilityMappingTGV", "MriResearchTools"])
     ```
 
 ## Example to run TGV
@@ -61,10 +93,6 @@ This project is an improvement of the [Python source code](http://www.neuroimagi
 
 The first execution might take some time to compile the kernels (~1min).
 
-## Command Line Interface
-
-Coming soon
-
 ## Settings
 
 The default settings were optimized for brain QSM and should lead to good results independently of the acquired resolution.
@@ -91,7 +119,7 @@ It uses the number of CPU threads julia was started with. You can use `julia --t
 
 This implementation doesn't support data with an oblique angle acquisition yet. For rotated data, it is recommended to use the [QSMxT pipeline](https://qsmxt.github.io/QSMxT/) for susceptibility mapping, which applies TGV internally
 
-## Self contained example to test if everything works
+## Self contained example to test if package works
 
 ```julia
 using QuantitativeSusceptibilityMappingTGV
@@ -120,7 +148,7 @@ The parallel CPU version is about twice as fast as the Cython version, the GPU v
 
 ## Run with other GPU types
 
-Other GPU types should be supported, however are only minimally tested.
+Other GPU types don't work with the command line script. They have to be accessed via Julia (or the command line script modified).
 
 ```julia
 using AMDGPU

src/tgv.jl

@@ -1,3 +1,12 @@
+"""
+qsm_tgv(phase, mask, res; TE, kwargs...)
+
+Optional arguments:
+
+Example:
+    chi = qsm_tgv(randn(20,20,20), trues(20,20,20), [1,1,2]; TE=0.025, fieldstrength=3)
+
+"""
 function qsm_tgv(phase, mask, res; TE, B0_dir=[0, 0, 1], fieldstrength=3, regularization=2, alpha=get_default_alpha(regularization), step_size=3, iterations=get_default_iterations(res, step_size), erosions=3, dedimensionalize=false, correct_laplacian=true, laplacian=get_laplace_phase_del, type=Float32, gpu=CUDA.functional(), nblocks=32)
     device, cu = select_device(gpu)
     phase, res, alpha, fieldstrength, mask = adjust_types(type, phase, res, alpha, fieldstrength, mask)
@@ -176,7 +185,7 @@ function select_device(library::Module)
         println("Using the GPU via Metal")
         return library.MetalKernels.MetalBackend(), library.MtlArray
     else
-        println("Using $(Threads.nthreads()) CPU cores")
+        println("Using $(Threads.nthreads()) CPU threads")
         return CPU(), identity
     end
 end
@@ -186,7 +195,7 @@ function select_device(gpu)
         println("Using the GPU")
         return CUDA.CUDAKernels.CUDABackend(), CUDA.cu
     else
-        println("Using $(Threads.nthreads()) CPU cores")
+        println("Using $(Threads.nthreads()) CPU threads")
         return CPU(), identity
     end
 end

tgv_qsm.jl

@@ -0,0 +1,37 @@
+#!/usr/bin/env -S julia --color=yes --startup-file=no --threads=auto
+
+## Usage
+
+# Call with: `<path-to-file>/tgv_qsm.jl ARGS`
+# On windows use: `julia --threads=auto <path-to-file>/tgv_qsm.jl ARGS`
+
+# Example call:
+# `./tgv_qsm.jl phase.nii.gz mask.nii.gz --TE 0.025 --output output.nii.gz --no-gpu
+
+import Pkg
+
+## Uncomment to use a local julia package directory instead of the global one
+# package_dir = joinpath(@__DIR__, ".tgv_cmd_packages")
+# mkpath(package_dir)
+# Pkg.activate(package_dir)
+
+try
+    using QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+catch
+    Pkg.add(["QuantitativeSusceptibilityMappingTGV", "MriResearchTools", "Comonicon"])
+    using QuantitativeSusceptibilityMappingTGV, MriResearchTools, Comonicon
+end
+
+version = Comonicon.get_version(QuantitativeSusceptibilityMappingTGV)
+Comonicon.get_version(::Module) = version
+
+@main function tgv_qsm(fn_phase, fn_mask; TE::Float64, output::String="output.nii.gz", fieldstrength::Float64=3.0, regularization::Float64=2.0, erosions::Int=3, B0_dir::Array{Int}=[0,0,1], dedimensionalize::Bool=false, no_laplacian_correction::Bool=false, step_size::Float64=3.0, no_gpu::Bool=false, type::DataType=Float32, nblocks::Int=32)
+    println("Starting calculation...")
+    phase = readphase(fn_phase)
+    mask = niread(fn_mask) .!= 0
+    res = header(phase).pixdim[2:4]
+    println("Resolution from NIfTI header [mm]: $(round.(Float64.(res); digits=2))")
+    chi = qsm_tgv(phase, mask, res; TE, B0_dir, fieldstrength, regularization, erosions, dedimensionalize, correct_laplacian=!no_laplacian_correction, gpu=!no_gpu, step_size, type, nblocks)
+    println("Writing output")
+    savenii(chi, output; header=header(phase))
+end