Minktensor2D

computeMinkowskiTensorEstimators.m

function [Phi000,Phi010,Phi020,Phi001,Phi011,Phi002,Phi100,Phi110,Phi120,Phi101,Phi111,Phi102] = computeMinkowskiTensorEstimators(radii,innerPoints,outerPoints,VoronoiIndices,VoronoiVertices,a)
% [Phi000,Phi010,Phi020,Phi001,Phi011,Phi002,Phi100,Phi110,Phi120,Phi101,Phi111,Phi102] = 
% computeMinkowskiTensorsEstimators(radii,innerPoints,outerPoints,VoronoiIndices,VoronoiVertices,a,Phi20,Phi21,Phi22)
% returns all computed estimators (not including the volume tensors) for 
% the Minkowski tensors of a set with positive reach given an array 
% containing three radii below the reach of the set, arrays with the inner 
% and outer points of the digitisation of the set, an array VoronoiVertices 
% with the Voronoi vertices of the outer points as well as a cell array 
% containing the indices of the vertices for each sampling point. Finally, 
% the lattice distance a must be specified.

% Necessary .m-files:
% computeVoronoiTensorMeasures.m

% Compute the Voronoi tensor measures
[VR00,VR20,VR02,VR01,VR10,VR11] = ...
    computeVoronoiTensorMeasures(radii,innerPoints,outerPoints,VoronoiIndices,VoronoiVertices);

Kappa = [1 2 pi 4*pi/3 pi^2/2]; % Volumes of the unit ball in Euclidean n-space where n equals the column number of Kappa minus 1

% Radii and Voronoi tensor measures for the non-scaled digitisation
radii = a*radii;
VR00 = a^(2)*VR00;
VR20 = a^(4)*VR20;
VR02 = a^(4)*VR02;
VR01 = a^(3)*VR01;
VR10 = a^(3)*VR10;
VR11 = a^(4)*VR11;

% Inverses of the constant matrices for s=0,1,2
M0 = [Kappa(1) Kappa(2)*radii(1) Kappa(3)*radii(1)^2;...
      Kappa(1) Kappa(2)*radii(2) Kappa(3)*radii(2)^2;...
      Kappa(1) Kappa(2)*radii(3) Kappa(3)*radii(3)^2];
  
M1 = [Kappa(2)*radii(1) Kappa(3)*radii(1)^2 Kappa(4)*radii(1)^3;...
      Kappa(2)*radii(2) Kappa(3)*radii(2)^2 Kappa(4)*radii(2)^3;...
      Kappa(2)*radii(3) Kappa(3)*radii(3)^2 Kappa(4)*radii(3)^3];
  
M2 = [Kappa(3)*radii(1)^2 Kappa(4)*radii(1)^3 Kappa(5)*radii(1)^4;...
      Kappa(3)*radii(2)^2 Kappa(4)*radii(2)^3 Kappa(5)*radii(2)^4;...
      Kappa(3)*radii(3)^2 Kappa(4)*radii(3)^3 Kappa(5)*radii(3)^4]; 

% Solving the systems of linear equations for the Minkowski tensors
% r = s = 0
b00 = [VR00(1) ; VR00(2) ; VR00(3)];
Phi00 = M0\b00;

Phi100 = Phi00(2)*(abs(Phi00(2))>1e-6);
Phi000 = Phi00(3)*(abs(Phi00(3))>1e-6);

% r = 2, s = 0
b201 = [VR20(1,1) ; VR20(1,2) ; VR20(1,3)];
EtPhi20 = M0\b201;
b202 = [VR20(2,1) ; VR20(2,2) ; VR20(2,3)];
ToPhi20 = M0\b202;
b203 = [VR20(3,1) ; VR20(3,2) ; VR20(3,3)];
TrePhi20 = M0\b203;

Phi120 = [EtPhi20(2)*(abs(EtPhi20(2))>1e-6) TrePhi20(2)*(abs(TrePhi20(2))>1e-6);TrePhi20(2)*(abs(TrePhi20(2))>1e-6) ToPhi20(2)*(abs(ToPhi20(2))>1e-6)]/2;
Phi020 = [EtPhi20(3)*(abs(EtPhi20(3))>1e-6) TrePhi20(3)*(abs(TrePhi20(3))>1e-6);TrePhi20(3)*(abs(TrePhi20(3))>1e-6) ToPhi20(3)*(abs(ToPhi20(3))>1e-6)]/2;

% r = 0, s = 2
b021 = [VR02(1,1) ; VR02(1,2) ; VR02(1,3)];
EtPhi02 = M2\b021;
b022 = [VR02(2,1) ; VR02(2,2) ; VR02(2,3)];
ToPhi02 = M2\b022;
b023 = [VR02(3,1) ; VR02(3,2) ; VR02(3,3)];
TrePhi02 = M2\b023;

Phi102 = [EtPhi02(2)*(abs(EtPhi02(2))>1e-6) TrePhi02(2)*(abs(TrePhi02(2))>1e-6);TrePhi02(2)*(abs(TrePhi02(2))>1e-6) ToPhi02(2)*(abs(ToPhi02(2))>1e-6)]/2;
Phi002 = [EtPhi02(3)*(abs(EtPhi02(3))>1e-6) TrePhi02(3)*(abs(TrePhi02(3))>1e-6);TrePhi02(3)*(abs(TrePhi02(3))>1e-6) ToPhi02(3)*(abs(ToPhi02(3))>1e-6)]/2;

% r = 0, s = 1
b011 = [VR01(1,1) ; VR01(1,2) ; VR01(1,3)];
EtPhi01 = M1\b011;
b012 = [VR01(2,1) ; VR01(2,2) ; VR01(2,3)];
ToPhi01 = M1\b011;

Phi101 = [EtPhi01(2)*(abs(EtPhi01(2))>1e-6) ; ToPhi01(2)*(abs(ToPhi01(2))>1e-6)];
Phi001 = [EtPhi01(3)*(abs(EtPhi01(3))>1e-6) ; ToPhi01(3)*(abs(ToPhi01(3))>1e-6)];

% r = 1, s = 0
b101 = [VR10(1,1) ; VR10(1,2) ; VR10(1,3)];
EtPhi10 = M0\b101;
b102 = [VR10(2,1) ; VR10(2,2) ; VR10(2,3)];
ToPhi10 = M0\b102;

Phi110 = [EtPhi10(2)*(abs(EtPhi10(2))>1e-6) ; ToPhi10(2)*(abs(ToPhi10(2))>1e-6)];
Phi010 = [EtPhi10(3)*(abs(EtPhi10(3))>1e-6) ; ToPhi10(3)*(abs(ToPhi10(3))>1e-6)];

% r = 1, s = 1
b111 = [VR11(1,1) ; VR11(1,2) ; VR11(1,3)];
EtPhi11 = M1\b111;
b112 = [VR11(2,1) ; VR11(2,2) ; VR11(2,3)];
ToPhi11 = M1\b112;
b113 = [VR11(3,1) ; VR11(3,2) ; VR11(3,3)];
TrePhi11 = M1\b113;

Phi111 = [EtPhi11(2)*(abs(EtPhi11(2))>1e-6) TrePhi11(2)*(abs(TrePhi11(2))>1e-6) ; TrePhi11(2)*(abs(TrePhi11(2))>1e-6) ToPhi11(2)*(abs(ToPhi11(2))>1e-6)];
Phi011 = [EtPhi11(3)*(abs(EtPhi11(3))>1e-6) TrePhi11(3)*(abs(TrePhi11(3))>1e-6) ; TrePhi11(3)*(abs(TrePhi11(3))>1e-6) ToPhi11(3)*(abs(ToPhi11(3))>1e-6)];

computeTensorMeasureComponents.m

function [componentVR00,componentVR20,componentVR02,componentVR01,componentVR10,componentVR11] = computeTensorMeasureComponents(R,x,VoronoiVertices)
% [ComponentVR00,ComponentVR20,ComponentVR02,ComponentVR01,ComponentVR10,ComponentVR11] = 
% computeTensorMeasureComponents(r,x,vertices) returns the contributions to
% the Voronoi tensors measures of x when x has R-bounded Voronoi vertices
% VoronoiVertices

vertexCounter = size(VoronoiVertices,1);
x1 = x(1);
x2 = x(2);

% Functions for integrals involved in the Minkowski tensor measures in
% polar coordinates
ints211Ice = @(tMin,tMax,rMax) (1/16)*(rMax.^4).*(2*(tMax-tMin)+sin(2*tMax)-sin(2*tMin)); % Primitive function for the coordinate (1,1) for tensors with s=2 ('ice' region)
ints222Ice = @(tMin,tMax,rMax) (1/16)*(rMax.^4).*(2*(tMax-tMin)-sin(2*tMax)+sin(2*tMin)); % Primitive function for the coordinate (2,2) for tensors with s=2 ('ice' region)
ints212Ice = @(tMin,tMax,rMax) ((sin(tMax)^2 - sin(tMin)^2)*(rMax^4))/8; % Primitive function for the coordinate (1,2)=(2,1) for tensors with s=2 ('ice' region)

ints11Ice = @(tMin,tMax,rMax) (rMax^3*(sin(tMax) - sin(tMin)))/3; % Primitive function for the coordinate (1,1) for tensors with s=1 ('ice' region)
ints12Ice = @(tMin,tMax,rMax) -(rMax^3*(cos(tMax) - cos(tMin)))/3; % Primitive function for the coordinate (2,1) for tensors with s=1 ('ice' region)

% If the R-bounded Voronoi cell is a disc, we simply integrate over this
if (~vertexCounter || vertexCounter == 1)
    
    volume = pi*R^2;
    
    % r = s = 0
    componentVR00 = volume;
    
    % r = 2, s = 0
    componentVR20 = [x1^2*volume ; x2^2*volume ; x1*x2*volume];
    
    % r = 0, s = 2
    componentVR02 = [ints211Ice(0,2*pi,R) ; ints222Ice(0,2*pi,R) ; ints212Ice(0,2*pi,R)];
    
    % r = 0, s = 1
    componentVR01 = [ints11Ice(0,2*pi,R) ; ints12Ice(0,2*pi,R)];
    
    % r = 1, s = 0
    componentVR10 = [x1*volume ; x2*volume];
    
    % r = 1, s = 1
    componentVR11 = [x1*componentVR01(1) ; x2*componentVR01(2) ; 0.5*(x1*componentVR01(2) + x2*componentVR01(1))];

else
    
    concatenatedVertices = [VoronoiVertices; VoronoiVertices(1,:)];
    
    % We prepare the components for the Voronoi tensor measures
    componentVR00 = 0;                  % r = s = 0 tensor
    componentVR20 = [0 0 0]';           % r = 2, s = 0 tensor
    componentVR02 = componentVR20;      % r = 0, s = 2 tensor
    componentVR01 = [0 0]';             % r = 0, s = 1 tensor
    componentVR10 = componentVR01;      % r = 1, s = 0 tensor
    componentVR11 = componentVR20;      % r = 1, s = 1 tensor
    
    ints0Ice = @(tMin,tMax,rMax) (1/2)*(rMax.^2).*(tMax-tMin); % Primitive function for tensors with s=0 ('ice' region)
    
    % Functions for integrals involved in the Minkowski tensor measures in
    % Cartesian coordinates of a rotated coordinate system
    ints0Triangle = @(s0,t0,s1) (s0*t0)/2 - (t0*(s0 - s1))/2; % Primitive function for tensors with s=0 ('triangle' region)
    
    ints2XXTriangle = @(s0,t0,s1,xs,xt) (s0*t0*(3*s0^2*xs^2 + 3*s0*t0*xs*xt + t0^2*xt^2))/12 -...
        (t0*(s0 - s1)*(3*s0^2*xs^2 + 2*s0*s1*xs^2 + 3*s0*t0*xs*xt + s1^2*xs^2 + s1*t0*xs*xt + t0^2*xt^2))/12; % Primitive function for the coordinates (1,1) and (2,2) for tensors with s=2 ('triangle' region)
    
    ints2XYTriangle = @(s0,t0,s1,xs,xt,ys,yt) (s0*t0*(6*s0^2*xs*ys + 2*t0^2*xt*yt + 3*s0*t0*xs*yt + 3*s0*t0*xt*ys))/24 - ...
        (t0*(s0 - s1)*(6*s0^2*xs*ys + 2*s1^2*xs*ys + 2*t0^2*xt*yt + 4*s0*s1*xs*ys + 3*s0*t0*xs*yt + 3*s0*t0*xt*ys + s1*t0*xs*yt + s1*t0*xt*ys))/24; % Primitive function for the coordinates (1,2)=(2,1) for tensors with s=2 ('triangle' region)
    
    ints1Triangle = @(s0,t0,s1,xs,xt) (s0*t0*(2*s0*xs + t0*xt))/6 - (t0*(s0 - s1)*(2*s0*xs + s1*xs + t0*xt))/6;  % Primitive function for the coordinates (1,1) and (2,1) for tensors with s=1 ('triangle' region)
    
    % Coordinates are translated to a coordinate system with x as origin.
    s = concatenatedVertices(1,1:2)-x;
    
    for jjj = 1:vertexCounter
        v = s;                                  % Current vertex
        s = concatenatedVertices(jjj+1,1:2)-x;  % Subsequent vertex
        
        vType = concatenatedVertices(jjj,3);    % Connection type of the current vertex
        
        % The vertex connection type now determines which type of integral we
        % need to compute
        
        % Connection type 1 (arc)
        if vType
            
            % Find the angles of the two points
            if v(1) == 0
                vAngle = norm(sign(v(2)))*(pi/2 + (mod(sign(v(2))+2,3))*pi);
            else
                vAngle = atan(v(2)/v(1)) + (v(1)<0)*pi + (v(1)>0 && v(2)<0)*2*pi;
            end
            
            if s(1) == 0
                sAngle = norm(sign(s(2)))*(pi/2 + (mod(sign(s(2))+2,3))*pi);
            else
                sAngle = atan(s(2)/s(1)) + (s(1)<0)*pi + (s(1)>0 && s(2)<0)*2*pi;
            end
            
            % Adjustment of angles for the integration
            if vAngle > sAngle
                sAngle = sAngle + 2*pi;
            end
            
            % Computation of the Voronoi tensor measure contributions from
            % the current R-bounded Voronoi vertex
            integrals0Ice = ints0Ice(vAngle,sAngle,R); % Integral calculated for all component terms where s = 0 for inner Voronoi vertices
            integrals11Ice = ints11Ice(vAngle,sAngle,R);
            integrals12Ice = ints12Ice(vAngle,sAngle,R);
            
            % r = s = 0
            componentVR00 = componentVR00 + integrals0Ice;
            
            % r = 2, s = 0
            componentVR20(1) = componentVR20(1) + x1^2*integrals0Ice;
            componentVR20(2) = componentVR20(2) + x2^2*integrals0Ice;
            componentVR20(3) = componentVR20(3) + x1*x2*integrals0Ice;
            
            % r = 0, s = 2
            componentVR02(1) = componentVR02(1) + ints211Ice(vAngle,sAngle,R);
            componentVR02(2) = componentVR02(2) + ints222Ice(vAngle,sAngle,R);
            componentVR02(3) = componentVR02(3) + ints212Ice(vAngle,sAngle,R);
            
            % r = 0, s = 1
            componentVR01(1) = componentVR01(1) + integrals11Ice;
            componentVR01(2) = componentVR01(2) + integrals12Ice;
            
            % r = 1, s = 0
            componentVR10(1) = componentVR10(1) + x1*integrals0Ice;
            componentVR10(2) = componentVR10(2) + x2*integrals0Ice;
            
            % r = 1, s = 1
            componentVR11(1) = componentVR11(1) + x1*integrals11Ice;
            componentVR11(2) = componentVR11(2) + x2*integrals12Ice;
            componentVR11(3) = componentVR11(3) + 0.5*(x1*integrals12Ice + x2*integrals11Ice);
            
        % Connection type 0 (line segment)    
        else
            
            % Constants for computing the Voronoi tensor measure
            % contribution from the current R-bounded Voronoi vertex
            normV = norm(v);
            u = v/normV;
            v = [-u(2) u(1)];
            s0 = dot(s,u);
            t0 = dot(s,v);
            xs = dot([1 0],u);
            xt = dot([1 0],v);
            ys = dot([0 1],u);
            yt = dot([0 1],v);
            
            integrals0Triangle = ints0Triangle(s0,t0,normV);
            integrals1xTriangle = ints1Triangle(s0,t0,normV,xs,xt);
            integrals2yTriangle = ints1Triangle(s0,t0,normV,ys,yt);
            
            % r = s = 0
            componentVR00 = componentVR00 + integrals0Triangle;
            
            % r = 2, s = 0
            componentVR20(1) = componentVR20(1) + x1^2*integrals0Triangle;
            componentVR20(2) = componentVR20(2) + x2^2*integrals0Triangle;
            componentVR20(3) = componentVR20(3) + x1*x2*integrals0Triangle;
            
            % r = 0, s = 2
            componentVR02(1) = componentVR02(1) + ints2XXTriangle(s0,t0,normV,xs,xt);
            componentVR02(2) = componentVR02(2) + ints2XXTriangle(s0,t0,normV,ys,yt);
            componentVR02(3) = componentVR02(3) + ints2XYTriangle(s0,t0,normV,xs,xt,ys,yt);
            
            % r = 0, s = 1
            componentVR01(1) = componentVR01(1) + integrals1xTriangle;
            componentVR01(2) = componentVR01(2) + integrals2yTriangle;
            
            % r = 1, s = 0
            componentVR10(1) = componentVR10(1) + x1*integrals0Triangle;
            componentVR10(2) = componentVR10(2) + x2*integrals0Triangle;
            
            % r = 1, s = 1
            componentVR11(1) = componentVR11(1) + x1*integrals1xTriangle;
            componentVR11(2) = componentVR11(2) + x2*integrals2yTriangle;
            componentVR11(3) = componentVR11(3) + 0.5*(x1*integrals2yTriangle + x2*integrals1xTriangle);
            
        end
    end
end

computeVolumeTensorEstimators.m

function [Phi200,Phi210,Phi220] = computeVolumeTensorEstimators(outerPoints,innerPoints,a)
% [Phi200,Phi210,Phi220] = 
% computeVolumeTensorEstimators(outerPoints,innerPoints,a) returns
% estimators for the non-trivial volume tensors of rank at most two in
% dimension two based on arrays outerPoints and innerPoints of coordinates 
% of all outer and inner sampling points of a ditigisation with lattice 
% distance a

numberOfPoints = size(outerPoints,1) + size(innerPoints,1);
allPoints = [outerPoints;innerPoints];

% (r=0)-tensor estimator
sumxr0 = numberOfPoints;
Phi20 = a^2*sumxr0;

% (r=1)-tensor estimator
sumxr1 = sum(allPoints,1);
Phi21 = a^3*sumxr1;

% (r=2)-tensor estimator
sumxr21 = sum(allPoints.^2,1);
sumxr22 = sum(allPoints(:,1).*allPoints(:,2),1);
Phi22a = (a^4/2)*sumxr21;
Phi22b = (a^4/2)*sumxr22;
Phi22 = [Phi22a,Phi22b];

% Tensors very close to zero are rounded off
Phi200 = Phi20*(Phi20>1e-6);
Phi210 = [Phi21(1)*(abs(Phi21(1))>1e-6) ; Phi21(2)*(abs(Phi21(2))>1e-6)];
Phi220 = [Phi22(1)*(abs(Phi22(1))>1e-6) Phi22(3)*(abs(Phi22(3))>1e-6);Phi22(3)*(abs(Phi22(3))>1e-6) Phi22(2)*(abs(Phi22(2))>1e-6)];

computeVoronoiCells.m

function [outerPoints,innerPoints,VoronoiIndices,VoronoiVertices] = computeVoronoiCells(digitisation,r,origin)
% [outerPoints,innerPoints,VoronoiIndices,VoronoiVertices] =
% computeVoronoiCells(digitisation,r,origin) returns an array outerPoints
% of coordinates of the sampling points in digitisation which have a
% non-trivial Voronoi cell; an array innerPoints of coordinates of inner
% points in digitisation, which have trivial Voronoi cell; an array
% VoronoiVertices of vertices of Voronoi cells of points in
% [outerPoints;innerPoints], and a cell array VoronoiIndices whose n'th row
% correspondonds to the indices of the Voronoi vertices in VoronoiVertices
% of the n'th point in [outerPoints;innerPoints]

% Necessary .m-files:
% convertDigitisationToCoordinates.m

% Convert digitisation to points in the plane
[innerPoints,outerPoints,VoronoiRelevant] = convertDigitisationToCoordinates(digitisation,origin);

if isnan(innerPoints)
    allPoints = outerPoints;
else
    allPoints = [innerPoints;outerPoints];
end

numberOfPoints = size(allPoints,1);

% Find the extremal coordinates of the digitisation
xMin = min(allPoints(1:numberOfPoints,1));
xMax = max(allPoints(1:numberOfPoints,1));
yMin = min(allPoints(1:numberOfPoints,2));
yMax = max(allPoints(1:numberOfPoints,2));

xMid = (xMax - xMin)/2 + xMin;
yMid = (yMax - yMin)/2 + yMin;

% We will insert 'border points' far away from the sampling points in order
% to obtain bounded Voronoi cells.
% The distance from digitisation to border points is set to three times the 
% maximal radius. Thus the relevant parts of the Voronoi cells of the 
% digitisation will not be affected
borderPointDistance = 3*ceil(r);

% Create border points
borderPoints = zeros(8,2);
borderPoints(1,:) = [xMin-borderPointDistance,yMid];
borderPoints(2,:) = [xMid,yMax+borderPointDistance];
borderPoints(3,:) = [xMax+borderPointDistance,yMid];
borderPoints(4,:) = [xMid,yMin-borderPointDistance];
borderPoints(5,:) = [xMin-borderPointDistance,yMin-borderPointDistance];
borderPoints(6,:) = [xMin-borderPointDistance,yMax+borderPointDistance];
borderPoints(7,:) = [xMax+borderPointDistance,yMin-borderPointDistance];
borderPoints(8,:) = [xMax+borderPointDistance,yMax+borderPointDistance];

% Create the Voronoi vertices and indices for all relevant sampling points
% of the digitisation
if isnan(VoronoiRelevant) 
    VoronoiRelevantPoints = []; 
    
% If there are Voronoi relevant inner points, we use them when computing 
% the Voronoi cells
else
   VoronoiRelevantPoints = VoronoiRelevant;
end

% The Voronoi diagram is now computed for all relevant points
[VoronoiVertices,excessVoronoiIndices] = voronoin([outerPoints;borderPoints;VoronoiRelevantPoints]);

% Clean up VoronoiIndices by removing indices correponding to artificial
% points
numberOfRelevantIndices = size(outerPoints,1);
VoronoiIndices = excessVoronoiIndices(1:numberOfRelevantIndices,:);

computeVoronoiTensorMeasures.m

function [VR00,VR20,VR02,VR01,VR10,VR11] = computeVoronoiTensorMeasures(radii,innerPoints,outerPoints,VoronoiIndices,VoronoiVertices)
% [VR00,VR20,VR02,VR01,VR10,VR11] =
% computeVoronoiTensorMeasures(radii,innerPoints,outerPoints,VoronoiIndices,VoronoiVertices)
% returns the Voronoi tensor measures of a set with positive reach given an
% array containing three radii below the reach of the set and greater than 
% 1/sqrt(2), arrays with the inner and outer points of the digitisation of 
% the set, a cell array VoronoiVertices with the Voronoi vertices of the 
% outer points as well as an array containing the indices of the vertices 
% for each sampling point.

% Necessary .m-files:
% verticesToCounterclockwiseOrder.m
% intersectLinesegmentCircle.m
% RBoundedVoronoiVertex.m
% computeTensorMeasureComponents.m

numberOfInnerPoints = size(innerPoints,1)*(~isnan(innerPoints(1)));
numberOfOuterPoints = size(outerPoints,1);

% Arrays containing the entries for the Voronoi tensor measures are
% prepared. The three columns correspond to the three radii.
VR00 = [0 0 0];     % r=0, s=0 (0-tensor)
VR20 = zeros(3,3);  % r=2, s=0 (2-tensor)
VR02 = VR20;        % r=0, s=2 (2-tensor)
VR01 = zeros(2,3);  % r=0, s=1 (1-tensor)
VR10 = VR01;        % r=1, s=0 (1-tensor)
VR11 = VR20;        % r=0, s=1 (1-tensor)

% First, contributions from the inner points to all Voronoi tensor measures
% are calculated
if numberOfInnerPoints
    for iii = 1:numberOfInnerPoints
        
        x = innerPoints(iii,:);
        x1 = x(1);
        x2 = x(2);
        
        for kkk = 1:3
            VR00(kkk) = VR00(kkk) + 1;
            VR20(:,kkk) = VR20(:,kkk) + [x1^2 ; x2^2 ; x1*x2];
            VR02(:,kkk) = VR02(:,kkk) + [1 ; 1 ; 0]*(1/12);
            VR10(:,kkk) = VR10(:,kkk) + [x1 ; x2];
        end
    end
end

% Next, contributions from the outer points to all Voronoi tensor measures
% are calculated
for iii = 1:numberOfOuterPoints
    x = outerPoints(iii,:);
    
    indices = VoronoiIndices{iii}; % Indices of the Voronoi vertices of x
    vertices = VoronoiVertices(indices,:); % Voronoi vertices of x
    vertices = verticesToCounterclockwiseOrder(vertices,x); % Make sure the Voronoi vertices are in counterclockwise order with respect to x
    nVertices = size(vertices,1);
    
    % Arrays containing the contributions to the entries of the Voronoi
    % tensor measures from x are prepared
    componentsVR00 = [0 0 0];
    componentsVR20 = zeros(3,3);
    componentsVR02 = componentsVR20;
    componentsVR01 = zeros(2,3);
    componentsVR10 = zeros(2,3);
    componentsVR11 = componentsVR20;
    
    for kkk = 1:3
        R = radii(kkk);
        
        column = ones(nVertices,1);
        coordinateMatrix = [x(1)*column x(2)*column];
        
        vTypes = (R < sqrt(sum((vertices-coordinateMatrix).^2,2)));  % Assign a 1 to vertices that are further away from the current coordinate than R; 0 otherwise
        
        ConcatenatedVertices = [vertices(nVertices,:); vertices; vertices(1,:)]; % Concatenated array with Vertices where the last row is copied to the first row and the first row is copied to the last row
        
        % We now adjust all Voronoi vertices to vertices lying on the
        % R-bounded Voronoi cell. Additionally, we determine whether a
        % vertex is connected to its subsequent neighbour by an arc
        % (type 1) or a line (type 0).
        
        adjustedVertices = zeros(2*nVertices,3); % For saving the adjusted vertices and their integration type. Each vertex is adjusted to at most 2 new vertices
        vertexCounter = 0;
        
        for jjj = 1:nVertices
            % The types of the vertex and its two neighbours decide how the
            % vertex should be adjusted
            v = vertices(jjj,:); % v is the jjj'th Voronoi vertex of x
            vType = vTypes(jjj); % vType is the type of v
            
            % If v is inside the R-ball with centre x, it is a vertex of
            % the R-bounded Voronoi cell of type 0
            if ~vType
                vertexCounter = vertexCounter + 1;
                adjustedVertices(vertexCounter,1:2) = v;
                adjustedVertices(vertexCounter,3) = 0;
                
            else
                p = ConcatenatedVertices(jjj,:);
                s = ConcatenatedVertices(jjj+2,:);
                
                % Check whether sv and vp intersect the R-ball
                intersectionsPV = intersectLinesegmentCircle(v,p,x,R);
                intersectionsVS = intersectLinesegmentCircle(v,s,x,R);
                
                if (~intersectionsPV && ~intersectionsVS)   % If none of the line segments intersect, the vertex is not relevant for our computations
                    continue
                    
                else
                    
                    if intersectionsPV
                        v1 = RBoundedVoronoiVertex(p,v,x,R); % Adjust v to the point where pv intersects the R-ball
                        
                        if isequal(v1,p)
                            continue
                        else
                            vertexCounter = vertexCounter + 1;
                            adjustedVertices(vertexCounter,1:2) = v1;
                            adjustedVertices(vertexCounter,3) = 1;
                        end
                    end
                    
                    if intersectionsVS
                        v2 = RBoundedVoronoiVertex(s,v,x,R); % Adjust v to the point where vs intersects the R-ball
                        
                        if isequal(v2,s)
                            continue
                        else
                            vertexCounter = vertexCounter + 1;
                            adjustedVertices(vertexCounter,1:2) = v2;
                            adjustedVertices(vertexCounter,3) = 0;
                        end
                    end
                end
            end
        end
        
        % Remove all zero-rows from adjustedVertices to obtain the vertices
        % of R-bounded Voronoi cell of x
        RboundedVertices = adjustedVertices(1:vertexCounter,:);
        
        % Add the terms for this coordinate and radius to the Voronoi
        % tensor measures
        [componentsVR00(kkk),componentsVR20(:,kkk),componentsVR02(:,kkk),componentsVR01(:,kkk),componentsVR10(:,kkk),componentsVR11(:,kkk)] = computeTensorMeasureComponents(R,x,RboundedVertices);
        
    end
    
    VR00(1) = VR00(1) + componentsVR00(1);
    VR00(2) = VR00(2) + componentsVR00(2);
    VR00(3) = VR00(3) + componentsVR00(3);
    
    VR20(:,1) = VR20(:,1) + componentsVR20(:,1);
    VR20(:,2) = VR20(:,2) + componentsVR20(:,2);
    VR20(:,3) = VR20(:,3) + componentsVR20(:,3);
    
    VR02(:,1) = VR02(:,1) + componentsVR02(:,1);
    VR02(:,2) = VR02(:,2) + componentsVR02(:,2);
    VR02(:,3) = VR02(:,3) + componentsVR02(:,3);
    
    VR01(:,1) = VR01(:,1) + componentsVR01(:,1);
    VR01(:,2) = VR01(:,2) + componentsVR01(:,2);
    VR01(:,3) = VR01(:,3) + componentsVR01(:,3);
    
    VR10(:,1) = VR10(:,1) + componentsVR10(:,1);
    VR10(:,2) = VR10(:,2) + componentsVR10(:,2);
    VR10(:,3) = VR10(:,3) + componentsVR10(:,3);
    
    VR11(:,1) = VR11(:,1) + componentsVR11(:,1);
    VR11(:,2) = VR11(:,2) + componentsVR11(:,2);
    VR11(:,3) = VR11(:,3) + componentsVR11(:,3);
    
end

convertDigitisationToCoordinates.m

function [innerPoints,outerPoints,VoronoiRelevant] = convertDigitisationToCoordinates(digitisation,origin)
% [innerPoints,outerPoints,VoronoiRelevant] =
% convertDigitisationToCoordinates(digitisation,origin) returns an array
% outerPoints of coordinates of the sampling points in digitisation which
% have a non-trivial Voronoi cell; an array innerPoints of coordinates of
% inner points in digitisation, which have trivial Voronoi cell, and an
% array of those points in digitisation which define the non-trivial
% Voronoi cells of the outer points

digitisation = flipud(digitisation); % Matrix indices are read from the
% upper left corner whereas planar coordinates are read from the bottom
% left corner, so we flip digitisation horisontally

numberOfRows = size(digitisation,1);
numberOfColumns = size(digitisation,2);

originx = origin(2) - 1;
originy = numberOfRows - origin(1);

% Prepare arrays for the three types of sampling points and allocate space
% for the maximal number of rows
arrayInnerPoints = zeros(numberOfRows*numberOfColumns,2);
arrayOuterPoints = arrayInnerPoints;
arrayVoronoiRelevant = arrayInnerPoints;

% Counters for the number of sampling points of each type
countInner = 0;
countOuter = 0;
countVoronoi = 0;

for iii = 1:(numberOfRows)
    for jjj = 1:(numberOfColumns)
        
        x = digitisation(iii,jjj);
        
        % Only 1's are converted to coordinates
        if x
            
            notBorder = (iii-1)*(jjj-1)*(iii-numberOfRows)*(jjj-numberOfColumns); % false if x lies in first or last column and/or row
            
            if notBorder
                
                xr = digitisation(iii+1,jjj); % point to the right of x
                xl = digitisation(iii-1,jjj); % point to the left of x
                xu = digitisation(iii,jjj+1); % point above x
                xd = digitisation(iii,jjj-1); % point below x
                
                inner = xr*xl*xu*xd; % true if x is an inner point
                
                if inner
                    
                    countInner = countInner + 1;
                    arrayInnerPoints(countInner,:) = [jjj-1,iii-1];
                    
                    innerBorder = ~((iii-2)*(jjj-2)*(iii-numberOfRows+1)*...
                        (jjj-numberOfColumns+1)); % true if x is a four-neighbour of a non-inner point. If so, x determines the Voronoi vertices of its neighbour
                    
                    if innerBorder
                        
                        countVoronoi = countVoronoi + 1;
                        arrayVoronoiRelevant(countVoronoi,:) = [jjj-1,iii-1];
                        
                    else
                        % Four-neighbours of the four-neighbours of x
                        xrr = digitisation(iii+2,jjj);
                        xru = digitisation(iii+1,jjj+1);
                        xrd = digitisation(iii+1,jjj-1);
                        xll = digitisation(iii-2,jjj);
                        xlu = digitisation(iii-1,jjj+1);
                        xld = digitisation(iii-1,jjj-1);
                        xuu = digitisation(iii,jjj+2);
                        xdd = digitisation(iii,jjj-2);
                        
                        % Check which of the four-neighbours of x are inner
                        % points
                        rightInner = xrr*xru*xrd;
                        leftInner = xll*xlu*xld;
                        upInner = xuu*xru*xlu;
                        downInner = xdd*xrd*xld;
                        
                        relevant = ~(rightInner*leftInner*upInner*downInner);
                        
                        if relevant
                            
                            countVoronoi = countVoronoi + 1;
                            arrayVoronoiRelevant(countVoronoi,:) = [jjj-1,iii-1];
                            
                        end
                        
                    end
                else
                    
                    countOuter = countOuter + 1;
                    arrayOuterPoints(countOuter,:) = [jjj-1,iii-1];
                    
                end
            else
                
                countOuter = countOuter + 1;
                arrayOuterPoints(countOuter,:) = [jjj-1,iii-1];
                
            end
        end
    end
end

if countInner
    
    % innerPoints trims arrayInnerPoints
    innerPoints = arrayInnerPoints(1:countInner,:);
    innerPoints(:,1) = innerPoints(:,1) - originx;
    innerPoints(:,2) = innerPoints(:,2) - originy;
    
else
    
    innerPoints = NaN;
    
end

if countVoronoi
    
    % VoronoiRelevant trims arrayVoronoiRelevant
    VoronoiRelevant = arrayVoronoiRelevant(1:countVoronoi,:);
    VoronoiRelevant(:,1) = VoronoiRelevant(:,1) - originx;
    VoronoiRelevant(:,2) = VoronoiRelevant(:,2) - originy;
    
else
    
    VoronoiRelevant = NaN;
    
end

% outerPoints trims arrayOuterPoints
outerPoints = arrayOuterPoints(1:countOuter,:);
outerPoints(:,1) = outerPoints(:,1) - originx;
outerPoints(:,2) = outerPoints(:,2) - originy;

intersectLineCircle.m

function [xout,yout] = intersectLineCircle(point1x,point1y,point2x,point2y,centrex,centrey,radius)
% [xout,yout] = intersectLineCircle(point1x,point1y,point2x,point2y,centrex,centrey,radius)
% returns arrays xout and yout with the x- and y-coordinates of the
% intersection points of a line with a circle in the plane. If the line is
% tangent to the circle, the point of intersection is stored as two
% identical points; if there are no points of intersection, the arrays 
% contain NaN as entries.
% Input arguments are coordinates (point1x,point1y) and (point2x,point2y)
% of two points on the line, the coordinates (centrex,centrey) of the 
% centre of the circle, and the radius of the circle. 

% Translate the coordinates to a system with the centre of the circle at
% the origin
x1 = point1x-centrex;
y1 = point1y-centrey;
x2 = point2x-centrex;
y2 = point2y-centrey;

% Compute the parameters of the system
dx = x2-x1;
dy = y2-y1;
dr = sqrt(dx^2+dy^2);
determinant = x1*y2-x2*y1;
discriminant = radius^2*dr^2-determinant^2;
epsilon = 1e-14; % Precision

% Compute the points of intersection (i.e. the solutions to the second
% order equation)
if discriminant > epsilon % Two points of intersection
    x = [(determinant*dy+sgn(dy)*dx*sqrt(discriminant))/dr^2 ; (determinant*dy-sgn(dy)*dx*sqrt(discriminant))/dr^2];
    y = [(-determinant*dx+abs(dy)*sqrt(discriminant))/dr^2 ; (-determinant*dx-abs(dy)*sqrt(discriminant))/dr^2];
elseif abs(discriminant) < epsilon % One point of intersection (the line is tangent to the circle)
    x = determinant*dy/dr^2*ones(2,1);
    y = -determinant*dx/dr^2*ones(2,1);
else % No points of intersection
    x = NaN*ones(2,1);
    y = NaN*ones(2,1);
end

% Translate the points of intersection back to a system with the centre of
% the circle at (centrex,centrey)
xout = x+centrex;
yout = y+centrey;

function s = sgn(x)
% sgn is a signum function which returns -1 for all negative numbers and 1
% otherwise
s = sign(sign(x)+0.5);

intersectLinesegmentCircle.m

function [intersection] = intersectLinesegmentCircle(p1,p2,c,r)
% [intersections] = intersectLinesegmentCircle(p1,p2,c,r) determines
% whether the line segment from p1 to p2 intersects the circle with radius
% r and centre c. If there are no intersection points, intersection is
% false; otherwise, intersection is true

if isequal(p1,p2)
    error('The two points are equal. No line segment can be determined.')
end

p = p2-p1; % Vector from p1 to p2
np = norm(p);
punit = p/np;
l = c - p1; % Vector from p1 to the centre of the circle
nProj = dot(l,punit); % Length of the projection of l onto p

% Find the closest point on the line segment to the centre c
if nProj <= 0
    closest = p1;
elseif nProj > np
    closest = p2;
else
    Proj = punit*nProj;
    closest = Proj + p1;
end

distance = c - closest; % Vector from the closest point to the centre of the circle
nDistance = norm(distance);

% Find out whether the line segment and the circle intersect
if nDistance > r
    intersection = 0;  % The line segment does not intersect the circle
elseif nDistance < r
    intersection = 2;  % The line segment intersects the circle in two points
else
    intersection = 1;  % The line segment is tangent to the circle
end

Minktensor2D.m

% Minktensor2D computes estimators for Minkowski tensors in dimension 2 for
% compact sets with positive reach given a digitisation of the set. 

% User must input a digitisation in form of a datafile containing an array
% of 0's and 1's, where a 1 indicates that the corresponding lattice point
% is contained within the set with positive reach; a 0 that it is not. The
% user must also specify the lattice distance as well an the origin, where
% the latter must be expressed as a row and column number corresponding to
% the 0-1-array (where these numbers need not be integers nor positive
% since the origin may not coincide with a lattice point). Finally, a lower 
% and an upper radius for use in the algorithm must be chosen. The radii 
% must not be lower than 1/sqrt(2) times the lattice distance, and neither 
% radius must exceed the reach of the digitised object.

% Primary (necessary) .m-files:
% computeMinkowskiTensorEstimators.m
% computeVolumeTensorEstimators.m
% computeVoronoiCells.m
% retrieveData.m
% retrieveRadii.m
% 
% Secondary .m-files used by the primary .m-files
% computeTensorMeasureComponents.m
% computeVoronoiTensorMeasures.m
% convertDigitisationToCoordinates.m
% intersectLineCircle.m
% intersectLinesegmentCircle.m
% RBoundedVoronoiVertex.m
% verticesToCounterclockwiseOrder.m

% Retrieve digitisation data, lattice distance, and origin from user
[digitisation,a,origin] = retrieveData();

% Retrieve radii from user
radii = retrieveRadii(a);

% Find planar coordinates of the sampling points in digitisation as well as 
% the vertices and indices of the Voronoi cells of the sampling points
[outerPoints,innerPoints,VoronoiIndices,VoronoiVertices] = computeVoronoiCells(digitisation,radii(3),origin);

% Compute estimators for the volume tensors
[Phi200,Phi210,Phi220] = computeVolumeTensorEstimators(outerPoints,innerPoints,a);

% Compute estimators for the Minkowski tensors of rank at most two
[Phi000,Phi010,Phi020,Phi001,Phi011,Phi002,Phi100,Phi110,Phi120,Phi101,Phi111,Phi102] = computeMinkowskiTensorEstimators(radii,innerPoints,outerPoints,VoronoiIndices,VoronoiVertices,a);

% Display all estimators
disp('The estimators for the Minkowski tensors are:')

disp('Phi000'),disp(Phi000);
disp('Phi010'),disp(Phi010);
disp('Phi020'),disp(Phi020);
disp('Phi001'),disp(Phi001);
disp('Phi011'),disp(Phi011);
disp('Phi002'),disp(Phi002);

disp('Phi100'),disp(Phi100);
disp('Phi110'),disp(Phi110);
disp('Phi120'),disp(Phi120);
disp('Phi101'),disp(Phi101);
disp('Phi111'),disp(Phi111);
disp('Phi102'),disp(Phi102);

disp('Phi200'),disp(Phi200);
disp('Phi210'),disp(Phi210);
disp('Phi220'),disp(Phi220);

RBoundedVoronoiVertex.m

function [correctedVertex] = RBoundedVoronoiVertex(neighbourVertex,vertex,c,R)
% [correctedVertex] = RBoundedVoronoiVertex(neighbourVertex,vertex,c,r)
% returns the R-bounded Voronoi vertex corresponding to vertex by
% translating it to the nearest point of intersection between the circle 
% with radius R and centre c and the line segment from neighbourVertex to vertex.

% Necessary .m-files:
% intersectLineCircle.m

[xIntersect,yIntersect] = intersectLineCircle(neighbourVertex(1),neighbourVertex(2),vertex(1),vertex(2),c(1),c(2),R); % Coordinates of intersection of the line and the circle with given radius and centre

[minimalDistance,index] = min([norm(vertex - [xIntersect(1),yIntersect(1)]) norm(vertex - [xIntersect(2),yIntersect(2)])]); % Find the index of the intersection point closest to vertex
correctedVertex = [xIntersect(index) yIntersect(index)];

retrieveData.m

function [digitisation,latDist,origin] = retrieveData()
% digitisation = retrieveData() returns an array corresponding to a loaded 
% data file specified by the user

% [dagitisation,latDist,origin] = retrieveData() also returns the lattice
% distance latDist and origin of the digitisation both specified by the
% user

% Ask user to choose a data file of type .dat
[DataFileName,DataPathName] = uigetfile('*.dat','Select data file (.dat) containing the digitisation');

% Load the specified data file
digitisation = load([DataPathName DataFileName]);

% Determine the size of the data file
DataColumns = size(digitisation,2);
DataRows = size(digitisation,1);

% Ask user to input the lattice distance of the digitisation
latDist = input('>> Please enter the lattice distance of your digitisation. ');

% If latDist is not a positive number, request a new lattice distance value
while (( ~isnumeric(latDist) ) || ( isempty(latDist) ) || (latDist <= 0))
    latDist = input('>> Your input is not a positive number. Please try again: ');
end

disp('Entries of the data file are given by coordinates [row number,column number]')
disp('read from top to bottom respectively left to right.')
disp('Coordinates of the origin are required. Since the origin may not be a data point,' )
disp('in the following, row and column numbers are allowed to be any real number.')

% Ask the user for the origin of the data file
originr = input('>> Please enter the row number of the origin in the data file. '); % 'Row number' of the origin

while (( originr < 1 ) || (originr > DataRows))
    originr = input('>> Your input is not a valid row number. Please try again: ');
end

originc = input('>> Please enter the column number of the origin in the data file. '); % 'Column number' of the origin

while (( originc < 1 ) || (originc > DataColumns))
    originc = input('>> Your input is not a valid column number. Please try again: ');
end

origin = [originr,originc];

retrieveRadii.m

function [radii] = retrieveRadii(a)
% radii = retrieveRadii(a) returns three radii [R0,R1,R2] (where the R0<R2 
% are specified by the user and R1=(R0+R2)/2) for use in the algorithm 
% Minktensor2D for a digitisation with lattice distance a

% Ask user for lower radius for use in the algorithm
lowerRadius = input('>> Please enter a lower radius for use in the algorithm. ');

% Continue to prompt for lowerRadius if user does not input a positive
% number greater than a/sqrt(2)
while ( ~isnumeric(lowerRadius) ) || ( isempty(lowerRadius) ) || ( lowerRadius <= (a/sqrt(2)) )
    disp('>> Your input is not a positive number greater than the lattice distance divided by the square root of two.')
    lowerRadius = input('Please try again: ');
end

% Ask user for upper radius for use in the algorithm
disp('>> Please enter the upper radius for use in the algorithm.')
upperRadius = input('(Number should not exceed the reach of the digitised object!) ');

% Continue to prompt for upperRadius if user does not input a positive
% number greater than lowerRadius
while ( ~isnumeric(upperRadius) ) || ( isempty(upperRadius) ) || ( upperRadius <= lowerRadius ) || ( isnan(upperRadius) ) || ( isinf(upperRadius) )
    upperRadius = input('>> Your input is not a positive number (greater than the lower radius chosen). Please try again: ');
end

% Scale radii with lattice distance
upperRadius = upperRadius/a;
lowerRadius = lowerRadius/a;

radii = linspace(lowerRadius,upperRadius,3);

verticesToCounterclockwiseOrder.m

function [counterclockwiseVertices] = verticesToCounterclockwiseOrder(vertices,referencePoint)
% [counterClockwiseVertices] = 
% verticesToCounterclockwiseOrder(vertices,referencePoint) returns an array
% of vertices of a polygon in counterclockwise order with respect to a
% reference point, referencePoint, given an array of vertices of the
% polygon. It is assumed that the polygon is closed, that the polygon is
% not self-intersecting or has holes, and that the vertices are not all
% collinear.

% Return an error if less than three vertices are provided
if all(size(vertices)<3)
    error('Three points are needed to define polygon.')
end

if size(vertices,1) < 3
    flip = 1;
    adjustedVertices = vertices';
    adjustedVertices(:,1) = adjustedVertices(:,1) - referencePoint(1);
    adjustedVertices(:,2) = adjustedVertices(:,2) - referencePoint(2);
else
    flip = 0;
    adjustedVertices = vertices;
    adjustedVertices(:,1) = adjustedVertices(:,1) - referencePoint(1);
    adjustedVertices(:,2) = adjustedVertices(:,2) - referencePoint(2);
end

count = 0;
n = size(adjustedVertices,1);

% Sum over the edges between all vertices
for iii = 1:(n-1)
    count = count + (adjustedVertices(iii+1,1)-adjustedVertices(iii,1))*(adjustedVertices(iii+1,2)+adjustedVertices(iii,2));
end
% Also sum the edge from the last to the first point
count = count + (adjustedVertices(1,1)-adjustedVertices(n,1))*(adjustedVertices(1,2)+adjustedVertices(n,2));

% Count is plus/minus twice the enclosed area, with a plus if the vertices
% are in clockwise order relative to the reference point. If so, we flip
% the sequence of the vertices
if count > 0
    adjustedVertices = flipud(adjustedVertices);
end

% Translate back to the original coordinate system
adjustedVertices(:,1) = adjustedVertices(:,1) + referencePoint(1);
adjustedVertices(:,2) = adjustedVertices(:,2) + referencePoint(2);

% If the vertices were transposed, transpose them back
if flip
    counterclockwiseVertices = adjustedVertices';
else
    counterclockwiseVertices = adjustedVertices;
end