image.hpp

Gaussian Mixture Algorithm

  		  	}
  		  	if (outsize >= 4)
  		  		outt[3] = *(background+3);
  		  }
	  	}}
	  }}
}

template<class T> void image_interpolate_bilin(T* background, T *pin,
		int indimi, int indimj, int inmodi, int inmodj, int ppi, int ppj,
		T* out, int outsize) {
	int li0, lj0;
	register int li1, lj1;
	int deltai, ndeltai;
	int deltaj, ndeltaj;
	register T *pin00;
	register T *v00, *v01, *v10, *v11;
	li0 = ppi >> 16;
	li1 = li0+1;
	deltai = ppi & 0x0000ffff;
	ndeltai = 0x00010000 - deltai;
	lj0 = ppj  >> 16;
	lj1 = lj0+1;
	deltaj = ppj & 0x0000ffff;
	ndeltaj = 0x00010000 - deltaj;
	pin00 = (T*)(pin) + inmodi * li0 + inmodj * lj0;
	if ((li1>0)&&(li1<indimi)) {
		if ((lj1>0)&&(lj1<indimj)) {
			v00 = (pin00);
			v01 = (pin00+inmodj);
			v11 = (pin00+inmodi+inmodj);
			v10 = (pin00+inmodi);
		} else if (lj1==0) {
			v00 = background;
			v01 = (pin00+inmodj);
			v11 = (pin00+inmodi+inmodj);
			v10 = background;
		} else if (lj1==indimj) {
			v00 = (pin00);
			v01 = background;
			v11 = background;
			v10 = (pin00+inmodi);
		} else {
			v00 = background;
			v01 = background;
			v11 = background;
			v10 = background;
		}
	} else if (li1==0) {
		if ((lj1>0)&&(lj1<indimj)) {
			v00 = background;
			v01 = background;
			v11 = (pin00+inmodi+inmodj);
			v10 = (pin00+inmodi);
		} else if (lj1==0) {
			v00 = background;
			v01 = background;
			v11 = (pin00+inmodi+inmodj);
			v10 = background;
		} else if (lj1==indimj) {
			v00 = background;
			v01 = background;
			v11 = background;
			v10 = (pin00+inmodi);
		} else {
			v00 = background;
			v01 = background;
			v11 = background;
			v10 = background;
		}
	} else if (li1==indimi) {
		if ((lj1>0)&&(lj1<indimj)) {
			v00 = (pin00);
			v01 = (pin00+inmodj);
			v11 = background;
			v10 = background;
		} else if (lj1==0) {
			v00 = background;
			v01 = (pin00+inmodj);
			v11 = background;
			v10 = background;
		} else if (lj1==indimj) {
			v00 = (pin00);
			v01 = background;
			v11 = background;
			v10 = background;
		} else {
			v00 = background;
			v01 = background;
			v11 = background;
			v10 = background;
		}
	} else {
		v00 = background;
		v01 = background;
		v11 = background;
		v10 = background;
	}
	// TODO: this does not work for ubyte (cf ubyte implementation in image.cpp),
	// make it generic to avoid code redondancy.
	double dd = 1.0 / 65536.0;
	T d = (T) dd;
	if (sizeof (T) <= 2)
		d = 0;
	if (outsize == 1)
	  *out = (ndeltaj * (( *v10*deltai + *v00*ndeltai )*d) +
		  	   deltaj  * (( *v11*deltai + *v01*ndeltai )*d))*d;
	else {
		if (outsize >= 3) {
		   *out = (ndeltaj * (( v10[0]*deltai + v00[0]*ndeltai )*d) +
	               deltaj  * (( v11[0]*deltai + v01[0]*ndeltai )*d))*d;
		   *(out + 1) = (ndeltaj * (( v10[1]*deltai + v00[1]*ndeltai )*d) +
	               deltaj  * (( v11[1]*deltai + v01[1]*ndeltai )*d))*d;
		   *(out + 2) = (ndeltaj * (( v10[2]*deltai + v00[2]*ndeltai )*d) +
	               deltaj  * (( v11[2]*deltai + v01[2]*ndeltai )*d))*d;
		}
		if (outsize >= 4) {
		   *(out + 3) = (ndeltaj * (( v10[3]*deltai + v00[3]*ndeltai )*d) +
	               deltaj  * (( v11[3]*deltai + v01[3]*ndeltai )*d))*d;
		}
	}
}

template<class T> void compute_bilin_transform(Idx<int> &dispi, Idx<int> &dispj,
		float x1, float y1, float x2, float y2, float x3, float y3,
		float x4, float y4, float p1, float q1, float p3, float q3) {
	// compute transformation matrix from coordinates
	// in target (rectangular) space to coordinates
	// in original (irregular quadrilateral) image
	// transformation matrix is in 16.16 fixed point format.
	float k = 65536 / ((p3 - p1) * (q3 - q1));
	float x41 = x4 - x1;
	float x21 = x2 - x1;
	float x43 = x4 - x3;
	float x23 = x2 - x3;
	int mx0 = k * ((q3 * x21) + (q1 * x43));
	int mx1 = k * ((p3 * x41) + (p1 * x23));
	int mx2 = (-k) * (x41 + x23);
	int mx3 = k * (((p3 * q3 * x1) + (p1 * q1 * x3)) -
			((p1 * q3 * x2) + (p3 * q1 * x4)));
	float y41 = y4 - y1;
	float y21 = y2 - y1;
	float y43 = y4 - y3;
	float y23 = y2 - y3;
	int my0 = k * ((q3 * y21) + (q1 * y43));
	int my1 = k * ((p3 * y41) + (p1 * y23));
	int my2 = (-k) * (y41 + y23);
	int my3 = k * (((p3 * q3 * y1) + (p1 * q1 * y3)) -
			((p1 * q3 * y2) + (p3 * q1 * y4)));
	int q = 0, p = 0;
	{ idx_bloop2(ispi, dispi, int, ispj, dispj, int) {
		p = 0;
		{ idx_bloop2(di, ispi, int, dj, ispj, int) {
			di->set(my0 * p + my1 * q + my2 * p * q + my3);
			dj->set(mx0 * p + mx1 * q + mx2 * p * q + mx3);
			p++;
		}}
		q++;
	}}
}

template<class T> void image_rotscale(Idx<T> &src, Idx<T> &out,
		double sx, double sy, double dx, double dy,
		double angle, double coeff, Idx<ubyte> &bg){
	double q = 1000;
	double coeff_inv = 1/coeff;
	double sa = q*sin(angle * 0.017453292);
	double ca = q*cos(angle * 0.017453292);
	double ca_plus_sa = coeff_inv * (sa + ca);
	double ca_minus_sa = coeff_inv * (ca - sa);
	float x1 = sx - ca_plus_sa;
	float y1 = sy - ca_minus_sa;
	float x2 = sx + ca_minus_sa;
	float y2 = sy - ca_plus_sa;
	float x3 = sx + ca_plus_sa;
	float y3 = sy + ca_minus_sa;
	float x4 = sx - ca_minus_sa;
	float y4 = sy + ca_plus_sa;
	float p1 = dx-q;
	float q1 = dy-q;
	float p3 = dx + q;
	float q3 = dy + q;
	image_warp_quad(src, out, bg, 1, x1, y1, x2, y2, x3, y3, x4, y4, p1, q1, p3, q3);
}



template<class T> void image_draw_box(Idx<T> &img, T val,
		unsigned int x, unsigned int y, unsigned int dx, unsigned int dy) {
	idx_checkorder1(img, 2);
  for (unsigned int i = x; i < x + dx; ++i) {
  	img.set(val, i, y);
  	img.set(val, i, y + dy);
  }
  for (unsigned int j = y; j < y + dy; ++j) {
  	img.set(val, x, j);
  	img.set(val, x + dx, j);
  }
}

template<class T> bool pnm_fread_into_rgbx(const char *fname, Idx<T> &out) {
	Idx<ubyte> tmp(1,1,1);
	bool ret = pnm_fread_into_rgbx(fname, tmp);
	out.resize(tmp.dim(0), tmp.dim(1), tmp.dim(2));
	idx_copy(tmp, out);
	return ret;
}

template<class T> bool image_read_rgbx(const char *fname, Idx<T> &out) {
	Idx<ubyte> tmp(1,1,1);
	bool ret = image_read_rgbx(fname, tmp);
	out.resize(tmp.dim(0), tmp.dim(1), tmp.dim(2));
	idx_copy(tmp, out);
	return ret;
}

////////////////////////////////////////////////////////////////
// Utilities

template<class T> void RGBtoYUV1D(Idx<T> &rgb, Idx<T> &yuv) {
	if (rgb.idx_ptr() == yuv.idx_ptr()) {
		ylerror("RGBtoYUV: dst must be different than src");
		return ;
	}
	yuv.set(0.257 * rgb.get(0) + 0.504 * rgb.get(1) + 0.098 * rgb.get(2) + 16, 0);
	yuv.set(0.439 * rgb.get(0) - (0.368 * rgb.get(1)) - (0.071 * rgb.get(2)) + 128, 1);
	yuv.set(-(0.148 * rgb.get(0)) - (0.291 * rgb.get(1)) + 0.439 * rgb.get(2) + 128, 2);
}

template<class T> void RGBtoYUV(Idx<T> &rgb, Idx<T> &yuv) {
	idx_checknelems2_all(rgb, yuv);
	switch (rgb.order()) {
	case 1: // process 1 pixel
		RGBtoYUV1D(rgb, yuv);
		return ;
	case 3: // process 2D image
		{ idx_bloop2(rg, rgb, T, yu, yuv, T) {
			{ idx_bloop2(r, rg, T, y, yu, T) {
				RGBtoYUV1D(r, y);
			}}
		}}
		return ;
	default:
		ylerror("RGBtoYUV dimension not implemented");
	}
}

template<class T> void YUVtoRGB1D(Idx<T> &yuv, Idx<T> &rgb) {
	if (rgb.idx_ptr() == yuv.idx_ptr()) {
		ylerror("YUVtoRGB: dst must be different than src");
		return ;
	}
	rgb.set(1.164 * (yuv.get(0) - 16) + 2.018 * (yuv.get(1) - 128), 0);
	rgb.set(1.164 * (yuv.get(0) - 16) - 0.813 * (yuv.get(2) - 128)
			- 0.391 * (yuv.get(1) - 128), 1);
	rgb.set(1.164 * (yuv.get(0) - 16) + 1.596 * (yuv.get(2) - 128), 2);
}

template<class T> void YUVtoRGB(Idx<T> &yuv, Idx<T> &rgb) {
	idx_checknelems2_all(rgb, yuv);
	switch (yuv.order()) {
	case 1: // process 1 pixel
		YUVtoRGB1D(yuv, rgb);
		return ;
	case 3: // process 2D image
		{ idx_bloop2(rg, rgb, T, yu, yuv, T) {
			{ idx_bloop2(r, rg, T, y, yu, T) {
				YUVtoRGB1D(y, r);
			}}
		}}
		return ;
	default:
		ylerror("YUVtoRGB dimension not implemented");
	}
}

} // end namespace ebl

#endif /* IMAGE_HPP_ */