Merge pull request #8 from jessdtate/master

Adding matlab inverse examples

MatlabLibrary/README.txt

@@ -0,0 +1,14 @@
+README
+======
+
+This directory contains a library of Matlab functions developed
+at the SCI institute and the CIBC, which are provided to extend the basic Matlab
+functions with code which helps in Biological Simulation Applications
+and with the integration into SCIRun.
+
+The functions in this library are organized in directories which correspond
+to different projects being pursued at the SCI institute and the CIBC.
+
+Add this directory and all the subdirectories to your matlab path to use these 
+functions in SCIRun and the Fwd/Inv Toolkit. 
+

MatlabLibrary/Total_Variation/BuildLocMassMat_LinearTet.m

@@ -0,0 +1,53 @@
+function [ LocMassMat, jcbInt ] = BuildLocMassMat_LinearTet( rhsType,vertX, vertY, vertZ, ws, phi )
+% This function build the local mass matrix system. Only works for linear tet element.
+% It hard codes the calculation of derivatives, quadratures are used only for integration
+%Copyright: Dafang Wang, SCI Institute, 2008-01-30
+
+%  rhsType: rright hand side type, 0 for laplace
+%  elmtOrder  the order of the base test polynomials,1 for linear, 2 for quadratic
+%  vertX: x-coordinates for 4 vertices of the tet
+%  vertY, vertZ: same as vertX
+%  zp, ws, the quadrature points and weights in x,y,z dimensions respectively. 
+% 
+%  IMPORTANT: 
+% Basis func phi are values defined over quaduatures. It doesn't matter whether '-1to1' or '0to1' space. 
+% The integration is done based on the '-1to1' tet. The weights have been set in that way.
+%
+%Output:
+
+ w_x = ws(:,1);
+  w_y = ws(:,2);
+w_z = ws(:,3);
+localMatDim = 4;
+
+    %calculate the jacobi from [-1,1] std tet to physical tet elements
+    [ jcbInt, jcbmatInt ] = JcbStdTet2GeneralTet( vertX, vertY, vertZ, '-1to1' );
+
+    LocMassMat = zeros( localMatDim, localMatDim);
+
+    for p = 1:localMatDim
+        for q = p:localMatDim
+            integrand = phi(:,:,:,p) .* phi(:,:,:,q);
+             LocMassMat(p,q) = IntegrationInTet(integrand, w_x, w_y, w_z, jcbInt );
+             LocMassMat(q,p) = LocMassMat(p,q);
+        end
+    end
+return;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% calculate the integration in the tet, follow the section 4.1.1.2 in spencer's book
+% jcb, the jacobian from [-1,1] tet to physical tet, wx,wy,wz are weights
+% derived from [-1,1] and have been adapted for tet integration
+function r = IntegrationInTet(f, wx, wy, wz, jcb )
+    nx = length(wx);  ny = length(wy);  nz = length(wz);
+    for i = 1:nz
+        f(:,:, i ) = f(:,:,i) * wz(i);
+    end
+    for i = 1:ny
+        f(:,i,:) = f(:,i,:) * wy(i);
+    end
+    for i = 1:1:nx
+        f(i,:,:) = f(i,:,:) * wx( i );
+    end
+    r = sum(sum(sum(f ) ) ) * abs(jcb );
+return;

MatlabLibrary/Total_Variation/BuildLocStiffMat_LinearTet.m

@@ -0,0 +1,80 @@
+function [ LocStMat, LocDerivMat, vol, jcbInt, localRHS ] = BuildLocStiffMat_LinearTet( rhsFlag, ...
+    conduct, vertX, vertY, vertZ, varargin)
+% This function build the local stiffness matrix system
+%   1. It hard codes the calculation of derivatives
+%   2. The integration in the stiffness matrix is hard coded, as the derivatives are all zero.
+%   3. If right-hand size source is needed, quadratures are used for its integration.
+%   4. It differs from BuildLocStMatTet() in that it doesn't use quadrature for integration
+%Copyright: Dafang Wang, SCI Institute, 2011-10-30
+
+%  rhsType: right hand side type, 0 for laplace
+%  rhsFlag:  1 if compute rhs
+%  conduct:  the conductivity within this element
+%  vertX: x-coordinates for 4 vertices of the tet
+%  vertY, vertZ: same as vertX
+% coefmatBase: each column is the coefficient of the basis in [0,1] canonical basis
+% zp, ws, the quadrature points and weights in x,y,z dimensions respectively. 
+%    Note that quadratures points are defined within a [-1,1] cube, but the
+%    weights have been adapted for integration in tetrahedra
+
+% Output:
+%   - vol: element volume
+
+localMatDim = 4;
+LocStMat = zeros(4,4);        LocDerivMat = zeros(3,4);
+
+%Calculate the coefficients of basis functions in the current tet, in physical space
+coefmatBase = CalPhyBasisTet( 1, [vertX, vertY, vertZ] );      
+for i = 1:4
+    LocDerivMat(:, i ) = coefmatBase( 2:4, i ); %each column contains the dx, dy, dz for one basis func
+end
+%grad(u) = localDerivMat * [u1;u2;u3; u4]; the field value at 4 vertices
+
+%calculate the jacobi from [-1,1] std tet to physical tet elements
+[ jcbInt, jcbmatInt ] = JcbStdTet2GeneralTet( vertX, vertY, vertZ, '-1to1' );
+vol = ComputeTetVolume( [vertX, vertY, vertZ] );
+
+LocStMat = zeros( localMatDim, localMatDim);
+for i =1:4
+    for j = 1: i
+        k = (conduct * LocDerivMat(:, j) )' * LocDerivMat(:, i)  * vol;
+        LocStMat( i, j ) = k;
+        LocStMat( j, i ) = k;
+    end
+end
+
+if rhsFlag == 1 %need to compute rhs
+    if ~exist('varargin')    disp('BuildLocStiffMat_LinearTet.m:  error \n');  return; end
+    rhsType = varargin{1};
+    if rhsType == 0
+        localRHS = zeros( localMatDim, 1);
+    else %use quadratuves
+        zp = varargin{2};        ws = varargin{3};
+        z_x = zp(:,1);      w_x = ws(:,1);
+        z_y = zp(:,2);      w_y = ws(:,2);
+        z_z = zp(:,3);      w_z = ws(:,3);
+        qdrNumX = length(z_x);       qdrNumY = length(z_y);  qdrNumZ = length(z_z);
+        elmtX = zeros(qdrNumX, qdrNumY, qdrNumZ);   
+        [elmtX, elmtY, elmtZ] = TransformStdQdr2PhySpace(z_x, z_y, z_z, vertX, vertY, vertZ, 'Tet');
+        vecXYZ = cat(4, elmtX, elmtY, elmtZ); %the last dimension is the xyz value
+
+        [phi, phidX, phidY, phidZ] = EvalBasisTet(elmtOrder, coefmatBase, vecXYZ );
+
+        %calculate the rhs source term
+        eF = Test_Right_Side_Function( 3, rhsType, conduct, elmtX, elmtY, elmtZ );
+        for p = 1:localMatDim
+            rhsTerm = phi(:,:,:,p).* eF;
+            localRHS( p,1 ) = IntegrationInTet(rhsTerm, w_x, w_y, w_z, jcbInt );
+        end
+    end
+end
+
+
+return;
+end
+
+function [vol] = ComputeTetVolume( vertices )
+
+	vol = elemVolumeCalculation([1 2 3 4], vertices);
+
+end

MatlabLibrary/Total_Variation/Build_FEmat4TotalVar.m

@@ -0,0 +1,81 @@
+%Build FE-based matrices to be used for total variation
+clear;
+filename = '../ts_ht1_mesh1_Tr1.mat';
+FileOut = 'tvFE_mesh1_Tr1.mat';
+load(filename);
+
+feDerivMat = zeros( eleNumH, 3, 4); % partial deriv in each element
+feStiffMat = zeros( eleNumH, 4, 4); %stiffness matrix in each cell
+
+% stH = sparse( ndNumH, ndNumH);
+conduct = eye(3,3);
+eleVol_H = zeros( eleNumH,1);
+
+%%%%%%%%%%%%% feDerivMat and feStiffMat
+
+%-----------------Build the stiffness matrix. -----------------
+for i = 1:eleNumH
+    ch = elmtH( i, 1:4);
+    vx = ndMatH( ch, 1);   vy = ndMatH(ch, 2);    vz = ndMatH(ch, 3);
+
+    [ LocStMat, LocDerivMat, vol ] = BuildLocStiffMat_LinearTet( 0, conduct, vx, vy, vz);
+    eleVol_H( i ) = vol;
+    feDerivMat( i, :, :) = LocDerivMat;
+    feStiffMat( i, :, :) = LocStMat;
+%     stH( ch, ch) = stH(ch, ch) + LocStMat;
+end
+% full( MatrixDiff( stH, stMatH_std)) %validation
+
+%-----------------Build the mass matrix. -----------------
+[z_x, w_x] = JacobiGLZW( 4, 0, 0);
+[z_y, w_y] = JacobiGRZW( 4, 1, 0 );     w_y = w_y / 2;
+[z_z, w_z] =  JacobiGRZW( 4, 2, 0 );     w_z = w_z / 4;
+zpTet = [ z_x, z_y, z_z];          wsTet = [w_x, w_y, w_z];
+qdTet = Transform2StdTetSpace( z_x, z_y, z_z, '0to1');
+coefmatBaseTet = LocalCanonicalTetBasis('Base', 1);
+phiTet = EvalBasisTet(1, coefmatBaseTet, qdTet);
+
+gMassMat = sparse( ndNum, ndNum);
+for i = 1:eleNum
+    cn = Elmt( i, :);
+    vx = NdMat(cn, 1);  vy = NdMat(cn,2);       vz = NdMat(cn, 3);
+     [ locMassMat, jcbInt ] = BuildLocMassMat_LinearTet( 0,vx, vy, vz, wsTet, phiTet );
+     gMassMat( cn, cn) = gMassMat(cn, cn) + locMassMat;
+end
+
+%-----------------Build the FE matrix for dual -----------------
+[z_x, w_x] = JacobiGLZW( 4, 0, 0);
+[z_y, w_y] = JacobiGRZW( 4, 1, 0 );     w_y = w_y / 2;
+[z_z, w_z] =  JacobiGRZW( 4, 2, 0 );     w_z = w_z / 4;
+zpTet = [ z_x, z_y, z_z];          wsTet = [w_x, w_y, w_z];
+qdTet = Transform2StdTetSpace( z_x, z_y, z_z, '0to1');
+coefmatBaseTet = LocalCanonicalTetBasis('Base', 1);
+phiTet = EvalBasisTet(1, coefmatBaseTet, qdTet);
+
+Ax = sparse( ndNumH, ndNumH);       Ay = Ax;        Az = Ax;
+for i = 1:eleNumH
+    ch = elmtH( i, 1:4);
+    vx = ndMatH( ch, 1);   vy = ndMatH(ch, 2);    vz = ndMatH(ch, 3);
+    locDerivMat = squeeze(feDerivMat( i, :, :));
+    [ locAx, locAy, locAz ] = BuildLocCrossDerivMat_LinearTet( vx, vy, vz, wsTet, phiTet, locDerivMat );
+
+    Ax( ch, ch) = Ax( ch, ch) + locAx;
+    Ay( ch, ch) = Ay( ch, ch) + locAy;
+    Az( ch, ch) = Az( ch, ch) + locAz;
+end
+
+massMatT = full(gMassMat( ndv_TsSurf, ndv_TsSurf));
+M1 = mR' * (gStMat \ mQ'); %M1 is full matrix
+M1A = M1 \ [ Ax, Ay, Az];     
+
+%-----------------------  End -------------------------------------------------
+clear i z_* w_* cn ch Loc*Mat loc*Mat vx vy vz opt jcb* locA* vol
+save( FileOut );
+
+% %Compute true TV value
+% feGradMag = zeros( eleNumH, 1);  %gradient magnitude
+% feGradMag0 = ComputeGradMagPerTet( elmtH, feDerivMat, vTMP0); %true grad mag
+% feTV0 = dot( eleVol_H, feGradMag0);
+
+
+

MatlabLibrary/Total_Variation/CalPhyBasisTet.m

@@ -0,0 +1,84 @@
+function coefmat = CalPhyBasisTet( order, vts )
+% Each column of coefmat is the coefficients of the basis function defined
+% in a prism in  physical domain. The vertices of the prism is given by vts 
+
+% Input:  each row of vts gives the x/y/z coordinates of a vertice, the
+% ordering of vertices must follow the right-hand rule:
+% cross((v2-v1),(v3-v1)) = v4-v1
+
+    switch order 
+        case 1
+            sz = 4;
+            Lb = zeros(sz, sz);
+                for i =1:sz
+                    x = vts(i,1);
+                    y = vts(i,2);
+                    z = vts(i,3);
+                    evec = [1, x, y, z];
+                    Lb( i, :) = evec;
+                end
+             rhs = eye(sz);
+             coefmat = Lb \ rhs;
+
+        case 2
+                sz = 10;
+                Lb = zeros(sz, sz);
+                Onode = zeros( sz, 3); %all vertices
+                Onode(1:4,:) = vts;
+                Onode( 5, :) = mean( vts([1,2], :) );
+                Onode( 6, :) = mean( vts([1,3], :) );
+                Onode( 7, :) = mean( vts([1,4], :) );
+                Onode( 8, :) = mean( vts([2,3], :) );
+                Onode( 9, :) = mean( vts([2,4], :) );
+                Onode( 10, :) = mean( vts([3,4], :) );
+                 for s = 1:sz
+                    x = Onode(s,1);
+                    y = Onode(s,2);
+                    z = Onode(s,3);
+                    evec = [ 1, x, y, z, x*x, y*y, z*z, x*y, x*z, y*z];
+                    Lb( s, :) = evec;
+                 end
+                rhs = eye(sz);
+                coefmat = Lb \ rhs;   
+
+        case 3
+%             disp('CalPhyBasisTet(): cubic case needed to be completed');
+                sz = 20;
+                Lb = zeros(sz, sz);
+                Onode = zeros( sz, 3); %all vertices
+                %node:
+                Onode(1:4,:) = vts;
+                %edge:
+                Onode( 5, :) = (2*vts(1,:) + vts(2,:) ) / 3; %x^2
+                Onode( 6, :) = (vts(1,:) + 2*vts(2,:) ) / 3; %x^3
+                Onode( 7, :) =  (2*vts(1,:) + vts(3,:) ) / 3; %y^2
+                Onode( 8, :) =  (vts(1,:) + 2*vts(3,:) ) / 3; %y^3
+                Onode( 9, :) =  (2*vts(1,:) + vts(4,:) ) / 3; %z^2
+                Onode( 10, :) =  (vts(1,:) + 2*vts(4,:) ) / 3; %z^3;
+
+                Onode( 11, :) = (2*vts(2,:) + vts(3,:) ) /  3; %x^2 * y
+                Onode( 12, :) = (vts(2,:) + 2*vts(3,:) ) /3; %x * y^2
+                Onode( 13, :) =  (2*vts(2,:) + vts(4,:) )/3; %x^2 * z
+                Onode( 14, :) =  (vts(2,:) + 2*vts(4,:) )/3; %x* z^2
+                Onode( 15, :) =  (2*vts(3,:) + vts(4,:) ) /3; %y^2 * z
+                Onode( 16, :) =  (vts(3,:) + 2*vts(4,:) ) /3; %y* z^2;
+
+                %face:
+                Onode(17, :) = mean( [Onode(1,:); Onode(2,:); Onode(3,:)] ); %xy
+                Onode(18, :) = mean( [Onode(1,:); Onode(2,:); Onode(4,:)] );  %xz
+                Onode(19, :) = mean( [Onode(1,:); Onode(3,:); Onode(4,:)] );  %yz
+                Onode(20, :) = mean( [Onode(2,:); Onode(3,:); Onode(4,:)] ); %xYz
+
+                 for s = 1:sz
+                    x = Onode(s,1);
+                    y = Onode(s,2);
+                    z = Onode(s,3);
+                    xx = x*x;   yy = y*y;  zz = z*z;
+                    xy = x*y; yz=y*z; xz = x*z;
+                    ev = [ 1, x, y, z, xx, xx*x, yy, yy*y,zz, zz*z, xx*y, x*yy, xx*z, x*zz, yy*z, y*zz, xy, xz, yz, x*y*z ];
+                    Lb( s, :) = ev;
+                 end
+                rhs = eye(sz);
+                coefmat = Lb \ rhs;   
+    end
+return

MatlabLibrary/Total_Variation/ComputeGradMagPerTet.m

@@ -0,0 +1,14 @@
+function feGradMag = ComputeGradMagPerTet( elmtH, feDerivMat, v)
+%compute gradient magnitude for each linear tetrahedral element
+% Input:
+%   - v: the field data located at nodes
+%   - feDerivMat:  a Nele*3*4 matrix, each 3*4 matric is the partial deriv for the 4 local basis functions in each ele
+%       in 2D triangles, each tri has an 2*3 matrix for the partial deriv
+eleNumH = size( elmtH, 1);
+feGradMag = zeros( eleNumH, 1);
+for i = 1:eleNumH
+    cv = v(  elmtH(i,:) );
+    feGradMag( i ) = norm( squeeze( feDerivMat(i,:,:) ) * cv, 2);
+end
+
+end