From fc7a046ad041978bf4071d1e15e6dc2e64176d4c Mon Sep 17 00:00:00 2001 From: Ivan Mahonin Date: Feb 26 2020 09:45:37 +0000 Subject: fix calculation of resolution --- diff --git a/synfig-core/src/modules/lyr_std/warp.cpp b/synfig-core/src/modules/lyr_std/warp.cpp index 1c93696..19c90f9 100644 --- a/synfig-core/src/modules/lyr_std/warp.cpp +++ b/synfig-core/src/modules/lyr_std/warp.cpp @@ -201,6 +201,135 @@ namespace { } + class OptimalResolutionSolver { + private: + Matrix matrix; + bool affine; + Vector affine_resolution; + Vector focus_a; + Vector focus_b; + Vector focus_m; + Vector fp_kw; + Vector dir; + Real len; + + public: + explicit OptimalResolutionSolver(const Matrix &matrix): + affine(), len() + { + this->matrix = matrix; + + // w-horizon line equatation is A[2]*x + B[2]*y + C[2] = w + const Vector3 &A = this->matrix.row_x(); + const Vector3 &B = this->matrix.row_y(); + const Vector3 &C = this->matrix.row_z(); + + const Real wsqr = A[2]*A[2] + B[2]*B[2]; + affine = (wsqr <= real_precision()*real_precision()); + affine_resolution = approximate_zero(C[2]) + ? Vector() + : rendering::TransformationAffine::calc_optimal_resolution(Matrix2( + this->matrix.row_x().to_2d()/C[2], + this->matrix.row_y().to_2d()/C[2] )); + const Real wsqr_div = !affine ? 1/wsqr : Real(0); + + // focus points + bool invertible = this->matrix.is_invertible(); + const Matrix back_matrix = this->matrix.get_inverted(); + const bool focus_a_exists = invertible && approximate_not_zero(back_matrix.m02); + const bool focus_b_exists = invertible && approximate_not_zero(back_matrix.m12); + 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().to_2d()/back_matrix.m02 : Vector(); + focus_b = focus_b_exists ? back_matrix.row_y().to_2d()/back_matrix.m12 : Vector(); + focus_m = focus_m_exists ? (focus_a + focus_b)*0.5 + : focus_a_exists ? focus_a : focus_b; + + const Vector dist = focus_m_exists ? focus_b - focus_a : Vector(); + len = dist.mag()*0.5; + dir = approximate_zero(len) ? Vector() : dist/(2*len); + + // projection of focus points to w-horizon line + fp_kw = Vector(A[2], B[2])*wsqr_div; + } + + private: + Real ratio_for_point(const Vector &point, Real w) const { + const Vector v = matrix.get_transformed(point); + const Vector ox( matrix.m00 - matrix.m02*v[0]*w, + matrix.m10 - matrix.m12*v[0]*w ); + const Vector oy( matrix.m01 - matrix.m02*v[1]*w, + matrix.m11 - matrix.m12*v[1]*w ); + return -ox.mag()-oy.mag(); + } + + Vector resolution_for_point(const Vector &point, Real w) const { + const Vector v = matrix.get_transformed(point); + const Matrix2 m( + Vector( (matrix.m00 - matrix.m02*v[0]*w)*w, + (matrix.m01 - matrix.m02*v[1]*w)*w ), + Vector( (matrix.m10 - matrix.m12*v[0]*w)*w, + (matrix.m11 - matrix.m12*v[1]*w)*w )); + return rendering::TransformationAffine::calc_optimal_resolution(m); + } + + // returns Vector(l, ratio) + Vector find_max(const Vector &point, const Vector &dir, Real maxl, Real w) const { + if (maxl <= 1 || maxl >= 1e+10) + return Vector(0, ratio_for_point(point, w)); + + Real l0 = 0; + Real l1 = maxl; + Real ll0 = (l0 + l1)*0.5; + 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; + } + } + + return Vector(ll0, vv0); + } + + public: + Vector solve(Real w) const { + if (affine) + return affine_resolution; + if (w < real_precision()) + return Vector(); + + Vector center; + const Vector offset_w = fp_kw/w; + if (len <= 1) { + center = focus_m + offset_w; + } else { + const Vector solution_a = find_max(focus_a + offset_w, dir, len, w); + const Vector solution_b = find_max(focus_b + offset_w, -dir, len, w); + center = solution_a[1] > solution_b[1] + ? focus_a + offset_w + dir*solution_a[0] + : focus_b + offset_w - dir*solution_b[0]; + } + + return resolution_for_point(center, w); + } + }; + + class TransformationPerspective: public rendering::Transformation { public: @@ -295,12 +424,14 @@ namespace { if (!bounds.rect.valid()) return Bounds(); - const Matrix norm_matrix = matrix.get_normalized_by_z(); + const Matrix norm_matrix = matrix.get_normalized_by_det(); const Real A = norm_matrix.m02; const Real B = norm_matrix.m12; const Real C = norm_matrix.m22; - if (A*A + B*B <= real_precision() * real_precision()) - return rendering::TransformationAffine::transform_bounds_affine(norm_matrix, bounds); + if (A*A + B*B <= real_precision() * real_precision()) { + const Real C_inv = approximate_zero(C) ? Real(0) : 1/C; + return rendering::TransformationAffine::transform_bounds_affine(norm_matrix*C_inv, bounds); + } const Real horizonw = real_precision(); @@ -343,11 +474,10 @@ namespace { } const Real midw = exp((log(minw) + log(maxw))*0.5); - const Real kwx = 1/(bounds.resolution[0] * midw); - const Real kwy = 1/(bounds.resolution[1] * midw); - out_bounds.resolution = rendering::TransformationAffine::calc_optimal_resolution( - Matrix2( norm_matrix.m00*kwx, norm_matrix.m01*kwx, - norm_matrix.m10*kwy, norm_matrix.m11*kwy )); + const Real krx = 1/bounds.resolution[0]; + const Real kry = 1/bounds.resolution[1]; + const OptimalResolutionSolver resolution_solver( norm_matrix*Matrix().set_scale(krx, kry) ); + out_bounds.resolution = resolution_solver.solve(midw); return out_bounds; } @@ -361,7 +491,7 @@ namespace { if (!bounds.is_valid()) return; - Matrix norm_matrix = matrix.get_normalized_by_z(); + const Matrix norm_matrix = matrix.get_normalized_by_det(); step = std::max(Real(1.1), step); // calc coefficient for equatation of "horizontal" line: A*x + B*y + C = 1/w (aka A*x + B*y = 1/w - C) @@ -381,19 +511,15 @@ namespace { // orthogonal projection, no perspective, no subdiviosions, just make single layer out_layers.push_back(Layer( bounds.rect, - rendering::TransformationAffine::transform_bounds_affine(norm_matrix, bounds), + rendering::TransformationAffine::transform_bounds_affine(norm_matrix*(1/C), bounds), Matrix3( 0, 0, 0, // alpha always one 0, 0, 0, 1, 1, 1 ) )); return; } - // calc base resolution - const Vector resolution = rendering::TransformationAffine::calc_optimal_resolution(Matrix2( - norm_matrix.m00*krx, norm_matrix.m01*krx, - norm_matrix.m10*kry, norm_matrix.m11*kry )); - if (resolution[0] <= real_precision() || resolution[1] <= real_precision()) - return; // cannot calc resolution, seems if matrix is not invertible + // resolution solver + const OptimalResolutionSolver resolution_solver( norm_matrix*Matrix().set_scale(krx, kry) ); // corners Vector3 corners[4] = { @@ -483,10 +609,17 @@ namespace { if (!layer_rect.valid() || !layer_rect_orig.valid()) continue; + // calc resolution + Vector resolution = resolution_solver.solve(w); + if (resolution[0] <= real_precision() || resolution[1] <= real_precision()) + continue; + // make layer out_layers.push_back( Layer( layer_rect_orig, - Bounds(layer_rect, resolution/w), + Bounds( + layer_rect, + resolution ), make_alpha_matrix( 1/aw0, 1/aw1, 1/bw0, 1/bw1, Vector3(norm_matrix.m02, norm_matrix.m12, norm_matrix.m22) ) )); @@ -524,9 +657,8 @@ namespace { const rendering::Transformation::Bounds bounds( source_rect, get_pixels_per_unit().multiply_coords(supersample)); const Matrix back_matrix = transformation->matrix - .get_normalized_by_z() - .get_inverted() - .get_normalized_by_z(); + .get_normalized_by_det() + .get_inverted(); sub_tasks.clear(); layers.clear(); @@ -632,9 +764,8 @@ namespace { const Matrix to_pixels_matrix = from_pixels_matrix.get_inverted(); const Matrix back_matrix = transformation->matrix - .get_normalized_by_z() - .get_inverted() - .get_normalized_by_z(); + .get_normalized_by_det() + .get_inverted(); const Matrix base_matrix = back_matrix * from_pixels_matrix; @@ -783,7 +914,7 @@ Warp::sync() * Matrix().set_scale( 1/src_dist[0], 1/src_dist[1] ) * Matrix().set_translate( -src_tl ); if (matrix.is_invertible()) { - back_matrix = matrix.get_inverted().get_normalized_by_z(); + back_matrix = matrix.get_normalized_by_det().get_inverted(); affine = approximate_zero(matrix.m02) && approximate_zero(matrix.m12); valid = true; } diff --git a/synfig-core/src/synfig/matrix.cpp b/synfig-core/src/synfig/matrix.cpp index a26751c..92d19df 100644 --- a/synfig-core/src/synfig/matrix.cpp +++ b/synfig-core/src/synfig/matrix.cpp @@ -321,16 +321,10 @@ Matrix3::get_inverted()const } Matrix3& -Matrix3::normalize_by_z() +Matrix3::normalize_by_det() { - Real k = m02*m02 + m12*m12 + m22*m22; - if (k > real_precision()*real_precision()) { - k = 1/sqrt(k); - if (m22 < 0) k = -k; - for(int i = 0; i < 3; ++i) - for(int j = 0; j < 3; ++j) - m[i][j] *= k; - } + Real d = det(); + if (approximate_not_zero(d)) *this *= 1/cbrt(fabs(d)); return *this; } diff --git a/synfig-core/src/synfig/matrix.h b/synfig-core/src/synfig/matrix.h index e051cbf..de7e4a6 100644 --- a/synfig-core/src/synfig/matrix.h +++ b/synfig-core/src/synfig/matrix.h @@ -382,13 +382,10 @@ public: Matrix3& invert() { return *this = get_inverted(); } - Matrix3 get_normalized_by_z() const - { return Matrix3(*this).normalize_by_z(); } + Matrix3 get_normalized_by_det() const + { return Matrix3(*this).normalize_by_det(); } - Matrix3& normalize_by_z(); - - Matrix2 to_2d() const - { return Matrix2(m00, m01, m10, m11); } + Matrix3& normalize_by_det(); //!Get the string of the Matrix //!@return String type. A string representation of the matrix diff --git a/synfig-core/src/synfig/rendering/primitive/transformationaffine.cpp b/synfig-core/src/synfig/rendering/primitive/transformationaffine.cpp index 884e5e5..15efbf1 100644 --- a/synfig-core/src/synfig/rendering/primitive/transformationaffine.cpp +++ b/synfig-core/src/synfig/rendering/primitive/transformationaffine.cpp @@ -63,8 +63,8 @@ TransformationAffine::derivative_vfunc(const Point&) const Vector TransformationAffine::calc_optimal_resolution(const Matrix2 &matrix) { const Real max_overscale_sqr = 1.0*4.0; - const Real max_overscale_sqrsqr = max_overscale_sqr*max_overscale_sqr; const Real real_precision_sqr = real_precision() * real_precision(); + const Real real_precision_sqrsqr = real_precision_sqr * real_precision_sqr; const Real a = matrix.m00 * matrix.m00; const Real b = matrix.m01 * matrix.m01; @@ -76,12 +76,19 @@ TransformationAffine::calc_optimal_resolution(const Matrix2 &matrix) { e = 1.0/e; const Real sum = a*d + b*c; + const Real ab2 = 2*a*b + real_precision_sqrsqr; + const Real cd2 = 2*c*d + real_precision_sqrsqr; + bool abgt = (ab2 >= sum); + bool cdgt = (cd2 >= sum); + if (abgt && cdgt) // when both is greater than sum select only lower of them + (ab2 > cd2 ? abgt : cdgt) = false; + Vector scale; - if (2*a*b + real_precision_sqr >= sum) { + if (abgt) { scale[0] = sqrt(2*b)*e; scale[1] = sqrt(2*a)*e; } else - if (2*c*d + real_precision_sqr >= sum) { + if (cdgt) { scale[0] = sqrt(2*d)*e; scale[1] = sqrt(2*c)*e; } else { @@ -90,11 +97,9 @@ TransformationAffine::calc_optimal_resolution(const Matrix2 &matrix) { scale[1] = sqrt(dif/(d - b))*e; } - const Real sx2 = scale[0]*scale[0]; - const Real sy2 = scale[1]*scale[1]; - const Real sqrsqr = (a*sx2 + b*sy2)*(c*sx2 + d*sy2); - if (sqrsqr > max_overscale_sqrsqr) - scale *= sqrt(sqrt(max_overscale_sqrsqr / sqrsqr)); + const Real sqr = scale[0]*scale[1]*(matrix.m00 + matrix.m10)*(matrix.m01 + matrix.m11); + if (sqr > max_overscale_sqr) + scale *= sqrt(max_overscale_sqr/sqr); return scale[0] <= real_precision() || scale[1] <= real_precision() ? Vector() : scale;