ale/d2/trans_single.h
2022-07-30 14:46:04 -03:00

841 lines
21 KiB
C++

// Copyright 2002, 2004 David Hilvert <dhilvert@auricle.dyndns.org>,
// <dhilvert@ugcs.caltech.edu>
/* This file is part of the Anti-Lamenessing Engine.
The Anti-Lamenessing Engine is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with the Anti-Lamenessing Engine; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
* trans_single.h: Represent transformations of the kind q = c(b^-1(p)),
* where p is a point in the source coordinate system, q is a point in the
* target coordinate system, b^-1 is a transformation correcting barrel
* distortion, and c is a transformation of projective or Euclidean type.
* (Note that ^-1 in this context indicates the function inverse rather than
* the exponential.)
*/
#ifndef __trans_single_h__
#define __trans_single_h__
#include "trans_abstract.h"
/*
* transformation: a structure to describe a transformation of kind q =
* c(b^-1(p)), where p is a point in the source coordinate system, q is a point
* in the target coordinate system, b^-1 is a transformation correcting barrel
* distortion, and c is a projective or Euclidean transformation. (Note that
* ^-1 in this case indicates a function inverse, not exponentiation.) Data
* elements are divided into those describing barrel distortion correction and
* those describing projective/Euclidean transformations.
*
* Barrel distortion correction estimates barrel distortion using polynomial
* functions of distance from the center of an image, following (roughly) the
* example set by Helmut Dersch in his PanoTools software:
*
* http://www.path.unimelb.edu.au/~dersch/barrel/barrel.html
*
* Projective transformation data member names roughly correspond to a typical
* treatment of projective transformations from:
*
* Heckbert, Paul. "Projective Mappings for Image Warping." Excerpted
* from his Master's Thesis (UC Berkeley, 1989). 1995.
*
* http://www.cs.cmu.edu/afs/cs/project/classes-ph/862.95/www/notes/proj.ps
*
* For convenience, Heckbert's 'x' and 'y' are noted here numerically by '0'
* and '1', respectively. 'x0' is denoted 'x[0][0]'; 'y0' is 'x[0][1]'.
*
* eu[i] are the parameters for euclidean transformations.
*
* We consider points to be transformed as homogeneous coordinate vectors
* multiplied on the right of the transformation matrix, and so we consider the
* transformation matrix as
*
* - -
* | a b c |
* | d e f |
* | g h i |
* - -
*
* where element i is always equal to 1.
*
*/
struct trans_single : public trans_abstract {
private:
point x[4];
ale_pos eu[3];
mutable ale_pos a, b, c, d, e, f, g, h; // matrix
mutable ale_pos _a, _b, _c, _d, _e, _f, _g, _h; // matrix inverse
int _is_projective;
pixel tonal_multiplier;
mutable int resultant_memo;
mutable int resultant_inverse_memo;
/*
* Calculate resultant matrix values.
*/
void resultant() const {
/*
* If we already know the answers, don't bother calculating
* them again.
*/
if (resultant_memo)
return;
int ale_pos_casting = ale_pos_casting_status();
ale_pos_enable_casting();
if (_is_projective) {
/*
* Calculate resultant matrix values for a general
* projective transformation given that we are mapping
* from the source domain of dimension input_height *
* input_width to a specified arbitrary quadrilateral.
* Follow the calculations outlined in the document by
* Paul Heckbert cited above for the case in which the
* source domain is a unit square and then divide to
* correct for the scale factor in each dimension.
*/
/*
* First, perform calculations as outlined in Heckbert.
*/
ale_pos delta_01 = x[1][0] - x[2][0];
ale_pos delta_02 = x[3][0] - x[2][0];
ale_pos sigma_0 = x[0][0] - x[1][0] + x[2][0] - x[3][0];
ale_pos delta_11 = x[1][1] - x[2][1];
ale_pos delta_12 = x[3][1] - x[2][1];
ale_pos sigma_1 = x[0][1] - x[1][1] + x[2][1] - x[3][1];
g = (sigma_0 * delta_12 - sigma_1 * delta_02)
/ (delta_01 * delta_12 - delta_11 * delta_02);
h = (delta_01 * sigma_1 - delta_11 * sigma_0 )
/ (delta_01 * delta_12 - delta_11 * delta_02);
a = (x[1][0] - x[0][0] + g * x[1][0]);
b = (x[3][0] - x[0][0] + h * x[3][0]);
c = x[0][0];
d = (x[1][1] - x[0][1] + g * x[1][1]);
e = (x[3][1] - x[0][1] + h * x[3][1]);
f = x[0][1];
/*
* Finish by scaling so that our transformation maps
* from a rectangle of width and height matching the
* width and height of the input image.
*/
a /= input_height;
b /= input_width;
d /= input_height;
e /= input_width;
g /= input_height;
h /= input_width;
} else {
/*
* Calculate matrix values for a euclidean
* transformation.
*
* We want to translate the image center by (eu[0],
* eu[1]) and rotate the image about the center by
* eu[2] degrees. This is equivalent to the following
* sequence of affine transformations applied to the
* point to be transformed:
*
* translate by (-h/2, -w/2)
* rotate by eu[2] degrees about the origin
* translate by (h/2, w/2)
* translate by (eu[0], eu[1])
*
* The matrix assigned below represents the result of
* combining all of these transformations. Matrix
* elements g and h are always zero in an affine
* transformation.
*/
ale_pos theta = (double) eu[2] * M_PI / 180;
a = (double) cos(theta) * (double) scale_factor;
b = (double) sin(theta) * (double) scale_factor;
c = 0.5 * ((double) input_height * ((double) scale_factor - (double) a)
- (double) input_width * (double) b) + (double) eu[0]
* (double) scale_factor;
d = -b;
e = a;
f = 0.5 * ((double) input_height * (double) b
+ (double) input_width * ((double) scale_factor
- (double) a))
+ (double) eu[1] * (double) scale_factor;
g = 0;
h = 0;
}
resultant_memo = 1;
if (!ale_pos_casting)
ale_pos_disable_casting();
}
/*
* Calculate the inverse transform matrix values.
*/
void resultant_inverse () const {
/*
* If we already know the answers, don't bother calculating
* them again.
*/
if (resultant_inverse_memo)
return;
resultant();
int ale_pos_casting = ale_pos_casting_status();
ale_pos_enable_casting();
/*
* For projective transformations, we calculate
* the inverse of the forward transformation
* matrix.
*/
double scale = (double) a * (double) e - (double) b * (double) d;
_a = ((double) e * 1 - (double) f * (double) h) / scale;
_b = ((double) h * (double) c - 1 * (double) b) / scale;
_c = ((double) b * (double) f - (double) c * (double) e) / scale;
_d = ((double) f * (double) g - (double) d * 1) / scale;
_e = (1 * (double) a - (double) g * (double) c) / scale;
_f = ((double) c * (double) d - (double) a * (double) f) / scale;
_g = ((double) d * (double) h - (double) e * (double) g) / scale;
_h = ((double) g * (double) b - (double) h * (double) a) / scale;
resultant_inverse_memo = 1;
if (!ale_pos_casting)
ale_pos_disable_casting();
}
public:
trans_single &operator=(const trans_single &ta) {
this->trans_abstract::operator=(*((trans_abstract *) &ta));
for (int i = 0; i < 4; i++) {
x[i] = ta.x[i];
}
for (int i = 0; i < 3; i++) {
eu[i] = ta.eu[i];
}
_is_projective = ta._is_projective;
tonal_multiplier = ta.tonal_multiplier;
resultant_memo = 0;
resultant_inverse_memo = 0;
return *this;
}
trans_single(const trans_single &ta) {
operator=(ta);
}
trans_single() {
}
/*
* Returns non-zero if the transformation might be non-Euclidean.
*/
int is_projective() const {
return _is_projective;
}
/*
* Projective/Euclidean transformation
*/
struct point pe(struct point p) const {
struct point result;
resultant();
result[0] = (a * p[0] + b * p[1] + c)
/ (g * p[0] + h * p[1] + 1);
result[1] = (d * p[0] + e * p[1] + f)
/ (g * p[0] + h * p[1] + 1);
return result;
}
/*
* Projective/Euclidean inverse
*/
struct point pei(struct point p) const {
struct point result;
resultant_inverse();
result[0] = (_a * p[0] + _b * p[1] + _c)
/ (_g * p[0] + _h * p[1] + 1);
result[1] = (_d * p[0] + _e * p[1] + _f)
/ (_g * p[0] + _h * p[1] + 1);
return result;
}
#if 0
/*
* operator() is the transformation operator.
*/
struct point operator()(struct point p) {
return transform(p);
}
#endif
/*
* Calculate projective transformation parameters from a euclidean
* transformation.
*/
void eu_to_gpt() {
assert(!_is_projective);
x[0] = transform_unscaled(point( 0 , 0 ) );
x[1] = transform_unscaled(point( input_height, 0 ) );
x[2] = transform_unscaled(point( input_height, input_width ) );
x[3] = transform_unscaled(point( 0 , input_width ) );
resultant_memo = 0;
resultant_inverse_memo = 0;
_is_projective = 1;
}
/*
* Calculate euclidean identity transform for a given image.
*/
static struct trans_single eu_identity(const image *i = NULL, ale_pos scale_factor = 1) {
struct trans_single r;
r.resultant_memo = 0;
r.resultant_inverse_memo = 0;
r.eu[0] = 0;
r.eu[1] = 0;
r.eu[2] = 0;
r.input_width = i ? i->width() : 2;
r.input_height = i ? i->height() : 2;
r.scale_factor = scale_factor;
r._is_projective = 0;
r.tonal_multiplier = pixel(1, 1, 1);
r.bd_set(0, (ale_pos *) NULL);
return r;
}
/*
* Calculate projective identity transform for a given image.
*/
static trans_single gpt_identity(const image *i, ale_pos scale_factor) {
struct trans_single r = eu_identity(i, scale_factor);
r.eu_to_gpt();
return r;
}
/*
* Set the tonal multiplier
*/
void set_tonal_multiplier(pixel p) {
tonal_multiplier = p;
}
pixel get_tonal_multiplier(struct point p) const {
return tonal_multiplier;
}
pixel get_inverse_tonal_multiplier(struct point p) const {
return tonal_multiplier;
}
/*
* Modify a euclidean transform in the indicated manner.
*/
void eu_modify(int i1, ale_pos diff) {
assert(!_is_projective);
resultant_memo = 0;
resultant_inverse_memo = 0;
if (i1 < 2)
eu[i1] += diff / scale_factor;
else
eu[i1] += diff;
}
/*
* Rotate about a given point in the original reference frame.
*/
void eu_rotate_about_scaled(point center, ale_pos diff) {
assert(center.defined());
point fixpoint = scaled_inverse_transform(center);
eu_modify(2, diff);
point offset = center - transform_scaled(fixpoint);
eu_modify(0, offset[0]);
eu_modify(1, offset[1]);
}
/*
* Modify all euclidean parameters at once.
*/
void eu_set(ale_pos eu[3]) {
resultant_memo = 0;
resultant_inverse_memo = 0;
this->eu[0] = eu[0] / scale_factor;
this->eu[1] = eu[1] / scale_factor;
this->eu[2] = eu[2];
if (_is_projective) {
_is_projective = 0;
eu_to_gpt();
}
}
/*
* Get the specified euclidean parameter
*/
ale_pos eu_get(int param) const {
assert (!_is_projective);
assert (param >= 0);
assert (param < 3);
if (param < 2)
return eu[param] * scale_factor;
else
return eu[param];
}
/*
* Modify a projective transform in the indicated manner.
*/
void gpt_modify(int i1, int i2, ale_pos diff) {
assert (_is_projective);
resultant_memo = 0;
resultant_inverse_memo = 0;
x[i2][i1] += diff;
}
/*
* Modify a projective transform according to the group operation.
*/
void gr_modify(int i1, int i2, ale_pos diff) {
assert (_is_projective);
assert (i1 == 0 || i1 == 1);
point diff_vector = (i1 == 0)
? point(diff, 0)
: point(0, diff);
trans_single t = *this;
t.resultant_memo = 0;
t.resultant_inverse_memo = 0;
t.input_height = (unsigned int) scaled_height();
t.input_width = (unsigned int) scaled_width();
t.scale_factor = 1;
t.bd_set(0, (ale_pos *) NULL);
resultant_memo = 0;
resultant_inverse_memo = 0;
x[i2] = t.transform_scaled(t.scaled_inverse_transform(x[i2]) + diff_vector);
}
/*
* Modify all projective parameters at once.
*/
void gpt_set(point x[4]) {
resultant_memo = 0;
resultant_inverse_memo = 0;
_is_projective = 1;
for (int i = 0; i < 4; i++)
this->x[i] = x[i];
}
void gpt_set(point x1, point x2, point x3, point x4) {
point x[4] = {x1, x2, x3, x4};
gpt_set(x);
}
void snap(ale_pos interval) {
for (int i = 0; i < 4; i++)
for (int j = 0; j < 2; j++)
x[i][j] = round(x[i][j] / interval) * interval;
interval /= scale();
for (int i = 0; i < 2; i++)
eu[i] = round(eu[i] / interval) * interval;
interval *= 2 * 180 / M_PI
/ sqrt(pow(unscaled_height(), 2)
+ pow(unscaled_width(), 2));
eu[2] = round(eu[2] / interval) * interval;
resultant_memo = 0;
resultant_inverse_memo = 0;
}
/*
* Get the specified projective parameter
*/
point gpt_get(int point) const {
assert (_is_projective);
assert (point >= 0);
assert (point < 4);
return x[point];
}
/*
* Get the specified projective parameter
*/
ale_pos gpt_get(int point, int dim) {
assert (_is_projective);
assert (dim >= 0);
assert (dim < 2);
return gpt_get(point)[dim];
}
/*
* Translate by a given amount
*/
void translate(point p) {
resultant_memo = 0;
resultant_inverse_memo = 0;
if (_is_projective)
for (int i = 0; i < 4; i++)
x[i] += p;
else {
eu[0] += p[0] / scale_factor;
eu[1] += p[1] / scale_factor;
}
}
/*
* Rotate by a given amount about a given point.
*/
void rotate(point center, ale_pos degrees) {
if (_is_projective)
for (int i = 0; i <= 4; i++) {
ale_pos radians = (double) degrees * M_PI / (double) 180;
x[i] -= center;
x[i] = point(
(double) x[i][0] * cos(radians) + (double) x[i][1] * sin(radians),
(double) x[i][1] * cos(radians) - (double) x[i][0] * sin(radians));
x[i] += center;
resultant_memo = 0;
resultant_inverse_memo = 0;
} else {
assert(center.defined());
point fixpoint = scaled_inverse_transform(center);
eu_modify(2, degrees);
point offset = center - transform_scaled(fixpoint);
eu_modify(0, offset[0]);
eu_modify(1, offset[1]);
}
}
void reset_memos() {
resultant_memo = 0;
resultant_inverse_memo = 0;
}
/*
* Rescale a transform with a given factor.
*/
void specific_rescale(ale_pos factor) {
resultant_memo = 0;
resultant_inverse_memo = 0;
if (_is_projective) {
for (int i = 0; i < 4; i++)
for (int j = 0; j < 2; j++)
x[i][j] *= factor;
} else {
#if 0
/*
* Euclidean scaling is handled in resultant().
*/
for (int i = 0; i < 2; i++)
eu[i] *= factor;
#endif
}
}
/*
* Set the dimensions of the image.
*/
void specific_set_dimensions(const image *im) {
int new_height = (int) im->height();
int new_width = (int) im->width();
if (_is_projective) {
/*
* If P(w, x, y, z) is a projective transform mapping
* the corners of the unit square to points w, x, y, z,
* and Q(w, x, y, z)(i, j) == P(w, x, y, z)(ai, bj),
* then we have:
*
* Q(w, x, y, z) == P( P(w, x, y, z)(0, 0),
* P(w, x, y, z)(a, 0),
* P(w, x, y, z)(a, b),
* P(w, x, y, z)(0, b) )
*
* If we note that P(w, x, y, z)(0, 0) == w, we can
* omit a calculation.
*
* We take 'a' as the ratio (new_height /
* old_height) and 'b' as the ratio (new_width /
* old_width) if we want the common upper left-hand
* region of both new and old images to map to the same
* area.
*
* Since we're not mapping from the unit square, we
* take 'a' as new_height and 'b' as new_width to
* accommodate the existing scale factor.
*/
point _x, _y, _z;
_x = transform_unscaled(point(new_height, 0 ));
_y = transform_unscaled(point(new_height, new_width));
_z = transform_unscaled(point( 0 , new_width));
x[1] = _x;
x[2] = _y;
x[3] = _z;
}
}
/*
* Modify all projective parameters at once. Accommodate bugs in the
* version 0 transformation file handler (ALE versions 0.4.0p1 and
* earlier). This code is only called when using a transformation data
* file created with an old version of ALE.
*/
void gpt_v0_set(point x[4]) {
_is_projective = 1;
/*
* This is slightly modified code from version
* 0.4.0p1.
*/
ale_pos delta_01 = x[1][0] - x[2][0];
ale_pos delta_02 = x[3][0] - x[2][0];
ale_pos sigma_0 = x[0][0] - x[1][0] + x[2][0] - x[3][0];
ale_pos delta_11 = x[1][1] - x[2][1];
ale_pos delta_12 = x[3][1] - x[2][1];
ale_pos sigma_1 = x[0][1] - x[1][1] + x[2][1] - x[3][1];
g = (sigma_0 * delta_12 - sigma_1 * delta_02)
/ (delta_01 * delta_12 - delta_11 * delta_02)
/ (input_width * scale_factor);
h = (delta_01 * sigma_1 - delta_11 * sigma_0 )
/ (delta_01 * delta_12 - delta_11 * delta_02)
/ (input_height * scale_factor);
a = (x[1][0] - x[0][0] + g * x[1][0])
/ (input_width * scale_factor);
b = (x[3][0] - x[0][0] + h * x[3][0])
/ (input_height * scale_factor);
c = x[0][0];
d = (x[1][1] - x[0][1] + g * x[1][1])
/ (input_width * scale_factor);
e = (x[3][1] - x[0][1] + h * x[3][1])
/ (input_height * scale_factor);
f = x[0][1];
resultant_memo = 1;
resultant_inverse_memo = 0;
this->x[0] = scaled_inverse_transform( point( 0 , 0 ) );
this->x[1] = scaled_inverse_transform( point( (input_height * scale_factor), 0 ) );
this->x[2] = scaled_inverse_transform( point( (input_height * scale_factor), (input_width * scale_factor) ) );
this->x[3] = scaled_inverse_transform( point( 0 , (input_width * scale_factor) ) );
resultant_memo = 0;
resultant_inverse_memo = 0;
}
/*
* Modify all euclidean parameters at once. Accommodate bugs in the
* version 0 transformation file handler (ALE versions 0.4.0p1 and
* earlier). This code is only called when using a transformation data
* file created with an old version of ALE.
*/
void eu_v0_set(ale_pos eu[3]) {
/*
* This is slightly modified code from version
* 0.4.0p1.
*/
int i;
x[0][0] = 0; x[0][1] = 0;
x[1][0] = (input_width * scale_factor); x[1][1] = 0;
x[2][0] = (input_width * scale_factor); x[2][1] = (input_height * scale_factor);
x[3][0] = 0; x[3][1] = (input_height * scale_factor);
/*
* Rotate
*/
ale_pos theta = (double) eu[2] * M_PI / 180;
for (i = 0; i < 4; i++) {
ale_pos _x[2];
_x[0] = ((double) x[i][0] - (double) (input_width * scale_factor)/2) * cos(theta)
+ ((double) x[i][1] - (double) (input_height * scale_factor)/2) * sin(theta)
+ (input_width * scale_factor)/2;
_x[1] = ((double) x[i][1] - (double) (input_height * scale_factor)/2) * cos(theta)
- ((double) x[i][0] - (double) (input_width * scale_factor)/2) * sin(theta)
+ (input_height * scale_factor)/2;
x[i][0] = _x[0];
x[i][1] = _x[1];
}
/*
* Translate
*/
for (i = 0; i < 4; i++) {
x[i][0] += eu[0];
x[i][1] += eu[1];
}
if (_is_projective) {
gpt_v0_set(x);
return;
}
/*
* Reconstruct euclidean parameters
*/
gpt_v0_set(x);
point center((input_height * scale_factor) / 2, (input_width * scale_factor) / 2);
point center_image = transform_scaled(center);
this->eu[0] = (center_image[0] - center[0]) / scale_factor;
this->eu[1] = (center_image[1] - center[1]) / scale_factor;
point center_left((input_height * scale_factor) / 2, 0);
point center_left_image = transform_scaled(center_left);
ale_pos displacement = center_image[0] - center_left_image[0];
this->eu[2] = asin(2 * displacement / (input_width * scale_factor)) / M_PI * 180;
if (center_left_image[1] > center_image[1])
this->eu[2] = this->eu[2] + 180;
resultant_memo = 0;
resultant_inverse_memo = 0;
_is_projective = 0;
}
void debug_output() {
fprintf(stderr, "[t.do ih=%u, iw=%d x=[[%f %f] [%f %f] [%f %f] [%f %f]] eu=[%f %f %f]\n"
" a-f=[%f %f %f %f %f %f %f %f] _a-_f=[%f %f %f %f %f %f %f %f]\n"
" bdcnm=%d ip=%d rm=%d rim=%d sf=%f]\n",
input_height, input_width,
(double) x[0][0], (double) x[0][1],
(double) x[1][0], (double) x[1][1],
(double) x[2][0], (double) x[2][1],
(double) x[3][0], (double) x[3][1],
(double) eu[0], (double) eu[1], (double) eu[2],
(double) a, (double) b, (double) c, (double) d,
(double) e, (double) f, (double) g, (double) h,
(double) _a, (double) _b, (double) _c, (double) _d,
(double) _e, (double) _f, (double) _g, (double) _h,
bd_count(),
_is_projective,
resultant_memo,
resultant_inverse_memo,
(double) scale_factor);
}
};
#endif