1 // Copyright 2002, 2004 David Hilvert <dhilvert@auricle.dyndns.org>,
2 // <dhilvert@ugcs.caltech.edu>
4 /* This file is part of the Anti-Lamenessing Engine.
6 The Anti-Lamenessing Engine is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 3 of the License, or
9 (at your option) any later version.
11 The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with the Anti-Lamenessing Engine; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 * trans_single.h: Represent transformations of the kind q = c(b^-1(p)),
23 * where p is a point in the source coordinate system, q is a point in the
24 * target coordinate system, b^-1 is a transformation correcting barrel
25 * distortion, and c is a transformation of projective or Euclidean type.
26 * (Note that ^-1 in this context indicates the function inverse rather than
30 #ifndef __trans_single_h__
31 #define __trans_single_h__
33 #include "trans_abstract.h"
36 * transformation: a structure to describe a transformation of kind q =
37 * c(b^-1(p)), where p is a point in the source coordinate system, q is a point
38 * in the target coordinate system, b^-1 is a transformation correcting barrel
39 * distortion, and c is a projective or Euclidean transformation. (Note that
40 * ^-1 in this case indicates a function inverse, not exponentiation.) Data
41 * elements are divided into those describing barrel distortion correction and
42 * those describing projective/Euclidean transformations.
44 * Barrel distortion correction estimates barrel distortion using polynomial
45 * functions of distance from the center of an image, following (roughly) the
46 * example set by Helmut Dersch in his PanoTools software:
48 * http://www.path.unimelb.edu.au/~dersch/barrel/barrel.html
50 * Projective transformation data member names roughly correspond to a typical
51 * treatment of projective transformations from:
53 * Heckbert, Paul. "Projective Mappings for Image Warping." Excerpted
54 * from his Master's Thesis (UC Berkeley, 1989). 1995.
56 * http://www.cs.cmu.edu/afs/cs/project/classes-ph/862.95/www/notes/proj.ps
58 * For convenience, Heckbert's 'x' and 'y' are noted here numerically by '0'
59 * and '1', respectively. 'x0' is denoted 'x[0][0]'; 'y0' is 'x[0][1]'.
61 * eu[i] are the parameters for euclidean transformations.
63 * We consider points to be transformed as homogeneous coordinate vectors
64 * multiplied on the right of the transformation matrix, and so we consider the
65 * transformation matrix as
73 * where element i is always equal to 1.
76 struct trans_single
: public trans_abstract
{
80 mutable ale_pos a
, b
, c
, d
, e
, f
, g
, h
; // matrix
81 mutable ale_pos _a
, _b
, _c
, _d
, _e
, _f
, _g
, _h
; // matrix inverse
83 pixel tonal_multiplier
;
84 mutable int resultant_memo
;
85 mutable int resultant_inverse_memo
;
88 * Calculate resultant matrix values.
90 void resultant() const {
93 * If we already know the answers, don't bother calculating
100 int ale_pos_casting
= ale_pos_casting_status();
101 ale_pos_enable_casting();
103 if (_is_projective
) {
106 * Calculate resultant matrix values for a general
107 * projective transformation given that we are mapping
108 * from the source domain of dimension input_height *
109 * input_width to a specified arbitrary quadrilateral.
110 * Follow the calculations outlined in the document by
111 * Paul Heckbert cited above for the case in which the
112 * source domain is a unit square and then divide to
113 * correct for the scale factor in each dimension.
117 * First, perform calculations as outlined in Heckbert.
120 ale_pos delta_01
= x
[1][0] - x
[2][0];
121 ale_pos delta_02
= x
[3][0] - x
[2][0];
122 ale_pos sigma_0
= x
[0][0] - x
[1][0] + x
[2][0] - x
[3][0];
123 ale_pos delta_11
= x
[1][1] - x
[2][1];
124 ale_pos delta_12
= x
[3][1] - x
[2][1];
125 ale_pos sigma_1
= x
[0][1] - x
[1][1] + x
[2][1] - x
[3][1];
127 g
= (sigma_0
* delta_12
- sigma_1
* delta_02
)
128 / (delta_01
* delta_12
- delta_11
* delta_02
);
130 h
= (delta_01
* sigma_1
- delta_11
* sigma_0
)
131 / (delta_01
* delta_12
- delta_11
* delta_02
);
133 a
= (x
[1][0] - x
[0][0] + g
* x
[1][0]);
134 b
= (x
[3][0] - x
[0][0] + h
* x
[3][0]);
137 d
= (x
[1][1] - x
[0][1] + g
* x
[1][1]);
138 e
= (x
[3][1] - x
[0][1] + h
* x
[3][1]);
142 * Finish by scaling so that our transformation maps
143 * from a rectangle of width and height matching the
144 * width and height of the input image.
157 * Calculate matrix values for a euclidean
160 * We want to translate the image center by (eu[0],
161 * eu[1]) and rotate the image about the center by
162 * eu[2] degrees. This is equivalent to the following
163 * sequence of affine transformations applied to the
164 * point to be transformed:
166 * translate by (-h/2, -w/2)
167 * rotate by eu[2] degrees about the origin
168 * translate by (h/2, w/2)
169 * translate by (eu[0], eu[1])
171 * The matrix assigned below represents the result of
172 * combining all of these transformations. Matrix
173 * elements g and h are always zero in an affine
177 ale_pos theta
= (double) eu
[2] * M_PI
/ 180;
179 a
= (double) cos(theta
) * (double) scale_factor
;
180 b
= (double) sin(theta
) * (double) scale_factor
;
181 c
= 0.5 * ((double) input_height
* ((double) scale_factor
- (double) a
)
182 - (double) input_width
* (double) b
) + (double) eu
[0]
183 * (double) scale_factor
;
186 f
= 0.5 * ((double) input_height
* (double) b
187 + (double) input_width
* ((double) scale_factor
189 + (double) eu
[1] * (double) scale_factor
;
196 if (!ale_pos_casting
)
197 ale_pos_disable_casting();
201 * Calculate the inverse transform matrix values.
203 void resultant_inverse () const {
206 * If we already know the answers, don't bother calculating
210 if (resultant_inverse_memo
)
215 int ale_pos_casting
= ale_pos_casting_status();
216 ale_pos_enable_casting();
219 * For projective transformations, we calculate
220 * the inverse of the forward transformation
224 double scale
= (double) a
* (double) e
- (double) b
* (double) d
;
226 _a
= ((double) e
* 1 - (double) f
* (double) h
) / scale
;
227 _b
= ((double) h
* (double) c
- 1 * (double) b
) / scale
;
228 _c
= ((double) b
* (double) f
- (double) c
* (double) e
) / scale
;
229 _d
= ((double) f
* (double) g
- (double) d
* 1) / scale
;
230 _e
= (1 * (double) a
- (double) g
* (double) c
) / scale
;
231 _f
= ((double) c
* (double) d
- (double) a
* (double) f
) / scale
;
232 _g
= ((double) d
* (double) h
- (double) e
* (double) g
) / scale
;
233 _h
= ((double) g
* (double) b
- (double) h
* (double) a
) / scale
;
235 resultant_inverse_memo
= 1;
237 if (!ale_pos_casting
)
238 ale_pos_disable_casting();
243 trans_single
&operator=(const trans_single
&ta
) {
245 this->trans_abstract::operator=(*((trans_abstract
*) &ta
));
247 for (int i
= 0; i
< 4; i
++) {
251 for (int i
= 0; i
< 3; i
++) {
255 _is_projective
= ta
._is_projective
;
257 tonal_multiplier
= ta
.tonal_multiplier
;
260 resultant_inverse_memo
= 0;
265 trans_single(const trans_single
&ta
) {
274 * Returns non-zero if the transformation might be non-Euclidean.
276 int is_projective() const {
277 return _is_projective
;
281 * Projective/Euclidean transformation
283 struct point
pe(struct point p
) const {
288 result
[0] = (a
* p
[0] + b
* p
[1] + c
)
289 / (g
* p
[0] + h
* p
[1] + 1);
290 result
[1] = (d
* p
[0] + e
* p
[1] + f
)
291 / (g
* p
[0] + h
* p
[1] + 1);
297 * Projective/Euclidean inverse
299 struct point
pei(struct point p
) const {
304 result
[0] = (_a
* p
[0] + _b
* p
[1] + _c
)
305 / (_g
* p
[0] + _h
* p
[1] + 1);
306 result
[1] = (_d
* p
[0] + _e
* p
[1] + _f
)
307 / (_g
* p
[0] + _h
* p
[1] + 1);
315 * operator() is the transformation operator.
317 struct point
operator()(struct point p
) {
323 * Calculate projective transformation parameters from a euclidean
328 assert(!_is_projective
);
330 x
[0] = transform_unscaled(point( 0 , 0 ) );
331 x
[1] = transform_unscaled(point( input_height
, 0 ) );
332 x
[2] = transform_unscaled(point( input_height
, input_width
) );
333 x
[3] = transform_unscaled(point( 0 , input_width
) );
336 resultant_inverse_memo
= 0;
342 * Identity transform for given dimensions
344 static struct trans_single
eu_identity(unsigned int height
= 2, unsigned int width
= 2, ale_pos scale_factor
= 1) {
345 struct trans_single r
;
347 r
.resultant_memo
= 0;
348 r
.resultant_inverse_memo
= 0;
353 r
.input_width
= width
;
354 r
.input_height
= height
;
355 r
.scale_factor
= scale_factor
;
357 r
._is_projective
= 0;
359 r
.tonal_multiplier
= pixel(1, 1, 1);
361 r
.bd_set(0, (ale_pos
*) NULL
);
367 * Calculate euclidean identity transform for a given image.
369 static struct trans_single
eu_identity(const image
*i
, ale_pos scale_factor
= 1) {
371 return eu_identity(2, 2, scale_factor
);
373 return eu_identity(i
->height(), i
->width(), scale_factor
);
377 * Projective identity transform for given dimensions.
379 static trans_single
gpt_identity(unsigned int height
, unsigned int width
, ale_pos scale_factor
) {
380 struct trans_single r
= eu_identity(height
, width
, scale_factor
);
388 * Calculate projective identity transform for a given image.
390 static trans_single
gpt_identity(const image
*i
, ale_pos scale_factor
) {
391 struct trans_single r
= eu_identity(i
, scale_factor
);
399 * Set the tonal multiplier
401 void set_tonal_multiplier(pixel p
) {
402 tonal_multiplier
= p
;
405 pixel
get_tonal_multiplier(struct point p
) const {
406 return tonal_multiplier
;
409 pixel
get_inverse_tonal_multiplier(struct point p
) const {
410 return tonal_multiplier
;
414 * Modify a euclidean transform in the indicated manner.
416 void eu_modify(int i1
, ale_pos diff
) {
418 assert(!_is_projective
);
421 resultant_inverse_memo
= 0;
424 eu
[i1
] += diff
/ scale_factor
;
431 * Rotate about a given point in the original reference frame.
433 void eu_rotate_about_scaled(point center
, ale_pos diff
) {
434 assert(center
.defined());
435 point fixpoint
= scaled_inverse_transform(center
);
437 point offset
= center
- transform_scaled(fixpoint
);
439 eu_modify(0, offset
[0]);
440 eu_modify(1, offset
[1]);
444 * Modify all euclidean parameters at once.
446 void eu_set(ale_pos eu
[3]) {
449 resultant_inverse_memo
= 0;
451 this->eu
[0] = eu
[0] / scale_factor
;
452 this->eu
[1] = eu
[1] / scale_factor
;
455 if (_is_projective
) {
463 * Get the specified euclidean parameter
465 ale_pos
eu_get(int param
) const {
466 assert (!_is_projective
);
471 return eu
[param
] * scale_factor
;
477 * Modify a projective transform in the indicated manner.
479 void gpt_modify(int i1
, int i2
, ale_pos diff
) {
481 assert (_is_projective
);
484 resultant_inverse_memo
= 0;
491 * Modify projective transform with control points at the corners of a
492 * specified bounded region.
495 void gpt_modify_bounded(int i1
, int i2
, ale_pos diff
, elem_bounds_int_t eb
) {
499 trans_single t
= gpt_identity(eb
.imax
- eb
.imin
, eb
.jmax
- eb
.jmin
, 1);
501 t
.gpt_modify(i1
, i2
, diff
);
503 for (int p
= 0; p
< 4; p
++)
504 new_x
[p
] = t
.transform_scaled(x
[p
] - point(eb
.imin
, eb
.jmin
)) + point(eb
.imin
, eb
.jmin
);
511 * Modify a projective transform according to the group operation.
513 void gr_modify(int i1
, int i2
, ale_pos diff
) {
514 assert (_is_projective
);
515 assert (i1
== 0 || i1
== 1);
517 point diff_vector
= (i1
== 0)
521 trans_single t
= *this;
523 t
.resultant_memo
= 0;
524 t
.resultant_inverse_memo
= 0;
525 t
.input_height
= (unsigned int) scaled_height();
526 t
.input_width
= (unsigned int) scaled_width();
528 t
.bd_set(0, (ale_pos
*) NULL
);
531 resultant_inverse_memo
= 0;
533 x
[i2
] = t
.transform_scaled(t
.scaled_inverse_transform(x
[i2
]) + diff_vector
);
537 * Modify all projective parameters at once.
539 void gpt_set(point x
[4]) {
542 resultant_inverse_memo
= 0;
546 for (int i
= 0; i
< 4; i
++)
551 void gpt_set(point x1
, point x2
, point x3
, point x4
) {
552 point x
[4] = {x1
, x2
, x3
, x4
};
556 void snap(ale_pos interval
) {
557 for (int i
= 0; i
< 4; i
++)
558 for (int j
= 0; j
< 2; j
++)
559 x
[i
][j
] = round(x
[i
][j
] / interval
) * interval
;
563 for (int i
= 0; i
< 2; i
++)
564 eu
[i
] = round(eu
[i
] / interval
) * interval
;
566 interval
*= 2 * 180 / M_PI
567 / sqrt(pow(unscaled_height(), 2)
568 + pow(unscaled_width(), 2));
570 eu
[2] = round(eu
[2] / interval
) * interval
;
573 resultant_inverse_memo
= 0;
577 * Get the specified projective parameter
579 point
gpt_get(int point
) const {
580 assert (_is_projective
);
588 * Get the specified projective parameter
590 ale_pos
gpt_get(int point
, int dim
) {
591 assert (_is_projective
);
595 return gpt_get(point
)[dim
];
599 * Translate by a given amount
601 void translate(point p
) {
604 resultant_inverse_memo
= 0;
607 for (int i
= 0; i
< 4; i
++)
610 eu
[0] += p
[0] / scale_factor
;
611 eu
[1] += p
[1] / scale_factor
;
616 * Rotate by a given amount about a given point.
618 void rotate(point center
, ale_pos degrees
) {
620 for (int i
= 0; i
<= 4; i
++) {
621 ale_pos radians
= (double) degrees
* M_PI
/ (double) 180;
625 (double) x
[i
][0] * cos(radians
) + (double) x
[i
][1] * sin(radians
),
626 (double) x
[i
][1] * cos(radians
) - (double) x
[i
][0] * sin(radians
));
630 resultant_inverse_memo
= 0;
632 assert(center
.defined());
633 point fixpoint
= scaled_inverse_transform(center
);
634 eu_modify(2, degrees
);
635 point offset
= center
- transform_scaled(fixpoint
);
637 eu_modify(0, offset
[0]);
638 eu_modify(1, offset
[1]);
644 resultant_inverse_memo
= 0;
649 * Rescale a transform with a given factor.
651 void specific_rescale(ale_pos factor
) {
654 resultant_inverse_memo
= 0;
656 if (_is_projective
) {
658 for (int i
= 0; i
< 4; i
++)
659 for (int j
= 0; j
< 2; j
++)
665 * Euclidean scaling is handled in resultant().
667 for (int i
= 0; i
< 2; i
++)
674 * Set the dimensions of the image.
676 void specific_set_dimensions(const image
*im
) {
678 int new_height
= (int) im
->height();
679 int new_width
= (int) im
->width();
681 if (_is_projective
) {
684 * If P(w, x, y, z) is a projective transform mapping
685 * the corners of the unit square to points w, x, y, z,
686 * and Q(w, x, y, z)(i, j) == P(w, x, y, z)(ai, bj),
689 * Q(w, x, y, z) == P( P(w, x, y, z)(0, 0),
690 * P(w, x, y, z)(a, 0),
691 * P(w, x, y, z)(a, b),
692 * P(w, x, y, z)(0, b) )
694 * If we note that P(w, x, y, z)(0, 0) == w, we can
695 * omit a calculation.
697 * We take 'a' as the ratio (new_height /
698 * old_height) and 'b' as the ratio (new_width /
699 * old_width) if we want the common upper left-hand
700 * region of both new and old images to map to the same
703 * Since we're not mapping from the unit square, we
704 * take 'a' as new_height and 'b' as new_width to
705 * accommodate the existing scale factor.
710 _x
= transform_unscaled(point(new_height
, 0 ));
711 _y
= transform_unscaled(point(new_height
, new_width
));
712 _z
= transform_unscaled(point( 0 , new_width
));
722 * Modify all projective parameters at once. Accommodate bugs in the
723 * version 0 transformation file handler (ALE versions 0.4.0p1 and
724 * earlier). This code is only called when using a transformation data
725 * file created with an old version of ALE.
727 void gpt_v0_set(point x
[4]) {
732 * This is slightly modified code from version
736 ale_pos delta_01
= x
[1][0] - x
[2][0];
737 ale_pos delta_02
= x
[3][0] - x
[2][0];
738 ale_pos sigma_0
= x
[0][0] - x
[1][0] + x
[2][0] - x
[3][0];
739 ale_pos delta_11
= x
[1][1] - x
[2][1];
740 ale_pos delta_12
= x
[3][1] - x
[2][1];
741 ale_pos sigma_1
= x
[0][1] - x
[1][1] + x
[2][1] - x
[3][1];
743 g
= (sigma_0
* delta_12
- sigma_1
* delta_02
)
744 / (delta_01
* delta_12
- delta_11
* delta_02
)
745 / (input_width
* scale_factor
);
747 h
= (delta_01
* sigma_1
- delta_11
* sigma_0
)
748 / (delta_01
* delta_12
- delta_11
* delta_02
)
749 / (input_height
* scale_factor
);
751 a
= (x
[1][0] - x
[0][0] + g
* x
[1][0])
752 / (input_width
* scale_factor
);
753 b
= (x
[3][0] - x
[0][0] + h
* x
[3][0])
754 / (input_height
* scale_factor
);
757 d
= (x
[1][1] - x
[0][1] + g
* x
[1][1])
758 / (input_width
* scale_factor
);
759 e
= (x
[3][1] - x
[0][1] + h
* x
[3][1])
760 / (input_height
* scale_factor
);
764 resultant_inverse_memo
= 0;
766 this->x
[0] = scaled_inverse_transform( point( 0 , 0 ) );
767 this->x
[1] = scaled_inverse_transform( point( (input_height
* scale_factor
), 0 ) );
768 this->x
[2] = scaled_inverse_transform( point( (input_height
* scale_factor
), (input_width
* scale_factor
) ) );
769 this->x
[3] = scaled_inverse_transform( point( 0 , (input_width
* scale_factor
) ) );
772 resultant_inverse_memo
= 0;
776 * Modify all euclidean parameters at once. Accommodate bugs in the
777 * version 0 transformation file handler (ALE versions 0.4.0p1 and
778 * earlier). This code is only called when using a transformation data
779 * file created with an old version of ALE.
781 void eu_v0_set(ale_pos eu
[3]) {
784 * This is slightly modified code from version
790 x
[0][0] = 0; x
[0][1] = 0;
791 x
[1][0] = (input_width
* scale_factor
); x
[1][1] = 0;
792 x
[2][0] = (input_width
* scale_factor
); x
[2][1] = (input_height
* scale_factor
);
793 x
[3][0] = 0; x
[3][1] = (input_height
* scale_factor
);
799 ale_pos theta
= (double) eu
[2] * M_PI
/ 180;
801 for (i
= 0; i
< 4; i
++) {
804 _x
[0] = ((double) x
[i
][0] - (double) (input_width
* scale_factor
)/2) * cos(theta
)
805 + ((double) x
[i
][1] - (double) (input_height
* scale_factor
)/2) * sin(theta
)
806 + (input_width
* scale_factor
)/2;
807 _x
[1] = ((double) x
[i
][1] - (double) (input_height
* scale_factor
)/2) * cos(theta
)
808 - ((double) x
[i
][0] - (double) (input_width
* scale_factor
)/2) * sin(theta
)
809 + (input_height
* scale_factor
)/2;
819 for (i
= 0; i
< 4; i
++) {
824 if (_is_projective
) {
830 * Reconstruct euclidean parameters
835 point
center((input_height
* scale_factor
) / 2, (input_width
* scale_factor
) / 2);
836 point center_image
= transform_scaled(center
);
838 this->eu
[0] = (center_image
[0] - center
[0]) / scale_factor
;
839 this->eu
[1] = (center_image
[1] - center
[1]) / scale_factor
;
841 point
center_left((input_height
* scale_factor
) / 2, 0);
842 point center_left_image
= transform_scaled(center_left
);
844 ale_pos displacement
= center_image
[0] - center_left_image
[0];
846 this->eu
[2] = asin(2 * displacement
/ (input_width
* scale_factor
)) / M_PI
* 180;
848 if (center_left_image
[1] > center_image
[1])
849 this->eu
[2] = this->eu
[2] + 180;
852 resultant_inverse_memo
= 0;
858 void debug_output() {
859 fprintf(stderr
, "[t.do ih=%u, iw=%d x=[[%f %f] [%f %f] [%f %f] [%f %f]] eu=[%f %f %f]\n"
860 " a-f=[%f %f %f %f %f %f %f %f] _a-_f=[%f %f %f %f %f %f %f %f]\n"
861 " bdcnm=%d ip=%d rm=%d rim=%d sf=%f]\n",
862 input_height
, input_width
,
863 (double) x
[0][0], (double) x
[0][1],
864 (double) x
[1][0], (double) x
[1][1],
865 (double) x
[2][0], (double) x
[2][1],
866 (double) x
[3][0], (double) x
[3][1],
867 (double) eu
[0], (double) eu
[1], (double) eu
[2],
868 (double) a
, (double) b
, (double) c
, (double) d
,
869 (double) e
, (double) f
, (double) g
, (double) h
,
870 (double) _a
, (double) _b
, (double) _c
, (double) _d
,
871 (double) _e
, (double) _f
, (double) _g
, (double) _h
,
875 resultant_inverse_memo
,
876 (double) scale_factor
);