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                 return p;
 167         }
 168
 169         /*
 170          * Barrel distortion inverse.
 171          */
 172         struct point bdi(struct point p) const {
 173                 if (bdcnum > 0) {
 174                         point half_diag = point(unscaled_height(), unscaled_width()) / 2;
 175
 176                         p -= half_diag;
 177
 178                         ale_pos r = p.norm() / half_diag.norm();
 179                         ale_pos s = r;
 180
 181                         while (fabs(r - bdr(s)) * half_diag.norm() > BARREL_INV_ERROR)
 182                                 s += (r - bdr(s)) / bdrd(s);
 183
 184                         if (r > 0.0001)
 185                                 p *= s / r;
 186
 187                         p += half_diag;
 188                 }
 189
 190                 assert (!isnan(p[0]) && !isnan(p[1]));
 191
 192                 return p;
 193         }
 194
 195         /*
 196          * Transformation sans barrel distortion
 197          */
 198         virtual struct point pe(struct point p) const = 0;
 199
 200         /*
 201          * Transformation inverse sans barrel distortion
 202          */
 203         virtual struct point pei(struct point p) const = 0;
 204
 205         /*
 206          * Map unscaled point p.
 207          */
 208         struct point transform_unscaled(struct point p) const {
 209                 return pe(bdi(p));
 210         }
 211
 212         /*
 213          * Transform point p.
 214          *
 215          * Barrel distortion correction followed by a projective/euclidean
 216          * transformation.
 217          */
 218         struct point transform_scaled(struct point p) const {
 219                 return transform_unscaled(p / scale_factor);
 220         }
 221
 222 #if 0
 223         /*
 224          * operator() is the transformation operator.
 225          */
 226         struct point operator()(struct point p) {
 227                 return transform(p);
 228         }
 229 #endif
 230
 231         /*
 232          * Map point p using the inverse of the transform into
 233          * the unscaled image space.
 234          */
 235         struct point unscaled_inverse_transform(struct point p) const {
 236                 return bd(pei(p));
 237         }
 238
 239         /*
 240          * Map point p using the inverse of the transform.
 241          *
 242          * Projective/euclidean inverse followed by barrel distortion.
 243          */
 244         struct point scaled_inverse_transform(struct point p) const {
 245                 assert (p.defined());
 246                 point q = unscaled_inverse_transform(p);
 247
 248                 q[0] *= scale_factor;
 249                 q[1] *= scale_factor;
 250
 251                 return q;
 252         }
 253
 254         /*
 255          * Calculate projective transformation parameters from a euclidean
 256          * transformation.
 257          */
 258         virtual void eu_to_gpt() = 0;
 259
 260         /*
 261          * Modify a euclidean transform in the indicated manner.
 262          */
 263         virtual void eu_modify(int i1, ale_pos diff) = 0;
 264
 265         /*
 266          * Rotate about a given point in the original reference frame.
 267          */
 268         virtual void eu_rotate_about_scaled(point center, ale_pos diff) = 0;
 269
 270         /*
 271          * Modify all euclidean parameters at once.
 272          */
 273         virtual void eu_set(ale_pos eu[3]) = 0;
 274
 275         /*
 276          * Get the specified euclidean parameter
 277          */
 278         virtual ale_pos eu_get(int param) const = 0;
 279
 280         /*
 281          * Modify a projective transform in the indicated manner.
 282          */
 283         virtual void gpt_modify(int i1, int i2, ale_pos diff) = 0;
 284
 285         /*
 286          * Modify a projective transform according to the group operation.
 287          */
 288         virtual void gr_modify(int i1, int i2, ale_pos diff) = 0;
 289
 290         /*
 291          * Modify all projective parameters at once.
 292          */
 293         virtual void gpt_set(point x[4]) = 0;
 294
 295         virtual void gpt_set(point x1, point x2, point x3, point x4) = 0;
 296
 297         /*
 298          * Snap positional parameters to the specified resolution.
 299          */
 300
 301         virtual void snap(ale_pos interval) = 0;
 302
 303         /*
 304          * Get the specified projective parameter
 305          */
 306         virtual point gpt_get(int point) const = 0;
 307
 308         /*
 309          * Get the specified projective parameter
 310          */
 311         virtual ale_pos gpt_get(int point, int dim) = 0;
 312
 313         /*
 314          * Check equality of transformation parameters.
 315          */
 316         virtual int operator==(const trans_abstract &t) const {
 317                 /*
 318                  * Small tolerances (< 10^-6?) can cause odd errors,
 319                  * possibly due to float<->double conversion issues.
 320                  */
 321                 double zero_tolerance = 0.01;
 322
 323                 if (scale() != t.scale())
 324                         return 0;
 325
 326                 if (is_projective() != t.is_projective())
 327                         return 0;
 328
 329                 if (is_projective()) {
 330                         assert (t.is_projective());
 331                         for (int i = 0; i < 4; i++)
 332                         for (int d = 0; d < 2; d++) {
 333                                 double abs_difference = fabs(gpt_get(i)[d] - t.gpt_get(i)[d]);
 334
 335                                 if (abs_difference > zero_tolerance)
 336                                         return 0;
 337                         }
 338                 } else {
 339                         assert (!t.is_projective());
 340                         for (int i = 0; i < 3; i++) {
 341                                 double abs_difference = fabs(eu_get(i) - t.eu_get(i));
 342
 343                                 if (abs_difference > zero_tolerance)
 344                                         return 0;
 345                         }
 346                 }
 347
 348                 return 1;
 349         }
 350
 351         virtual int operator!=(const trans_abstract &t) const {
 352                 return !(operator==(t));
 353         }
 354
 355
 356         /*
 357          * Translate by a given amount
 358          */
 359         virtual void translate(point p) = 0;
 360
 361         /*
 362          * Rotate by a given amount about a given point.
 363          */
 364         virtual void rotate(point p, ale_pos degrees) = 0;
 365
 366         /*
 367          * Set the specified barrel distortion parameter.
 368          */
 369         void bd_set(unsigned int degree, ale_pos value) {
 370                 assert (degree < bdcnum);
 371                 bdc[degree] = value;
 372         }
 373
 374         /*
 375          * Set all barrel distortion parameters.
 376          */
 377         void bd_set(unsigned int degree, ale_pos values[BARREL_DEGREE]) {
 378                 assert (degree <= BARREL_DEGREE);
 379                 bdcnum = degree;
 380                 for (unsigned int d = 0; d < degree; d++)
 381                         bdc[d] = values[d];
 382         }
 383
 384         /*
 385          * Get all barrel distortion parameters.
 386          */
 387         void bd_get(ale_pos result[BARREL_DEGREE]) {
 388                 for (unsigned int d = 0; d < bdcnum; d++)
 389                         result[d] = bdc[d];
 390         }
 391
 392         /*
 393          * Get the specified barrel distortion parameter.
 394          */
 395         ale_pos bd_get(unsigned int degree) {
 396                 assert (degree < bdcnum);
 397                 return bdc[degree];
 398         }
 399
 400         /*
 401          * Get the number of barrel distortion parameters.
 402          */
 403         unsigned int bd_count() {
 404                 return bdcnum;
 405         }
 406
 407         /*
 408          * Get the maximum allowable number of barrel distortion parameters.
 409          */
 410         unsigned int bd_max() {
 411                 return BARREL_DEGREE;
 412         }
 413
 414         /*
 415          * Modify the specified barrel distortion parameter.
 416          */
 417         void bd_modify(unsigned int degree, ale_pos diff) {
 418                 assert (degree < bdcnum);
 419                 bd_set(degree, bd_get(degree) + diff);
 420         }
 421
 422         /*
 423          * Rescale a transform with a given factor.
 424          */
 425         void rescale(ale_pos factor) {
 426                 specific_rescale(factor);
 427                 scale_factor *= factor;
 428         }
 429
 430         /*
 431          * Set a new domain.
 432          */
 433
 434         void set_domain(unsigned int new_height, unsigned int new_width) {
 435                 reset_memos();
 436                 input_width = new_width;
 437                 input_height = new_height;
 438         }
 439
 440         /*
 441          * Set the dimensions of the image.
 442          */
 443         void set_dimensions(const image *im) {
 444
 445                 int new_height = (int) im->height();
 446                 int new_width  = (int) im->width();
 447
 448                 reset_memos();
 449                 specific_set_dimensions(im);
 450                 input_height = new_height;
 451                 input_width  = new_width;
 452         }
 453
 454         /*
 455          * Get the position and dimensions of a pixel P mapped from one
 456          * coordinate system to another, using the forward transformation.
 457          * This function uses scaled input coordinates.
 458          */
 459         void map_area(point p, point *q, ale_pos d[2]) {
 460
 461                 /*
 462                  * Determine the coordinates in the target frame for the source
 463                  * image pixel P and two adjacent source pixels.
 464                  */
 465
 466                     (*q) = transform_scaled(p);
 467                 point q0 = transform_scaled(point(p[0] + 1, p[1]));
 468                 point q1 = transform_scaled(point(p[0], p[1] + 1));
 469
 470                 /*
 471                  * Calculate the distance between source image pixel and
 472                  * adjacent source pixels, measured in the coordinate system of
 473                  * the target frame.
 474                  */
 475
 476                 ale_pos ui = fabs(q0[0] - (*q)[0]);
 477                 ale_pos uj = fabs(q0[1] - (*q)[1]);
 478                 ale_pos vi = fabs(q1[0] - (*q)[0]);
 479                 ale_pos vj = fabs(q1[1] - (*q)[1]);
 480
 481                 /*
 482                  * We map the area of the source image pixel P onto the target
 483                  * frame as a rectangular area oriented on the target frame's
 484                  * axes.  Note that this results in an area that may be the
 485                  * wrong shape or orientation.
 486                  *
 487                  * We define two estimates of the rectangle's dimensions below.
 488                  * For rotations of 0, 90, 180, or 270 degrees, max and sum are
 489                  * identical.  For other orientations, sum is too large and max
 490                  * is too small.  We use the mean of max and sum, which we then
 491                  * divide by two to obtain the distance between the center and
 492                  * the edge.
 493                  */
 494
 495                 ale_pos maxi = (ui > vi) ? ui : vi;
 496                 ale_pos maxj = (uj > vj) ? uj : vj;
 497                 ale_pos sumi = ui + vi;
 498                 ale_pos sumj = uj + vj;
 499
 500                 d[0] = (maxi + sumi) / 4;
 501                 d[1] = (maxj + sumj) / 4;
 502         }
 503
 504         /*
 505          * Get the position and dimensions of a pixel P mapped from one
 506          * coordinate system to another, using the forward transformation.
 507          * This function uses unscaled input coordinates.
 508          */
 509         void map_area_unscaled(point p, point *q, ale_pos d[2]) {
 510
 511                 /*
 512                  * Determine the coordinates in the target frame for the source
 513                  * image pixel P and two adjacent source pixels.
 514                  */
 515
 516                     (*q) = transform_unscaled(p);
 517                 point q0 = transform_unscaled(point(p[0] + 1, p[1]));
 518                 point q1 = transform_unscaled(point(p[0], p[1] + 1));
 519
 520                 /*
 521                  * Calculate the distance between source image pixel and
 522                  * adjacent source pixels, measured in the coordinate system of
 523                  * the target frame.
 524                  */
 525
 526                 ale_pos ui = fabs(q0[0] - (*q)[0]);
 527                 ale_pos uj = fabs(q0[1] - (*q)[1]);
 528                 ale_pos vi = fabs(q1[0] - (*q)[0]);
 529                 ale_pos vj = fabs(q1[1] - (*q)[1]);
 530
 531                 /*
 532                  * We map the area of the source image pixel P onto the target
 533                  * frame as a rectangular area oriented on the target frame's
 534                  * axes.  Note that this results in an area that may be the
 535                  * wrong shape or orientation.
 536                  *
 537                  * We define two estimates of the rectangle's dimensions below.
 538                  * For rotations of 0, 90, 180, or 270 degrees, max and sum are
 539                  * identical.  For other orientations, sum is too large and max
 540                  * is too small.  We use the mean of max and sum, which we then
 541                  * divide by two to obtain the distance between the center and
 542                  * the edge.
 543                  */
 544
 545                 ale_pos maxi = (ui > vi) ? ui : vi;
 546                 ale_pos maxj = (uj > vj) ? uj : vj;
 547                 ale_pos sumi = ui + vi;
 548                 ale_pos sumj = uj + vj;
 549
 550                 d[0] = (maxi + sumi) / 4;
 551                 d[1] = (maxj + sumj) / 4;
 552         }
 553
 554         /*
 555          * Get the position and dimensions of a pixel P mapped from one
 556          * coordinate system to another, using the inverse transformation.  If
 557          * SCALE_FACTOR is not equal to one, divide out the scale factor to
 558          * obtain unscaled coordinates.  This method is very similar to the
 559          * map_area method above.
 560          */
 561         void unscaled_map_area_inverse(point p, point *q, ale_pos d[2]) {
 562
 563                 /*
 564                  * Determine the coordinates in the target frame for the source
 565                  * image pixel P and two adjacent source pixels.
 566                  */
 567
 568                     (*q) = scaled_inverse_transform(p);
 569                 point q0 = scaled_inverse_transform(point(p[0] + 1, p[1]));
 570                 point q1 = scaled_inverse_transform(point(p[0], p[1] + 1));
 571
 572
 573                 /*
 574                  * Calculate the distance between source image pixel and
 575                  * adjacent source pixels, measured in the coordinate system of
 576                  * the target frame.
 577                  */
 578
 579                 ale_pos ui = fabs(q0[0] - (*q)[0]);
 580                 ale_pos uj = fabs(q0[1] - (*q)[1]);
 581                 ale_pos vi = fabs(q1[0] - (*q)[0]);
 582                 ale_pos vj = fabs(q1[1] - (*q)[1]);
 583
 584                 /*
 585                  * We map the area of the source image pixel P onto the target
 586                  * frame as a rectangular area oriented on the target frame's
 587                  * axes.  Note that this results in an area that may be the
 588                  * wrong shape or orientation.
 589                  *
 590                  * We define two estimates of the rectangle's dimensions below.
 591                  * For rotations of 0, 90, 180, or 270 degrees, max and sum are
 592                  * identical.  For other orientations, sum is too large and max
 593                  * is too small.  We use the mean of max and sum, which we then
 594                  * divide by two to obtain the distance between the center and
 595                  * the edge.
 596                  */
 597
 598                 ale_pos maxi = (ui > vi) ? ui : vi;
 599                 ale_pos maxj = (uj > vj) ? uj : vj;
 600                 ale_pos sumi = ui + vi;
 601                 ale_pos sumj = uj + vj;
 602
 603                 d[0] = (maxi + sumi) / 4;
 604                 d[1] = (maxj + sumj) / 4;
 605
 606                 if (scale_factor != 1) {
 607                         d[0] /= scale_factor;
 608                         d[1] /= scale_factor;
 609                         (*q)[0] /= scale_factor;
 610                         (*q)[1] /= scale_factor;
 611                 }
 612         }
 613
 614         /*
 615          * Modify all projective parameters at once.  Accommodate bugs in the
 616          * version 0 transformation file handler (ALE versions 0.4.0p1 and
 617          * earlier).  This code is only called when using a transformation data
 618          * file created with an old version of ALE.
 619          */
 620         virtual void gpt_v0_set(point x[4]) = 0;
 621
 622         /*
 623          * Modify all euclidean parameters at once.  Accommodate bugs in the
 624          * version 0 transformation file handler (ALE versions 0.4.0p1 and
 625          * earlier).  This code is only called when using a transformation data
 626          * file created with an old version of ALE.
 627          */
 628         virtual void eu_v0_set(ale_pos eu[3]) = 0;
 629
 630         virtual void debug_output() = 0;
 631
 632         virtual ~trans_abstract() {
 633         }
 634 };
 635
 636 #endif