How to discretize a 3D mesh to generate a 3D binary matrix I apply the algorithm proposed by TLousky discretize a mesh to discretize the entire mesh. I get this result

import numpy as np
import bpy, bmesh
from mathutils import Vector
print('Start')
def check_raycast(ray_origin, ray_destination, obj):
    ''' This function casts a virtual ray from "ray_origin" to
        "ray_destination", and finds any intersections along
        the ray's path with the mesh object referenced in the "obj" param.
        If there are any intersections, it will return True, else False.
    '''
    mat = obj.matrix_local.inverted()
    f   = obj.ray_cast(mat * ray_origin, mat * ray_destination)
    loc, normal, face_idx = f

    if face_idx == -1:
        return 0

    return 1
matrix=np.zeros((76,76,88));
#the size of the cube (voxel) is 2*2*2 cm
i=-1
for x in np.arange(-.75,.75,.02):
    i=i+1
    j=-1 
    for y in np.arange(-.75,.75,.02):
        j=j+1;
        k=-1
        for z in np.arange (0,1.75,.02):
            k=k+1
            c = bpy.data.objects['Cube']
            o = bpy.data.objects['modelPerson:Body']
            bpy.context.object.location[0]=x
            bpy.context.object.location[1]=y
            bpy.context.object.location[2]=z

            # Generate bmesh object from cube mesh data
            bm = bmesh.new()
            bm.from_mesh( c.data )
            bm.edges.ensure_lookup_table() # Generates edge   index table
            bm.verts.ensure_lookup_table() # Generates vertex index table
            intersectsMesh = 0
            for e in bm.edges:
                # Find the global coordinates of each edge's two vertices
                coos = [ c.matrix_world * v.co for v in e.verts ]
                # Set these verts as ray casting origin and destination
                ray_origin, ray_destination = coos
                if check_raycast(ray_origin, ray_destination, o):
                    intersectsMesh = 1
                    break
            insideMesh = check_raycast( c.location, Vector( (0,0,1000) ),  o )
            print( "Intersects: ", intersectsMesh )
            print( "Inside Mesh: ", insideMesh )
            matrix[i][j][k]=   insideMesh|intersectsMesh
            print("value",  matrix[i][j][k])
print("max value of matrix before reshape ", np.max(matrix))
matrix2=matrix.reshape(1,16*18*16);
print("max value of matrix ", np.amax(matrix2))
'''I reshape the matrix for using savetxt() after I can reshape the matrix using Mathematica'''
np.savetxt('3DMatrix.csv', matrix.reshape(1,76*76*88),delimiter=',')
print('finish') 






here is my code
You're getting this because of the matrix[i][j][k]=   insideMesh|intersectsMesh line. 

The "insideMesh" test is currently not accurate enough: If the cube intersects with the mesh at a certain Z height, the "insideMesh" test will also be true for all the Z (k) values below that height in that matrix column (that have the same X and Y values), since the test is based on a ray cast straight up towards (0,0,1000).

You can fix this either by replaceing the or | with and &amp; (which might require smaller cubes but nor necessarily), or by making the "insideMesh" test more robust by checking for intersections with all 6 possible directions:

directions = [
    (0,0,1000),
    (0,0,-1000),
    (0,1000,0),
    (0,-1000,0),
    (1000,0,0),
    (-1000,0,0),
]

insideMesh = len( [ i for i in range(6) if check_raycast( c.location, directions[i], o ) ] ) == 6