This project has moved. For the latest updates, please go here.

New frequency domain button on the fluence/PhD maps - MATLAB yet?

Aug 11, 2013 at 12:16 AM
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 ...

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 carriage-return 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.70195812655052E-05    4.9729648657582E-05  ...
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.
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); 
                fs =  Vts.Modeling.ForwardSolvers.PointSourceSDAForwardSolver();
            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);
            % 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;
            r = complex(fReal, fImag);
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); % s-d separation, in mm
zs = linspace(0.1,19.9,100); % z range in mm
fts = linspace(0,1,2); % frequency range in GHz

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!
Aug 11, 2013 at 5:10 PM
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