From 4dc49ca1178f28ff6a65a949733671c0a6779cc2 Mon Sep 17 00:00:00 2001
From: Ivan Mahonin <bh@icystar.com>
Date: Feb 14 2020 12:36:54 +0000
Subject: perspective: Perspective class


---

diff --git a/c++/perspective/src/SConstruct b/c++/perspective/src/SConstruct
index 216a321..22c9433 100644
--- a/c++/perspective/src/SConstruct
+++ b/c++/perspective/src/SConstruct
@@ -21,6 +21,7 @@ sources = [
 	'main.cpp',
 	'mainwindow.cpp',
 	'matrix.cpp',
+	'perspective.cpp',
 	'perspectivesandboxview.cpp',
 	'surface.cpp',
 	'view.cpp' ]
diff --git a/c++/perspective/src/common.h b/c++/perspective/src/common.h
index 029e22a..dedca83 100644
--- a/c++/perspective/src/common.h
+++ b/c++/perspective/src/common.h
@@ -51,6 +51,36 @@ public:
 	
 	inline explicit Color(Channel r = 0.0, Channel g = 0.0, Channel b = 0.0, Channel a = 0.0):
 		r(r), g(g), b(b), a(a) { }
+	
+	inline Color operator+ (const Color &other) const
+		{ return Color(r + other.r, g + other.g, b + other.b, a + other.a); }
+	inline Color operator- (const Color &other) const
+		{ return Color(r - other.r, g - other.g, b - other.b, a - other.a); }
+	inline Color operator* (Channel k) const
+		{ return Color(r*k, g*k, b*k, a*k); }
+	
+	inline Color& operator+= (const Color &other)
+		{ r += other.r, g += other.g, b += other.b, a += other.a; return *this; }
+	inline Color& operator-= (const Color &other)
+		{ r -= other.r, g -= other.g, b -= other.b, a -= other.a; return *this; }
+	inline Color& operator*= (Channel k)
+		{ r *= k, g *= k, b *= k, a *= k; return *this; }
+	
+	inline Color get_mult_alpha() const
+		{ return Color(r*a, g*a, b*a, a); }
+	inline Color get_demult_alpha() const  {
+		if (fabs(a) < real_precision) return *this;
+		Channel k = 1/a;
+		return Color(r*k, g*k, b*k, a);
+	}
+	
+	inline void mult_alpha()
+		{ r*=a, g*=a, b*=a; }
+	inline void demult_alpha() {
+		if (fabs(a) < real_precision) return;
+		Real k = 1/a;
+		r*=k, g*=k, b*=k;
+	}
 };
 
 
diff --git a/c++/perspective/src/generator.cpp b/c++/perspective/src/generator.cpp
index 67a34d8..d80baec 100644
--- a/c++/perspective/src/generator.cpp
+++ b/c++/perspective/src/generator.cpp
@@ -117,15 +117,15 @@ Generator::generate(const Vector2 &coord, const Vector2 &precision) const {
 }
 
 void
-Generator::generate(Surface &target, const Vector2 &lt, const Vector2 &rb) const {
+Generator::generate(Surface &target, const Pair2 &bounds) const {
 	if (target.empty()) return;
-	Vector2 d = rb - lt;
+	Vector2 d = bounds.distance();
 	d.x /= target.width();
 	d.y /= target.height();
 	Vector2 precision = d/2.0;
 	
-	Vector2 p = lt;
-	for(int r = 0; r < target.height(); ++r, p.x = lt.x, p.y += d.y) {
+	Vector2 p = bounds.p0;
+	for(int r = 0; r < target.height(); ++r, p.x = bounds.p0.x, p.y += d.y) {
 		Color *row = target[r];
 		for(int c = 0; c < target.width(); ++c, p.x += d.x) {
 			Vector4 v = generate(p, precision);
diff --git a/c++/perspective/src/generator.h b/c++/perspective/src/generator.h
index e53476a..3d44432 100644
--- a/c++/perspective/src/generator.h
+++ b/c++/perspective/src/generator.h
@@ -68,7 +68,7 @@ public:
 	Real random(const ULongIntVector2 &coord, int index) const;
 	Vector4 generate(const ULongIntVector2 &coord, bool shrink_cache = true) const;
 	Vector4 generate(const Vector2 &coord, const Vector2 &precision) const;
-	void generate(Surface &target, const Vector2 &lt, const Vector2 &rb) const;
+	void generate(Surface &target, const Pair2 &bounds) const;
 };
 
 
diff --git a/c++/perspective/src/mainwindow.cpp b/c++/perspective/src/mainwindow.cpp
index 0935197..3057964 100644
--- a/c++/perspective/src/mainwindow.cpp
+++ b/c++/perspective/src/mainwindow.cpp
@@ -8,15 +8,262 @@
 #include "mainwindow.h"
 
 
+// it's not a math cross product
+// just speculative name for function: dot(a, Vector2(-y, x))
+static Real
+cross_product2(const Vector2 &a, const Vector2 &b)
+	{ return a*b.perp(); }
+
+static bool
+line_cross(const Pair2 &line_a, const Pair2 &line_b, Vector2 &out) {
+	Vector2 da = line_a.distance();
+	Vector2 db = line_b.distance();
+	Real denominator = cross_product2(da, db);
+	if (real_equal(denominator, 0.0))
+		return false;
+	Real numerator = cross_product2(da, line_a.p0 - line_b.p0);
+	out = line_a.p0 + db*(numerator/denominator);
+	return true;
+}
+
+static Vector4
+make_matrix_row(const Vector2 &p0, const Vector2 &p1)
+	{ return Vector4(p1.x - p0.x, p1.y - p0.y, 0.0, 0.0); }
+
+static Vector4
+make_matrix_row(const Vector2 &p0, const Vector2 &p1, const Vector2 &p2) {
+	Real l2 = (p2 - p1).length();
+	if (real_less_or_equal(l2, 0.0))
+		return Vector4();
+	Real l1 = (p1 - p0).length();
+	Real w = l1/l2;
+	return Vector4(p2.x*w, p2.y*w, 0.0, w);
+}
+
+static Matrix4
+make_perspective_matrix(const Vector2 &p0, const Vector2 &px, const Vector2 &py, const Vector2 &p1) {
+	Vector2 focus;
+	Matrix4 m;
+	m[0] = line_cross(Pair2(p0, px), Pair2(py, p1), focus)
+		 ? make_matrix_row(p0, px, focus)
+		 : make_matrix_row(p0, px);
+	m[1] = line_cross(Pair2(p0, py), Pair2(px, p1), focus)
+		 ? make_matrix_row(p0, py, focus)
+		 : make_matrix_row(p0, py);
+	m[2] = Vector4(0.0, 0.0, 1.0, 0.0);
+	m[3] = Vector4(p0.x, p0.y, 0.0, 1.0);
+	return m;
+}
+
+static Matrix4
+make_plain_matrix(const Pair2 &bounds) {
+	const Vector2 &p = bounds.p0;
+	const Vector2 d = bounds.distance();
+	return Matrix4(
+		Vector4(d.x, 0.0, 0.0, 0.0),
+		Vector4(0.0, d.y, 0.0, 0.0),
+		Vector4(0.0, 0.0, 1.0, 0.0),
+		Vector4(p.x, p.y, 0.0, 1.0) );
+}
+
+static void
+make_layers(const Pair2 &bounds, const Matrix4 &matrix, MainWindow::LayerList &out_layers) {
+	const Real min_distance_to_horizon = 1.0; // one pixel
+	
+	if (bounds.empty() || real_equal(matrix.row_w().w, 0.0))
+		return;
+	
+	// get focuses
+	Vector2 focus_x, focus_y;
+	bool focus_x_exists = real_not_equal(matrix.row_x().w, 0.0);
+	bool focus_y_exists = real_not_equal(matrix.row_y().w, 0.0);
+	if (focus_x_exists)
+		focus_x = matrix.row_x().vec2() / matrix.row_x().w;
+	if (focus_y_exists)
+		focus_y = matrix.row_y().vec2() / matrix.row_y().w;
+
+	// get horizon
+	Vector2 h, dh;
+	if (focus_x_exists && focus_x_exists) {
+		h = focus_x;
+		dh = (focus_y - focus_x).normalized();
+	} else
+	if (focus_x_exists) {
+		h = focus_x;
+		dh = matrix.row_y().vec2().normalized();
+	} else
+	if (focus_y_exists) {
+		h = focus_y;
+		dh = matrix.row_x().vec2().normalized();
+	} else {
+		// TODO
+		return;
+	}
+	
+	// get horison perpendicular
+	Vector2 hp = dh.perp().normalized();
+	
+	// get perspective bounds
+	Vector2 corners[] = {
+		Vector2(bounds.p0.x, bounds.p0.y),
+		Vector2(bounds.p0.x, bounds.p1.y),
+		Vector2(bounds.p1.x, bounds.p1.y),
+		Vector2(bounds.p1.x, bounds.p0.y),
+		Vector2(bounds.p0.x, bounds.p0.y) }; // close the ring
+	Real lmin = 1.0/real_precision;
+	Real lmax = 0.0;
+	for(int i = 0; i < 4; ++i) {
+		Real p = hp*(corners[i] - h);
+		if (p < lmin) lmin = p;
+		if (p > lmax) lmax = p;
+	}
+	lmin = std::max(min_distance_to_horizon, lmin);
+	lmax = std::max(min_distance_to_horizon, lmax);
+	if (real_less_or_equal(lmax, min_distance_to_horizon))
+		return;
+	
+	// layers
+	Real scale_step = 4.0;
+	Real l = lmin;
+	Matrix4 back_matrix = matrix.inverted();
+	while(real_less(l, lmax)) {
+		Vector2 p0 = h + hp*l;
+		l *= scale_step;
+		Vector2 p1 = h + hp*l;
+		
+		Vector4 points[8];
+		int points_count = 0;
+		for(int i = 0; i < 4; ++i) {
+			const Vector2 &cp0 = corners[i];
+			const Vector2 &cp1 = corners[i+1];
+			Vector2 dc = cp1 - cp0;
+			Real kk = hp*dc;
+			if (real_not_equal(kk, 0.0)) {
+				kk = 1.0/kk;
+				Real cl0 = clamp( hp*(p0 - cp0)*kk, 0.0, 1.0 );
+				Real cl1 = clamp( hp*(p1 - cp0)*kk, 0.0, 1.0 );
+				if (real_not_equal(cl1 - cl0, 0.0)) {
+					points[points_count++] = Vector4(cp0 + dc*cl0, 0.0, 1.0);
+					points[points_count++] = Vector4(cp0 + dc*cl1, 0.0, 1.0);
+				}
+			}
+		}
+
+		Vector2 src_res;
+		{
+			Vector2 c = (p0 + p1)*0.5;
+			Vector2 px = c*0.75 + p1*0.25;
+			Vector2 dx = px - c;
+			Vector2 dy = dx.perp();
+			Vector2 py = c + dy;
+			Real len = dx.length();
+			
+			Vector2 cc  = (back_matrix*Vector4(c , 0, 1)).persp_divide().vec2();
+			Vector2 cpx = (back_matrix*Vector4(px, 0, 1)).persp_divide().vec2();
+			Vector2 cpy = (back_matrix*Vector4(py, 0, 1)).persp_divide().vec2();
+			Vector2 ddx = cpx - cc;
+			Vector2 ddy = cpy - cc;
+			src_res = Vector2(
+				len/std::max(fabs(ddx.x), fabs(ddy.x)),
+				len/std::max(fabs(ddx.y), fabs(ddy.y)) );
+		}
+		
+
+		Pair2 src_rect;
+		Pair2 dst_rect;
+		for(int i = 0; i < points_count; ++i) {
+			const Vector4 &p = points[i];
+			Vector4 pp = back_matrix*p;
+			assert(real_not_equal(pp.w, 0.0));
+			pp /= pp.w;
+			if (i) {
+				if (src_rect.p0.x > pp.x) src_rect.p0.x = pp.x;
+				if (src_rect.p0.y > pp.y) src_rect.p0.y = pp.y;
+				if (src_rect.p1.x < pp.x) src_rect.p1.x = pp.x;
+				if (src_rect.p1.y < pp.y) src_rect.p1.y = pp.y;
+
+				if (dst_rect.p0.x > p.x) dst_rect.p0.x = p.x;
+				if (dst_rect.p0.y > p.y) dst_rect.p0.y = p.y;
+				if (dst_rect.p1.x < p.x) dst_rect.p1.x = p.x;
+				if (dst_rect.p1.y < p.y) dst_rect.p1.y = p.y;
+			} else {
+				src_rect.p0 = src_rect.p1 = pp.vec2();
+				dst_rect.p0 = dst_rect.p1 = points[i].vec2();
+			}
+		}
+		
+		if (!src_rect.empty() && !dst_rect.empty())
+		{
+			Vector2 src_size_r = src_rect.size();
+			IntVector2 src_size(
+				(int)ceil(src_size_r.x * src_res.x),
+				(int)ceil(src_size_r.y * src_res.y) );
+			if (src_size.x > 0 && src_size.y > 0 && src_size.x*src_size.x < 100000000) {
+				out_layers.push_back(MainWindow::Layer());
+				MainWindow::Layer &layer = out_layers.back();
+				layer.src_bounds = src_rect;
+				layer.src_size = src_size;
+				layer.dst_bounds.p0.x = (int)floor(dst_rect.p0.x) - 2;
+				layer.dst_bounds.p0.y = (int)floor(dst_rect.p0.y) - 2;
+				layer.dst_bounds.p1.x = (int)ceil(dst_rect.p1.x) + 2;
+				layer.dst_bounds.p1.y = (int)ceil(dst_rect.p1.y) + 2;
+				layer.back_transform = back_matrix;
+			}
+		}
+	}
+}
+
+static void
+make_surface(const Generator &generator, MainWindow::Layer &layer) {
+	if (!layer.src_surface) {
+		layer.dst_surface.reset();
+		layer.cairo_surface.clear();
+		
+		if (layer.src_size.x > 0 && layer.src_size.y > 0) {
+			Surface *surface = new DataSurface(layer.src_size.x, layer.src_size.y);
+			generator.generate(*surface, layer.src_bounds);
+			layer.src_surface = RefPtr<Surface>(surface);
+		}
+	}
+	
+	if (layer.src_surface && !layer.dst_surface) {
+		layer.cairo_surface.clear();
+		
+		if ( layer.dst_bounds.p0.x < layer.dst_bounds.p1.x
+		  && layer.dst_bounds.p0.y < layer.dst_bounds.p1.y )
+		{
+			int w = layer.dst_bounds.p1.x - layer.dst_bounds.p0.x;
+			int h = layer.dst_bounds.p1.y - layer.dst_bounds.p0.y;
+			Surface *surface = new DataSurface(w, h);
+			for(int y = 0; y < h; ++y) {
+				for(int x = 0; x < w; ++x) {
+					Vector4 pos = layer.back_transform * Vector4(x + layer.dst_bounds.p0.x, y + layer.dst_bounds.p0.y, 0, 1);
+					if (fabs(pos.w) > real_precision)
+						pos /= pos.w;
+					surface->row(y)[x] = layer.src_surface->get_pixel(pos.x, pos.y);
+				}
+			}
+			layer.dst_surface = RefPtr<Surface>(surface);
+		}
+	}
+	
+	if (layer.dst_surface && !layer.cairo_surface) {
+		layer.cairo_surface = layer.dst_surface->to_cairo_surface();
+	}
+}
+
+
+
+
 MainWindow::MainWindow():
-	src_rect0   (new View::Point( Vector2(-1, -1) )),
-	src_rect1   (new View::Point( Vector2( 1,  1) )),
-	dst_rect0   (new View::Point( Vector2(-1, -1) )),
-	dst_rect1   (new View::Point( Vector2( 1, -1) )),
-	dst_rect2   (new View::Point( Vector2( 1,  1) )),
-	dst_rect3   (new View::Point( Vector2(-1,  1) )),
-	dst_bounds0 (new View::Point( Vector2(-2, -2) )),
-	dst_bounds1 (new View::Point( Vector2( 2,  2) ))
+	src_rect_p0   (new View::Point( Vector2(-1, -1) )),
+	src_rect_p1   (new View::Point( Vector2( 1,  1) )),
+	dst_rect_p0   (new View::Point( Vector2(-1, -1) )),
+	dst_rect_px   (new View::Point( Vector2( 1, -1) )),
+	dst_rect_py   (new View::Point( Vector2(-1,  1) )),
+	dst_rect_p1   (new View::Point( Vector2( 1,  1) )),
+	dst_bounds_p0 (new View::Point( Vector2(-2, -2) )),
+	dst_bounds_p1 (new View::Point( Vector2( 2,  2) ))
 {
 	std::cout << "generator seed: " << generator.seed() << std::endl;
 	set_default_size(750, 500);
@@ -25,8 +272,8 @@ MainWindow::MainWindow():
 	hbox->set_homogeneous(true);
 	
 	src_view.transform.scale(50, 50);
-	src_view.points.push_back(src_rect0);
-	src_view.points.push_back(src_rect1);
+	src_view.points.push_back(src_rect_p0);
+	src_view.points.push_back(src_rect_p1);
 	src_view.signal_point_motion.connect(
 		sigc::mem_fun(*this, &MainWindow::on_point_motion) );
 	src_view.signal_point_changed.connect(
@@ -40,12 +287,12 @@ MainWindow::MainWindow():
 	hbox->add(src_view);
 	
 	dst_view.transform.scale(50, 50);
-	dst_view.points.push_back(dst_rect0);
-	dst_view.points.push_back(dst_rect1);
-	dst_view.points.push_back(dst_rect2);
-	dst_view.points.push_back(dst_rect3);
-	dst_view.points.push_back(dst_bounds0);
-	dst_view.points.push_back(dst_bounds1);
+	dst_view.points.push_back(dst_rect_p0);
+	dst_view.points.push_back(dst_rect_px);
+	dst_view.points.push_back(dst_rect_py);
+	dst_view.points.push_back(dst_rect_p1);
+	dst_view.points.push_back(dst_bounds_p0);
+	dst_view.points.push_back(dst_bounds_p1);
 	dst_view.signal_point_motion.connect(
 		sigc::mem_fun(*this, &MainWindow::on_point_motion) );
 	dst_view.signal_point_changed.connect(
@@ -64,13 +311,14 @@ MainWindow::MainWindow():
 	show_all();
 	
 	update_view_surface();
+	update_layers();
 }
 
 Gtk::Widget*
 MainWindow::generator_demo() {
 	std::cout << "generator demo: generate" << std::endl;
 	DataSurface surface(256, 256);
-	generator.generate(surface, Vector2(0, 0), Vector2(16, 16));
+	generator.generate(surface, Pair2(Vector2(0, 0), Vector2(16, 16)));
 	
 	std::cout << "generator demo: convert" << std::endl;
 	Glib::RefPtr<Gdk::Pixbuf> pixbuf = Gdk::Pixbuf::create(
@@ -96,8 +344,9 @@ MainWindow::update_view_surface() {
 	int w = src_view.get_width();
 	int h = src_view.get_height();
 	
-	Vector2 rect0(transform * Vector4(0, 0, 0, 1));
-	Vector2 rect1(transform * Vector4(w, h, 0, 1));
+	Pair2 rect(
+		(transform * Vector4(0, 0, 0, 1)).vec2(),
+		(transform * Vector4(w, h, 0, 1)).vec2() );
 	
 	if (w <= 0 || h <= 0) {
 		view_surface.clear();
@@ -105,14 +354,59 @@ MainWindow::update_view_surface() {
 	if ( !view_surface
 	  || view_surface->get_width() != w
 	  || view_surface->get_height() != h
-	  || view_surface_rect0 != rect0
-	  || view_surface_rect1 != rect1 )
+	  || view_surface_rect != rect )
 	{
 		DataSurface surface(w, h);
-		generator.generate(surface, rect0, rect1);
+		generator.generate(surface, rect);
 		view_surface = surface.to_cairo_surface();
-		view_surface_rect0 = rect0;
-		view_surface_rect1 = rect1;
+		view_surface_rect = rect;
+	}
+}
+
+void
+MainWindow::update_layers() {
+	Matrix dst_transform = dst_view.transform_to_pixels();
+	Vector2 dst_rect_p0_pixels(dst_transform * Vector4(dst_rect_p0->position, 0.0, 1.0));
+	Vector2 dst_rect_px_pixels(dst_transform * Vector4(dst_rect_px->position, 0.0, 1.0));
+	Vector2 dst_rect_py_pixels(dst_transform * Vector4(dst_rect_py->position, 0.0, 1.0));
+	Vector2 dst_rect_p1_pixels(dst_transform * Vector4(dst_rect_p1->position, 0.0, 1.0));
+	Vector2 dst_bounds_p0_pixels(dst_transform * Vector4(dst_bounds_p0->position, 0.0, 1.0));
+	Vector2 dst_bounds_p1_pixels(dst_transform * Vector4(dst_bounds_p1->position, 0.0, 1.0));
+
+	Matrix in_matrix = make_plain_matrix( Pair2(
+		src_rect_p0->position,
+		src_rect_p1->position ));
+	Matrix out_matrix = make_perspective_matrix(
+		dst_rect_p0_pixels,
+		dst_rect_px_pixels,
+		dst_rect_py_pixels,
+		dst_rect_p1_pixels );
+	Matrix matrix = out_matrix * in_matrix.inverted();
+
+	Pair2 bounds(
+		dst_bounds_p0_pixels,
+		dst_bounds_p1_pixels );
+	bounds.sort();
+
+	LayerList layers;
+	make_layers(bounds, matrix, layers);
+	
+	for(LayerList::iterator i = layers.begin(); i != layers.end(); ++i)
+		make_surface(generator, *i);
+	
+	target_surface.clear();
+	int w = dst_view.get_width();
+	int h = dst_view.get_height();
+	if (w > 0 && h > 0) {
+		target_surface = Cairo::ImageSurface::create(Cairo::FORMAT_ARGB32, w, h);
+		Cairo::RefPtr<Cairo::Context> context = Cairo::Context::create(target_surface);
+		context->set_operator(Cairo::OPERATOR_ADD);
+		for(LayerList::const_iterator i = layers.begin(); i != layers.end(); ++i) {
+			if (i->cairo_surface) {
+				context->set_source(i->cairo_surface, i->dst_bounds.p0.x, i->dst_bounds.p0.y);
+				context->paint();
+			}
+		}
 	}
 }
 
@@ -123,6 +417,7 @@ MainWindow::on_point_motion(const View::PointPtr &/*point*/) {
 
 void
 MainWindow::on_point_changed(const View::PointPtr &/*point*/) {
+	update_layers();
 }
 
 void
@@ -145,16 +440,16 @@ MainWindow::on_draw_src_view(const Cairo::RefPtr<Cairo::Context> &context) {
 		update_view_surface();
 	if (view_surface) {
 		context->save();
-		src_view.back_transform_context(context);
+		context->transform( src_view.transform_from_pixels().to_cairo() );
 		context->set_source(view_surface, 0, 0);
 		context->paint();
 		context->restore();
 	}
 	
-	context->move_to(src_rect0->position.x, src_rect0->position.y);
-	context->line_to(src_rect0->position.x, src_rect1->position.y);
-	context->line_to(src_rect1->position.x, src_rect1->position.y);
-	context->line_to(src_rect1->position.x, src_rect0->position.y);
+	context->move_to(src_rect_p0->position.x, src_rect_p0->position.y);
+	context->line_to(src_rect_p0->position.x, src_rect_p1->position.y);
+	context->line_to(src_rect_p1->position.x, src_rect_p1->position.y);
+	context->line_to(src_rect_p1->position.x, src_rect_p0->position.y);
 	context->close_path();
 	context->set_line_width(ps);
 	context->set_source_rgba(1, 1, 0, 0.5);
@@ -168,20 +463,27 @@ MainWindow::on_draw_dst_view(const Cairo::RefPtr<Cairo::Context> &context) {
 	const Real ps = src_view.get_pixel_size();
 
 	context->save();
+	if (target_surface) {
+		context->save();
+		context->transform( dst_view.transform.inverted().to_cairo() );
+		context->set_source(target_surface, 0, 0);
+		context->paint();
+		context->restore();
+	}
 	
-	context->move_to(dst_rect0->position.x, dst_rect0->position.y);
-	context->line_to(dst_rect1->position.x, dst_rect1->position.y);
-	context->line_to(dst_rect2->position.x, dst_rect2->position.y);
-	context->line_to(dst_rect3->position.x, dst_rect3->position.y);
+	context->move_to(dst_rect_p0->position.x, dst_rect_p0->position.y);
+	context->line_to(dst_rect_px->position.x, dst_rect_px->position.y);
+	context->line_to(dst_rect_p1->position.x, dst_rect_p1->position.y);
+	context->line_to(dst_rect_py->position.x, dst_rect_py->position.y);
 	context->close_path();
 	context->set_line_width(ps);
 	context->set_source_rgba(1, 1, 0, 0.5);
 	context->stroke();
 
-	context->move_to(dst_bounds0->position.x, dst_bounds0->position.y);
-	context->line_to(dst_bounds0->position.x, dst_bounds1->position.y);
-	context->line_to(dst_bounds1->position.x, dst_bounds1->position.y);
-	context->line_to(dst_bounds1->position.x, dst_bounds0->position.y);
+	context->move_to(dst_bounds_p0->position.x, dst_bounds_p0->position.y);
+	context->line_to(dst_bounds_p0->position.x, dst_bounds_p1->position.y);
+	context->line_to(dst_bounds_p1->position.x, dst_bounds_p1->position.y);
+	context->line_to(dst_bounds_p1->position.x, dst_bounds_p0->position.y);
 	context->close_path();
 	context->set_line_width(ps);
 	context->set_source_rgba(1, 0, 0, 0.5);
diff --git a/c++/perspective/src/mainwindow.h b/c++/perspective/src/mainwindow.h
index 16ee403..8e8ac63 100644
--- a/c++/perspective/src/mainwindow.h
+++ b/c++/perspective/src/mainwindow.h
@@ -2,6 +2,8 @@
 #define MAINWINDOW_H
 
 
+#include <vector>
+
 #include <gtkmm/window.h>
 
 #include "view.h"
@@ -11,24 +13,42 @@
 
 class MainWindow: public Gtk::Window {
 public:
+	class Layer {
+	public:
+		IntVector2 src_size;
+		Pair2 src_bounds;
+		RefPtr<Surface> src_surface;
+
+		IntPair2 dst_bounds;
+		Matrix back_transform;
+		RefPtr<Surface> dst_surface;
+		Cairo::RefPtr<Cairo::ImageSurface> cairo_surface;
+
+		void process(const RefPtr<Surface> &dst_surface) const;
+	};
+
+	typedef std::vector<Layer> LayerList;
+
 	Generator generator;
 
 	View src_view;
 	View dst_view;
 	
 	View::PointPtr
-		src_rect0, src_rect1,
-		dst_rect0, dst_rect1, dst_rect2, dst_rect3,
-		dst_bounds0, dst_bounds1;
+		src_rect_p0, src_rect_p1,
+		dst_rect_p0, dst_rect_px, dst_rect_py, dst_rect_p1,
+		dst_bounds_p0, dst_bounds_p1;
 	
-	Vector2 view_surface_rect0;
-	Vector2 view_surface_rect1;
+	Pair2 view_surface_rect;
 	Cairo::RefPtr<Cairo::ImageSurface> view_surface;
-	
+
+	Cairo::RefPtr<Cairo::ImageSurface> target_surface;
+
 	MainWindow();
 	
 	Gtk::Widget* generator_demo();
 	
+	void update_layers();
 	void update_view_surface();
 	
 	void on_point_motion(const View::PointPtr &point);
diff --git a/c++/perspective/src/matrix.cpp b/c++/perspective/src/matrix.cpp
index 048f9ad..29cbd3c 100644
--- a/c++/perspective/src/matrix.cpp
+++ b/c++/perspective/src/matrix.cpp
@@ -2,6 +2,42 @@
 
 
 Real
+Matrix3::det() const {
+    return m00*(m11*m22 - m12*m21)
+         + m01*(m12*m20 - m10*m22)
+         + m02*(m10*m21 - m11*m20);
+}
+
+Matrix3
+Matrix3::inverted(bool *success) const {
+	const Real rm00 = m11*m22 - m12*m21;
+	const Real rm10 = m12*m20 - m10*m22;
+	const Real rm20 = m10*m21 - m11*m20;
+	const Real d = m00*rm00 + m01*rm10 + m02*rm20;
+	
+	Real k;
+	if (fabs(d) > real_precision) {
+		k = 1/d;
+		if (success) *success = true;
+	} else {
+		k = 0;
+		if (success) *success = false;
+	}
+	
+	return Matrix3(
+		Vector3( m11*m22 - m12*m21,
+				 m02*m21 - m01*m22,
+				 m01*m12 - m02*m11 )*k,
+		Vector3( m12*m20 - m10*m22,
+				 m00*m22 - m02*m20,
+				 m02*m10 - m00*m12 )*k,
+		Vector3( m10*m21 - m11*m20,
+				 m01*m20 - m00*m21,
+				 m00*m11 - m01*m10 )*k );
+}
+
+
+Real
 Matrix4::det() const {
     return m00 * (m11*(m22*m33 - m23*m32) + m12*(m23*m31 - m21*m33) + m13*(m21*m32 - m22*m31))
          + m01 * (m01*(m23*m32 - m22*m33) + m02*(m21*m33 - m23*m31) + m03*(m22*m31 - m21*m32))
diff --git a/c++/perspective/src/matrix.h b/c++/perspective/src/matrix.h
index 9c49992..a7d962c 100644
--- a/c++/perspective/src/matrix.h
+++ b/c++/perspective/src/matrix.h
@@ -7,6 +7,89 @@
 #include "vector.h"
 
 
+
+class Matrix3 {
+public:
+	union {
+		struct {
+			Real m00, m01, m02,
+				 m10, m11, m12,
+				 m20, m21, m22;
+		};
+		struct { Real m[3][3]; };
+		struct { Real a[9]; };
+	};
+	
+	inline explicit Matrix3(
+		const Vector3 &x = Vector3(1, 0, 0),
+		const Vector3 &y = Vector3(0, 1, 0),
+		const Vector3 &z = Vector3(0, 0, 1)
+	):
+		m00(x.x), m01(x.y), m02(x.z),
+		m10(y.x), m11(y.y), m12(y.z),
+		m20(z.x), m21(z.y), m22(z.z) { }
+	
+	inline Vector3& operator[] (int index) { return Vector3::cast(m[index]); }
+	inline const Vector3& operator[] (int index) const { return Vector3::cast(m[index]); }
+	
+	inline Vector3& row_x() { return (*this)[0]; }
+	inline Vector3& row_y() { return (*this)[1]; }
+	inline Vector3& row_z() { return (*this)[2]; }
+
+	inline const Vector3& row_x() const { return (*this)[0]; }
+	inline const Vector3& row_y() const { return (*this)[1]; }
+	inline const Vector3& row_z() const { return (*this)[2]; }
+	
+	inline const Vector3 get_col(int index) const
+		{ return Vector3( m[0][index], m[1][index], m[2][index] ); }
+	
+	inline Vector3 operator* (const Vector3 &v) const
+		{ return row_x()*v.x + row_y()*v.y + row_z()*v.z; }
+	inline Matrix3 operator* (const Matrix3 &other) const {
+		return Matrix3( *this * other.row_x(),
+						*this * other.row_y(),
+						*this * other.row_z() );
+	}
+	
+	inline Matrix3& operator*= (const Matrix3 &other)
+		{ return *this = *this * other; }
+	
+	Real det() const;
+	Matrix3 inverted(bool *success = NULL) const;
+	inline Matrix3& invert(bool *success = NULL)
+		{ return *this = inverted(); }
+	
+	static Real det(const Matrix3 matrix, const Matrix3 preinverted_matrix);
+	static inline Matrix3 zero() { return Matrix3(Vector3(), Vector3(), Vector3()); }
+
+	static inline Matrix3 identity() { return Matrix3(); }
+	static inline Matrix3 translation(const Vector2 &v) {
+		return Matrix3(
+			Vector3(1, 0, 0),
+			Vector3(0, 1, 0),
+			Vector3(v   , 1) );
+	}
+	static inline Matrix3 scaling(const Vector2 &v) {
+		return Matrix3(
+			Vector3(v.x,   0, 0),
+			Vector3(  0, v.y, 0),
+			Vector3(  0,   0, 1) );
+	}
+	
+	std::string to_string() const {
+		return "(" + row_x().to_string() + ", "
+			       + row_y().to_string() + ", "
+			       + row_z().to_string() + ")";
+	}
+	
+	Cairo::Matrix to_cairo() const {
+		return Cairo::Matrix(
+			m00, m10, m01,
+			m11, m20, m21 );
+	}
+};
+
+
 class Matrix4 {
 public:
 	union {
diff --git a/c++/perspective/src/perspective.cpp b/c++/perspective/src/perspective.cpp
new file mode 100644
index 0000000..7af2f54
--- /dev/null
+++ b/c++/perspective/src/perspective.cpp
@@ -0,0 +1,374 @@
+
+#include "perspective.h"
+
+
+namespace {
+	const int border_width  = 2;
+	const int max_width     = 16384 - 2*border_width;
+	const int max_height    = 16384 - 2*border_width;
+
+	const int max_area_width = 4096 - 2*border_width;
+	const int max_area       = max_area_width * max_area_width;
+
+	const Real max_scale_sqr = 4.0 * 4.0;
+}
+
+
+
+int
+Perspective::truncate_line(
+	Vector2 *out_points,
+	const Pair2 &bounds,
+	Real a,
+	Real b,
+	Real c )
+{
+	// equatation of line is a*x + b*y + c = 0
+	
+	int count = 0;
+	
+	if (fabs(a) > real_precision) {
+		Real x0 = -(c + bounds.p0.y*b)/a;
+		if ( x0 + real_precision >= bounds.p0.x
+		  && x0 - real_precision <= bounds.p1.x )
+		{
+			if (out_points) out_points[count] = Vector2(x0, bounds.p0.y);
+			if (++count >= 2) return count;
+		}
+		
+		Real x1 = -(c + bounds.p1.y*b)/a;
+		if ( x1 + real_precision >= bounds.p0.x
+		  && x1 - real_precision <= bounds.p1.x )
+		{
+			if (out_points) out_points[count] = Vector2(x1, bounds.p1.y);
+			if (++count >= 2) return count;
+		}
+	}
+
+	if (fabs(b) > real_precision) {
+		Real y0 = -(c + bounds.p0.x*a)/b;
+		if ( y0 + real_precision >= bounds.p0.y
+		  && y0 - real_precision <= bounds.p1.y )
+		{
+			if (out_points) out_points[count] = Vector2(bounds.p0.x, y0);
+			if (++count >= 2) return count;
+		}
+		
+		Real y1 = -(c + bounds.p1.x*a)/b;
+		if ( y1 + real_precision >= bounds.p0.y
+		  && y1 - real_precision <= bounds.p1.y )
+		{
+			if (out_points) out_points[count] = Vector2(bounds.p1.x, y1);
+			if (++count >= 2) return count;
+		}
+	}
+
+	return count;
+}
+
+
+Vector2
+Perspective::calc_optimal_resolution(
+	const Vector2 &ox,
+	const Vector2 &oy )
+{
+	const Real a = ox.x * ox.x;
+	const Real b = ox.y * ox.y;
+	const Real c = oy.x * oy.x;
+	const Real d = oy.y * oy.y;
+	Real e = fabs(ox.x*oy.y - ox.y*oy.x);
+	if (e < real_precision)
+		return Vector2(); // paranoid check, e must be non-zero when matrix is invertible
+	e = 1.0/e;
+	
+	const Real sum = a*d + b*c;
+	Vector2 scale;
+	if (2*a*b + real_precision >= sum) {
+		scale.x = sqrt(2*b)*e;
+		scale.y = sqrt(2*a)*e;
+	} else
+	if (2*c*d + real_precision >= sum) {
+		scale.x = sqrt(2*d)*e;
+		scale.y = sqrt(2*c)*e;
+	} else {
+		const Real dif = a*d - b*c;
+		scale.x = sqrt(dif/(a - c))*e;
+		scale.y = sqrt(dif/(d - b))*e;
+	}
+
+	return scale.x <= real_precision || scale.y <= real_precision
+	     ? Vector2() : scale;
+}
+
+
+Vector3
+Perspective::make_alpha_matrix_col(
+	Real w0,
+	Real w1,
+	const Vector3 &w_col )
+{
+	Real k = w1 - w0;
+	if (fabs(k) <= real_precision)
+		return w_col;
+
+	k = w1/k;
+	return Vector3(
+		k*w_col.x,
+		k*w_col.y,
+		k*(w_col.z - w0) );
+}
+
+
+Matrix3
+Perspective::make_alpha_matrix(
+	Real aw0, Real aw1,
+	Real bw0, Real bw1,
+	const Vector3 &w_col )
+{
+	const Vector3 a_col = make_alpha_matrix_col(aw0, aw1, w_col);
+	const Vector3 b_col = make_alpha_matrix_col(bw0, bw1, w_col);
+	return Matrix3(
+		Vector3( a_col.x, b_col.x, w_col.x ),
+		Vector3( a_col.y, b_col.y, w_col.y ),
+		Vector3( a_col.z, b_col.z, w_col.z ) );
+}
+
+
+void
+Perspective::calc_raster_size(
+	Pair2 &out_bounds,
+	IntVector2 &out_raster_size,
+	Matrix3 &out_raster_matrix,
+	const Vector2 &resolution,
+	const Pair2 &bounds )
+{
+	const Vector2 offset = bounds.p0;
+	const Vector2 raster_size_orig(
+		(bounds.p1.x - bounds.p0.x)*resolution.x,
+		(bounds.p1.y - bounds.p0.y)*resolution.y );
+	
+	Vector2 raster_size_float = raster_size_orig;
+	if (raster_size_float.x > max_width)
+		raster_size_float.x = max_width;
+	if (raster_size_float.y > max_height)
+		raster_size_float.y = max_height;
+	IntVector2 raster_size(
+		(int)ceil( raster_size_float.x - real_precision ),
+		(int)ceil( raster_size_float.y - real_precision ) );
+	if (raster_size.x * raster_size.y > max_area) {
+		const Real k = sqrt(Real(max_area)/(raster_size.x * raster_size.y));
+		raster_size.x = std::max(1, (int)floor(raster_size.x*k + real_precision));
+		raster_size.y = std::max(1, (int)floor(raster_size.y*k + real_precision));
+	}
+	
+	const Vector2 new_resolution(
+		raster_size_orig.x > real_precision ? resolution.x * raster_size.x / raster_size_orig.x : resolution.x,
+		raster_size_orig.y > real_precision ? resolution.y * raster_size.y / raster_size_orig.y : resolution.y );
+	
+	out_bounds = bounds.inflated( new_resolution * border_width );
+	out_raster_size = raster_size + IntVector2(border_width, border_width)*2;
+	out_raster_matrix = Matrix3::translation(Vector2(border_width, border_width))
+					  * Matrix3::scaling(new_resolution)
+					  * Matrix3::translation(-offset);
+}
+
+
+void
+Perspective::create_layers(
+	LayerList &out_layers,
+	const Matrix3 &matrix,
+	const IntPair2 &dst_bounds,
+	const Real step )
+{
+	bool is_invertible = false;
+	const Matrix3 back_matrix = matrix.inverted(&is_invertible);
+	if (!is_invertible)
+		return; // matrix is collapsed
+
+	// calc src resolution
+	const Vector2 resolution = calc_optimal_resolution(
+		back_matrix.row_x().vec2(),
+		back_matrix.row_y().vec2() );
+	if (resolution.x <= real_precision || resolution.y <= real_precision)
+		return; // 
+
+	// corners
+	Vector3 dst_corners[4] = {
+		Vector3(dst_bounds.p0.x, dst_bounds.p0.y, 1),
+		Vector3(dst_bounds.p1.x, dst_bounds.p0.y, 1),
+		Vector3(dst_bounds.p1.x, dst_bounds.p1.y, 1),
+		Vector3(dst_bounds.p0.x, dst_bounds.p1.y, 1) };
+	
+	// calc coefficient for equatation of "horizontal" line: A*x + B*y + C = 1/w (aka A*x + B*y = 1/w - C)
+	//                     equatation of line of horizon is: A*x + B*y + C = 0   (aka A*x + B*y = -C)
+	const Real A = back_matrix.m02;
+	const Real B = back_matrix.m12;
+	const Real C = back_matrix.m22;
+	
+	Real hd = sqrt(A*A + B*B);
+	if (hd <= real_precision) {
+		// orthogonal projection - no perspective - no subdiviosions
+		
+		if (fabs(C) < real_precision)
+			return; // only when matrix was not invertible (additional check)
+		
+		// force w coord to one to avoid division to w
+		Real k = 1/C;
+		Matrix3 bm = back_matrix;
+		for(int i = 0; i < 9; ++i) bm.a[i] *= k;
+
+		// calc bounds
+		Pair2 layer_src_bounds(Vector2(INFINITY, INFINITY), Vector2(-INFINITY, -INFINITY));
+		for(int i = 0; i < 4; ++i)
+			layer_src_bounds.expand( (bm*dst_corners[i]).vec2() );
+		if (layer_src_bounds.empty())
+			return;
+		
+		// make layer
+		Layer layer;
+		layer.dst_bounds = dst_bounds;
+		calc_raster_size(
+			layer.src_bounds,
+			layer.src_size,
+			layer.back_matrix,
+			resolution * k,
+			layer_src_bounds );
+		layer.back_matrix *= bm;
+		layer.back_alpha_matrix = Matrix3(
+			Vector3(0, 0, 0),
+			Vector3(0, 0, 0),
+			Vector3(1, 1, 1) );
+		out_layers.push_back(layer);
+
+		return; 
+	}
+
+	// find visible w range
+	hd = 1/hd;
+	const Real horizonw1 = hd;
+	const Real horizonw2 = hd/std::min(Real(2), step);
+	const Real horizonw3 = hd/step;
+	Real maxw = -INFINITY, minw = INFINITY;
+	Vector3 src_corners[4];
+	for(int i = 0; i < 4; ++i) {
+		Vector3 src = back_matrix * dst_corners[i];
+		if (fabs(src.z) > real_precision) {
+			Real w = 1/src.z;
+			if (w > 0 && w < horizonw1)
+				src_corners[i] = Vector3(src.x*w, src.y*w, w);
+			else
+				w = horizonw1;
+			if (minw > w) minw = w;
+			if (maxw < w) maxw = w;
+		}
+	}
+	if (minw >= maxw - real_precision)
+		return; // all bounds too thin
+	const Real maxw3 = std::min(maxw, horizonw3);
+
+	// steps
+	const Real stepLog = log(step);
+	int minlog = (int)floor(log(minw)/stepLog + real_precision);
+	int maxlog = (int)ceil(log(maxw3)/stepLog - real_precision);
+	if (maxlog < minlog) maxlog = minlog;
+	
+	Real w = pow(step, Real(minlog));
+	for(int i = minlog; i <= maxlog; ++i, w *= step) {
+		// w range
+		const Real w0  = w/step;
+		const Real w1  = std::min(w*step, horizonw1);
+		
+		// alpha ranges
+		const Real aw0 = w0;
+		const Real aw1 = w;
+		const Real bw0 = i == maxlog ? horizonw1 : w1;
+		const Real bw1 = i == maxlog ? horizonw2 : w;
+		
+		// calc bounds
+		Pair2 layer_src_bounds(Vector2(INFINITY, INFINITY), Vector2(-INFINITY, -INFINITY));
+		Pair2 layer_dst_bounds(Vector2(INFINITY, INFINITY), Vector2(-INFINITY, -INFINITY));
+		for(int j = 0; j < 4; ++j) {
+			if (src_corners[j].z && src_corners[j].z > w0 && src_corners[j].z < w1) {
+				layer_src_bounds.expand(src_corners[j].vec2());
+				layer_dst_bounds.expand(dst_corners[j].vec2());
+			}
+		}
+		Vector2 line[2];
+		int line_count = truncate_line(line, Pair2(dst_bounds), A, B, 1/w0 - C);
+		for(int j = 0; j < line_count; ++j) {
+			layer_src_bounds.expand( (back_matrix * Vector3(line[j], 1)).vec2() * w0 );
+			layer_dst_bounds.expand( line[j] );
+		}
+		
+		line_count = truncate_line(line, Pair2(dst_bounds), A, B, 1/w1 - C);
+		for(int j = 0; j < line_count; ++j) {
+			layer_src_bounds.expand( (back_matrix * Vector3(line[j], 1)).vec2() * w1 );
+			layer_dst_bounds.expand( line[j] );
+		}
+		
+		if (layer_src_bounds.empty() || layer_dst_bounds.empty())
+			continue;
+		
+		// make layer
+		Layer layer;
+		layer.dst_bounds = IntPair2(
+			IntVector2( std::max(dst_bounds.p0.x, (int)floor(layer_dst_bounds.p0.x + real_precision)),
+						std::max(dst_bounds.p0.y, (int)floor(layer_dst_bounds.p0.y + real_precision)) ),
+			IntVector2( std::min(dst_bounds.p1.x, (int)ceil (layer_dst_bounds.p1.x - real_precision)),
+						std::min(dst_bounds.p1.y, (int)ceil (layer_dst_bounds.p1.y - real_precision)) ));
+		calc_raster_size(
+			layer.src_bounds,
+			layer.src_size,
+			layer.back_matrix,
+			resolution * w,
+			layer_src_bounds );
+		layer.back_matrix *= back_matrix;
+		layer.back_alpha_matrix = make_alpha_matrix(
+			aw0, aw1, bw0, bw1,
+			layer.back_matrix.get_col(2) );
+		out_layers.push_back(layer);
+	}
+}
+
+
+void
+Perspective::add_premulted(
+	const Layer &layer,
+	const RefPtr<Surface> &src_surface,
+	const RefPtr<Surface> &dst_surface )
+{
+	if (!src_surface || src_surface->empty()) return;
+	if (!dst_surface || dst_surface->empty()) return;
+	
+	const int x0 = std::max(layer.dst_bounds.p0.x, 0);
+	const int y0 = std::max(layer.dst_bounds.p0.y, 0);
+	const int x1 = std::min(layer.dst_bounds.p1.x, dst_surface->width());
+	const int y1 = std::min(layer.dst_bounds.p1.y, dst_surface->width());
+	if (x0 >= x1 || y0 >= y1) return;
+	
+	const int w = x1 - x0;
+	const int h = y1 - y0;
+	if (w <= 0 || h <= 0) return;
+	
+	const Vector3 coord_dx = layer.back_matrix.row_x();
+	const Vector3 coord_dy = layer.back_matrix.row_y() + coord_dx*(x0 - x1);
+	Vector3 coord = layer.back_matrix * Vector3(x0, y0, 1);
+
+	const Vector2 alpha_dx = layer.back_alpha_matrix.row_x().vec2();
+	const Vector2 alpha_dy = layer.back_alpha_matrix.row_y().vec2() + alpha_dx*(x0 - x1);
+	Vector2 alpha = (layer.back_alpha_matrix * Vector3(x0, y0, 1)).vec2();
+	
+	const int dr = dst_surface->pitch() - x1 + x0;
+	Color *c = &dst_surface->row(y0)[x0];
+	for(int r = y0; r < y1; ++r, c += dr, coord += coord_dy, alpha += alpha_dy) {
+		for(Color *end = c + w; c != end; ++c, coord += coord_dx, alpha += alpha_dx) {
+			if (coord.z > real_precision) {
+				const Real w = 1/coord.z;
+				const Real a = clamp(alpha.x*w, 0, 1) * clamp(alpha.y*w, 0, 1);
+				if (a > real_precision)
+					*c += src_surface->get_pixel_premulted(coord.x*w, coord.y*w) * a;
+			}
+		}
+	}
+}
+
diff --git a/c++/perspective/src/perspective.h b/c++/perspective/src/perspective.h
new file mode 100644
index 0000000..7c96821
--- /dev/null
+++ b/c++/perspective/src/perspective.h
@@ -0,0 +1,71 @@
+#ifndef PERSPECTIVE_H
+#define PERSPECTIVE_H
+
+#include <vector>
+
+#include "vector.h"
+#include "matrix.h"
+#include "surface.h"
+
+
+class Perspective {
+public:
+	class Layer {
+	public:
+		IntPair2   dst_bounds;
+		Pair2      src_bounds;
+		IntVector2 src_size;
+		Matrix3    back_matrix;
+		Matrix3    back_alpha_matrix;
+	};
+	
+	typedef std::vector<Layer> LayerList;
+	
+	
+	static int truncate_line(
+		Vector2 *out_points,
+		const Pair2 &bounds,
+		Real a,
+		Real b,
+		Real c );
+
+	static Vector2 calc_optimal_resolution(
+		const Vector2 &ox,
+		const Vector2 &oy );
+
+	static Vector3 make_alpha_matrix_col(
+		Real w0,
+		Real w1,
+		const Vector3 &w_col );
+
+	static Matrix3 make_alpha_matrix(
+		Real aw0, Real aw1,
+		Real bw0, Real bw1,
+		const Vector3 &w_col );
+	
+	static void calc_raster_size(
+		Pair2 &out_bounds,
+		IntVector2 &out_raster_size,
+		Matrix3 &out_raster_matrix,
+		const Vector2 &resolution,
+		const Pair2 &bounds );
+
+	static void create_layers(
+		LayerList &out_layers,
+		const Matrix3 &matrix,
+		const IntPair2 &dst_bounds,
+		const Real step = 4.0 );
+	
+	static void add_premulted(
+		const Layer &layer,
+		const RefPtr<Surface> &src_surface,
+		const RefPtr<Surface> &dst_surface );
+};
+
+
+
+
+
+
+
+#endif
diff --git a/c++/perspective/src/perspectivesandboxview.cpp b/c++/perspective/src/perspectivesandboxview.cpp
index b317084..aee0140 100644
--- a/c++/perspective/src/perspectivesandboxview.cpp
+++ b/c++/perspective/src/perspectivesandboxview.cpp
@@ -1,4 +1,6 @@
 
+#include <iostream>
+
 #include "perspectivesandboxview.h"
 
 
@@ -71,7 +73,108 @@ PerspectiveSandBoxView::draw_line(
 	Real c,
 	const Color &color )
 {
+	// equatation of line is a*x + b*y = c
+
+	int count = 0;
+	Vector2 points[4];
 	
+	if (fabs(a) > real_precision) {
+		Real x0 = (c - bounds.p0.y*b)/a;
+		if ( x0 + real_precision >= bounds.p0.x
+		  && x0 - real_precision <= bounds.p1.x )
+			points[count++] = Vector2(x0, bounds.p0.y);
+		
+		Real x1 = (c - bounds.p1.y*b)/a;
+		if ( x1 + real_precision >= bounds.p0.x
+		  && x1 - real_precision <= bounds.p1.x )
+			points[count++] = Vector2(x1, bounds.p1.y);
+	}
+
+	if (fabs(b) > real_precision) {
+		Real y0 = (c - bounds.p0.x*a)/b;
+		if ( y0 + real_precision >= bounds.p0.y
+		  && y0 - real_precision <= bounds.p1.y )
+			points[count++] = Vector2(bounds.p0.x, y0);
+		
+		Real y1 = (c - bounds.p1.x*a)/b;
+		if ( y1 + real_precision >= bounds.p0.y
+		  && y1 - real_precision <= bounds.p1.y )
+			points[count++] = Vector2(bounds.p1.x, y1);
+	}
+
+	if (count >= 2) {
+		context->set_source_rgba(color.r, color.g, color.b, color.a);
+		context->move_to(points[0].x, points[0].y);
+		context->line_to(points[1].x, points[1].y);
+		context->stroke();
+	}
+}
+
+void
+PerspectiveSandBoxView::draw_subdivisions(
+	const Cairo::RefPtr<Cairo::Context> &context,
+	const Matrix &matrix,
+	const Pair2 &bounds,
+	const Color &color )
+{
+	const Real ps = get_pixel_size();
+	
+	Vector4 a = matrix.row_x();
+	Vector4 b = matrix.row_y();
+	Vector4 c = matrix.row_w();
+
+	// calc coefficient for equatation of "horizontal" line: A*x + B*y = C + D/w
+	//                     equatation of line of horizon is: A*x + B*y = C
+	Real A = a.y*b.w - a.w*b.y;
+	Real B = a.w*b.x - a.x*b.w;
+	Real C = a.y*b.x - a.x*b.y;
+	Real D = a.w*(b.x*c.y - b.y*c.x) + b.w*(a.y*c.x - a.x*c.y) - C*c.w;
+	
+	if (D < 0) {
+		A = -A;
+		B = -B;
+		C = -C;
+		D = -D;
+	}
+
+	Real hd = ps*sqrt(A*A + B*B);
+	if (hd <= real_precision)
+		return; // orthogonal projection - no prespective - no subdiviosions
+	hd = D/hd;
+	Real horizonw  = hd;
+	Real horizonw2 = hd/4;
+
+	Real maxw = -INFINITY, minw = INFINITY;
+	Vector2 corners[] = {
+		Vector2(bounds.p0.x, bounds.p0.y),
+		Vector2(bounds.p1.x, bounds.p0.y),
+		Vector2(bounds.p1.x, bounds.p1.y),
+		Vector2(bounds.p0.x, bounds.p1.y) };
+	for(int i = 0; i < 4; ++i) {
+		Real w = A*corners[i].x + B*corners[i].y - C;
+		if (fabs(w) > real_precision) {
+			w = D/w;
+			if (w < 0 || w > horizonw) w = horizonw;
+			if (minw > w) minw = w;
+			if (maxw < w) maxw = w;
+		}
+	}
+
+	if (minw >= maxw - real_precision)
+		return; // all bounds too thin
+	Real maxw2 = maxw < horizonw2 ? maxw : horizonw2;
+
+	// draw limits
+	Color limits_color = color;
+	limits_color.a *= 0.5;
+	draw_line(context, get_bounds(), A, B, C + D/minw, limits_color);
+	draw_line(context, get_bounds(), A, B, C + D/maxw, limits_color);
+
+	// split
+	int minlog = (int)ceil (log2(minw ) + real_precision);
+	int maxlog = (int)floor(log2(maxw2) - real_precision);
+	for(int i = minlog; i <= maxlog; ++i)
+		draw_line(context, bounds, A, B, C + D/exp2(i), color);
 }
 
 void
@@ -100,6 +203,12 @@ PerspectiveSandBoxView::on_draw_view(const Cairo::RefPtr<Cairo::Context> &contex
 		Real aw = B.x*C.y - B.y*C.x;
 		Real bw = A.y*C.x - A.x*C.y;
 		
+		if (cw < 0) {
+			aw = -aw;
+			bw = -bw;
+			cw = -cw;
+		}
+		
 		Vector2 c = p0*cw;
 		Vector2 a = px*(cw + aw) - c;
 		Vector2 b = py*(cw + bw) - c;
@@ -110,6 +219,7 @@ PerspectiveSandBoxView::on_draw_view(const Cairo::RefPtr<Cairo::Context> &contex
 	}
 
 	draw_grid(context, matrix, Color(0, 1, 0, 1));
+	draw_subdivisions(context, matrix, bounds, Color(0, 0, 1, 1));
 
 	// draw frames
 	context->set_source_rgba(0, 0, 1, 0.75);
diff --git a/c++/perspective/src/perspectivesandboxview.h b/c++/perspective/src/perspectivesandboxview.h
index 2a05038..45afde1 100644
--- a/c++/perspective/src/perspectivesandboxview.h
+++ b/c++/perspective/src/perspectivesandboxview.h
@@ -29,6 +29,12 @@ public:
 		Real c,
 		const Color &color );
 
+	void draw_subdivisions(
+		const Cairo::RefPtr<Cairo::Context> &context,
+		const Matrix &matrix,
+		const Pair2 &bounds,
+		const Color &color );
+
 	void on_draw_view(const Cairo::RefPtr<Cairo::Context> &context);
 };
 
diff --git a/c++/perspective/src/surface.cpp b/c++/perspective/src/surface.cpp
index 55124a1..90bf4cd 100644
--- a/c++/perspective/src/surface.cpp
+++ b/c++/perspective/src/surface.cpp
@@ -14,6 +14,95 @@ Surface::reset() {
 	m_origin = NULL;
 }
 
+void
+Surface::mult_alpha() {
+	if (empty()) return;
+	const int h = height();
+	const int w = width();
+	const int dr = pitch() - width();
+	Color *c = m_origin;
+	for(int r = 0; r < h; ++r, c += dr)
+		for(Color *end = c + w; c != end; ++c)
+			c->mult_alpha();
+}
+
+void
+Surface::demult_alpha() {
+	if (empty()) return;
+	const int h = height();
+	const int w = width();
+	const int dr = pitch() - width();
+	Color *c = m_origin;
+	for(int r = 0; r < h; ++r, c += dr)
+		for(Color *end = c + w; c != end; ++c)
+			c->demult_alpha();
+}
+
+Color
+Surface::get_pixel(Real x, Real y) const {
+	// linear interpolation
+	if ( empty()
+	  || x < 0 || x > width() - 1
+	  || y < 0 || y > height() - 1 )
+		return Color();
+	
+	int x0 = (int)floor(x), x1 = x0 + 1;
+	int y0 = (int)floor(y), y1 = y0 + 1;
+	Real px = x - x0, ppx = 1.0 - px;
+	Real py = y - y0, ppy = 1.0 - py;
+	
+	x0 = iclamp(x0, 0, width() - 1);
+	x1 = iclamp(x1, 0, width() - 1);
+	y0 = iclamp(y0, 0, height() - 1);
+	y1 = iclamp(y1, 0, height() - 1);
+
+	Real k[] = { ppx*ppy, px*ppy, ppx*py, px*py };
+	const Color *corners[] = { &row(y0)[x0], &row(y0)[x1], &row(y1)[x0], &row(y1)[x1] };
+	Color sum;
+	for(int i = 0; i < 4; ++i) {
+		const Color &c = *corners[i];
+		Real a = c.a * k[i];
+		for(int j = 0; j < 3; ++j)
+			sum.channels[j] += c.channels[j]*a;
+		sum.a += a;
+	}
+	
+	if (sum.a <= real_precision)
+		return Color();
+
+	Real aa = 1.0/sum.a;
+	for(int j = 0; j < 3; ++j)
+		sum.channels[j] *= aa;
+	return sum;
+}
+
+Color
+Surface::get_pixel_premulted(Real x, Real y) const {
+	// linear interpolation
+	if ( empty()
+	  || x < 0 || x > width() - 1
+	  || y < 0 || y > height() - 1 )
+		return Color();
+	
+	int x0 = (int)floor(x), x1 = x0 + 1;
+	int y0 = (int)floor(y), y1 = y0 + 1;
+	Real px = x - x0, ppx = 1.0 - px;
+	Real py = y - y0, ppy = 1.0 - py;
+	
+	x0 = iclamp(x0, 0, width() - 1);
+	x1 = iclamp(x1, 0, width() - 1);
+	y0 = iclamp(y0, 0, height() - 1);
+	y1 = iclamp(y1, 0, height() - 1);
+
+	Color sum;
+	Real k[] = { ppx*ppy, px*ppy, ppx*py, px*py };
+	const Color *corners[] = { &row(y0)[x0], &row(y0)[x1], &row(y1)[x0], &row(y1)[x1] };
+	for(int i = 0; i < 4; ++i)
+		for(int j = 0; j < 4; ++j)
+			sum.channels[j] += corners[i]->channels[j]*k[i];
+	return sum;
+}
+
 void Surface::copy(int width, int height, const Color *src_origin, Color *dest_origin, int src_pitch, int dest_pitch) {
 #ifndef NDEBUG
 	{ // check source and destination ranges
@@ -68,7 +157,7 @@ Surface::to_pixbuf() const {
 }
 
 Cairo::RefPtr<Cairo::ImageSurface>
-Surface::to_cairo_surface() const {
+Surface::to_cairo_surface(bool premulted) const {
 	if (empty())
 		return Cairo::RefPtr<Cairo::ImageSurface>();
 	Cairo::RefPtr<Cairo::ImageSurface> cairo_surface =
@@ -81,7 +170,8 @@ Surface::to_cairo_surface() const {
 		unsigned char *pixel = pixels;
 		pixels += stride;
 		for(int c = 0; c < width(); ++c) {
-			const Color &color = src_row[c];
+			Color color = src_row[c];
+			if (!premulted) color.mult_alpha();
 			*(pixel++) = (guint8)(color.r*255.9);
 			*(pixel++) = (guint8)(color.g*255.9);
 			*(pixel++) = (guint8)(color.b*255.9);
diff --git a/c++/perspective/src/surface.h b/c++/perspective/src/surface.h
index c149f2f..e1f36db 100644
--- a/c++/perspective/src/surface.h
+++ b/c++/perspective/src/surface.h
@@ -22,7 +22,7 @@ protected:
 	Surface(const Surface &);
 	Surface& operator= (const Surface &) { return *this; }
 	void reset();
-
+	
 public:
 	virtual ~Surface() { }
 	
@@ -34,7 +34,7 @@ public:
 	
 	inline Color* origin() { return m_origin; }
 	inline const Color* origin() const { return m_origin; }
-
+	
 	inline int data_count()   const { return data_count  (width(), height(), pitch()); }
 	inline size_t data_size() const { return data_size   (width(), height(), pitch()); }
 	inline int begin_offset() const { return begin_offset(width(), height(), pitch()); }
@@ -52,7 +52,12 @@ public:
 	
 	inline Color* operator[] (int index) { return row(index); }
 	inline const Color* operator[] (int index) const { return row(index); }
-
+	
+	void mult_alpha();
+	void demult_alpha();
+	Color get_pixel(Real x, Real y) const;
+	Color get_pixel_premulted(Real x, Real y) const;
+	
 	static inline int data_count(int width, int height, int pitch)
 		{ return width + abs(pitch)*(height - 1); }
 	static inline size_t data_size(int width, int height, int pitch)
@@ -63,9 +68,9 @@ public:
 		{ return pitch > 0 ? data_count(width, height, pitch) : width; }
 
 	static void copy(int width, int height, const Color *src_origin, Color *dest_origin, int src_pitch = 0, int dest_pitch = 0);
-
+	
 	Glib::RefPtr<Gdk::Pixbuf> to_pixbuf() const;
-	Cairo::RefPtr<Cairo::ImageSurface> to_cairo_surface() const;
+	Cairo::RefPtr<Cairo::ImageSurface> to_cairo_surface(bool premulted = false) const;
 };
 
 
diff --git a/c++/perspective/src/vector.h b/c++/perspective/src/vector.h
index 6e2b563..fd19895 100644
--- a/c++/perspective/src/vector.h
+++ b/c++/perspective/src/vector.h
@@ -270,7 +270,18 @@ public:
 		p0(p), p1(p) { }
 	inline PairT(const Vector &p0, const Vector &p1):
 		p0(p0), p1(p1) { }
-	
+
+	template<typename TT>
+	inline explicit PairT(const VectorT<TT> &p):
+		p0(p), p1(p) { }
+	template<typename TT>
+	inline explicit PairT(const VectorT<TT> &p0, const VectorT<TT> &p1):
+		p0(p0), p1(p1) { }
+
+	template<typename TT>
+	inline explicit PairT(const PairT<TT> &other):
+		PairT(other.p0, other.p1) { }
+		
 	inline bool operator< (const PairT &other) const {
 		return p0 < other.p0 ? true
 			 : other.p0 < p0 ? false
@@ -293,13 +304,34 @@ public:
 				return true;
 		return false;
 	}
-	inline void sort() {
+	inline PairT& sort() {
 		for(int i = 0; i < Vector::Count; ++i)
-			if (!Vector::coord_less(p0[i], p1[i]))
+			if (p0[i] < p1[i])
 				std::swap(p0[i], p1[i]);
+		return *this;
 	}
 	inline PairT sorted() const
-		{ PairT x(*this); x.sort(); return x; }
+		{ return PairT(*this).sort(); }
+
+	inline PairT& expand(const Vector &p) {
+		if (empty())
+			return *this = PairT(p);
+		for(int i = 0; i < Vector::Count; ++i) {
+			if (p0[i] < p[i]) p0[i] = p[i];
+			if (p[i] < p1[i]) p1[i] = p[i];
+		}
+		return *this;
+	}
+	inline PairT expanded(const Vector &p) const
+		{ return PairT(*this).expand(p); }
+
+	inline PairT& inflate(const Vector &p) {
+		if (empty()) return *this;
+		p0 -= p; p1 += p;
+		return *this;
+	}
+	inline PairT inflated(const Vector &p) const
+		{ return PairT(*this).inflate(p); }
 };
 
 
diff --git a/c++/perspective/src/view.cpp b/c++/perspective/src/view.cpp
index c477133..cb05f0b 100644
--- a/c++/perspective/src/view.cpp
+++ b/c++/perspective/src/view.cpp
@@ -37,9 +37,17 @@ View::View():
 View::~View() { }
 
 Real
-View::get_pixel_size()
+View::get_pixel_size() const
 	{ return 1.0/sqrt(fabs(transform.m00 * transform.m11)); }
 
+Pair2
+View::get_bounds() const {
+	Matrix m = transform_from_pixels();
+	return Pair2(
+		(m*Vector4(0          , 0           , 0, 1)).vec2(),
+		(m*Vector4(get_width(), get_height(), 0, 1)).vec2() );
+}
+
 Matrix
 View::transform_to_pixels() const {
 	Matrix matrix = transform;
diff --git a/c++/perspective/src/view.h b/c++/perspective/src/view.h
index 198087e..a5070f8 100644
--- a/c++/perspective/src/view.h
+++ b/c++/perspective/src/view.h
@@ -64,7 +64,8 @@ public:
 	void point_changed(const PointPtr &point)
 		{ point_motion(point); signal_point_changed(point); }
 
-	Real get_pixel_size();
+	Real get_pixel_size() const;
+	Pair2 get_bounds() const;
 	
 	Matrix transform_to_pixels() const;
 	Matrix transform_from_pixels() const;