Merge branch 'master' into 'main'

Streamline packing density computation;

See merge request sapme1/hexagonalsampling!2

ComputeMostHexagonalSampling.m

-function [packingDensity, L_optimum, isDegenerate, packingDensityCAIPI, L_CAIPI] = ...
+function [packingDensity, L_optimum, isDegenerate, packingDensityCAIPI, L_CAIPI, fh] = ...
     ComputeMostHexagonalSampling(R_inplane, FOV_inplane, sliceSpacing, nSlicesSimultaneously, ...
-    nSlicesTotal, sliceThickness, sliceGap, doPlot, doDebug, gifName)
+    nSlicesTotal, sliceThickness, sliceGap, doPlot, doDebug, gifName, Lincrement)
 % This function computes how to optimally sample the phase encoding plane for 3D/Multiband
 % EPI/Spiral to achieve the maximum packing density of circles in the image domain
 %
@@ -34,6 +34,9 @@ function [packingDensity, L_optimum, isDegenerate, packingDensityCAIPI, L_CAIPI]
 %
 %   gifName                 name of gif in case debug plot should be saved as gif
 %
+%   Lincrement              A number of Ls is tested, Lincrement is the resolution, i.e. how precise
+%                           is the determination of the optimum L {0.01}
+%
 %
 % OUT
 %   packingDensity          circle packing density [%]
@@ -43,6 +46,9 @@ function [packingDensity, L_optimum, isDegenerate, packingDensityCAIPI, L_CAIPI]
 %
 %   isDegenerate            flag whether packing density metric is degenerate
 %
+%   fh                      cell figure handles to grid figure (fh{1}) and packing density
+%                           figure (fh{2})
+%
 %
 % Example: [packingDensity, L_optimum] = ComputeMostHexagonalSampling(1, 192, 19.2, 3);
 %          [packingDensity, L_optimum] = ComputeMostHexagonalSampling(1, 192, [], 3, 3, 4, 15.2)
@@ -54,11 +60,11 @@ if nargin == 7
     sliceSpacing = nSlicesTotal*(sliceThickness+sliceGap)/nSlicesSimultaneously;
 end
 
-if nargin < 8
+if nargin < 8 || isempty(doPlot)
     doPlot = false;
 end
 
-if nargin < 9
+if nargin < 9 || isempty(doDebug)
     doDebug = false;
 end
 
@@ -74,9 +80,17 @@ kzSlabThickness = 2*pi/sliceSpacing;
 % the next point after a down-blip, the grids are going to get only worse = less homogeneous, i.e.
 % 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)
+    Lmax = max(Lmax,nSlicesSimultaneously);
+end
 % For L<2 there is an L>2 which produces a rotationally symmetric grid i.e. same packing density
 Lmin = 2;
-Lincrement = 0.01;%0.001;
+
+if nargin < 11
+    Lincrement = 0.01;
+end
 
 LArray = Lmin:Lincrement:Lmax;
 
@@ -88,7 +102,7 @@ DELTA = 1:nSlicesSimultaneously;
 LNarsudeUnique = nSlicesSimultaneously./DELTA;
 LNarsudeUnique = LNarsudeUnique(mod(LNarsudeUnique,1)~=0);
 LNarsudeUnique = LNarsudeUnique(LNarsudeUnique>2);
-LArray = sort([LArray LNarsudeUnique]);
+LArray = unique(sort([LArray LNarsudeUnique]));
 
 howFar2Look = 10; % TODO: Find out what to put here!
 d_min = zeros(size(LArray));
@@ -110,7 +124,7 @@ for iL = 1:numel(LArray)
     d = sqrt(kyArray.^2 + kzArray.^2);
     d_min(iL) = min(d(d>0));
 
-    packingDensityArray(iL) = pi/4*d_min(iL)^2/(2*pi*deltak_inplane/sliceSpacing);
+    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);
@@ -120,18 +134,19 @@ for iL = 1:numel(LArray)
         PlotGrid(kyArray, kzArray, sliceSpacing, ah);
 
         subplot(1,2,2); hold on;
-        plot(LArray(iL),packingDensityArray(iL)*100,'.','Color','k');
+        plot(LArray(iL),packingDensityArray(iL),'.','Color','k');
 
         if isBlippedCAIPI(iL) || isNarsude(iL)
-            plot(LArray(iL),packingDensityArray(iL)*100,'x','Color','r','MarkerSize',10,'LineWidth',3);
+            plot(LArray(iL),packingDensityArray(iL),'x','Color','r','MarkerSize',10,'LineWidth',3);
         end
         
         xlabel('L'); ylabel('Packing density [%]');
-        xlim([Lmin Lmax]); ylim([min(packingDensityArray*100) 90.7]);
+        xlim([Lmin Lmax]); ylim([min(packingDensityArray) 90.7]);
         sgtitle(sprintf('L = %.2f',LArray(iL)));
 
         %% Make GIF
         if ~isempty(gifName)
+            set(gcf, 'color', 'white');
             frame = getframe(1);
             im = frame2im(frame);
             [imind,cm] = rgb2ind(im,256);
@@ -158,17 +173,38 @@ 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);
 
 
 %% Plot grid and resulting sampling density for all options
 if doPlot
-    % For plotting get grid points on a bit larger range
-    [kyArray, kzArray] = GetGridPoints(2*ceil(Lmax), deltak_inplane, sliceSpacing, L_optimum);
-    PlotGrid(kyArray, kzArray, sliceSpacing);
+    fh{1} = figure('Name','Sampling grid');
+    [Loptions, indOptions] = sort([L_optimum L_CAIPI]);
+    PackingDensityoptions = round([packingDensity packingDensityCAIPI]);
+    for iL = 1:numel(Loptions)
+        bh = subplot(1,numel(L_CAIPI)+1,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))))
+        if iL>1, ylabel([]); yticks([]); end
+    end
 
-    figure('Name','Packing density'); hold on;
-    plot(LArray, packingDensityArray*100);
-    plot(L_CAIPI, packingDensityCAIPI*100, 'x', 'MarkerSize', 10, 'LineWidth', 3);
+    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');
+    legend;
     xlabel('L');
     ylabel('Packing density [%]');
 end
@@ -194,13 +230,16 @@ scatter(ky, kz, 'MarkerEdgeColor', [0.5 0.5 0.5], 'MarkerFaceColor', [0.5 0.5 0.
 
 hold on;
 
-yline(0,'r','LineWidth',3); 
-yline(2*pi/sliceSpacing,'r','LineWidth',3);
+slabBorderColour = [100 149 237]/255; % light blue
+yline(0,'r','LineWidth',3, 'Color', slabBorderColour); 
+yline(2*pi/sliceSpacing, 'r', 'LineWidth', 3, 'Color', slabBorderColour);
+
+plot(ky(1:ceil(numel(ky)/3)), kz(1:ceil(numel(ky)/3)), 'LineWidth', 3, 'Color', [253 218 13]/255); % yellow
 
 axis equal;
-ylim(2*pi/sliceSpacing*[-1 2]);
+ylim(2*pi/sliceSpacing*[-0.2 1.2]);
 xlim([min(ky) max(ky)]);
-xlabel('k_y'); ylabel('k_z');
+xlabel('k_y [rad/m]'); ylabel('k_z [rad/m]');
 hold off;
 end
 
@@ -208,5 +247,5 @@ end
 function packingDensity = ComputePackingDensity(d, deltak_inplane, sliceSpacing)
 A_circle = pi/4*d.^2;
 A_unitcell = 2*pi*deltak_inplane/sliceSpacing;
-packingDensity = A_circle/A_unitcell;
+packingDensity = A_circle/A_unitcell*100;
 end
\ No newline at end of file

MakeExampleGridsPlot.m

+clear all;
+
+baseFolder = pwd;
+
+doSave = true;
+saveFolder = pwd;
+if doSave
+    dfws = get(0,'DefaultFigureWindowStyle');
+    set(0,'DefaultFigureWindowStyle','normal');
+end
+
+Nsim = 5;
+R = 4.4/Nsim;
+FOVinplane = 220;
+
+sliceSpacing = FOVinplane/2/Nsim;
+
+[packingDensity, L_optimum, isDegenerate, packingDensityCAIPItmp, L_CAIPItmp, fh] = ...
+    ComputeMostHexagonalSampling(R, FOVinplane, sliceSpacing, Nsim, [], [], [], true, 0, [], 0.001);
+
+if doSave
+    set(fh{1},'Position',[29 174 1479 588])
+    export_fig(fullfile(saveFolder, 'ExampleGrids'),'-png','-transparent','-a1','-r300',fh{1});
+
+    set(fh{2},'Position',[845.8  418.6  372.8  244.8])
+    hLegend = findobj(gcf, 'Type', 'Legend');
+    hLegend.Position = [0.2517    0.2370    0.5359    0.2708];
+    export_fig(fullfile(saveFolder, 'ExampleGridsPackingDensity'),'-png','-transparent','-a1','-r300',fh{2});
+    set(0,'DefaultFigureWindowStyle',dfws);
+end
\ No newline at end of file

MakeOptDvsRPlot.m

@@ -44,11 +44,17 @@ for iNSim = 1:numel(NSimArray)
             name = sprintf(', N_{sim} = %i',Nsim);
         end
         nexttile;  hold on;
-        plot(RArray(1:numel(packingDensity))*Nsim,packingDensity*100,...
-            'DisplayName', ['As hexagonal as possible' name],'LineWidth',2,'Color','k');
-        plot(RArray(1:numel(packingDensity))*Nsim,packingDensityCAIPI*100,...
-            'DisplayName',['Blipped CAIPIRINHA' 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%');
+        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,...
+            '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);
+
+        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]);
         if iNSim==numel(NSimArray) || doSave==false
@@ -57,6 +63,10 @@ for iNSim = 1:numel(NSimArray)
         else
             xlabel(''); xticks('');
         end
+
+        if NSimArray(iNSim)==5
+            xline(4.4, ':', 'LineWidth', 2); % Mark example cases from other figure
+        end
     end
 end