% Radon transform example. This version uses
% an image containing smooth features
% Note: The inversion of the Radom Transform
% requires working with 2D Fourier transforms,
% which is nitty-gritty stuff and requires
% extensive knowledge of their properties
% image mtrue(i,j) is N by N, in (x=i,y=j) space,
% range of both x, y is +/- 1
% test image is a sum of several Normal curves
mtrue = exp( -( ((x-xbar1).^2)*ones(1,Ny) + ones(Nx,1)*((y'-ybar1).^2))/(2*s1*s1));
mtrue = mtrue+exp( -( ((x-xbar2).^2)*ones(1,Ny) + ones(Nx,1)*((y'-ybar2).^2))/(2*s2*s2));
mtrue = mtrue+exp( -( ((x-xbar3).^2)*ones(1,Ny) + ones(Nx,1)*((y'-ybar3).^2))/(2*s3*s3));
% some statistics, to compare with reconstructed image
fprintf('m has minimum value %f \n', min(min(mtrue)) );
m has minimum value 0.000000
fprintf('m has maximum value %f \n', max(max(mtrue)) );
m has maximum value 1.009415
fprintf('Caption: True image\n');
% radon transform R(i,j) is in (u=i, q=theta=j) space
% range of u is +/- 1 (note corners of box cur off)
Dq = (qmax-qmin)/(Nq-1); % note includes both -pi/2 and pi/2
q = qmin + Dq*[0:Nq-1]'; % for ease of interpolation
% integration over s, range +/- 1
% brute force Radon transform via ray sums
% object here is top build Nu by Ns array of image, m,
% evaluated at proper values of (u, s)
% (x,y) from (u,s); these are Nu by Ns matrices
xp = u*sin(q(i))*ones(1, Ns) - ones(Nu,1)*s'*cos(q(i))';
yp = u*cos(q(i))*ones(1, Ns) + ones(Nu,1)*s'*sin(q(i));
% indices corresponding to (x, y)
ixp = floor((xp-xmin)/Dx) + 1;
iyp = floor((yp-ymin)/Dy) + 1;
% handle out of range x indices
% handle out of range y indices
% nasty array manipulation to sample image, mtrue,
% at a series of (u, s) points, putting the values
% in an Nu by Ns matrix, S
% Step 1: convert Nu by Ns arrays to vectors
Lixp = ixp( [1:Nu*Ns]' );
Liyp = iyp( [1:Nu*Ns]' );
% Step 2: generate linear index into image
L = sub2ind( [Nx, Ny], Lixp, Liyp);
% Step 3: sample image, mtrue, and rearrange it back into Nu by Ns
S = reshape(mtrue(L),Nu,Ns);
% now handle points that were off the edge of (x,y) space by
% integrate (sum) S(u,s) over s to get R(u) for fixed q
% plot the radon transform
%title('Radon transform');
fprintf('Caption: Radon transform of image.\n');
Caption: Radon transform of image.
% statistics of the Radon transform
fprintf('R has minimum value %f \n', min(min(R)) );
R has minimum value 0.000000
fprintf('R has maximum value %f \n', max(max(R)) );
R has maximum value 0.375994
% Inverse Radon Transform via Fourier transforms
% reorder array so that s=0, currenty at row N/2, is at row 1.
R = circshift(R, [-Nx/2+1,0] );
% now take 1D fourier transform over columns of R
% invert Radon transform via central slice theorem
% Fourier transform of image
% fourier transform u -> ku
ku=Dku*[0:Nu/2,-Nu/2+1:-1]';
Nku2 = Nu/2+1; % non-negative frequuencies
% fourier transform x -> kx
kx=Dkx*[0:Nx/2,-Nx/2+1:-1]';
Nkx2 = Nx/2+1; % non-negative frequuencies
% fourier transform y -> ky
ky=Dky*[0:Ny/2,-Ny/2+1:-1]';
Nky2 = Ny/2+1; % non-negative frequuencies
% Now use the Fourier Slice theorem. All the
% work is in sampling the radon transform at
% a series of points that correspond to wavenumber
% pairs on a rectangular grid. That requires
% the Radon tranform to be interpolated. I
% provide two ways of performing the interpolation.
% Note: Method 1 seems to work better. I think
% that's because my complex interpolator works
% better than using MatLab's TriScattered Interp
% on the real and imaginary parts, separately
if( 1 ) % Method 1: element by element
for j = (1:Nky2) % only need left side, will rebuild right using symmetry
theta = atan2( kx(i), ky(j) );
if( theta > (pi/2) ) % map large thetas onto negative thetas
kr = sqrt(kx(i)^2 + ky(j)^2);
% index into randon transform array
continue; % ignore points that are out of bounds
% treat indices are real numbers
krindex = (kr-kumin)/Dku+1;
thindex = (theta-qmin)/Dq+1;
% bracket real number indices with integers
% Lagrangian polynomial interpoaltion
% note that the indices are separated by unity
CKR = ((B-krindex)*Rt(A,C)+(krindex-A)*Rt(B,C));
DKR = ((B-krindex)*Rt(A,D)+(krindex-A)*Rt(B,D));
else % Method 2: all at once
% Rt is defined art these ku's and q's
% Rtinterpr = TriScatteredInterp( kkuu(:), qq(:), real(Rt(:)), 'natural' );
% Rtinterpi = TriScatteredInterp( kkuu(:), qq(:), imag(Rt(:)), 'natural' );
% MatLab's TriScatteredInterp() interpolator can't handle complex
% numbers, so interpolate real and imaginary parts separately. The
% Delaunay triangulation needs to be performed only once, so do it
triangles = DelaunayTri( kkuu(:), qq(:) );
Rtinterpr = TriScatteredInterp( triangles, real(Rt(:)), 'natural' );
Rtinterpi = TriScatteredInterp( triangles, imag(Rt(:)), 'natural' );
% need all these elements of of mtt(i,j)
% ii = [1:Nx]'*ones(1,Nky2);
% jj = ones(Nx,1)*[1:Nky2];
% these indices have kr and theta
kkrr = sqrt( (kx(1:Nx)*ones(1,Nky2)).^2 + (ones(Nx,1)*ky(1:Nky2)').^2 );
qq = atan2( kx(1:Nx)*ones(1,Nky2), ones(Nx,1)*ky(1:Nky2)' );
% interpolate, set NaN's to zero, and populate mtt
tmp = complex( Rtinterpr( kkrr(:), qq(:) ), Rtinterpi( kkrr(:), qq(:) ) );
mtt(1:Nx, 1:Nky2) = reshape( tmp, Nx, Nky2 );
% Fourier transform of image needs to be phase shifted
% else origin is in the wrong place
phase = 2*pi*(kx*ones(1,Ny)*xshift + ones(Nx,1)*(ky')*yshift);
mtt = mtt .* complex( cos(phase), sin(phase) );
% impose all necessary symmetries required for a real image
% these four elements must be real
mtt(Nx/2+1,1)=real(mtt(Nx/2+1,1));
mtt(1,Ny/2+1)=real(mtt(1,Ny/2+1,1));
mtt(Nx/2+1,Ny/2+1)=real(mtt(Nx/2+1,Ny/2+1));
% bottoms of these two columns must be conjugates of top
mtt(Nx:-1:Nx/2+2,1)=conj(mtt(2:Nx/2,1));
mtt(Nx:-1:Nx/2+2,Ny/2+1)=conj(mtt(2:Nx/2,Ny/2+1));
% right hand side flipped conjugate of left hand side
mtt(1,mp) = conj(mtt(1,m));
mtt(2:Nx,mp) = flipud(conj(mtt(2:Nx,m)));
% now take inverse fourier transform
% Note: normalization not necessarily correct for non-square images
mest = real(mest); % throw away imaginary part (which should be zero)
% plot reconstructed image
% title('reconstructed image');
fprintf('Caption: Image reconstructed from Radon transform of image.\n');
Caption: Image reconstructed from Radon transform of image.
% compare amplitudes of reconstructed image to original,
% to be sure that the normalization is correct
fprintf('mest has minimum value %f \n', min(min(mest)) );
mest has minimum value -0.008493
fprintf('mest has maximum value %f \n', max(max(mest)) );
mest has maximum value 0.996005
% Radon transform example. This version uses
% Note: The inversion of the Radom Transform
% requires working with 2D Fourier transforms,
% which is nitty-gritty stuff and requires
% extensive knowledge of their properties
% image mtrue(i,j) is N by N, in (x=i,y=j) space,
% range of both x, y is +/- 1
% test image is a tiff file with a hypothetical depiction of a magma
% chamber; Image must be 256 by 256 in dimension to work.
mtrue=double(imread('../data/magma.tif','tif'));
mtrue = -(mtrue-mean(mean(mtrue)));
% some statistics, to compare with reconstructed image
fprintf('m has minimum value %f \n', min(min(mtrue)) );
m has minimum value -121.647095
fprintf('m has maximum value %f \n', max(max(mtrue)) );
m has maximum value 89.352905
imagesc([-1, 1], [-1, 1], mtrue);
fprintf('Caption: True image.\n');
% radon transform R(i,j) is in (u=i, q=theta=j) space
% range of u is +/- 1 (note corners of box cur off)
Dq = (qmax-qmin)/(Nq-1); % note includes both -pi/2 and pi/2
q = qmin + Dq*[0:Nq-1]'; % for ease of interpolation
% integration over s, range +/- 1
% radon transform R(u,q) tvia ray sums
% object here is top build Nu by Ns array of image, m,
% evaluated at proper values of (u, s)
% (x,y) from (u,s); these are Nu by Ns matrices
xp = u*sin(q(i))*ones(1, Ns) - ones(Nu,1)*s'*cos(q(i))';
yp = u*cos(q(i))*ones(1, Ns) + ones(Nu,1)*s'*sin(q(i));
% indices corresponding to (x, y)
ixp = floor((xp-xmin)/Dx) + 1;
iyp = floor((yp-ymin)/Dy) + 1;
% handle out of range x indices
% handle out of range y indices
% nasty array manipulation to sample image, mtrue,
% at a series of (u, s) points, putting the values
% in an Nu by Ns matrix, S
% Step 1: convert Nu by Ns arrays to vectors
Lixp = ixp( [1:Nu*Ns]' );
Liyp = iyp( [1:Nu*Ns]' );
% Step 2: generate linear index into image
L = sub2ind( [Nx, Ny], Lixp, Liyp);
% Step 3: sample image, mtrue, and rearrange it back into Nu by Ns
S = reshape(mtrue(L),Nu,Ns);
% now handle points that were off the edge of (x,y) space by
% integrate (sum) S(u,s) over s to get R(u) for fixed q
% plot the radon transform
imagesc([-1, 1], [-1, 1], R);
% title('Radon transform');
fprintf('Caption: Radon transform of true image.\n');
Caption: Radon transform of true image.
% statistics of the Radon transform
fprintf('R has minimum value %f \n', min(min(R)) );
R has minimum value -46.525462
fprintf('R has maximum value %f \n', max(max(R)) );
R has maximum value 76.184901
% Inverse Radon Transform via Fourier transforms
% reorder array so that s=0, currenty at row N/2, is at row 1.
R = circshift(R, [-Nx/2+1,0] );
% now take 1D fourier transform over columns of R
% invert radon transform via central slice theorem
% fourier transform of image
% fourier transform u -> ku
ku=Dku*[0:Nu/2,-Nu/2+1:-1]';
Nku2 = Nu/2+1; % non-negative frequuencies
% fourier transform x -> kx
kx=Dkx*[0:Nx/2,-Nx/2+1:-1]';
Nkx2 = Nx/2+1; % non-negative frequuencies
% fourier transform y -> ky
ky=Dky*[0:Ny/2,-Ny/2+1:-1]';
Nky2 = Ny/2+1; % non-negative frequuencies
% Now use the Fourier Slice theorem. All the
% work is in sampling the radon transform at
% a series of points that correspond to wavenumber
% pairs on a rectangular grid. That requires
% the Radon tranform to be interpolated. I
% provide two ways of performing the interpolation.
% Note: Method 1 seems to work better. I think
% that's because my complex interpolator works
% better than using MatLab's TriScattered Interp
% on the real and imaginary parts, separately
if( 0 ) % Method 1: element by element
for j = (1:Nky2) % only need left side, will rebuild right using symmetry
theta = atan2( kx(i), ky(j) );
if( theta > (pi/2) ) % map large thetas onto negative thetas
kr = sqrt(kx(i)^2 + ky(j)^2);
% index into randon transform array
continue; % ignore points that are out of bounds
% treat indices are real numbers
krindex = (kr-kumin)/Dku+1;
thindex = (theta-qmin)/Dq+1;
% bracket real number indices with integers
% Lagrangian polynomial interpoaltion
% note that the indices are separated by unity
CKR = ((B-krindex)*Rt(A,C)+(krindex-A)*Rt(B,C));
DKR = ((B-krindex)*Rt(A,D)+(krindex-A)*Rt(B,D));
else % Method 2: all at once
% Rt is defined art these ku's and q's
% Rtinterpr = TriScatteredInterp( kkuu(:), qq(:), real(Rt(:)), 'natural' );
% Rtinterpi = TriScatteredInterp( kkuu(:), qq(:), imag(Rt(:)), 'natural' );
% MatLab's TriScatteredInterp() interpolator can't handle complex
% numbers, so interpolate real and imaginary parts separately. The
% Delaunay triangulation needs to be performed only once, so do it
triangles = DelaunayTri( kkuu(:), qq(:) );
Rtinterpr = TriScatteredInterp( triangles, real(Rt(:)), 'natural' );
Rtinterpi = TriScatteredInterp( triangles, imag(Rt(:)), 'natural' );
% need all these elements of of mtt(i,j)
% ii = [1:Nx]'*ones(1,Nky2);
% jj = ones(Nx,1)*[1:Nky2];
% these indices have kr and theta
kkrr = sqrt( (kx(1:Nx)*ones(1,Nky2)).^2 + (ones(Nx,1)*ky(1:Nky2)').^2 );
qq = atan2( kx(1:Nx)*ones(1,Nky2), ones(Nx,1)*ky(1:Nky2)' );
% interpolate, set NaN's to zero, and populate mtt
tmp = complex( Rtinterpr( kkrr(:), qq(:) ), Rtinterpi( kkrr(:), qq(:) ) );
mtt(1:Nx, 1:Nky2) = reshape( tmp, Nx, Nky2 );
% Fourier transform of image needs to be phase shifted
% else origin is in the wrong place
phase = 2*pi*(kx*ones(1,Ny)*xshift + ones(Nx,1)*(ky')*yshift);
mtt = mtt .* complex( cos(phase), sin(phase) );
% impose all necessary symmetries required for a real image
% these four elements must be real
mtt(Nx/2+1,1)=real(mtt(Nx/2+1,1));
mtt(1,Ny/2+1)=real(mtt(1,Ny/2+1,1));
mtt(Nx/2+1,Ny/2+1)=real(mtt(Nx/2+1,Ny/2+1));
% bottoms of these two columns must be conjugates of top
mtt(Nx:-1:Nx/2+2,1)=conj(mtt(2:Nx/2,1));
mtt(Nx:-1:Nx/2+2,Ny/2+1)=conj(mtt(2:Nx/2,Ny/2+1));
% right hand side flipped conjugate of left hand side
mtt(1,mp) = conj(mtt(1,m));
mtt(2:Nx,mp) = flipud(conj(mtt(2:Nx,m)));
% now take inverse fourier transform
% Note: normalization not necessarily correct for non-square images
imagesc([-1, 1], [-1, 1], mest);
title('reconstructed image');
fprintf('Caption: Image reconstructed from Radon trsnsform.\n');
Caption: Image reconstructed from Radon trsnsform.
fprintf('mest has minimum value %f \n', min(min(mest)) );
mest has minimum value -99.766491
fprintf('mest has maximum value %f \n', max(max(mest)) );
mest has maximum value 94.227659
% gradient descent inversion of d(x)=Lm(x),
% where L is a linear operator. The gradient of
% error dE/dm is calculated using the adjoint method
% The linear operator is 0.1*(D + 2*I) where D is a
% first derivative and I is an integral
% a true model function m(x)that is zero at the boundaries
mtrue = sin(4*pi*x/xmax)+0.5*sin(7*pi*x/xmax)+0.25*sin(13*pi*x/xmax);
% derivative and integral operators
D = (1/Dx)*toeplitz( [1, -1, zeros(1,M-2)]', [1, zeros(1,M-1)] );
I = Dx*toeplitz( ones(M,1), [1, zeros(1,M-1)] );
% data is a linear operator L on the model parameters, d=Lm
% where the linear operator contains both a derivative and an intergal:
% d = 0.1*( dm/dx + 2 * (integral 0 x) m dm )
% brute force dE/dm, to check
% Eo = (dtrue-do)'*(dtrue-do);
% E1 = (dtrue-d1)'*(dtrue-d1);
% dEdmA = -2*(L'*(dtrue-do));
% fprintf('%d %f %f\n', i, dEdm(i), dEdmA(i) );
Ehist = zeros(Niter+1,1);
% error and its gradient at the trial solution
Ego = (dtrue-dgo)'*(dtrue-dgo);
dEdmo = -2*Dx*(L'*(dtrue-dgo));
% gradient descent method
v = -dEdmo / sqrt(dEdmo'*dEdmo);
Eg = (dtrue-dg)'*(dtrue-dg);
dEdm= -2*(L'*(dtrue-dg));
if( (Eg<=(Ego + c1*alpha*v'*dEdmo)) )
Dmg = sqrt( (mg-mgo)'*(mg-mgo) );
% plot one intermediate result
% terminate on small change
fprintf('%d iterations\n',k);
% plot the true model parameters
plot(x,mtrue,'k-','LineWidth',6);
plot(x,mgo1,'g:','LineWidth',3);
plot(x,mgo,'g-','LineWidth',3);
plot(x,mgo2,'g--','LineWidth',3);
fprintf('Caption: True solution (black), starting solution (green dotted line),\n');
Caption: True solution (black), starting solution (green dotted line),
fprintf('solution in progress (green dashed curve), estimated solution (green curve)\n');
solution in progress (green dashed curve), estimated solution (green curve)
plot(x,dtrue,'r-','LineWidth',3);
plot(x,dgo,'k-','LineWidth',1);
fprintf('Caption: True data (black curve), predicted data (red curve).\n');
Caption: True data (black curve), predicted data (red curve).
axis( [1, 101, -4, 0 ] );
plot([1:101],log10(Ehist(1:101)/Ehist(1)),'k-','LineWidth',3);
ylabel('log10 error, E(i)');
fprintf('Caption: Prediction error E Eas a function of iteration number i.\n');
Caption: Prediction error E Eas a function of iteration number i.
% simple 1D backprojection example
% the relationship d(x) = L m(x) is "inveted"
% using the back-propagation approximation m=L'd
% (where L' is the adjoint of L). In this
% example, L is the integral from 0 to x, Linv
% is the first derivative, and its adjoint is the
% integral from infinity to x
mtrue = sin(f1*x).*exp(-((x-xbar1).^2)/(2*sigma2)) + sin(f2*x).*exp(-((x-xbar2).^2)/(2*sigma2));
% data is integral of model
axis( [0, M, min(mtrue), max(mtrue) ] );
plot( x, mtrue, 'k-', 'LineWidth', 2 );
fprintf('Caption: True solution\n');
% exact solution is m = Linv d
% where Linv the derivative operator d/dx
axis( [0, M, min(mest1), max(mest1) ] );
plot( x(1:end-1), mest1, 'k-', 'LineWidth', 2 );
fprintf('Caption: Estimated solution using exact inversion.\n');
Caption: Estimated solution using exact inversion.
% solution via backprojection
mest2 = flipud(Dx*cumsum(flipud(dobs)));
% plot back-projected solution
axis( [0, M, min(mest2), max(mest2) ] );
plot( x, mest2, 'k-', 'LineWidth', 2 );
% title('m back-projected');
fprintf('Caption: Estimated solution using back-projection.\n');
Caption: Estimated solution using back-projection.
% backprojection in a a simple tomography example
% The data are "travel times", that is integrals
% of m(x,y) along straight line rays between sources
% and receivers on the edges of a square region.
clear gdaGsparse gdaepsilon;
global gdaGsparse gdaepsilon;
% true model m(x,y) is read from a file
Sbig=double(imread('../data/magma.tif','tif'));
[Nybig, Nxbig] = size(Sbig);
S=Sbig(1:idec:Nybig-1, 1:idec:Nxbig-1);
% rearrage model into vector
% order of model parameters in vector: index = (ix-1)*Ny + iy
rowindex = zeros(Nx*Ny,1);
colindex = zeros(Nx*Ny,1);
set(gca, 'YDir', 'reverse' );
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
axis( [xmin, xmax, ymin, ymax] );
image( [xmin, xmax], [ymin, ymax], S, 'CDataMapping', 'scaled');
fprintf('Caption: True model.\n');
% define sources and receivers on the edges of the box
xs(1:Nside,1)=xmin+dx/10;
xs(1:Nside,2)=ymin+Dy*[1:Nside]'/(Nside+1);
xs(Nside+1:2*Nside,1)=xmax-dx/10;
xs(Nside+1:2*Nside,2)=ymin+Dy*[1:Nside]'/(Nside+1);
xs(2*Nside+1:3*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(2*Nside+1:3*Nside,2)=ymin+dy/10;
xs(3*Nside+1:4*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(3*Nside+1:4*Nside,2)=ymax-dy/10;
R = sqrt( (xs(i,1)-xs(j,1))^2 + (xs(i,2)-xs(j,2))^2 );
if (R<Rmin) % exclude really short rays
% order of model parameters in vector: index = (ix-1)*Ny + iy
% use row-column-value method, instead of gdaGsparse=spalloc(N,M,2*N*Nx);
% also use build non-sparse row Gk first method
r = sqrt( (x2-x1)^2 + (y2-y1)^2 );
% I use a sloppy way of computing the length of the ray
% in each pixel. I subdivide the ray into Nr pieces, and
% assign each piece to exactly one pixel, the one its
if( 0 ) % slow but sure way
ix = 1+floor( (Nx-1)*(x-xmin)/Dx );
iy = 1+floor( (Ny-1)*(y-ymin)/Dy );
% gdaGsparse(k,q) = gdaGsparse(k,q) + dr;
else % faster way, or so we hope
% calculate the array indices of all the ray pieces
xv = x1 + (x2-x1)*[1:Nr]'/Nr;
yv = y1 + (y2-y1)*[1:Nr]'/Nr;
ixv = 1+floor( (Nx-1)*(xv-xmin)/Dx );
iyv = 1+floor( (Ny-1)*(yv-ymin)/Dy );
% now count of the ray segments in each pixel of the
% image, and use the count to increment the appropriate
% element of G. The use of the hist() function to do
% the counting is a bit weird, but it seems to work
icount = find( count~=0 );
% gdaGsparse(k,icount) = gdaGsparse(k,icount) + count(icount)*dr;
Gk(icount)=Gk(icount)+count(icount)'*dr;
% now accumulate Gk into element arrays
gdaGsparse = sparse(Gr,Gc,Gv,N,M);
% dont clear element arrays yet
% compute true travel times
% inversion of data by damped least squares
gdaepsilon = 1.0e-2*max(max(abs(gdaGsparse)));
mest1=bicg(@gda_dlsmul,GTd,1e-5,3*M);
bicg converged at iteration 40 to a solution with relative residual 9.5e-06.
% inversion of data by back projection
% normalize each row of the data kernel to unity
% with compensatory normalization of the data
rowsum = sum(gdaGsparse,2);
rowsumr(find(~isfinite(rowsumr)))=0;
% manipulaste row-col-value arrays instead of
% GP(:,i) = gdaGsparse(:,i) .* rowsumr;
Gv(el) = Gv(el) * rowsumr(k);
GP = sparse(Gr,Gc,Gv,N,M);
clear Gr Gc Gv; % now clear them
% rearrange m back into image
% plot estimated image (damped least squares)
set(gca, 'YDir', 'reverse' );
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
axis( [xmin, xmax, ymin, ymax] );
image( [xmin, xmax], [ymin, ymax], Sest1, 'CDataMapping', 'scaled');
% title( 'LS estimate' );
fprintf('Caption: Estimated image, using damped least squares\n');
Caption: Estimated image, using damped least squares
% plot estimated image (back projection)
set(gca, 'YDir', 'reverse' );
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
axis( [xmin, xmax, ymin, ymax] );
Sest2 = Sest2-mean(Sest2(:));
caxis( [zmin, zmax] ); % color scale fudged here
image( [xmin, xmax], [ymin, ymax], Sest2, 'CDataMapping', 'scaled');
% title( 'BP estimate' );
fprintf('Caption: Estimated image, using back projection\n');
Caption: Estimated image, using back projection
% temperature / chemical reaction example
% adjoint method used to build data kernel
% and the problem is then solved usign least squares
% kernel for differential equation
% differential equation (d/dt + c)u = f
% has green function G(t,t')=H(t-t')exp(-c(t-t'))
% create a heat function m
% allow two options, complex and simple
% complicated m built out of Gaussians
m = exp(-0.5*((t-t1).^2)/(sigma1^2));
m = m+2*exp(-0.5*((t-t1).^2)/(sigma1^2));
m = m+0.5*exp(-0.5*((t-t1).^2)/(sigma1^2));
m = m+exp(-0.5*((t-t1).^2)/(sigma1^2));
% very simple m for test purposes
% perform the Green function integral to predict u
% use convolution function; faster than coding the integral
u = Dt*conv(Ho,m); % Green function integral is a convolution
% data is P with random noise
dobs = dtrue+random('Normal',0,sigmad,N,1);
% plot Green function for t'=10;
H = (t>tp).*exp(-c*(t-tp));
plot( t, H, 'k-', 'LineWidth', 3 );
axis( [0, tmax, 0, 1.2*max(m)] );
plot( t, m, 'k-', 'LineWidth', 3 );
axis( [0, tmax, 0, max(u)] );
plot( t, u, 'k-', 'LineWidth', 3 );
axis( [0, tmax, 0, 1.1*max(dobs)] );
plot( t, dobs, 'k-', 'LineWidth', 3 );
fprintf('Caption: (Top row) Green function, (2nd row) true solution,\n');
Caption: (Top row) Green function, (2nd row) true solution,
fprintf('(3rd row) field u, (bottom row) true data.\n');
(3rd row) field u, (bottom row) true data.
% adjoint differential equation (-d/dt + c)u = h
% has green function G(t,t')=H(t'-t)exp(c(t-t'))
Gd(i,j)=Dt*(-b/c)*(exp(-c*(ti-tee))-1);
% solve with dampled least squares
mest = (Gd'*Gd + e2*eye(M,M))\(Gd'*dobs);
% plot estimsted solution
plot( t, mest, 'r--', 'LineWidth', 3 );
plot( t, m, 'k-', 'LineWidth', 2 );
% plot green function for t'=10;
H = Dt*(t<tp).*exp(c*(t-tp));
plot( t, H, 'k-', 'LineWidth', 3 );
ylabel('Q(t,tau=100)','FontSize',10);
% plot a few rowa of the data kernel
axis( [0, tmax, 0, 1.1*max(Gd(i,:)) ] );
plot( t, Gd(i,:)', 'k-', 'LineWidth', 3 );
ylabel('g(i=125,t)','FontSize',10);
axis( [0, tmax, 0, 1.1*max(Gd(i,:)) ] );
plot( t, Gd(i,:)', 'k-', 'LineWidth', 3 );
ylabel('g(i=250,t)','FontSize',10);
axis( [0, tmax, 0, 1.1*max(Gd(i,:)) ] );
plot( t, Gd(i,:)', 'k-', 'LineWidth', 3 );
ylabel('g(i=375,t)','FontSize',10);
fprintf('Caption: (top row) Adjoint Green function. (2nd-4th row)\n');
Caption: (top row) Adjoint Green function. (2nd-4th row)
fprintf('selected rows of the data kernel.\n');
selected rows of the data kernel.
fprintf('Caption: Data kernel G\n');
% Uses heat flow problem to illustrate the use of
% the adjoint method to compute the data kernel dd/dm.
% This deriavtive is calculated two ways, using finite
% difference and using the adjoint method (achieving the
% same result, of course).
% reference solution is analytic
u0 = (t>0) .* exp( -c0*t );
% plot the reference solution
axis( [t(1), t(end), -1.1, 1.1] );
% the data is the b times the integral of the solution
% predicted data associated with the reference field
d0 = b*Dt*cumsum( (t>0) .* u0 );
% point scatterer at t0, that is
% c(t) = c0 + dm * delta(t-t0)
t0 = 0.5; % time of scatterer
it0 = find(t>=t0,1); % index of scatterer
dm = 0.3; % amplitude of scatterer
% field perturbation according the the Born approximation
du = -dm * u0(it0) * (t>t0) .* exp( -c0*(t-t0) );
u = u0 + du; % total field
d = b*Dt*cumsum( (t>0) .* u ); % predicted data
% plot the reference and total field
plot( t, u, 'k-', 'LineWidth', 3 );
plot( t, u0, 'r-', 'LineWidth', 2 );
% plot the predicted data associated with u and u0
axis( [t(1), t(end), -1.1, 1.1] );
plot( t, d, 'k-', 'LineWidth', 3 );
plot( t, d0, 'r-', 'LineWidth', 2 );
% derivative by finite difference
axis( [t(1), t(end), -yr, yr] );
plot( t, dddm, 'k-', 'LineWidth', 3 );
% derivative by adjoint formula
% which should (and does) exactly match the adjoint formula
dddm2 = -(b/c0)*(t>0).*(t>t0).*(1-exp(-c0*(t-t0))).*exp(-c0*t0);
plot( t, dddm2, 'r-', 'LineWidth', 2 );
fprintf('Caption: (Top) Field u (black curve) and reference field u0 (red curve),\n');
Caption: (Top) Field u (black curve) and reference field u0 (red curve),
fprintf('(middle) data predicted by the field (black) and by the reference field (red),\n');
(middle) data predicted by the field (black) and by the reference field (red),
fprintf('(Bottom) derivative dd/dm by finite differences (black) and the adjoint method\n');
(Bottom) derivative dd/dm by finite differences (black) and the adjoint method
% data kernel for heat flow problem with time-variable thermal constant c(t)
ctrue = c0+c01*sin(pi*t/tmax).*(t>0.0);
% true solution and correponding data
utrue(i) = utrue(i-1) - ctrue(i)*Dt*utrue(i-1);
dtrue = b*Dt*cumsum(utrue);
% reference solution and corresponding data
u0(i) = u0(i-1) - c0(i)*Dt*u0(i-1);
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, utrue, 'k-', 'LineWidth', 3);
plot( t, u0, 'r-', 'LineWidth', 2);
fprintf('Caption: True field u (black curve) and reference field (red curve)\n');
Caption: True field u (black curve) and reference field (red curve)
axis( [t(1), t(N), -1.1, 1.1] );
plot( t, dtrue, 'k-', 'LineWidth', 3);
plot( t, d0, 'r-', 'LineWidth', 2);
fprintf('Caption: True data d (black curve) and data predicted by the reference field (red curve)\n');
Caption: True data d (black curve) and data predicted by the reference field (red curve)
for i = (1:N) % loop over data points
% ( -d/dt + c0 ) lami = hi where hi = b H(ti-t)
hi = b*(t(i)*ones(N,1)>=t);
% descrete version of adjoint equation
% (-lamj(k+1)+lamj(k))/Dt + c0 lamj(k+1) = hi
% lamj(k) = lamj(k+1) - c0 = Dt lamj(k+1) + Dt hi
lami(k) = lami(k+1) - c0(k+1)*Dt*lami(k+1) + Dt*hi(k+1);
% perturbed equation is L u = 0 with
% L = ( d/dt + c0 + mj delta(t-tj) )
% ddidmj = ( lami, dLdmi u0 ) = -( lami, delta(t-tj) u0 ) =
% ddidmj = - lami(j) u0(j)
ddidmj(i,1:N) = - ( lami .* u0 )';
G = ddidmj; % data kernel
fprintf('Caption: Data kernel\n');
% solve inverse problem for the function c(t) in
% the cooling problem using adjoint methods
% the equation of cooling is
% ( d/dt + c(t) ) u = 0 with initial condition u(0)=1
% the true c(t) is a constant plus a sinusoid
% the problem is solved iteratively with Newton's
% methid using the linearized data kernel Gij = ddi/dmj
% calculated by adjoint methods.
% both the original equation and the adjoint equation
% are solved numerically with a very simple recursion based on
% Euler's method. A more accurate method, such as Runga
% Kutta could be substituted.
% At each iteration, the function c(t) is parameteried
% c = c0 + (sum over i) mi delta(t-ti)
% which implies that the update is c0 = c0 + m/Dt
% where the factor of Dt arises from the delta funcion
% that is, the integral of the mi delta(t-ti) over one
% pixel is mi, and the integral of mi/Dt over one pixel
% syntheic data are created by adding a small amount of
% random noise to the true data.
% The problem solved via generalized least squares, with
% both flatness and smallness prior information. A good
% solution is achieved in about five iterations.
ctrue = c0+c01*sin(pi*t/tmax).*(t>0.0);
% true solution and correponding data
utrue(i) = utrue(i-1) - ctrue(i)*Dt*utrue(i-1);
dtrue = b*Dt*cumsum(utrue);
dobs = dtrue + random('Normal',0.0,sd,N,1);
% flatness matrix for generalized least squares
r = [-1.0; zeros(N-2,1)];
c = [-1.0, 1.0, zeros(1,N-2) ];
h1 = zeros(N-1,1); % right hand side
% reference solution and corresponding data
u0(i) = u0(i-1) - c0(i)*Dt*u0(i-1);
for i = (1:N) % loop over data points
% ( -d/dt + c0 ) lami = hi where hi = b H(ti-t)
hi = b*(t(i)*ones(N,1)>=t);
% descrete version of adjoint equation
% (-lamj(k+1)+lamj(k))/Dt + c0 lamj(k+1) = hi
% lamj(k) = lamj(k+1) - c0 = Dt lamj(k+1) + Dt hi
lami(k) = lami(k+1) - c0(k+1)*Dt*lami(k+1) + Dt*hi(k+1);
% perturbed equation is L u = 0 with
% L = ( d/dt + c0 + mj delta(t-tj) )
% ddidmj = ( lami, dLdmi u0 ) = -( lami, delta(t-tj) u0 ) =
% ddidmj = - lami(j) u0(j)
ddidmj(i,1:N) = - ( lami .* u0 )';
G = ddidmj; % data kernel
% H Dm = h1 - H c0 or Hm = h1
Dh = h1 - H*c0; % apply prior info to solution
Dh = h1; % apply to perturbation, only
% generalized least squares with prior informaiton of flatness and smallness
% flatness can be applied to either solution or perturbation
% smallness is applied to perturbation
F = [G; epsi1*H; epsi2*eye(N) ];
f = [Dd; epsi1*Dh; epsi2*h0 ];
c0 = c0 + m/Dt; % factor of Dt from delta function
axis( [t(1), t(N), 0.0, 1.5] );
plot( t, ctrue, 'k-', 'LineWidth', 3);
plot( t, c0, 'r-', 'LineWidth', 2);
plot( t, cstart, 'g--', 'LineWidth', 2);
fprintf('Caprion: True model parameter c (black curve), reference value\n');
Caprion: True model parameter c (black curve), reference value
fprintf('(green curve) and estimated (red curve)\n');
(green curve) and estimated (red curve)
% Data Assimilation example
% A insulated rod of length Lx cools because its ends
% are held at a temperature u(0)=u(Lx)=0. Heat moves
% according to the thermal diffusion equation
% second derivative matrix
r = [-2.0; 1.0; zeros(Lx-2,1) ];
c = [-2.0, 1.0, zeros(1,Lx-2) ];
D2 = toeplitz( r,c) / (Dx^2);
D = Dt*kappa*D2 + eye(Lx);
m0 = exp( -0.5 * ((x-xmid).^2) / sx2 );
% true solution, by stepping the equation forward in time
utrue(1:Lx,i) = D*utrue(1:Lx,i-1);
% data positions and time are randomly assigned
Ndt = 35; % data per time
N = (Lt-1)*Ndt; % total number of data
% no data at zero time and no data at boundaries
% index array kd gives locations of data
kd(1:Ndt,i-1) = sort(j(1:Ndt));
% although noise is defined everywhere, only values at
% data positions-times are used
nse = random('Normal',0.0,sd,Lx,Lt);
sx = 2.0; % half-width of Gaussian source
mg = exp( -0.5 * ((x-xmid).^2) / sx2 );
% guessed solution, treat source as initial condition
% and use recurion to propagate solution forward in time
u(1:Lx,i) = D * u(1:Lx,i-1);
% error = observed temperature - predicted temperature
% observed temperature = true temperature + noise
e = zeros(Lx,Lt); % although the error is defined
% for all poitions-times, it is non-zero only at observation points
e(k,i) = (utrue(k,i)+nse(k,i))-u(k,i);
% simplified gradent method
if( (iter>1) && (emax>emaxlast) )
% adjoint field, step backward in time, source is error
lam(1:Lx,i-1) = D*lam(1:Lx,i) + Dt*e(1:Lx,i);
% dE/dmi is adjoint field at t=0
% change in -dEdm direction
axis( [0, maxiter, 0, 1.1*max(Emax) ] );
plot([0:maxiter-1]',Emax,'k-','LineWidth',2);
fprintf('Caption: Error E as a function of iteration number i\n');
Caption: Error E as a function of iteration number i
% plot images of solution and etc.
axis( [0, tmax, 0, xmax] );
axis( [0, tmax, 0, xmax] );
axis( [0, tmax, 0, xmax] );
e0absmax = max(abs(e0(:)));
caxis( [-e0absmax, e0absmax ]);
axis( [0, tmax, 0, xmax] );
lam0absmax = max(abs(lam0(:)));
caxis( [-lam0absmax, lam0absmax] );
axis( [0, tmax, 0, xmax] );
axis( [0, tmax, 0, xmax] );
axis( [0, tmax, 0, xmax] );
caxis( [-e0absmax, e0absmax ]);
axis( [0, tmax, 0, xmax] );
caxis( [-lam0absmax, lam0absmax ]);
fprintf('Caption: Data assimilation problem. (Top row) (left) True field u,\n');
Caption: Data assimilation problem. (Top row) (left) True field u,
fprintf('(2nd) initial field u0, (3rd) initial error e0 (right) initial adjoint\n');
(2nd) initial field u0, (3rd) initial error e0 (right) initial adjoint
fprintf('field lambda0. (Bottom) (left) True field u, (2nd) final predicted u,\n');
field lambda0. (Bottom) (left) True field u, (2nd) final predicted u,
fprintf('(3rd) final error e, (right) final adjoint field lambda\n');
(3rd) final error e, (right) final adjoint field lambda
% In order for this script to actually write the datafile
% you must set dowrite to 1
dm = a1*exp(-((t-t01).^2)/(2*s1*s1))+a2*exp(-((t-t02).^2)/(2*s2*s2));
axis( [t(1), t(end), -mr, mr] );
plot( t, dm, 'k-', 'LineWidth', 3 );
u0 = (t>0) .* exp( -c0*t );
axis( [t(1), t(end), -ur, ur] );
plot( t, u0, 'k-', 'LineWidth', 3 );
% sum up point scatterers to get perturbation in solution
dui = -dm(i) * u0(i) * (t>t0i) .* exp( -c0*(t-t0i) );
plot( t, u, 'r-', 'LineWidth', 2 );
% data predicted by reference solution
d0 = b*Dt*cumsum( (t>0) .* u0 );
% data predicted by heterogenous solution
d = b*Dt*cumsum( (t>0) .* u );
axis( [t(1), t(end), -dr, dr] );
plot( t, d0, 'k-', 'LineWidth', 3 );
plot( t, d, 'r-', 'LineWidth', 3 );
fprintf('Caption: (Top) Heterogeneities in the model parsmeters. (Middle) True field u (red curve) and\n');
Caption: (Top) Heterogeneities in the model parsmeters. (Middle) True field u (red curve) and
fprintf('reference field u0 (black curve), (bottom) true data (red curve) and reference data (black curve)\n');
reference field u0 (black curve), (bottom) true data (red curve) and reference data (black curve)
dlmwrite('gda1106_data.txt', [t,dm,d], 'delimiter', '\t' );
dlmwrite('../data/gda1106_data.txt', [t,dm,d], 'delimiter', '\t' );