d2/trans_abstract.h

   1 // Copyright 2002, 2004, 2007 David Hilvert <dhilvert@auricle.dyndns.org>,
   2 //                                          <dhilvert@ugcs.caltech.edu>
   3
   4 /*  This file is part of the Anti-Lamenessing Engine.
   5
   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.
  10
  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.
  15
  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
  19 */
  20
  21 /*
  22  * trans_abstract.h: Abstract transformation superclass.
  23  */
  24
  25 #ifndef __trans_abstract_h__
  26 #define __trans_abstract_h__
  27
  28 #include "image.h"
  29 #include "point.h"
  30
  31 #ifndef M_PI
  32 #define M_PI 3.14159265358979323846
  33 #endif
  34
  35 /*
  36  * Number of coefficients used in correcting barrel distortion.
  37  */
  38
  39 #define BARREL_DEGREE 5
  40
  41 /*
  42  * Acceptable error for inverse barrel distortion, measured in scaled output
  43  * pixels.
  44  */
  45
  46 #define BARREL_INV_ERROR 0.01
  47
  48 struct trans_abstract {
  49 private:
  50         ale_pos bdc[BARREL_DEGREE];                             // barrel-dist. coeffs.
  51         unsigned int bdcnum;                                    // number of bdcs
  52
  53 protected:
  54         ale_pos scale_factor;
  55         unsigned int input_height, input_width;
  56
  57         virtual void specific_rescale(ale_pos factor) = 0;
  58         virtual void reset_memos() = 0;
  59         virtual void specific_set_dimensions(const image *im) = 0;
  60
  61 public:
  62
  63         trans_abstract() {
  64                 bdcnum = 0;
  65         }
  66
  67         trans_abstract &operator=(const trans_abstract &ta) {
  68                 scale_factor = ta.scale_factor;
  69                 input_width = ta.input_width;
  70                 input_height = ta.input_height;
  71
  72                 bdcnum = ta.bdcnum;
  73
  74                 assert (bdcnum < BARREL_DEGREE);
  75
  76                 for (unsigned int d = 0; d < bdcnum; d++)
  77                         bdc[d] = ta.bdc[d];
  78
  79                 return *this;
  80         }
  81
  82         trans_abstract (const trans_abstract &ta) {
  83                 operator=(ta);
  84         }
  85
  86         /*
  87          * Returns non-zero if the transformation might be non-Euclidean.
  88          */
  89         virtual int is_projective() const = 0;
  90
  91         /*
  92          * Get scale factor.
  93          */
  94
  95         ale_pos scale() const {
  96                 return scale_factor;
  97         }
  98
  99         /*
 100          * Get width of input image.
 101          */
 102         ale_pos scaled_width() const {
 103                 return (input_width * scale_factor);
 104         }
 105
 106         /*
 107          * Get unscaled width of input image.
 108          */
 109         unsigned int unscaled_width() const {
 110                 return (unsigned int) input_width;
 111         }
 112
 113         /*
 114          * Get height of input image;
 115          */
 116         ale_pos scaled_height() const {
 117                return (input_height * scale_factor);
 118         }
 119
 120         /*
 121          * Get unscaled height of input image.
 122          */
 123         unsigned int unscaled_height() const {
 124                 return (unsigned int) input_height;
 125         }
 126
 127         /*
 128          * Barrel distortion radial component.
 129          */
 130         ale_pos bdr(ale_pos r) const {
 131                 assert (bdcnum < BARREL_DEGREE);
 132                 ale_pos s = r;
 133                 for (unsigned int d = 0; d < bdcnum; d++)
 134                         s += bdc[d] * (pow(r, d + 2) - r);
 135                 return s;
 136         }
 137
 138         /*
 139          * Derivative of the barrel distortion radial component.
 140          */
 141         ale_pos bdrd(ale_pos r) const {
 142                 assert (bdcnum < BARREL_DEGREE);
 143                 ale_pos s = 1;
 144                 for (unsigned int d = 0; d < bdcnum; d++)
 145                         s += bdc[d] * (pow(r, d + 1) - 1);
 146                 return s;
 147         }
 148
 149         /*
 150          * Barrel distortion.
 151          */
 152         struct point bd(struct point p) const {
 153                 if (bdcnum > 0) {
 154                         point half_diag = point(unscaled_height(), unscaled_width()) / 2;
 155
 156                         p -= half_diag;
 157
 158                         ale_pos r = p.norm() / half_diag.norm();
 159
 160                         if (r > 0.00001)
 161                                 p *= bdr(r)/r;
 162
 163                         p += half_diag;
 164                 }
 165
 166                 assert (!isnan(p[0]) && !isnan(p[1]));
 167
 168                 return p;
 169         }
 170
 171         /*
 172          * Barrel distortion inverse.
 173          */
 174         struct point bdi(struct point p) const {
 175                 if (bdcnum > 0) {
 176                         point half_diag = point(unscaled_height(), unscaled_width()) / 2;
 177
 178                         p -= half_diag;
 179
 180                         ale_pos r = p.norm() / half_diag.norm();
 181                         ale_pos s = r;
 182
 183                         while (fabs(r - bdr(s)) * half_diag.norm() > BARREL_INV_ERROR)
 184                                 s += (r - bdr(s)) / bdrd(s);
 185
 186                         if (r > 0.0001)
 187                                 p *= s / r;
 188
 189                         p += half_diag;
 190                 }
 191
 192                 assert (!isnan(p[0]) && !isnan(p[1]));
 193
 194                 return p;
 195         }
 196
 197         /*
 198          * Transformation sans barrel distortion
 199          */
 200         virtual struct point pe(struct point p) const = 0;
 201
 202         /*
 203          * Transformation inverse sans barrel distortion
 204          */
 205         virtual struct point pei(struct point p) const = 0;
 206
 207         /*
 208          * Map unscaled point p.
 209          */
 210         struct point transform_unscaled(struct point p) const {
 211                 return pe(bdi(p));
 212         }
 213
 214         /*
 215          * Transform point p.
 216          *
 217          * Barrel distortion correction followed by a projective/euclidean
 218          * transformation.
 219          */
 220         struct point transform_scaled(struct point p) const {
 221                 return transform_unscaled(p / scale_factor);
 222         }
 223
 224 #if 0
 225         /*
 226          * operator() is the transformation operator.
 227          */
 228         struct point operator()(struct point p) {
 229                 return transform(p);
 230         }
 231 #endif
 232
 233         /*
 234          * Map point p using the inverse of the transform into
 235          * the unscaled image space.
 236          */
 237         struct point unscaled_inverse_transform(struct point p) const {
 238                 return bd(pei(p));
 239         }
 240
 241         /*
 242          * Map point p using the inverse of the transform.
 243          *
 244          * Projective/euclidean inverse followed by barrel distortion.
 245          */
 246         struct point scaled_inverse_transform(struct point p) const {
 247                 assert (p.defined());
 248                 point q = unscaled_inverse_transform(p);
 249
 250                 q[0] *= scale_factor;
 251                 q[1] *= scale_factor;
 252
 253                 return q;
 254         }
 255
 256         /*
 257          * Calculate projective transformation parameters from a euclidean
 258          * transformation.
 259          */
 260         virtual void eu_to_gpt() = 0;
 261
 262         /*
 263          * Modify a euclidean transform in the indicated manner.
 264          */
 265         virtual void eu_modify(int i1, ale_pos diff) = 0;
 266
 267         /*
 268          * Rotate about a given point in the original reference frame.
 269          */
 270         virtual void eu_rotate_about_scaled(point center, ale_pos diff) = 0;
 271
 272         /*
 273          * Modify all euclidean parameters at once.
 274          */
 275         virtual void eu_set(ale_pos eu[3]) = 0;
 276
 277         /*
 278          * Get the specified euclidean parameter
 279          */
 280         virtual ale_pos eu_get(int param) const = 0;
 281
 282         /*
 283          * Modify a projective transform in the indicated manner.
 284          */
 285         virtual void gpt_modify(int i1, int i2, ale_pos diff) = 0;
 286
 287         /*
 288          * Modify a projective transform according to the group operation.
 289          */
 290         virtual void gr_modify(int i1, int i2, ale_pos diff) = 0;
 291
 292         /*
 293          * Modify all projective parameters at once.
 294          */
 295         virtual void gpt_set(point x[4]) = 0;
 296
 297         virtual void gpt_set(point x1, point x2, point x3, point x4) = 0;
 298
 299         /*
 300          * Snap positional parameters to the specified resolution.
 301          */
 302
 303         virtual void snap(ale_pos interval) = 0;
 304
 305         /*
 306          * Get the specified projective parameter
 307          */
 308         virtual point gpt_get(int point) const = 0;
 309
 310         /*
 311          * Get the specified projective parameter
 312          */
 313         virtual ale_pos gpt_get(int point, int dim) = 0;
 314
 315         /*
 316          * Check equality of transformation parameters.
 317          */
 318         virtual int operator==(const trans_abstract &t) const {
 319                 /*
 320                  * Small tolerances (< 10^-6?) can cause odd errors,
 321                  * possibly due to float<->double conversion issues.
 322                  */
 323                 double zero_tolerance = 0.01;
 324
 325                 if (scale() != t.scale())
 326                         return 0;
 327
 328                 if (is_projective() != t.is_projective())
 329                         return 0;
 330
 331                 if (is_projective()) {
 332                         assert (t.is_projective());
 333                         for (int i = 0; i < 4; i++)
 334                         for (int d = 0; d < 2; d++) {
 335                                 double abs_difference = fabs(gpt_get(i)[d] - t.gpt_get(i)[d]);
 336
 337                                 if (abs_difference > zero_tolerance)
 338                                         return 0;
 339                         }
 340                 } else {
 341                         assert (!t.is_projective());
 342                         for (int i = 0; i < 3; i++) {
 343                                 double abs_difference = fabs(eu_get(i) - t.eu_get(i));
 344
 345                                 if (abs_difference > zero_tolerance)
 346                                         return 0;
 347                         }
 348                 }
 349
 350                 return 1;
 351         }
 352
 353         virtual int operator!=(const trans_abstract &t) const {
 354                 return !(operator==(t));
 355         }
 356
 357
 358         /*
 359          * Translate by a given amount
 360          */
 361         virtual void translate(point p) = 0;
 362
 363         /*
 364          * Rotate by a given amount about a given point.
 365          */
 366         virtual void rotate(point p, ale_pos degrees) = 0;
 367
 368         /*
 369          * Set the specified barrel distortion parameter.
 370          */
 371         void bd_set(unsigned int degree, ale_pos value) {
 372                 assert (degree < bdcnum);
 373                 bdc[degree] = value;
 374         }
 375
 376         /*
 377          * Set all barrel distortion parameters.
 378          */
 379         void bd_set(unsigned int degree, ale_pos values[BARREL_DEGREE]) {
 380                 assert (degree <= BARREL_DEGREE);
 381                 bdcnum = degree;
 382                 for (unsigned int d = 0; d < degree; d++)
 383                         bdc[d] = values[d];
 384         }
 385
 386         /*
 387          * Get all barrel distortion parameters.
 388          */
 389         void bd_get(ale_pos result[BARREL_DEGREE]) {
 390                 for (unsigned int d = 0; d < bdcnum; d++)
 391                         result[d] = bdc[d];
 392         }
 393
 394         /*
 395          * Get the specified barrel distortion parameter.
 396          */
 397         ale_pos bd_get(unsigned int degree) {
 398                 assert (degree < bdcnum);
 399                 return bdc[degree];
 400         }
 401
 402         /*
 403          * Get the number of barrel distortion parameters.
 404          */
 405         unsigned int bd_count() {
 406                 return bdcnum;
 407         }
 408
 409         /*
 410          * Get the maximum allowable number of barrel distortion parameters.
 411          */
 412         unsigned int bd_max() {
 413                 return BARREL_DEGREE;
 414         }
 415
 416         /*
 417          * Modify the specified barrel distortion parameter.
 418          */
 419         void bd_modify(unsigned int degree, ale_pos diff) {
 420                 assert (degree < bdcnum);
 421                 bd_set(degree, bd_get(degree) + diff);
 422         }
 423
 424         /*
 425          * Rescale a transform with a given factor.
 426          */
 427         void rescale(ale_pos factor) {
 428                 specific_rescale(factor);
 429                 scale_factor *= factor;
 430         }
 431
 432         /*
 433          * Set a new domain.
 434          */
 435
 436         void set_domain(unsigned int new_height, unsigned int new_width) {
 437                 reset_memos();
 438                 input_width = new_width;
 439                 input_height = new_height;
 440         }
 441
 442         /*
 443          * Set the dimensions of the image.
 444          */
 445         void set_dimensions(const image *im) {
 446
 447                 int new_height = (int) im->height();
 448                 int new_width  = (int) im->width();
 449
 450                 reset_memos();
 451                 specific_set_dimensions(im);
 452                 input_height = new_height;
 453                 input_width  = new_width;
 454         }
 455
 456         /*
 457          * Get the position and dimensions of a pixel P mapped from one
 458          * coordinate system to another, using the forward transformation.
 459          * This function uses scaled input coordinates.
 460          */
 461         void map_area(point p, point *q, ale_pos d[2]) {
 462
 463                 /*
 464                  * Determine the coordinates in the target frame for the source
 465                  * image pixel P and two adjacent source pixels.
 466                  */
 467
 468                     (*q) = transform_scaled(p);
 469                 point q0 = transform_scaled(point(p[0] + 1, p[1]));
 470                 point q1 = transform_scaled(point(p[0], p[1] + 1));
 471
 472                 /*
 473                  * Calculate the distance between source image pixel and
 474                  * adjacent source pixels, measured in the coordinate system of
 475                  * the target frame.
 476                  */
 477
 478                 ale_pos ui = fabs(q0[0] - (*q)[0]);
 479                 ale_pos uj = fabs(q0[1] - (*q)[1]);
 480                 ale_pos vi = fabs(q1[0] - (*q)[0]);
 481                 ale_pos vj = fabs(q1[1] - (*q)[1]);
 482
 483                 /*
 484                  * We map the area of the source image pixel P onto the target
 485                  * frame as a rectangular area oriented on the target frame's
 486                  * axes.  Note that this results in an area that may be the
 487                  * wrong shape or orientation.
 488                  *
 489                  * We define two estimates of the rectangle's dimensions below.
 490                  * For rotations of 0, 90, 180, or 270 degrees, max and sum are
 491                  * identical.  For other orientations, sum is too large and max
 492                  * is too small.  We use the mean of max and sum, which we then
 493                  * divide by two to obtain the distance between the center and
 494                  * the edge.
 495                  */
 496
 497                 ale_pos maxi = (ui > vi) ? ui : vi;
 498                 ale_pos maxj = (uj > vj) ? uj : vj;
 499                 ale_pos sumi = ui + vi;
 500                 ale_pos sumj = uj + vj;
 501
 502                 d[0] = (maxi + sumi) / 4;
 503                 d[1] = (maxj + sumj) / 4;
 504         }
 505
 506         /*
 507          * Get the position and dimensions of a pixel P mapped from one
 508          * coordinate system to another, using the forward transformation.
 509          * This function uses unscaled input coordinates.
 510          */
 511         void map_area_unscaled(point p, point *q, ale_pos d[2]) {
 512
 513                 /*
 514                  * Determine the coordinates in the target frame for the source
 515                  * image pixel P and two adjacent source pixels.
 516                  */
 517
 518                     (*q) = transform_unscaled(p);
 519                 point q0 = transform_unscaled(point(p[0] + 1, p[1]));
 520                 point q1 = transform_unscaled(point(p[0], p[1] + 1));
 521
 522                 /*
 523                  * Calculate the distance between source image pixel and
 524                  * adjacent source pixels, measured in the coordinate system of
 525                  * the target frame.
 526                  */
 527
 528                 ale_pos ui = fabs(q0[0] - (*q)[0]);
 529                 ale_pos uj = fabs(q0[1] - (*q)[1]);
 530                 ale_pos vi = fabs(q1[0] - (*q)[0]);
 531                 ale_pos vj = fabs(q1[1] - (*q)[1]);
 532
 533                 /*
 534                  * We map the area of the source image pixel P onto the target
 535                  * frame as a rectangular area oriented on the target frame's
 536                  * axes.  Note that this results in an area that may be the
 537                  * wrong shape or orientation.
 538                  *
 539                  * We define two estimates of the rectangle's dimensions below.
 540                  * For rotations of 0, 90, 180, or 270 degrees, max and sum are
 541                  * identical.  For other orientations, sum is too large and max
 542                  * is too small.  We use the mean of max and sum, which we then
 543                  * divide by two to obtain the distance between the center and
 544                  * the edge.
 545                  */
 546
 547                 ale_pos maxi = (ui > vi) ? ui : vi;
 548                 ale_pos maxj = (uj > vj) ? uj : vj;
 549                 ale_pos sumi = ui + vi;
 550                 ale_pos sumj = uj + vj;
 551
 552                 d[0] = (maxi + sumi) / 4;
 553                 d[1] = (maxj + sumj) / 4;
 554         }
 555
 556         /*
 557          * Get the position and dimensions of a pixel P mapped from one
 558          * coordinate system to another, using the inverse transformation.  If
 559          * SCALE_FACTOR is not equal to one, divide out the scale factor to
 560          * obtain unscaled coordinates.  This method is very similar to the
 561          * map_area method above.
 562          */
 563         void unscaled_map_area_inverse(point p, point *q, ale_pos d[2]) {
 564
 565                 /*
 566                  * Determine the coordinates in the target frame for the source
 567                  * image pixel P and two adjacent source pixels.
 568                  */
 569
 570                     (*q) = scaled_inverse_transform(p);
 571                 point q0 = scaled_inverse_transform(point(p[0] + 1, p[1]));
 572                 point q1 = scaled_inverse_transform(point(p[0], p[1] + 1));
 573
 574
 575                 /*
 576                  * Calculate the distance between source image pixel and
 577                  * adjacent source pixels, measured in the coordinate system of
 578                  * the target frame.
 579                  */
 580
 581                 ale_pos ui = fabs(q0[0] - (*q)[0]);
 582                 ale_pos uj = fabs(q0[1] - (*q)[1]);
 583                 ale_pos vi = fabs(q1[0] - (*q)[0]);
 584                 ale_pos vj = fabs(q1[1] - (*q)[1]);
 585
 586                 /*
 587                  * We map the area of the source image pixel P onto the target
 588                  * frame as a rectangular area oriented on the target frame's
 589                  * axes.  Note that this results in an area that may be the
 590                  * wrong shape or orientation.
 591                  *
 592                  * We define two estimates of the rectangle's dimensions below.
 593                  * For rotations of 0, 90, 180, or 270 degrees, max and sum are
 594                  * identical.  For other orientations, sum is too large and max
 595                  * is too small.  We use the mean of max and sum, which we then
 596                  * divide by two to obtain the distance between the center and
 597                  * the edge.
 598                  */
 599
 600                 ale_pos maxi = (ui > vi) ? ui : vi;
 601                 ale_pos maxj = (uj > vj) ? uj : vj;
 602                 ale_pos sumi = ui + vi;
 603                 ale_pos sumj = uj + vj;
 604
 605                 d[0] = (maxi + sumi) / 4;
 606                 d[1] = (maxj + sumj) / 4;
 607
 608                 if (scale_factor != 1) {
 609                         d[0] /= scale_factor;
 610                         d[1] /= scale_factor;
 611                         (*q)[0] /= scale_factor;
 612                         (*q)[1] /= scale_factor;
 613                 }
 614         }
 615
 616         /*
 617          * Modify all projective parameters at once.  Accommodate bugs in the
 618          * version 0 transformation file handler (ALE versions 0.4.0p1 and
 619          * earlier).  This code is only called when using a transformation data
 620          * file created with an old version of ALE.
 621          */
 622         virtual void gpt_v0_set(point x[4]) = 0;
 623
 624         /*
 625          * Modify all euclidean parameters at once.  Accommodate bugs in the
 626          * version 0 transformation file handler (ALE versions 0.4.0p1 and
 627          * earlier).  This code is only called when using a transformation data
 628          * file created with an old version of ALE.
 629          */
 630         virtual void eu_v0_set(ale_pos eu[3]) = 0;
 631
 632         virtual void debug_output() = 0;
 633
 634         virtual ~trans_abstract() {
 635         }
 636 };
 637
 638 #endif