deconv_rrrl2

#!/usr/bin/env python3

# Copyright (C) 2017 Vladimir Nadvornik
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.


import tifffile
import numpy as np
import cv2

from ctypes import *
cv_lib = CDLL("libopencv_core.so.4.2.0")  
print(cv_lib._ZN2cv5utils7logging11setLogLevelENS1_8LogLevelE(6))
print(cv_lib._ZN2cv5utils7logging11getLogLevelEv(6))


import sys
import os
import argparse
import pfor
import gc
import pickle
import time

from astro_utils import noise_level, poly_bg, poly_res, poly_array,  extrapolate_transp
from centroid import centroid, sym_center
from detect_lines import detect_lines, draw_lines
from scipy.optimize import curve_fit

cv2.setNumThreads(-1)

def crop_param(s):
	try:
		if s is None:
			return None
		x1,x2,y1,y2 = s.split(',')
		return (slice(int(y1),int(y2)), slice(int(x1),int(x2)))
	except:
    		raise argparse.ArgumentTypeError("Filter param must be r,s")


parser = argparse.ArgumentParser()
parser.add_argument("outfile",
                    help="output tiff file")
parser.add_argument("infile", nargs='+',
                    help="input tiff file")

parser.add_argument("--psffile",
                    help="psf tiff file if different from infile")

parser.add_argument("--iter", type=int, default=10,
                    help="number of iteration")
parser.add_argument("--diameter", type=int, default=15,
                    help="psf diameter")
parser.add_argument("--tiles", type=int, default=4,
                    help="number of tiles NxN")
parser.add_argument('--var-psf', action='store_true',
                    help="variable psf")

parser.add_argument("--lp-depth-min", type=int, default=2,
                    help="min depth for setup")

parser.add_argument('--update-scales', action='store_true',
                    help="update scales")

parser.add_argument("--psf-filter-sigma-res", type=float, default=1.2,
                    help="psf filter")
parser.add_argument("--psf-filter-sigma", type=float, default=0.5,
                    help="psf filter final sigma")
parser.add_argument("--psf-filter-moffat-scale", type=float, default=1,
                    help="psf filter moffat scale")
parser.add_argument("--psf-filter-order", type=int, default=3,
                    help="psf filter order")

parser.add_argument("--update-iter", type=int, default=0,
                    help="start updating at iteration i")
parser.add_argument("--update-accel", type=float, default=20,
                    help="update accel")
parser.add_argument("--update-thr", type=float, default=0,
                    help="update thr")
parser.add_argument("--update-psf-robustness", type=float, default=2,
                    help="update thr")

parser.add_argument("--blur-iter", type=int, default=0,
                    help="start updating at iteration i")
parser.add_argument("--blur-sigma", type=float, default=2,
                    help="update accel")

parser.add_argument("--reg2", type=float, default=[0.1], nargs = '+',
                    help="regularization2")
parser.add_argument("--reg-lambda2", type=float, default=[2.0], nargs = '+',
                    help="lambda2")
parser.add_argument("--reg4", type=float, default=[0.1], nargs = '+',
                    help="regularization4")
parser.add_argument("--reg-lambda4", type=float, default=[2.0], nargs = '+',
                    help="lambda4")
parser.add_argument("--regf", type=float, default=[0.1], nargs = '+',
                    help="regularization flatten")
parser.add_argument("--reg-lambdaf", type=float, default=[2.0], nargs = '+',
                    help="lambda f")
parser.add_argument("--lambda-eval-iter", type=int, default=1000,
                    help="stop updating lambda at iteration i")
parser.add_argument("--reg-plus", type=float, default=1,
                    help="reg plus")

parser.add_argument("--reg-mode", default='openclose',
                    help="openclose open close dilate erode median normal")


parser.add_argument("--reg-morph-radius", type=int, default=2,
                    help="morph radius")

parser.add_argument("--top-thres", type=float, default=0.5,
                    help="top-thres")

parser.add_argument("--robustness", type=float, default=2,
                    help="robustness")

parser.add_argument("--cut-out", type=float, default=1,
                    help="robustness")

parser.add_argument("--dering", type=float, default=0.9,
                    help="deringing")
parser.add_argument("--dering2", type=float, default=1,
                    help="deringing")
parser.add_argument("--dering2-radius", type=int, default=13,
                    help="deringing")
parser.add_argument("--dering2-thr", type=float, default=3,
                    help="deringing")


parser.add_argument("--gaussian", type=float, default=0.0,
                    help="blind with gaussian")

parser.add_argument("--gradient-radius", type=int, default=300,
                    help="gradient radius")

parser.add_argument("--crop", type=crop_param, default=None,
                    help="crop x1,x2,y1,y2")

parser.add_argument("--scale", type=float, default=1,
                    help="output scale")

parser.add_argument("--zero", type=float, default=1024,
                    help="zero level")

parser.add_argument("--gamma", type=float, default=1.0,
                    help="process gamma")

parser.add_argument("--wb", type=float, default=[1.0], nargs = '+',
                    help="white balance")

parser.add_argument("--channel-weights", type=float, default=[], nargs = '+',
                    help="channel weights")

parser.add_argument("--show-gamma", type=float, default=0.1,
                    help="display gamma")

parser.add_argument("--overexp", type=float, default=0.8,
                    help="overexp level")

parser.add_argument("--tile-overlap", type=float, default=0.25,
                    help="tile overlap")

parser.add_argument("--weight-dilate", type=int, default=2,
                    help="weight dilate")

parser.add_argument("--delta-x", type=float, default=0,
                    help="delta x")
parser.add_argument("--delta-y", type=float, default=0,
                    help="delta y")


parser.add_argument("--save-psf",
                    help="save psf")
parser.add_argument("--load-psf",
                    help="load psf")

parser.add_argument("--test-ptlist", action='store_true',
                    help="test ptlist")


parser.add_argument('--poly-bg', type=int, default = 0,
                   help='poly_bg_order')
parser.add_argument('--poly-bg-iter', type=int, default = 10,
                   help='poly_bg_iter')
parser.add_argument('--poly-bg-kappa', type=float, default = 2,
                   help='poly_bg_kappa')

parser.add_argument('--stars-per-tile', type=int, default = 4,
                   help='min number of stars per tile')

args = parser.parse_args()

while len(args.reg2) < len(args.reg_lambda2):
	args.reg2.append(args.reg2[-1])
while len(args.reg_lambda2) < len(args.reg2):
	args.reg_lambda2.append(args.reg_lambda2[-1])
while len(args.reg4) < len(args.reg_lambda4):
	args.reg4.append(args.reg4[-1])
while len(args.reg_lambda4) < len(args.reg4):
	args.reg_lambda4.append(args.reg_lambda4[-1])
while len(args.regf) < len(args.reg_lambdaf):
	args.regf.append(args.regf[-1])
while len(args.reg_lambdaf) < len(args.regf):
	args.reg_lambdaf.append(args.reg_lambdaf[-1])


def normalize(img):
        dst = np.empty_like(img)
        return cv2.normalize(img, dst, alpha = 0, beta = 255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8UC1)


hfr_size = None
hfr_mat_mask = None
hfr_mat = None
xsize_mat = None
ysize_mat = None
d1size_mat = None
d2size_mat = None

def set_htr_size(diameter):
	global hfr_size
	global hfr_mat_mask
	global hfr_mat
	global xsize_mat
	global ysize_mat
	global d1size_mat
	global d2size_mat
	
	hfr_size = diameter
	hfr_mat_mask = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (hfr_size * 2 + 1, hfr_size * 2 + 1))
	hfr_mat = cv2.multiply(np.array([[(x**2 + y**2)**0.5 for x in range(-hfr_size, hfr_size + 1) ] for y in range(-hfr_size, hfr_size + 1) ], dtype=np.float), hfr_mat_mask, dtype=cv2.CV_32FC1)
	xsize_mat = cv2.multiply(np.array([[ abs(x) for x in range(-hfr_size, hfr_size + 1) ] for y in range(-hfr_size, hfr_size + 1) ], dtype=np.float), hfr_mat_mask, dtype=cv2.CV_32FC1)
	ysize_mat = cv2.multiply(np.array([[ abs(y) for x in range(-hfr_size, hfr_size + 1) ] for y in range(-hfr_size, hfr_size + 1) ], dtype=np.float), hfr_mat_mask, dtype=cv2.CV_32FC1)
	d1size_mat = cv2.multiply(np.array([[ abs(x + y) for x in range(-hfr_size, hfr_size + 1) ] for y in range(-hfr_size, hfr_size + 1) ], dtype=np.float), hfr_mat_mask, dtype=cv2.CV_32FC1)
	d2size_mat = cv2.multiply(np.array([[ abs(x - y) for x in range(-hfr_size, hfr_size + 1) ] for y in range(-hfr_size, hfr_size + 1) ], dtype=np.float), hfr_mat_mask, dtype=cv2.CV_32FC1)

def hfr(a , mat = None):
	if mat is None:
		mat = hfr_mat
	s = cv2.sumElems(cv2.multiply(a,  hfr_mat_mask, dtype=cv2.CV_32FC1))[0]
	if s == 0.0:
		return hfr_size
	r = cv2.sumElems(cv2.multiply(a,  mat, dtype=cv2.CV_32FC1))[0] / s
	return r

def poly_array2(X, Y):
	res = np.empty([X.shape[0], 4])
	res[:, 0] = 1
	res[:, 1] = X ** 2
	res[:, 2] = Y ** 2
	res[:, 3] = Y * X

	return res



def find_ptlist(img, zero = 0.0):
	kernel = hfr_mat_mask #np.ones((args.diameter,args.diameter),np.uint8)
	img = np.array(img, dtype=np.float32, copy=True)
	bl = cv2.GaussianBlur(img,(5, 5),0)
	dil = cv2.dilate(bl, kernel)

	cmpmax = cv2.compare(bl, dil, cv2.CMP_GE)
	
	

	#er = cv2.erode(bl, kernel)
	#bg = cv2.GaussianBlur(er,(19, 19),0)
	#img -= bg
	#img[np.where(img < 0)] = 0

	sigma = noise_level(img)
	cmpmax[np.where(bl <= zero + sigma * 10)] = 0
	
	cv2t.imshow("cmpmax", normalize(cmpmax))

	ptlist = []
	for (y, x) in zip(*cmpmax.nonzero()):
		psf = cv2.getRectSubPix(img, (args.diameter * 2 + 1, args.diameter * 2 + 1), (x, y))
		
		(xc, yc) = centroid(psf)
		ptlist.append((y + yc, x + xc))

	return ptlist

def psf_bg_w(psf):

	bg = np.median(psf)
	mask = np.zeros_like(psf, dtype=np.uint8)
	mask[psf <= bg] = 255

	for i in range(10):
		bg, stddev = cv2.meanStdDev(psf, mask = mask)
		mask = np.zeros_like(psf, dtype=np.uint8)
		mask[psf < bg + 2 * stddev] = 255

    
	psf2  = psf - bg
	w = np.sum(psf2)
	
	return bg, w


def fix_hfr(v, hfr, it):
	A0 = np.ones((len(v), 3))
	A0[:,0] = v
	A0[:,2] = 1.0 / v
	W0 = v ** 0.5
	hfr0 = hfr
	A = A0
	W = W0
	print()
	for i in range(0, it):
		c = np.linalg.lstsq(A * W[:, np.newaxis], hfr * W)[0]
		print(c)
		fit =  np.dot(A, c)
		d2 = (hfr - fit) ** 2
		s2 = np.average(d2)
		
		keep=np.where((d2 < s2 * 10) | (hfr < fit))
		A = A0[keep]
		hfr = hfr0[keep]
		W = W0[keep]
	return c[0:2]

def filter_ptlist(img, i_ptlist, debug=False):
	minv = 64 / 65535.0

	maxval = np.max(img)

	(height,width) = img.shape
	hfrlist = []
	ptlist = []

	for (y, x) in i_ptlist:
		psf = cv2.getRectSubPix(img, (args.diameter * 2 + 1, args.diameter * 2 + 1), (x, y))
		
		bl = cv2.GaussianBlur(cv2.medianBlur(psf, 5), (7, 7), 0)
		minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(bl)
		
		if maxVal > maxval * 0.9:
			continue

		(xc, yc) = sym_center(psf)
		#(xc, yc) = (0, 0)
		
		#bg, w = psf_bg_w(psf)

		#if maxVal - bg < minv:
		#	continue

		psf2 = cv2.getRectSubPix(img, (args.diameter * 2 + 1, args.diameter * 2 + 1), (x + xc, y + yc))
		(xc2, yc2) = sym_center(psf2)
		xc += xc2
		yc += yc2
		psf2 = cv2.getRectSubPix(img, (args.diameter * 2 + 1, args.diameter * 2 + 1), (x + xc, y + yc))
		
		bg, w = psf_bg_w(psf2)
		if w <= 0:
			continue
		
		psf2 -= bg
		#psf2 = np.abs(psf2)
		
		hf = []
		skip = False
		for mat in [hfr_mat, xsize_mat, ysize_mat, d1size_mat, d2size_mat]:
			h = hfr(psf2, mat)
			if h <= 0:
				skip = True
				break
			hf.append(h)
		if skip or hf[0] > args.diameter * 0.5:
			continue
		hf.append(xc)
		hf.append(yc)
		#w = max(0, maxVal - bg)

		hfrlist.append([y, x, w] + hf)
		ptlist.append((y, x))

	hfrlist = np.array(hfrlist)
	
	tile_n = np.empty((len(hfrlist),), dtype=np.int)
	tiles = args.tiles ** 2
	tile_h = float(height) / args.tiles
	tile_w = float(width) / args.tiles


	if debug:
		show = normalize(img ** 0.4)
		for (y, x) in hfrlist[:, 0:2]:
			cv2.circle(show, (int(x),int(y)), 10, (20), 1)

	
	for i in range(0, args.tiles):
		for j in range(0, args.tiles):
			tile_n[(hfrlist[:, 0] >= i * tile_h) & (hfrlist[:, 0] < (i + 1) * tile_h) 
			     & (hfrlist[:, 1] >= j * tile_w) & (hfrlist[:, 1] < (j + 1) * tile_w)] = i * args.tiles + j 
	
	#sigma = 10
	keep_all = None

	for t in range(0, tiles):
		
		tile_mask = (tile_n == t)
		if not np.any(tile_mask):
			continue
		avg = np.median(hfrlist[tile_mask,2])
		print("tile avg", t, np.log(avg))

		keep = (hfrlist[:,2] >= avg) & tile_mask
		if not np.any(tile_mask):
			continue
		
		for i in range(0, 10):
			med = np.average(hfrlist[keep, 3], weights = hfrlist[keep, 2])
			sigma2 = np.average((hfrlist[keep, 3] - med)**2, weights = hfrlist[keep, 2])
			sigma2 = max(sigma2, 0.1)
			print("m s", med, sigma2)
			keep = (hfrlist[:, 3] <= med + sigma2**0.5 * 3)  & tile_mask
		#if sigma2 > 6:
		#	continue

		if keep_all is None:
			keep_all = keep
		else:
			keep_all = keep_all | keep

	keep_all = np.where(keep_all)
	hfrlist = hfrlist[keep_all]
	ptlist = [ptlist[i]  for i in keep_all[0]]


	if debug:
		print("plot")
		import matplotlib.pyplot as plt
		
		p = hfrlist#[hfrlist[:, 2] > 20]
		x = p[:, 2]
		y = p[:, 3]# / p[:, 3]
		
		xm = np.amax(x)
		for i in range(100, 100):
			m, c = fix_hfr(x, y, i)
			print(m,c, xm)
			plt.plot([0, xm], [c, c + xm * m], 'b-')

		plt.plot(np.log(x), y, 'ro')
		plt.show()
#	print(hfrlist)
#	print(ptlist)

	Y = hfrlist[:,0]
	X = hfrlist[:,1]
	
	achar0 = np.array(hfrlist[:, 3 :])
	
	#achar0[:, 0 : 4] /= hfrlist[:, 3][:, np.newaxis]
	
	Y = (Y / height * 2.0) - 1.0
	X = (X / width * 2.0) - 1.0
	
	A0 = poly_array(X, Y, args.psf_filter_order)
	
	W = np.array(hfrlist[:,2])


	W /= np.amax(W)
	W = W * (1.0 - W)
	W **= 0.5
	
	Ws = np.sqrt(W)
	
	A0w = A0 * Ws[:,np.newaxis]
	
	achar0w = achar0 * Ws[:,np.newaxis]
	
	achar = achar0
	A = A0
	
	acharw = achar0w
	Aw = A0w
	
	nchar = len(achar0[0])
	
	char_sigma2 = [10000.0] * nchar
	char_coef = [None] * nchar
	char_fit0 = [None] * nchar
	

	tile_n0 = np.empty((len(hfrlist),), dtype=np.int)
	for i in range(0, args.tiles):
		for j in range(0, args.tiles):
			tile_n0[(hfrlist[:, 0] >= i * tile_h) & (hfrlist[:, 0] < (i + 1) * tile_h) 
			     & (hfrlist[:, 1] >= j * tile_w) & (hfrlist[:, 1] < (j + 1) * tile_w)] = i * args.tiles + j 
	
	tile_n = tile_n0

	for i in range(0, 30):
		for j in range(0, nchar):
			char_coef[j] = np.linalg.lstsq(Aw, acharw[:, j])[0]
			char_fit = np.dot(A, char_coef[j])
		
			#hfr_sigma2 = np.average((ahfr - hfr_fit)**2, weights = WW)
			
			kappa = 4

			char_fit0[j] = np.dot(A0, char_coef[j])
			
			keep = None
			for t in range(0, tiles):
		
				tile_mask0 = (tile_n0 == t)
				if not np.any(tile_mask0):
					continue
				
				tile_mask = (tile_n == t)
			
				char_sigma2 = max(np.average((achar[tile_mask, j] - char_fit[tile_mask])**2), 0.1)
				print("i{} j{} t{} sigma2 {}".format(i, j, t, char_sigma2))
			
				t_keep = ((achar0[:, j] - char_fit0[j])**2 <= char_sigma2 * kappa) & tile_mask0
				
				if keep is None:
					keep = t_keep
				else:
					keep = keep | t_keep
				
				
				
			if j == 0:
				gkeep = keep
			else:
				gkeep = gkeep & keep
		
		gkeep = np.where(gkeep)

		#if i == 20:
		#	A0 = poly_array(X, Y, 5)
		#	A0w = A0 * Ws[:,np.newaxis]
		
		A = A0[gkeep]
		achar = achar0[gkeep]

		Aw = A0w[gkeep]
		acharw = achar0w[gkeep]
		
		tile_n = np.empty((len(hfrlist[gkeep]),), dtype=np.int)
		for i in range(0, args.tiles):
			for j in range(0, args.tiles):
				idx = ((hfrlist[gkeep[0], 0] >= i * tile_h) & (hfrlist[gkeep[0], 0] < (i + 1) * tile_h) 
				     & (hfrlist[gkeep[0], 1] >= j * tile_w) & (hfrlist[gkeep[0], 1] < (j + 1) * tile_w))
				
				tile_n[idx] = i * args.tiles + j 

		
		print('len', len(achar), char_coef)
	

	if debug:
		print("plot")
		import matplotlib.pyplot as plt
		
		p = hfrlist[gkeep]#[hfrlist[:, 2] > 20]
		x = p[:, 2]
		y = p[:, 3]# / p[:, 3]
		
		xm = np.amax(x)
		for i in range(100, 100):
			m, c = fix_hfr(x, y, i)
			print(m,c, xm)
			plt.plot([0, xm], [c, c + xm * m], 'b-')

		plt.plot(np.log(x), y, 'ro')
		plt.show()

	if debug:
		for (y, x) in ptlist:
			cv2.circle(show, (int(x),int(y)), 12, (100), 1)

	print(gkeep)
	ptlist = [ptlist[i]  for i in gkeep[0]]


	if debug:
		for (y, x) in ptlist:
			cv2.circle(show, (int(x),int(y)), 14, (255), 1)
		cv2.imwrite("pts_t.tif", show)

	return ptlist
	

def get_psf_list(img, ptlist, over):
	img = np.array(img, dtype=np.float32)
	
	psflist = []
	for (y, x) in ptlist:
		psf = cv2.getRectSubPix(img, (args.diameter * 2 + 1, args.diameter * 2 + 1), (x, y))
		if np.amax(psf) > over:
			continue

		psflist.append(psf)
	return psflist


def rl_simple(src, psf, it):
	src = np.array(src, dtype=np.float64)
	res = np.ones_like(src)
	for i in range(0, it):
		b_res = psf.apply(res)
		cor = cv2.divide(src, b_res)
		cor[np.where(cor < 0.00001)] = 0.00001
		cor = psf.apply(cor, flip = True)
		res *= cor
	return res


def rrrl_simple(src, psf, it):
	psf_zero = 0.001
	scales = np.ones((len(src),), dtype=np.float64)
	weights = np.ones((len(src),), dtype=np.float64)
	for i in range(len(src)):
		bg, w = psf_bg_w(src[i])
		#bl = cv2.GaussianBlur(cv2.medianBlur(src[i], 5), (7, 7), 0)
		#minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(bl)
		
		#src[i] -= minVal - psf_zero
		#scales[i] = 1 / (maxVal + 0.00000000001 - minVal) 
		#weights[i] = (maxVal - minVal) #** 2
		
		if w < 0:
			w = 0
		src[i] -= bg - psf_zero
		scales[i] = 1 / (w + 0.00000000001)
		weights[i] = w #** 2
		
		
	src = np.atleast_3d(src)
	src[np.where(src < 0)] = 0
	
	res = np.mean(src, axis = 0)
	res = cv2.GaussianBlur(res, (5,5), 0)
	
	
	for i in range(0, it):
		kappa = (it - i) / it * 10 + 3
		b_res = psf.apply(res)
		b_res[np.where(b_res < 0.0000001)] = 0.0000001
		
		sum_w = np.zeros_like(res)
		sum_corw = np.zeros_like(res)
		
		scaled_src = src * scales[:, None, None] + psf_zero * (1.0 - scales[:, None, None])
		
		cor = scaled_src / b_res[None, :, :]
		cor[np.where(cor < 0.00001)] = 0.00001
		
		dat = b_res[None, :, :] - scaled_src# * (1.0 - np.log(cor))
		dat = dat * dat
		
		var = np.mean(dat, axis = 0)
		clip = np.where(dat > (var * kappa**2)[None, :, :])
		
		var = np.mean(var)
		
		dat = (dat + var) ** -0.25 * weights[:, None, None]
		
		dat[clip] = 0.0
		
		for c in range(0, len(src)):
			#cor_mean = np.average(cor[c], weights = dat[c] + 0.0000001)
			#scales[c] /= cor_mean

			corw = psf.apply(cv2.multiply(dat[c], cor[c]), flip = True)
			w =  psf.apply(dat[c], flip = True)

			corw[np.where(corw < 0.0)] = 0.0
			w[np.where(w < 0)] = 0

			sum_w += w
			sum_corw += corw
		
		sum_w[np.where(sum_w < 1e-8)] = 1e-8
		cor = cv2.divide(sum_corw, sum_w)
		
		bl = cv2.GaussianBlur(res, (3, 3), 0)
		lap = res - bl
		g = (1.0 - np.clip(res / 0.1, 0, 1)) ** 2
		lap *= g
		lap[res < psf_zero] = (res - psf_zero)[res < psf_zero]
		
		res *= cor
		res -= lap
		
	#print(scales)
	res -= psf_zero
	res[np.where(res < 0)] = 0
	
	for c in range(0,len(src)):
		weights[c] = np.mean(dat[c])
	
	target = np.average(scaled_src, axis = 0, weights = weights)
	
	return res, target
		

optimal_dft_size_cache = {}
def optimal_dft_size(s):
	if s in optimal_dft_size_cache:
		return optimal_dft_size_cache[s]
	r = int(cv2.getOptimalDFTSize(s))
	while r in [1875, 2025, 2160, 2304, 2400, 2430, 2500, 2560, 2592, 2700, 2880]:
		r = int(cv2.getOptimalDFTSize(r + 1))
	optimal_dft_size_cache[s] = r
	print("optimal_dft_size_cache %d: %d" % (s, r))
	return r

class PsfGauss:
	def __init__(self, size, sigma):
		self.size = size
		self.sigma = sigma
	
	def apply(self, img, flip = False):
		return cv2.GaussianBlur(img, (self.size, self.size), self.sigma)

class PsfMoffat:
	def __init__(self, size, alpha = 1, beta=4.765, psf = None, residual = None):
		self.size = size
		self.mx = (size - 1) / 2.0
		self.my = (size - 1) / 2.0
		self.shape = [size, size]
		
		if psf is not None:
			self.fit(psf, residual)
		else:
			self.set(alpha, beta)

	def show(self, name, psf, residual, moffat):
		h = 4
		w = len(psf)
		img = np.zeros((h * self.size, w * self.size), dtype = np.uint8)
		def show1(x, y, i1):
			img[y * self.size : (y+1) * self.size, x * self.size : (x+1) * self.size] = normalize(i1)
		
		for i in range(0, w):
			show1(i, 0, residual[i] ** 0.4)
			show1(i, 1, psf[i] ** 0.4)
			show1(i, 2, moffat[i] ** 0.4)
			show1(i, 3, psf[i] - moffat[i]) 
		img = cv2.resize(img, (img.shape[1] * 2, img.shape[0] * 2), interpolation=cv2.INTER_NEAREST)
		cv2t.imshow(name, img)

	def fit(self, psf, residual = None):
		r2 = np.array([[ (x - self.mx) ** 2 + (y - self.my) ** 2 for x in range(0, self.size)] for y in range(0, self.size)])
		psf = np.array(psf)
		
		A = np.ones((psf.shape[1] * psf.shape[2], 2))
		
		def moffat2d(syx, alpha, beta):
			#s, y, x = syx
			#print(alpha, beta, mag, shift)
			moffat = (1 + (r2 / (alpha * alpha))) ** -beta
			ret = np.empty_like(psf)
			for i in range(0, len(psf)):
				#if residual is not None:
				moffat_i = cv2.filter2D(residual[i], -1, moffat)
				A[:, 0] = moffat_i.ravel()
				m, c = np.linalg.lstsq(A, psf[i].ravel())[0]
				ret[i] = c + m * moffat_i
			#print(ret)
			ret *= 1000.0
			return ret.ravel()#[(s, y, x)]

		slices = len(psf)

		x0 = np.arange(self.size)
		y0 = np.arange(self.size)
		s0 = np.arange(slices)
		s, y, x = np.meshgrid(s0, y0, x0, indexing='ij')
		syx = (s.ravel(), y.ravel(), x.ravel())
		s, y, x = syx
		z = psf[syx] * 1000.0

		min_alpha = (args.psf_filter_sigma - args.psf_filter_sigma_res)
		
		try:
			popt, pcov = curve_fit(moffat2d, syx, z, p0 = [3, 5], bounds = ([min_alpha, 1], [6.0, 10]), loss='soft_l1', verbose = 1, xtol=1e-12, ftol=1e-12, gtol=1e-12, jac = '3-point')
			print ("fit1a", popt)
			self.show("fit1a", psf * 1000.0, residual, moffat2d(syx, *popt).reshape(psf.shape))
			
			for i in range(0,3):
				sigma = z - moffat2d(syx, *popt)
				ss = np.median(sigma ** 2) ** 0.5 * 0.1
				sigma[sigma < ss] = ss
				sigma *= 2
				popt, pcov = curve_fit(moffat2d, syx, z, p0 = popt, bounds = ([0.01, 1], [8.0, 10]), loss='soft_l1', verbose = 1, xtol=1e-12, ftol=1e-12, gtol=1e-12, jac = '3-point', sigma = sigma)
			
			print ("fit1b", popt)
			
			self.show("fit1b", psf * 1000.0, residual, moffat2d(syx, *popt).reshape(psf.shape))
		except  Exception as e:
			print(e)
			popt = [1.0, 4.765]
		self.set(popt[0], popt[1])




	def set(self, alpha, beta):
		self.alpha = alpha
		self.beta = beta
		
		r2 = np.array([[((x - self.mx) **2 + (y - self.my)**2) for x in range(self.size) ] for y in range(self.size) ], dtype=np.float)
		
		self.psf = (1 + r2 / (alpha **2 )) ** -beta

		cv2t.imshow("moffat", normalize(self.psf ** 0.4))
	
		self.psf /= self.psf.sum()
		self.psf = np.array(self.psf, dtype=np.float32)
		
	def interpolate(self, s):
		if s < 0.01:
			s = 0.01
		return PsfMoffat(self.size, self.alpha * s, self.beta)
	
	def apply(self, img, flip = False, img_shape = None):
		if img_shape is None:
			img_shape = img.shape
			to_mat = True
		else:
			to_mat = False
			
		dft_shape = (optimal_dft_size(img_shape[0] + self.shape[0]), optimal_dft_size(img_shape[1] + self.shape[1]))

		eimg = cv2.copyMakeBorder(cv2.UMat(img), top=0, bottom=dft_shape[0]-img_shape[0], left=0, right=dft_shape[1]-img_shape[1], borderType=cv2.BORDER_CONSTANT, value=0)
		epsf = cv2.copyMakeBorder(cv2.UMat(self.psf), top=0, bottom=dft_shape[0]-self.shape[0], left=0, right=dft_shape[1]-self.shape[1], borderType=cv2.BORDER_CONSTANT, value=0)

		IMG = cv2.dft(eimg, flags=cv2.DFT_COMPLEX_OUTPUT, nonzeroRows = img_shape[0])
		PSF = cv2.dft(epsf, flags=cv2.DFT_COMPLEX_OUTPUT, nonzeroRows = self.shape[0])
		RES = cv2.mulSpectrums(IMG, PSF, 0)

		res = cv2.idft(RES, flags=cv2.DFT_SCALE | cv2.DFT_REAL_OUTPUT )
		
		anchor = (self.shape[1] // 2, self.shape[0] // 2)
		res = cv2.UMat(res, (anchor[0], anchor[0] + img_shape[0]), (anchor[1], anchor[1] + img_shape[1]))
		if to_mat:
			return res.get()
		return res

def lowblur(img, thr, n):
	for i in range(n):
		bl = cv2.GaussianBlur(img, (3, 3), 0)
		lap = img - bl
		g = (1.0 - np.clip(img / thr, 0, 1)) ** 2
		lap *= g
#		lap -= cv2.GaussianBlur(lap, (3, 3), 0)
		img -= lap
	return img

class Psf:
	def set(self, psf):
		psf = np.array(psf, dtype=np.float32)
		psf[np.where(psf < 0)] = 0
		psf /= psf.sum()
		self.psf = psf
		self.psf_flip = self.psf[::-1, ::-1]
		self.updated = False
		self.center = (0.0, 0.0)
		self.shape = psf.shape
		

	def extract(self, img, ptlist, psf_filter_sigma):
		psflist = get_psf_list(img, ptlist)
		self.from_psflist(psflist, psf_filter_sigma)

	def from_psflist(self, psflist, psf_filter_sigma):
		if len(psflist) == 0:
			return self.gaussian(None, 1.0)
		#psf = rrrl_simple(psflist, PsfMoffat(args.diameter * 2 + 1, 3, 2), 30)
		psf, target = rrrl_simple(psflist, PsfGauss(args.diameter * 2 + 1, psf_filter_sigma), 60)
		
#		psfgrad = (cv2.Sobel(psf, -1, 1,0,ksize=3)**2 + cv2.Sobel(psf, -1, 0,1,ksize=3)**2) ** 0.5
		#psfgrad2 = (cv2.Sobel(psfgrad, -1, 1,0,ksize=3)**2 + cv2.Sobel(psfgrad, -1, 0,1,ksize=3)**2) ** 0.5
		#cv2t.imshow("psfgrad2" , normalize(psfgrad2))
		
#		med_grad = float(np.mean(psfgrad))
		#print(med_grad)
		
		
		
#		psfgrad[np.where(psfgrad < med_grad * 3)] = 0
#		cv2t.imshow("psfgrad" , normalize(psfgrad))
		
#		floodmask = np.zeros((psf.shape[0] + 2, psf.shape[1] + 2), dtype=np.uint8)
#		cv2.floodFill(psfgrad, floodmask, (0,0), 255, loDiff = 0, upDiff = med_grad * 3, flags = 256 + 8 + cv2.FLOODFILL_MASK_ONLY)
#		mask = 1 - floodmask[1:-1, 1:-1]
#		cv2t.imshow("psfgradmask", mask * 255)



		


		
		floodmask = np.zeros((psf.shape[0] + 2, psf.shape[1] + 2), dtype=np.uint8)
		bl_psf = cv2.medianBlur(psf, 3)
		bl_psf = cv2.GaussianBlur(bl_psf, (3, 3), 0)
		
		y, x = np.ogrid[-args.diameter:args.diameter + 1, -args.diameter:args.diameter + 1]
		r = np.array((x * x + y * y) ** 0.5 / args.diameter, dtype = np.float32) ** 2 / 4
		
		minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(bl_psf)
		
		bl_psf_r = bl_psf + r * (maxVal - minVal)
		thr = (maxVal - minVal) / 2.0 + minVal
		
		while True:
			cv2.floodFill(bl_psf_r, floodmask, maxLoc, 2, loDiff = 1000.0, upDiff = 0, flags = 256 + 8 + cv2.FLOODFILL_MASK_ONLY)
			mask = floodmask[1:-1, 1:-1]
			bl_psf[np.where(mask > 0)] = 0
			minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(bl_psf)
			if maxVal < thr:
				break
		
		mask[0,:] = 0
		mask[-1,:] = 0
		mask[:, 0] = 0
		mask[:, -1] = 0
		
		cv2t.imshow("psf1" , normalize(psf ** 0.4))
		
#		mask = cv2.erode(mask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))
		#border = cv2.dilate(mask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3))) - mask


		low = np.median(psf)
		psf -= low
		psf[psf < 0] = 0
		psf[mask == 0] = 0

		thr = np.mean(psf[mask != 0]) / 5
		mask = (psf / (psf + thr)) ** 2
		mask = cv2.medianBlur(mask, 5)
		bask = cv2.GaussianBlur(mask, (3, 3), 0)
		cv2t.imshow("mask", mask)
		
		psf *= mask
		psf = lowblur(psf, thr, 3)
		
		if psf.sum() == 0:
			psf[args.diameter, args.diameter] += 1
		cv2t.imshow("psf2" , normalize(psf ** 0.4))
			
		self.target = target
		#target = cv2.erode(target, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))
		self.residual = cv2.GaussianBlur(psf, (args.diameter * 2 + 1, args.diameter * 2 + 1), args.psf_filter_sigma_res)
		#residual = cv2.dilate(residual, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))
		
		cv2t.imshow("srcfit" , normalize(psf ** 0.4))
		
		#self.target = np.average(psflist, axis = 0, weights = weights)
		
		self.set(psf)
		#cv2t.imshow("psf%d" % id(self) , normalize(self.psf ** 0.5))
	
	def gaussian(self, img, sigma):
		psf = np.zeros((args.diameter * 2 + 1, args.diameter * 2 + 1), dtype=np.float64)
		psf[args.diameter, args.diameter] = 1
		psf = cv2.GaussianBlur(psf, (args.diameter * 2 + 1, args.diameter * 2 + 1), sigma)
		self.set(psf)


	def update(self, cor_psf, dat, res, shape, moffat = None, scale_only = False):
		if not self.updated and not scale_only:
			self.updated = True
			
			hfr = self.hfr()
			lo_psf = cv2.erode(self.psf, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (int(hfr * 2), int(hfr * 2))))
			self.lo_thr = np.amax(lo_psf)
			
			#lo_psf = cv2.dilate(lo_psf, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (8, 8)))
			lo_psf = cv2.GaussianBlur(lo_psf, (args.diameter * 2 + 1, args.diameter * 2 + 1), args.diameter / 6)
			self.psf = cv2.max(self.psf, lo_psf)
			
			self.psf = cv2.GaussianBlur(self.psf, (args.diameter * 2 + 1, args.diameter * 2 + 1), 0.5)
			
			#self.psf += cv2.GaussianBlur(self.psf, (args.diameter * 2 + 1, args.diameter * 2 + 1), args.diameter / 3) * 1e-2
			
			if moffat:
				self.psf = moffat.apply(self.psf)
				
			self.accel_mask = np.zeros((args.diameter * 2 + 1, args.diameter * 2 + 1), dtype=np.float64)
			self.accel_mask[args.diameter, args.diameter] = 1
			self.accel_mask = cv2.GaussianBlur(self.accel_mask, (args.diameter * 2 + 1, args.diameter * 2 + 1), args.diameter / 2.0)
			self.accel_mask *= args.update_accel / np.amax(self.accel_mask)
			
			self.thr_mask = np.zeros((args.diameter * 2 + 1, args.diameter * 2 + 1), dtype=np.uint8)
			self.thr_mask[args.diameter, args.diameter] = 1
			self.thr_mask = cv2.dilate(self.thr_mask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (int(args.diameter * 1.6), int(args.diameter * 1.6))))
		
			
		#rs = cv2.sumElems(res)[0]
		#if rs == 0:
		#	return 1
			
		#res = cv2.divide(res, rs)
		
		
		dft_shape = (optimal_dft_size(shape[0] + self.shape[0]//2 - int(self.center[1]) + 1), 
		             optimal_dft_size(shape[1] + self.shape[1]//2 - int(self.center[0]) + 1))


		cordat = cv2.multiply(cor_psf, dat)

		res = cv2.flip(res, -1)
		eres = cv2.copyMakeBorder(res, top=0, bottom=dft_shape[0]-shape[0], left=0, right=dft_shape[1]-shape[1], borderType=cv2.BORDER_CONSTANT, value=0)
		RES = cv2.dft(eres, flags=cv2.DFT_COMPLEX_OUTPUT, nonzeroRows = shape[0])

		edat = cv2.copyMakeBorder(dat, top=0, bottom=dft_shape[0]-shape[0], left=0, right=dft_shape[1]-shape[1], borderType=cv2.BORDER_CONSTANT, value=0)
		DAT = cv2.dft(edat, flags=cv2.DFT_COMPLEX_OUTPUT, nonzeroRows = shape[0])
		cv2.mulSpectrums(DAT, RES, 0, DAT)
		cv2.idft(DAT, edat, flags=cv2.DFT_SCALE | cv2.DFT_REAL_OUTPUT )


		ecordat = cv2.copyMakeBorder(cordat, top=0, bottom=dft_shape[0]-shape[0], left=0, right=dft_shape[1]-shape[1], borderType=cv2.BORDER_CONSTANT, value=0)
		cv2.dft(ecordat, DAT, flags=cv2.DFT_COMPLEX_OUTPUT, nonzeroRows = shape[0])
		cv2.mulSpectrums(DAT, RES, 0, DAT)
		cv2.idft(DAT, ecordat, flags=cv2.DFT_SCALE | cv2.DFT_REAL_OUTPUT )
		
		
#res = cv2.UMat(res, (psf_shape[0] // 2, psf_shape[0] // 2 + img_shape[0]), (psf_shape[1] // 2, psf_shape[1] // 2 + img_shape[1]))

		#cor_psf = cv2.divide(cv2.filter2D(cv2.multiply(cor_psf, dat), -1, res, anchor = anchor), cv2.filter2D(dat, -1, res, anchor = anchor))
		#cor_psf = cv2.divide(cv2.filter2D(cv2.multiply(cor_psf, dat), -1, res), cv2.filter2D(dat, -1, res))
		#cor_psf = cv2.filter2D(cor_psf, -1, res, anchor = anchor)

		#cor_psf = cv2.divide(ecordat, edat)

		#cv2t.imshow("cor_psf1", normalize(cor_psf.get()))
		#cor_psf = np.array(cor_psf, dtype = np.float32)
		#cor_psf = cv2.getRectSubPix(cor_psf, (self.shape[1], self.shape[0]), (shape[1] - self.center[0], shape[0] - self.center[1]))
		
		cor_psf = cv2.divide(cv2.getRectSubPix(ecordat, (self.shape[1], self.shape[0]), (shape[1] - self.center[0], shape[0] - self.center[1])),
		                     cv2.getRectSubPix(edat   , (self.shape[1], self.shape[0]), (shape[1] - self.center[0], shape[0] - self.center[1])))
		
		cv2.patchNaNs(cor_psf, 1.0)
		cor_psf = cor_psf.get()
		
		#cor_psf /= np.mean(cor_psf)
		#print cor_psf
		#print "mean", np.mean(cor_psf)
		cor_psf = cv2.medianBlur(cor_psf, 5)
		cor_psf = cv2.GaussianBlur(cor_psf, (3,3), 0)
		
		m = np.sum(cor_psf * self.psf)
		
		if np.isnan(m):
			return 1.0
		
		if scale_only:
			return m
		
		if m > 0:
			cor_psf = np.array(cor_psf, dtype=np.float64) / m
		
		cor_psf = np.clip(cor_psf, 0.0001, 10000)
		
		sc = 0
#		while np.amin(cor_psf) < 1.0 - 0.6 / args.update_accel or np.amax(cor_psf) > 1.0 + 0.8 / args.update_accel:
		while np.amin(cor_psf + 0.6 / self.accel_mask) < 1.0 or np.amax(cor_psf - 0.8 / self.accel_mask) > 1.0:
			sc += 1
			cor_psf **= 0.8
		cv2t.imshow("cor_psf", normalize(cor_psf))
		print("cor_psf mean", m, "scale", sc)
		
		cor_psf = (cor_psf - 1) * self.accel_mask + 1
		
		if np.all(cor_psf >= 0) and not np.isnan(cor_psf).any():
			bl_dif = lowblur(self.psf, self.lo_thr, 1) - self.psf
			#bl_dif = cv2.GaussianBlur(self.psf, (3, 3), 0) - self.psf
		
			self.psf *= cor_psf
			
			update_thr = np.mean(self.psf[self.thr_mask == 0] ** 2)**0.5 * 2
			#print("update_thr", update_thr)
			if update_thr > 1e-8:
				self.psf *= ((self.psf + args.update_thr * update_thr) / (self.psf + update_thr)) #** 0.5

			self.psf += bl_dif

		cv2t.imshow("psf", normalize(self.psf ** 0.1))
		
		if self.psf.sum() == 0:
			self.psf[args.diameter, args.diameter] += 1

		self.psf /= self.psf.sum()
		self.psf_flip = self.psf[::-1, ::-1]
		return m

	def apply(self, img, flip = False, img_shape=None):
		if img_shape is None:
			img_shape = img.shape
			to_mat = True
		else:
			to_mat = False
		x = int(self.center[0])
		y = int(self.center[1])
		if not flip:
			x = -x
			y = -y
		
		anchor = (self.shape[0] // 2 + max(-y, 0), self.shape[1] // 2 + max(-x, 0))
		psf_shape = (self.shape[0] + max(y, 0), self.shape[1] + max(x, 0))
			
		dft_shape = (optimal_dft_size(img_shape[0] + psf_shape[0] + max(-y, 0)), optimal_dft_size(img_shape[1] + psf_shape[1] + max(-x, 0)))

		eimg = cv2.copyMakeBorder(cv2.UMat(img), top=0, bottom=dft_shape[0]-img_shape[0], left=0, right=dft_shape[1]-img_shape[1], borderType=cv2.BORDER_CONSTANT, value=0)
		epsf = cv2.UMat(dft_shape[0], dft_shape[1], cv2.CV_32F, 0.0)
		cv2.flip((self.psf_flip, self.psf)[flip], -1, cv2.UMat(epsf, (max(y, 0), psf_shape[0]), (max(x, 0), psf_shape[1])))

		IMG = cv2.dft(eimg, flags=cv2.DFT_COMPLEX_OUTPUT, nonzeroRows = img_shape[0])
		PSF = cv2.dft(epsf, flags=cv2.DFT_COMPLEX_OUTPUT, nonzeroRows = psf_shape[0])
		cv2.mulSpectrums(IMG, PSF, 0, IMG)

		cv2.idft(IMG, eimg, flags=cv2.DFT_SCALE | cv2.DFT_REAL_OUTPUT )
		#print("xy", x, y, "psf_shape", psf_shape, "anchor", anchor, "img_shape", img_shape, "dft_shape", dft_shape)
		res = cv2.UMat(eimg, (anchor[0], anchor[0] + img_shape[0]), (anchor[1], anchor[1] + img_shape[1]))
		if to_mat:
			return res.get()
		return res
		#return cv2.filter2D(img, -1, psf, anchor = anchor)
	
	def hfr(self):
		
		(x, y) = centroid(self.psf)
		x += args.diameter
		y += args.diameter
		
		print(self.psf.shape)
		psf = cv2.getRectSubPix(self.psf, (args.diameter * 2 + 1, args.diameter * 2 + 1), (x, y))
		return hfr(psf)
	
class VarPsf:
	def extract(self, img, ptlist, psf_filter_sigma):

		flatpsflist = []
		psflist = get_psf_list(img, ptlist)
		flatpsflist = np.array(psflist).reshape((len(psflist), -1))
	
		print("flatpsflist", flatpsflist.shape)

	
		cov = np.cov(flatpsflist)
		print("cov", cov)
		w, v = np.linalg.eig(cov)
		print("eig", w)
		print(v)
	
		print("eig shape:", v.shape)
		print("len psflist", len(psflist))
		print("cov shape", cov.shape)
	
		num_comp = min(20, w.shape[0])
	
		psfcomp = []
		self.psfcomplist = []
		for i, ev in enumerate(v):
			psf_s = np.zeros((args.diameter * 2 + 1, args.diameter * 2 + 1), np.float64)
			for j, c in enumerate(ev):
				psf_s += psflist[j] * float(c)
		
			psfcomp.append(psf_s.flatten())
			self.psfcomplist.append(psf_s)
		
			cv2.imwrite("psf%d.tif" % i, normalize(psf_s))
			if (i >= num_comp):
				break
	
		psfcomp = np.array(psfcomp).T
		
		psfcoef = []
		for psf in flatpsflist:
			psfcoef.append(np.linalg.lstsq(psfcomp, psf)[0])
		
		psfcoef = np.array(psfcoef).T
	
		Y = np.array([p[0] for p in ptlist])
		X = np.array([p[1] for p in ptlist])
	
		print("X", X)
		A = poly_array(X / 1000.0, Y / 1000.0, 3)
		print("coef fit")
	
		self.coef_fit = []
		for coefval in psfcoef:
			self.coef_fit.append(np.linalg.lstsq(A, coefval)[0])
	

		print("comp")
		self.comp_scale = []
		for c in self.coef_fit:
			v = poly_res(img.shape, c, 3)
			print(cv2.minMaxLoc(v))
			self.comp_scale.append(v)
		
		self.psfcomplist_flip = [psf[::-1, ::-1] for psf in self.psfcomplist]
		
		self.scale = None
		self.scale = self.apply(np.ones_like(img, dtype = np.float32))
		
		print("scale")
		print(cv2.minMaxLoc(self.scale))

	def apply(self, img, flip = False):
		res = np.zeros_like(img)
		for s, psf in zip(self.comp_scale, (self.psfcomplist_flip, self.psfcomplist)[flip]):
			res += cv2.filter2D(img * s, -1, psf)
		if self.scale is not None:
			res = res / self.scale
		return res


class TilePsf:
	def extract(self, img, ptlist, psf_filter_sigma, mask = None):
		self.shape = img.shape
		
		print("max img", np.amax(img))
		
		if mask is not None and np.array(mask).size <= 1:
			mask = None
		
		stars_tile = len(ptlist) / args.tiles ** 2 / 8
		print("stars per tile", stars_tile)
		stars_tile = max(stars_tile, args.stars_per_tile)
		h, w = self.shape

		self.tile_size_h = int((h + args.tiles - 1) / args.tiles)
		self.tile_size_w = int((w + args.tiles - 1) / args.tiles)
		tile_overlap = int(self.tile_size_h * args.tile_overlap)
		self.tile_overlap = tile_overlap
		
		self.tiles = []
		self.hfrlist = []
		self.sparse = False
		for y in range(0, h, self.tile_size_h):
			for x in range(0, w, self.tile_size_w):
				ul = (y, x)
				lr = (min(y + self.tile_size_h, h), min(x + self.tile_size_w, w))
				
				tile_overlap_e = tile_overlap
				eul_mo = (max(0, y - tile_overlap), max(0, x - tile_overlap))
				elr_mo = (min(y + self.tile_size_h + tile_overlap, h), min(x + self.tile_size_w + tile_overlap, w))
				if mask is not None and np.count_nonzero(mask[ul[0] : lr[0], ul[1]: lr[1]]) == 0:
					self.tiles.append((None, ul, lr, eul_mo, elr_mo))
					self.sparse = True
					continue
				while True:
					eul = (max(0, y - tile_overlap_e), max(0, x - tile_overlap_e))
					elr = (min(y + self.tile_size_h + tile_overlap_e, h), min(x + self.tile_size_w + tile_overlap_e, w))
				
					t_ptlist = []
					for  py, px in ptlist:
						if (py >= eul[0] + args.diameter and px >= eul[1] + args.diameter and py < elr[0] - args.diameter -1 and px < elr[1] - args.diameter -1
						        and  ( mask is None or
						          np.count_nonzero(mask[int(py) - args.diameter: int(py) + args.diameter + 1, int(px) - args.diameter: int(px) + args.diameter + 1]) == (args.diameter * 2 + 1)**2 )):
							t_ptlist.append((py - eul[0], px - eul[1]))
				
					psflist = get_psf_list(img[eul[0] : elr[0], eul[1]: elr[1]], t_ptlist, np.amax(img) * 0.95)
					if len(psflist) > stars_tile:
						break
					if eul == (0,0) and elr == (h,w):
						break
					tile_overlap_e *= 2
				
				
				if len(psflist) == 0:
					self.tiles.append((None, ul, lr, eul_mo, elr_mo))
					self.sparse = True
					self.hfrlist.append(0)
					continue
				
				
				if args.var_psf:
					psf = VarPsf()
					psf.extract(img[eul[0] : elr[0], eul[1]: elr[1]], t_ptlist, psf_filter_sigma)
					self.tiles.append((psf, ul, lr, eul, elr))
				else:
					psf = Psf()
					
					psf.from_psflist(psflist, psf_filter_sigma)
					self.tiles.append((psf, ul, lr, eul_mo, elr_mo))
					self.hfrlist.append(psf.hfr())

					

		print(self.tiles)
		
		self.update_hfr_field(np.clip(args.tiles, 1, 3))
		
		#target_psf = [t[0].target for t in self.tiles if t[0] is not None]
		#residual = [t[0].residual for t in self.tiles if t[0] is not None]

		#self.moffat_base = PsfMoffat(args.diameter * 2 + 1, psf = target_psf, residual = residual)
		#self.moffat = self.moffat_base
		self.moffat = None
		self.moffat_base = None
		

	def gaussian(self, img, sigma):
		
		self.shape = img.shape
		h, w = self.shape

		self.tile_size_h = int((h + args.tiles - 1) / args.tiles)
		self.tile_size_w = int((w + args.tiles - 1) / args.tiles)
		tile_overlap = int(self.tile_size_h * args.tile_overlap)
		self.tile_overlap = tile_overlap
		
		self.tiles = []
		self.hfrlist = []
		self.sparse = False
		for y in range(0, h, self.tile_size_h):
			for x in range(0, w, self.tile_size_w):
				ul = (y, x)
				lr = (min(y + self.tile_size_h, h), min(x + self.tile_size_w, w))
				
				tile_overlap_e = tile_overlap
				eul = (max(0, y - tile_overlap_e), max(0, x - tile_overlap_e))
				elr = (min(y + self.tile_size_h + tile_overlap_e, h), min(x + self.tile_size_w + tile_overlap_e, w))
				
				psf = Psf()
				psf.gaussian(img[eul[0] : elr[0], eul[1]: elr[1]], sigma)
				self.tiles.append((psf, ul, lr, eul, elr))
				self.hfrlist.append(psf.hfr())

		self.update_hfr_field(np.clip(args.tiles, 1, 3))
		self.moffat_base = None
		self.moffat = None
	
	def add_delta(self, delta):
		for psf, ul, lr, eul, elr in self.tiles:
			psf.center = delta


	_tile_cache = {}
	def _get_tile_mask(self, ul, lr, eul, elr):
		eh = elr[0] - eul[0]
		ew = elr[1] - eul[1]
		y1 = ul[0] - eul[0]
		y2 = lr[0] - eul[0]
		x1 = ul[1] - eul[1]
		x2 = lr[1] - eul[1]
		
		key = (eh, ew, y1, y2, x1, x2, self.tile_overlap)
		if key in self._tile_cache:
			return self._tile_cache[key]
		
		t_mask = np.zeros((eh, ew), dtype=np.float32)
		t_mask[y1:y2, x1:x2] = 1
		t_mask = cv2.blur(t_mask, (self.tile_overlap - 1, self.tile_overlap - 1))
		t_mask = cv2.UMat(t_mask)
		self._tile_cache[key] = t_mask
		return t_mask
		
		
	
	
	def apply(self, img, flip = False):
		# if self.shape != img.shape 
		if len(self.tiles) != args.tiles**2:
			raise ValueError("psf shape differs")
	
	
		#res = np.zeros_like(img)
		res = cv2.UMat(self.shape[0], self.shape[1], cv2.CV_32F, 0.0)
		if self.sparse:
			#mask = np.zeros_like(img)
			mask = cv2.UMat(self.shape[0], self.shape[1], cv2.CV_32F, 0.0)
		
		for psf, ul, lr, eul, elr in self.tiles:
			if psf is None:
				continue
			
			#tile = img[eul[0] : elr[0], eul[1]: elr[1]]
			tile = cv2.UMat(img, (eul[0], elr[0]), (eul[1], elr[1]))
			
			t_res = psf.apply(tile, flip, img_shape=(elr[0] - eul[0], elr[1] - eul[1]))

			t_mask = self._get_tile_mask(ul, lr, eul, elr)
			
			dest = cv2.UMat(res, (eul[0], elr[0]), (eul[1], elr[1]))
			cv2.add(dest, cv2.multiply(t_res, t_mask), dest)
			#res[eul[0] : elr[0], eul[1]: elr[1]] += t_res * t_mask
			if self.sparse:
				mdest = cv2.UMat(mask, (eul[0], elr[0]), (eul[1], elr[1]))
				cv2.add(mdest, t_mask, mdest)
				#mask[eul[0] : elr[0], eul[1]: elr[1]] += t_mask
			#cv2t.imshow("tile", mask)
			#cv2t.waitKey(0)

		if self.sparse:
			mask[mask == 0] = 100.0
			cv2.divide(res, mask, res)
			#res /= mask
		
		
		if self.moffat is not None:
			res = self.moffat.apply(res, img_shape=self.shape)
		return res
		
	def update(self, cor_psf, dat, res, scale_only = False):
		m = []
		for psf, ul, lr, eul, elr in self.tiles:
			#t_cor_psf = cor_psf[eul[0] : elr[0], eul[1]: elr[1]]
			t_cor_psf = cv2.UMat(cor_psf, (eul[0], elr[0]), (eul[1], elr[1]))
			#t_dat = dat[eul[0] : elr[0], eul[1]: elr[1]]
			t_dat = cv2.UMat(dat, (eul[0], elr[0]), (eul[1], elr[1]))
			#t_res = res[eul[0] : elr[0], eul[1]: elr[1]]
			t_res = cv2.UMat(res, (eul[0], elr[0]), (eul[1], elr[1]))
			t_mask = self._get_tile_mask(ul, lr, eul, elr)

			t_dat = cv2.multiply(t_dat, t_mask)
			m.append(psf.update(t_cor_psf, t_dat, t_res, (elr[0] - eul[0], elr[1] - eul[1]), self.moffat, scale_only))
		self.moffat_base = None
		self.moffat = None
		return np.mean(m)
	
	def interpolate(self, s):
		if self.moffat_base is not None:
			self.moffat = self.moffat_base.interpolate(s)

	def hfr(self):
		hfrlist = [h for h in self.hfrlist if h is not None]
		return np.median(hfrlist)
		
	def update_hfr_field(self, order):
		YX = np.indices((args.tiles, args.tiles), dtype = np.float64) / args.tiles
		X = YX[1]
		Y = YX[0]
	
		print("hfrlist", self.hfrlist)
		a = np.array(self.hfrlist)
	
		Xf = X.ravel()
		Yf = Y.ravel()
		
		af = a[a != 0]
		Xf = Xf[a != 0]
		Yf = Yf[a != 0]
	
		if len(af) == 0:
			self.hfr_field = np.ones((args.tiles, args.tiles), dtype=np.float32)
			return
	
		A = poly_array(Xf, Yf, order = order)
	
		coef = np.linalg.lstsq(A, af)[0]
	
		self.hfr_coef = coef
		
		A = poly_array(X, Y, order)
	
		self.hfr_field = np.dot(A, coef)
		
		self.hfr_field[self.hfr_field < 0.5] = 0.5
		self.hfr_field = np.array(self.hfr_field, dtype=np.float32)
		
		
		print(a.reshape((args.tiles, args.tiles)))
		print(self.hfr_field)

	def get_hfr_field(self, scale):
		hfr_field = self.hfr_field * scale
		print(hfr_field)
		f = cv2.UMat(hfr_field)
		return cv2.resize(f, (self.shape[1], self.shape[0]), interpolation=cv2.INTER_LINEAR)

	
	def save_debug(self, fn):
		vlist = []
		hlist = []
		for psf, ul, lr, eul, elr in self.tiles:
			if psf is None:
				continue
			p = np.zeros_like(psf.psf)
			p[args.diameter, args.diameter] = 1
			p = psf.apply(p)
			if self.moffat is not None:
				p = self.moffat.apply(p)

			hlist.append(p)
			if len(hlist) == args.tiles:
				vlist.append(np.concatenate(hlist, axis = 1))
				hlist = []
		res = np.concatenate(vlist, axis = 0)
		res = normalize(res ** 0.1)
		cv2.imwrite(fn, res)

def testPsf(shape, psf):
	test = np.zeros(shape, dtype = np.float64)
	test[:: args.diameter * 2, :: args.diameter * 2] = 1
	test = psf.apply(test)
	return test


def get_hfr(img, ptlist):
	hf = 0.0
	img = np.array(img, dtype = np.float32)
	for (y, x) in ptlist:
		psf = cv2.getRectSubPix(img, (args.diameter * 2 + 1, args.diameter * 2 + 1), (x, y))
		bg, w = psf_bg_w(psf)
		psf -= bg
		psf = np.abs(psf)
		
		hf += hfr(psf)
	return hf / len(ptlist)

def up_erode(img, r):
	size = (img.shape[1], img.shape[0])
	up = cv2.pyrUp(img)
	up = cv2.erode(up, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (r,r)))
	return cv2.pyrDown(up, dstsize = size)

def up_dilate(img, r):
	size = (img.shape[1], img.shape[0])
	up = cv2.pyrUp(img)
	up = cv2.dilate(up, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (r,r)))
	return cv2.pyrDown(up, dstsize = size)

def pyrBlur(img):
#	kernel = np.array([
#		1.0, 4.0, 6.0, 4.0, 1.0,
#		4.0,16.0,24.0,16.0, 4.0,
#		6.0,24.0,36.0,24.0, 6.0,
#		4.0,16.0,24.0,16.0, 4.0,
#		1.0, 4.0, 6.0, 4.0, 1.0]).reshape((5,5)) / 256.0
#	return cv2.filter2D(img, -1, kernel)
#1	return cv2.GaussianBlur(img, (5,5), 0.9) * 0.9
#	return cv2.GaussianBlur(img, (5,5), 0.6)
	size = (img.shape[1], img.shape[0])
	img = cv2.pyrDown(img)
	return cv2.pyrUp(img, dstsize = size)


def mean_clip2(src2):
	u_src2 = cv2.UMat(src2)
	avg = cv2.mean(u_src2)[0]
	for i in range(0, 8):
		print("avg", avg)
		mask = cv2.inRange(u_src2, 1e-10, float(avg) * 4)
		avg = cv2.mean(u_src2, mask=mask)[0]
	
	avg = float(avg)
	print("mean_{}: {}".format(src2.shape, avg ** 0.5))
	return avg


def mean_clip(src):

	avg2 = mean_clip2(src ** 2)
	return avg2 ** 0.5

class reg_pm_2:
	def step1(self, res_a, lm, m = None):
		g = []
		
		res = np.mean(res_a, axis = 0)
		p_res = np.pad(res, ((1, 1), (1, 1)), 'edge')
		
		#p_res -= cv2.GaussianBlur(p_res, (11,11), 0)
		
		grad_x = np.diff(p_res, axis = 1)
		grad_xm = np.mean([grad_x[0:-1, :], grad_x[1:, :]], axis = 0)
		grad_y = np.diff(p_res, axis = 0)
		grad_ym = np.mean([grad_y[:, 0:-1], grad_y[:, 1:]], axis = 0)
	
		lap = cv2.GaussianBlur(res, (7,7), 0) - res
		lap = np.pad(lap, ((1, 1), (1, 1)), 'edge')
		lap = np.mean([lap[:-1,:-1], lap[:-1,1:], lap[1:,:-1], lap[1:,1:]], axis = 0)
		lap[np.where(lap > 0)] = 0
			
			
		g = grad_xm ** 2 + grad_ym ** 2 + lap ** 2 * 2

		self.m = m
		if self.m is None:
			self.m = mean_clip2(g)

		l = lm ** 2 * self.m
		g = 1. / (1. + (g / l)) ** 2
		g = cv2.GaussianBlur(g, (3,3), 0)
		
		cv2t.imshow("g2_%d" % g.shape[1], normalize(g))
		self.gx = np.mean([g[1:, 1:-1], g[:-1, 1:-1]], axis = 0)
		self.gy = np.mean([g[1:-1, 1:], g[1:-1, :-1]], axis = 0)

	def step2(self, c):
		grad_x = np.diff(c, axis = 1)
		grad_y = np.diff(c, axis = 0)
		grad_xg = grad_x * self.gx
		grad_yg = grad_y * self.gy


		del grad_x
		del grad_y
	
		grad_xg2 = np.diff(grad_xg, axis = 1)
	
		grad_yg2 = np.diff(grad_yg, axis = 0)
		del grad_xg
		del grad_yg
	
		reg = np.zeros_like(c)
		reg[:, 1:-1] +=  grad_xg2
		reg[1:-1, :] +=  grad_yg2

		del grad_xg2
		del grad_yg2
		return reg


class Pyramid(object):
	def __init__(self, depth = 3, reg = None, l = None, mode = 'normal'):
		if reg is not None:
			depth = len(reg)
			self.l_mul = l
			self.reg = reg
		else:
			self.l_mul = [1.0] * depth
			self.reg = [1.0] * depth
		
		self.depth = depth
		self.mode = mode


	def gaussian_pyr(self, img, mode = 'normal', weight = None, wp = None):
		d = args.reg_morph_radius * 2 + 1
		#kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (d, d))
		if weight is not None:
			wp = self.gaussian_pyr(weight)
		img = np.array(img)
		gp = [img]
		for i in range(self.depth):
			try:
				if mode == 'openclose':
					op = up_erode(img, d)
					op = up_dilate(op, d)
					cl = up_dilate(img, d)
					cl = up_erode(cl, d)
					img = (op + cl) / 2.0
				elif mode == 'openclose2':
					i1 = up_erode(img, d)
					i1 = up_dilate(i1, d)
					i1 = up_dilate(i1, d)
					i1 = up_erode(i1, d)
					i2 = up_dilate(img, d)
					i2 = up_erode(i2, d)
					i2 = up_erode(i2, d)
					i2 = up_dilate(i2, d)
					img = (i1 + i2) / 2.0
				elif mode == 'open':
					img = up_erode(img, d)
					img = up_dilate(img, d)
				elif mode == 'close':
					img = up_dilate(img, d)
					img = up_erode(img, d)
				elif mode == 'median':
					img = cv2.medianBlur(img, d)
				elif mode == 'erode':
					img = up_erode(img, d)
				elif mode == 'dilate':
					img = up_dilate(img, d)
				elif mode != 'normal':
					raise Exception('bad mode', mode)
				
				if wp is not None:
					img = img * wp[i]
				img = cv2.pyrDown(img)
				if wp is not None:
					img = cv2.divide(img, wp[i + 1])
					img[wp[i + 1] == 0.0] = 0
			except:
				pass
			gp.append(img)
		return gp

	def laplacian_pyr(self, img, weight = None):
		wp = None
		if weight is not None:
			wp = self.gaussian_pyr(weight)

		gp = self.gaussian_pyr(img, self.mode, wp = wp)
		
		lp = [np.array(gp[self.depth])]
		for i in range(self.depth,0,-1):
			size = (gp[i - 1].shape[1], gp[i - 1].shape[0])
			up = cv2.pyrUp(gp[i], dstsize = size)
			lap = cv2.subtract(gp[i-1], up)
			if wp is not None:
				lap[wp[i-1] == 0] = 0
			lp.append(lap)
		return lp[::-1]

	def collapse(self, pyr):
		depth = len(pyr)
		up = pyr[depth - 1]
		for i in range(depth - 2, -1, -1):
			size = (pyr[i].shape[1], pyr[i].shape[0])
			up = cv2.pyrUp(up, dstsize = size)
			up = up + pyr[i]
		return up
	
	def upscale_to(self, src, shape):
		shapes = []
		while shape[0] - src.shape[0] > 1 and shape[1] - src.shape[1] > 1:
			shapes.append(shape)
			shape = ((shape[0] + 1) // 2, (shape[1] + 1) // 2)
		for shape in shapes[::-1]:
			size = (shape[1], shape[0])
			src = cv2.pyrUp(src, dstsize = size)
		return src

	def downscale(self, src, n):
		for i in range(n):
			src = cv2.pyrDown(src)
		return src

class reg_pm_2_pyr(Pyramid):
	def __init__(self, depth = 3, reg = None, l = None):
		super(reg_pm_2_pyr, self).__init__(depth, reg, l)
		self.m = [1] * self.depth
	
	def step1(self, res_a, eval_noise):
		
		gp = []
		for res in res_a:
			gp.append(self.gaussian_pyr(res))
		
		self.levels = []
		for i, (img, lm) in enumerate(zip(zip(*gp), self.l_mul)):
			pm = reg_pm_2()
			
			if eval_noise:
				pm.step1(img, lm)
				self.m[i] = pm.m
			else:
				pm.step1(img, lm, self.m[i])
			self.levels.append(pm)
	
	def step2(self, c):
		gp = self.gaussian_pyr(c)
		resp = []
		for pm, img, r in zip(self.levels, gp, self.reg):
			print(pm.gx.shape, img.shape)
			resp.append(pm.step2(img) * r) 
		
		
		return self.collapse(resp)
		


class reg_pyr(Pyramid):
	def __init__(self, depth = 3, reg = None, l = None):
		super(reg_pyr, self).__init__(depth, reg, l)
		self.m = [1] * self.depth
	
	def step1(self, res_a, eval_noise):
		
		lps = []
		res = np.mean(res_a, axis = 0)
		lps.append(self.laplacian_pyr(res))
			
		self.g_lp = []
		for i in range(self.depth):
			g = []
			for j in range(len(lps)):
				lap = lps[j][i]
				
				g.append(lap * lap)
				
			
			g = np.mean(g, axis = 0)
			if eval_noise:
				self.m[i] = mean_clip2(g)
			lm = self.l_mul[i] ** 2 * self.m[i]

			g = 1. / (1. + (g / lm) ** 1)
			#if i == self.depth - 1:
			#	g = cv2.medianBlur(g, 3)
			#g = cv2.erode(g, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))
			#g = cv2.GaussianBlur(g, (3,3), 0)

			#g = cv2.dilate(g, cv2.getStructuringElement(cv2.MORPH_CROSS, (3,3)))
		
			cv2t.imshow("g_pyr_%d" % len(self.g_lp), g)
			self.g_lp.append(g)
		

	def step2(self, c):
		lp = self.laplacian_pyr(c)
		reg_lp = []
		for i in range(self.depth):
			lp[i] *= self.g_lp[i]
			reg_lp.append(-(lp[i] - cv2.GaussianBlur(lp[i], (3,3), 0))  * self.reg[i])
			#-cv2.Laplacian(lp[i], -1, ksize = 3))
		
		
		return self.collapse(reg_lp)



class reg_pyr_flatten(Pyramid):
	def __init__(self, depth = 3, reg = None, l = None):
		super(reg_pyr_flatten, self).__init__(depth, reg, l, mode = args.reg_mode)
		self.m = [1] * self.depth
	
	def step1(self, res_a, eval_noise, channel_weights = None, weights_a = None):
		
		if channel_weights is None:
			channel_weights = [1.0] * len(res_a)
		lps = []
		#res = np.mean(res_a, axis = 0)
		for i, res_i in enumerate(res_a):
			if weights_a is not None:
				lps.append(self.laplacian_pyr(res_i, weight = weights_a[i]))
			else:
				lps.append(self.laplacian_pyr(res_i))
			
		self.g_lp = []
		for i in range(self.depth):
			g = []
			g2 = []
			for j in range(len(lps)):
				lap = lps[j][i]
				lap = lap - pyrBlur(lap)
				
				g2.append(lap * lap * channel_weights[j]**2)
				g.append(lap * channel_weights[j])
				
			
			g2 = np.mean(g2, axis = 0)
			g = np.mean(g, axis = 0)
			
			if eval_noise:
				print(g2.shape)
				self.m[i] = mean_clip2(g2)
			
			#g = g2
			g = g * g
			
			lm = self.l_mul[i] ** 2 * self.m[i]

			g = (g / lm)# ** 0.9 # ** 0.5
			
			g = cv2.GaussianBlur(g, (5,5), 0.7)# 0.6
			
			g = g / (g + 1)
			#g = cv2.GaussianBlur(g, (5,5), 0.7)
			
			
			#g = cv2.medianBlur(g, 3)

			#g = up_erode(g, 2)
			##g = up_dilate(g, 3)
			##g = cv2.GaussianBlur(g, (5,5), 2)# 0.6
			
			g = 1.0 - g ** self.reg[i]
			
			cv2t.imshow("g_flat_%d" % len(self.g_lp), g)
			self.g_lp.append(g)
		

	def step2(self, c, weight = None):
		lp = self.laplacian_pyr(c, weight = weight)
		reg_lp = []
		for i in range(self.depth):
			lap = lp[i]
			lap = lap * self.g_lp[i]
			lap = lap - pyrBlur(lap)
			
			reg_lp.append(-lap)
		
		
		return self.collapse(reg_lp)


def pyr_noise_level(img):
	pyr = Pyramid(3)
	img = np.array(img, dtype=np.float)
	lp = pyr.laplacian_pyr(img)
	
	sigma2 = []
	for l in lp[:-1]:
		sigma2.append(mean_clip2(l * l))
	return np.sum(sigma2) ** 0.5


def signal_level(img):
	img = np.array(img, dtype=np.float)
	avg, sigma = cv2.meanStdDev(cv2.blur(img, (5, 5)) - cv2.blur(img, (100, 100)))
	sigma = float(sigma)
	print("signal_level", sigma)
	return sigma

def bg_level(src):
	print("1")
	res = cv2.UMat(src)
	res = cv2.GaussianBlur(res, (5, 5), 0)
	res = cv2.erode(res, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5, 5)), iterations = args.diameter)
	
	print("1")
	#res = cv2.blur(res, (3, 3))
	#print("2")
	#res = cv2.erode(res, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (args.diameter * 4 + 1,args.diameter * 4 + 1)))
	#print("3")
	res = cv2.blur(res, (args.diameter * 4 + 1,args.diameter * 4 + 1))
	print("4")
	res = cv2.blur(res, (args.diameter * 4 + 1,args.diameter * 4 + 1))
	print("5")
	res = cv2.blur(res, (args.diameter * 4 + 1,args.diameter * 4 + 1))
	print("6")
	res = cv2.blur(res, (args.diameter * 4 + 1,args.diameter * 4 + 1))
	print("7")
	return res.get()

def average_src(src, white, weights = None):
	sum_src = np.zeros_like(src[0], dtype = np.float32)
	sum_weights = np.zeros_like(src[0], dtype = np.float32)
	
	if weights is None:
		weights = [1.0] * len(src)
	
	
	def average1(i, lock):
		
		s = np.array(src[i], dtype = np.float32)
		w = np.array(weights[i], dtype = np.float32)
		gr = int(np.count_nonzero(w) ** 0.5 / 2) + 1
		w2 = cv2.blur(w, (gr, gr))
		w2 = cv2.blur(w2, (gr, gr))
		w2 = cv2.blur(w2, (gr, gr))
		w2 **= 6
		w *= w2
		del w2
		cv2t.imshow("w", normalize(w))
		
		s = s / white[i] * w
		
		with lock:
			sum_src[:] += s
			sum_weights[:] += w

	pfor.pfor(average1, list(range(0, len(src)))) 

	res = cv2.divide(sum_src, sum_weights)
	res[np.where(sum_weights <= 0.001)] = args.zero / white[0]
	return res, sum_weights, args.zero / white[0]

def cluster_exps(exps):
	minv = np.amin(exps)
	maxv = np.amax(exps)
	hist, bins = np.histogram(exps, bins = np.arange(minv - 1.0, maxv + 2.0, 1.0), density=False)
	
	print(hist, bins)
	
	res = []
	while True:
		max_i = np.argmax(hist)
		if hist[max_i] == 0:
			break
		
		low = max(max_i - 1, 0)
		high = min(max_i + 1, len(hist) - 1)
		
		w = np.where((exps >= bins[low]) & (exps < bins[high + 1]))[0]
		res.append(w)
		
		hist[low] = 0
		hist[high] = 0
		hist[max_i] = 0
		print(hist)
	print(res)
	return res
		
def mask_overlaps(m1, m2):
	nz1 = np.count_nonzero(m1)
	nz2 = np.count_nonzero(m2)
	cross = np.count_nonzero(m1 * m2)
	ret = cross > min(nz1, nz2) * 0.8
	print('mask_overlap', nz1, nz2, cross, ret)
	return ret

def cluster_overlap(exp_clusters, masks):
	res = []
	for cl in exp_clusters:
		res_clusters = []
		res_masks = []
		for c in cl:
			added = False
			for i, m in enumerate(res_masks):
				if mask_overlaps(m, masks[c]):
					res_clusters[i].append(c)
					added = True
					break
			if not added:
				res_clusters.append([c])
				res_masks.append(masks[c])
		res += res_clusters
	return res


def save_output(it = None):
	res_ch = []
	for outc in range(0, outchannels):
		if background is not None:
			res_wb = res[outc] - background[outc] + args.zero / 65535.0
		else:
			res_wb = res[outc]
		res_wb = res_wb * white_balance[outc] + args.zero / 65535.0 * (1.0 - white_balance[outc])
		#res_ch.append(cv2.multiply(res_wb, 65535.0, dtype = cv2.CV_16UC1))
		res_ch.append(res_wb)

	outidx = 0
	for i, fn in enumerate(outfiles):
		res_i = []
		for j in range(0, outfiles_c[i]):
			res_i.append(res_ch[outidx])
			outidx += 1
		kwargs = {}
                
		if np.asarray(transp[i]).shape == res[0].shape:
			if len(res_i) == 1:
				kwargs = {'planarconfig': 'contig', 'photometric' : 'minisblack'}
			#res_i.append(transp[i])
		
		res_i = cv2.merge(res_i)
		if it is not None:
			fn = "i%04d_%s" % (it, fn)
		tifffile.imsave(fn, res_i, **kwargs)
		del res_i
	del res_ch

def weighted_lin_fit(a, b, w):

	ma = np.average(a, weights = w)
	mb = np.average(b, weights = w)
	print("weighted_lin_fit average ", ma, mb)
	
	sa = np.average((a - ma)**2, weights = w) ** 0.5
	sb = np.average((b - mb)**2, weights = w) ** 0.5
	cor = np.average((a - ma) * (b - mb), weights = w) / (sa * sb)
	
	m = cor * sb / sa
	c = mb - m * ma
	return m, c

def histogram_weight(src, n = 50):
	isrc = np.array(src * n, dtype=np.int)
	isrc[isrc < 0] = 0
	isrc[isrc >= n] = n - 1
	hist, edges = np.histogram(isrc, bins=np.arange(n+1))
	print("hist", hist)
	m = max(np.median(hist), 1)
	print(1.0 / (np.array(hist, dtype=np.float) + m))
	histval = np.take(hist, isrc)
	ret = 1.0 / (np.array(histval, dtype=np.float) + m)
	return ret


def AvgDev(img):
	return np.mean(np.abs(img - np.median(img)))

def MAD(img):
	return np.median(np.abs(img - np.median(img)))

def BVMV(img):
	m = np.median(img)
	mdif = np.abs(img - m)
	Y = mdif / ( 9 * np.median(mdif))
	mask = (Y < 1.0)
	
	Ym2 = Y[mask]**2
	
	ret = img.size * np.sum(mdif[mask]**2 * (1.0 - Ym2)**4) / np.abs(np.sum((1.0 - Ym2)*(1.0 - 5 * Ym2)))**2
	return ret


def gradient_dif_filter(dif, w, r):
	difw = dif * w
	r1 = r // 4
	difw = cv2.blur(difw, (r1, r1))
	mean_w = cv2.blur(w, (r1, r1))
	difw = cv2.blur(difw, (r1, r1))
	mean_w = cv2.blur(mean_w, (r1, r1))
	mean_w = cv2.max(mean_w, 1e-6)

	if np.amin(dif.shape) > r * 5:
		dif = cv2.divide(difw, mean_w)
		cv2.patchNaNs(dif, 0)
		dif = cv2.min(cv2.max(dif, -1), 1)
		print("check difi1", np.isnan(dif).any(), cv2.minMaxLoc(dif))

		shape = dif.shape
		resize_w = int((shape[1] + r - 1) / r)
		resize_h = int((shape[0] + r - 1) / r)
		dif = cv2.resize(dif, (resize_w, resize_h), interpolation=cv2.INTER_AREA)
		dif = cv2.medianBlur(dif, 5)
		print("check difi2", np.isnan(dif).any(), cv2.minMaxLoc(dif))
		dif = cv2.resize(dif, (shape[1], shape[0]), interpolation=cv2.INTER_LINEAR)
		print("check difi3", np.isnan(dif).any(), cv2.minMaxLoc(dif))
		dif = cv2.blur(dif, (r, r))
		dif = cv2.blur(dif, (r, r))
		dif = cv2.blur(dif, (r, r))
		dif = cv2.blur(dif, (r, r))
		dif = cv2.blur(dif, (r, r))
	else:
		difw = cv2.blur(difw, (r, r))
		mean_w = cv2.blur(mean_w, (r, r))
		difw = cv2.blur(difw, (r, r))
		mean_w = cv2.blur(mean_w, (r, r))
		difw = cv2.blur(difw, (r, r))
		mean_w = cv2.blur(mean_w, (r, r))
		dif = cv2.divide(difw, mean_w)
		cv2.patchNaNs(dif, 0)

	dif = np.array(dif, dtype=np.float32)

	print("check difi4", np.isnan(dif).any(), cv2.minMaxLoc(dif))

	return dif


def gradient_dif_filter2(dif, weight, r):
	h, w = dif.shape
	resize_w = int((w + r - 1) / r)
	resize_h = int((w + r - 1) / r)
	res = np.empty((resize_h, resize_w), dtype = dif.dtype)
	for i in range(0, resize_h):
		for j in range(0, resize_w):
			for s in range(1, w):
				y0 = np.clip(int((i - s) * h / resize_h), 0, h - 1)
				y1 = np.clip(int((i + 1 + s) * h / resize_h), 0, h - 1)
				x0 = np.clip(int((j - s) * w / resize_w), 0, w - 1)
				x1 = np.clip(int((j + 1 + s) * w / resize_w), 0, w - 1)
				tile = dif[y0:y1, x0:x1]
				tile_w = weight[y0:y1, x0:x1]
				tile_m = tile[tile_w > 0]
				if len(tile_m) > 0:
					break
			res[i, j] = np.median(tile_m)
	res = cv2.resize(res, (w, h), interpolation=cv2.INTER_LINEAR)
	res = cv2.blur(res, (r, r))
	res = cv2.blur(res, (r, r))
	res = cv2.blur(res, (r, r))


	return res

def mean_sigma_clip(a0, kappa=2, n=10, weights = None):
	a0 = np.array(a0)
	a = a0
	weights0 = weights
	for i in range(0, n):
		m = np.average(a, weights=weights)
		d2 = (a - m)**2
		s2 = np.average(d2, weights=weights)
		d02 = (a0 - m)**2
		where = d02 <= s2 * kappa**2
		a = a0[where]
		if weights0 is not None:
			weights = np.array(weights0)[where]
	return m



#########################################################################################################


np.set_printoptions(precision=6, linewidth=160)
set_htr_size(args.diameter)

if args.test_ptlist:
	order = args.psf_filter_order
	f = open(args.infile[0], 'rb')
	args, img, ptlist = pickle.load(f)
	f.close()
	args.psf_filter_order = order
	set_htr_size(args.diameter)
	f_ptlist = filter_ptlist(img, ptlist, debug=True)
	sys.exit(0)

from cv2t import cv2t


src = []
src_white = []
src_sigma = []
src_median = []
src_signal_level = []

psfsrc = []
psfsrc_white = []
imgidx = []
transp = []
weights = []
maxval = []

exps = []
gradient_weight = []
time_mov = []

outmap = []

outfiles = []
outfiles_c = []
outchannels = 0


nc = 1
for f in args.outfile.split(','):
	try:
		nc, fn = f.split(':')
		nc = int(nc)
	except:
		fn = f
	outfiles.append(fn)
	outfiles_c.append(nc)
	outchannels += nc

outc_num = [0] * outchannels


outc = 0
for f in args.infile:
	fn = f
	exp = 1.0
	grad_weight = 1.0
	t_mov = 0.0
	if not os.path.isfile(fn):
		try:
			outc, fn = fn.split(':')
			outc = int(outc)
		except:
			pass

	for fn in fn.split(','):
		if not os.path.isfile(fn):
			try:
				fn, t_mov = fn.split('#')
				t_mov = float(t_mov)
			except:
				pass
			try:
				fn, exp_s, grad_weight_s = fn.split('@')
				exp = float(exp_s)
				grad_weight = float(grad_weight_s)
			except: 
				try:
					fn, exp = fn.split('@')
					exp = float(exp)
				except:
					pass
		if not os.path.isfile(fn):
			print(fn, "not found")
			sys.exit(1)
			


		print("Read ", fn)
		img = tifffile.TiffFile(fn).asarray(memmap=True)
		if img is None:
			print('Failed to load fn:', fn)
			sys.exit(1)

		white = np.iinfo(img.dtype).max
		img = np.atleast_3d(img)
		if args.crop is not None:
			img = img[args.crop[0], args.crop[1]]
		h, w, channels = img.shape
		col = [1,1,1,3,3][channels]
	
		transp_c = white
		if channels > col:
			transp_c = img[:,:,col]
		
		#mask = np.amin(img, axis = 2)
		#img = extrapolate_transp(img[:,:,0:col], mask, add = True)
	
		outi = outc
		cidx = []
		for c in range(0, col):
			cidx.append(len(src))
			src.append(img[:,:, c])
			src_white.append(white)
			weights.append(transp_c)
			src_sigma.append(pyr_noise_level(img[:,:, c]))
			if channels > col:
				src_median.append(np.median(img[:,:, c][transp_c > 0]))
			else:
				src_median.append(np.median(img[:,:, c]))
			src_signal_level.append(signal_level(img[:,:, c]))
			
			mv = np.amax(cv2.medianBlur(img[:,:, c], 3))
			maxval.append(mv)
		
			outmap.append(outi)
			outc_num[outi] += 1
			outi += 1
		
			exps.append(exp)
			gradient_weight.append(grad_weight)
			time_mov.append(t_mov)
		imgidx.append(cidx)
		transp.append(transp_c)

outmap = np.array(outmap)

srcmap = [[] for outc in range(outchannels)]
for c in range(0, len(src)):
	outc = outmap[c]
	srcmap[outc].append(c)

exps = np.array(exps, dtype = np.float64)

scales = 2.0 ** exps

scales /= np.amax(scales)
scales *= args.scale

white_balance = [1.0] * outchannels
for i, wb in enumerate(args.wb):
	white_balance[i] = wb
print("white balance", white_balance)

#channel_scales = [1.0] * outchannels
#for c in range(0, len(scales)):
#	outc = outmap[c]
#	channel_scales[outc] *= scales[c]
#	print("scale", c, outc, scales[c], channel_scales[outc])
#for outc in range(0, outchannels):
#	channel_scales[outc] **= 1.0 / outc_num[outc]

#print("channel_scales", channel_scales)
#print("outc_num", outc_num)

channel_weights = [0.0] * outchannels


for c in range(0, len(src)):
	print("scale sigma", scales[c], src_sigma[c], scales[c] * src_sigma[c], src_signal_level[c])
	outc = outmap[c]
	channel_weights[outc] += (src_signal_level[c] / (scales[c] * src_sigma[c])) ** 2 / outc_num[outc]

#channel_weights = np.array([1.0, 0.0, 0.0, 0.0])

if args.channel_weights:
	channel_weights = [0.0] * outchannels
	for i, cw in enumerate(args.channel_weights):
		channel_weights[i] = cw

channel_weights = np.array(channel_weights)#** 0.5
print("channel_weights", channel_weights)

exp_clusters = cluster_exps(exps)
exp_clusters = cluster_overlap(exp_clusters, weights)
print('exp_clusters', exp_clusters)

sum_src, sumweights, sum_zero = average_src(src, src_white, weights)

print("sumweights", sumweights)
print("sum_src", sum_src)

if args.psffile:
	for f in [args.psffile]:
		img = tifffile.TiffFile(f).asarray(memmap=True)
		if img is None:
			print('Failed to load fn1:', f)
			sys.exit(1)

		white = np.iinfo(img.dtype).max
		img = np.atleast_3d(img)
		h, w, channels = img.shape
		col = [1,1,1,3,3][channels]

		for c in range(0, col):
			psfsrc.append(img[:,:, c])
			psfsrc_white.append(white)

	sum_psfsrc, sum_psfweights, sum_zero = average_src(psfsrc, psfsrc_white)
else:
	psfsrc = src
	psfsrc_white = src_white
	sum_psfsrc = sum_src



col = len(src)



#s_res = cv2.GaussianBlur(sum_src,(5,5),0)
#s_res = cv2.erode(s_res, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (7,7)))
s_res = sum_src

if args.gaussian == 0.0:
	ptlist = find_ptlist(sum_psfsrc, zero = sum_zero)

	show = normalize(sum_psfsrc ** 0.4)
	for (y, x) in ptlist:
		cv2.circle(show, (int(x),int(y)), 10, (255), 1)
	cv2.imwrite("pts_all.tif", show)

	del sum_psfsrc

	cluster_idx = [None] * col
	cluster_ptlists = []
	for cl_i, cl in enumerate(exp_clusters):
		cl_psfsrc = []
		cl_psfsrc_white = []
		cl_weights = []
		for i in cl:
			cl_psfsrc.append(psfsrc[i])
			cl_psfsrc_white.append(psfsrc_white[i])
			cl_weights.append(weights[i])
			cluster_idx[i] = cl_i

		sum_psfsrc, sum_psfweights, sum_zero = average_src(cl_psfsrc, cl_psfsrc_white, cl_weights)
		f = open("pts%d.data" % cl_i, 'wb')
		pickle.dump((args, sum_psfsrc, ptlist), f)
		f.close()
		cl_ptlist = filter_ptlist(sum_psfsrc, ptlist)

		show = normalize(sum_psfsrc ** 0.4)
		for (y, x) in cl_ptlist:
			cv2.circle(show, (int(x),int(y)), 10, (255), 1)
		cv2.imwrite("pts%d.tif" % cl_i, show)

		del sum_psfsrc
		cluster_ptlists.append(cl_ptlist)


psf = [None] * col
c_ptlists = [None] * col

skipped = [False] * col

res = [None] * outchannels
var = [1.0] * outchannels
reg_prev = [None] * outchannels
#hist_weight = [None] * outchannels
kappa = 10.0

glob_var = [1.0] * outchannels
glob_var_sumw = [0.0] * outchannels
glob_var_w = [0.0] * outchannels


cor_scales = np.array([ 1.0 ] * col)

try:
	if args.load_psf:
		f = open(args.load_psf, 'rb')
		psf = pickle.load(f)
		f.close()
		for c in range(0, col):
			psf[c].save_debug("testpsf%d.tif" % (c))

	else:
		raise IOError
except:
	psf = [None] * col

	def step1(c, lock):
		print("start0  %d" % c)
		psf[c] = TilePsf()
		psfs = np.array(psfsrc[c], dtype = np.float32) ** args.gamma / psfsrc_white[c]
		if args.gaussian == 0.0:
			deltax = args.delta_x * (time_mov[c] - 0.5)
			deltay = args.delta_y * (time_mov[c] - 0.5)
			print("delta %f %f" % (deltax, deltay))
			c_ptlists[c] = [(y, x) for (y, x) in cluster_ptlists[cluster_idx[i]]]
			psf[c].extract(psfs, c_ptlists[c], args.psf_filter_sigma, mask = weights[c])
			psf[c].add_delta((deltax, deltay))
			psf[c].save_debug("testpsf%d.tif" % (c))

		else:
			psf[c].gaussian(psfs, args.gaussian)
		print("end0  %d" % c)


	pfor.pfor(step1, list(range(0, col)))

	if args.save_psf:
		f = open(args.save_psf, 'wb')
		pickle.dump(psf, f)
		f.close()

hfr_list = []
for c in range(0, col):
	hfr_list.append(psf[c].hfr())

def get_scaled_src_w(c, eval_noise = 1.0, med_bg = True):
	scaled_src = (np.array(src[c], dtype=np.float32) ** args.gamma * scales[c] + args.zero * (1.0 - scales[c])) / src_white[c]

	weight = np.empty_like(src[c], dtype=np.float32)
	weight[:, :] = weights[c]

	if med_bg:
		scaled_src[weight == 0.0] = (src_median[c] * scales[c] + args.zero * (1.0 - scales[c])) / src_white[c]
	else:
		scaled_src[weight == 0.0] = 0.0

	gr = args.diameter
	w2 = cv2.blur(weight, (gr,gr)) / src_white[c]
	w2 = cv2.blur(w2, (gr,gr))
	w2 *= w2
	weight *= w2

	print("noise_level", c, scales[c] * src_sigma[c])
	weight /= src_white[c] * (scales[c] * src_sigma[c])**(2 * eval_noise)

	return scaled_src, weight

def lin_fit(X, Y, weight, size = 500):
	wm = cv2.blur(weight, (size, size))
	Xm = cv2.divide(cv2.blur(cv2.multiply(X, weight), (size, size)), wm)
	cv2.patchNaNs(Xm, 0)
	Ym = cv2.divide(cv2.blur(cv2.multiply(Y, weight), (size, size)), wm)
	cv2.patchNaNs(Ym, 0)

	Xc = cv2.subtract(X, Xm)
	Yc = cv2.subtract(Y, Ym)
	
	Xcw = cv2.multiply(Xc, weight)
	Ycw = cv2.multiply(Yc, weight)
	
	cov = cv2.mean(cv2.multiply(Xcw, Ycw))[0]
	var = cv2.mean(cv2.pow(Xcw, 2))[0]
	b = cov / var
	dif = cv2.subtract(Y, cv2.multiply(X, b))
	dif = cv2.divide(cv2.blur(cv2.multiply(dif, weight), (size, size)), wm)
	return dif, b
	
def blur_dif(X, Y, weight, size = 3000):
	wm = cv2.blur(weight, (size, size))
	wm = cv2.blur(wm, (size, size))
	wm = cv2.blur(wm, (size, size))
	wm = cv2.blur(wm, (size, size))

	dif = cv2.subtract(Y, X)
	cv2.multiply(dif, weight, dif)
	dif = cv2.blur(dif, (size, size))
	dif = cv2.blur(dif, (size, size))
	dif = cv2.blur(dif, (size, size))
	dif = cv2.blur(dif, (size, size))
	cv2.divide(dif, wm, dif)
	cv2.patchNaNs(dif, 0)
	return dif


src_lp_cor = [None] * col
partial = [True] * outchannels

ref0 = [None] * outchannels

lp_depth_min = args.lp_depth_min
lp_depth_max = int(np.log2(np.amin(src[0].shape)) - 2)
print("depth max", lp_depth_max)
pyr = Pyramid(depth = lp_depth_max - 1, mode = 'median')
for outc in range(0, outchannels):
	srcpyr = [[None] * outc_num[outc] for i in range(lp_depth_min, lp_depth_max)]
#	hp_pyr = [None] * lp_depth_min
#	hp_pyr_stack = [[] for i in range(lp_depth_min)]
	weightpyr = [[None] * outc_num[outc] for i in range(lp_depth_min, lp_depth_max)]
	for si in range([1, 3][args.update_scales]):
		for ic in range(outc_num[outc]):
			c = srcmap[outc][ic]
			scaled_src, weight = get_scaled_src_w(c, med_bg = False)
			

			#mask = np.ones_like(weight)
			#mask[weight == 0.0] = 0.0
		
			slp = pyr.laplacian_pyr(scaled_src, weight = weight)
		
			weight *= gradient_weight[c]
			
			if np.count_nonzero(weight == 0) == 0:
				print("partial False", outc)
				partial[outc] = False
			else:
				if not partial[outc]:
					weight *= 0
			
			wgp = pyr.gaussian_pyr(weight)
			
			for i in range(lp_depth_min, lp_depth_max):
				if ic == 0:
					cv2.imwrite("p%d.tif" % (i), normalize(slp[i]))
					cv2.imwrite("w%d.tif" % (i), normalize(wgp[i]))

		
				srcpyr[i - lp_depth_min][ic] = slp[i]
			
				weight = cv2.erode(wgp[i], cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (11,11)))
				weight = cv2.GaussianBlur(weight, (7,7), 0)
				weight[weight == 0.0] = 1e-12
				weightpyr[i - lp_depth_min][ic] = weight
		
#			for i in range(lp_depth_min):
#				hp_pyr_stack[i].append(slp[i])
#				if len(hp_pyr_stack[i]) >= 5 or ic == outc_num[outc] - 1:
#					if hp_pyr[i] is None:
#						hp_pyr[i] = np.median(hp_pyr_stack[i], axis = 0) / ((outc_num[outc] + 4) // 5)
#					else:
#						hp_pyr[i] += np.median(hp_pyr_stack[i], axis = 0) / ((outc_num[outc] + 4) // 5)
#					hp_pyr_stack[i] = []

		
	
		res_pyr = []
		for i in range(lp_depth_max - lp_depth_min):
			weightpyr[i] = np.array(weightpyr[i])
			for it in range(0, 5):
				avg = np.average(srcpyr[i], axis = 0, weights = weightpyr[i])
				d2 = (np.array(srcpyr[i]) - avg[np.newaxis, :, :]) ** 2
				s2 = np.average(d2, axis = 0, weights = weightpyr[i])
				weightpyr[i][d2 > s2 * 4] = 0
			
			res_pyr.append(avg)

		up_scales = []
		res_hp = pyr.collapse(res_pyr[0: (lp_depth_max - lp_depth_min) - 2])
		for ic in range(outc_num[outc]):
#			ratios = []
#			c = srcmap[outc][ic]
#			for i in range((lp_depth_max - lp_depth_min) // 2 + 1):
#				r = 1.0
#				for j in range(0, 10):
#					dif2 = (srcpyr[i][ic] * r - res_pyr[i]) ** 2
#					dvar = np.mean(dif2)
#					where = dif2 < dvar * 4
#					sl = srcpyr[i][ic][where]
#					rl = res_pyr[i][where]
#					w = rl / (np.abs(rl) + np.mean(rl*rl)**0.5 * 10.0)
#					s_m = np.mean(sl * w)
#					r_m = np.mean(rl * w)
#					r = r_m / s_m
#					print(j, r)
#				ratios.append(r)
#			r = mean_sigma_clip(ratios)
#			print('ratios', r, ratios)
#			up_scales.append(r)
			
			#c = srcmap[outc][ic]
#			src_hp = []
#			for i in range((lp_depth_max - lp_depth_min) - 2):
#				src_hp.append(srcpyr[i][ic])
#			src_hp = pyr.collapse(src_hp)
#			r = 1.0
#			for j in range(0, 6):
#				dif2 = (src_hp * r - res_hp) ** 2
#				dvar = np.mean(dif2)
#				where = dif2 < dvar * 4
#				sl = src_hp[where]
#				rl = res_hp[where]
#				rl2 = rl * rl
#				w = rl / (rl2 + np.mean(rl2) * 10.0)
#				s_m = np.mean(sl * w)
#				r_m = np.mean(rl * w)
#				r = r_m / s_m
#				print("ratios", j, r)
#			up_scales.append(r)
			

			ratios = []
			c = srcmap[outc][ic]
			for i in range(1, lp_depth_max - lp_depth_min - 1):
				r = 1.0
				
				src2 = srcpyr[i][ic] ** 2
				res2 = res_pyr[i] ** 2

				src2 = cv2.erode(src2, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))
				res2 = cv2.erode(res2, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))

				var_src = np.mean(src2)  
				var_res = np.mean(res2) 
				where = srcpyr[i][ic] * res_pyr[i] > 0 
				where = np.logical_and(where, src2 > var_src / 2)
				where = np.logical_and(where, res2 > var_res / 2)

				ld = np.log(res2[where]) - np.log(src2[where])
				r = np.exp(mean_sigma_clip(ld) / 2)
				ratios.append(r)
			r = mean_sigma_clip(ratios)
			print('ratios', r, ratios)
			up_scales.append(r)

			
			
		up_scales = np.array(up_scales) / np.exp(np.mean(np.log(up_scales)))
		if args.update_scales:
			for ic in range(outc_num[outc]):
				c = srcmap[outc][ic]
				scales[c] *= up_scales[ic]
		print("up_scales", up_scales)
			
	
	del srcpyr
	del weightpyr
	
	mean = pyr.collapse(res_pyr)
	cv2.imwrite("pyr%d.tif" % (outc), normalize(mean**0.1))
	mean = pyr.upscale_to(mean, src[0].shape)
	
	ref0[outc] = np.array(mean, copy=True)
	
	for ic in range(outc_num[outc]):
		c = srcmap[outc][ic]
		scaled_src, weight = get_scaled_src_w(c)
		#scaled_src = pyr.downscale(scaled_src, lp_depth_min)
		#weight = pyr.downscale(weight, lp_depth_min)
		#dif = scaled_src - mean
		#bl = int(args.gradient_radius / 2**lp_depth_min) + 1
		#dif = gradient_dif_filter(dif, weight, bl)


		dif = scaled_src - mean
		print("check dif", np.isnan(dif).any(), cv2.minMaxLoc(dif))
		dif = gradient_dif_filter(dif, weight, args.gradient_radius)
		print("check dif2", np.isnan(dif).any(), cv2.minMaxLoc(dif))
		dif = pyr.downscale(dif, lp_depth_min)

			
		src_lp_cor[c] = dif
		cv2.imwrite("dif%d.tif" % (c), normalize(dif))

	
	#hp_pyr.append(mean)
	#mean = pyr.upscale_to(mean, src[0].shape)
	res[outc] = mean
	ref0[outc] = np.array(mean, copy=True)

print("scales",scales)
res = np.array(res)
offset = np.empty_like(res)

noise = noise_level(sum_src) + 0.00000001
print("noise:", noise * 65535.0)

#reg1 = reg_pm_2_pyr(reg = args.reg2, l = args.reg_lambda2)
#reg2 = reg_pyr(reg = args.reg4, l = args.reg_lambda4)
reg3 = reg_pyr_flatten(reg = args.regf, l = args.reg_lambdaf)

#for c in range(0, col):
#	psf[c].interpolate(args.psf_filter_moffat_scale)

#for c in range(0,col):
#	h = psf[c].hfr() ** 2
#	src_sigma[c] *= h
#	print("hfr cor", c, h, src_sigma[c])

background = None

bright_w = [0] * outchannels
ref_bg = [None] * outchannels

lines = [None] * col

it = args.iter

for i in range(0, it):

	print(i)
	if i <= it / 2:
		for c in range(0, col):
			psf[c].interpolate(args.psf_filter_moffat_scale * np.clip((i - 5) / (min(it - 3, args.update_iter) - 5), 0, 1))
			#psf[c].interpolate(args.psf_filter_moffat_scale)


	#s_res = np.mean(res, axis = 0)
	#

	#if i < 10:
	#noise = noise_level(s_res) + 0.00000001
	#print("noise:", noise * 65535.0)
	

	setup = (i <= 3)

	if args.poly_bg > 0 and (i == 6 or i == 15):
		background, bg_stddev = poly_bg(res, order = args.poly_bg, scale = 2, erode = 1, kappa = args.poly_bg_kappa, it = args.poly_bg_iter, save_mask = "mask", get_stddev = True)
		background = np.array(background, dtype=np.float32)


	sum_w = np.zeros_like(res)
	sum_corw = np.zeros_like(res)
	var_next = np.zeros_like(res)

	#reg_weight = 1 + (float(i) / it * args.robustness ** 0.5)**2
	#dat_eps = 1.0 / args.robustness**4
	med_hfr = np.median(hfr_list)
	
	
	u_res = [None] * outchannels
	for outc in range(0, outchannels):
		if background is not None:
			res[outc] = np.maximum(res[outc], background[outc] - bg_stddev[outc])
			
			bright_w[outc] = np.maximum(res[outc] - background[outc], 0)
			bright_w[outc] = 4 * bright_w[outc] / (bright_w[outc] + 3*glob_var[outc]**0.5)
			bright_w[outc] = cv2.UMat(bright_w[outc])
			bright_w[outc] = cv2.dilate(bright_w[outc], cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (7,7)))
			bright_w[outc] = cv2.GaussianBlur(bright_w[outc], (11, 11), 0)
	
		#hist_weight[outc] = histogram_weight(res[outc])
	
		#offset[outc] = np.amin(res[outc]) * args.dering
		if setup:
			offset[outc] = 0.0
		else:
			offset[outc] = bg_level(res[outc]) * args.dering
		
		ref_bg[outc] = ref0[outc] - offset[outc]
		
		res[outc] -= offset[outc]
		res[outc][res[outc] <=0] = 0.0000000001
		
		
		u_res[outc] = cv2.UMat(res[outc])
	
	def step2(c, lock):
		if skipped[c]:
			return
		print("start2 %d %d" %(i,c))
		kappa = 3 #np.clip(3 - i, 0, 5) + 3
		outc = outmap[c]
		scaled_src, weight = get_scaled_src_w(c, eval_noise = np.clip(i / 5.0, 0.0, 1))

		if setup:
			b_res = u_res[outc]
		else:
			b_res = psf[c].apply(u_res[outc])

		scaled_src -= offset[outc]
		
		
		if src_lp_cor[c] is not None:
			scaled_src -= pyr.upscale_to(src_lp_cor[c], scaled_src.shape)
		else:
			dif = scaled_src - b_res.get()
			dif = gradient_dif_filter(dif, weight, args.gradient_radius)
			
			scaled_src -= dif
			src_lp_cor[c] = pyr.downscale(dif, lp_depth_min)

		scaled_src = cv2.UMat(scaled_src)

		
		dif = cv2.subtract(b_res, scaled_src)

		if i == 3:
			lines[c], mask = detect_lines(dif)
			draw_lines(mask, lines[c], 127)
			print("lines%d_%d.tif" % (i, c))
			cv2.imwrite("lines%d_%d.tif" % (i, c), mask)
			
		
		#if i < 8:
		#	kwargs = {'planarconfig': 'contig', 'photometric' : 'minisblack'}
		#	fn = "dif%02d_%04d.tif" % (i, c)
		#	tifffile.imsave(fn, 0.5 - dif.get(), **kwargs)


		print("maxval", c, maxval[c])
		#where_over = np.where((src[c] > maxval[c] * args.overexp) & (scaled_src < b_res.get()))
		where_over = np.where(src[c] > maxval[c] * args.overexp)
		#where_over = np.where((src[c] > maxval[c] * args.overexp) & (scaled_src < res[outc]))
		
		w_over = np.clip((src[c][where_over] - maxval[c] * args.overexp) / src_white[c] + cv2.blur(dif, (3, 3)).get()[where_over] / scales[c], 0, 1)
		print("w_over", np.amax(w_over), np.amax(cv2.blur(dif, (3, 3)).get()[where_over] / scales[c]))
		
		weight[where_over] /= 1.0 + (w_over * 100) ** 2
		
		if lines[c]:
			draw_lines(weight, lines[c], 0)

		#d = blur_dif(scaled_src, ref_bg[outc], cv2.divide(weight, ref_bg[outc]))
		#cv2.add(scaled_src, d, scaled_src)
		#cv2t.imshow("lin_fit", normalize(d.get()))
		

		cor = cv2.max(cv2.divide(scaled_src, b_res), 0.00001)
		
		#del b_res
		
		#cv2.patchNaNs(cor, 0.00001)
		
		#dif_mean = np.average(dif, weights = weight)
		##dif_mean = cv2.mean(dif)[0]
		#dif_mean = np.median(dif)
		##cv2.subtract(dif, dif_mean, dif)
		#dif -= dif_mean
		##print(dif_mean)
		##print("dif_mean %d %d %f" %(i,c,dif_mean))
		
		dat = cv2.pow(dif, 2)
		
		dat_m = dat.get()
		with lock:
			var_next[outc] += dat_m
		
		if not setup:
			#clip = np.where(np.logical_and(dat_m > var[outc] * kappa * kappa, dif.get() < 0))
			clip = np.where(dat_m > var[outc] * kappa * kappa)
			weight[clip] = 0
		#del dif
		
		

		var_sumw = np.sum(dat_m * weight)
		var_w = np.sum(weight)
		
		with lock:
			glob_var_sumw[outc] += var_sumw
			glob_var_w[outc] += var_w

		mean_dat = var_sumw / var_w
		print("dat_mean %d %d %f glob %f" %(i,c,mean_dat**0.5, glob_var[outc]**0.5))
		
		if not setup:
			mean_dat = glob_var[outc]

		if args.update_iter > 0 and i > args.update_iter:
			dat_psf = cv2.pow(cv2.add(dat, mean_dat / args.update_psf_robustness**2), -0.25)
			cv2.multiply(dat_psf, weight, dat_psf)
			#dat_psf *= weight

		cv2.add(dat, mean_dat / args.robustness**2, dat)
		cv2.pow(dat, -0.25 * args.cut_out, dat) #-0.3)
		#dat = (dat + mean_dat / args.robustness**2) ** -0.25
		cv2.multiply(dat, weight, dat)
		#del weight
		
		hfr_f = psf[c].get_hfr_field(1.0 / med_hfr)
		
		#print("check hfr_f1", cv2.minMaxLoc(hfr_f.get()))
		
		cv2.log(hfr_f, hfr_f)
		#print("check hfr_f2", cv2.minMaxLoc(hfr_f.get()))
		cv2.multiply(hfr_f, bright_w[outc], hfr_f)
		#print("check hfr_f3", cv2.minMaxLoc(hfr_f.get()))
		cv2.exp(hfr_f, hfr_f)
		#print("check hfr_f4", cv2.minMaxLoc(hfr_f.get()))

		
		cv2.divide(dat, hfr_f, dat)
		
		#cv2t.imshow("hfr_f", normalize(hfr_f.get()))
		#del hfr_f

		dat = cv2.erode(dat, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))

			

		#dat[clip] = 0
		
		
		cor_psf = cv2.pow(cor, 2)
		#cor_psf = cv2.medianBlur(cor, 5)

		
		if setup:
			corw = cv2.multiply(dat, cor)
			w = dat
		else:
			corw = psf[c].apply(cv2.multiply(dat, cor), flip = True)
			w =  psf[c].apply(dat, flip = True)
		#cor = psf.apply(cor, flip = True)
		
		cv2.max(corw, 0, corw)
		cv2.max(w, 0, w)
		
		corw = corw.get()
		w = w.get()
		#corw[np.where(corw < 0.0)] = 0.0
		#w[np.where(w < 0.0)] = 0.0



#		if i >= 4 and i <= 8:
#			dif = gradient_dif_filter(dif, dat, args.gradient_radius)
#			dif = pyr.downscale(dif, lp_depth_min)

#			
#			src_lp_cor[c] -= dif
#			
#			with lock:
#				cv2.imwrite("dif_up%d.tif" % (c), normalize(dif))
#				cv2.imwrite("dif%d.tif" % (c), normalize(src_lp_cor[c]))

#				print("check dif", np.isnan(dif).any(), cv2.minMaxLoc(dif))



		
		if c % (col // 5) == 0:
			cv2t.imshow("dat%d" % c, normalize(-dat.get()**0.5))
		
				
		if args.update_iter > 0 and i > args.update_iter:
			m = psf[c].update(cor_psf, dat_psf, u_res[outc])
			
			
			if i == args.update_iter + 1:
				for k in range(3):
					b_res = psf[c].apply(u_res[outc])
					cor_psf = cv2.max(cv2.divide(scaled_src, b_res), 0.00001)
					cor_psf = cv2.pow(cor_psf, 2)
#					cor_psf = cv2.medianBlur(cor_psf, 5)

		
					m = psf[c].update(cor_psf, dat_psf, u_res[outc])
			
			cor_scales[c] = m
			#if args.update_scales:
			#	scales[c] /= m**0.2
#		else:
#			cor_scales[c] = psf[c].update(cor, dat, u_res[outc], scale_only = True)

		with lock:
			#cv2.imwrite("reg_p%d_%d.tif" % (c, i), normalize(reg_p))
#			print("check res", np.isnan(res[c]).any(), cv2.minMaxLoc(res[c]))
		
#			print("check reg_p", np.isnan(reg_p).any(), cv2.minMaxLoc(reg_p))
#			print("check reg_m", np.isnan(reg_m).any(), cv2.minMaxLoc(reg_m))
			#print("check dat", np.isnan(dat.get()).any(), cv2.minMaxLoc(dat))
			#print("check w", np.isnan(w).any(), cv2.minMaxLoc(w))
			#print("check corw", np.isnan(corw).any(), cv2.minMaxLoc(corw))
			#print("scales_c ", c, scales[c])

			sum_w[outc] += w
			sum_corw[outc] += corw

		print("end2 %d %d" %(i,c))
	pfor.pfor(step2, list(range(0, col)) if i % 2 == 0 else reversed(list(range(0, col))))
	
	del u_res

	var = var_next


	norm_w = []
	for sw in sum_w:
		norm_w.append(sw / np.median(sw))
	norm_w = 1.0 / np.mean(norm_w, axis = 0)
	
	norm_w = cv2.dilate(norm_w, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (args.weight_dilate * 2 + 1, args.weight_dilate * 2 + 1)))
	norm_w = np.clip(norm_w, 1.0, 10.0)
	norm_w = cv2.GaussianBlur(norm_w, (args.weight_dilate * 2 + 1, args.weight_dilate * 2 + 1), 0)
		
	cv2t.imshow("norm_w", norm_w / 10.0)

	#reg1.step1(res, i < args.lambda_eval_iter)
	#reg2.step1(res, i < args.lambda_eval_iter)
	reg3.step1(res + offset, i < args.lambda_eval_iter, channel_weights = channel_weights, weights_a = sum_w)

	print("cor_scales", cor_scales)
	def step3(outc, lock):
		sc_var = np.mean(np.log(cor_scales[outmap == outc]) ** 2)
		kappa = 2
		while True:
			sc_where = (np.log(cor_scales) ** 2 > sc_var * kappa) & (outmap == outc)
			if np.count_nonzero(sc_where) <= outc_num[outc] / 4 + 1:
				break
			kappa *= 1.1
		print("cor_scales", outc, sc_var, kappa)
		print(sc_where)
		print(scales)
		for c in range(col):
			if sc_where[c]:
				mul = 1.0 / np.clip(cor_scales[c] ** 0.3, 0.8, 1.3)
				#scales[c] *= mul
				#src_lp_cor[c] *= mul
				src_lp_cor[c] = None
		print(scales)
	
	
	
	
	
		var[outc] /= outc_num[outc]
		
		if i < 5 or glob_var[outc] > glob_var_sumw[outc] / glob_var_w[outc]:
			glob_var[outc] = glob_var_sumw[outc] / glob_var_w[outc]
		glob_var_sumw[outc] = 0.0
		glob_var_w[outc] = 0.0
		

		
		#print("check res_pre", np.isnan(res[outc]).any(), cv2.minMaxLoc(res[outc]))

		#cor = cv2.divide(sum_corw[outc], sum_w[outc])

		#res[outc] *= cor
	
	
		#reg_weight = args.reg * min(float(i) / (it * 0.3), 1)
		reg = reg3.step2(res[outc], weight = sum_w[outc])
		#print("check reg1", np.isnan(reg).any(), cv2.minMaxLoc(reg))
		#print("check reg2", np.isnan(reg).any(), cv2.minMaxLoc(reg))

		missing = np.where(sum_w[outc] < 1e-6)
		sum_w[outc][missing] = 1e-6

		#reg = cv2.divide(reg, sum_w[outc])
		#reg *= np.mean(sum_w[outc])
		#print("check reg3", np.isnan(reg).any(), cv2.minMaxLoc(reg))

		#mean, stddev = cv2.meanStdDev(reg)
		#stddev = float(stddev)
		#reg = np.clip(reg, -stddev * 3, stddev * 3)

		
		#hp = cv2.Laplacian(res[outc], -1, ksize = 7)
		#coefs = np.linalg.lstsq(np.array([reg.ravel()]).T, hp.ravel().T)[0]
		#print("reg fit ", coefs)
		
		
		#reg = cv2.divide(reg, np.clip(res[outc], 1024.0/65535.0, 1))

		#print("check reg4", np.isnan(reg).any(), cv2.minMaxLoc(reg))

		reg[missing] = 0
		#print("check reg5", np.isnan(reg).any(), cv2.minMaxLoc(reg))


		#reg[np.where(src[c] > 0.95)] /= 1000
	
		#reg = cv2.GaussianBlur(reg, (3,3), 0.6)
		#if i > 0:
		#	reg = reg * 0.1 + reg_prev[outc] * 0.9
		#	reg_prev[outc] = reg
		#else:
		#	reg_prev[outc] = (reg * 0.1)
		
		#print reg
		#reg *= (1.0 + args.reg_plus * (np.clip(norm_w, 1 , 3) - 1.0) / 2.0)
		reg_p = cv2.max(reg, 0)
		reg_m = cv2.min(reg, 0)
	
		#res[outc] *= cv2.divide(1 + reg_p, 1 - reg_m)
		if True: #not setup:
			ures = cv2.UMat(res[outc])
			rings = cv2.dilate(ures, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (9, 9)))
			rings = cv2.erode(rings, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (9, 9)))
			cv2.subtract(rings, res[outc], rings)
			rings = cv2.max(rings, 0)
		
		
			bg = cv2.erode(ures, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (7,7)), 4)
			bg = cv2.GaussianBlur(bg, (args.diameter * 2 + 1, args.diameter * 2 + 1), 0)
			
			ringmask = cv2.max(cv2.subtract(ures, bg), 0)
			
			rm_2 = cv2.pow(ringmask, 2)
			rm_var = cv2.GaussianBlur(rm_2, (args.diameter * 2 + 1, args.diameter * 2 + 1), 0)
			rm_var_glob = float(cv2.mean(rm_2)[0])

			cv2.divide(ringmask, cv2.add(cv2.add(ringmask, cv2.multiply(cv2.pow(rm_var, 0.5), args.dering2_thr)), rm_var_glob**0.5), ringmask)

			print("check ringmask", cv2.minMaxLoc(ringmask))
			#cv2.threshold(ringmask, 0.1, 1, cv2.THRESH_TOZERO, ringmask)
			
			ringmask = cv2.dilate(ringmask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (args.dering2_radius * 2 + 1,args.dering2_radius * 2 + 1)))
			
			ringmask_e = cv2.erode(ringmask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (args.dering2_radius * 2 + 1, args.dering2_radius * 2 + 1)))
			thr = float(cv2.mean(cv2.pow(ringmask_e, 2))[0]) ** 0.5
			cv2.threshold(ringmask_e, thr, 1, cv2.THRESH_BINARY, ringmask_e)
			
			cv2.subtract(ringmask, ringmask_e, ringmask)
			cv2.max(ringmask, 0, ringmask)

			print("check ringmask2", cv2.minMaxLoc(ringmask), thr)
			ringmask = cv2.GaussianBlur(ringmask, (5, 5), 0)

			cv2.threshold(ringmask, thr * 3, 1, cv2.THRESH_TOZERO, ringmask)

			#ringmask = cv2.pow(ringmask, 2)
		
			cv2.multiply(rings, ringmask, rings)
			rings = cv2.GaussianBlur(rings, (3, 3), 0)
			print("check rings", cv2.minMaxLoc(rings))
			
			rings = rings.get()

			reg += rings * args.dering2
		
			rings = np.clip(rings, 0, np.mean(rings) * 20)
			cv2t.imshow("rings%d" % outc, normalize(rings))
			cv2t.imshow("ringmask%d" % outc, normalize(ringmask.get()))


		#cor = cv2.divide(sum_corw[outc] + reg_p, sum_w[outc] - reg_m)
		cor = cv2.divide(sum_corw[outc], sum_w[outc])
		cor[missing] = 1.0
		cor[np.where(cor < 0.01)] = 0.01
		cor[np.where(cor > 100.)] = 100.
		
		res[outc] *= cor
		if setup:
			res[outc] = cv2.GaussianBlur(res[outc], (5, 5), 0)
		else:
			res[outc] += reg
		
		minVal, maxVal, minLoc, maxLoc = cv2.minMaxLoc(res[outc])
		maxVal = min(maxVal, 1.0)
		avg, stddev = cv2.meanStdDev(res[outc])
		
		floodmask = np.zeros((res[outc].shape[0] + 2, res[outc].shape[1] + 2), dtype=np.uint8)
		print(minVal, maxVal, minLoc, maxLoc)
		cv2.floodFill(res[outc], floodmask, minLoc, 2, loDiff = 0, upDiff = (maxVal - minVal) * args.top_thres, flags = 256 + 8 + cv2.FLOODFILL_MASK_ONLY + cv2.FLOODFILL_FIXED_RANGE)
		floodmask = floodmask[1:-1, 1:-1]
		#cv2t.imshow("fm%d" % outc, normalize(floodmask))
		top = cv2.dilate(res[outc], cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (21,21)))
		top = cv2.erode(top, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (21,21)))
		top[floodmask > 0] = 0
		cv2t.imshow("top%d" % outc, normalize(top))
		res[outc] = np.maximum(res[outc], top)

		if i < args.blur_iter:
			res[outc] = cv2.GaussianBlur(res[outc], (17, 17), args.blur_sigma * (1.0 - i / args.blur_iter))

		#print("check res", np.isnan(res[outc]).any(), cv2.minMaxLoc(res[outc]))

		#print("check reg_p", np.isnan(reg_p).any(), cv2.minMaxLoc(reg_p))
		#print("check reg_m", np.isnan(reg_m).any(), cv2.minMaxLoc(reg_m))
		

		reg_p = np.clip(reg_p, 0, np.mean(reg_p) * 20)
		reg_m = np.clip(reg_m, np.mean(reg_m) * 20, 0)
		
		cv2t.imshow("reg_p%d" % outc, normalize(reg_p))
		cv2t.imshow("reg_m%d" % outc, normalize(-reg_m))
		
		res[outc] += offset[outc]
		cv2t.imshow("res%d" % outc, normalize(np.clip(res[outc] - 1010/65535.0,10/65535.0,1)** args.show_gamma))
		cv2t.imshow("off%d" % outc, normalize(offset[outc]))

#		if args.gaussian == 0.0:
#			curhfr = get_hfr(res[outc], ptlist)
#			print("hfr", curhfr)
			#if i < 3 and curhfr > 2:
			#	res = cv2.erode(res, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3)))
	pfor.pfor(step3, list(range(0, outchannels)))


	if args.update_iter > 0 and i > args.update_iter:
#		hfr_list = []
#		for c in range(0, col):
#			hfr_list.append(psf[c].hfr())

#		print("res_hfr", i, hfr_list)

		for c in range(0, col):
			psf[c].save_debug("testpsf%d_%d.tif" % (i, c))
		
#		if i == 150:
#			s_hfrlist = sorted(hfr_list)
#			thr = s_hfrlist[col / 8]
#			for c in range(0, col):
#				skipped[c] = hfr_list[c] > thr
#				cv2.imwrite("testpsf%d.tif" % (c), normalize(testPsf(psfsrc[c].shape, psf[c])))
#			print "hfr thr", thr
#			print skipped
			
	save_output(i)
	

	cv2t.waitKey(1)
	gc.collect()
	
#for outc in range(0, outchannels):
#	if args.gaussian == 0.0:
#		psf = TilePsf()
#		pts = set_psf_level(res[outc], ptlist)
#		psf.extract(res[outc], pts, 0.01)
#		cv2.imwrite("testpsf%d_res.tif" % (outc), normalize(testPsf(src[outc].shape, psf)))
	

save_output()

cv2t.waitKey(1000)