
Hi  thanks to the VP team for adding the magic frequency domain button to the fluence maps page. I am finding this very helpful.
Question: Can this feature be replicated in the MATLAB version?
I ask not only because this would be cool, but to export numerical data in the VTS web page doesnt seem to work well. you can export the image fine, but trying to read that data for plotting etc fails. Having all of this in MATLAB makes for easy analysis ...
Thanks!


Coordinator
Aug 11, 2013 at 12:58 AM
Edited Aug 11, 2013 at 12:59 AM

Sorry you're having trouble, Albert.
Regarding data export in ASCII format, it looks like we added the line "% Map Values: " at the beginning of the data. Originally, you could just load this file as ascii data, with the x and z data points in the header (mainly for your reference).
If you open the text file, and hit a carriagereturn after "% Map Values: ", then you should have the 2D data available for matlab import:
% X Values: 19.9 19.7 19.5 ...
% Y Values: 0.1 0.3 0.5 0.7 0.9 ...
% Map Values:
4.70195812655052E05 4.9729648657582E05 ...
Then, in Matlab, you can write:
x = 19.9:0.1:20;
z = 0.1:0.2:19.9;
temp = load('test.txt', 'ascii');
fluence = reshape(temp,[length(x)/2 length(z)])';
figure; imagesc(log(fluence)); axis equal; axis off; colormap hot; colorbar;
This should let you plot something (with the default frequency domain settings) that looks like this:
Regarding direct Matlab computations  you're right; it looks like FluenceOfRhoAndZAndFt not in there yet. This should be pretty straightforward  would be happy to work on this.


Coordinator
Aug 11, 2013 at 1:34 AM

Albert, I've just added the following to the source code, which will go out in the next release.
To vts\VtsSolvers.m, add the following function, say at line 254 (after FluenceOfRhoAndZ):
function r = FluenceOfRhoAndZAndFt(op, rhos, zs, fts)
%% FluenceOfRhoAndZ
% FluenceOfRhoAndZ(OP, RHOS, ZS)
%
% OP is an N x 4 matrix of optical properties
% eg. OP = [[mua1, mus'1, g1, n1]; [mua2, mus'2, g2, n2]; ...];
% RHO is an 1 x M array of detector locations (in mm)
% eg. RHO = [1:10];
% Z is a 1 x M array of z values (in mm)
% eg. Z = linspace(0.1,19.9,100);
nop = size(op,1);
nrho = length(rhos);
nft = length(fts);
nz = length(zs);
if strfind(VtsSolvers.Options.SolverType, 'SDA')
fs = Vts.Factories.SolverFactory.GetForwardSolver(VtsSolvers.Options.SolverType);
else
fs = Vts.Modeling.ForwardSolvers.PointSourceSDAForwardSolver();
end
op_net = NET.createArray('Vts.OpticalProperties', nop);
for i=1:nop
op_net(i) = Vts.OpticalProperties;
op_net(i).Mua = op(i,1);
op_net(i).Musp = op(i,2);
op_net(i).G = op(i,3);
op_net(i).N = op(i,4);
end;
% call the solver, which returns an array of (.NET) Complex structs
fComplexNET = fs.FluenceOfRhoAndZAndFt(op_net,rhos,zs,fts);%,[length(fts) length(zs) length(rhos) nop]);
% create a native Matlab array
fReal = zeros([nft nz nrho nop]);
fImag = zeros([nft nz nrho nop]);
ci = 1;
for i=1:nop
for j=1:nrho
for k=1:nz
for el=1:nft
fReal(el,k,j,i) = fComplexNET(ci).Real;
fImag(el,k,j,i) = fComplexNET(ci).Imaginary;
ci=ci+1;
end
end
end
end
r = complex(fReal, fImag);
end
Then, in 'vts_solver_demo.m' (or in your own scripting code), here's the demo of the functionality:
%% Example FluenceOfRhoAndZAndFt
% Evaluate fluence as a function of rho and z using one set of optical
% properties and a distributed gaussian source SDA solver type.
op = [0.01 1 0.8 1.4];
rhos = linspace(0.1,19.9,100); % sd separation, in mm
zs = linspace(0.1,19.9,100); % z range in mm
fts = linspace(0,1,2); % frequency range in GHz
VtsSolvers.SetSolverType('DistributedPointSourceSDA');
test = VtsSolvers.FluenceOfRhoAndZAndFt(op, rhos, zs, fts);
xs = [fliplr(rhos(2:end)),rhos];
% xs = [rhos(end:1:2), rhos];
% f = figure; imagesc(log(test));
f = figure; imagesc(xs,zs,log([fliplr(squeeze(test(1,:,2:end))),squeeze(test(1,:,:))]));
axis image
title('Fluence of \rho and z and ft (ft=0GHz)');
xlabel('\rho [mm]')
ylabel('z [mm]')
set(f,'Name','Fluence of Rho and z and ft (ft=0GHz)');
f = figure; imagesc(xs,zs,log([fliplr(squeeze(test(2,:,2:end))),squeeze(test(2,:,:))]));
axis image
title('Fluence of \rho and z and ft (ft=1GHz)');
xlabel('\rho [mm]')
ylabel('z [mm]')
set(f,'Name','Fluence of Rho and z and ft (ft=1GHz)');
modulation = squeeze(test(2,:,:)./test(1,:,:));
figure; imagesc(xs,zs,[fliplr(modulation(:,2:end)), modulation]);
axis image
title('Modulation of fluence (AC/DC) of \rho and z and ft (ft=1GHz)');
xlabel('\rho [mm]')
ylabel('z [mm]')
set(f,'Name','Modulation of fluence (AC/DC) of Rho and z and ft (ft=1GHz)');
You should see a few images like this:
...and this:
Let me know how that goes for you!



thanks David  I think that did the trick.
hopefully sometime there will be an equivalent for the time domain (to explore time gating strategies).
*******************************************
Albert Cerussi, PhD
Director, Diffuse Optical Spectroscopic Imaging (DOSI) Lab
Beckman Laser Institute and Medical Clinic
University of California, Irvine
*******************************************

