How to discretize a 3D mesh to generate a 3D binary matrix
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