% ##########
% NOTES:
% Before clicking Run, please make sure to enter values for the following parameters
% into the 'Program input' field on the top-right:
% Conductivity [S/m]
% Permittivity
% Frequency [MHz]
% Volume [ml]
% This is a version of the phantom recipe generator found at:
% https: / / amri.ninds.nih.gov / phantomrecipe.html
%
% CREATED:
% 2025-03-27..04-21, JAdZ
%
% LICENSE:
% Copyright (C) <= 2025 NIH, NYU and others
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, see <http://[Log in to view URL]
%
% REFERENCE:
% C Ianniello, JA de Zwart, Q Duan, CM Deniz, L Alon, JS Lee, R Lattanzi, R Brown
% Synthesized tissue-equivalent dielectric phantoms using salt and polyvinylpyrrolidone
% solutions. Magn Reson Med 2018, 80:413 (doi: 10.1002/mrm.27005)
%
% MODIFICATION HISTORY:
% 2025-04-22, JAdZ
% Implemented the lowfrq code, but the correct transition frequency is unknown, I
% emailed Ryan about this.
% ##########
% ##### make this a script file #####
% see https://[Log in to view URL]
1;
% ##### PVP_fitting() #####
function [M_NaCl,M_PVP, M_water] = PVP_fitting(sig,epi_r, f, M_water)
if length(M_water(:)) == 1
M_water = M_water * ones(size(sig));
end
M_NaCl = zeros(size(sig));
M_PVP = M_NaCl;
for k = 1:length(sig(:))
[PVP_ratio, Salt_ratio] = solve_PVP_poly(sig(k),epi_r(k),f(k));
M_NaCl(k) = Salt_ratio*M_water(k);
M_PVP(k) = PVP_ratio*M_water(k);
end
endfunction
% ##### solve_PVP_poly() #####
function [PVP_ratio, Salt_ratio] = solve_PVP_poly(si,ei,fi)
tmp = fminsearch(@(x) err_fun(x,fi,si,ei), [0.65 0.0125]);
PVP_ratio = tmp(1);
Salt_ratio =tmp(2);
end
function errs = err_fun(x,fi,si,ei)
[s_fit,e_fit] = PVP_get_property(x(1), x(2), fi);
x1 = [ei,si];
x2 = [e_fit,s_fit];
diffs = (x1-x2)./x1;
errs = sqrt(sum(diffs.^2));
endfunction
% ##### ND_polyval() #####
function [Y,EM,Porder] = ND_polyval(p,X,varargin)
% 2016/08/31, Jacco A. de Zwart
%
% Copied from Qi Duan's SANBA_get_recipe_togo.m, since it is also needed in
% PVP_get_recipe.m
if(~isvector(p))
error('p must be a vector');
end
p=p(:);
Porder=[];
PM=[];
Nterm = length(p);
[Nsample,Ndim]=size(X);
switch nargin
case 2
%first, we need to figure out the order
Porder = ND_polyFindPorder(Nterm,Ndim);
%calculate the power matrix
PM = flipud(perms_sum([0 Porder],Ndim));
case 3
PM = varargin{1};
if (length(p)~=size(PM,1))
error('p and PM are not compatible');
end
if (size(X,2)~=size(PM,2))
error('X and PM are not compatible');
end
otherwise
error('Incorrect number of input variables');
end
buff = ones(Nterm,Ndim,Nsample);
%alternatively, presumable faster since usually Nsample>Nterm>Ndim
for l = 1:Ndim
dmatrix = X(:,l);
for k = 1:Nterm
pwr = dmatrix.^PM(k,l);
buff(k,l,:) = pwr;
end
end
%Y = sum(squeeze(prod(buff,2)).*repmat(p(:),[1,Nsample]),1);
%compute encoding matrix
EM = squeeze(prod(buff,2)); %Nterm x Nsample
EM = EM'; %Nsample x Nterm
Y = EM*p;
Y = Y(:);
endfunction
% ##### PVP_get_property() #####
function [sig,epi_r] = PVP_get_property(PVP_Water_Ratio, Salt_Water_Ratio, Freq)
if (length(PVP_Water_Ratio)~=length(Salt_Water_Ratio))||(length(PVP_Water_Ratio)~=length(Freq))
error(' PVP_get_property: Three inputs have to be the same size');
end
% Cannot load S from anywhere, hard-code locally instead
%S = load('PVP_poly');
S.PVP_scale=1.;
S.f_scale=1.e-9;
S.salt_scale=100;
S.PM_e=[2 2 3; 2 1 4; 2 0 5; 1 2 4; 1 1 5; 1 0 6; 0 2 5; 0 1 6; 0 0 7; 2 2 2; ...
2 1 3; 2 0 4; 1 2 3; 1 1 4; 1 0 5; 0 2 4; 0 1 5; 0 0 6; 2 2 1; 2 1 2; ...
2 0 3; 1 2 2; 1 1 3; 1 0 4; 0 2 3; 0 1 4; 0 0 5; 2 2 0; 2 1 1; 2 0 2; ...
1 2 1; 1 1 2; 1 0 3; 0 2 2; 0 1 3; 0 0 4; 2 1 0; 2 0 1; 1 2 0; 1 1 1; ...
1 0 2; 0 2 1; 0 1 2; 0 0 3; 2 0 0; 1 1 0; 1 0 1; 0 2 0; 0 1 1; 0 0 2; ...
1 0 0; 0 1 0; 0 0 1; 0 0 0];
S.PM_s=[2 2 3; 2 1 4; 2 0 5; 1 2 4; 1 1 5; 1 0 6; 0 2 5; 0 1 6; 0 0 7; 2 2 2; ...
2 1 3; 2 0 4; 1 2 3; 1 1 4; 1 0 5; 0 2 4; 0 1 5; 0 0 6; 2 2 1; 2 1 2; ...
2 0 3; 1 2 2; 1 1 3; 1 0 4; 0 2 3; 0 1 4; 0 0 5; 2 2 0; 2 1 1; 2 0 2; ...
1 2 1; 1 1 2; 1 0 3; 0 2 2; 0 1 3; 0 0 4; 2 1 0; 2 0 1; 1 2 0; 1 1 1; ...
1 0 2; 0 2 1; 0 1 2; 0 0 3; 2 0 0; 1 1 0; 1 0 1; 0 2 0; 0 1 1; 0 0 2; ...
1 0 0; 0 1 0; 0 0 1; 0 0 0];
if (Freq < 10.e9)
S.p_e=transpose([ 0.013557011 ; -0.11805195 ; 0.0076131663 ; 0.017256696 ; -0.0025769092 ; 0.0043416256 ; ...
0.0019328785 ; 0.012710021 ; -0.052021413 ; -0.11585241 ; 1.2706060 ; -0.065454174 ; -0.20849495 ; ...
0.10610378 ; -0.089433297 ; -0.041223019 ; -0.21115615 ; 0.87362259 ; 0.30538797 ; -4.6591597 ; ...
0.12849309 ; 0.87737280 ; -0.94482957 ; 0.70015815 ; 0.30637670 ; 1.4282038 ; -5.9062002 ; ...
-0.81108291 ; 8.2017696 ; -0.046788800 ; -1.6194214 ; 3.1846617 ; -2.7062164 ; -1.0086707 ; ...
-5.0432326 ; 20.553815 ; -3.1016459 ; 1.0955239 ; 1.9014745 ; -5.3279383 ; 6.3525175 ; 1.5715749 ; ...
9.9090984 ; -38.998115 ; 13.771371 ; 1.2597579 ; -13.529784 ; -1.2140807 ; -10.590180 ; 38.859382 ; ...
-45.968063 ; 4.1233002 ; -18.486451 ; 80.341140 ]);
S.p_s=transpose([ 0.00036017521 ; 0.0050755034 ; 0.00014456315 ; -0.00050220931 ; -0.00023334900 ; -0.00016861715 ; ...
5.2075024e-06 ; -3.2204460e-05 ; -0.00068067540 ; -0.00051555143 ; -0.036620046 ; -0.0059069587 ; ...
0.0029268847 ; -0.0028363596 ; 0.0019522635 ; 0.00043826712 ; 0.00062842420 ; 0.011035856 ; ...
-0.0010011897 ; 0.032691044 ; 0.049412973 ; -0.0026304978 ; 0.034005829 ; 0.0032538226 ; ...
-0.0033452600 ; -0.0046115161 ; -0.070069727 ; -0.15545854 ; -0.21589380 ; -0.074497826 ; ...
0.018787587 ; -0.052190288 ; -0.11182866 ; 0.0061203598 ; 0.018290554 ; 0.21801830 ; 1.8702903 ; ...
-0.0053568649 ; 0.25684365 ; 0.19235538 ; 0.19846533 ; -0.021200856 ; -0.065384737 ; -0.34914141 ; ...
-0.35021451 ; -3.3460019 ; 0.25250292 ; -0.10960254 ; 0.11795664 ; 0.53461132 ; 0.14660646 ; ...
1.6373667 ; -0.075738217 ; 0.10935135 ]);
else
S.p_e=transpose([ 0.011805335 ; -0.14464843 ; 0.0069357918 ; 0.013468636 ; -0.0059084284 ; 0.010146790 ; ...
0.00078412937 ; 0.0040864598 ; -0.015136267 ; -0.095121764 ; 1.5438998 ; -0.035230146 ; ...
-0.16863082 ; 0.22784519 ; -0.17019150 ; -0.025323176 ; -0.069588966 ; 0.25411302 ; 0.28599051 ; ...
-5.5439930 ; -0.14294256 ; 0.72280651 ; -1.9631792 ; 1.0517561 ; 0.22518252 ; 0.48922931 ; ...
-1.7269313 ; -0.082226131 ; 8.3682848 ; 0.85457320 ; -1.3976375 ; 6.3753179 ; -2.9454405 ; ...
-0.81507715 ; -1.8610658 ; 6.1196438 ; -4.7693457 ; 0.34965290 ; 0.51793472 ; -8.2129545 ; ...
4.5307713 ; 1.3414729 ; 4.2276048 ; -12.137641 ; 15.144375 ; 5.5387159 ; -10.556753 ; -0.53304948 ; ...
-5.8848859 ; 13.082945 ; -50.563847 ; 1.5125963 ; -7.1994831 ; 79.442701 ]);
S.p_s=transpose([ 0.00050292956 ; 0.0044025354 ; 0.00081702467 ; -0.00035627814 ; -0.00040857246 ; -0.00058451886 ; ...
8.6110845e-06 ; -1.5125634e-05 ; -0.00041641555 ; -0.0042955366 ; -0.033546825 ; -0.012998685 ; ...
0.0017412356 ; -0.00044430190 ; 0.0054871013 ; 0.00031705808 ; 0.00059590591 ; 0.0079475010 ; ...
0.00098822067 ; 0.027044634 ; 0.077991047 ; 0.0052494236 ; 0.024306980 ; -0.0031549297 ; ...
-0.0025991520 ; -0.0060862995 ; -0.058010820 ; -0.093563581 ; -0.11462494 ; -0.12628383 ; ...
0.0086047419 ; -0.026771431 ; -0.13466063 ; 0.0019014825 ; 0.027826888 ; 0.20578417 ; 1.7840338 ; ...
-0.0081733085 ; 0.15066118 ; 0.064693028 ; 0.29797058 ; -0.012923270 ; -0.083016310 ; -0.37726294 ; ...
-0.38131759 ; -3.2915352 ; 0.22245810 ; -0.055920294 ; 0.14949451 ; 0.59238066 ; 0.31376929 ; ...
1.5957381 ; -0.15722930 ; 0.050235394 ]);
end
sig = zeros(size(Freq));
epi_r = sig;
x=[PVP_Water_Ratio(:),Salt_Water_Ratio(:)];
sig(:) = ND_polyval(S.p_s,[x(:,1)*S.PVP_scale,x(:,2)*S.salt_scale, Freq(:)*S.f_scale],S.PM_s);
epi_r(:) = ND_polyval(S.p_e,[x(:,1)*S.PVP_scale,x(:,2)*S.salt_scale, Freq(:)*S.f_scale],S.PM_e);
endfunction
% ##### MAIN #####
fprintf(1,"Make sure to enter values for the following in the 'Program input' box\n");
fprintf(1,"directly above this text (newline separated):\n");
fprintf(1," Conductivity [S/m]\n");
fprintf(1," Permittivity\n");
fprintf(1," Frequency [MHz]\n");
fprintf(1," Volume [ml]\n");
fprintf(1,"\n");
cond=input("Conductivity [S/m]: ")
perm=input("Permittivity : ")
freq=input("Frequency [MHz] : ")
volu=input("Volume [ml] : ")
%[nacl,pvp,water] = PVP_fitting(0.6,50.,297400000.,1000.)
[nacl,pvp,water] = PVP_fitting(cond,perm,freq*1.e6,volu);
fprintf(1,"\n");
fprintf(1,"REQUIRED INGREDIENTS:\n")
fprintf(1,"NaCl: %f g\n",nacl);
fprintf(1,"PVP: %f g\n",pvp);
fprintf(1,"Water: %f g\n",water);
fprintf(1,"\n");
To embed this project on your website, copy the following code and paste it into your website's HTML: