From 43cb045bf87eab046cf339584608d842c088e65e Mon Sep 17 00:00:00 2001 From: Ivan Mahonin Date: Feb 25 2020 17:43:00 +0000 Subject: perspective: proper calculation of resolution --- diff --git a/c++/perspective/src/mainwindow.cpp b/c++/perspective/src/mainwindow.cpp index 236c4f7..85632e0 100644 --- a/c++/perspective/src/mainwindow.cpp +++ b/c++/perspective/src/mainwindow.cpp @@ -10,8 +10,8 @@ #include "mainwindow.h" -const bool checks = false; -const bool paint_outside_bounds = false; +const bool checks = true; +const bool paint_outside_bounds = true; MainWindow::MainWindow(): @@ -192,8 +192,8 @@ MainWindow::update_dst_surface() { Perspective::print_layers(layers, Log::tab(2)); std::cout << Log::tab() << "---------------" << std::endl; - std::cout << Log::tab() << "optimal step: " - << Perspective::find_optimal_step(matrix, dst_bounds_pixels) << std::endl; + //std::cout << Log::tab() << "optimal step: " + // << Perspective::find_optimal_step(matrix, dst_bounds_pixels) << std::endl; // override dst_bounds if (paint_outside_bounds) @@ -209,6 +209,8 @@ MainWindow::update_dst_surface() { layer_surface.mult_alpha(); Perspective::add_premulted(*i, layer_surface, surface); } + for(Perspective::LayerList::const_iterator i = layers.begin(); i != layers.end(); ++i) + Perspective::paint_cross(surface, i->center); dst_surface = surface.to_cairo_surface(true); std::cout << "update_dst_surface() - end" << std::endl; diff --git a/c++/perspective/src/perspective.cpp b/c++/perspective/src/perspective.cpp index 0d59e20..a2b1ecc 100644 --- a/c++/perspective/src/perspective.cpp +++ b/c++/perspective/src/perspective.cpp @@ -17,6 +17,146 @@ namespace { const Real max_overscale = 2.0; Real max_overscale_sqr = max_overscale*max_overscale; + + + class OptimalResolutionSolver { + private: + Matrix3 matrix; + bool affine; + Vector2 affine_resolution; + Vector2 focus_a; + Vector2 focus_b; + Vector2 focus_m; + Vector2 fp_kw; + Vector2 dir; + Real len; + + public: + explicit OptimalResolutionSolver(const Matrix3 &matrix): + matrix(matrix), affine(), len() + { + // w-horizon line equatation is A.z*x + B.z*y + C.z = w + const Vector3 &A = matrix.row_x(); + const Vector3 &B = matrix.row_y(); + const Vector3 &C = matrix.row_z(); + + const Real wsqr = A.z*A.z + B.z*B.z; + affine = wsqr <= real_precision_sqr; + affine_resolution = fabs(C.z) > real_precision + ? Perspective::calc_optimal_resolution( + matrix.row_x().vec2()/C.z, + matrix.row_y().vec2()/C.z ) + : Vector2(); + const Real wsqr_div = !affine ? 1/wsqr : Real(0); + + // focus points + bool invertible = false; + const Matrix3 back_matrix = matrix.inverted(&invertible); + const bool focus_a_exists = invertible && fabs(back_matrix.m02) > real_precision; + const bool focus_b_exists = invertible && fabs(back_matrix.m12) > real_precision; + const bool focus_m_exists = focus_a_exists && focus_b_exists; + assert(focus_a_exists || focus_b_exists); + focus_a = focus_a_exists ? back_matrix.row_x().vec2()/back_matrix.m02 : Vector2(); + focus_b = focus_b_exists ? back_matrix.row_y().vec2()/back_matrix.m12 : Vector2(); + focus_m = focus_m_exists ? (focus_a + focus_b)*0.5 + : focus_a_exists ? focus_a : focus_b; + + const Vector2 dist = focus_m_exists ? focus_b - focus_a : Vector2(); + len = dist.length()*0.5; + dir = fabs(len) > real_precision ? dist/(2*len) : Vector2(); + + // projection of focus points to w-horizon line + fp_kw = Vector2(A.z, B.z)*wsqr_div; + } + + private: + Real ratio_for_point(const Vector2 &point, Real w) const { + const Vector3 v = matrix*Vector3(point, 1); + const Vector2 ox( matrix.m00 - matrix.m02*v.x*w, + matrix.m10 - matrix.m12*v.x*w ); + const Vector2 oy( matrix.m01 - matrix.m02*v.y*w, + matrix.m11 - matrix.m12*v.y*w ); + const Real ratio = -ox.length()-oy.length(); + std::cout << Log::to_string(point, 8) << ": " << ratio << std::endl; + return ratio; + } + + Vector2 resolution_for_point(const Vector2 &point, Real w) const { + const Vector3 v = matrix*Vector3(point, 1); + const Vector2 ox( (matrix.m00 - matrix.m02*v.x*w)*w, + (matrix.m01 - matrix.m02*v.y*w)*w ); + const Vector2 oy( (matrix.m10 - matrix.m12*v.x*w)*w, + (matrix.m11 - matrix.m12*v.y*w)*w ); + return Perspective::calc_optimal_resolution(ox, oy); + } + + // returns (l, ratio) + Vector2 find_max(const Vector2 &point, const Vector2 &dir, Real maxl, Real w) { + if (maxl <= 1 || maxl >= 1e+10) + return Vector2(0, ratio_for_point(point, w)); + + Real l0 = 0; + Real l1 = maxl; + Real ll0 = (l0 + l1)*0.5; + std::cout << "begin" << std::endl; + Real vv0 = ratio_for_point(point + dir*ll0, w); + while(l1 - l0 > 1) { + Real ll1, vv1; + if (ll0 - l0 < l1 - ll0) { + ll1 = (ll0 + l1)*0.5; + vv1 = ratio_for_point(point + dir*ll1, w); + } else { + ll1 = ll0; + vv1 = vv0; + ll0 = (l0 + ll0)*0.5; + vv0 = ratio_for_point(point + dir*ll0, w); + } + + if (vv0 > vv1) { + l1 = ll1; + } else { + l0 = ll0; + ll0 = ll1; + vv0 = vv1; + } + } + + std::cout << "end" << std::endl; + return Vector2(ll0, vv0); + } + + public: + Vector2 solve(Real w, Vector2 *out_center = nullptr) { + if (out_center) *out_center = Vector2(); + + if (affine) + return affine_resolution; + if (w < real_precision) + return Vector2(); + + Vector2 center; + const Vector2 offset_w = fp_kw/w; + if (len <= 1) { + center = focus_m + offset_w; + std::cout << "focus_m: " << Log::to_string(focus_m, 8) << std::endl; + } else { + const Vector2 solution_a = find_max(focus_a + offset_w, dir, len, w); + const Vector2 solution_b = find_max(focus_b + offset_w, -dir, len, w); + center = solution_a.y > solution_b.y + ? focus_a + offset_w + dir*solution_a.x + : focus_b + offset_w - dir*solution_b.x; + std::cout << "solution_a: " << Log::to_string(solution_a, 8) << ", len: " << len << std::endl; + std::cout << "solution_b: " << Log::to_string(solution_b, 8) << std::endl; + std::cout << "focus_a: " << Log::to_string(focus_a, 8) << std::endl; + std::cout << "focus_b: " << Log::to_string(focus_b, 8) << std::endl; + } + std::cout << "offset_w: " << Log::to_string(offset_w, 8) << std::endl; + std::cout << "center: " << Log::to_string(center, 8) << std::endl; + + if (out_center) *out_center = center; + return resolution_for_point(center, w); + } + }; } @@ -54,7 +194,7 @@ Perspective::make_matrix( Matrix3 -Perspective::normalize_matrix( +Perspective::normalize_matrix_by_w( const Matrix3 &matrix ) { Real k = matrix.get_col(2).length(); @@ -131,7 +271,7 @@ Perspective::calc_optimal_resolution( const Real d = oy.y * oy.y; Real e = fabs(ox.x*oy.y - ox.y*oy.x); if (e < real_precision_sqr) - return Vector2(); // paranoid check, e must be non-zero when matrix is invertible + return Vector2(); // matrix 2x2 ox oy is not invertible e = 1.0/e; const Real sum = a*d + b*c; @@ -243,17 +383,11 @@ Perspective::create_layers( const Real step ) { bool is_invertible = false; - Matrix3 back_matrix = normalize_matrix(matrix).inverted(&is_invertible); + Matrix3 norm_matrix = normalize_matrix_by_w(matrix); + Matrix3 back_matrix = norm_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; // cannot calc resolution, this can happen if matrix is (almost) not invertible - // corners Vector3 dst_corners[4] = { Vector3(dst_bounds.p0.x, dst_bounds.p0.y, 1), @@ -274,6 +408,13 @@ Perspective::create_layers( if (fabs(C) < real_precision) return; // only when matrix was not invertible (additional check) + // 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; // cannot calc resolution, this can happen if matrix is (almost) not invertible + // calc bounds Pair2 layer_src_bounds(Vector2(INFINITY, INFINITY), Vector2(-INFINITY, -INFINITY)); for(int i = 0; i < 4; ++i) @@ -333,6 +474,8 @@ Perspective::create_layers( int maxlog = (int)ceil(log(maxw3)/stepLog - real_precision); if (maxlog < minlog) maxlog = minlog; + OptimalResolutionSolver resolution_solver(back_matrix); + Real w = pow(step, Real(minlog)); for(int i = minlog; i <= maxlog; ++i, w *= step) { // w range @@ -365,6 +508,11 @@ Perspective::create_layers( 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] ); + + Vector3 v = back_matrix*Vector3(line[j], 1); + Vector3 vv = norm_matrix*(v*w1); + assert(real_equal(v.z, 1/w1)); + assert(real_equal(vv.z, w1)); } if (layer_src_bounds.empty() || layer_dst_bounds.empty()) @@ -381,7 +529,7 @@ Perspective::create_layers( layer.src_bounds, layer.src_size, layer.back_matrix, - resolution / w, + resolution_solver.solve(w, &layer.center), layer_src_bounds, Vector2(layer.dst_bounds.size()) ); layer.back_matrix *= back_matrix; @@ -463,6 +611,29 @@ Perspective::add_premulted( } } +void +Perspective::paint_cross( + Surface &dst_surface, + const Vector2 &point ) +{ + const Color black(0, 0, 0, 1); + const Color white(1, 1, 1, 1); + + int cx = (int)floor(point.x); + int cy = (int)floor(point.y); + for(int i = 0; i < 10; ++i) { + dst_surface.put_pixel(cx-i, cy, white); + dst_surface.put_pixel(cx, cy-i, white); + dst_surface.put_pixel(cx+1+i, cy+1, white); + dst_surface.put_pixel(cx+1, cy+1+i, white); + + dst_surface.put_pixel(cx+1+i, cy, black); + dst_surface.put_pixel(cx+1, cy-i, black); + dst_surface.put_pixel(cx-i, cy+1, black); + dst_surface.put_pixel(cx, cy+1+i, black); + } +} + void Perspective::print_layer(const Layer &layer, const std::string &prefix) { @@ -474,6 +645,7 @@ Perspective::print_layer(const Layer &layer, const std::string &prefix) { << 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()); + std::cout << prefix << "cetner: " << Log::to_string(layer.center, w) << std::endl; } diff --git a/c++/perspective/src/perspective.h b/c++/perspective/src/perspective.h index 6d75a62..d9eced0 100644 --- a/c++/perspective/src/perspective.h +++ b/c++/perspective/src/perspective.h @@ -18,6 +18,7 @@ public: IntVector2 src_size; Matrix3 back_matrix; Matrix3 back_alpha_matrix; + Vector2 center; }; typedef std::vector LayerList; @@ -29,9 +30,9 @@ public: const Vector2 &py, const Vector2 &p1 ); - static Matrix3 normalize_matrix( + static Matrix3 normalize_matrix_by_w( const Matrix3 &matrix ); - + static int truncate_line( Vector2 *out_points, const Pair2 &bounds, @@ -75,6 +76,10 @@ public: const Layer &layer, Surface &src_surface, Surface &dst_surface ); + + static void paint_cross( + Surface &dst_surface, + const Vector2 &point ); static void print_layer(const Layer &layer, const std::string &prefix = std::string()); diff --git a/c++/perspective/src/perspectivesandboxview.cpp b/c++/perspective/src/perspectivesandboxview.cpp index 650b403..4687875 100644 --- a/c++/perspective/src/perspectivesandboxview.cpp +++ b/c++/perspective/src/perspectivesandboxview.cpp @@ -36,13 +36,22 @@ static Vector2 calc_optimal_resolution( return Vector2(); // paranoid check, e must be non-zero when matrix is invertible e = 1.0/e; - const Real sum = a*d + b*c; + const Real sum = a*d + b*c - real_precision_sqr; Vector2 scale; - if (2*a*b + real_precision_sqr >= sum) { + if (2*a*b >= sum && 2*c*d >= sum) { + if (a*b < c*d) { + scale.x = sqrt(2*b)*e; + scale.y = sqrt(2*a)*e; + } else { + scale.x = sqrt(2*d)*e; + scale.y = sqrt(2*c)*e; + } + } else + if (2*a*b >= sum) { scale.x = sqrt(2*b)*e; scale.y = sqrt(2*a)*e; } else - if (2*c*d + real_precision_sqr >= sum) { + if (2*c*d >= sum) { scale.x = sqrt(2*d)*e; scale.y = sqrt(2*c)*e; } else { @@ -89,6 +98,92 @@ void PerspectiveSandBoxView::update_background( background = surface.to_cairo_surface(); } + +void PerspectiveSandBoxView::update_background_best_perimeter( + const Matrix3 &matrix, + int width, + int height ) +{ + const int log_min = -10000; + const int log_max = 10000; + const Real log_base = log(1.001); + + std::vector map_pk_max(log_max - log_min + 2, 0); + std::vector map_pk_min(log_max - log_min + 2, 1); + std::vector map_pk(width*height, 0); + std::vector map_log_index(width*height, 0); + + + // fill map + for(int r = 0; r < height; ++r) { + for(int c = 0; c < width; ++c) { + Vector3 v = matrix*Vector3(c, r, 1); + if (v.z <= real_precision) + continue; + + Real l = round(log(v.z)/log_base); + if (l-1 < log_min || l+1 > log_max) + continue; + int log_index = (int)l - log_min + 1; + assert(log_index > 0 && log_index < (int)map_pk_max.size() && log_index < (int)map_pk_min.size()); + + const Real k = 1/(v.z*v.z); + const Vector2 ox = Vector2( v.z*matrix.m00 - matrix.m02*v.x, + v.z*matrix.m10 - matrix.m12*v.x )*k; + const Vector2 oy = Vector2( v.z*matrix.m01 - matrix.m02*v.y, + v.z*matrix.m11 - matrix.m12*v.y )*k; + + const Real area = fabs(ox.x*oy.y - ox.y*oy.x); + const Real optimal_side = sqrt(area); + const Real optimal_perimeter = 4*optimal_side; + const Real perimeter = (ox.length() + oy.length())*2; + + const Real pk = optimal_perimeter/perimeter; + + map_pk[r*width + c] = pk; + map_log_index[r*width + c] = log_index; + if (map_pk_max[log_index] < pk) + map_pk_max[log_index] = pk; + if (map_pk_min[log_index] > pk) + map_pk_min[log_index] = pk; + } + } + + // paint + DataSurface surface(width, height); + for(int r = 0; r < height; ++r) { + for(int c = 0; c < width; ++c) { + const int log_index = map_log_index[r*width + c]; + assert(log_index >= 0 && log_index < (int)map_pk_max.size() && log_index < (int)map_pk_min.size()); + if (log_index) { + const Real v = map_pk[r*width + c]; + const Real v0 = map_pk_min[log_index]; + const Real v1 = map_pk_max[log_index]; + if (v0 + real_precision < v1) { + Real k = (v - v0)/(v1 - v0); + const Real p = 0.5; + k -= p; + if (k > 0) { + k /= 1 - p; + k = pow(k, 20); + surface[r][c] = Color(k, k, k, 1); + } + } + } + } + } + + // paint centers + Perspective::LayerList layers; + Perspective::create_layers(layers, matrix.inverted(), IntPair2(IntVector2(), IntVector2(width, height)), 2); + for(Perspective::LayerList::const_iterator i = layers.begin(); i != layers.end(); ++i) + Perspective::paint_cross(surface, i->center); + std::cout << "layers count: " << layers.size() << std::endl; + + background = surface.to_cairo_surface(); +} + + void PerspectiveSandBoxView::draw_grid( const Cairo::RefPtr &context, @@ -284,8 +379,8 @@ PerspectiveSandBoxView::on_draw_view(const Cairo::RefPtr &contex Vector3( m4.m00, m4.m01, m4.m03 ), Vector3( m4.m10, m4.m11, m4.m13 ), Vector3( m4.m30, m4.m31, m4.m33 ) ); - Matrix3 back_matrix = Perspective::normalize_matrix( m3.inverted() ); - update_background(back_matrix, get_width(), get_height()); + Matrix3 back_matrix = Perspective::normalize_matrix_by_w( m3.inverted() ); + update_background_best_perimeter(back_matrix, get_width(), get_height()); } if (background) { diff --git a/c++/perspective/src/perspectivesandboxview.h b/c++/perspective/src/perspectivesandboxview.h index 115923b..6c9a3df 100644 --- a/c++/perspective/src/perspectivesandboxview.h +++ b/c++/perspective/src/perspectivesandboxview.h @@ -23,6 +23,11 @@ public: int width, int height ); + void update_background_best_perimeter( + const Matrix3 &matrix, + int width, + int height ); + void draw_grid( const Cairo::RefPtr &context, const Matrix &matrix, diff --git a/c++/perspective/src/surface.h b/c++/perspective/src/surface.h index e1f36db..5d734da 100644 --- a/c++/perspective/src/surface.h +++ b/c++/perspective/src/surface.h @@ -53,6 +53,11 @@ public: inline Color* operator[] (int index) { return row(index); } inline const Color* operator[] (int index) const { return row(index); } + inline void put_pixel(int x, int y, const Color &color) { + if (x >= 0 && y >= 0 && x < width() && y < height()) + row(y)[x] = color; + } + void mult_alpha(); void demult_alpha(); Color get_pixel(Real x, Real y) const;