image.hpp

Gaussian Mixture Algorithm

/***************************************************************************
 *   Copyright (C) 2008 by Yann LeCun and Pierre Sermanet *
 *   [e-mail], [e-mail] *
 *
 * 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);