Gradient waveform design with projection algorithm: the example of Rosette trajectories

This document shows how to use the algorithm of trajectory projection to design feasible gradient waveforms. It corresponds to Fig. 5 of the paper [Chauffert et al., Gradient Waveform Design for variable density sampling in Magnetic Resonance Imaging]

Contents

close all
clear all
clc
addpath('MRI/')
addpath('functions/')
addpath('operators/')

Enter the Gradient constraints

Parameters of the scanner (here use in [Lustig et al, IEEE TMI 2008])

Gmax = 40e-3;  % T/m
Smax = 150e-3; % T/m/ms
Kmax = 600;     % m^-1

gamma = 42.576*1e3; % kHz/T

alpha = gamma*Gmax;  % in m^-1 ms^-1
beta  = gamma*Smax;  % in m^-1 ms^-2
Dt    = .004;        % sampling time in ms

Choose an input trajectory for the algorithm

Give an input trajectory

w1 = 14.7*2*pi*Gmax;
w2 = 8.7/1.02*2*pi*Gmax;
T = .17/Gmax;
t = 0e-3:Dt:T;
C = Kmax*sin(w1*t').*exp(1i*w2*t');

x=[real(C)';imag(C)'];
s0=parameterize_maximum_speed(x,.9*alpha,Dt)';
sub=1; % subsampling of the curve for visualization
figure, plot(s0(1:sub:end,1),s0(1:sub:end,2),'b.','linewidth',2)
axis equal, axis off
set(gcf,'Color',[1 1 1])
legend('input trajectory')


Gmx = 4;
Smx = 15;
T = 17/Gmx;
Kmx = 6;
w1 = 0.147*2*pi*Gmx;
w2 = 0.087/1.02*2*pi*Gmx;
t = 0e-3:4e-3:T;
C = Kmx*sin(w1*t').*exp(1i*w2*t');
C = [real(C) imag(C) 0*C];
tic
[C_riv, time_riv, g_riv, s_riv, k_riv] = minTimeGradient(C,0);          % Rotationally invariant solution
toc
tic
[C_rv, time_rv, g_rv, s_rv, k_rv]= minTimeGradient(C,1, 0);  % Rotationally variant solution
toc

Rotationally Invariant Solution
Const arc-length parametrization
Compute geometry dependent constraints
Solve ODE forward
Solve ODE backwards
Final Interpolations
Done
Elapsed time is 19.105257 seconds.
Rotationally Variant Solution
Const arc-length parametrization
Compute geometry dependent constraints
Solve ODE forward
Solve ODE backwards
Final Interpolations
Done
Elapsed time is 39.906054 seconds.

Specify constraints

CRV=set_MRI_constraints_RV(alpha,beta,Dt);
CRIV=set_MRI_constraints_RIV(alpha,beta,Dt);

% No additional affine constraint:
C_linear=set_Linear_constraints(size(s0,1),size(s0,2));

Algo_param.nit = 80000;  % number of iterations
Algo_param.L=16;    % Lipschitz constant of the gradient
Algo_param.discretization_step=Dt;

% optional paramteres
Algo_param.show_progression = 0;
Algo_param.display_results = 0;

Project curve with Rotation-Invariant Constraints

tic
s1=Project_Curve_Affine_Constraints(s0,CRIV,C_linear,Algo_param);
toc
% Compute gradients
T=size(s1,1);
t=1:T;
g1=Prime(s1,Dt)/gamma;
gg1=Second(s1,Dt)/gamma;

Elapsed time is 76.880312 seconds.

Project curve with Rotation-Variant Constraints

tic
s2=Project_Curve_Affine_Constraints(s0,CRV,C_linear,Algo_param);
toc

% Compute gradients
g2=Prime(s2,Dt)/gamma;
gg2=Second(s2,Dt)/gamma;

Elapsed time is 55.987054 seconds.

Display projected curves

figure, plot(s0(1:sub:end,1),s0(1:sub:end,2),'--','color',[.5 .5 .5],'linewidth',3)
axis equal, axis off
set(gcf,'Color',[1 1 1])
hold on,
plot(s1(:,1),s1(:,2),'k','linewidth',3)
axis equal, axis off
set(gcf,'Color',[1 1 1])
hold off
title('Projection with RIV constraints')
h=legend('input trajectory', 'projected trajectory','location','sw');
set(h,'FontSize',20,'interpreter','latex');

figure, plot(s0(1:sub:end,1),s0(1:sub:end,2),'--','color',[.5 .5 .5],'linewidth',3)
axis equal, axis off
set(gcf,'Color',[1 1 1])
hold on
plot(s2(:,1),s2(:,2),'k','linewidth',3)
axis equal, axis off
set(gcf,'Color',[1 1 1])
hold off
title('Projection with RV constraints')
h=legend('input trajectory', 'projected trajectory','location','sw');
set(h,'FontSize',20,'interpreter','latex');

Display gradients for rotation variant constraints

T_rv=1.1*time_rv/Dt;
t_rv=1:T_rv;
figure
plot(g_rv(:,1)*1e-2,'r','linewidth',3), axis([0 T_rv -Gmax*1.1 Gmax*1.1])
hold on, plot(g_rv(:,2)*1e-2,'b','linewidth',3), axis([0 T_rv -Gmax*1.1 Gmax*1.1])
set(gcf,'Color',[1 1 1])
hold on, plot(t_rv,0*t_rv, '--k','lineWidth',3)
hold on, plot(t_rv,0*t_rv+Gmax, '--k','lineWidth',2)
hold on, plot(t_rv,0*t_rv-Gmax, '--k','lineWidth',2)
set(gca,'XTick',[0,T_rv/1.1])
set(gca,'XTickLabel',{'0',time_rv})
set(gca,'YTick',[-Gmax,Gmax])
set(gca,'YTickLabel',{-Gmax*1e2,Gmax*1e2})
set(gca,'FontSize',15)
ylabel('G/cm')
xlabel('ms')
hold off
h_legend=legend('$ g_x(t)$','$ g_y(t)$');
set(h_legend,'FontSize',20,'interpreter','latex');
title('reparameterization')

% display gradients
figure,
plot(g2(:,1),'r','linewidth',3), axis([0 T_rv -Gmax*1.1 Gmax*1.1])
hold on, plot(g2(:,2),'b','linewidth',3), axis([0 T_rv -Gmax*1.1 Gmax*1.1])
set(gcf,'Color',[1 1 1])
hold on, plot(t_rv,0*t_rv, '--k','lineWidth',3)
hold on, plot(t_rv,0*t_rv+Gmax, '--k','lineWidth',2)
hold on, plot(t_rv,0*t_rv-Gmax, '--k','lineWidth',2)
set(gca,'XTick',[0,T])
set(gca,'XTickLabel',{'0',T*Dt})
set(gca,'YTick',[-Gmax,Gmax])
set(gca,'YTickLabel',{-Gmax*1e2,Gmax*1e2})
set(gca,'FontSize',15)
ylabel('G/cm')
xlabel('ms')
hold off
h_legend=legend('$ g_x(t)$','$ g_y(t)$');
set(h_legend,'FontSize',20,'interpreter','latex');
title('projection')

Display gradients for rotation invariant constraints

T_riv=1.1*time_riv/Dt;
t_riv=1:T_riv;
figure, plot(g_riv(:,1)*1e-2,'r','linewidth',3), axis([0 T_riv -Gmax*1.1 Gmax*1.1])
hold on, plot(g_riv(:,2)*1e-2,'b','linewidth',3), axis([0 T_riv -Gmax*1.1 Gmax*1.1])
g2n=sqrt(g_riv(:,1).^2+g_riv(:,2).^2)*1e-2;
hold on, plot(g2n,'k','linewidth',3), axis([0 T_riv -Gmax*1.1 Gmax*1.1])
set(gcf,'Color',[1 1 1])
hold on, plot(t_riv,0*t_riv, '--k','lineWidth',3)
hold on, plot(t_riv,0*t_riv+Gmax, '--k','lineWidth',2)
hold on, plot(t_riv,0*t_riv-Gmax, '--k','lineWidth',2)
set(gca,'XTick',[0,T_riv/1.1])
set(gca,'XTickLabel',{'0',time_riv})
set(gca,'YTick',[-Gmax,Gmax])
set(gca,'YTickLabel',{-Gmax*1e2,Gmax*1e2})
set(gca,'FontSize',15)
ylabel('G/cm')
xlabel('ms')
hold off
h_legend=legend('$ g_x(t)$','$ g_y(t)$','$|g(t)|$');
set(h_legend,'FontSize',20,'interpreter','latex');
title('reparameterization')

% display gradients
figure, plot(g1(:,1),'r','linewidth',3), axis([0 T_riv -Gmax*1.1 Gmax*1.1])
hold on, plot(g1(:,2),'b','linewidth',3), axis([0 T_riv -Gmax*1.1 Gmax*1.1])
g1n=sqrt(g1(:,1).^2+g1(:,2).^2);
hold on, plot(g1n,'k','linewidth',3), axis([0 T_riv -Gmax*1.1 Gmax*1.1])
set(gcf,'Color',[1 1 1])
hold on, plot(t_riv,0*t_riv, '--k','lineWidth',3)
hold on, plot(t_riv,0*t_riv+Gmax, '--k','lineWidth',2)
hold on, plot(t_riv,0*t_riv-Gmax, '--k','lineWidth',2)
set(gca,'XTick',[0,T])
set(gca,'XTickLabel',{'0',T*Dt})
set(gca,'YTick',[-Gmax,Gmax])
set(gca,'YTickLabel',{-Gmax*1e2,Gmax*1e2})
set(gca,'FontSize',15)
ylabel('G/cm')
xlabel('ms')
hold off
h_legend=legend('$ g_x(t)$','$ g_y(t)$','$|g(t)|$');
set(h_legend,'FontSize',20,'interpreter','latex');
title('projection')

Published with MATLAB® R2014a