computeMixedVolumeTensor

PURPOSE ^

calculate the coefficients for measure calculation in S-representation

SYNOPSIS ^

function tensor = computeMixedVolumeTensor(hMatrix, handleH, dim,incidenceMatrix, verbose)

DESCRIPTION ^

 calculate the coefficients for measure calculation in S-representation

 IMPORTANT: this function uses the sum of hC-vectors instead of the
   Minkoski sum. This is valid as long as the tensor of mixed volume is
   applied to unified partitions. It CANNOT be applied to arbitrary
   combinations of structuring elements.

 THIS FUNCTION IS NO USER FUNCTION.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function tensor = computeMixedVolumeTensor(hMatrix, handleH, dim, ...
0002     incidenceMatrix, verbose)
0003 % calculate the coefficients for measure calculation in S-representation
0004 %
0005 % IMPORTANT: this function uses the sum of hC-vectors instead of the
0006 %   Minkoski sum. This is valid as long as the tensor of mixed volume is
0007 %   applied to unified partitions. It CANNOT be applied to arbitrary
0008 %   combinations of structuring elements.
0009 %
0010 % THIS FUNCTION IS NO USER FUNCTION.
0011 
0012 % The elk-library: convex geometry applied to crystallization modeling.
0013 %   Copyright (C) 2012 Alexander Reinhold
0014 %
0015 % This program is free software: you can redistribute it and/or modify it
0016 %   under the terms of the GNU General Public License as published by the
0017 %   Free Software Foundation, either version 3 of the License, or (at your
0018 %   option) any later version.
0019 %
0020 % This program is distributed in the hope that it will be useful, but
0021 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0022 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0023 %   General Public License for more details.
0024 %
0025 % You should have received a copy of the GNU General Public License along
0026 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0027 if ~exist('verbose', 'var') || isempty(verbose)
0028     verbose = 0;
0029 end
0030 
0031 if verbose
0032     disp('Computing mixed volumes..');
0033     fprintf('......');
0034 end
0035 switch dim
0036     case 1
0037         tensor = zeros([1 size(hMatrix, 1)]);
0038         for iSe = 1:size(hMatrix, 1)
0039             tensor(iSe) = handleH(hMatrix(iSe, :)');
0040         end
0041         
0042     case 2
0043         tensor = nan([size(hMatrix, 1) size(hMatrix, 1)]);
0044         for iSe = 1:size(hMatrix, 1)
0045             tensor(iSe,iSe) = handleH(hMatrix(iSe, :)');
0046         end
0047         % Measure for all combinations of two polytopes
0048         % V(lA*A + lB*B) = lA^2*V(A) + lB^2*V(B) + 2*lA*lB V(A, B)
0049         %
0050         % V(A,B) = (V_AB - V(A) - V(B)) / 2
0051         % but as V(A,B) is applied it would be multiplied by 2, hence, we
0052         % omit this.
0053         for iSe = 1:size(hMatrix, 1)
0054             for jSe = (iSe+1):size(hMatrix, 1)
0055                 if incidenceMatrix(iSe, jSe) == 0
0056                     continue
0057                 end
0058                 tensor(iSe,jSe) = (handleH(...
0059                                      hMatrix(iSe, :)' + hMatrix(jSe, :)')...
0060                                    - tensor(iSe,iSe) ...
0061                                    - tensor(jSe,jSe));
0062             end
0063         end
0064     case 3
0065         % initialise 3D matrix
0066         tensor = nan([size(hMatrix, 1), ...
0067                         size(hMatrix, 1), ...
0068                         size(hMatrix, 1)]);
0069         
0070         % Measure for all structuring elements, which can be calculated
0071         % directly.
0072         for iSe = 1:size(hMatrix, 1)
0073             if verbose
0074                 fprintf('\b\b\b\b\b\b%5.1f%%', ...
0075                     (iSe/size(hMatrix, 1))*100);
0076             end
0077             tensor(iSe, iSe, iSe) = handleH(hMatrix(iSe, :)');
0078         end
0079         
0080         % Measure for all combinations of two polytopes
0081         % V(lA A + lB B) = lA^3 * V(A) + 3*lA^2*lB V(A, A, B) + ...
0082         %                  lB^3 * V(B) + 3*lA*lB^2 V(A, B, B)
0083         %
0084         % V(A + B)  = V(A) +   V(B) + 3*V(A, A, B) +  3*V(A, B, B)
0085         % V(A + 2B) = V(A) + 8*V(B) + 6*V(A, A, B) + 12*V(A, B, B)
0086         %
0087         % V(A + 2B) - 2*V(A + B) = -V(A) + 6*V(B) + 6*V(A, B, B)
0088         % V(A, B, B) = (V(A + 2B) - 2*V(A + B) + V(A) - 6 V(B)) / 6
0089         %
0090         % V(A, A, B) = (V(A + B) - V(A) - V(B) - 3*V(A, B, B)) / 3
0091         for iSe = 1:size(hMatrix, 1)
0092             if verbose
0093                 fprintf('\b\b\b\b\b\b%5.1f%%', ...
0094                        (iSe/size(hMatrix, 1))*100);
0095             end
0096             for jSe = (iSe+1):size(hMatrix, 1)
0097                 if incidenceMatrix(iSe, jSe) == 0
0098                     continue
0099                 end
0100                 V_AB = handleH(hMatrix(iSe,:)' + hMatrix(jSe,:)');
0101                 V_A2B = handleH(hMatrix(iSe,:)' + 2*hMatrix(jSe,:)');
0102                 tensor(iSe,jSe,jSe) = (V_A2B - 2*V_AB ...
0103                                  + tensor(iSe,iSe,iSe) ...
0104                                  -6*tensor(jSe,jSe,jSe) ...
0105                                 )/2;
0106                 tensor(iSe,iSe,jSe) = (V_AB ...
0107                                  - tensor(iSe,iSe,iSe) ...
0108                                  - tensor(jSe,jSe,jSe) ...
0109                                  - tensor(iSe,jSe,jSe) ...
0110                                 );
0111             end
0112         end
0113         
0114         % Measure for all combinations of three polytopes
0115         % V(A+B+C) = V(A) + V(B) + V(C) + ...
0116         %            3*V(A,B,B) + 3*V(A,C,C) + 3*V(B,C,C) +
0117         %            6*V(A,B,C)
0118         % V(A,B,C) = (V(A+B+C) - V(A) - V(B) - V(C) ...
0119         %             -3*V(A,B,B) -3*V(A,C,C) -3*V(B,C,C)) / 6
0120         for iSe = 1:size(hMatrix, 1)
0121             if verbose
0122                 fprintf('\b\b\b\b\b\b%5.1f%%', ...
0123                     (iSe/size(hMatrix, 1))*100);
0124             end
0125             for jSe = (iSe+1):size(hMatrix, 1)
0126                 if incidenceMatrix(iSe, jSe) == 0
0127                     continue
0128                 end
0129                 for c = (jSe+1):size(hMatrix, 1)
0130                     if incidenceMatrix(iSe, c) == 0 || ...
0131                        incidenceMatrix(jSe, c) == 0
0132                         continue
0133                     end
0134                     V_ABC = handleH(hMatrix(iSe,:)' + hMatrix(jSe,:)' + hMatrix(c,:)');
0135                     tensor(iSe,jSe,c) = (V_ABC ...
0136                                      -tensor(iSe,iSe,iSe) ...
0137                                      -tensor(jSe,jSe,jSe) ...
0138                                      -tensor(c,c,c) ...
0139                                      -tensor(iSe,jSe,jSe) ...
0140                                      -tensor(iSe,c,c) ...
0141                                      -tensor(jSe,c,c) ...
0142                                      -tensor(iSe,iSe,jSe) ...
0143                                      -tensor(iSe,iSe,c) ...
0144                                      -tensor(jSe,jSe,c) ...
0145                                     );
0146                 end
0147             end
0148         end
0149 end
0150 if verbose, fprintf('\b\b\b\b\b\b');end
0151 end
Generated on Sat 18-Jul-2015 16:45:31 by m2html © 2005