|
|
4dc49c |
|
|
|
3a6f9b |
#include <iostream></iostream>
|
|
|
3a6f9b |
#include <iomanip></iomanip>
|
|
|
3a6f9b |
|
|
|
3a6f9b |
#include "log.h"
|
|
|
3a6f9b |
|
|
|
4dc49c |
#include "perspective.h"
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
4dc49c |
namespace {
|
|
|
3a6f9b |
const int border_width = 2;
|
|
|
3a6f9b |
const int max_width = 16384 - 2*border_width;
|
|
|
3a6f9b |
const int max_height = 16384 - 2*border_width;
|
|
|
4dc49c |
|
|
|
4dc49c |
const int max_area_width = 4096 - 2*border_width;
|
|
|
4dc49c |
const int max_area = max_area_width * max_area_width;
|
|
|
5f2fc2 |
|
|
|
57974f |
const Real max_overscale = 2.0;
|
|
|
57974f |
Real max_overscale_sqr = max_overscale*max_overscale;
|
|
|
43cb04 |
|
|
|
43cb04 |
|
|
|
43cb04 |
class OptimalResolutionSolver {
|
|
|
43cb04 |
private:
|
|
|
43cb04 |
Matrix3 matrix;
|
|
|
43cb04 |
bool affine;
|
|
|
43cb04 |
Vector2 affine_resolution;
|
|
|
43cb04 |
Vector2 focus_a;
|
|
|
43cb04 |
Vector2 focus_b;
|
|
|
43cb04 |
Vector2 focus_m;
|
|
|
43cb04 |
Vector2 fp_kw;
|
|
|
43cb04 |
Vector2 dir;
|
|
|
43cb04 |
Real len;
|
|
|
43cb04 |
|
|
|
43cb04 |
public:
|
|
|
43cb04 |
explicit OptimalResolutionSolver(const Matrix3 &matrix):
|
|
|
43cb04 |
matrix(matrix), affine(), len()
|
|
|
43cb04 |
{
|
|
|
43cb04 |
// w-horizon line equatation is A.z*x + B.z*y + C.z = w
|
|
|
43cb04 |
const Vector3 &A = matrix.row_x();
|
|
|
43cb04 |
const Vector3 &B = matrix.row_y();
|
|
|
43cb04 |
const Vector3 &C = matrix.row_z();
|
|
|
43cb04 |
|
|
|
43cb04 |
const Real wsqr = A.z*A.z + B.z*B.z;
|
|
|
43cb04 |
affine = wsqr <= real_precision_sqr;
|
|
|
43cb04 |
affine_resolution = fabs(C.z) > real_precision
|
|
|
43cb04 |
? Perspective::calc_optimal_resolution(
|
|
|
43cb04 |
matrix.row_x().vec2()/C.z,
|
|
|
43cb04 |
matrix.row_y().vec2()/C.z )
|
|
|
43cb04 |
: Vector2();
|
|
|
43cb04 |
const Real wsqr_div = !affine ? 1/wsqr : Real(0);
|
|
|
43cb04 |
|
|
|
43cb04 |
// focus points
|
|
|
43cb04 |
bool invertible = false;
|
|
|
43cb04 |
const Matrix3 back_matrix = matrix.inverted(&invertible);
|
|
|
43cb04 |
const bool focus_a_exists = invertible && fabs(back_matrix.m02) > real_precision;
|
|
|
43cb04 |
const bool focus_b_exists = invertible && fabs(back_matrix.m12) > real_precision;
|
|
|
43cb04 |
const bool focus_m_exists = focus_a_exists && focus_b_exists;
|
|
|
43cb04 |
assert(focus_a_exists || focus_b_exists);
|
|
|
43cb04 |
focus_a = focus_a_exists ? back_matrix.row_x().vec2()/back_matrix.m02 : Vector2();
|
|
|
43cb04 |
focus_b = focus_b_exists ? back_matrix.row_y().vec2()/back_matrix.m12 : Vector2();
|
|
|
43cb04 |
focus_m = focus_m_exists ? (focus_a + focus_b)*0.5
|
|
|
43cb04 |
: focus_a_exists ? focus_a : focus_b;
|
|
|
43cb04 |
|
|
|
43cb04 |
const Vector2 dist = focus_m_exists ? focus_b - focus_a : Vector2();
|
|
|
43cb04 |
len = dist.length()*0.5;
|
|
|
43cb04 |
dir = fabs(len) > real_precision ? dist/(2*len) : Vector2();
|
|
|
43cb04 |
|
|
|
43cb04 |
// projection of focus points to w-horizon line
|
|
|
43cb04 |
fp_kw = Vector2(A.z, B.z)*wsqr_div;
|
|
|
43cb04 |
}
|
|
|
43cb04 |
|
|
|
43cb04 |
private:
|
|
|
43cb04 |
Real ratio_for_point(const Vector2 &point, Real w) const {
|
|
|
43cb04 |
const Vector3 v = matrix*Vector3(point, 1);
|
|
|
43cb04 |
const Vector2 ox( matrix.m00 - matrix.m02*v.x*w,
|
|
|
43cb04 |
matrix.m10 - matrix.m12*v.x*w );
|
|
|
43cb04 |
const Vector2 oy( matrix.m01 - matrix.m02*v.y*w,
|
|
|
43cb04 |
matrix.m11 - matrix.m12*v.y*w );
|
|
|
43cb04 |
const Real ratio = -ox.length()-oy.length();
|
|
|
43cb04 |
std::cout << Log::to_string(point, 8) << ": " << ratio << std::endl;
|
|
|
43cb04 |
return ratio;
|
|
|
43cb04 |
}
|
|
|
43cb04 |
|
|
|
43cb04 |
Vector2 resolution_for_point(const Vector2 &point, Real w) const {
|
|
|
43cb04 |
const Vector3 v = matrix*Vector3(point, 1);
|
|
|
43cb04 |
const Vector2 ox( (matrix.m00 - matrix.m02*v.x*w)*w,
|
|
|
43cb04 |
(matrix.m01 - matrix.m02*v.y*w)*w );
|
|
|
43cb04 |
const Vector2 oy( (matrix.m10 - matrix.m12*v.x*w)*w,
|
|
|
43cb04 |
(matrix.m11 - matrix.m12*v.y*w)*w );
|
|
|
43cb04 |
return Perspective::calc_optimal_resolution(ox, oy);
|
|
|
43cb04 |
}
|
|
|
43cb04 |
|
|
|
43cb04 |
// returns (l, ratio)
|
|
|
43cb04 |
Vector2 find_max(const Vector2 &point, const Vector2 &dir, Real maxl, Real w) {
|
|
|
43cb04 |
if (maxl <= 1 || maxl >= 1e+10)
|
|
|
43cb04 |
return Vector2(0, ratio_for_point(point, w));
|
|
|
43cb04 |
|
|
|
43cb04 |
Real l0 = 0;
|
|
|
43cb04 |
Real l1 = maxl;
|
|
|
43cb04 |
Real ll0 = (l0 + l1)*0.5;
|
|
|
43cb04 |
std::cout << "begin" << std::endl;
|
|
|
43cb04 |
Real vv0 = ratio_for_point(point + dir*ll0, w);
|
|
|
43cb04 |
while(l1 - l0 > 1) {
|
|
|
43cb04 |
Real ll1, vv1;
|
|
|
43cb04 |
if (ll0 - l0 < l1 - ll0) {
|
|
|
43cb04 |
ll1 = (ll0 + l1)*0.5;
|
|
|
43cb04 |
vv1 = ratio_for_point(point + dir*ll1, w);
|
|
|
43cb04 |
} else {
|
|
|
43cb04 |
ll1 = ll0;
|
|
|
43cb04 |
vv1 = vv0;
|
|
|
43cb04 |
ll0 = (l0 + ll0)*0.5;
|
|
|
43cb04 |
vv0 = ratio_for_point(point + dir*ll0, w);
|
|
|
43cb04 |
}
|
|
|
43cb04 |
|
|
|
43cb04 |
if (vv0 > vv1) {
|
|
|
43cb04 |
l1 = ll1;
|
|
|
43cb04 |
} else {
|
|
|
43cb04 |
l0 = ll0;
|
|
|
43cb04 |
ll0 = ll1;
|
|
|
43cb04 |
vv0 = vv1;
|
|
|
43cb04 |
}
|
|
|
43cb04 |
}
|
|
|
43cb04 |
|
|
|
43cb04 |
std::cout << "end" << std::endl;
|
|
|
43cb04 |
return Vector2(ll0, vv0);
|
|
|
43cb04 |
}
|
|
|
43cb04 |
|
|
|
43cb04 |
public:
|
|
|
43cb04 |
Vector2 solve(Real w, Vector2 *out_center = nullptr) {
|
|
|
43cb04 |
if (out_center) *out_center = Vector2();
|
|
|
43cb04 |
|
|
|
43cb04 |
if (affine)
|
|
|
43cb04 |
return affine_resolution;
|
|
|
43cb04 |
if (w < real_precision)
|
|
|
43cb04 |
return Vector2();
|
|
|
43cb04 |
|
|
|
43cb04 |
Vector2 center;
|
|
|
43cb04 |
const Vector2 offset_w = fp_kw/w;
|
|
|
43cb04 |
if (len <= 1) {
|
|
|
43cb04 |
center = focus_m + offset_w;
|
|
|
43cb04 |
std::cout << "focus_m: " << Log::to_string(focus_m, 8) << std::endl;
|
|
|
43cb04 |
} else {
|
|
|
43cb04 |
const Vector2 solution_a = find_max(focus_a + offset_w, dir, len, w);
|
|
|
43cb04 |
const Vector2 solution_b = find_max(focus_b + offset_w, -dir, len, w);
|
|
|
43cb04 |
center = solution_a.y > solution_b.y
|
|
|
43cb04 |
? focus_a + offset_w + dir*solution_a.x
|
|
|
43cb04 |
: focus_b + offset_w - dir*solution_b.x;
|
|
|
43cb04 |
std::cout << "solution_a: " << Log::to_string(solution_a, 8) << ", len: " << len << std::endl;
|
|
|
43cb04 |
std::cout << "solution_b: " << Log::to_string(solution_b, 8) << std::endl;
|
|
|
43cb04 |
std::cout << "focus_a: " << Log::to_string(focus_a, 8) << std::endl;
|
|
|
43cb04 |
std::cout << "focus_b: " << Log::to_string(focus_b, 8) << std::endl;
|
|
|
43cb04 |
}
|
|
|
43cb04 |
std::cout << "offset_w: " << Log::to_string(offset_w, 8) << std::endl;
|
|
|
43cb04 |
std::cout << "center: " << Log::to_string(center, 8) << std::endl;
|
|
|
43cb04 |
|
|
|
43cb04 |
if (out_center) *out_center = center;
|
|
|
43cb04 |
return resolution_for_point(center, w);
|
|
|
43cb04 |
}
|
|
|
43cb04 |
};
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
3a6f9b |
Matrix3
|
|
|
3a6f9b |
Perspective::make_matrix(
|
|
|
3a6f9b |
const Vector2 &p0,
|
|
|
3a6f9b |
const Vector2 &px,
|
|
|
3a6f9b |
const Vector2 &py,
|
|
|
3a6f9b |
const Vector2 &p1 )
|
|
|
3a6f9b |
{
|
|
|
3a6f9b |
Vector2 A = px - p1;
|
|
|
3a6f9b |
Vector2 B = py - p1;
|
|
|
3a6f9b |
Vector2 C = p0 + p1 - px - py;
|
|
|
3a6f9b |
|
|
|
3a6f9b |
Real cw = A.y*B.x - A.x*B.y;
|
|
|
3a6f9b |
Real aw = B.x*C.y - B.y*C.x;
|
|
|
3a6f9b |
Real bw = A.y*C.x - A.x*C.y;
|
|
|
3a6f9b |
|
|
|
3a6f9b |
// normalize and force cw to be positive
|
|
|
3a6f9b |
Real k = aw*aw + bw*bw + cw*cw;
|
|
|
3a6f9b |
k = k > real_precision ? 1/sqrt(k) : 1;
|
|
|
3a6f9b |
if (cw < 0) k = -k;
|
|
|
3a6f9b |
aw *= k;
|
|
|
3a6f9b |
bw *= k;
|
|
|
3a6f9b |
cw *= k;
|
|
|
3a6f9b |
|
|
|
3a6f9b |
Vector2 c = p0*cw;
|
|
|
3a6f9b |
Vector2 a = px*(cw + aw) - c;
|
|
|
3a6f9b |
Vector2 b = py*(cw + bw) - c;
|
|
|
3a6f9b |
|
|
|
3a6f9b |
return Matrix3( Vector3(a, aw),
|
|
|
3a6f9b |
Vector3(b, bw),
|
|
|
3a6f9b |
Vector3(c, cw) );
|
|
|
3a6f9b |
}
|
|
|
3a6f9b |
|
|
|
4dc49c |
|
|
|
5f2fc2 |
Matrix3
|
|
|
db8a1f |
Perspective::matrix_mult(
|
|
|
db8a1f |
const Matrix3 &matrix,
|
|
|
db8a1f |
Real value )
|
|
|
5f2fc2 |
{
|
|
|
5f2fc2 |
Matrix3 m = matrix;
|
|
|
5f2fc2 |
for(int i = 0; i < 9; ++i)
|
|
|
db8a1f |
m.a[i] *= value;
|
|
|
5f2fc2 |
return m;
|
|
|
5f2fc2 |
}
|
|
|
5f2fc2 |
|
|
|
5f2fc2 |
|
|
|
db8a1f |
Matrix3
|
|
|
db8a1f |
Perspective::normalize_matrix_by_det(
|
|
|
db8a1f |
const Matrix3 &matrix )
|
|
|
db8a1f |
{
|
|
|
db8a1f |
const Real d = matrix.det();
|
|
|
db8a1f |
return fabs(d) <= real_precision ? matrix : matrix_mult(matrix, fabs(1/cbrt(d)));
|
|
|
db8a1f |
}
|
|
|
db8a1f |
|
|
|
db8a1f |
|
|
|
4dc49c |
int
|
|
|
4dc49c |
Perspective::truncate_line(
|
|
|
4dc49c |
Vector2 *out_points,
|
|
|
4dc49c |
const Pair2 &bounds,
|
|
|
4dc49c |
Real a,
|
|
|
4dc49c |
Real b,
|
|
|
4dc49c |
Real c )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
// equatation of line is a*x + b*y + c = 0
|
|
|
4dc49c |
|
|
|
4dc49c |
int count = 0;
|
|
|
4dc49c |
|
|
|
4dc49c |
if (fabs(a) > real_precision) {
|
|
|
4dc49c |
Real x0 = -(c + bounds.p0.y*b)/a;
|
|
|
4dc49c |
if ( x0 + real_precision >= bounds.p0.x
|
|
|
4dc49c |
&& x0 - real_precision <= bounds.p1.x )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
if (out_points) out_points[count] = Vector2(x0, bounds.p0.y);
|
|
|
4dc49c |
if (++count >= 2) return count;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
Real x1 = -(c + bounds.p1.y*b)/a;
|
|
|
4dc49c |
if ( x1 + real_precision >= bounds.p0.x
|
|
|
4dc49c |
&& x1 - real_precision <= bounds.p1.x )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
if (out_points) out_points[count] = Vector2(x1, bounds.p1.y);
|
|
|
4dc49c |
if (++count >= 2) return count;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
if (fabs(b) > real_precision) {
|
|
|
4dc49c |
Real y0 = -(c + bounds.p0.x*a)/b;
|
|
|
4dc49c |
if ( y0 + real_precision >= bounds.p0.y
|
|
|
4dc49c |
&& y0 - real_precision <= bounds.p1.y )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
if (out_points) out_points[count] = Vector2(bounds.p0.x, y0);
|
|
|
4dc49c |
if (++count >= 2) return count;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
Real y1 = -(c + bounds.p1.x*a)/b;
|
|
|
4dc49c |
if ( y1 + real_precision >= bounds.p0.y
|
|
|
4dc49c |
&& y1 - real_precision <= bounds.p1.y )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
if (out_points) out_points[count] = Vector2(bounds.p1.x, y1);
|
|
|
4dc49c |
if (++count >= 2) return count;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
return count;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
4dc49c |
Vector2
|
|
|
4dc49c |
Perspective::calc_optimal_resolution(
|
|
|
4dc49c |
const Vector2 &ox,
|
|
|
4dc49c |
const Vector2 &oy )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
const Real a = ox.x * ox.x;
|
|
|
4dc49c |
const Real b = ox.y * ox.y;
|
|
|
4dc49c |
const Real c = oy.x * oy.x;
|
|
|
4dc49c |
const Real d = oy.y * oy.y;
|
|
|
4dc49c |
Real e = fabs(ox.x*oy.y - ox.y*oy.x);
|
|
|
5f2fc2 |
if (e < real_precision_sqr)
|
|
|
43cb04 |
return Vector2(); // matrix 2x2 ox oy is not invertible
|
|
|
4dc49c |
e = 1.0/e;
|
|
|
4dc49c |
|
|
|
4dc49c |
const Real sum = a*d + b*c;
|
|
|
4dc49c |
Vector2 scale;
|
|
|
db8a1f |
bool abgt = (2*a*b + real_precision_sqr >= sum);
|
|
|
db8a1f |
bool cdgt = (2*c*d + real_precision_sqr >= sum);
|
|
|
db8a1f |
|
|
|
db8a1f |
if (abgt && cdgt) {
|
|
|
db8a1f |
if (a*b < c*d) {
|
|
|
db8a1f |
scale.x = sqrt(2*b)*e;
|
|
|
db8a1f |
scale.y = sqrt(2*a)*e;
|
|
|
db8a1f |
} else {
|
|
|
db8a1f |
scale.x = sqrt(2*d)*e;
|
|
|
db8a1f |
scale.y = sqrt(2*c)*e;
|
|
|
db8a1f |
}
|
|
|
db8a1f |
} else
|
|
|
db8a1f |
if (abgt) {
|
|
|
4dc49c |
scale.x = sqrt(2*b)*e;
|
|
|
4dc49c |
scale.y = sqrt(2*a)*e;
|
|
|
4dc49c |
} else
|
|
|
db8a1f |
if (cdgt) {
|
|
|
4dc49c |
scale.x = sqrt(2*d)*e;
|
|
|
4dc49c |
scale.y = sqrt(2*c)*e;
|
|
|
4dc49c |
} else {
|
|
|
4dc49c |
const Real dif = a*d - b*c;
|
|
|
4dc49c |
scale.x = sqrt(dif/(a - c))*e;
|
|
|
4dc49c |
scale.y = sqrt(dif/(d - b))*e;
|
|
|
4dc49c |
}
|
|
|
5f2fc2 |
|
|
|
4dc49c |
return scale.x <= real_precision || scale.y <= real_precision
|
|
|
4dc49c |
? Vector2() : scale;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
4dc49c |
Vector3
|
|
|
4dc49c |
Perspective::make_alpha_matrix_col(
|
|
|
4dc49c |
Real w0,
|
|
|
4dc49c |
Real w1,
|
|
|
4dc49c |
const Vector3 &w_col )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
Real k = w1 - w0;
|
|
|
4dc49c |
if (fabs(k) <= real_precision)
|
|
|
4dc49c |
return w_col;
|
|
|
4dc49c |
|
|
|
4dc49c |
k = w1/k;
|
|
|
4dc49c |
return Vector3(
|
|
|
4dc49c |
k*w_col.x,
|
|
|
4dc49c |
k*w_col.y,
|
|
|
4dc49c |
k*(w_col.z - w0) );
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
4dc49c |
Matrix3
|
|
|
4dc49c |
Perspective::make_alpha_matrix(
|
|
|
4dc49c |
Real aw0, Real aw1,
|
|
|
4dc49c |
Real bw0, Real bw1,
|
|
|
4dc49c |
const Vector3 &w_col )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
const Vector3 a_col = make_alpha_matrix_col(aw0, aw1, w_col);
|
|
|
4dc49c |
const Vector3 b_col = make_alpha_matrix_col(bw0, bw1, w_col);
|
|
|
4dc49c |
return Matrix3(
|
|
|
4dc49c |
Vector3( a_col.x, b_col.x, w_col.x ),
|
|
|
4dc49c |
Vector3( a_col.y, b_col.y, w_col.y ),
|
|
|
4dc49c |
Vector3( a_col.z, b_col.z, w_col.z ) );
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
4dc49c |
void
|
|
|
4dc49c |
Perspective::calc_raster_size(
|
|
|
4dc49c |
Pair2 &out_bounds,
|
|
|
4dc49c |
IntVector2 &out_raster_size,
|
|
|
4dc49c |
Matrix3 &out_raster_matrix,
|
|
|
4dc49c |
const Vector2 &resolution,
|
|
|
5f2fc2 |
const Pair2 &bounds,
|
|
|
5f2fc2 |
const Vector2 &dst_size )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
const Vector2 offset = bounds.p0;
|
|
|
4dc49c |
const Vector2 raster_size_orig(
|
|
|
4dc49c |
(bounds.p1.x - bounds.p0.x)*resolution.x,
|
|
|
4dc49c |
(bounds.p1.y - bounds.p0.y)*resolution.y );
|
|
|
4dc49c |
|
|
|
4dc49c |
Vector2 raster_size_float = raster_size_orig;
|
|
|
5f2fc2 |
Real sqr = raster_size_float.square();
|
|
|
5f2fc2 |
Real sqr_max = dst_size.square() * max_overscale_sqr;
|
|
|
5f2fc2 |
if (sqr_max > real_precision && sqr > sqr_max)
|
|
|
5f2fc2 |
raster_size_float *= sqrt(sqr_max/sqr);
|
|
|
4dc49c |
if (raster_size_float.x > max_width)
|
|
|
4dc49c |
raster_size_float.x = max_width;
|
|
|
4dc49c |
if (raster_size_float.y > max_height)
|
|
|
4dc49c |
raster_size_float.y = max_height;
|
|
|
5f2fc2 |
|
|
|
4dc49c |
IntVector2 raster_size(
|
|
|
5f2fc2 |
std::max(1, (int)ceil( raster_size_float.x - real_precision )),
|
|
|
5f2fc2 |
std::max(1, (int)ceil( raster_size_float.y - real_precision )) );
|
|
|
4dc49c |
if (raster_size.x * raster_size.y > max_area) {
|
|
|
4dc49c |
const Real k = sqrt(Real(max_area)/(raster_size.x * raster_size.y));
|
|
|
4dc49c |
raster_size.x = std::max(1, (int)floor(raster_size.x*k + real_precision));
|
|
|
4dc49c |
raster_size.y = std::max(1, (int)floor(raster_size.y*k + real_precision));
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
const Vector2 new_resolution(
|
|
|
4dc49c |
raster_size_orig.x > real_precision ? resolution.x * raster_size.x / raster_size_orig.x : resolution.x,
|
|
|
4dc49c |
raster_size_orig.y > real_precision ? resolution.y * raster_size.y / raster_size_orig.y : resolution.y );
|
|
|
3a6f9b |
const Vector2 new_border(
|
|
|
3a6f9b |
new_resolution.x > real_precision ? border_width/new_resolution.x : Real(0),
|
|
|
3a6f9b |
new_resolution.y > real_precision ? border_width/new_resolution.y : Real(0) );
|
|
|
4dc49c |
|
|
|
3a6f9b |
out_bounds = bounds.inflated( new_border );
|
|
|
4dc49c |
out_raster_size = raster_size + IntVector2(border_width, border_width)*2;
|
|
|
4dc49c |
out_raster_matrix = Matrix3::translation(Vector2(border_width, border_width))
|
|
|
4dc49c |
* Matrix3::scaling(new_resolution)
|
|
|
4dc49c |
* Matrix3::translation(-offset);
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
4dc49c |
void
|
|
|
4dc49c |
Perspective::create_layers(
|
|
|
4dc49c |
LayerList &out_layers,
|
|
|
4dc49c |
const Matrix3 &matrix,
|
|
|
4dc49c |
const IntPair2 &dst_bounds,
|
|
|
4dc49c |
const Real step )
|
|
|
4dc49c |
{
|
|
|
4dc49c |
bool is_invertible = false;
|
|
|
db8a1f |
Matrix3 norm_matrix = normalize_matrix_by_det(matrix);
|
|
|
db8a1f |
Matrix3 back_matrix = normalize_matrix_by_det( norm_matrix.inverted(&is_invertible) );
|
|
|
4dc49c |
if (!is_invertible)
|
|
|
4dc49c |
return; // matrix is collapsed
|
|
|
3a6f9b |
|
|
|
4dc49c |
// corners
|
|
|
4dc49c |
Vector3 dst_corners[4] = {
|
|
|
4dc49c |
Vector3(dst_bounds.p0.x, dst_bounds.p0.y, 1),
|
|
|
4dc49c |
Vector3(dst_bounds.p1.x, dst_bounds.p0.y, 1),
|
|
|
4dc49c |
Vector3(dst_bounds.p1.x, dst_bounds.p1.y, 1),
|
|
|
4dc49c |
Vector3(dst_bounds.p0.x, dst_bounds.p1.y, 1) };
|
|
|
4dc49c |
|
|
|
4dc49c |
// calc coefficient for equatation of "horizontal" line: A*x + B*y + C = 1/w (aka A*x + B*y = 1/w - C)
|
|
|
4dc49c |
// equatation of line of horizon is: A*x + B*y + C = 0 (aka A*x + B*y = -C)
|
|
|
4dc49c |
const Real A = back_matrix.m02;
|
|
|
4dc49c |
const Real B = back_matrix.m12;
|
|
|
4dc49c |
const Real C = back_matrix.m22;
|
|
|
4dc49c |
|
|
|
4dc49c |
Real hd = sqrt(A*A + B*B);
|
|
|
4dc49c |
if (hd <= real_precision) {
|
|
|
4dc49c |
// orthogonal projection - no perspective - no subdiviosions
|
|
|
4dc49c |
|
|
|
4dc49c |
if (fabs(C) < real_precision)
|
|
|
4dc49c |
return; // only when matrix was not invertible (additional check)
|
|
|
db8a1f |
back_matrix = matrix_mult(back_matrix, 1/C);
|
|
|
4dc49c |
|
|
|
43cb04 |
// calc src resolution
|
|
|
43cb04 |
const Vector2 resolution = calc_optimal_resolution(
|
|
|
43cb04 |
back_matrix.row_x().vec2(),
|
|
|
43cb04 |
back_matrix.row_y().vec2() );
|
|
|
43cb04 |
if (resolution.x <= real_precision || resolution.y <= real_precision)
|
|
|
43cb04 |
return; // cannot calc resolution, this can happen if matrix is (almost) not invertible
|
|
|
43cb04 |
|
|
|
4dc49c |
// calc bounds
|
|
|
4dc49c |
Pair2 layer_src_bounds(Vector2(INFINITY, INFINITY), Vector2(-INFINITY, -INFINITY));
|
|
|
4dc49c |
for(int i = 0; i < 4; ++i)
|
|
|
3a6f9b |
layer_src_bounds.expand( (back_matrix*dst_corners[i]).vec2() );
|
|
|
4dc49c |
if (layer_src_bounds.empty())
|
|
|
4dc49c |
return;
|
|
|
4dc49c |
|
|
|
3a6f9b |
std::cout << "debug: layer_src_bounds: " << Log::to_string(layer_src_bounds) << std::endl;
|
|
|
3a6f9b |
std::cout << "debug: back_matrix:"<< std::endl << Log::to_string(back_matrix, 10);
|
|
|
3a6f9b |
|
|
|
4dc49c |
// make layer
|
|
|
4dc49c |
Layer layer;
|
|
|
4dc49c |
layer.dst_bounds = dst_bounds;
|
|
|
4dc49c |
calc_raster_size(
|
|
|
4dc49c |
layer.src_bounds,
|
|
|
4dc49c |
layer.src_size,
|
|
|
4dc49c |
layer.back_matrix,
|
|
|
3a6f9b |
resolution,
|
|
|
5f2fc2 |
layer_src_bounds,
|
|
|
5f2fc2 |
Vector2(layer.dst_bounds.size()) );
|
|
|
3a6f9b |
layer.back_matrix *= back_matrix;
|
|
|
4dc49c |
layer.back_alpha_matrix = Matrix3(
|
|
|
4dc49c |
Vector3(0, 0, 0),
|
|
|
4dc49c |
Vector3(0, 0, 0),
|
|
|
4dc49c |
Vector3(1, 1, 1) );
|
|
|
4dc49c |
out_layers.push_back(layer);
|
|
|
4dc49c |
|
|
|
4dc49c |
return;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
// find visible w range
|
|
|
4dc49c |
hd = 1/hd;
|
|
|
4dc49c |
const Real horizonw1 = hd;
|
|
|
4dc49c |
const Real horizonw2 = hd/std::min(Real(2), step);
|
|
|
4dc49c |
const Real horizonw3 = hd/step;
|
|
|
4dc49c |
Real maxw = -INFINITY, minw = INFINITY;
|
|
|
4dc49c |
Vector3 src_corners[4];
|
|
|
4dc49c |
for(int i = 0; i < 4; ++i) {
|
|
|
4dc49c |
Vector3 src = back_matrix * dst_corners[i];
|
|
|
4dc49c |
if (fabs(src.z) > real_precision) {
|
|
|
4dc49c |
Real w = 1/src.z;
|
|
|
4dc49c |
if (w > 0 && w < horizonw1)
|
|
|
4dc49c |
src_corners[i] = Vector3(src.x*w, src.y*w, w);
|
|
|
4dc49c |
else
|
|
|
4dc49c |
w = horizonw1;
|
|
|
4dc49c |
if (minw > w) minw = w;
|
|
|
4dc49c |
if (maxw < w) maxw = w;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
if (minw >= maxw - real_precision)
|
|
|
4dc49c |
return; // all bounds too thin
|
|
|
4dc49c |
const Real maxw3 = std::min(maxw, horizonw3);
|
|
|
4dc49c |
|
|
|
4dc49c |
// steps
|
|
|
4dc49c |
const Real stepLog = log(step);
|
|
|
4dc49c |
int minlog = (int)floor(log(minw)/stepLog + real_precision);
|
|
|
4dc49c |
int maxlog = (int)ceil(log(maxw3)/stepLog - real_precision);
|
|
|
4dc49c |
if (maxlog < minlog) maxlog = minlog;
|
|
|
4dc49c |
|
|
|
43cb04 |
OptimalResolutionSolver resolution_solver(back_matrix);
|
|
|
43cb04 |
|
|
|
4dc49c |
Real w = pow(step, Real(minlog));
|
|
|
4dc49c |
for(int i = minlog; i <= maxlog; ++i, w *= step) {
|
|
|
4dc49c |
// w range
|
|
|
4dc49c |
const Real w0 = w/step;
|
|
|
4dc49c |
const Real w1 = std::min(w*step, horizonw1);
|
|
|
4dc49c |
|
|
|
4dc49c |
// alpha ranges
|
|
|
4dc49c |
const Real aw0 = w0;
|
|
|
4dc49c |
const Real aw1 = w;
|
|
|
4dc49c |
const Real bw0 = i == maxlog ? horizonw1 : w1;
|
|
|
4dc49c |
const Real bw1 = i == maxlog ? horizonw2 : w;
|
|
|
4dc49c |
|
|
|
4dc49c |
// calc bounds
|
|
|
4dc49c |
Pair2 layer_src_bounds(Vector2(INFINITY, INFINITY), Vector2(-INFINITY, -INFINITY));
|
|
|
4dc49c |
Pair2 layer_dst_bounds(Vector2(INFINITY, INFINITY), Vector2(-INFINITY, -INFINITY));
|
|
|
4dc49c |
for(int j = 0; j < 4; ++j) {
|
|
|
4dc49c |
if (src_corners[j].z && src_corners[j].z > w0 && src_corners[j].z < w1) {
|
|
|
4dc49c |
layer_src_bounds.expand(src_corners[j].vec2());
|
|
|
4dc49c |
layer_dst_bounds.expand(dst_corners[j].vec2());
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
Vector2 line[2];
|
|
|
00b5ca |
int line_count = truncate_line(line, Pair2(dst_bounds), A, B, C - 1/w0);
|
|
|
4dc49c |
for(int j = 0; j < line_count; ++j) {
|
|
|
4dc49c |
layer_src_bounds.expand( (back_matrix * Vector3(line[j], 1)).vec2() * w0 );
|
|
|
4dc49c |
layer_dst_bounds.expand( line[j] );
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
00b5ca |
line_count = truncate_line(line, Pair2(dst_bounds), A, B, C - 1/w1);
|
|
|
4dc49c |
for(int j = 0; j < line_count; ++j) {
|
|
|
4dc49c |
layer_src_bounds.expand( (back_matrix * Vector3(line[j], 1)).vec2() * w1 );
|
|
|
4dc49c |
layer_dst_bounds.expand( line[j] );
|
|
|
43cb04 |
|
|
|
43cb04 |
Vector3 v = back_matrix*Vector3(line[j], 1);
|
|
|
43cb04 |
Vector3 vv = norm_matrix*(v*w1);
|
|
|
43cb04 |
assert(real_equal(v.z, 1/w1));
|
|
|
43cb04 |
assert(real_equal(vv.z, w1));
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
if (layer_src_bounds.empty() || layer_dst_bounds.empty())
|
|
|
4dc49c |
continue;
|
|
|
4dc49c |
|
|
|
4dc49c |
// make layer
|
|
|
4dc49c |
Layer layer;
|
|
|
4dc49c |
layer.dst_bounds = IntPair2(
|
|
|
4dc49c |
IntVector2( std::max(dst_bounds.p0.x, (int)floor(layer_dst_bounds.p0.x + real_precision)),
|
|
|
4dc49c |
std::max(dst_bounds.p0.y, (int)floor(layer_dst_bounds.p0.y + real_precision)) ),
|
|
|
4dc49c |
IntVector2( std::min(dst_bounds.p1.x, (int)ceil (layer_dst_bounds.p1.x - real_precision)),
|
|
|
4dc49c |
std::min(dst_bounds.p1.y, (int)ceil (layer_dst_bounds.p1.y - real_precision)) ));
|
|
|
4dc49c |
calc_raster_size(
|
|
|
4dc49c |
layer.src_bounds,
|
|
|
4dc49c |
layer.src_size,
|
|
|
4dc49c |
layer.back_matrix,
|
|
|
43cb04 |
resolution_solver.solve(w, &layer.center),
|
|
|
5f2fc2 |
layer_src_bounds,
|
|
|
5f2fc2 |
Vector2(layer.dst_bounds.size()) );
|
|
|
4dc49c |
layer.back_matrix *= back_matrix;
|
|
|
4dc49c |
layer.back_alpha_matrix = make_alpha_matrix(
|
|
|
00b5ca |
1/aw0, 1/aw1, 1/bw0, 1/bw1,
|
|
|
4dc49c |
layer.back_matrix.get_col(2) );
|
|
|
4dc49c |
out_layers.push_back(layer);
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
4dc49c |
|
|
|
57974f |
Real
|
|
|
57974f |
Perspective::find_optimal_step(
|
|
|
57974f |
const Matrix3 &matrix,
|
|
|
57974f |
const IntPair2 &dst_bounds )
|
|
|
57974f |
{
|
|
|
57974f |
max_overscale_sqr = 1e10;
|
|
|
57974f |
|
|
|
57974f |
ULongInt min_area = -1ull;
|
|
|
57974f |
Real optimal_step = 0;
|
|
|
57974f |
LayerList layers;
|
|
|
57974f |
for(Real step = 2; step < 16; step *= 1.01) {
|
|
|
57974f |
create_layers(layers, matrix, dst_bounds, step);
|
|
|
57974f |
ULongInt area = 0;
|
|
|
57974f |
for(LayerList::const_iterator i = layers.begin(); i != layers.end(); ++i)
|
|
|
57974f |
area += i->src_size.square();
|
|
|
57974f |
if (area < min_area) {
|
|
|
57974f |
min_area = area;
|
|
|
57974f |
optimal_step = step;
|
|
|
57974f |
}
|
|
|
57974f |
layers.clear();
|
|
|
57974f |
}
|
|
|
57974f |
|
|
|
57974f |
max_overscale_sqr = max_overscale*max_overscale;
|
|
|
57974f |
|
|
|
57974f |
return optimal_step;
|
|
|
57974f |
}
|
|
|
57974f |
|
|
|
57974f |
|
|
|
57974f |
|
|
|
4dc49c |
void
|
|
|
4dc49c |
Perspective::add_premulted(
|
|
|
4dc49c |
const Layer &layer,
|
|
|
3a6f9b |
Surface &src_surface,
|
|
|
3a6f9b |
Surface &dst_surface )
|
|
|
4dc49c |
{
|
|
|
3a6f9b |
if (src_surface.empty()) return;
|
|
|
3a6f9b |
if (dst_surface.empty()) return;
|
|
|
4dc49c |
|
|
|
4dc49c |
const int x0 = std::max(layer.dst_bounds.p0.x, 0);
|
|
|
4dc49c |
const int y0 = std::max(layer.dst_bounds.p0.y, 0);
|
|
|
3a6f9b |
const int x1 = std::min(layer.dst_bounds.p1.x, dst_surface.width());
|
|
|
3a6f9b |
const int y1 = std::min(layer.dst_bounds.p1.y, dst_surface.height());
|
|
|
4dc49c |
if (x0 >= x1 || y0 >= y1) return;
|
|
|
4dc49c |
|
|
|
4dc49c |
const int w = x1 - x0;
|
|
|
4dc49c |
const int h = y1 - y0;
|
|
|
4dc49c |
if (w <= 0 || h <= 0) return;
|
|
|
4dc49c |
|
|
|
4dc49c |
const Vector3 coord_dx = layer.back_matrix.row_x();
|
|
|
4dc49c |
const Vector3 coord_dy = layer.back_matrix.row_y() + coord_dx*(x0 - x1);
|
|
|
4dc49c |
Vector3 coord = layer.back_matrix * Vector3(x0, y0, 1);
|
|
|
4dc49c |
|
|
|
4dc49c |
const Vector2 alpha_dx = layer.back_alpha_matrix.row_x().vec2();
|
|
|
4dc49c |
const Vector2 alpha_dy = layer.back_alpha_matrix.row_y().vec2() + alpha_dx*(x0 - x1);
|
|
|
4dc49c |
Vector2 alpha = (layer.back_alpha_matrix * Vector3(x0, y0, 1)).vec2();
|
|
|
4dc49c |
|
|
|
3a6f9b |
const int dr = dst_surface.pitch() - x1 + x0;
|
|
|
3a6f9b |
Color *c = &dst_surface[y0][x0];
|
|
|
4dc49c |
for(int r = y0; r < y1; ++r, c += dr, coord += coord_dy, alpha += alpha_dy) {
|
|
|
4dc49c |
for(Color *end = c + w; c != end; ++c, coord += coord_dx, alpha += alpha_dx) {
|
|
|
4dc49c |
if (coord.z > real_precision) {
|
|
|
4dc49c |
const Real w = 1/coord.z;
|
|
|
4dc49c |
const Real a = clamp(alpha.x*w, 0, 1) * clamp(alpha.y*w, 0, 1);
|
|
|
4dc49c |
if (a > real_precision)
|
|
|
3a6f9b |
*c += src_surface.get_pixel_premulted(coord.x*w, coord.y*w) * a;
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
}
|
|
|
4dc49c |
|
|
|
43cb04 |
void
|
|
|
43cb04 |
Perspective::paint_cross(
|
|
|
43cb04 |
Surface &dst_surface,
|
|
|
43cb04 |
const Vector2 &point )
|
|
|
43cb04 |
{
|
|
|
43cb04 |
const Color black(0, 0, 0, 1);
|
|
|
43cb04 |
const Color white(1, 1, 1, 1);
|
|
|
43cb04 |
|
|
|
43cb04 |
int cx = (int)floor(point.x);
|
|
|
43cb04 |
int cy = (int)floor(point.y);
|
|
|
43cb04 |
for(int i = 0; i < 10; ++i) {
|
|
|
43cb04 |
dst_surface.put_pixel(cx-i, cy, white);
|
|
|
43cb04 |
dst_surface.put_pixel(cx, cy-i, white);
|
|
|
43cb04 |
dst_surface.put_pixel(cx+1+i, cy+1, white);
|
|
|
43cb04 |
dst_surface.put_pixel(cx+1, cy+1+i, white);
|
|
|
43cb04 |
|
|
|
43cb04 |
dst_surface.put_pixel(cx+1+i, cy, black);
|
|
|
43cb04 |
dst_surface.put_pixel(cx+1, cy-i, black);
|
|
|
43cb04 |
dst_surface.put_pixel(cx-i, cy+1, black);
|
|
|
43cb04 |
dst_surface.put_pixel(cx, cy+1+i, black);
|
|
|
43cb04 |
}
|
|
|
43cb04 |
}
|
|
|
43cb04 |
|
|
|
3a6f9b |
|
|
|
3a6f9b |
void
|
|
|
3a6f9b |
Perspective::print_layer(const Layer &layer, const std::string &prefix) {
|
|
|
3a6f9b |
const int w = 10;
|
|
|
3a6f9b |
std::cout << prefix << "dst_bounds: " << Log::to_string(layer.dst_bounds, w) << std::endl;
|
|
|
3a6f9b |
std::cout << prefix << "src_bounds: " << Log::to_string(layer.dst_bounds, w) << std::endl;
|
|
|
3a6f9b |
std::cout << prefix << "src_size: " << Log::to_string(layer.src_size, w) << std::endl;
|
|
|
3a6f9b |
std::cout << prefix << "back_matrix:" << std::endl
|
|
|
3a6f9b |
<< Log::to_string(layer.back_matrix, w, prefix + Log::tab());
|
|
|
3a6f9b |
std::cout << prefix << "back_alpha_matrix:" << std::endl
|
|
|
3a6f9b |
<< Log::to_string(layer.back_alpha_matrix, w, prefix + Log::tab());
|
|
|
43cb04 |
std::cout << prefix << "cetner: " << Log::to_string(layer.center, w) << std::endl;
|
|
|
3a6f9b |
}
|
|
|
3a6f9b |
|
|
|
3a6f9b |
|
|
|
3a6f9b |
void
|
|
|
3a6f9b |
Perspective::print_layers(const LayerList &layers, const std::string &prefix) {
|
|
|
3a6f9b |
int index = 0;
|
|
|
3a6f9b |
for(LayerList::const_iterator i = layers.begin(); i < layers.end(); ++i, ++index) {
|
|
|
3a6f9b |
std::cout << prefix << "layer #" << index << std::endl;
|
|
|
3a6f9b |
print_layer(*i, prefix + Log::tab());
|
|
|
3a6f9b |
}
|
|
|
3a6f9b |
}
|