import json,numpy as np
from scipy.optimize import least_squares
mods=json.load(open('/workspace/scratch/colourmodels.json'));ref=np.array(json.load(open('/workspace/input/chart_reference.json'))['srgb8'])/255
lin=lambda v:np.where(v<=.04045,v/12.92,((v+.055)/1.055)**2.4)
enc=lambda v:np.where(v<=.0031308,12.92*v,1.055*np.maximum(v,0)**(1/2.4)-.055)
out={}
for s,m in mods.items():
 obs=np.array(m['observed']);x=lin(obs);tar=lin(ref)
 p0=np.r_[np.ones(3),np.linalg.lstsq(x,tar,rcond=None)[0].ravel(),np.zeros(3)]
 def fun(p):return (enc(np.maximum(x-p[12:15],1e-10)**p[:3]@p[3:12].reshape(3,3))-ref).ravel()
 fit=least_squares(fun,p0,max_nfev=300)
 corrected=np.maximum(x-fit.x[12:15],1e-10)**fit.x[:3]@fit.x[3:12].reshape(3,3)
 anchor=float(x[21,1]/corrected[21,1])
 print(s,fit.x[:3].round(3),'offset',fit.x[12:].round(4),'mean/max',round(abs(fun(fit.x)).mean()*255,3),round(abs(fun(fit.x)).max()*255,3))
 out[s]={'gamma':fit.x[:3].tolist(),'matrix':fit.x[3:12].reshape(3,3).tolist(),'input_offset':fit.x[12:].tolist(),'anchor':anchor}
json.dump(out,open('/workspace/scratch/finalcolour.json','w'),indent=2)
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)
{
 "chart": "ColorChecker Classic, 24 patches, row-major from the top-left (dark skin) to the bottom-right (black)",
 "values": "BabelColor average sRGB (D65) 8-bit; patches 19-24 are the neutral series",
 "srgb8": [
  [
   115,
   82,
   68
  ],
  [
   194,
   150,
   130
  ],
  [
   98,
   122,
   157
  ],
  [
   87,
   108,
   67
  ],
  [
   133,
   128,
   177
  ],
  [
   103,
   189,
   170
  ],
  [
   214,
   126,
   44
  ],
  [
   80,
   91,
   166
  ],
  [
   193,
   90,
   99
  ],
  [
