https://github.com/wagenadl/sbemalign
Tip revision: d76dcc55e7dad3e7bca91de24d20d201696a5339 authored by Daniel A. Wagenaar on 07 March 2020, 05:53:12 UTC
Cleaned up repo for paper submission
Cleaned up repo for paper submission
Tip revision: d76dcc5
renderq5rigid.py
#!/usr/bin/python3
import aligndb
import sys
import os
import numpy as np
import rawimage
import pyqplot as qp
import factory
import swiftir
import cv2
import config
monttbl = 'solveq5mont'
rigidtbl = 'solveq5rigidtile'
X = Y = 684*5 # Full tile size in q5 space!
nz = 200
db = aligndb.DB()
ri = db.runinfo()
def imagespace(z0):
x0, y0, x1, y1 = db.sel(f'''select
min(mo.x+ri.x), min(mo.y+ri.y),
max(mo.x+ri.x), max(mo.y+ri.y)
from {monttbl} as mo
inner join {rigidtbl} as ri on mo.z0=ri.z0 and mo.r=ri.r and mo.m=ri.m
where mo.z0={z0}''')[0]
x0 = int(x0)
y0 = int(y0)
W = int(x1+X+1)-x0
H = int(y1+Y+1)-y0
return (x0,y0,W,H)
def renderslice(z0, r, s):
# Returns an image and x,y coordinates in model space of its top left
# Positions of tiles
x0, y0 = db.vsel(f'''select x,y from {monttbl}
where z0={z0} and r={r} order by m''')
x1, y1 = db.vsel(f'''select x,y from {rigidtbl}
where z0={z0} and r={r} and s={s} order by m''')
# Full extent of all tiles
xxl = x0 + x1
xxr = xxl + X
yyt = y0 + y1
yyb = yyt + Y
R = ri.nrows(r)
C = ri.ncolumns(r)
(x0,y0,W,H) = imagespace(z0)
img = np.zeros((H,W), dtype=np.uint8) + 255
if C==1:
yjoin = np.zeros(R+1, dtype=int)
yjoin[0] = int(yyt[0]+1)
yjoin[-1] = int(yyb[-1])
for row in range(R-1):
yjoin[row+1] = int((yyb[row] + yyt[row+1])/2 + .5)
for row in range(R):
m = row
print(f' Loading fullq5img {r} {m} {s}')
im1 = rawimage.fullq5img(r, m, s)
if im1 is None:
print(f'Skipping tile R{r} M{m} S{s} - no image')
continue
yt = yjoin[row] # in model coordinates
yb = yjoin[row+1]
xl = int(xxl[m]+1)
xr = xl + X - 2
yt1 = yjoin[row] - yyt[m] # in tile coordinates
yb1 = yt1 + yb-yt
xl1 = xl - xxl[m]
xr1 = xl1 + xr-xl
cen = ((xl1+xr1)/2, (yt1+yb1)/2)
siz = (xr-xl, yb-yt)
im1 = swiftir.extractStraightWindow(im1, cen, siz)
img[yt-y0:yb-y0,xl-x0:xr-x0] = im1
elif C==2:
yjoin = np.zeros(R-1, dtype=int)
xjoin = np.zeros(R, dtype=int)
for row in range(R-1):
yjoin[row] = (yyb[row*C] + yyb[row*C+1]
+ yyt[(row+1)*C] + yyt[(row+1)*C+1]) / 4
for row in range(R):
xjoin[row] = (xxr[row*C] + xxl[row*C+1])/2
for col in range(C):
for row in range(R):
m = col + C*row
print(f' Loading fullq5img {r} {m} {s}')
im1 = rawimage.fullq5img(r, m, s)
if im1 is None:
print(f'Skipping tile R{r} M{m} S{s} - no image')
continue
if row==0:
yt = int(yyt[m]+1)
else:
yt = yjoin[row-1]
if row==R-1:
yb = int(yyb[m])
else:
yb = yjoin[row]
if col==0:
xl = int(xxl[m]+1)
else:
xl = xjoin[row]
if col==C-1:
xr = int(xxr[m])
else:
xr = xjoin[row]
yt1 = yt - yyt[m]
yb1 = yt1 + yb-yt
xl1 = xl - xxl[m]
xr1 = xl1 + xr-xl
cen = ((xl1+xr1)/2, (yt1+yb1)/2)
siz = (xr-xl, yb-yt)
im1 = swiftir.extractStraightWindow(im1, cen, siz)
img[yt-y0:yb-y0,xl-x0:xr-x0] = im1
else:
raise ValueError('Bad column count')
return img, x0, y0
argc = len(sys.argv)
if argc>=2:
z0 = int(sys.argv[1])
else:
z0 = None
if argc>=3:
r = int(sys.argv[2])
else:
r = None
if argc>=4:
s = int(sys.argv[3])
else:
s = None
if s is None:
fac = factory.Factory(12)
def rendersvslice(z0, r, s):
odir = f'{config.root}/rigidq5/Z{z0}'
if not os.path.exists(odir):
os.mkdir(odir)
ofn = odir + f'/R{r}S{s}.jpg'
if not os.path.exists(ofn):
print(f'Rendering R{r}S{s} in z0={z0}')
img, x0, y0 = renderslice(z0, r, s)
rawimage.saveimage(img, ofn)
def rendersvrun(z0, r, s):
if s is None:
sv = ri.subvolume(z0, nz)
for rl in sv:
if rl.r==r:
ss = range(rl.s0, rl.s1)
for s in ss:
fac.request(rendersvslice, z0, r, s)
else:
rendersvslice(z0, r, s)
def rendersubvol(z0, r, s):
if r is None:
sv = ri.subvolume(z0, nz)
for rl in sv:
rendersvrun(z0, rl.r, None)
else:
rendersvrun(z0, r, s)
if z0 is None:
for z0 in range(4,9500,200):
rendersubvol(z0, None, None)
else:
rendersubvol(z0, r, s)
if s is None:
fac.shutdown()