i=1;  %one test element
p=0;  %initial value for Voxel Nr. inside the element
outside=0;
boundary=0;

for j=1:4
    
Tetras(i,j) = j;

end
     %Coordinates of one Tetrahedron Element with 4 Nodes 
  
     Nodes(Tetras(i,1),1)= 0;   %Tetras(Element Nr, Node Nr),Coordinate X)
     Nodes(Tetras(i,1),2)= 0;   %Tetras(Element Nr, Node Nr),Coordinate Y)
     Nodes(Tetras(i,1),3)= 0;   %Tetras(Element Nr, Node Nr),Coordinate Z)
     Nodes(Tetras(i,2),1)= 1;   %Tetras(Element Nr, Node Nr),Coordinate X)
     Nodes(Tetras(i,2),2)= 0;   %Tetras(Element Nr, Node Nr),Coordinate Y)
     Nodes(Tetras(i,2),3)= 0;   %Tetras(Element Nr, Node Nr),Coordinate Z)
     Nodes(Tetras(i,3),1)= 1;   %Tetras(Element Nr, Node Nr),Coordinate X)
     Nodes(Tetras(i,3),2)= 1;   %Tetras(Element Nr, Node Nr),Coordinate Y)
     Nodes(Tetras(i,3),3)= 0;   %Tetras(Element Nr, Node Nr),Coordinate Z)
     Nodes(Tetras(i,4),1)= 1;   %Tetras(Element Nr, Node Nr),Coordinate X)
     Nodes(Tetras(i,4),2)= 0;   %Tetras(Element Nr, Node Nr),Coordinate Y)
     Nodes(Tetras(i,4),3)= 1;   %Tetras(Element Nr, Node Nr),Coordinate Z)

     V1 = [Nodes(Tetras(i,1),1) Nodes(Tetras(i,1),2) Nodes(Tetras(i,1),3)];   %Vectors including X,Y and Z coordinate of every Node
     V2 = [Nodes(Tetras(i,2),1) Nodes(Tetras(i,2),2) Nodes(Tetras(i,2),3)];
     V3 = [Nodes(Tetras(i,3),1) Nodes(Tetras(i,3),2) Nodes(Tetras(i,3),3)];
     V4 = [Nodes(Tetras(i,4),1) Nodes(Tetras(i,4),2) Nodes(Tetras(i,4),3)];
     
                %Position of the center of one test Voxel (inside or outside the element 
                
                X = 0.5 ;
                Y = 0 ;
                Z = 0 ;
                V5 = [X Y Z];
                
                D0 = [V1 1; V2 1; V3 1; V4 1];    %Matrices containing Node Coordinate Vectors
                D1 = [V5 1; V2 1; V3 1; V4 1];
                D2 = [V1 1; V5 1; V3 1; V4 1];
                D3 = [V1 1; V2 1; V5 1; V4 1];
                D4 = [V1 1; V2 1; V3 1; V5 1];
                
                Dsum = D1 + D2 + D3 + D4;        
               
                DetD0 = det(D0);  %Determinants important for test
                DetD1 = det(D1); 
                DetD2 = det(D2);
                DetD3 = det(D3);
                DetD4 = det(D4);
                
                VD0 = sign(DetD0); %determination of determinants sign
                VD1 = sign(DetD1);
                VD2 = sign(DetD2);
                VD3 = sign(DetD3);
                VD4 = sign(DetD4);
                 
                %dicission about the Voxel Position (inside/outside/border)
                
                if (VD1 == 0) | (VD2 == 0) | (VD3 == 0) | (VD4 == 0)
                    
                    boundary = boundary +1;
                    
                elseif (VD0 == VD1) & (VD0 == VD2) & (VD0 == VD3) & (VD0 == VD4) 
                                        
                    p=p+1;
                    HUwert(p) = 100;
                    
                else
                    
                    outside = outside +1 ;
                end