Sample program of Poisson surface reconstruction. The y-axis of surface normal should be the upper direction. The surface normal should be the numpy file with 3D vector in each pixel. The image file of object region is also needed. The object region should be white, and the background should be black in this file. Black should be 0 and white should be 255, and 1-254 should not be used. Open Google Colab, create new file, copy&paste the text below. Run, push upload button, select "normal.npy" and "inside.bmp". Output obj faile is downloaded.
from google.colab import files
import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
import cv2
files.upload()
normal=np.load('normal.npy')
inside=cv2.imread('inside.bmp',cv2.IMREAD_GRAYSCALE)
inside=inside.astype(dtype=bool)
rows,cols,_=normal.shape
print('Gradient...',end='',flush=True)
GRADIENTMAX=40
gradp=np.zeros((rows,cols),dtype=np.float64)
gradq=np.zeros((rows,cols),dtype=np.float64)
for y in range(0,rows):
for x in range(0,cols):
nx=normal[y,x,0]
ny=normal[y,x,1]
nz=normal[y,x,2]
if nz<1e-15:
nz=1e-15
gradp[y,x]=np.clip(-nx/nz,-GRADIENTMAX,GRADIENTMAX)
gradq[y,x]=np.clip(-ny/nz,-GRADIENTMAX,GRADIENTMAX)
gradq[y,x]=-gradq[y,x]
print('done')
print('Matrix...',end='',flush=True)
size=rows*cols
left=lil_matrix((size,size),dtype=np.float64)
right=lil_matrix((size,1),dtype=np.float64)
for y in range(0,rows):
for x in range(0,cols):
i=y*cols+x
iy0=(y-1)*cols+x
iy1=(y+1)*cols+x
ix0=y*cols+(x-1)
ix1=y*cols+(x+1)
if inside[y,x]:
w=0
g=0
if y-1>=0:
if inside[y-1,x]:
left[i,iy0]=-1
w+=1
g+=(gradq[y-1,x]+gradq[y,x])/2
if x-1>=0:
if inside[y,x-1]:
left[i,ix0]=-1
w+=1
g+=(gradp[y,x-1]+gradp[y,x])/2
if y+1<rows:
if inside[y+1,x]:
left[i,iy1]=-1
w+=1
g+=-(gradq[y+1,x]+gradq[y,x])/2
if x+1<cols:
if inside[y,x+1]:
left[i,ix1]=-1
w+=1
g+=-(gradp[y,x+1]+gradp[y,x])/2
left[i,i]=w
right[i,0]=g
left[i,i]+=1e-15*1
right[i,0]+=1e-15*0
else:
left[i,i]=1
right[i,0]=0
print('done')
print('Integrate...',end='',flush=True)
height=spsolve(left.tocsr(),right.tocsr())
height=height.reshape((rows,cols))
height-=np.min(height)
print('done')
print('Save...',end='',flush=True)
OBJSCALE=0.01
vertex=np.zeros((rows,cols),dtype=np.int32)
with open('shape.obj',mode='w') as f:
f.write('o shape\n')
num = 1
for y in range(rows-1,-1,-1):
for x in range(0,cols):
if inside[y,x]:
dx=x-cols/2
dy=y-rows/2
dy=-dy
dz=height[y,x]
dx*=OBJSCALE
dy*=OBJSCALE
dz*=OBJSCALE
f.write('v %1.7f %1.7f
%1.7f\n'%(dx,dy,dz))
vertex[y,x]=num
num+=1
else:
vertex[y,x]=0
for y in range(rows-1,0,-1):
for x in range(0,cols-1):
i00=vertex[y,x]
i10=vertex[y,x+1]
i11=vertex[y-1,x+1]
i01=vertex[y-1,x]
if i00>0 and i11>0 and i01>0:
f.write('f %d %d %d\n'%(i00,i11,i01))
if i11>0 and i00>0 and i10>0:
f.write('f %d %d %d\n'%(i11,i00,i10))
print('done')
files.download('shape.obj')
