/*	parsek.cpp

	Probabilistic Alignment Raw Sticher Experiment from Kentucky

	A parsec is a unit of distance corresponding to a "parallax of one
	second" -- i.e., the distance at which 1 AU subtends an angle of
	one arcsecond, or approximately 3.26 light years. That's big. And
	that's what Parsek is supposed to do: make big images by looking
	at small angle offsets.

	See Usage line for options.

	2023 by H. Dietz
*/

#include <ctype.h>
#include <errno.h>
#include <fcntl.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <setjmp.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <time.h>
#include <sys/types.h>
#include <stdint.h>
#include <sys/mman.h>

// C++ and OpenCV stuff
#include <opencv2/opencv.hpp>
#include <iostream>
#include <cstdio>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
 
using namespace cv;
using namespace std;

// OpenMP
#include <omp.h>

// use 16 bit instad of 8 bit for all processing...
// except alignment, which is done on 8 bit grayscale
#define	CV_TYP	CV_16UC3
#define	VecTYP	Vec3w
#define	pixel_t	ushort

#define	NO_CFA	"unknown"		// CFA pattern is unknown

char	*prog;				// program name for error messages
int	xin, yin;			// x,y input dimensions
int	bpp = 0;			// set to number of initial bits per pixel
int	bppl, bppr;			// left and right shifts to get to 16 bits
int	aligniter = 50;			// alignment iterations, generally 10-50
const	char *cfa;			// Color Filter Array pattern as 2x2, e.g., "RGGB"
int	dcraw = 1;			// use dcraw, else use LibRaw
double	epsilon = 1e-12;		// for TermCriteria, was 1e-10
char	*emfile = 0;			// error model file name
float	first = 0;			// filter distance from first with this weighting, 0 disables
float	ingamma = 2.0;			// default gamma for linearizing 8-bit input pixels
ushort	gtable[256];			// gamma table for 8-bit stuff
float	histexp = 1;			// histogram noise model exponent
int	interp = 1;			// interpolate between pixels in inputs
double	maxcon = 4.0;			// maximum confidence from 1 to 1000000
double	invmaxcon = 1.0/4.0;		// 1.0 / maxcon
int	neigh = 0;			// do neighborhood filtering?
const	char *outname = "parsek.png";	// default output file name
float	outliers = 100.0;		// filter outliers with this weighting, 0 disables
double	prune = 0.01;			// pruning cutoff on alignment matrix
int	rawin = 0;			// input files in a raw format?
int	smooth = 0;			// smooth 24bpp to raw conversion
int	verbose = 0;			// how verbose messages?
int	xout = 0;			// pixel columns in output file
int	yout = 0;			// pixel rows in output file
float	distlim;			// compute as 1.0 / scale factor

pixel_t	*org;			// original image

// color offsets from CFA
// note that input PGM file is in CFA order, NOT OpenCV order
#define	CFAX(POS)	((POS)&1)
#define	CFAY(POS)	(((POS)>>1)&1)
int	cfar, cfag, cfaG, cfab;	// encoded offsets
int	cof[4];			// color of, translates local offset to BGR 0,1,2

// basic contribution from one image to one pixel color channel
#define	bval	val[0]
#define	gval	val[1]
#define	rval	val[2]
#define	bcon	con[0]
#define	gcon	con[1]
#define	rcon	con[2]
typedef	struct {
	float	val[3];	// values
	float	con[3];	// confidences
} vc_t;

#define	MAXIMGS	32	// better not exceed this!
char	*inname[MAXIMGS];
Mat	imgin[MAXIMGS];
Mat	imgout;
Mat	warpin[MAXIMGS];
int	imgcnt;
float	xscale, yscale;
int	fpruned[MAXIMGS];
int	opruned;

void
mkgamma(float gamma)
{
	float scale = 0xffff / powf(255.0, gamma);
	for (int i=0; i<256; ++i) {
		gtable[i] = (scale * powf(((float) i), gamma)) + 0.5;
	}
}

inline ushort
cvt8to16(uchar c)
{
	// convert 8 to 16 bit
	ushort s = c;
	return(s * s);
}

inline ushort
cvt16(ushort s)
{
	// convert "fake" 16 bit to real 16 bit
	// assumes that bpp is number of valid bits
	return((s << bppl) | (s >> bppr));
}

inline int
cfasafe(int v, int max)
{
	if (v < 0) return(v & 1);
	if (v >= max) return((max-1) - (((max-1) & 1) != (v & 1)));
	return(v);
}

inline float
condist(float wantx, float wanty, float havex, float havey)
{
	// basically returns 1/(dist*dist), catching errors
	float difx = wantx - havex;
	float dify = wanty - havey;
	difx *= difx;
	dify *= dify;
	float dist = difx + dify;
	if (dist >= 4.0) return(0.0);
	if (dist < distlim) return(1.0);
	return(distlim / dist);
}

#define	HBITS	8
#define	HDIM	(1<<(HBITS))
#define	HVAL(V)	(((V)>>(16-HBITS))&(HDIM-1))

float	hist[HDIM][HDIM][3];

float
hcon(int ref, int val, int c)
{
	// return histogram-based confidence
	return(hist[HVAL(ref)][HVAL(val)][c]);
}

void
beginhist()
{

	// initialize hist to 0, but 1 on the edges
	for (int j=0; j<HDIM; ++j) {
		for (int i=0; i<HDIM; ++i) {
			for (int c=0; c<3; ++c) {
				hist[j][i][c] = ((i == 0) || (i == (HDIM-1)));
			}
		}
	}
}

void
mkhistrange(int v0, int v1, int c)
{
	// outline bounding box of range v0..v1 in color c
	// double counts if v0==v1, but that's OK

	int i0 = HVAL(v0);
	int i1 = HVAL(v1);
	for (int i=i0; i<=i1; ++i) {
#pragma omp atomic
		++hist[i][i0][c];
#pragma omp atomic
		++hist[i][i1][c];
	}
}

void
mkhist()
{
	// create histogram using similarity of pixel values
	// across all supplied images

	// fill histogram
	switch (imgin[0].elemSize()) {
	case 3:
#pragma omp parallel for schedule(guided)
		for (int y=0; y<yin; ++y) {
			for (int x=0; x<xin; ++x) {
				int rgb[MAXIMGS][3];
				uchar *p = imgin[0].ptr(y, x);
				rgb[0][0] = cvt8to16(p[0]);
				rgb[0][1] = cvt8to16(p[1]);
				rgb[0][2] = cvt8to16(p[2]);
				int k = 1;
				for (int i=1; i<imgcnt; ++i) {
					const float *wp = warpin[i].ptr<float>(0);
					int ix = ((int) ((wp[0] * x + wp[1] * y + wp[2]) + 0.5));
					int iy = ((int) ((wp[3] * x + wp[4] * y + wp[5]) + 0.5));
					if ((iy >= 0) &&
					    (iy < yin) &&
					    (ix >= 0) &&
					    (ix < xin)) {
						uchar *p = imgin[i].ptr(iy, ix);
						rgb[k][0] = cvt8to16(p[0]);
						rgb[k][1] = cvt8to16(p[1]);
						rgb[k][2] = cvt8to16(p[2]);
						++k;
					}
				}

				// all combinations
				for (int j=0; j<k; ++j) {
					for (int i=j+1; i<k; ++i) {
						mkhistrange(rgb[j][0], rgb[i][0], 0);
						mkhistrange(rgb[j][1], rgb[i][1], 0);
						mkhistrange(rgb[j][2], rgb[i][2], 0);
					}
				}
			}
		}
		break;
	case 6:
#pragma omp parallel for schedule(guided)
		for (int y=0; y<yin; ++y) {
			for (int x=0; x<xin; ++x) {
				int rgb[MAXIMGS][3];
				ushort *p = ((ushort *) imgin[0].ptr(y, x));
				rgb[0][0] = p[0];
				rgb[0][1] = p[1];
				rgb[0][2] = p[2];
				int k = 1;
				for (int i=1; i<imgcnt; ++i) {
					const float *wp = warpin[i].ptr<float>(0);
					int ix = ((int) ((wp[0] * x + wp[1] * y + wp[2]) + 0.5));
					int iy = ((int) ((wp[3] * x + wp[4] * y + wp[5]) + 0.5));
					if ((iy >= 0) &&
					    (iy < yin) &&
					    (ix >= 0) &&
					    (ix < xin)) {
						ushort *p = ((ushort *) imgin[i].ptr(y, x));
						rgb[k][0] = p[0];
						rgb[k][1] = p[1];
						rgb[k][2] = p[2];
						++k;
					}
				}

				// all combinations
				for (int j=0; j<k; ++j) {
					for (int i=j+1; i<k; ++i) {
						mkhistrange(rgb[j][0], rgb[i][0], 0);
						mkhistrange(rgb[j][1], rgb[i][1], 0);
						mkhistrange(rgb[j][2], rgb[i][2], 0);
					}
				}
			}
		}
		break;
	default:
#pragma omp parallel for schedule(guided)
		for (int y=0; y<yin; ++y) {
			for (int x=0; x<xin; ++x) {
				int rgb[2*MAXIMGS];
				int c = cof[((y & 1) << 1) | (x & 1)];
				ushort *p = ((ushort *) imgin[0].ptr(y, x));
				rgb[0] = *p;
				int k = 1;
				for (int i=1; i<imgcnt; ++i) {
					// takes nearest of same color channel...
					// but that could be TWO values for green
					const float *wp = warpin[i].ptr<float>(0);
					int ix = ((int) ((wp[0] * x + wp[1] * y + wp[2]) + 0.5));
					int iy = ((int) ((wp[3] * x + wp[4] * y + wp[5]) + 0.5));
					if ((iy >= 0) &&
					    (iy < yin) &&
					    (ix >= 0) &&
					    (ix < xin) &&
					    (c == cof[((iy & 1) << 1) | (ix & 1)])) {
						ushort *p = ((ushort *) imgin[i].ptr(y, x));
						rgb[k++] = *p;
					}
					++ix;
					if ((iy >= 0) &&
					    (iy < yin) &&
					    (ix >= 0) &&
					    (ix < xin) &&
					    (c == cof[((iy & 1) << 1) | (ix & 1)])) {
						ushort *p = ((ushort *) imgin[i].ptr(y, x));
						rgb[k++] = *p;
					}
					--ix; ++iy;
					if ((iy >= 0) &&
					    (iy < yin) &&
					    (ix >= 0) &&
					    (ix < xin) &&
					    (c == cof[((iy & 1) << 1) | (ix & 1)])) {
						ushort *p = ((ushort *) imgin[i].ptr(y, x));
						rgb[k++] = *p;
					}
					++ix;
					if ((iy >= 0) &&
					    (iy < yin) &&
					    (ix >= 0) &&
					    (ix < xin) &&
					    (c == cof[((iy & 1) << 1) | (ix & 1)])) {
						ushort *p = ((ushort *) imgin[i].ptr(y, x));
						rgb[k++] = *p;
					}
				}

				// all combinations
				for (int j=0; j<k; ++j) {
					for (int i=j+1; i<k; ++i) {
						mkhistrange(rgb[j], rgb[i], c);
					}
				}
			}
		}
	}
}

void
endhist()
{

	// make histogram monotonic
#pragma omp parallel for schedule(guided)
	for (int j=0; j<HDIM; ++j) {
		float v[3];
		v[2] = v[1] = v[0] = 0;
		for (int i=0; i<=j; ++i) {
			for (int c=0; c<3; ++c) {
				v[c] = (hist[j][i][c] += v[c]);
			}
		}
	}

#pragma omp parallel for schedule(guided)
	for (int j=0; j<HDIM; ++j) {
		float v[3];
		v[2] = v[1] = v[0] = 0;
		for (int i=HDIM-1; i>=j; --i) {
			for (int c=0; c<3; ++c) {
				v[c] = (hist[j][i][c] += v[c]);
			}
		}
	}

	// normalize histogram so 100% confidence
	// of value being itself
#pragma omp parallel for schedule(guided)
	for (int j=0; j<HDIM; ++j) {
		float v[3];
		v[0] = hist[j][j][0];
		v[1] = hist[j][j][1];
		v[2] = hist[j][j][2];
		for (int i=0; i<HDIM; ++i) {
			for (int c=0; c<3; ++c) {
				hist[j][i][c] = powf(hist[j][i][c] / v[c], histexp);
			}
		}
	}

	if (emfile) {
		// write out noise model as an image
		Mat errmod = Mat(HDIM, HDIM, CV_8UC3);
#pragma omp parallel for schedule(guided)
		for (int j=0; j<HDIM; ++j) {
			float *fp = &(hist[j][0][0]);
			uchar *ucp = errmod.ptr(j);
			for (int i=0; i<HDIM; ++i) {
				*(ucp++) = ((uchar) (255 * *(fp++)));
				*(ucp++) = ((uchar) (255 * *(fp++)));
				*(ucp++) = ((uchar) (255 * *(fp++)));
			}
		}
		imwrite(emfile, errmod);
		if (verbose) {
			fprintf(stderr,
				"writing pixel error model to file %s\n",
				emfile);
		}
	}
}

vc_t
getvc(int x, int y, int imgno)
{
	// synthesize red, green, blue value & confidence for pixel at output x,y
	// uses weighted combination of nearest of same color
	// note that val is multiplied by the con
	vc_t v;

	// convert to local image coordinates
	const float *wp = warpin[imgno].ptr<float>(0);
	float sx = (x * xscale);
	float sy = (y * yscale);
	float fx = wp[0] * sx + wp[1] * sy + wp[2];
	float fy = wp[3] * sx + wp[4] * sy + wp[5];

	// convert to nearest pixel coordinates
	int nx = ((int) fx);
	int ny = ((int) fy);

	// find each color in neighborhood
	v.rval = v.gval = v.bval = 0;
	v.rcon = v.gcon = v.bcon = 0;

	switch (imgin[imgno].elemSize()) {
	case 3:
		// three 8-bit colors per pixel
		if ((ny>=0) && (ny<yin) && (nx>=0) && (nx<xin)) {
			v.rcon = v.gcon = v.bcon = condist(fx, fy, nx, ny);
			uchar *p = imgin[imgno].ptr(ny, nx);
			v.bval = cvt8to16(p[0]) * v.bcon;
			v.gval = cvt8to16(p[1]) * v.gcon;
			v.rval = cvt8to16(p[2]) * v.rcon;
		}
		break;
	case 6:
		// three 16-bit colors per pixel
		if ((ny>=0) && (ny<yin) && (nx>=0) && (nx<xin)) {
			v.rcon = v.gcon = v.bcon = condist(fx, fy, nx, ny);
			ushort *p = ((ushort *) imgin[imgno].ptr(ny, nx));
			v.bval = p[0] * v.bcon;
			v.gval = p[1] * v.gcon;
			v.rval = p[2] * v.rcon;
		}
		break;
	default:
		// CFA pattern on raw
		if (smooth) {
			for (int j=ny-1; j<=ny+1; ++j) if ((j>=0) && (j<yin)) {
				int patj = ((j & 1) << 1);
				for (int i=nx-1; i<=nx+1; ++i) if ((i>=0) && (i<xin)) {
					int pat = patj + (i & 1);
					int c = cof[pat];
					float d = condist(fx, fy, i, j);
					v.con[c] += d;
					v.val[c] += (d * imgin[imgno].at<pixel_t>(j, i));
				}
			}
		} else {
			for (int j=ny-1; j<=ny+1; ++j) if ((j>=0) && (j<yin)) {
				int patj = ((j & 1) << 1);
				for (int i=nx-1; i<=nx+1; ++i) if ((i>=0) && (i<xin)) {
					int pat = patj + (i & 1);
					int c = cof[pat];
					float d = condist(fx, fy, i, j);
					if (d > v.con[c]) {
						v.con[c] = d;
						v.val[c] = (d * imgin[imgno].at<pixel_t>(j, i));
					}
				}
			}
		}
	}

	return(v);
}

VecTYP
getv(int x, int y)
{
	// get vector of color channel values in OpenCV BGR order
	// values are synthesized using value & confidence metrics
	vc_t valcon[MAXIMGS];
	vc_t uvalcon[MAXIMGS]; // unscaled value version
	double rgbdist[MAXIMGS];

	// get value, confidence from each image
	// and compute weighted average
	vc_t va;
	for (int i=0; i<imgcnt; ++i) {
		uvalcon[i] = (valcon[i] = getvc(x, y, i));

		// undo scaling val by con, but don't divide by 0
		if (uvalcon[i].bcon) uvalcon[i].bval /= valcon[i].bcon;
		if (uvalcon[i].gcon) uvalcon[i].gval /= valcon[i].gcon;
		if (uvalcon[i].rcon) uvalcon[i].rval /= valcon[i].rcon;
	}

	if (first > 0) {
		// first is weight for being close to first image value
		for (int i=0; i<imgcnt; ++i) {
			if (uvalcon[i].bcon) {
				uvalcon[i].bcon += first * hcon(((int) uvalcon[i].bval), ((int) uvalcon[0].bval), 0);
			}
			if (uvalcon[i].gcon) {
				uvalcon[i].gcon += first * hcon(((int) uvalcon[i].gval), ((int) uvalcon[0].gval), 1);
			}
			if (uvalcon[i].rcon) {
				uvalcon[i].rcon += first * hcon(((int) uvalcon[i].rval), ((int) uvalcon[0].rval), 2);
			}
		}
	} else if (first < 0) {
		// first is pruning threshhold
		// first of -v means remove any channel that isn't
		// confidence weighting that isn't at least v
		for (int i=0; i<imgcnt; ++i) {
			if (uvalcon[i].bcon) {
				if (hcon(((int) uvalcon[i].bval), ((int) uvalcon[0].bval), 0) < -first) {
					uvalcon[i].bcon = 0;
#pragma omp atomic
					++fpruned[i];
				}
			}
			if (uvalcon[i].gcon) {
				if (hcon(((int) uvalcon[i].gval), ((int) uvalcon[0].gval), 1) < -first) {
					uvalcon[i].gcon = 0;
#pragma omp atomic
					++fpruned[i];
				}
			}
			if (uvalcon[i].rcon) {
				if (hcon(((int) uvalcon[i].rval), ((int) uvalcon[0].rval), 2) < -first) {
					uvalcon[i].rcon = 0;
#pragma omp atomic
					++fpruned[i];
				}
			}
		}
	}

	if (!outliers) {
		// don't filter outliers
		va.rval = va.gval = va.bval = 0;
		va.rcon = va.gcon = va.bcon = 0;
		for (int i=0; i<imgcnt; ++i) {
			// add impacts of scaled values
			// which may have been altered by first filter
			va.rval += uvalcon[i].rval * uvalcon[i].rcon;
			va.rcon += uvalcon[i].rcon;
			va.gval += uvalcon[i].gval * uvalcon[i].gcon;
			va.gcon += uvalcon[i].gcon;
			va.bval += uvalcon[i].bval * uvalcon[i].bcon;
			va.bcon += uvalcon[i].bcon;
		}
	} else {
		// try to remove outlier values...
		// outlier confidences are weighted by outliers factor

		// bubble sort into decreasing value
		// with 0 confidences at the bottom
		for (int k=0; k<imgcnt; ++k) {
			for (int j=k+1; j<imgcnt; ++j) {
				if ((uvalcon[k].rcon == 0) ||
				    (uvalcon[k].rval < uvalcon[j].rval)) {
					float c = uvalcon[k].rcon;
					float v = uvalcon[k].rval;
					uvalcon[k].rcon = uvalcon[j].rcon;
					uvalcon[k].rval = uvalcon[j].rval;
					uvalcon[j].rcon = c;
					uvalcon[j].rval = v;
				}
				if ((uvalcon[k].gcon == 0) ||
				    (uvalcon[k].gval < uvalcon[j].gval)) {
					float c = uvalcon[k].gcon;
					float v = uvalcon[k].gval;
					uvalcon[k].gcon = uvalcon[j].gcon;
					uvalcon[k].gval = uvalcon[j].gval;
					uvalcon[j].gcon = c;
					uvalcon[j].gval = v;
				}
				if ((uvalcon[k].bcon == 0) ||
				    (uvalcon[k].bval < uvalcon[j].bval)) {
					float c = uvalcon[k].bcon;
					float v = uvalcon[k].bval;
					uvalcon[k].bcon = uvalcon[j].bcon;
					uvalcon[k].bval = uvalcon[j].bval;
					uvalcon[j].bcon = c;
					uvalcon[j].bval = v;
				}
			}
		}

		// find first 0 confidence entry
		int bgrend[3];
		bgrend[2] = bgrend[1] = bgrend[0] = 0;
		for (int i=0; i<imgcnt; ++i) {
			if (uvalcon[i].bcon > 0) bgrend[0] = i;
			if (uvalcon[i].gcon > 0) bgrend[1] = i;
			if (uvalcon[i].rcon > 0) bgrend[2] = i;
		}

		// let's assume that the middle entry is our best guess
		// if so, we should be weighting relative to that
		float bv = uvalcon[bgrend[0] / 2].bval;
		if (!(bgrend[0] & 1)) bv = (bv + uvalcon[bgrend[0] / 2 + 1].bval) / 2.0;
		float gv = uvalcon[bgrend[1] / 2].gval;
		if (!(bgrend[1] & 1)) gv = (gv + uvalcon[bgrend[1] / 2 + 1].gval) / 2.0;
		float rv = uvalcon[bgrend[2] / 2].rval;
		if (!(bgrend[2] & 1)) rv = (rv + uvalcon[bgrend[2] / 2 + 1].rval) / 2.0;

		// perform confidence-weighted sum
		va.rval = va.gval = va.bval = 0;
		va.rcon = va.gcon = va.bcon = 0;
		for (int j=0; j<imgcnt; ++j) {
			float con;
			con = hcon(((int) uvalcon[j].bval), ((int) bv), 0);
			if (outliers < 0) {
				if (con >= -outliers) {
					// leave it in if it clears outliers threshold
					va.bval += (uvalcon[j].bval * uvalcon[j].bcon);
					va.bcon += uvalcon[j].bcon;
				} else {
#pragma omp atomic
					++opruned;
				}
			} else if (uvalcon[j].bcon) {
				// weight it by outliers factor
				con = (con * outliers) + uvalcon[j].bcon;
				va.bval += (uvalcon[j].bval * con);
				va.bcon += con;
			}
			con = hcon(((int) uvalcon[j].gval), ((int) gv), 1);
			if (outliers < 0) {
				if (con >= -outliers) {
					// leave it in if it clears outliers threshold
					va.gval += (uvalcon[j].gval * uvalcon[j].gcon);
					va.gcon += uvalcon[j].gcon;
				} else {
#pragma omp atomic
					++opruned;
				}
			} else if (uvalcon[j].gcon) {
				// weight it by outliers factor
				con = (con * outliers) + uvalcon[j].gcon;
				va.gval += (uvalcon[j].gval * con);
				va.gcon += con;
			}
			con = hcon(((int) uvalcon[j].rval), ((int) rv), 2);
			if (outliers < 0) {
				if (con >= -outliers) {
					// leave it in if it clears outliers threshold
					va.rval += (uvalcon[j].rval * uvalcon[j].rcon);
					va.rcon += uvalcon[j].rcon;
				} else {
#pragma omp atomic
					++opruned;
				}
			} else if (uvalcon[j].rcon) {
				// weight it by outliers factor
				con = (con * outliers) + uvalcon[j].rcon;
				va.rval += (uvalcon[j].rval * con);
				va.rcon += con;
			}
		}
	}

	VecTYP res;
	res[0] = ((int) ((va.bval / va.bcon) + 0.5));
	res[1] = ((int) ((va.gval / va.gcon) + 0.5));
	res[2] = ((int) ((va.rval / va.rcon) + 0.5));
	return(res);
}

VecTYP
getneigh(int x, int y, float bavg, float gavg, float ravg)
{
	// second pass of getv that works a lot like OUTLIERS, but
	// the second pass looks for the colors most similar to
	// the input averages rarther than the middle value
	// values are synthesized using value & confidence metrics
	vc_t valcon[MAXIMGS];

	// get value, confidence from each image
	// and compute weighted average
	vc_t va;
	valcon[0] = (va = getvc(x, y, 0));
	for (int i=1; i<imgcnt; ++i) {
		valcon[i] = getvc(x, y, i);
		va.rval += valcon[i].rval;
		va.rcon += valcon[i].rcon;
		va.gval += valcon[i].gval;
		va.gcon += valcon[i].gcon;
		va.bval += valcon[i].bval;
		va.bcon += valcon[i].bcon;
	}

	// try to remove values too different from average...
	// this is really about reducing motion artifacts
	vc_t uvalcon[MAXIMGS];

	// make temporary copy of valcon[] with vals unscaled
	for (int i=0; i<imgcnt; ++i) {
		uvalcon[i] = valcon[i];
		if (uvalcon[i].rcon) uvalcon[i].rval /= valcon[i].rcon;
		if (uvalcon[i].gcon) uvalcon[i].gval /= valcon[i].gcon;
		if (uvalcon[i].bcon) uvalcon[i].bval /= valcon[i].bcon;
	}
printf("Bubbling\n");
	// bubble sort 'em into decreasing value
	// with 0 confidences at the bottom
	for (int k=0; k<imgcnt; ++k) {
		for (int j=k+1; j<imgcnt; ++j) {
			if ((uvalcon[k].rcon == 0) ||
			    (uvalcon[k].rval < uvalcon[j].rval)) {
				float c = uvalcon[k].rcon;
				float v = uvalcon[k].rval;
				uvalcon[k].rcon = uvalcon[j].rcon;
				uvalcon[k].rval = uvalcon[j].rval;
				uvalcon[j].rcon = c;
				uvalcon[j].rval = v;
			}
			if ((uvalcon[k].gcon == 0) ||
			    (uvalcon[k].gval < uvalcon[j].gval)) {
				float c = uvalcon[k].gcon;
				float v = uvalcon[k].gval;
				uvalcon[k].gcon = uvalcon[j].gcon;
				uvalcon[k].gval = uvalcon[j].gval;
				uvalcon[j].gcon = c;
				uvalcon[j].gval = v;
			}
			if ((uvalcon[k].bcon == 0) ||
			    (uvalcon[k].bval < uvalcon[j].bval)) {
				float c = uvalcon[k].bcon;
				float v = uvalcon[k].bval;
				uvalcon[k].bcon = uvalcon[j].bcon;
				uvalcon[k].bval = uvalcon[j].bval;
				uvalcon[j].bcon = c;
				uvalcon[j].bval = v;
			}
		}
	}
printf("Sorted\n");
	// find first 0 confidence entry
	int rend = 0, gend = 0, bend = 0;
	for (int i=0; i<imgcnt; ++i) {
		rend += (uvalcon[i].rcon > 0);
		gend += (uvalcon[i].gcon > 0);
		bend += (uvalcon[i].bcon > 0);
	}

	// find middle value for each as value
	// closest to average, resolve ties by
	// taking the entry with highest confidence
	int rmid = 0, gmid = 0, bmid = 0;
	float absdif;
	absdif = fabsf(uvalcon[rmid].rval - ravg);
	for (int i=1; i<rend; ++i) {
		float t = fabsf(uvalcon[i].rval - ravg);
		if ((absdif > t) ||
		    ((absdif == t) && (uvalcon[rmid].rcon < uvalcon[i].rcon))) {
			rmid = i;
			absdif = t;
		}
	}
	absdif = fabsf(uvalcon[gmid].gval - gavg);
	for (int i=1; i<gend; ++i) {
		float t = fabsf(uvalcon[i].gval - gavg);
		if ((absdif > t) ||
		    ((absdif == t) && (uvalcon[gmid].gcon < uvalcon[i].gcon))) {
			gmid = i;
			absdif = t;
		}
	}
	absdif = fabsf(uvalcon[bmid].bval - bavg);
	for (int i=1; i<bend; ++i) {
		float t = fabsf(uvalcon[i].bval - bavg);
		if ((absdif > t) ||
		    ((absdif == t) && (uvalcon[bmid].bcon < uvalcon[i].bcon))) {
			bmid = i;
			absdif = t;
		}
	}
printf("Pruning...\n");
	// for each color, keep adding most-similar samples
	// to mid until confidence > 1/2 full confidence
	// then recompute va using only lo..hi entries
	int hi, lo;
	float consum;
	// red channel
	hi = (lo = rmid);
	consum = uvalcon[rmid].rcon;
	while ((2.0 * consum) < va.rcon) {
		if (lo > 0) {
			if (hi < (rend-1)) {
				if (fabsf(uvalcon[lo].rval - uvalcon[lo-1].rval) <
				    fabsf(uvalcon[hi].rval - uvalcon[hi+1].rval)) {
					consum += uvalcon[--lo].rcon;
				} else {
					consum += uvalcon[++hi].rcon;
				}
			} else consum += uvalcon[--lo].rcon;
		} else if (hi < (rend-1)) consum += uvalcon[++hi].rcon;
		else break;
	}
	va.rval = (va.rcon = 0);
	for (int i=lo; i<=hi; ++i) {
		va.rval += (uvalcon[i].rval * uvalcon[i].rcon);
		va.rcon += uvalcon[i].rcon;
	}
	// green channel
	hi = (lo = gmid);
	consum = uvalcon[gmid].gcon;
	while ((2.0 * consum) < va.gcon) {
		if (lo > 0) {
			if (hi < (gend-1)) {
				if (fabsf(uvalcon[lo].gval - uvalcon[lo-1].gval) <
				    fabsf(uvalcon[hi].gval - uvalcon[hi+1].gval)) {
					consum += uvalcon[--lo].gcon;
				} else {
					consum += uvalcon[++hi].gcon;
				}
			} else consum += uvalcon[--lo].gcon;
		} else if (hi < (gend-1)) consum += uvalcon[++hi].gcon;
		else break;
	}
	va.gval = (va.gcon = 0);
	for (int i=lo; i<=hi; ++i) {
		va.gval += (uvalcon[i].gval * uvalcon[i].gcon);
		va.gcon += uvalcon[i].gcon;
	}
	// blue channel
	hi = (lo = bmid);
	consum = uvalcon[bmid].bcon;
	while ((2.0 * consum) < va.bcon) {
		if (lo > 0) {
			if (hi < (bend-1)) {
				if (fabsf(uvalcon[lo].bval - uvalcon[lo-1].bval) <
				    fabsf(uvalcon[hi].bval - uvalcon[hi+1].bval)) {
					consum += uvalcon[--lo].bcon;
				} else {
					consum += uvalcon[++hi].bcon;
				}
			} else consum += uvalcon[--lo].bcon;
		} else if (hi < (bend-1)) consum += uvalcon[++hi].bcon;
		else break;
	}
	va.bval = (va.bcon = 0);
	for (int i=lo; i<=hi; ++i) {
		va.bval += (uvalcon[i].bval * uvalcon[i].bcon);
		va.bcon += uvalcon[i].bcon;
	}

	VecTYP res;
	res[0] = ((int) ((va.bval / va.bcon) + 0.5));
	res[1] = ((int) ((va.gval / va.gcon) + 0.5));
	res[2] = ((int) ((va.rval / va.rcon) + 0.5));
	return(res);
}

int
has_cfa()
{
	return(strcmp(cfa, NO_CFA) != 0);
}

void
decode_cfa(const char *ref)
{
	// set global variables for color filter array pattern
	// the character string
	cfa = ref;

	// decode CFA offsets
	cfar = cfag = cfaG = cfab = 4;
	for (int i=0; i<4; ++i) {
		switch (cfa[i]) {
		case 'R': case 'r': if (cfar != 4) goto badcfa; cfar = i; break;
		case 'B': case 'b': if (cfab != 4) goto badcfa; cfab = i; break;
		case 'G': case 'g': if (cfag == 4) cfag = i; else {
					if (cfaG != 4) goto badcfa;
					cfaG = i; }
				    break;
		default:
badcfa:
			fprintf(stderr,
				"%s: invalid CFA specification '%s'\n",
				prog, cfa);
			exit(2);
		}
	}
	cof[cfab] = 0;
	cof[cfag] = cof[cfaG] = 1;
	cof[cfar] = 2;
}

void
usage()
{
	fprintf(stderr,
"Usage: %s -options input_files\n"
"-a int\tset number of alignment iterations (%d)\n"
"-c str\t2x2 RGB CFA pattern is %s\n"
"-d    \tuse DCRAW instead of LibRaw (%s)\n"
"-e flt\talignment termination epsilon (%g)\n"
"-E str\twrite pixel error model file (%s)\n"
"-f flt\tfilter distance from first with this weighting (%g)\n"
"      \tpositive amplifies confidence, negative prunes\n"
"-g flt\tset gamma for 8-bit input pixels (%g)\n"
"-h flt\thistogram noise model exponent, 0 disables (%g)\n"
"-i    \ttoggle interpolation within input images (%s)\n"
"-m flt\tset maximum confidence (%g)\n"
"-n    \ttoggle neighborhood filtering (%s)\n"
"-o str\tset output file name (%s)\n"
"-O flt\tfilter outliers with this weighting (%g)\n"
"      \tpositive amplifies confidence, negative prunes\n"
"-p flt\tprune image if alignment rotation/scaling exceeds (%lf)\n"
"-r    \ttoggle input images are raw (%s)\n"
"-s    \ttoggle smoothing in 24bpp to raw conversion (%s)\n"
"-v    \tincrement verbosity of messages\n"
"-X int\tset final image int times input image x (columns); 2 by default\n"
"-x int\tset final image x as int (columns)\n"
"-Y int\tset final image int times input image y (rows); 2 by default\n"
"-y int\tset final image y as int (rows)\n",
		prog,
		aligniter,
		cfa,
		(dcraw ? "ON" : "OFF"),
		epsilon,
		(emfile ? emfile : "not written"),
		first,
		ingamma,
		histexp,
		(interp ? "ON" : "OFF"),
		maxcon,
		(neigh ? "ON" : "OFF"),
		outname,
		outliers,
		prune,
		(rawin ? "ON" : "OFF"),
		(smooth ? "ON" : "OFF"));
}


int
main(int argc, char **argv)
{
	// set argument defaults
	prog = argv[0];
	cfa = NO_CFA;
	mkgamma(ingamma);	// default gamma for 8-bit input pixels

	// parse arguments
	for (int i=1; i<argc; ++i) {
		char *ap = argv[i];
		switch (*ap) {
		case '-':
			switch (ap[1]) {
			case 'a': case 'A':
				// set alignment iteration count
				aligniter = atoi(argv[i+1]);
				++i;
				break;
			case 'c': case 'C':
				// set CFA pattern
				decode_cfa(argv[i+1]);
				++i;
				break;
			case 'd': case 'D':
				dcraw = !dcraw;
				break;
			case 'e':
				// set alignment epsilon
				epsilon = atof(argv[i+1]);
				++i;
				break;
			case 'E':
				// set error model file name
				emfile = argv[i+1];
				++i;
				break;
			case 'f': case 'F':
				// do first image filtering?
				first = atof(argv[i+1]);
				++i;
				break;
			case 'g': case 'G':
				// set gamma for 8-bit input
				ingamma = atof(argv[i+1]);
				mkgamma(ingamma);
				++i;
				break;
			case 'h': case 'H':
				// histogram noise model exponent, 0 to disable
				histexp = atof(argv[i+1]);
				++i;
				break;
			case 'i': case 'I':
				interp = !interp;
				break;
			case 'm': case 'M':
				// maximum confidence
				maxcon = atof(argv[i+1]);
				invmaxcon = 1.0 / maxcon;
				++i;
				break;
			case 'n': case 'N':
				// do neighborhood filtering?
				neigh = !neigh;
				break;
			case 'o':
				// set output file name
				outname = argv[i+1];
				++i;
				break;
			case 'O':
				// do outlier filtering?
				outliers = atof(argv[i+1]);
				++i;
				break;
			case 'p': case 'P':
				// threshold for pruning rotation/scaling
				prune = atof(argv[i+1]);
				++i;
				break;
			case 'r': case 'R':
				// inputs are raw images
				rawin = !rawin;
				break;
			case 's': case 'S':
				// smooth while creating raw from 24bpp
				smooth = !smooth;
				break;
			case 'v': case 'V':
				// verbose?
				++verbose;
				break;
			case 'x':
				// set final X dimension
				xout = atoi(argv[i+1]);
				++i;
				break;
			case 'X':
				// set final X dimension multiplier
				xout = -atoi(argv[i+1]);
				++i;
				break;
			case 'y':
				// set final X dimension
				yout = atoi(argv[i+1]);
				++i;
				break;
			case 'Y':
				// set final X dimension multiplier
				yout = -atoi(argv[i+1]);
				++i;
				break;
			default:
				usage();
				exit(1);
			}
			break;
		default:
			// should be an input image
			if (imgcnt >= MAXIMGS) {
				fprintf(stderr, "%s: too many input images (max %d)\n", prog, MAXIMGS);
				exit(2);
			}

			// read the file
			inname[imgcnt] = ap;

			// note: this reads 16-bit pixels into 8-bit values
			// and color into grayscale...
			// but the alignment computation needs it that way
			if (rawin) {
				// make a very raw 16bpp PGM version
re_read:
				if (!has_cfa()) decode_cfa("RGGB");
				char buf[4096];
				if (dcraw) {
					sprintf(buf, "dcraw -4 -D -c %s >%s.pgm", inname[imgcnt], inname[imgcnt]);
				} else {
					sprintf(buf, "unprocessed_raw -q %s", inname[imgcnt]);
				}
				system(buf);
				sprintf(buf, "%s.pgm", inname[imgcnt]);
				imgin[imgcnt] = imread(buf, IMREAD_GRAYSCALE);
				// get rid of raw PGM file
				unlink(buf);
			} else {
				// read in the image as native depth and colors
				imgin[imgcnt] = imread(inname[imgcnt], IMREAD_GRAYSCALE);
				if ((imgcnt == 0) && imgin[imgcnt].empty()) {
					// maybe it should be raw?
					if (verbose) {
						fprintf(stderr,
							"trying to read %s as raw\n",
							inname[imgcnt]);
					}
					rawin = 1;
					goto re_read;
				}
			}
			if (imgin[imgcnt].empty()) {
				fprintf(stderr,
					"%s: could not read %s\n",
					prog, inname[imgcnt]);
				exit(3);
			}
			if (verbose) {
				fprintf(stderr,
					"%s: %dx%d\n",
					inname[imgcnt],
					imgin[imgcnt].cols, imgin[imgcnt].rows);
			}

			// all input images must be same size
			if ((imgin[imgcnt].cols != imgin[0].cols) ||
			    (imgin[imgcnt].rows != imgin[0].rows)) {
				fprintf(stderr,
					"%s: dimensions of %s and %s differ\n",
					prog, inname[imgcnt], inname[0]);
				exit(4);
			}

			// compute alignment with first input image
			warpin[imgcnt] = Mat::eye(2, 3, CV_32F);
			if (imgcnt == 0) {
				if (verbose) {
					printf("%s is the reference image for alignment\n",
					       inname[imgcnt]);
				}
				++imgcnt;
			} else {
				const double start_time  = (double) getTickCount();
				const int warp_mode = MOTION_EUCLIDEAN;
				TermCriteria criteria (TermCriteria::COUNT+TermCriteria::EPS, aligniter, epsilon);
				double cc = findTransformECC(imgin[0], imgin[imgcnt], warpin[imgcnt], warp_mode, criteria);
				const double end_time  = (double) getTickCount();
				const double total_time = (end_time - start_time) / getTickFrequency();
				imgin[imgcnt].release(); // don't need this anymore...
				if (cc == -1) {
					fprintf(stderr,
						"%s: image %s pruned because alignment with %s failed\n",
						prog, inname[imgcnt], inname[0]);
				} else {
					// are the alignments found OK?
					const float* wp = warpin[imgcnt].ptr<float>(0);
					const double t = (fabsf(wp[0]-1)+fabsf(wp[1])+fabsf(wp[3])+fabsf(wp[4]-1));

					if (verbose) {
						printf("alignment of %s is translation %g,%g with %g correlation in %gs\n",
						       inname[imgcnt], wp[2], wp[5], cc, total_time);
						if (verbose > 1) {
							printf("complete alignment matrix is:\n"
							       "%g\t%g\t%g\n"
							       "%g\t%g\t%g\n",
							       wp[0], wp[1], wp[2],
							       wp[3], wp[4], wp[5]);
						}
					}

					if (t > prune) {
						// we will apply a bit more than just pure translation,
						// but empirically quality drops sharply if this is too much
						// so it is better to just skip the image
						fprintf(stderr,
							"%s: image %s pruned due to %g difference from pure translation\n",
							prog, inname[imgcnt], t);
					} else {
						// keep this image
						++imgcnt;
					}
				}
			}
		}
	}

	if (imgcnt < 1) {
		fprintf(stderr, "%s: no usable input images\n", prog);
		usage();
		exit(2);
	}
	imgin[0].release(); // don't need this anymore

	// now re-read everything as native type,
	// which could be 16 bits per pixel monochrome
	for (int i=0; i<imgcnt; ++i) {
		if (rawin) {
			// make a very raw 16bpp PGM version
			char buf[4096];
			if (dcraw) {
				sprintf(buf, "dcraw -4 -D -c %s >%s.pgm", inname[i], inname[i]);
			} else {
				sprintf(buf, "unprocessed_raw -q %s", inname[i]);
			}
			system(buf);
			sprintf(buf, "%s.pgm", inname[i]);
			imgin[i] = imread(buf, IMREAD_UNCHANGED);
			// get rid of raw PGM file
			unlink(buf);
		} else {
			// read in the image as native depth and colors
			imgin[i] = imread(inname[i], IMREAD_UNCHANGED);
		}
		if (verbose) {
			printf("%s read as %dx%d %ldbpp\n",
			       inname[i], imgin[i].cols, imgin[i].rows, 8*imgin[i].elemSize());
		}

		// make it into a 16-bit monochrome image with all bits active
		switch (imgin[i].elemSize()) {
		case 1: {
				// 8-bit monochrome input converted to 16-bit
				Mat t = Mat(imgin[i].rows, imgin[i].cols, CV_16UC1);
				bpp = 8;
#pragma omp parallel for schedule(guided)
				for (int k=0; k<imgin[i].rows; ++k) {
					for (int j=0; j<imgin[i].cols; ++j) {
						pixel_t v = imgin[i].at<uchar>(k, j);
						imgin[i].at<pixel_t>(k, j) = cvt8to16(v);
					}
				}
				imgin[i] = t;
			} break;
		case 2: {
				// 16-bit, but how many bits are really valid?
				if (!bpp) {
					// initial guess is 12 valid bits per pixel
					bpp = 12;
newbpp:
					bppl = (16 - bpp);		// left shift to 16 bit
					bppr = (bpp - (16 - bpp));	// right shift to 16 bit
					for (int k=0; k<imgin[i].rows; ++k) {
						for (int j=0; j<imgin[i].cols; ++j) {
							pixel_t v = imgin[i].at<pixel_t>(k, j);
							if (v > ((1<<bpp)-1)) {
								do ++bpp; while (v > ((1<<bpp)-1));
								goto newbpp;
							}
						}
					}
					if (verbose) {
						printf("original data is %d bits per pixel\n", bpp);
					}
				}

				// now fix 'em
#pragma omp parallel for schedule(guided)
				for (int k=0; k<imgin[i].rows; ++k) {
					for (int j=0; j<imgin[i].cols; ++j) {
						pixel_t v = imgin[i].at<pixel_t>(k, j);
						imgin[i].at<pixel_t>(k, j) = cvt16(v);
					}
				}
			} break;
		case 3: if (has_cfa()) {
				// 24-bit color input converted to 16-bit monochrome
				// using the CFA pattern
				Mat t = Mat(imgin[i].rows, imgin[i].cols, CV_16UC1);

				if (verbose) {
					printf("convertng %s data from 24bpp color to 16bpp raw\n",
					       inname[i]);
				}

				if (smooth) {
					// offset for row above or below
					int ystep = imgin[i].cols * 3;
#pragma omp parallel for schedule(guided)
					for (int k=0; k<imgin[i].rows; ++k) {
						uchar *p = imgin[i].ptr(k);
						pixel_t kb = ((k & 1) << 1);
						for (int j=0; j<imgin[i].cols; ++j) {
							int v = p[ cof[kb | (j & 1)] ];
							v = 8 * cvt8to16(v);
							int w = 8;
							if (k > 0) {
								int va = p[cof[kb | (j & 1)] - ystep];
								v += cvt8to16(va);
								++w;
							}
							if (j > 0) {
								int va = p[cof[kb | (j & 1)] - 3];
								v += cvt8to16(va);
								++w;
							}
							if (j < imgin[i].cols-1) {
								int va = p[cof[kb | (j & 1)] + 3];
								v += cvt8to16(va);
								++w;
							}
							if (k < imgin[i].rows-1) {
								int va = p[cof[kb | (j & 1)] + ystep];
								v += cvt8to16(va);
								++w;
							}
							t.at<pixel_t>(k, j) = v / w;
							p += 3;
						}
					}
					imgin[i] = t;
				} else {
#pragma omp parallel for schedule(guided)
					for (int k=0; k<imgin[i].rows; ++k) {
						uchar *p = imgin[i].ptr(k);
						pixel_t kb = ((k & 1) << 1);
						for (int j=0; j<imgin[i].cols; ++j) {
							pixel_t v = p[ cof[kb | (j & 1)] ];
							t.at<pixel_t>(k, j) = cvt8to16(v);
							p += 3;
						}
					}
					imgin[i] = t;
				}
			} break;
		case 6: if (has_cfa()) {
				// 48-bit color input converted to 16-bit monochrome
				// using the CFA pattern
				Mat t = Mat(imgin[i].rows, imgin[i].cols, CV_16UC1);

				if (verbose) {
					printf("convertng %s data from 48bpp color to 16bpp raw\n",
					       inname[i]);
				}

#pragma omp parallel for schedule(guided)
				for (int k=0; k<imgin[i].rows; ++k) {
					ushort *p = ((ushort *) (imgin[i].ptr(k)));
					pixel_t kb = ((k & 1) << 1);
					for (int j=0; j<imgin[i].cols; ++j) {
						pixel_t v = p[ cof[kb | (j & 1)] ];
						t.at<pixel_t>(k, j) = cvt8to16(v);
						p += 3;
					}
				}
				imgin[i] = t;
			} break;
		default:
			fprintf(stderr,
				"could not handle %s as %ldbpp\n",
				inname[i], 8*imgin[i].elemSize());
			exit(6);
		}
	}

	// compute scaling
	xin = imgin[0].cols;
	yin = imgin[0].rows;
	if (xout < 1) {
		xout = (xout ? -xout * xin : 2 * xin);
	}
	if (yout < 1) {
		yout = (yout ? -yout * yin : 2 * yin);
	}
	xscale = xin / ((float) xout);
	yscale = yin / ((float) yout);
	if (verbose) {
		printf("output image %s will be %dx%d pixels\n",
		       outname, xout, yout);
	}
	distlim = (1.0 / ((xscale > yscale) ? xscale : yscale));
	distlim *= distlim;

	if (imgcnt < 2) {
		// turn off everything that needs more than  one image
		first = 0;
		histexp = 0;
		outliers = 0;
	}

	if (verbose > 1) {
		fprintf(stderr,
"All option settings used:\n"
"%d alignment iterations\n"
"CFA pattern is %s\n"
"Use %s for raw decoding\n"
"Alignment termination epsilon is %g\n"
"Pixel error model file is %s\n"
"Filter distance from first with this weighting (%g)\n"
"Gamma for input 8-bit pixel data is %g\n"
"Histogram noise model using %g exponent\n"
"Interpolation within input images is %s\n"
"Maximum confidence is %g\n"
"Neighborhood filtering is %s\n"
"Output file name is %s\n"
"Filter outliers with this weighting (%g)\n"
"Prune image if alignment rotation/scaling exceeds %lf\n"
"Input images are %sraw\n"
"Smoothing is %s\n"
"Verbosity level for messages is %d\n"
"Final image dimensions (columns, rows) are %dx%d\n",
		aligniter,
		cfa,
		(dcraw ? "DCRAW" : "LibRaw"),
		epsilon,
		(emfile ? emfile : "not written"),
		first,
		ingamma,
		histexp,
		(interp ? "ON" : "OFF"),
		maxcon,
		(neigh ? "ON" : "OFF"),
		outname,
		outliers,
		prune,
		(rawin ? "" : "NOT "),
		(smooth ? "ON" : "OFF"),
		verbose,
		xout,
		yout);
	}

	// compute the histogram noise model
	if (histexp) {
		if (verbose) {
			fprintf(stderr, "making pixel error model histogram\n");
		}
		beginhist();
		mkhist();
		endhist();
	}

	// make the output image
	imgout = Mat(yout, xout, CV_16UC3);

	// compute the superresolution image
	for (int i=0; i<imgcnt; ++i) fpruned[i] = 0;
	opruned = 0;
#pragma omp parallel for schedule(guided)
	for (int j=0; j<yout; ++j) {
		for (int i=0; i<xout; ++i) {
			imgout.at<VecTYP>(j, i) = getv(i, j);
		}
	}

	// release all input images because
	// the upcoming imwrite() needs memory space...
	if (verbose) {
		fprintf(stderr, "superresolution image built; freeing input images\n");
	}
	for (int i=0; i<imgcnt; ++i) {
		imgin[i].release();
	}

	if (verbose) {
		for (int i=0; i<imgcnt; ++i) {
			if (fpruned[0]) {
				fprintf(stderr,
"%0.2f%% of pixel values pruned from %s by first\n",
					(100.0 * fpruned[i]) / (3.0 * xout * yout),
					inname[i]);
			}
		}

		if (opruned) {
			fprintf(stderr,
"%0.2f%% of pixel values pruned as outliers\n",
				(100.0 * opruned) / (3.0 * xout * yout * imgcnt));
		}
	}

	if (verbose) {
		fprintf(stderr, "writing output file %s\n", outname);
	}
	fflush(stderr);
	fflush(stdout);

	imwrite(outname, imgout);

	return(0);
}