diff --git a/+b0shim/+compute/siemens_basis.m b/+b0shim/+compute/siemens_basis.m new file mode 100644 index 00000000..f908be4c --- /dev/null +++ b/+b0shim/+compute/siemens_basis.m @@ -0,0 +1,175 @@ +function [ basis ] = siemens_basis( ~, X, Y, Z ) +%siemens_basis Return 1st & 2nd order *ideal* Siemens shim fields +% +% SYNTAX +% +% basis = siemens_basis( orders, X, Y, Z ) +% +% DESCRIPTION +% +% The function first wraps `b0shim.compute.spherical_harmonics()` to generate +% 1st and 2nd order spherical harmonic `basis` fields at the grid positions +% given by arrays `X,Y,Z`. *Following Siemens convention*, `basis` is then: +% +% - Reordered along the 4th dimension as *X, Y, Z, Z2, ZX, ZY, X2-Y2, XY* +% +% - Rescaled to Hz/unit-shim, where "unit-shim" refers to the measure displayed +% in the Adjustments card of the Syngo console UI, namely: +% +% - 1 micro-T/m for *X,Y,Z* gradients (= 0.042576 Hz/mm) +% - 1 micro-T/m^2 for 2nd order terms (= 0.000042576 Hz/mm^2) +% +% The returned `basis` is thereby in the form of ideal "shim reference maps", +% ready for optimization. +% +% INPUTS +% +% orders (uint row-vector) +% Degrees of the desired terms in the series expansion, specified as a +% vector of non-negative integers (`[0:1:n]` yields harmonics up to n-th +% order) +% +% X (numeric 2- or 3-d array) +% "Right->Left" grid coordinates in the patient coordinate system +% (i.e. DICOM reference, units of mm) +% +% Y (numeric 2- or 3-d array) +% "Anterior->Posterior" grid coordinates in the patient coordinate system +% (i.e. DICOM reference, units of mm) +% +% Z (numeric 2- or 3-d array) +% "Inferior->Superior" grid coordinates in the patient coordinate system +% (i.e. DICOM reference, units of mm) +% +% OUTPUTS +% +% basis (double 4-D array) +% Spherical harmonic basis fields +% +% NOTES/TODO +% +% 1. For now, `orders` is, in fact, ignored: fixed as [1:2]—which is suitable +% for the Prisma (presumably other Siemens systems as well)—however, the +% 3rd-order shims of the Terra should ultimately be accommodated too. (Requires +% checking the Adjustments/Shim card to see what the corresponding terms and +% values actually are). So, for now, `basis` will always be returned with 8 terms +% along the 4th dim. +% +% See also +% [ b0shim.compute.spherical_harmonics ] + arguments + ~ ;%orders(1,:) {mustBeNumeric,mustBeNonnegative}; + X(:,:,:,1) {mustBeNumeric}; + Y(:,:,:,1) {mustBeNumeric}; + Z(:,:,:,1) {mustBeNumeric}; + end + +sh = b0shim.compute.spherical_harmonics( [1:2], X, Y, Z ) ; + +%% Reorder terms along 4th array dim. in line with Siemens shims: +% X, Y, Z, Z2, ZX, ZY, X2-Y2, XY +sh = reordertosiemens( sh ) ; + +scalingFactors = computenormalizationfactors() ; +basis = zeros( size( sh ) ) ; + +for iCh = 1 : size( sh, 4 ) + basis(:,:,:,iCh) = scalingFactors(iCh) * sh(:,:,:,iCh) ; +end + +end + +% ---------------- +%% Local functions + +function [ sh1 ] = reordertosiemens( sh0 ) +%REORDERTOSIEMENS +% +% sh1 = reordertosiemens( sh0 ) +% +% Reorder 1st-2nd order basis terms along 4th dim. from +% +% 1. Y, Z, X, XY, ZY, Z2, ZX, X2-Y2 (output by b0shim.compute.spherical_harmonics), to +% +% 2. X, Y, Z, Z2, ZX, ZY, X2-Y2, XY (in line with Siemens shims) + +assert( ( nargin == 1 ) && ( size( sh0, 4 ) == 8 ) ) + +sh1(:,:,:,1) = sh0(:,:,:,3) ; +sh1(:,:,:,2) = sh0(:,:,:,1) ; +sh1(:,:,:,3) = sh0(:,:,:,2) ; +sh1(:,:,:,4) = sh0(:,:,:,6) ; +sh1(:,:,:,5) = sh0(:,:,:,7) ; +sh1(:,:,:,6) = sh0(:,:,:,5) ; +sh1(:,:,:,7) = sh0(:,:,:,8) ; +sh1(:,:,:,8) = sh0(:,:,:,4) ; + +end %reordertosiemens() + +function [ scalingFactors ] = computenormalizationfactors() +%COMPUTENORMALIZATIONFACTORS +% +% scalingFactors = computenormalizationfactors() +% +% Returns a vector of `scalingFactors` to apply to the (Siemens-reordered) +% 1st+2nd order spherical harmonic fields for rescaling field terms as +% "shim reference maps" in units of Hz/unit-shim: +% +% Gx, Gy, and Gz should yield 1 micro-T of field shift per metre: +% equivalently, 0.042576 Hz/mm +% +% 2nd order terms should yield 1 micro-T of field shift per metre-squared: +% equivalently, 0.000042576 Hz/mm^2 +% +% Gist: given the stated nominal values, we can pick several arbitrary +% reference positions around the origin/isocenter at which we know what the +% field *should* be, and use that to calculate the appropriate scaling factor. +% +% NOTE: The method has been worked out empirically and has only been tested for +% 2 Siemens Prisma systems. E.g. re: Y, Z, ZX, and XY terms, it was noted that +% their polarity had to be flipped relative to the form given by +% b0shim.compute.spherical_harmonics(). To adapt the code to other systems, +% arbitrary (?) changes along these lines will likely be needed. + +% TODO: For concision/readability, consider defining+solving equations directly +% using MATLAB's symbolic toolbox + +%% create basis on small 3x3x3 mm^3 isotropic grid +[XIso, YIso, ZIso] = meshgrid( [-1:1], [-1:1], [-1:1] ) ; + +sh = b0shim.compute.spherical_harmonics( [1:2], XIso, YIso, ZIso ) ; +sh = reordertosiemens( sh ) ; + +nChannels = size( sh, 4) ; % = 8 +scalingFactors = zeros( nChannels, 1 ) ; + +% indices of reference positions for normalization: +iX1 = find( ( XIso == 1 ) & ( YIso == 0 ) & ( ZIso == 0 ) ) ; +iY1 = find( ( XIso == 0 ) & ( YIso == 1 ) & ( ZIso == 0 ) ) ; +iZ1 = find( ( XIso == 0 ) & ( YIso == 0 ) & ( ZIso == 1 ) ) ; + +iX1Z1 = find( ( XIso == 1 ) & ( YIso == 0 ) & ( ZIso == 1 ) ) ; +iY1Z1 = find( ( XIso == 0 ) & ( YIso == 1 ) & ( ZIso == 1 ) ) ; +iX1Y1 = find( ( XIso == 1 ) & ( YIso == 1 ) & ( ZIso == 0 ) ) ; + +% order the reference indices like the sh field terms +iRef = [iX1 iY1 iZ1 iZ1 iX1Z1 iY1Z1 iX1 iX1Y1]' ; + +% distance from iso/origin to adopted reference point [units: mm] +r = [1 1 1 1 sqrt(2) sqrt(2) 1 sqrt(2)] ; + +%% -------------- +% invert polarity +% Y, Z, ZX, and XY terms only (determined empirically) +sh(:,:,:,[2,3,5,8]) = -sh(:,:,:,[2,3,5,8] ) ; + +%% ------ +% scaling: +orders = [1 1 1 2 2 2 2 2] ; + +for iCh = 1 : nChannels + field = sh(:,:,:,iCh) ; + scalingFactors(iCh) = 42.576*( ( r(iCh) * 0.001 )^orders(iCh) )/field( iRef( iCh ) ) ; +end + +end %computenormalizationfactors() diff --git a/+b0shim/+compute/spherical_harmonics.m b/+b0shim/+compute/spherical_harmonics.m new file mode 100644 index 00000000..cea21b1b --- /dev/null +++ b/+b0shim/+compute/spherical_harmonics.m @@ -0,0 +1,206 @@ +function [ basis ]= spherical_harmonics( orders, X, Y, Z ) +%spherical_harmonics Return orthonormal spherical harmonic basis set +% +% SYNTAX +% +% basis = spherical_harmonics( orders, X, Y, Z ) +% +% DESCRIPTION +% +% Returns an array of spherical harmonic basis fields with the order/degree +% index along the 4th dimension. +% +% INPUTS +% +% orders (uint row-vector) +% Degrees of the desired terms in the series expansion, specified as a +% vector of non-negative integers (`[0:1:n]` yields harmonics up to n-th +% order) +% +% X (numeric 2- or 3-d array) +% 2- or 3-D arrays of grid coordinates +% +% Y (numeric 2- or 3-d array) +% 2- or 3-D arrays of grid coordinates +% +% Z (numeric 2- or 3-d array) +% 2- or 3-D arrays of grid coordinates +% +% X,Y,Z must be identically sized. +% +% OUTPUTS +% +% basis (double 4-D array) +% Spherical harmonic basis fields +% +% EXAMPLE +% +% ``` +% % Initialize grid positions +% [X,Y,Z] = ndgrid([-10:10],[-10:10],[-10:10]); +% +% % 0th-to-2nd order terms inclusive +% orders = [0:2]; +% +% basis = spherical_harmonics(orders, X, Y, Z); +% ``` +% +% - `basis(:,:,:,1)` corresponds to the 0th-order constant term (globally=unity) +% +% - `basis(:,:,:,2:4)` to 1st-order linear terms +% - 2: *y* +% - 3: *z* +% - 4: *x* +% +% - `basis(:,:,:,5:9)` to 2nd-order terms +% - 5: *xy* +% - 6: *zy* +% - 7: *z2* +% - 8: *zx* +% - 9: *x2y2* +% +% NOTES +% +% Based on calc_spherical_harmonics_arb_points_cz.m by jaystock@nmr.mgh.harvard.edu + arguments + orders(1,:) {mustBeNumeric,mustBeNonnegative}; + X(:,:,:,1) {mustBeNumeric}; + Y(:,:,:,1) {mustBeNumeric}; + Z(:,:,:,1) {mustBeNumeric}; + end + +%% Check inputs +assert( isequal(ndims(X), ndims(Y), ndims(Z)), ... + 'Input arrays X, Y, and Z must be identically sized') ; + +switch ndims(X) + case 2 + gridSize = [ size(X) 1 ] ; + case 3 + gridSize = size(X) ; + otherwise + error('Input arrays X, Y, and Z must have 2 or 3 dimensions') ; +end + +%% Initialize variables +nVoxels = numel(X); +nOrders = numel(orders) ; +harm_all = zeros( nVoxels, 1 ) ; + +ii=0; + +%% Compute basis +for iOrder = 1 : nOrders + + n = orders(iOrder); + m = -orders(iOrder):1:orders(iOrder); + + for mm = 1 : numel(m) + + ii = ii+1; + + % The first 2*n(1)+1 columns of the output correspond to harmonics of + % order n(1), and the next 2*n(2)+1 columns correspond to harmonics of + % order n(2), etc. + harm_all(:,ii) = leg_rec_harmonic_cz( n, m(mm), X(:), Y(:), Z(:)); + end + +end + +%% Reshape to initial gridSize +nBasis = size( harm_all, 2) ; +basis = zeros( [gridSize nBasis] ) ; + +for iBasis = 1 : nBasis + basis(:,:,:, iBasis) = reshape( harm_all(:, iBasis), gridSize ) ; +end + +% ---------------- +%% Local functions + +function out = leg_rec_harmonic_cz(n, m, pos_x, pos_y, pos_z) +% returns harmonic field for the required solid harmonic addressed by +% n, m based on the iterative Legendre polynomial calculation. +% +% Positive m values correspond to the cosine component, negative to the sine. +% +% returned fields will eventually follow RRI's convention +% pos_... can be both value and vector/matrix + + r2 = pos_x.^2+pos_y.^2+pos_z.^2; + r = r2.^0.5; + phi = atan2(pos_y, pos_x); + cos_theta = cos(atan2((pos_x.^2+pos_y.^2).^0.5, pos_z)); + %cos_theta=pos_z./r; + + if m>=0, + c=1; + else + c=0; + m=-m; + end + + Ymn = leg_rec(n,m,cos_theta); + + rri_norm = factorial(n+m+1)/factorial(n-m)/ffactorial(2*m); + + out = (n+m+1)*r.^(n).*(cos(m*phi)*c+sin(m*phi)*(1-c)).*Ymn/rri_norm; + + function out = ffactorial(n) + %FFACTORIAL FFactorial (double factorial) function. + + N = n(:); + if any(fix(N) ~= N) || any(N < 0) || ~isa(N,'double') || ~isreal(N) + error('N must be a matrix of non-negative integers.') + end + + if n==0 || n==1 + out=1; + else + out=n*ffactorial(n-2); + end + + end %ffactorial + + function out=leg_rec(n, m, u) + % compute legendre polynomial values for dat using recursive relations + + if m>0 + p_mm=(-1)^m*ffactorial(2*m-1)*(1-u.^2).^(m/2); + else + p_mm=1; + end + + if (n==m) + out=p_mm; + else + p_mm1=(2*m+1)*u.*p_mm; + + if (n==m+1) + out=p_mm1; + else + % recursive calculation needed + a=m+2; + p_ma_2=p_mm; + p_ma_1=p_mm1; + + while 1 + p_ma=((2*a-1)*u.*p_ma_1-(a+m-1)*p_ma_2)/(a-m); + + if a==n + break; + end + % prepare next iteration + p_ma_2=p_ma_1; + p_ma_1=p_ma; + a=a+1; + end + + out=p_ma; + end + end + end %leg_rec + +end %leg_rec_harmonic_cz + +end % main diff --git a/Coils/Shim_Sim/ShimOpt_SphericalHarmonics/ShimOpt_SphericalHarmonics.m b/Coils/Shim_Sim/ShimOpt_SphericalHarmonics/ShimOpt_SphericalHarmonics.m index 9a637aca..a3a3325b 100644 --- a/Coils/Shim_Sim/ShimOpt_SphericalHarmonics/ShimOpt_SphericalHarmonics.m +++ b/Coils/Shim_Sim/ShimOpt_SphericalHarmonics/ShimOpt_SphericalHarmonics.m @@ -171,294 +171,14 @@ methods(Static) % ========================================================================= function [ basisFields ]= generatebasisfields( orders, X, Y, Z ) -%GENERATEBASISFIELDS -% -% Generates orthonormal spherical harmonic (SH) basis set -% -% ....... -% -% Usage -% -% [ basisFields ] = GENERATEBASISFIELDS( orders, X, Y, Z ) -% -% Returns array of SH basis fields where the 4th dimension is the order/degree index -% e.g. if orders = [0:2], -% -% basisFields(:,:,:,1) corresponds to the 0th order term, -% -% basisFields(:,:,:,2:4) to 1st order terms -% 2 -> (y) -% 3 -> (z) -% 4 -> (x) -% -% basisFields(:,:,:,5:9) to 2nd orders -% 5 -> (xy) -% 6 -> (zy) -% 7 -> (z2) -% 8 -> (zx) -% 9 -> (x2y2) -% -% etc. -% -% Input -% -% orders -% vector of non-negative integers orders to calculate (e.g. [0:3]). -% Typically orders is specified as 0:1:N to obtain spherical harmonics -% up to order N -% -% X, Y, Z -% 3d arrays specifying voxel coordinates at which to calculate the harmonics -% -% ....... -% -% Based on calc_spherical_harmonics_arb_points_cz.m by jaystock@nmr.mgh.harvard.edu - - -% set the orders_to_calculate as a vector of integers n with 0 < n < N. This -% calculates accompanying spherical harmonic orders. The first 2*n(1)+1 columns -% of the output correspond to harmonics of order n(1), and the next -% 2*n(2)+1 columns correspond to harmonics of order n(2), etc. % -% note that Hetherington refers to n as "degree" instead of "order" -% - -assert( all( orders >=0 ) ) ; - -switch ndims(X) - case 2 - gridSize = [ size(X) 1 ] ; - case 3 - gridSize = size(X) ; - otherwise - error('Input arrays X,Y, and Z must have 2 or 3 dimensions') ; -end - -nVoxels = numel(X); -nOrders = numel(orders) ; - -harm_all = zeros( nVoxels, 1 ) ; - -ii=0; - -for iOrder = 1 : nOrders - - n = orders(iOrder); - m = -orders(iOrder):1:orders(iOrder); - - for mm = 1 : numel(m) - - ii = ii+1; - - harm_all(:,ii) = leg_rec_harmonic_cz( n, m(mm), X(:), Y(:), Z(:)); - end -end - -nBasisFields = size( harm_all, 2) ; -basisFields = zeros( [gridSize nBasisFields] ) ; - -for iBasisField = 1 : nBasisFields - basisFields(:,:,:, iBasisField) = reshape( harm_all(:, iBasisField), gridSize ) ; -end - -function out = leg_rec_harmonic_cz(n, m, pos_x, pos_y, pos_z) -% returns harmonic field for the required solid harmonic addressed by -% n, m based on the iterative Legendre polynomial calculation -% Positive m values correspond to cosine component and negative to sine -% -% returned fields will eventually follow RRI's convention -% pos_... can be both value and vector/matrix - - r2=pos_x.^2+pos_y.^2+pos_z.^2; - r=r2.^0.5; - phi=atan2(pos_y, pos_x); - cos_theta=cos(atan2((pos_x.^2+pos_y.^2).^0.5, pos_z)); - %cos_theta=pos_z./r; - - if m>=0, - c=1; - else - c=0; - m=-m; - end - - Ymn=leg_rec(n,m,cos_theta); - - rri_norm=factorial(n+m+1)/factorial(n-m)/ffactorial(2*m); - - out=(n+m+1)*r.^(n).*(cos(m*phi)*c+sin(m*phi)*(1-c)).*Ymn/rri_norm; - -function out = ffactorial(n) -%FFACTORIAL FFactorial (double factorial) function. - - N = n(:); - if any(fix(N) ~= N) || any(N < 0) || ~isa(N,'double') || ~isreal(N) - error('MATLAB:factorial:NNegativeInt', ... - 'N must be a matrix of non-negative integers.') - end - - if n==0 || n==1 - out=1; - else - out=n*ffactorial(n-2); - end -end - -function out=leg_rec(n, m, u) -% compute legendre polynomial values for dat using recursive relations - - if m>0 - p_mm=(-1)^m*ffactorial(2*m-1)*(1-u.^2).^(m/2); - else - p_mm=1; - end - - if (n==m) - out=p_mm; - else - p_mm1=(2*m+1)*u.*p_mm; - - if (n==m+1) - out=p_mm1; - else - % recursive calculation needed - a=m+2; - p_ma_2=p_mm; - p_ma_1=p_mm1; - - while 1 - p_ma=((2*a-1)*u.*p_ma_1-(a+m-1)*p_ma_2)/(a-m); - - if a==n - break; - end - % prepare next iteration - p_ma_2=p_ma_1; - p_ma_1=p_ma; - a=a+1; - end - - out=p_ma; - end - end -end - -end +basisFields = b0shim.compute.spherical_harmonics( orders, X, Y, Z); end % ========================================================================= function [ basisFields ] = generatebasisfields_siemens( X, Y, Z ) -%GENERATEBASISFIELDS_SIEMENS -% -% Wraps to ShimOpt_SphericalHarmonics.generatebasisfields(), reorders, and -% rescales the basis set to return ideal "shim reference maps" (in units of -% Hz/unit-shim) for the 1st and 2nd order spherical harmonic shims of the -% Siemens Prisma -% -% Usage -% -% [ basisFields ] = GENERATEBASISFIELDS_SIEMENS( X, Y, Z ) -% -% X, Y, Z are 3d arrays of the voxel positions [units: mm] at which to generate -% the basis functions (e.g. X,Y,Z are generally returned from MaRdI.getvoxelpositions() ) - -assert( nargin == 3 ) - -sh = ShimOpt_SphericalHarmonics.generatebasisfields( [1:2], X, Y, Z ) ; -% Reorder terms along 4th array dim. in line with Siemens shims: X, Y, Z, Z2, ZX, ZY, X2-Y2, XY -sh = reordertosiemens( sh ) ; - -scalingFactors = computenormalizationfactors() ; -basisFields = zeros( size( sh ) ) ; - -for iCh = 1 : size( sh, 4 ) - basisFields(:,:,:,iCh) = scalingFactors(iCh) * sh(:,:,:,iCh) ; -end - -return; - -function [ sh1 ] = reordertosiemens( sh0 ) -%REORDERTOSIEMENS -% -% basisFields1 = REORDERTOSIEMENS( basisFields0 ) -% -% basisFields returned from .GENERATEBASISFIELDS() are ordered (along the 4th -% array dimension) like: -% Y, Z, X, XY, ZY, Z2, ZX, X2-Y2 -% -% REORDERTOSIEMENS orders them in line with Siemens shims: -% X, Y, Z, Z2, ZX, ZY, X2-Y2, XY - -assert( ( nargin == 1 ) && ( size( sh0, 4 ) == 8 ) ) - -sh1(:,:,:,1) = sh0(:,:,:,3) ; -sh1(:,:,:,2) = sh0(:,:,:,1) ; -sh1(:,:,:,3) = sh0(:,:,:,2) ; -sh1(:,:,:,4) = sh0(:,:,:,6) ; -sh1(:,:,:,5) = sh0(:,:,:,7) ; -sh1(:,:,:,6) = sh0(:,:,:,5) ; -sh1(:,:,:,7) = sh0(:,:,:,8) ; -sh1(:,:,:,8) = sh0(:,:,:,4) ; - -end %reordertosiemens() - -function [ scalingFactors ] = computenormalizationfactors() -%COMPUTENORMALIZATIONFACTORS -% -% scalingFactors = computenormalizationfactors() -% -% returns a vector of scalingFactors to apply to the (properly reordered) -% ideal 1st+2nd order spherical harmonic basis returned from GENERATEBASISFIELD -% to scale the terms as "shim reference maps" in units of Hz/unit-shim -% ----- -% Gx, Gy, and Gz should yield 1 micro-T of field shift per metre -% equivalently, 0.042576 Hz/mm -% -% 2nd order terms should yield 1 micro-T of field shift per metre-squared -% equivalently, 0.000042576 Hz/mm^2 - -%% ------ -% create basis on small 3x3x3 mm^3 isotropic grid -[XIso, YIso, ZIso] = meshgrid( [-1:1], [-1:1], [-1:1] ) ; - -sh = ShimOpt_SphericalHarmonics.generatebasisfields( [1:2], XIso, YIso, ZIso ) ; -% Reorder terms along 4th array dim. in line with Siemens shims: X, Y, Z, Z2, ZX, ZY, X2-Y2, XY -sh = reordertosiemens( sh ) ; - -nChannels = size( sh, 4) ; % = 8 -scalingFactors = zeros( nChannels, 1 ) ; - -% indices of reference positions for normalization: -iX1 = find( ( XIso == 1 ) & ( YIso == 0 ) & ( ZIso == 0 ) ) ; -iY1 = find( ( XIso == 0 ) & ( YIso == 1 ) & ( ZIso == 0 ) ) ; -iZ1 = find( ( XIso == 0 ) & ( YIso == 0 ) & ( ZIso == 1 ) ) ; - -iX1Z1 = find( ( XIso == 1 ) & ( YIso == 0 ) & ( ZIso == 1 ) ) ; -iY1Z1 = find( ( XIso == 0 ) & ( YIso == 1 ) & ( ZIso == 1 ) ) ; -iX1Y1 = find( ( XIso == 1 ) & ( YIso == 1 ) & ( ZIso == 0 ) ) ; - -% order the reference indices like the sh field terms -iRef = [iX1 iY1 iZ1 iZ1 iX1Z1 iY1Z1 iX1 iX1Y1]' ; - -%% ------ -% scaling: -% 1st order terms yield 1 micro-T of field shift per m (i.e 0.042576 Hz/mm ) -% 2nd order terms yield 1 micro-T of field shift per m^2 (i.e 0.000042576 Hz/mm^2 ) - -% distance from iso/origin to adopted reference point [units: mm] -r = [1 1 1 1 sqrt(2) sqrt(2) 1 sqrt(2)] ; - -% invert polarity of certain terms: -sh(:,:,:,[2,3,5,8]) = -sh(:,:,:,[2,3,5,8] ) ; - -orders = [1 1 1 2 2 2 2 2] ; - -for iCh = 1 : nChannels - field = sh(:,:,:,iCh) ; - scalingFactors(iCh) = 42.576*( ( r(iCh) * 0.001 )^orders(iCh) )/field( iRef( iCh ) ) ; -end -end %computenormalizationfactors() +basisFields = b0shim.compute.siemens_basis([1:2],X,Y,Z); end % =========================================================================