法線から高さを計算するPoisson表面復元のサンプルプログラムです．法線のy軸は上向きにしてください．法線は，各画素に３次元ベクトルを格納したnumpyファイルにしてください．また，物体領域を表す画像ファイルも用意してください．画像中で，物体が存在している領域を白，背景を黒で指定した画像を物体領域ファイルとします．黒の値は0とし，白の値は255とし，画像中に1～254の値が存在ないようにしてください．Google Colabを開いて新規作成し，以下のソースコードのテキスト情報をコピペしてください．実行したら，ファイルアップロードのボタンが出てくるので，ボタンを押して"normal.npy"と"inside.bmp"を指定してください．計算結果のobjファイルがダウンロードされます．
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')
