/*************************************************************************** *   Copyright (C) 2008 by Yann LeCun and Pierre Sermanet * *   yann@cs.nyu.edu, pierre.sermanet@gmail.com * * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: *     * Redistributions of source code must retain the above copyright *       notice, this list of conditions and the following disclaimer. *     * Redistributions in binary form must reproduce the above copyright *       notice, this list of conditions and the following disclaimer in the *       documentation and/or other materials provided with the distribution. *     * Redistribution under a license not approved by the Open Source *       Initiative (http://www.opensource.org) must display the *       following acknowledgement in all advertising material: *        This product includes software developed at the Courant *        Institute of Mathematical Sciences (http://cims.nyu.edu). *     * The names of the authors may not be used to endorse or promote products *       derived from this software without specific prior written permission. * * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * DISCLAIMED. IN NO EVENT SHALL ThE AUTHORS BE LIABLE FOR ANY * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ***************************************************************************/#ifndef IMAGE_HPP_#define IMAGE_HPP_#include <algorithm>#include <math.h>using namespace std;#define BLK_AVRG(nlin, ncol) {\  int k,l; int norm = ncol * nlin;\  int acc0=0, acc1=0, acc2=0;\  for (k=0; k<nlin; k++) {\  register T *pinptr = pin+k*in_mod0;\  for (l=0; l<ncol; l++) {\  acc0 += pinptr[0];acc1 += pinptr[1];acc2 += pinptr[2];\  pinptr += in_mod1; }}\  pout[0] = acc0 / norm;pout[1] = acc1 / norm;pout[2] = acc2 / norm;}namespace ebl {template<class T> Idx<T> image_crop(Idx<T> &in, int x, int y, int w, int h){	Idx<T> bla = in.narrow(0, h, y).narrow(1, w, x);	Idx<T> bla3(bla.dim(0), bla.dim(1), bla.order() < 3 ? -1 : bla.dim(2));	idx_copy(bla, bla3);	return bla3;}template<class T> Idx<T> image_resize(Idx<T> &image, double w, double h, int mode) {	Idx<T> im(image);	if (im.order() < 2) ylerror("image must have at least an order of 2.");	intg imw = im.dim(1);	intg imh = im.dim(0);	int rw = 0;	int rh = 0;	// determine actual size of output image	if ((0 == w) || (0 == h)) {		if (0 == w) {		  if (0 == h) {		  	ylerror("desired width and height cannot be both zero");		  } else {		  	w = max(1, (int) (imw * (h / imh)));		  }		} else h = max(1, (int) (imh * (w / imw)));	}	else if (mode == 0) {		double r = min(w / imw, h / imh);		w = max(1, (int) (r * imw));    h = max(1, (int) (r * imh));	}	else if (mode == 1) {		w = max(1, (int) w);		h = max(1, (int) h);	}	else if (mode == 2) {		w = max(1, (int) (w * imw));		h = max(1, (int) (h * imh));	}	else ylerror("image_resize: illegal mode or desired dimensions");	// compute closest integer subsampling ratio	rw = (int) imw / w;	rh = (int) imh / h;	// subsample by integer ratio if necessary	if ((0 != rh) || (0 != rw)) {		im = image_subsample(im, rh, rw);		imw = im.dim(1);		imh = im.dim(0);	}	// resample from subsampled image with bilinear interpolation	Idx<T> rez(h, w, (im.order() == 3) ? im.dim(2) : 1);	Idx<T> bg(4);	idx_clear(bg);	// the 0.5 thingies are necessary because warp-bilin interprets	// integer coordinates as beiing at the center of each pixel.	float x1 = -0.5, y1 = -0.5, x3 = imw - 0.5, y3 = imh - 0.5;	float p1 = -0.5, q1 = -0.5, p3 = w - 0.5, q3 = h - 0.5;	image_warp_quad(im, rez, bg, 1, x1, y1, x3, y1, x3, y3, x1, y3, p1, q1, p3, q3);	if (im.order() == 2)		return rez.select(2, 0);	return rez;}template<class T> Idx<T> image_subsample_grayscale(Idx<T> &in, int nlin, int ncol) {	intg h = in.dim(0);	intg w = in.dim(1);	intg nh = h / nlin;	intg nw = w / ncol;	Idx<T> out(nh, nw);	if ((nlin == 1) && (ncol == 1)) {		idx_copy(in, out);		return out;	}	Idx<T> inp = in.narrow(0, nlin * nh, 0);	inp = inp.narrow(1, ncol * nw, 0);	T *_idx2loopc1, *pin;	T *_idx2loopc2, *pout;	int i, _imax = out.dim(0);	int j, _jmax = out.dim(1);	int _imat1_m0 = inp.mod(0);	int _imat1_m1 = inp.mod(1);	int _imat2_m0 = out.mod(0);	int _imat2_m1 = out.mod(1);	int pin_incr = ncol * _imat1_m1;	int norm = ncol * nlin;	register int acc0;	register int k,l;	register T *pinptr;	register int pinptr_incr = _imat1_m0 - ncol* _imat1_m1;	_idx2loopc1 = inp.idx_ptr();	_idx2loopc2 = out.idx_ptr();	for (i = 0; i < _imax; i++) {		pin = _idx2loopc1;		pout = _idx2loopc2;		for (j = 0; j < _jmax; j++) {			acc0 =0;			pinptr = pin;			for (k=0; k<nlin; k++) {				for (l=0; l<ncol; l++) {					acc0 += pinptr[0];					pinptr += _imat1_m1;				}				pinptr += pinptr_incr;			}			pout[0] = acc0/norm;			pin += pin_incr;			pout += _imat2_m1;		}		_idx2loopc1 += _imat1_m0 * nlin;		_idx2loopc2 += _imat2_m0;	}	return out;}template<class T> Idx<T> image_subsample_rgb(Idx<T> &in, int nlin, int ncol) {	intg h = in.dim(0);	intg w = in.dim(1);	intg nh = h / nlin;	intg nw = w / ncol;	Idx<T> out(nh, nw, in.dim(2));	if ((nlin == 1) && (ncol == 1)) {		idx_copy(in, out);		return out;	}	Idx<T> inp = in.narrow(0, nlin * nh, 0).narrow(1, ncol * nw, 0);	T *in_line, *pin;	T *out_line, *pout;	int i, _imax = out.dim(0);	int j, _jmax = out.dim(1);	int in_mod0 = inp.mod(0);	int in_mod1 = inp.mod(1);	int out_mod0 = out.mod(0);	int out_mod1 = out.mod(1);	int pin_incr = ncol * in_mod1;	in_line = inp.idx_ptr();	out_line = out.idx_ptr();	for (i = 0; i < _imax; i++) {		pin = in_line;		pout = out_line;		for (j = 0; j < _jmax; j++) {			BLK_AVRG(nlin, ncol);			pin += pin_incr;			pout += out_mod1;		}		in_line += in_mod0 * nlin;		out_line += out_mod0;	}	return out;}template<class T> Idx<T> image_subsample(Idx<T> &in, int nlin, int ncol) {	switch (in.order()) {	case 2:		return image_subsample_grayscale(in, nlin, ncol);  case 3:		return image_subsample_rgb(in, nlin, ncol);  default:  	ylerror("image must have at least an order of 2.");  	return in;	}}template<class T> void image_warp_quad(Idx<T> &in, Idx<T> &out,		Idx<T> &background, int mode,		float x1, float y1, float x2, float y2, float x3, float y3,		float x4, float y4, float p1, float q1, float p3, float q3) {  intg outi = out.dim(0);  intg outj = out.dim(1);  Idx<int> dispi(outi, outj);  Idx<int> dispj(outi, outj);  compute_bilin_transform<float>(dispi, dispj, x1, y1, x2, y2, x3, y3,  		x4, y4, p1, q1, p3, q3);  if (0 == mode)	  image_warp_fast(in, out, background.idx_ptr(), dispi, dispj);  else	  image_warp(in, out, background, dispi, dispj);}template<class T> void image_warp(Idx<T> &in, Idx<T> &out, Idx<T> &background,		Idx<int> &pi, Idx<int> &pj) {  T* pin = in.idx_ptr();  int indimi = in.dim(0);  int indimj = in.dim(1);  int inmodi = in.mod(0);  int inmodj = in.mod(1);  int ppi, ppj;  { idx_bloop3(lout, out, T, lpi, pi, int, lpj, pj, int) {  	{ idx_bloop3(llout, lout, T, llpi, lpi, int, llpj, lpj, int) {  		ppi = llpi.get();  		ppj = llpj.get();      image_interpolate_bilin(background.idx_ptr(), pin, indimi, indimj,      		inmodi, inmodj, ppi, ppj, llout->idx_ptr(), (int) out.dim(2));  	}}  }}}template<class T> void image_warp_fast(Idx<T> &in, Idx<T> &out, T *background,		Idx<int> &pi, Idx<int> &pj) {	  T* pin = in.idx_ptr();	  int indimi = in.dim(0);	  int indimj = in.dim(1);	  int inmodi = in.mod(0);	  int inmodj = in.mod(1);	  int ppi, ppj;	  register int li, lj;	  register T *inn, *outt;	  int outsize = out.dim(2);	  { idx_bloop3(lout, out, T, lpi, pi, int, lpj, pj, int) {	  	{ idx_bloop3(llout, lout, T, llpi, lpi, int, llpj, lpj, int) {	  		outt = llout->idx_ptr();	  		ppi = llpi.get();	  		ppj = llpj.get();  		  li = (ppi+0x7f) >> 16;  		  lj = (ppj+0x7f) >> 16;  		  if ((li>=0) && (li<indimi) && (lj>=0) && (lj<indimj)) {  		  	inn = (T*)(pin) + inmodi * li + inmodj * lj;  		  	outt[0] = inn[0];  		  	if (outsize >= 3) {  		  		outt[1] = inn[1];  		  		outt[2] = inn[2];  		  	}  		  	if (outsize >= 4)  		  		outt[3] = inn[3];  		  } else {  		  	outt[0] = *(background);  		  	if (outsize >= 3) {  		  		outt[1] = *(background+1);  		  		outt[2] = *(background+2);