cat > /workspace/scratch/colourfit.py <<'PY'
import json,cv2,numpy as np,tifffile
from scipy.optimize import least_squares
from scipy.interpolate import PchipInterpolator
from PIL import Image
src=json.load(open('/workspace/scratch/sources.json'));boxes=json.load(open('/workspace/scratch/chartboxes.json'));ref=np.array(json.load(open('/workspace/input/chart_reference.json'))['srgb8'])/255
out={}
def lin(v):return np.where(v<=.04045,v/12.92,((v+.055)/1.055)**2.4)
for s in src:
 a=tifffile.imread('/workspace/input/charts/'+s+'_chart.tif').astype('float64')/65535
 x0,y0,x1,y1=boxes[s];step=(x1-x0)/5.87;pw=step*.87;dx=(x1-x0-pw)/5;dy=(y1-y0-pw)/3
 obs=[]
 for row in range(4):
  for col in range(6):
   x=round(x0+pw/2+dx*col);y=round(y0+pw/2+dy*row);r=max(3,round(pw*.25));obs.append(np.median(a[y-r:y+r,x-r:x+r],axis=(0,1)))
 obs=np.array(obs)
 # inverse neutral tone curve: a power law on each delivered channel
 pars=[];corr=[]
 for c in range(3):
  fit=least_squares(lambda p:p[0]*obs[18:,c]**p[1]-ref[18:,c],[1,1],bounds=([.1,.3],[4,3]))
  pars.append(fit.x.tolist())
 cc=np.stack([pars[c][0]*obs[:,c]**pars[c][1] for c in range(3)],1)
 M=np.linalg.lstsq(cc,ref,rcond=None)[0]
 # joint fine fit of exponents and colour matrix
 def fun(p):return (obs**p[:3] @ p[3:].reshape(3,3)-ref).ravel()
 p0=np.r_[np.array(pars)[:,1],(np.diag(np.array(pars)[:,0])@M).ravel()]
 fit=least_squares(fun,p0,max_nfev=300)
 gamma=fit.x[:3];matrix=fit.x[3:].reshape(3,3)
 error=np.abs(fun(fit.x)).max()*255
 anchor=float(lin(obs[21,1])/lin(ref[21,1]))
 print(s,'G22',round(obs[21,1]*255,2),'gamma',gamma.round(4),'max residual',round(error,3),'brightness',round(anchor,4))
 out[s]={'gamma':gamma.tolist(),'matrix':matrix.tolist(),'anchor_linear':anchor,'observed':obs.tolist()}
json.dump(out,open('/workspace/scratch/colourmodels.json','w'),indent=2)
PY
OPENBLAS_NUM_THREADS=1 python /workspace/scratch/colourfit.py