ShiftP

Buddy of mine wants to straighten some images automatically, kinda like ShiftN :



Oh, you'll want Python version 2, I said:

#makefile
PIP2=pip
PYTHON2=python

setup2:
$(PIP2) install --user virtualenv
$(PYTHON2) -m virtualenv v2
source v2/bin/activate && pip install numpy Pillow pylsd

And Python 3, I said, also in a VirtualENVironment sandbox:

#makefile
PIP3=pip3
PYTHON3=python3

setup3:
#$(PYTHON3) -m venv v3
source v3/bin/activate && pip install numpy Pillow scipy

setup: setup2 setup3
-mkdir temp

run:
source v2/bin/activate && python FindLines.py Source.jpg

source v3/bin/activate && python Warp.py

Start by finding all your lines:

#FindLines.py
import json
import numpy
import sys

from PIL import Image
import pylsd.lsd

def ExportMeta(fileName,outName):
meta={'fileName':fileName}
image=Image.open(fileName)
meta['width']=image.width
meta['height']=image.height
grayScale=numpy.asarray(image.convert('L'))
lines=pylsd.lsd(grayScale)
lineArray=[]
for row in lines:
lineArray.append(list(row))
meta['lineArray']=lineArray

with open(outName,'w') as f:
f.write(json.dumps(meta,sort_keys=True,separators=(',', ': '),indent=4))


ExportMeta(sys.argv[1],'temp/meta.json')

Always prefilter your inputs, I nagged:

#Warp.py
def FindWeightedLines(lineArray,meta):
linesHorizontal=[]
linesVertical=[]
for (x0,y0,x1,y1,width) in lineArray:
if width<2:
continue
h0=RemapXY(x0,y0,meta)
h1=RemapXY(x1,y1,meta)
dx=abs(h0[0]-h1[0])
dy=abs(h0[1]-h1[1])
if max(dx,dy)<min(dx,dy)*4:
continue
magnitude=dx*dx+dy*dy
if dx<dy:
linesHorizontal.append([magnitude,h0,h1])
else:
linesVertical.append([magnitude,h0,h1])

return sorted(linesHorizontal)[-30:]+sorted(linesVertical)[-30:]
#Always prefilter your inputs!!!

Now here's the tricksy bit, in *TWO* parts, setup a perspective transform:

#Warp.py
def RemapXY(x,y,meta):
xx=(x-meta['width']/2)/meta['scale']
yy=(y-meta['height']/2)/meta['scale']
return (xx,yy,1)

def UnmapXYZ(xx,yy,zz,meta):
rx=xx/zz
ry=yy/zz
x=rx*meta['scale']+meta['width']/2
y=ry*meta['scale']+meta['height']/2
return (x,y)

*And* a non-linear warp. We don't need the full power of Chebyshev Polynomials here, I reminded him. We can just use 0, 1, x, y, x², y² and xy. Why? Because all spanning basis-es are equivalent in low dimensions!

#Warp.py
def ApplyTransform(transform,x,y,z):
rx=x*transform[0]+y*transform[1]+z*transform[2]
ry=x*transform[3]+y*transform[4]+z*transform[5]
rz=x*transform[6]+y*transform[7]+z*transform[8]
nonLinear=True
if nonLinear:
rx+=x*y*transform[9]
ry+=x*y*transform[10]
rx+=x*x*transform[11]
ry+=x*x*transform[12]
rx+=y*y*transform[13]
ry+=y*y*transform[14]
return (rx,ry,rz)

def ApplyTransformhomogenous(transform,x,y,z):
(hx,hy,hz)=ApplyTransform(transform,x,y,z)
return (hx/hz,hy/hz)

Next you'll need a loss function, weakly constrain your transform matrix, then setup a sum-of-squares for your error term:

#Warp.py
def loss(transform,meta):

result=0

for i in range(9):

t=transform[i]

if i==0 or i == 4 or i == 8:

t=transform[i]-1

result += t*t

for(weight,h0,h1) in meta['weightedLineArray']:

(x2,y2)=ApplyTransformHomogenous(transform,*h0)

(x3,y3)=ApplyTransformHomogenous(transform,*h1)

dx=abs(x3-x2)

dy=abs(y3-y2)

if dx<dy:
(dx,dy)=(dy,dx)
q=math.sqrt(dx*dx + dy*dy)
dx=dx/q
dy=dy/q

result += dy*dy

Are we done yet? Oh, a driver...:

#Warp.py
def Main():
with open('temp/meta.json','r') as f:
meta=json.loads(f.read())
meta['scale']=math.sqrt(meta['width']*meta['height'])/2

meta['weightedLineArray']=FindWeightedLines(meta['lineArray'],meta)

m=minimize(loss,[1,0,0,0,1,0,0,0,1,0,0,0,0,0,0],args=meta)

transform=m.x

(x0,y0,x1,y1) = FindClipRectangle(...)

source=Image.open(meta['fileName'])

image=Image.new('RGB',(x1-x0,y1-y0),(0,0,0))

draw=ImageDraw.Draw(image)

meta['splitCount']=64

for x in range(meta['splitCount']):

for y in range(meta['splitCount']):

p00=SquareMap(x,y,meta)

p01=SquareMap(x,y+1,meta)

p10=SquareMap(x+1,y,meta)

p11=SquareMap(x+1,y+1,meta)

r00=ApplyTransform(transform,*RemapXY(*p00,meta))

r01=ApplyTransform(transform,*RemapXY(*p01,meta))

r10=ApplyTransform(transform,*RemapXY(*p10,meta))

r11=ApplyTransform(transform,*RemapXY(*p11,meta))

s00=UnmapXYZ(*r00,meta)

s01=UnmapXYZ(*r01,meta)

s10=UnmapXYZ(*r10,meta)

s11=UnmapXYZ(*r11,meta)

TextureMapTriangle(draw,x0,y0,x1,y1,source,s00,s01,s10,p00,p01,p10)

TextureMapTriangle(draw,x0,y0,x1,y1,source,s10,s01,s11,p10,p01,p11)

print('Progress %d/%d'%(x,meta['splitCount']),flush=True)

image.save('temp/warped.jpg')

Oh, and you need texture mapped triangles? Python is terrible for that, there's no way to make it run fast.... Fine, here's one of those, just to get you started, but don't blame me if it's slow, this needs to be in OpenGL or something so you can run it on the GPU and apply proper gamma correction.

#TextureMapTriangle.py
def Left(p,a,b):
cross=(p[0]-a[0])*(b[1]-a[1])-(p[1]-a[1])*(b[0]-a[0])
return cross<0

def SampleMap(source,x,y,dx,dy):
if x<0:
x=0
if y<0:
y=0
if x>=source.width:
x=source.width-1
if y>=source.height:
y=source.height-1
return source.getpixel((x,y))


def TextureMapTriangle(draw,x0,y0,x1,y1,source,p0,p1,p2,uv0,uv1,uv2):
xy0=list(map(min,zip(p0,p1,p2)))
xy1=list(map(max,zip(p0,p1,p2)))

dx1=p1[0]-p0[0]
dy1=p1[1]-p0[1]
dx2=p2[0]-p0[0]
dy2=p2[1]-p0[1]
det=dx1*dy2-dx2*dy1
if xy0[0]<x0:
xy0[0]=x0
if xy0[1]<y0:
xy0[1]=y0
if xy1[0]>x1:
xy1[0]=x1
if xy1[1]>y1:
xy1[1]=y1

for x in range(math.floor(xy0[0]),math.ceil(xy1[0])):
for y in range(math.floor(xy0[1]),math.ceil(xy1[1])):
p=(x,y)
if Left(p,p0,p1):
continue
if Left(p,p1,p2):
continue
if Left(p,p2,p0):
continue

dx=x-p0[0]
dy=y-p0[1]
u=(dx*dy2-dy*dx2)/det
v=(-dx*dy1+dy*dx1)/det

uu=uv0[0]+u*(uv1[0]-uv0[0])+v*(uv2[0]-uv0[0])
vv=uv0[1]+u*(uv1[1]-uv0[1])+v*(uv2[1]-uv0[1])
c=SampleMap(source,uu,vv,1,1)
draw.point((x-x0,y-y0),tuple(c))

And then you could be like me, and license all of the above code under CC0. Yay!