Combined programs

CMakeLists.txt

@@ -0,0 +1,41 @@
+cmake_minimum_required( VERSION 2.6 FATAL_ERROR )
+project( toolkid )
+
+find_package( ITK REQUIRED )
+include( ${ITK_USE_FILE} )
+#include( UseITK.cmake )
+
+# cern ROOT requirements
+set( USE_ROOT ON CACHE BOOL "Use ROOT support" )
+if( USE_ROOT )
+    set( ROOT_HOME "" CACHE FILEPATH "ROOT home" )
+    include( FindROOT.cmake )
+endif( USE_ROOT )
+
+# GSL requirements
+set( USE_GSL ON CACHE BOOL "Use GSL support" )
+if( USE_GSL )
+    set( GSL_HOME "" CACHE FILEPATH "GSL home" )
+    include( FindGSL.cmake )
+endif( USE_GSL )
+
+# Boost
+set( Boost_USE_STATIC_LIBS ON )
+set( Boost_USE_MULTITHREAD ON )
+
+find_package( Boost COMPONENTS program_options system graph regex REQUIRED )
+
+message( "Boost program_options library: " ${Boost_PROGRAM_OPTIONS_LIBRARY} )
+message( "Boost system library: " ${Boost_SYSTEM_LIBRARY} )
+message( "Boost graph library: " ${Boost_GRAPH_LIBRARY} )
+message( "Boost regex library: " ${Boost_REGEX_LIBRARY} )
+
+message( "Boost include directories: " ${Boost_INCLUDE_DIRS} )
+
+set( USE_BLAS ON CACHE BOOL "Build with BLAS support" )
+mark_as_advanced( USE_BLAS )
+
+set( LIBRARY_OUTPUT_PATH ${toolkid_BINARY_DIR}/lib CACHE INTERNAL "Single library output directory" )
+set( EXECUTABLE_OUTPUT_PATH ${toolkid_BINARY_DIR}/bin CACHE INTERNAL "Single binary output directory" )
+
+subdirs( src )

FindGSL.cmake

@@ -0,0 +1,136 @@
+# Try to find gnu scientific library GSL
+# See
+# http://www.gnu.org/software/gsl/  and
+# http://gnuwin32.sourceforge.net/packages/gsl.htm
+#
+# Based on a script of Felix Woelk and Jan Woetzel
+# (www.mip.informatik.uni-kiel.de)
+#
+# It defines the following variables:
+#  GSL_FOUND - system has GSL lib
+#  GSL_INCLUDE_DIRS - where to find headers
+#  GSL_LIBRARIES - full path to the libraries
+#  GSL_LIBRARY_DIRS, the directory where the PLplot library is found.
+
+#  CMAKE_GSL_CXX_FLAGS  = Unix compiler flags for GSL, essentially "`gsl-config --cxxflags`"
+#  GSL_LINK_DIRECTORIES = link directories, useful for rpath on Unix
+#  GSL_EXE_LINKER_FLAGS = rpath on Unix
+
+set( GSL_FOUND OFF )
+set( GSL_CBLAS_FOUND OFF )
+
+# Windows, but not for Cygwin and MSys where gsl-config is available
+if( WIN32 AND NOT CYGWIN AND NOT MSYS )
+# look for headers
+    find_path( GSL_INCLUDE_DIR
+            NAMES gsl/gsl_cdf.h gsl/gsl_randist.h
+            )
+if( GSL_INCLUDE_DIR )
+# look for gsl library
+    find_library( GSL_LIBRARY
+            NAMES gsl
+            )
+    if( GSL_LIBRARY )
+    set( GSL_INCLUDE_DIRS ${GSL_INCLUDE_DIR} )
+    get_filename_component( GSL_LIBRARY_DIRS ${GSL_LIBRARY} PATH )
+    set( GSL_FOUND ON )
+endif( GSL_LIBRARY )
+
+# look for gsl cblas library
+    find_library( GSL_CBLAS_LIBRARY
+            NAMES gslcblas
+            )
+    if( GSL_CBLAS_LIBRARY )
+    set( GSL_CBLAS_FOUND ON )
+endif( GSL_CBLAS_LIBRARY )
+
+    set( GSL_LIBRARIES ${GSL_LIBRARY} ${GSL_CBLAS_LIBRARY} )
+endif( GSL_INCLUDE_DIR )
+
+    mark_as_advanced(
+            GSL_INCLUDE_DIR
+            GSL_LIBRARY
+            GSL_CBLAS_LIBRARY
+            )
+    else( WIN32 AND NOT CYGWIN AND NOT MSYS )
+if( UNIX OR MSYS )
+    find_program( GSL_CONFIG_EXECUTABLE gsl-config
+            /usr/bin/
+            /usr/local/bin
+            )
+
+    if( GSL_CONFIG_EXECUTABLE )
+set( GSL_FOUND ON )
+
+# run the gsl-config program to get cxxflags
+    execute_process(
+            COMMAND sh "${GSL_CONFIG_EXECUTABLE}" --cflags
+            OUTPUT_VARIABLE GSL_CFLAGS
+            RESULT_VARIABLE RET
+            ERROR_QUIET
+            )
+if( RET EQUAL 0 )
+    string( STRIP "${GSL_CFLAGS}" GSL_CFLAGS )
+separate_arguments( GSL_CFLAGS )
+
+# parse definitions from cflags; drop -D* from CFLAGS
+    string( REGEX MATCHALL "-D[^;]+"
+            GSL_DEFINITIONS  "${GSL_CFLAGS}" )
+    string( REGEX REPLACE "-D[^;]+;" ""
+            GSL_CFLAGS "${GSL_CFLAGS}" )
+
+# parse include dirs from cflags; drop -I prefix
+    string( REGEX MATCHALL "-I[^;]+"
+            GSL_INCLUDE_DIRS "${GSL_CFLAGS}" )
+    string( REPLACE "-I" ""
+            GSL_INCLUDE_DIRS "${GSL_INCLUDE_DIRS}")
+    string( REGEX REPLACE "-I[^;]+;" ""
+            GSL_CFLAGS "${GSL_CFLAGS}")
+
+    message("GSL_DEFINITIONS=${GSL_DEFINITIONS}")
+    message("GSL_INCLUDE_DIRS=${GSL_INCLUDE_DIRS}")
+    message("GSL_CFLAGS=${GSL_CFLAGS}")
+    else( RET EQUAL 0 )
+    set( GSL_FOUND FALSE )
+endif( RET EQUAL 0 )
+
+# run the gsl-config program to get the libs
+    execute_process(
+            COMMAND sh "${GSL_CONFIG_EXECUTABLE}" --libs
+            OUTPUT_VARIABLE GSL_LIBRARIES
+            RESULT_VARIABLE RET
+            ERROR_QUIET
+            )
+if( RET EQUAL 0 )
+    string(STRIP "${GSL_LIBRARIES}" GSL_LIBRARIES )
+separate_arguments( GSL_LIBRARIES )
+
+# extract linkdirs (-L) for rpath (i.e., LINK_DIRECTORIES)
+    string( REGEX MATCHALL "-L[^;]+"
+            GSL_LIBRARY_DIRS "${GSL_LIBRARIES}" )
+    string( REPLACE "-L" ""
+            GSL_LIBRARY_DIRS "${GSL_LIBRARY_DIRS}" )
+    else( RET EQUAL 0 )
+    set( GSL_FOUND FALSE )
+endif( RET EQUAL 0 )
+
+    MARK_AS_ADVANCED(
+            GSL_CFLAGS
+            )
+    message( STATUS "Using GSL from ${GSL_PREFIX}" )
+else( GSL_CONFIG_EXECUTABLE )
+    message( STATUS "FindGSL: gsl-config not found.")
+    endif( GSL_CONFIG_EXECUTABLE )
+    endif( UNIX OR MSYS )
+endif( WIN32 AND NOT CYGWIN AND NOT MSYS )
+
+    if( GSL_FOUND )
+if( NOT GSL_FIND_QUIETLY )
+    message( STATUS "FindGSL: Found both GSL headers and library" )
+    endif( NOT GSL_FIND_QUIETLY )
+    else( GSL_FOUND )
+if( GSL_FIND_REQUIRED )
+    message( FATAL_ERROR "FindGSL: Could not find GSL headers or library" )
+    endif( GSL_FIND_REQUIRED )
+endif( GSL_FOUND )
+

matlab/backup.looklocker.m

@@ -0,0 +1,36 @@
+% Name : procsems_ll
+% Optimized by Wim (26-05-2009)
+%
+% RUN as: looklocker( 'gems_01', 'processed/gems', 1, 1 )
+% -------------------------------------------------------
+function [im] = looklocker( filename, exportname, zffactor1, zffactor2 );
+
+fid                 = read_fid( filename );                   % size( fid ) -> 64 38400
+ne                  = read_procpar( filename, 'ne' );         % 24
+nv                  = read_procpar( filename, 'nv' );         % 64      phase-encoding steps
+np                  = read_procpar( filename, 'np' ) * 0.5;   % 64      readout steps: is always RO * 2 ! (so devide by 2).
+ns                  = read_procpar( filename, 'ns' );         % 25      slices
+pss                 = read_procpar( filename, 'pss' );        % -0.15 -0.05 0.05 0.15 0.25   [...] 1    (25 values...)
+[pssord, sortindex] = sort( pss );                            % sortindex ->  1 14 2 15 3 16 [...] 13   (25 values...)
+
+% reshape...
+fid = reshape( fid, nv, ns, ne, np );         % -> 64    25    24    64
+fid = permute( fid, [1,4,2,3]      );         % -> 64    64    25    24
+fid = reshape( fid, nv, np, ns, ne );         % -> 64    64    25    24
+
+% zerofill...
+zf1 = pow2( ceil( log2( nv ) ) ) * zffactor1;
+zf2 = pow2( ceil( log2( np ) ) ) * zffactor2; 
+
+im  = zeros( zf1, zf2, ns, ne );
+
+% fft...
+for slicecount=1:ns;       % 1:25
+    for arraycount=1:ne;   % 1:24
+      im( :, :, slicecount, arraycount ) = fft2d( fid ( :, :, sortindex( slicecount ), arraycount ), zf1, zf2 );
+    end;
+end;
+
+% save...
+exportbfloat( abs( im ), exportname );
+

matlab/epiphasemap.m

@@ -0,0 +1,45 @@
+% phasemap = epiphasemap( trace, zeroFillRO )
+%   Build a phase map for EPI phase correction
+%
+%   trace : FID trace (RO x PE matrix)
+%   zeroFillRO : output vector length in RO-direction
+%
+function phasemap = epiphasemap( trace, zeroFillRO ),
+
+numberOfRO = size( trace, 1 );
+numberOfPE = size( trace, 2 );
+
+% transform in RO direction
+transform = fft1d( trace, zeroFillRO );
+
+% mask background
+mask = mkmask( transform, 0.95 );
+
+% calculate fitting weights
+weight = abs( transform ).^2;
+weight = weight / max( weight );
+
+% determine full x-axis
+x = [ 1 : zeroFillRO ]';
+
+% initialize phasemap
+phasemap = zeros( zeroFillRO, numberOfPE );
+
+for i = 1 : numberOfPE,
+  % calculate current phase line's phase
+  phase = phasecalc( transform( :, i ) );
+
+  % fit a second-order polynomial a*x^2 + b*x + c
+  maskX = find( mask( :, i ) > 0 );
+  [ a, b, c ] = fitsvd( x( maskX ), phase( maskX ), weight( maskX ) );
+  phasemap( :, i ) = mkpoly( a, b, c, x );
+
+  % correct with mean phase
+  meanphase = phasecalc( phasewithmap( transform( :, i ), phasemap( :, i ) ) );
+
+  % store in phasemap
+  phasemap( :, i ) = phasemap( :, i ) - meanphase( round( zeroFillRO / 2 ) );
+end
+
+return
+

matlab/episortro.m

@@ -0,0 +1,35 @@
+% episortro( fid, filename )
+%   Re-order traces and points in an EPI FID
+%
+%   fid : FID (np x ntraces x nblocks)
+%   filename : path to the FID (for procpar)
+%
+function [ fid, np, nv ] = episortro( fid, filename ),
+
+% get header
+header = fid.header;
+dim = fid.dim;
+fid = fid.fid;
+
+% load RO-pattern table
+patternRO = procpar( filename, 'ropat' );
+load( [ '/home1/maurits/vnmrsys/epi_indices/' patternRO '.mat' ], 'EPI_indices', 'nv', 'np' );
+
+numberOfRO = np; % number of points per trace
+numberOfPE = nv; % number of phase encoding steps
+numberOfBlocks = header.nblocks; % number of blocks
+numberOfTraces = header.ntraces; % number of traces
+
+% order FID
+orderFid = zeros( numberOfRO, numberOfPE, numberOfTraces, numberOfBlocks );
+
+for i = 1 : numberOfRO,
+  for j = 1 : numberOfPE,
+    orderFid( i, j, :, : ) = fid( EPI_indices( i, j ), :, : );
+  end
+end
+
+fid = struct( 'header', header, 'dim', dim, 'fid', orderFid );
+
+return
+

matlab/fft1d.m

@@ -0,0 +1,36 @@
+% transform = fft1d( fid, N )
+%   Fourier transform (ifft) in dimension 1
+%
+%   fid : input spectrum (1D, 2D, or 3D)
+%   N   : length of the spectrum in dimension 1
+%         N > size( fid, 1 ): zero-filling the spectrum
+%         N < size( fid, 1 ): center of the spectrum is retained
+%
+function transform = fft1d( fid, N );
+
+inputSize = size( fid );
+outputSize = size( fid );
+
+% output size in the first dimension
+outputSize( 1 ) = N;
+
+% prepare output
+transform = zeros( outputSize );
+
+if N < inputSize( 1 ),
+  % keep center of k-space
+  length = inputSize( 1 ) - N;
+  startIndex = round( length / 2 );
+  endIndex = startIndex + N - 1;
+
+  transform( :, :, : ) = fid( startIndex : endIndex, :, : );
+else
+  % zero-filling / copying
+  transform( 1 : inputSize( 1 ), :, : ) = fid;
+end
+
+% shift, transform, center
+transform = circshift( ifft( circshift( transform, round( -1 * inputSize( 1 ) / 2 ) ) ), N / 2 );
+
+return
+

matlab/fft2d.m

@@ -0,0 +1,9 @@
+% transform = fft2d( fid, N, M )
+%   2-D Fourier transform
+%
+function transform = fft2d( fid, N, M ),
+
+transform = fft1d( fft1d( fid, N )', M )';
+
+return
+

matlab/fft3d.m

@@ -0,0 +1,21 @@
+% transform = fft3d( fid, N, M, P )
+%   3-D Fourier transform
+%
+function transform = fft3d( fid, N, M, P ),
+
+inputSize = size( fid );
+outputSize = [ N M P ];
+
+% prepare output
+transform = zeros( outputSize );
+
+% for each PE2 line, 2-D transform
+for i = 1 : inputSize( 3 ),
+  transform( :, :, i ) = fft2d( fid( :, :, i ), N, M );
+end
+
+% permute, transform in PE2-direction, permute back
+transform = permute( fft1d( permute( transform, [ 3 1 2 ] ), P ), [ 2 3 1 ] );
+
+return
+