Example1.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2004, Jonathan Stott
%
% 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 2
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tell matlab where to find the PMI files.  This should be done only
% once per session; I recommend putting it into your startup.m file.

pmipath('/homes/monte/1/home/jstott/matlab/PMI');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define a few variables for convenience
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

doscat = 1;		% flags, 1 or 0
doabs  = 1;
thick  = 6;		% Slab thickness

Method = 'Rytov';       % Born or Rytov

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate SD structure with source/detector/measurement information

SD.Lambda  = [690 800]; % Define 2 Wavelengths
SD.ModFreq = [  100  ]; % Single-frequency RF imager

% Source and Detector Positions

SD.SrcPos = SetOptode([-2:2:2], [-2:2:2], 0, 1);
SD.DetPos = SetOptode([-3:2:3], [-3:2:3], thick, 1);

% Source and Detector Amplitudes

SD.SrcAmp = 1e-3*ones(size(SD.SrcPos,1), ...
                      length(SD.Lambda), length(SD.ModFreq));
SD.DetAmp = 1e-3*ones(size(SD.DetPos,1), ...
                      length(SD.Lambda),length(SD.ModFreq));

% Generate a measurement list with all possible measurements

SD.MeasList = genMeasList(SD, 'all');

nWvl = length(SD.Lambda);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate structure describing Medium to be imaged

% Optical properties

Medium.idxRefr = [ 1.37 1.36 ];
Medium.Muao    = [ 0.02 0.04 ];
Medium.Muspo   = [ 9.50 8.35 ];

% Imaging volume

Medium.Geometry       = 'slab';
Medium.Slab_Thickness = thick;

% Volume to reconstruct

Medium.CompVol.Type  = 'computed';
Medium.CompVol.XStep = 0.5;
Medium.CompVol.X     = [ -4.0 4.0 ];
Medium.CompVol.YStep = 0.5;
Medium.CompVol.Y     = [ -4.0 4.0 ];
Medium.CompVol.ZStep = 0.5;	% Avoid the problematic z=0 case
Medium.CompVol.Z     = [ 0.25 thick ];

nVox = length(sampleVolume(Medium.CompVol));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate forward matrix ('Born'/'Rytov' solve the homogeneous problem)

% Two wavelengths were defined, but I'm only going to reconstruct data
%  from one of them, as an example of using a subset of the measurement
%  list (and to save time).  Also drop the long-distance imaging pairs.

dR = calcSep(SD, SD.MeasList);

ml1 = find(SD.MeasList(:,4) == 1 & dR < 9);
MeasList = SD.MeasList(ml1, :);
clear ml1 dR;

disp('Generating forward matrix - please wait');

[Phi0, A] = genBornMat(SD, Medium, MeasList, Method, [doabs doscat]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define perturbations - NOTE: abs/scat perturbations supported if and
%                         only if they're supported by the forward
%                         problem (i.e., they're included in A)

if (doabs)
   Medium.Object{1}.Type   = 'Sphere';         % Spherical scatterer
   Medium.Object{1}.Pos    = [ 1 1.25 2.9 ];
   Medium.Object{1}.Radius = 1.2;
   Medium.Object{1}.Mua    = Medium.Muao;
   Medium.Object{1}.Musp   = Medium.Muspo - 1;
end

if (doscat)
   if (doabs)
      n = 2;
   else
      n = 1;
   end
   
   Medium.Object{n}.Type   = 'Sphere';          % Spherical absorber
   Medium.Object{n}.Pos    = [ -2.25 -1.25 1.75 ];
   Medium.Object{n}.Radius = 1.2;
   Medium.Object{n}.Mua    = Medium.Muao + 0.02;
   Medium.Object{n}.Musp   = Medium.Muspo;
   
   clear n;
end

disp('Generating simulated data');

Phi = genBornData(SD, Medium, MeasList, Method, Phi0, A, [doabs doscat]);

disp('Adding noise');

Phi = addShotNoise(SD,MeasList,Phi,10,1);    	% approx 1000:1 SNR
Phi = addElecNoise(SD,MeasList,Phi,1000);    	% average 1000:1 SNR

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Invert the simulated data

disp('Inverting simulated data');

% A gives me d\mu_a and dD/D, convert to d\mu_a/\mu_a and dD/D to make
% both halves of order 1

if (doabs & doscat)
   for k = 1:nWvl
      % First half at each wavelength is absorbtion, second scattering
      
      A(:,2*(k-1)*nVox + [1:nVox]) = ...
	  A(:,2*(k-1)*nVox + [1:nVox]) * Medium.Muao(k);
   end
end

% Repack as real and imaginary for inversion

A1 = [ real(A); imag(A) ];

if (strcmpi(Method,'Rytov'))
   Y = [ real(log(Phi./Phi0)); imag(log(Phi./Phi0)) ];
else
   Y = [ real(Phi-Phi0); imag(Phi-Phi0) ];
end

clear A Phi0;

% Invert the data using filtered back-projection with Tikhonov
% regularization.  The regularization parameter alpha was adjusted 
% manually to get good images.

alpha = 1e-6 * normest(A1).^2;

X = fbp(A1, Y, alpha);

if (doabs & doscat)
   X = reshape(X, nVox, 2*nWvl);

   for k = 1:nWvl
      X(:,2*k-1) = X(:,2*k-1) *  Medium.Muao(k);
      X(:, 2*k ) = X(:, 2*k ) * -Medium.Muspo(k) / 3;
   end
else
   X = reshape(X, nVox, nWvl);
   
   if (doscat)
      for k = 1:nWvl
	 X(:,k) = X(:,k) * -Medium.Muspo(k) / 3;
      end
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display the reconstructed data [d\mu_a and d\mu_s']

if (doabs & doscat)
   dmua = X(:, 1:2:end);
   dmus = X(:, 2:2:end);
elseif (doabs)
   dmua = X;
   dmus = zeros(0,nWvl);
else
   dmua = zeros(0,nWvl);
   dmus = X;
end

clear X;

for k = 1:nWvl
   if (any(MeasList(:,4) == k))
      showImage(Medium, dmua(:,k), dmus(:,k));
   end
end

% Display the actual scattering perturbation(s), just as a reference

for k = 1:nWvl
   % Only display wavelengths actually used
   
   if (any(MeasList(:,4) == k))
      dmua0 =  calcDelMuA(Medium, [], k);
      dmus0 = calcDelMuSp(Medium, [], k);

      if (~isempty(dmua0) & all(dmua0(:) == 0))
	 dmua0 = [];
      end

      if (~isempty(dmus0) & all(dmus0(:) == 0))
	 dmus0 = [];
      end

      showImage(Medium, dmua0(:), dmus0(:));
   end
end

% Clear some variables before exiting to save memor

clear dmu0 k n nWvl nVox thick k alpha;

Original Perturbations
Reconstructed Perturbations