Skip to content
Snippets Groups Projects
Commit 05a1fb73 authored by Maria Engel's avatar Maria Engel
Browse files

Consequently rename Narsude methof to blipped-CAIPIRINHA B, and the other one A, as in abstract;

make sure to include optimum blipped CAIPIRINHA patterns for high L, even if they cannot be globally optimal anymore;
parent 07c48630
No related branches found
No related tags found
1 merge request!4Consequently rename Narsude methof to blipped-CAIPIRINHA B, and the other one A, as in abstract;
function [packingDensity, L_optimum, isDegenerate, packingDensityCAIPI, L_CAIPI, fh] = ...
function [packingDensity, L_optimum, isDegenerate, L_CAIPI_A, packingDensityCAIPI_A, ...
L_CAIPI_B, packingDensityCAIPI_B, fh] = ...
ComputeMostHexagonalSampling(R_inplane, FOV_inplane, sliceSpacing, nSlicesSimultaneously, ...
nSlicesTotal, sliceThickness, sliceGap, doPlot, doDebug, gifName, Lincrement)
% This function computes how to optimally sample the phase encoding plane for 3D/Multiband
......@@ -46,6 +47,16 @@ function [packingDensity, L_optimum, isDegenerate, packingDensityCAIPI, L_CAIPI,
%
% isDegenerate flag whether packing density metric is degenerate
%
% L_CAIPI_A options of L offered by blipped-CAIPRINHA using the formalism described
% in Zahneisen et al. MRM 2014 ("A" in ISMRM abstract Engel et al. 2024)
%
% packingDensityCAIPI_A packing density belonging to L_CAIPI_A
%
% L_CAIPI_B options of L offered by blipped-CAIPRINHA using the formalism described
% in Narsude et al. MRM 2016 ("B" in ISMRM abstract Engel et al. 2024)
%
% packingDensityCAIPI_B packing density belonging to L_CAIPI_B
%
% fh cell figure handles to grid figure (fh{1}) and packing density
% figure (fh{2})
%
......@@ -81,8 +92,8 @@ kzSlabThickness = 2*pi/sliceSpacing;
% no point in checking their packing density
Lmax = ceil(sqrt(0.5+sqrt(0.25+(kzSlabThickness/deltak_inplane)^2)));
% However for demonstration purposes, we still want to take into account all blipped-CAIPIRINHA
% patterns for plotting
if (doPlot || doDebug) && ~isempty(nSlicesSimultaneously)
% patterns for plotting or when we are interested in the optimum blipped-CAIPIRINHA pattern
if (doPlot || doDebug || nargout>3) && ~isempty(nSlicesSimultaneously)
Lmax = max(Lmax,nSlicesSimultaneously);
end
% For L<2 there is an L>2 which produces a rotationally symmetric grid i.e. same packing density
......@@ -99,10 +110,17 @@ LArray = Lmin:Lincrement:Lmax;
DELTA = 1:nSlicesSimultaneously;
% Translate these DELTA's to L's and find which ones present unique patterns (i.e. not yet covered
% by the previous L's)
LNarsudeUnique = nSlicesSimultaneously./DELTA;
LNarsudeUnique = LNarsudeUnique(mod(LNarsudeUnique,1)~=0);
LNarsudeUnique = LNarsudeUnique(LNarsudeUnique>2);
LArray = unique(sort([LArray LNarsudeUnique]));
L_CAIPI_B = nSlicesSimultaneously./DELTA;
L_CAIPI_B = L_CAIPI_B(L_CAIPI_B>=2);
L_CAIPI_B_Unique = L_CAIPI_B(mod(L_CAIPI_B,1)~=0);
LArray = unique(sort([LArray L_CAIPI_B_Unique]));
isBlippedCAIPI_A = ((ceil(LArray)==LArray) + (LArray<=nSlicesSimultaneously)) > 1; % This is being kind to blipped CAIPIRINHA!
isBlippedCAIPI_B = ismember(LArray, L_CAIPI_B);
colourCAIPI_A = [217, 83, 25]/255; % red
colourCAIPI_B = [148, 0, 211]/255; % purple
colourOptimum = [34, 139, 34]/255; % green
howFar2Look = 10; % TODO: Find out what to put here!
d_min = zeros(size(LArray));
......@@ -112,10 +130,8 @@ if doDebug
ah = subplot(1,2,1);
end
%% Initiate arrays
%% Initiate array
packingDensityArray = zeros(numel(LArray),1);
isBlippedCAIPI = zeros(numel(LArray),1);
isNarsude = zeros(numel(LArray),1);
%% Go through range of L's, compute the minimum distance between grid points and the resulting packing density
for iL = 1:numel(LArray)
......@@ -126,9 +142,6 @@ for iL = 1:numel(LArray)
packingDensityArray(iL) = ComputePackingDensity(d_min(iL), deltak_inplane, sliceSpacing);
isBlippedCAIPI(iL) = ceil(LArray(iL))==LArray(iL);% && LArray(iL)<=nSlicesSimultaneously; % This is being kind to blipped CAIPIRINHA!
isNarsude(iL) = ismember(LArray(iL), LNarsudeUnique);
if doDebug
[kyArray, kzArray] = GetGridPoints(3*Lmax, deltak_inplane, sliceSpacing, LArray(iL));
PlotGrid(kyArray, kzArray, sliceSpacing, ah);
......@@ -136,8 +149,10 @@ for iL = 1:numel(LArray)
subplot(1,2,2); hold on;
plot(LArray(iL),packingDensityArray(iL),'.','Color','k');
if isBlippedCAIPI(iL) || isNarsude(iL)
plot(LArray(iL),packingDensityArray(iL),'x','Color','r','MarkerSize',10,'LineWidth',3);
if isBlippedCAIPI_A(iL)
plot(LArray(iL),packingDensityArray(iL),'x','Color',colourCAIPI_A,'MarkerSize',10,'LineWidth',3);
elseif isBlippedCAIPI_B(iL)
plot(LArray(iL),packingDensityArray(iL),'x','Color',colourCAIPI_B,'MarkerSize',10,'LineWidth',3);
end
xlabel('L'); ylabel('Packing density [%]');
......@@ -169,41 +184,41 @@ isDegenerate = d_optimum > kzSlabThickness;
packingDensity = ComputePackingDensity(d_optimum, deltak_inplane, sliceSpacing);
%% Compute options that blipped-CAIPIRINHA would offer
indCAIPI = find(isBlippedCAIPI + isNarsude);
d_CAIPI = d_min(indCAIPI);
L_CAIPI = LArray(indCAIPI);
packingDensityCAIPI = ComputePackingDensity(d_CAIPI,deltak_inplane,sliceSpacing);
indNarsude = find(isNarsude);
d_Narsude = d_min(indNarsude);
L_Narsude = LArray(indNarsude);
packingDensityNarsude = ComputePackingDensity(d_Narsude,deltak_inplane,sliceSpacing);
indCAIPI_A = find(isBlippedCAIPI_A);
d_CAIPI_A = d_min(indCAIPI_A);
L_CAIPI_A = LArray(indCAIPI_A);
packingDensityCAIPI_A = ComputePackingDensity(d_CAIPI_A,deltak_inplane,sliceSpacing);
indCAIPI_B = find(isBlippedCAIPI_B);
d_CAIPI_B = d_min(indCAIPI_B);
L_CAIPI_B = LArray(indCAIPI_B);
packingDensityCAIPI_B = ComputePackingDensity(d_CAIPI_B,deltak_inplane,sliceSpacing);
%% Plot grid and resulting sampling density for all options
if doPlot
fh{1} = figure('Name','Sampling grid');
[Loptions, indOptions] = sort([L_optimum L_CAIPI]);
PackingDensityoptions = round([packingDensity packingDensityCAIPI]);
[Loptions, indOptions] = (unique([L_optimum L_CAIPI_A L_CAIPI_B]));
PackingDensityoptions = round([packingDensity packingDensityCAIPI_A packingDensityCAIPI_B]*10)/10;
for iL = 1:numel(Loptions)
bh = subplot(1,numel(L_CAIPI)+1,iL);
bh = subplot(1,numel(Loptions),iL);
% For plotting get grid points on a bit larger range
[kyArray, kzArray] = GetGridPoints(3*ceil(Lmax), deltak_inplane, sliceSpacing, Loptions(iL));
PlotGrid(kyArray, kzArray, sliceSpacing, bh);
title(sprintf('L = %G \n Packing density = %i%%',Loptions(iL),PackingDensityoptions(indOptions(iL))))
title(sprintf('L = %G %s = %G \n Packing density = %.1f%%', Loptions(iL), char(hex2dec('0394')), ...
nSlicesSimultaneously/Loptions(iL), PackingDensityoptions(indOptions(iL))))
if iL>1, ylabel([]); yticks([]); end
end
fh{2} = figure('Name','Packing density'); hold on;
plot(LArray, packingDensityArray, 'LineWidth', 3, 'Color', 'k', ...
'DisplayName', 'All feasible sampling patterns');
plot(L_CAIPI, packingDensityCAIPI, 'x', 'MarkerSize', 10, 'LineWidth', 3, ...
'DisplayName', 'Blipped-CAIPIRINHA - A');
if ~isempty(L_Narsude)
plot(L_Narsude, packingDensityNarsude, 'x', 'MarkerSize', 10, 'LineWidth', 3, ...
'Color', [148,0,211]/255, 'DisplayName', 'Blipped-CAIPIRINHA - B');
end
plot(L_optimum, packingDensity, 'x', 'MarkerSize', 10, 'LineWidth', 3, 'Color', ...
[34, 139, 34]/255, 'DisplayName', 'Optimum sampling');
plot(L_CAIPI_A, packingDensityCAIPI_A, 'x', 'MarkerSize', 10, 'LineWidth', 3, ...
'Color', colourCAIPI_A, 'DisplayName', 'Blipped-CAIPIRINHA - A');
plot(L_CAIPI_B, packingDensityCAIPI_B, 'x', 'MarkerSize', 10, 'LineWidth', 3, ...
'Color', colourCAIPI_B, 'DisplayName', 'Blipped-CAIPIRINHA - B');
plot(L_optimum, packingDensity, 'x', 'MarkerSize', 10, 'LineWidth', 3, ...
'Color', colourOptimum, 'DisplayName', 'Optimum sampling');
legend;
xlabel('L');
ylabel('Packing density [%]');
......
......@@ -15,7 +15,7 @@ FOVinplane = 220;
sliceSpacing = FOVinplane/2/Nsim;
[packingDensity, L_optimum, isDegenerate, packingDensityCAIPItmp, L_CAIPItmp, fh] = ...
[~, ~, ~, ~, ~, ~, ~, fh] = ...
ComputeMostHexagonalSampling(R, FOVinplane, sliceSpacing, Nsim, [], [], [], true, 0, [], 0.001);
if doSave
......
......@@ -23,14 +23,21 @@ t = tiledlayout(numel(NSimArray),1,'TileSpacing','Compact');
for iNSim = 1:numel(NSimArray)
Nsim = NSimArray(iNSim);
clear RArray packingDensity L_optimum packingDensityCAIPI L_CAIPI
clear RArray packingDensity L_optimum packingDensityCAIPI_A L_CAIPI_A packingDensityCAIPI_B L_CAIPI_B
RArray = (Rtot_start/Nsim):Rtot_increment:ceil(Rtot_end/Nsim);
sliceSpacing = FOVinplane/2/Nsim;
for iR = 1:numel(RArray)
[packingDensity(iR), L_optimum(iR), isDegenerate, packingDensityCAIPItmp, L_CAIPItmp] = ...
ComputeMostHexagonalSampling(RArray(iR), FOVinplane, sliceSpacing, Nsim);
[packingDensityCAIPI(iR), indCAIPI] = max(packingDensityCAIPItmp);
L_CAIPI(iR) = L_CAIPItmp(indCAIPI);
[packingDensity(iR), L_optimum(iR), isDegenerate, L_CAIPItmp_A, packingDensityCAIPItmp_A, ...
L_CAIPItmp_B, packingDensityCAIPItmp_B] = ...
ComputeMostHexagonalSampling(RArray(iR), FOVinplane, sliceSpacing, Nsim, [], [], [], false);
[packingDensityCAIPI_A(iR), indCAIPI] = max(packingDensityCAIPItmp_A);
L_CAIPI_A(iR) = L_CAIPItmp_A(indCAIPI);
if ~isempty(packingDensityCAIPItmp_B)
[packingDensityCAIPI_B(iR), indCAIPI] = max(packingDensityCAIPItmp_B);
L_CAIPI_B(iR) = L_CAIPItmp_B(indCAIPI);
end
if isDegenerate
break
end
......@@ -46,17 +53,15 @@ for iNSim = 1:numel(NSimArray)
nexttile; hold on;
plot(RArray(1:numel(packingDensity))*Nsim, packingDensity,...
'DisplayName', ['As hexagonal as possible' name], 'LineWidth', 2, 'Color', [34, 139, 34]/255);
plot(RArray(1:numel(packingDensity))*Nsim, packingDensityCAIPI,...
plot(RArray(1:numel(packingDensity))*Nsim, packingDensityCAIPI_A, 'Color', [217, 83, 25]/255,...
'DisplayName', ['Blipped-CAIPIRINHA - A' name], 'LineWidth', 2);
isNarsude = mod(L_CAIPI,1)>0;
plot(RArray(isNarsude)*Nsim, packingDensityCAIPI(isNarsude), '.',...
'DisplayName', ['Blipped-CAIPIRINHA - B' name], 'MarkerSize', 5, 'Color', [148,0,211]/255);
plot(RArray(1:min(numel(packingDensity),numel(packingDensityCAIPI_B)))*Nsim, packingDensityCAIPI_B, 'Color', [148, 0, 211]/255,...
'DisplayName', ['Blipped-CAIPIRINHA - B' name], 'LineWidth', 2);
yline(pi/2/sqrt(3)*100, 'Color', [0.5 0.5 0.5], 'LineWidth', 2, 'DisplayName',...
'Maximum circle packing density, 90.7%');
xlim([Rtot_start Rtot_end]);
ylim([60 92]);
ylim([50 92]);
if iNSim==numel(NSimArray) || doSave==false
xlabel('R_{total}');
legend;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment