L=TV_operator_3D( msh ) construct Total Variation operator for a 3D mesh INPUTS: msh = a 3D scaip msh structure with .msh.vtx_c_c, .elem_c defined OUTPUTS: L = TV operator (generally used for regularisation) Copyright 2004 Andrea Borsic, SC-AIP s.a.s. Scientific Computing & Applied Inverse Problems www.sc-aip.com Modifications (C) 2005 Andy Adler. License: GPL version 2 or version 3 $Id: TV_operator_3D.m 3039 2012-06-06 14:58:10Z bgrychtol $
0001 function L=TV_operator_3D( msh ) 0002 % L=TV_operator_3D( msh ) 0003 % construct Total Variation operator for a 3D mesh 0004 % 0005 % INPUTS: 0006 % msh = a 3D scaip msh structure with .msh.vtx_c_c, .elem_c defined 0007 % OUTPUTS: 0008 % L = TV operator (generally used for regularisation) 0009 % 0010 % Copyright 2004 Andrea Borsic, SC-AIP s.a.s. 0011 % Scientific Computing & Applied Inverse Problems www.sc-aip.com 0012 % Modifications (C) 2005 Andy Adler. 0013 % License: GPL version 2 or version 3 0014 % $Id: TV_operator_3D.m 3039 2012-06-06 14:58:10Z bgrychtol $ 0015 0016 %global SP3D % Sparse 3D matrix used in the computations 0017 %SP3D=[]; 0018 0019 num_vtx=size(msh.vtx_c,1); 0020 num_tet=size(msh.elem_c,1); 0021 0022 list=[]; % List is an auxiliary variable which will hold for each row the two facing thetrahydra and the shared face area 0023 0024 %We build a selection array, to index the T matrix 0025 0026 SEL=[ 1 1 1 2; 2 3 4 4; 3 4 2 3]; 0027 0028 % allocate_3Dsparse(num_vtx,num_vtx,num_vtx,num_tet*4); 0029 SP3D=spalloc(num_vtx,num_vtx^2,num_tet*4); 0030 0031 for k=1:num_tet 0032 0033 for j=1:4 % cycle on the faces of each thetrahydra 0034 0035 face_vtx(1)=msh.elem_c(k,SEL(1,j)); % face_vtx are the varticies on a face of thetrahydra k 0036 face_vtx(2)=msh.elem_c(k,SEL(2,j)); 0037 face_vtx(3)=msh.elem_c(k,SEL(3,j)); 0038 0039 face_vtx=sort(face_vtx); % faces must be unique 0040 0041 % simplex=read_from_3Dsparse(face_vtx,num_vtx,num_vtx,num_vtx); 0042 %function val=read_from_3Dsparse(mnp,M,N,P); 0043 % 0044 %global SP3D 0045 %val=SP3D(mnp(1),(mnp(2)-1)*N+mnp(3)); 0046 simplex=SP3D(face_vtx(1),(face_vtx(2)-1)*num_vtx+face_vtx(3)); 0047 0048 if (simplex==0) 0049 0050 % write_to_3Dsparse(mnp,M,N,P,val); 0051 % SP3D(mnp(1),(mnp(2)-1)*N+mnp(3))=val; 0052 % 0053 % write_to_3Dsparse(face_vtx,num_vtx,num_vtx,num_vtx,k); 0054 SP3D(face_vtx(1),(face_vtx(2)-1)*num_vtx+face_vtx(3))=k; 0055 0056 else 0057 0058 vec1=msh.vtx_c(face_vtx(1),:)-msh.vtx_c(face_vtx(2),:); % vec1,vec2 are vectors along edgs, used for the area calculation as cross prod. 0059 vec2=msh.vtx_c(face_vtx(1),:)-msh.vtx_c(face_vtx(3),:); 0060 facearea=0.5*norm(cross(vec1,vec2)); 0061 list=[list;[k,simplex,facearea]]; % Triangles and length of the shared edge are written into the list 0062 0063 end % if then else 0064 0065 end % for j 0066 0067 end % for k 0068 0069 L=spalloc(length(list),num_tet,2*length(list)); 0070 0071 for i=1:length(list) 0072 0073 L(i,list(i,1))=list(i,3); 0074 L(i,list(i,2))=-list(i,3); 0075 0076 end % for 0077 0078 %clear SP3D; 0079 0080 0081 %%%%%%%%%% Auxiliary functions for handling the 3D sparse matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0082 0083 %%%%%%%%%% mnp is a vector of 3 indexes into the matrix, and M,N,P the matrix size along the three dimensions 0084 % - removed these functions for inclusion with EIDORS -aa Dec05 0085 % 0086 %function allocate_3Dsparse(M,N,P,max_elements); 0087 % 0088 %global SP3D 0089 %SP3D=spalloc(M,N*P,max_elements); 0090 % 0091 %function write_to_3Dsparse(mnp,M,N,P,val); 0092 % 0093 %global SP3D 0094 %SP3D(mnp(1),(mnp(2)-1)*N+mnp(3))=val; 0095 % 0096 %function val=read_from_3Dsparse(mnp,M,N,P); 0097 % 0098 %global SP3D 0099 %val=SP3D(mnp(1),(mnp(2)-1)*N+mnp(3)); 0100 %