diff --git a/c++/perspective/src/generator.cpp b/c++/perspective/src/generator.cpp
index 48bc1f7..1e62fd4 100644
--- a/c++/perspective/src/generator.cpp
+++ b/c++/perspective/src/generator.cpp
@@ -1,5 +1,7 @@
 #include <ctime>
 
+#include <iostream>
+
 #include "generator.h"
 
 
@@ -7,7 +9,7 @@ Generator::Generator(): Generator(random(time(NULL))) { }
 
 Generator::Generator(ULongInt seed):
 	one(1ull << 32),
-	max_cache_count(8*1024*1024)
+	max_cache_count(4ll*1024*1024*1024/sizeof(CacheEntry))
 {
 	set_seed(seed);
 	LongIntVector2 a(1, 2);
@@ -128,8 +130,14 @@ Generator::generate(Surface &target, const Pair2 &bounds) const {
 	d.y /= target.height();
 	Vector2 precision = d/2.0;
 	
+	int percent = 0;
+	Real pk = Real(100)/target.height();
+	
 	Vector2 p = bounds.p0;
 	for(int r = 0; r < target.height(); ++r, p.x = bounds.p0.x, p.y += d.y) {
+		Real pr = r*pk;
+		while(percent < pr)
+			std::cout << (++percent) << "%" << std::endl;
 		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/main.cpp b/c++/perspective/src/main.cpp
index ab69fdd..6ce8f08 100644
--- a/c++/perspective/src/main.cpp
+++ b/c++/perspective/src/main.cpp
@@ -13,9 +13,11 @@ int main(int argc, char **argv) {
 		argc, argv, "org.coolbug.lab.perspective");
 	
 	std::cout << "create window" << std::endl;
-	Gtk::Window window;
-	window.add(*Gtk::manage(new PerspectiveSandBoxView()));
-	window.show_all();
+	//Gtk::Window window;
+	//window.add(*Gtk::manage(new PerspectiveSandBoxView()));
+	//window.add(*Gtk::manage(new PerspectiveSandBoxView()));
+	//window.show_all();
+	MainWindow window;
 	
 	std::cout << "run" << std::endl;
 	int result = application->run(window);
diff --git a/c++/perspective/src/mainwindow.cpp b/c++/perspective/src/mainwindow.cpp
index 3057964..80afd6e 100644
--- a/c++/perspective/src/mainwindow.cpp
+++ b/c++/perspective/src/mainwindow.cpp
@@ -1,258 +1,13 @@
 #include <iostream>
-#include <utility>
 
 #include <gdkmm/pixbuf.h>
 #include <gtkmm/hvbox.h>
 #include <gtkmm/image.h>
 
-#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();
-	}
-}
-
+#include "log.h"
+#include "perspective.h"
 
+#include "mainwindow.h"
 
 
 MainWindow::MainWindow():
@@ -310,8 +65,8 @@ MainWindow::MainWindow():
 	add(*hbox);
 	show_all();
 	
-	update_view_surface();
-	update_layers();
+	update_src_surface();
+	update_dst_surface();
 }
 
 Gtk::Widget*
@@ -339,7 +94,9 @@ MainWindow::generator_demo() {
 }
 
 void
-MainWindow::update_view_surface() {
+MainWindow::update_src_surface() {
+	std::cout << "update_src_surface() - begin" << std::endl;
+
 	Matrix4 transform = src_view.transform.inverted();
 	int w = src_view.get_width();
 	int h = src_view.get_height();
@@ -349,65 +106,92 @@ MainWindow::update_view_surface() {
 		(transform * Vector4(w, h, 0, 1)).vec2() );
 	
 	if (w <= 0 || h <= 0) {
-		view_surface.clear();
-	} else
-	if ( !view_surface
-	  || view_surface->get_width() != w
-	  || view_surface->get_height() != h
-	  || view_surface_rect != rect )
-	{
+		src_surface.clear();
+	} else {
 		DataSurface surface(w, h);
 		generator.generate(surface, rect);
-		view_surface = surface.to_cairo_surface();
-		view_surface_rect = rect;
+		src_surface = surface.to_cairo_surface();
 	}
+
+	std::cout << "update_src_surface() - end" << std::endl;
 }
 
 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));
+MainWindow::update_dst_surface() {
+	std::cout << "update_dst_surface() - begin" << std::endl;
+
+	dst_surface.clear();
 
-	Matrix in_matrix = make_plain_matrix( Pair2(
+	// src matrix
+	const Pair2 src_bounds(
 		src_rect_p0->position,
-		src_rect_p1->position ));
-	Matrix out_matrix = make_perspective_matrix(
+		src_rect_p1->position );
+	const Vector2 src_dist = src_bounds.distance();
+	if (fabs(src_dist.x) < real_precision || fabs(src_dist.y) < real_precision)
+		return;
+	
+	Real kx = 1/src_dist.x;
+	Real ky = 1/src_dist.x;
+	Real x0 = -src_bounds.p0.x*kx;
+	Real x1 = -src_bounds.p0.y*kx;
+	Matrix3 src_matrix(
+		Vector3(kx,  0, 0),
+		Vector3( 0, ky, 0),
+		Vector3(x0, x1, 1) );
+	
+	// dst matrix
+	const Matrix dst_transform = dst_view.transform_to_pixels();
+	const Vector2 dst_rect_p0_pixels(dst_transform * Vector4(dst_rect_p0->position, 0.0, 1.0));
+	const Vector2 dst_rect_px_pixels(dst_transform * Vector4(dst_rect_px->position, 0.0, 1.0));
+	const Vector2 dst_rect_py_pixels(dst_transform * Vector4(dst_rect_py->position, 0.0, 1.0));
+	const Vector2 dst_rect_p1_pixels(dst_transform * Vector4(dst_rect_p1->position, 0.0, 1.0));
+	Matrix3 dst_matrix = Perspective::make_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();
+	// final perspective matrix
+	Matrix3 matrix = dst_matrix * src_matrix;
+
+	// dst bounds
+	const Pair2 dst_bounds_pixels_float(
+		(dst_transform * Vector4(dst_bounds_p0->position, 0.0, 1.0)).vec2(),
+		(dst_transform * Vector4(dst_bounds_p1->position, 0.0, 1.0)).vec2() );
+	const IntPair2 dst_bounds_pixels(
+		IntVector2( (int)floor(dst_bounds_pixels_float.p0.x + real_precision),
+					(int)floor(dst_bounds_pixels_float.p0.y + real_precision) ),
+		IntVector2( (int)ceil (dst_bounds_pixels_float.p1.x - real_precision),
+					(int)ceil (dst_bounds_pixels_float.p1.y - real_precision) ));
+	if (dst_bounds_pixels.empty())
+		return;
 
-	LayerList layers;
-	make_layers(bounds, matrix, layers);
+	std::cout << Log::tab() << "dst_bounds_pixels: " << Log::to_string(dst_bounds_pixels) << std::endl;
+	std::cout << Log::tab() << "matrix:" << std::endl
+			  << Log::to_string(matrix, 10, Log::tab(2));
+
+	// make layers
+	Perspective::LayerList layers;
+	Perspective::create_layers(
+		layers,
+		matrix,
+		dst_bounds_pixels );
 	
-	for(LayerList::iterator i = layers.begin(); i != layers.end(); ++i)
-		make_surface(generator, *i);
+	std::cout << Log::tab() << "created layers:" << std::endl;
+	Perspective::print_layers(layers, Log::tab(2));
+	std::cout << Log::tab() << "---------------" << std::endl;
 	
-	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();
-			}
-		}
+	// make surface
+	DataSurface surface(dst_bounds_pixels.p1.x, dst_bounds_pixels.p1.y);
+	for(Perspective::LayerList::const_iterator i = layers.begin(); i != layers.end(); ++i) {
+		DataSurface layer_surface(i->src_size.x, i->src_size.y);
+		generator.generate(layer_surface, i->src_bounds);
+		layer_surface.mult_alpha();
+		Perspective::add_premulted(*i, layer_surface, surface);
 	}
+	dst_surface = surface.to_cairo_surface(true);
+
+	std::cout << "update_dst_surface() - end" << std::endl;
 }
 
 void
@@ -417,13 +201,13 @@ MainWindow::on_point_motion(const View::PointPtr &/*point*/) {
 
 void
 MainWindow::on_point_changed(const View::PointPtr &/*point*/) {
-	update_layers();
+	update_dst_surface();
 }
 
 void
 MainWindow::on_view_transform_changed(View *view) {
 	if (view == &src_view) {
-		update_view_surface();
+		update_src_surface();
 		view->queue_draw();
 	}
 }
@@ -434,14 +218,14 @@ MainWindow::on_draw_src_view(const Cairo::RefPtr<Cairo::Context> &context) {
 
 	context->save();
 
-	if ( !view_surface
-	  || src_view.get_width() != view_surface->get_width()
-	  || src_view.get_height() != view_surface->get_height() )
-		update_view_surface();
-	if (view_surface) {
+	if ( !src_surface
+	  || src_view.get_width() != src_surface->get_width()
+	  || src_view.get_height() != src_surface->get_height() )
+		update_src_surface();
+	if (src_surface) {
 		context->save();
 		context->transform( src_view.transform_from_pixels().to_cairo() );
-		context->set_source(view_surface, 0, 0);
+		context->set_source(src_surface, 0, 0);
 		context->paint();
 		context->restore();
 	}
@@ -463,10 +247,10 @@ MainWindow::on_draw_dst_view(const Cairo::RefPtr<Cairo::Context> &context) {
 	const Real ps = src_view.get_pixel_size();
 
 	context->save();
-	if (target_surface) {
+	if (dst_surface) {
 		context->save();
-		context->transform( dst_view.transform.inverted().to_cairo() );
-		context->set_source(target_surface, 0, 0);
+		context->transform( dst_view.transform_from_pixels().to_cairo() );
+		context->set_source(dst_surface, 0, 0);
 		context->paint();
 		context->restore();
 	}
diff --git a/c++/perspective/src/mainwindow.h b/c++/perspective/src/mainwindow.h
index 8e8ac63..127406c 100644
--- a/c++/perspective/src/mainwindow.h
+++ b/c++/perspective/src/mainwindow.h
@@ -13,22 +13,6 @@
 
 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;
@@ -39,17 +23,15 @@ public:
 		dst_rect_p0, dst_rect_px, dst_rect_py, dst_rect_p1,
 		dst_bounds_p0, dst_bounds_p1;
 	
-	Pair2 view_surface_rect;
-	Cairo::RefPtr<Cairo::ImageSurface> view_surface;
-
-	Cairo::RefPtr<Cairo::ImageSurface> target_surface;
+	Cairo::RefPtr<Cairo::ImageSurface> src_surface;
+	Cairo::RefPtr<Cairo::ImageSurface> dst_surface;
 
 	MainWindow();
 	
 	Gtk::Widget* generator_demo();
 	
-	void update_layers();
-	void update_view_surface();
+	void update_src_surface();
+	void update_dst_surface();
 	
 	void on_point_motion(const View::PointPtr &point);
 	void on_point_changed(const View::PointPtr &point);
diff --git a/c++/perspective/src/perspective.cpp b/c++/perspective/src/perspective.cpp
index 7af2f54..ad7f7d0 100644
--- a/c++/perspective/src/perspective.cpp
+++ b/c++/perspective/src/perspective.cpp
@@ -1,19 +1,54 @@
 
+#include <iostream>
+#include <iomanip>
+
+#include "log.h"
+
 #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 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;
 }
 
 
+Matrix3
+Perspective::make_matrix(
+	const Vector2 &p0,
+	const Vector2 &px,
+	const Vector2 &py,
+	const Vector2 &p1 )
+{
+	Vector2 A = px - p1;
+	Vector2 B = py - p1;
+	Vector2 C = p0 + p1 - px - py;
+	
+	Real cw = A.y*B.x - A.x*B.y;
+	Real aw = B.x*C.y - B.y*C.x;
+	Real bw = A.y*C.x - A.x*C.y;
+	
+	// normalize and force cw to be positive
+	Real k = aw*aw + bw*bw + cw*cw;
+	k = k > real_precision ? 1/sqrt(k) : 1;
+	if (cw < 0) k = -k;
+	aw *= k;
+	bw *= k;
+	cw *= k;
+	
+	Vector2 c = p0*cw;
+	Vector2 a = px*(cw + aw) - c;
+	Vector2 b = py*(cw + bw) - c;
+
+	return Matrix3( Vector3(a, aw),
+					Vector3(b, bw),
+					Vector3(c, cw) );
+}
+
 
 int
 Perspective::truncate_line(
@@ -164,8 +199,11 @@ Perspective::calc_raster_size(
 	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 );
+	const Vector2 new_border(
+		new_resolution.x > real_precision ? border_width/new_resolution.x : Real(0),
+		new_resolution.y > real_precision ? border_width/new_resolution.y : Real(0) );
 	
-	out_bounds = bounds.inflated( new_resolution * border_width );
+	out_bounds = bounds.inflated( new_border );
 	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)
@@ -181,16 +219,24 @@ Perspective::create_layers(
 	const Real step )
 {
 	bool is_invertible = false;
-	const Matrix3 back_matrix = matrix.inverted(&is_invertible);
+	Matrix3 back_matrix = matrix.inverted(&is_invertible);
 	if (!is_invertible)
 		return; // matrix is collapsed
+	
+	{ // normalize matrix by w-column
+		Real k = back_matrix.get_col(2).length();
+		if (k <= real_precision)
+			return; // matrix is (almost) not invertible
+		k = 1/k;
+		for(int i = 0; i < 9; ++i) back_matrix.a[i] *= k;
+	}
 
 	// 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; // 
+		return; // cannot calc resolution, this can happen if matrix is (almost) not invertible
 
 	// corners
 	Vector3 dst_corners[4] = {
@@ -212,18 +258,16 @@ Perspective::create_layers(
 		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() );
+			layer_src_bounds.expand( (back_matrix*dst_corners[i]).vec2() );
 		if (layer_src_bounds.empty())
 			return;
 		
+		std::cout << "debug: layer_src_bounds: " << Log::to_string(layer_src_bounds) << std::endl;
+		std::cout << "debug: back_matrix:"<< std::endl << Log::to_string(back_matrix, 10);
+		
 		// make layer
 		Layer layer;
 		layer.dst_bounds = dst_bounds;
@@ -231,9 +275,9 @@ Perspective::create_layers(
 			layer.src_bounds,
 			layer.src_size,
 			layer.back_matrix,
-			resolution * k,
+			resolution,
 			layer_src_bounds );
-		layer.back_matrix *= bm;
+		layer.back_matrix *= back_matrix;
 		layer.back_alpha_matrix = Matrix3(
 			Vector3(0, 0, 0),
 			Vector3(0, 0, 0),
@@ -334,16 +378,16 @@ Perspective::create_layers(
 void
 Perspective::add_premulted(
 	const Layer &layer,
-	const RefPtr<Surface> &src_surface,
-	const RefPtr<Surface> &dst_surface )
+	Surface &src_surface,
+	Surface &dst_surface )
 {
-	if (!src_surface || src_surface->empty()) return;
-	if (!dst_surface || dst_surface->empty()) return;
+	if (src_surface.empty()) return;
+	if (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());
+	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.height());
 	if (x0 >= x1 || y0 >= y1) return;
 	
 	const int w = x1 - x0;
@@ -358,17 +402,39 @@ Perspective::add_premulted(
 	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];
+	const int dr = dst_surface.pitch() - x1 + x0;
+	Color *c = &dst_surface[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;
+					*c += src_surface.get_pixel_premulted(coord.x*w, coord.y*w) * a;
 			}
 		}
 	}
 }
 
+
+void
+Perspective::print_layer(const Layer &layer, const std::string &prefix) {
+	const int w = 10;
+	std::cout << prefix << "dst_bounds: " << Log::to_string(layer.dst_bounds, w) << std::endl;
+	std::cout << prefix << "src_bounds: " << Log::to_string(layer.dst_bounds, w) << std::endl;
+	std::cout << prefix << "src_size:   " << Log::to_string(layer.src_size,   w) << std::endl;
+	std::cout << prefix << "back_matrix:" << std::endl
+			  << Log::to_string(layer.back_matrix, w, prefix + Log::tab()); 
+	std::cout << prefix << "back_alpha_matrix:" << std::endl
+			  << Log::to_string(layer.back_alpha_matrix, w, prefix + Log::tab()); 
+}
+
+
+void
+Perspective::print_layers(const LayerList &layers, const std::string &prefix) {
+	int index = 0;
+	for(LayerList::const_iterator i = layers.begin(); i < layers.end(); ++i, ++index) {
+		std::cout << prefix << "layer #" << index << std::endl;
+		print_layer(*i, prefix + Log::tab());
+	}
+}
diff --git a/c++/perspective/src/perspective.h b/c++/perspective/src/perspective.h
index 7c96821..64c4fc5 100644
--- a/c++/perspective/src/perspective.h
+++ b/c++/perspective/src/perspective.h
@@ -2,6 +2,7 @@
 #define PERSPECTIVE_H
 
 #include <vector>
+#include <string>
 
 #include "vector.h"
 #include "matrix.h"
@@ -22,6 +23,12 @@ public:
 	typedef std::vector<Layer> LayerList;
 	
 	
+	static Matrix3 make_matrix(
+		const Vector2 &p0,
+		const Vector2 &px,
+		const Vector2 &py,
+		const Vector2 &p1 );
+	
 	static int truncate_line(
 		Vector2 *out_points,
 		const Pair2 &bounds,
@@ -58,8 +65,12 @@ public:
 	
 	static void add_premulted(
 		const Layer &layer,
-		const RefPtr<Surface> &src_surface,
-		const RefPtr<Surface> &dst_surface );
+		Surface &src_surface,
+		Surface &dst_surface );
+	
+	
+	static void print_layer(const Layer &layer, const std::string &prefix = std::string());
+	static void print_layers(const LayerList &layers, const std::string &prefix = std::string());
 };
 
 
diff --git a/c++/perspective/src/vector.h b/c++/perspective/src/vector.h
index fd19895..9352fae 100644
--- a/c++/perspective/src/vector.h
+++ b/c++/perspective/src/vector.h
@@ -314,11 +314,9 @@ public:
 		{ 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];
+			if (p[i] < p0[i]) p0[i] = p[i];
+			if (p1[i] < p[i]) p1[i] = p[i];
 		}
 		return *this;
 	}