Changeset 11616 for ImageMagick/trunk

Show
Ignore:
Timestamp:
07/13/08 06:00:38 (7 weeks ago)
Author:
anthony
Message:
 
Files:
1 modified

Legend:

Unmodified
Added
Removed

  • ImageMagick/trunk/magick/distort.c

    r11612 r11616  
    197197} 
    198198 
     199#if 0   
     200   For some reason compiler errors on % in this! 
     201static inline double MagickFraction(double x) 
     202{ 
     203  /* value mod 1 in a  -.5 to .5 range */ 
     204  return ((((x%1.0)+1.5)%1.0)-0.5); 
     205} 
     206#endif 
    199207static inline double MagickRound(double x) 
    200208{ 
    201   assert(x >= (1.0*LONG_MIN-0.5)); 
    202   assert(x <= (1.0*LONG_MAX+0.5)); 
     209  /* round the fraction to nearest integer */ 
    203210  if (x >= 0.0) 
    204211    return((double) ((long) (x+0.5))); 
    205212  return((double) ((long) (x-0.5))); 
     213#if 0 
     214  return(x - MagickFraction(x)); 
     215#endif 
    206216} 
    207217 
     218static unsigned long poly_number_terms(double order) 
     219{ 
     220  /* Return the number of terms for a 2d polynomial 
     221     Order must either be an integer, or 1.5 to produce 
     222     the 2 dimentional polyminal function... 
     223        affine     1 (3)      u = c0 + c1*x + c2*y 
     224        bilinear  1.5 (4)     u = '' + c3*x*y 
     225        quadratic  2 (6)      u = '' + c4*x*x + c5*y*y 
     226        cubic      3 (10)     u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3 
     227        quartic    4 (15)     u = '' + c10*x^4 + ... + c14*y^4 
     228        quintic    5 (21)     u = '' + c15*x^5 + ... + c20*y^5 
     229     number in parenthesis is minimum number of points needed. 
     230     Anything beyond quintic, has not been implemented until 
     231     a more automated way of determined terms is found. 
     232   */ 
     233  if ( order < 1 || order > 5 || 
     234       ( order != floor(order) && (order-1.5) > MagickEpsilon) ) 
     235    return 0; /* invalid polynomial order */ 
     236  return ( (long)floor((order+1)*(order+2)/2) ); 
     237} 
     238 
     239static double poly_term(long n, double x, double y) 
     240{ 
     241  /* return the result for this polynomial term */ 
     242  switch(n) { 
     243    case  0:  return( 1.0 ); /* constant */ 
     244    case  1:  return(  x  ); 
     245    case  2:  return(  y  ); /* affine      order = 1   terms = 3 */ 
     246    case  3:  return( x*y ); /* bilinear    order = 1.5 terms = 4 */ 
     247    case  4:  return( x*x ); 
     248    case  5:  return( y*y ); /* quadratic   order = 2   terms = 6 */ 
     249    case  6:  return( x*x*x ); 
     250    case  7:  return( x*x*y ); 
     251    case  8:  return( x*y*y ); 
     252    case  9:  return( y*y*y ); /* cubic       order = 3   terms = 10 */ 
     253    case 10:  return( x*x*x*x ); 
     254    case 11:  return( x*x*x*y ); 
     255    case 12:  return( x*x*y*y ); 
     256    case 13:  return( x*y*y*y ); 
     257    case 14:  return( y*y*y*y ); /* quartic     order = 4   terms = 15 */ 
     258    case 15:  return( x*x*x*x*x ); 
     259    case 16:  return( x*x*x*x*y ); 
     260    case 17:  return( x*x*x*y*y ); 
     261    case 18:  return( x*x*y*y*y ); 
     262    case 19:  return( x*y*y*y*y ); 
     263    case 20:  return( y*y*y*y*y ); /* quintic     order = 5   terms = 21 */ 
     264  } 
     265  return( 0 ); /* should never happen */ 
     266} 
     267static const char *poly_term_str(long n) 
     268{ 
     269  /* return the result for this polynomial term */ 
     270  switch(n) { 
     271    case  0:  return(""); /* constant */ 
     272    case  1:  return("*xx"); 
     273    case  2:  return("*yy"); /* affine      order = 1   terms = 3 */ 
     274    case  3:  return("*xx*yy"); /* bilinear    order = 1.5 terms = 4 */ 
     275    case  4:  return("*xx*xx"); 
     276    case  5:  return("*yy*yy"); /* quadratic   order = 2   terms = 6 */ 
     277    case  6:  return("*xx*xx*xx"); 
     278    case  7:  return("*xx*xx*yy"); 
     279    case  8:  return("*xx*yy*yy"); 
     280    case  9:  return("*yy*yy*yy"); /* cubic       order = 3   terms = 10 */ 
     281    case 10:  return("*xx*xx*xx*xx"); 
     282    case 11:  return("*xx*xx*xx*yy"); 
     283    case 12:  return("*xx*xx*yy*yy"); 
     284    case 13:  return("*xx*yy*yy*yy"); 
     285    case 14:  return("*yy*yy*yy*yy"); /* quartic     order = 4   terms = 15 */ 
     286    case 15:  return("*xx*xx*xx*xx*xx"); 
     287    case 16:  return("*xx*xx*xx*xx*yy"); 
     288    case 17:  return("*xx*xx*xx*yy*yy"); 
     289    case 18:  return("*xx*xx*yy*yy*yy"); 
     290    case 19:  return("*xx*yy*yy*yy*yy"); 
     291    case 20:  return("*yy*yy*yy*yy*yy"); /* quintic     order = 5   terms = 21 */ 
     292  } 
     293  return( "UNKNOWN" ); /* should never happen */ 
     294} 
     295static double poly_term_dx(long n, double x, double y) 
     296{ 
     297  /* polynomial term for x derivative */ 
     298  switch(n) { 
     299    case  0:  return( 0.0 ); /* constant */ 
     300    case  1:  return( 1.0 ); 
     301    case  2:  return( 0.0 ); /* affine      order = 1   terms = 3 */ 
     302    case  3:  return(  y  ); /* bilinear    order = 1.5 terms = 4 */ 
     303    case  4:  return(  x  ); 
     304    case  5:  return( 0.0 ); /* quadratic   order = 2   terms = 6 */ 
     305    case  6:  return( x*x ); 
     306    case  7:  return( x*y ); 
     307    case  8:  return( y*y ); 
     308    case  9:  return( 0.0 ); /* cubic       order = 3   terms = 10 */ 
     309    case 10:  return( x*x*x ); 
     310    case 11:  return( x*x*y ); 
     311    case 12:  return( x*y*y ); 
     312    case 13:  return( y*y*y ); 
     313    case 14:  return( 0.0 ); /* quartic     order = 4   terms = 15 */ 
     314    case 15:  return( x*x*x*x ); 
     315    case 16:  return( x*x*x*y ); 
     316    case 17:  return( x*x*y*y ); 
     317    case 18:  return( x*y*y*y ); 
     318    case 19:  return( y*y*y*y ); 
     319    case 20:  return( 0.0 ); /* quintic     order = 5   terms = 21 */ 
     320  } 
     321  return( 0.0 ); /* should never happen */ 
     322} 
     323static double poly_term_dy(long n, double x, double y) 
     324{ 
     325  /* polynomial term for y derivative */ 
     326  switch(n) { 
     327    case  0:  return( 0.0 ); /* constant */ 
     328    case  1:  return( 0.0 ); 
     329    case  2:  return( 1.0 ); /* affine      order = 1   terms = 3 */ 
     330    case  3:  return(  x  ); /* bilinear    order = 1.5 terms = 4 */ 
     331    case  4:  return( 0.0 ); 
     332    case  5:  return(  y  ); /* quadratic   order = 2   terms = 6 */ 
     333    default:  return( poly_term_dx(n-1,x,y) ); /* weird but true */ 
     334  } 
     335  /* NOTE: the only reason that last is not true for 'quadtratic' 
     336     is due to the re-arrangement of terms to allow for 'bilinear' 
     337  */ 
     338} 
    208339 
    209340static double *DistortCoefficents(Image *image, DistortImageMethod *method, 
    … …  
    230361    i; 
    231362 
    232   MagickBooleanType 
    233     status; 
    234  
    235   /* currently the array is always 9 doubles, that may change in the future */ 
    236   coefficients = AcquireQuantumMemory(9,sizeof(*coefficients)); 
    237   if (coefficients == (double *) NULL) { 
    238     (void) ThrowMagickException(exception,GetMagickModule(), 
    239                   ResourceLimitError,"MemoryAllocationFailed", 
    240                   "%s", "DistortCoefficients"); 
    241     return((double *) NULL); 
    242   } 
    243   /* zero out coeffiecents array */ 
    244   for (i=0; i < 9L; i++) 
    245     coefficients[i] = 0.0; 
     363  unsigned long 
     364    number_coefficients;   /* number of coefficients needed */ 
    246365 
    247366  /* If not enough control point pairs are found for specific distortions 
    … …  
    256375  switch (*method) 
    257376  { 
     377    case ScaleRotateTranslateDistortion: 
     378    case AffineDistortion: 
     379    case AffineProjectionDistortion: 
     380      number_coefficients=6; 
     381      break; 
     382    case PerspectiveDistortion: 
     383    case PerspectiveProjectionDistortion: 
     384      number_coefficients=9; 
     385      break; 
     386    case BilinearDistortion: 
     387      number_coefficients=8; 
     388      break; 
     389    case PolynomialDistortion: 
     390      if ( number_arguments < 5 && (number_arguments-1)%4 != 0) 
     391      { 
     392        (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 
     393                     "InvalidArgument","%s : '%s'","distort Polynomial", 
     394                     "Invalid number of coord-pairs: order [sX,sY,dX,dY]*"); 
     395        return((double *) NULL); 
     396      } 
     397      number_coefficients = poly_number_terms(arguments[0]); 
     398      if ( number_coefficients == 0 ) 
     399      { 
     400        (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 
     401                     "InvalidArgument","%s : '%s'","distort Polynomial", 
     402                     "Invalid order, should be 1,1.5,2,3,4 or 5"); 
     403        return((double *) NULL); 
     404      } 
     405      if ( number_arguments < number_coefficients*4+1) 
     406      { 
     407        (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 
     408                     "InvalidArgument", "%s : '%s %ld'","distort Polynomial", 
     409                     "Not enough coord-pairs, minimum", number_coefficients); 
     410        return((double *) NULL); 
     411      } 
     412      number_coefficients = number_coefficients*2+1; 
     413      break; 
     414    case ArcDistortion: 
     415      number_coefficients=5; 
     416      break; 
     417    case UndefinedDistortion: 
     418    default: 
     419      number_coefficients=0; /* should never happen */ 
     420      break; 
     421  } 
     422 
     423  /* allocate the array of coefficients needed */ 
     424  coefficients = AcquireQuantumMemory(number_coefficients, 
     425                        sizeof(*coefficients)); 
     426  if (coefficients == (double *) NULL) { 
     427    (void) ThrowMagickException(exception,GetMagickModule(), 
     428                  ResourceLimitError,"MemoryAllocationFailed", 
     429                  "%s", "DistortCoefficients"); 
     430    return((double *) NULL); 
     431  } 
     432  /* zero out coeffiecents array */ 
     433  for (i=0; i < (long)number_coefficients; i++) 
     434    coefficients[i] = 0.0; 
     435 
     436  switch (*method) 
     437  { 
    258438    case AffineDistortion: 
    259439    { 
    260       /* Input Arguments are any number of pairs of distorted coodinates 
     440      /* Affine Distortion 
     441           u =  c0*x + c2*y + c4 
     442           v =  c1*x + c3*y + c5 
     443 
     444         Input Arguments are any number of pairs of distorted coodinates 
    261445         Which will be used to generate coefficients for Distortion. 
    262446            u1,v1, x1,y1,  u2,v2, x2,y2,  u3,v3, x3,y3,  u4,v4, x4,y4 ... 
    … …  
    289473      } 
    290474      else { 
    291         /* 
    292            Affine Distortion 
    293               u =  c0*x + c2*y + c4 
    294               v =  c1*x + c3*y + c5 
     475        /* 2 or more points (usally 3) given. 
     476           Solve a least squares simultanious equation for coefficients. 
    295477        */ 
    296478        double 
    … …  
    298480          **vectors, 
    299481          terms[3]; 
     482 
     483        MagickBooleanType 
     484          status; 
    300485 
    301486        *method = AffineDistortion; 
    … …  
    321506        if ( number_arguments == 8 ) { 
    322507          /* Only two pairs were given, but we need 3 to solve the affine. 
    323              Fake a 3rd control point pair by p1 around p0 by 90 degrees 
     508             Fake a pair by rotating p1 around p0 by 90 degrees. 
    324509             That is given    u0,v0 x0,y0   u1,v1 x1,y1 
    325510             Then   u2 = u0 - (v1-v0)   v2 = v0 + (u1-u0) 
    … …  
    330515 
    331516          uv2[0] = arguments[0] - arguments[5] + arguments[1];   /* u2 */ 
    332           uv2[1] = arguments[1] + arguments[4] + arguments[0];   /* v2 */ 
     517          uv2[1] = arguments[1] + arguments[4] - arguments[0];   /* v2 */ 
    333518          terms[0] = arguments[2] - arguments[7] + arguments[3]; /* x2 */ 
    334           terms[1] = arguments[3] + arguments[6] + arguments[2]; /* y2 */ 
     519          terms[1] = arguments[3] + arguments[6] - arguments[2]; /* y2 */ 
    335520          terms[2] = 1;                                          /* 1 */ 
    336521          LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL); 
    337522        } 
     523      /* Solve for LeastSquares Coefficients */ 
    338524        status=GaussJordanElimination(matrix,vectors,3UL,2UL); 
    339525        if ( status == MagickTrue ) { 
    … …  
    373559        } 
    374560        InvertAffineCoefficients(coefficients, inverse); 
    375         fprintf(stderr, "Affine ReverseMap\n"); 
    376         for (i=0; i<6; i++) 
    377           fprintf(stderr, "  %lf\n", coefficients[i]); 
    378         fprintf(stderr, "Affine Projection\n"); 
    379         for (i=0; i<6; i++) 
    380           fprintf(stderr, "  %lf\n", inverse[i]); 
     561        fprintf(stderr, "Affine Reverse Map\n"); 
     562        fprintf(stderr, "  -fx 'xx=%+lf*i %+lf*j %+lf;\n", 
     563            coefficients[0], coefficients[2], coefficients[4]); 
     564        fprintf(stderr, "       yy=%+lf*i %+lf*j %+lf; p{xx,yy}'\n", 
     565            coefficients[1], coefficients[3], coefficients[5]); 
     566        fprintf(stderr, "Affine Forward Map\n"); 
     567        fprintf(stderr, "  -distort AffineProjection \\\n  '"); 
     568        for (i=0; i<5; i++) 
     569          fprintf(stderr, "%lg,", inverse[i]); 
     570        fprintf(stderr, "%lg'\n", inverse[5]); 
    381571        inverse = RelinquishMagickMemory(inverse); 
    382572      } 
    … …  
    400590      InvertAffineCoefficients(arguments, coefficients); 
    401591      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) { 
    402         fprintf(stderr, "Affine ReverseMap\n"); 
    403         for (i=0; i<6; i++) 
    404           fprintf(stderr, "  %lf\n", coefficients[i]); 
    405         fprintf(stderr, "AffineProjection\n"); 
    406         for (i=0; i<6; i++) 
    407           fprintf(stderr, "  %lf\n", arguments[i]); 
     592        fprintf(stderr, "Affine Reverse Map\n"); 
     593        fprintf(stderr, "  -fx 'xx=%+lf*i %+lf*j %+lf;\n", 
     594            coefficients[0], coefficients[2], coefficients[4]); 
     595        fprintf(stderr, "       yy=%+lf*i %+lf*j %+lf; p{xx,yy}'\n", 
     596            coefficients[1], coefficients[3], coefficients[5]); 
    408597      } 
    409598      return(coefficients); 
    … …  
    434623        *vectors[1], 
    435624        terms[8]; 
     625 
     626      MagickBooleanType 
     627        status; 
    436628 
    437629      *method = PerspectiveDistortion; 
    … …  
    470662            8UL,1UL); 
    471663      } 
     664      /* Solve for LeastSquares Coefficients */ 
    472665      status=GaussJordanElimination(matrix,vectors,8UL,1UL); 
    473666      matrix = RelinquishMagickMatrix(matrix, 8UL); 
    … …  
    499692        } 
    500693        InvertPerspectiveCoefficients(coefficients, inverse); 
    501         fprintf(stderr, "Perspective ReverseMap\n"); 
    502         for (i=0; i<8; i++) 
    503           fprintf(stderr, "  %lf\n", coefficients[i]); 
    504         fprintf(stderr, "PerspectiveProjection\n"); 
    505         for (i=0; i<8; i++) 
    506           fprintf(stderr, "  %lf\n", inverse[i]); 
     694        fprintf(stderr, "Perspective Reverse Map\n"); 
     695        fprintf(stderr, "  -fx 'xx=%+lf*i %+lf*j %+lf;\n", 
     696            coefficients[0], coefficients[1], coefficients[2]); 
     697        fprintf(stderr, "       yy=%+lf*i %+lf*j %+lf;\n", 
     698            coefficients[3], coefficients[4], coefficients[5]); 
     699        fprintf(stderr, "       rr=%+lf*i %+lf*j + 1; p{xx/rr,yy/rr}'\n", 
     700            coefficients[6], coefficients[7]); 
     701        fprintf(stderr, "Perspective Forward Map\n"); 
     702        fprintf(stderr, "  -distort PerspectiveProjection \\\n      '"); 
     703        for (i=0; i<4; i++) 
     704          fprintf(stderr, "%lg,", inverse[i]); 
     705        fprintf(stderr, "\n       "); 
     706        for (; i<7; i++) 
     707          fprintf(stderr, "%lg,", inverse[i]); 
     708        fprintf(stderr, "%lg'\n", inverse[7]); 
    507709        inverse = RelinquishMagickMemory(inverse); 
    508710      } 
    … …  
    522724      InvertPerspectiveCoefficients(arguments, coefficients); 
    523725      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) { 
    524         fprintf(stderr, "Perspective ReverseMap\n"); 
    525         for (i=0; i<8; i++) 
    526           fprintf(stderr, "  %lf\n", coefficients[i]); 
    527         fprintf(stderr, "PerspectiveProjection\n"); 
    528         for (i=0; i<8; i++) 
    529           fprintf(stderr, "  %lf\n", arguments[i]); 
     726        fprintf(stderr, "Perspective Reverse Map\n"); 
     727        fprintf(stderr, "  -fx 'xx=%+lf*i %+lf*j %+lf;\n", 
     728            coefficients[0], coefficients[1], coefficients[2]); 
     729        fprintf(stderr, "       yy=%+lf*i %+lf*j %+lf;\n", 
     730            coefficients[3], coefficients[4], coefficients[5]); 
     731        fprintf(stderr, "       rr=%+lf*i %+lf*j + 1; p{xx/rr,yy/rr}'\n", 
     732            coefficients[6], coefficients[7]); 
    530733      } 
    531734      /* 
    … …  
    548751        terms[4]; 
    549752 
     753      MagickBooleanType 
     754        status; 
     755 
    550756      /* Bilinear Distortion (currently reversed -- this needs to be fixed) 
    551757            u = c0*x + c1*y + c2*x*y + c3; 
    … …  
    582788        LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i]),4UL,2UL); 
    583789      } 
     790      /* Solve for LeastSquares Coefficients */ 
    584791      status=GaussJordanElimination(matrix,vectors,4UL,2UL); 
    585792      matrix = RelinquishMagickMatrix(matrix, 4UL); 
    … …  
    591798        return((double *) NULL); 
    592799      } 
     800      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) { 
     801        fprintf(stderr, "Bilinear Reverse Map\n"); 
     802        fprintf(stderr, "  -fx 'xx=%+lf*i %+lf*j %+lf*i*j %+lf;\n", 
     803            coefficients[0], coefficients[1], coefficients[2], coefficients[3]); 
     804        fprintf(stderr, "       yy=%+lf*i %+lf*j %+lf*i*j %+lf;\n", 
     805            coefficients[4], coefficients[5], coefficients[6], coefficients[7]); 
     806        fprintf(stderr, "       p{xx,yy}'\n"); 
     807      } 
     808      return(coefficients); 
     809    } 
     810    case PolynomialDistortion: 
     811    { 
     812      /* Polynomial Distortion 
     813           c0 = number of terms in one polynomial equation 
     814 
     815            u = c1 + c1*x + c2*y + c3*x*y + ... + c8*x^3 + c9*y^3 
     816            v = .... 
     817 
     818         Input Arguments are pairs of distorted coodinates 
     819         Which will be used to generate the coefficients of the above. 
     820            u1,v1, x1,y1,  u2,v2, x2,y2,  u3,v3, x3,y3,  u4,v4, x4,y4 ... 
     821      */ 
     822      double 
     823        **matrix, 
     824        *vectors[2], 
     825        *terms; 
     826 
     827      long 
     828        nterms; 
     829 
     830      register long 
     831        j; 
     832 
     833      MagickBooleanType 
     834        status; 
     835 
     836      nterms = coefficients[0] = poly_number_terms(arguments[0]); 
     837 
     838      terms = AcquireQuantumMemory(nterms, sizeof(*terms)); 
     839      if (terms == (double *) NULL) { 
     840        coefficients = RelinquishMagickMemory(coefficients); 
     841        (void) ThrowMagickException(exception,GetMagickModule(), 
     842                  ResourceLimitError,"MemoryAllocationFailed", 
     843                  "%s", "DistortCoefficients"); 
     844        return((double *) NULL); 
     845      } 
     846      matrix = AcquireMagickMatrix(nterms,nterms); 
     847      if (matrix == (double **) NULL) { 
     848        terms = RelinquishMagickMemory(terms); 
     849        coefficients = RelinquishMagickMemory(coefficients); 
     850        (void) ThrowMagickException(exception,GetMagickModule(), 
     851                  ResourceLimitError,"MemoryAllocationFailed", 
     852                  "%s", "DistortCoefficients"); 
     853        return((double *) NULL); 
     854      } 
     855      /* fake a 2xN vectors matrix from rest of coefficients array */ 
     856      vectors[0] = &(coefficients[1]); 
     857      vectors[1] = &(coefficients[nterms+1]); 
     858      /* Add control points for least squares solving */ 
     859      for (i=2; i < (long) number_arguments; i+=4) { 
     860        for (j=0; j < nterms; j++) 
     861          terms[j] = poly_term(j,arguments[i],arguments[i+1]); 
     862        LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i]),nterms,2UL); 
     863      } 
     864      terms = RelinquishMagickMemory(terms); 
     865      /* Solve for LeastSquares Coefficients */ 
     866      status=GaussJordanElimination(matrix,vectors,nterms,2UL); 
     867      matrix = RelinquishMagickMatrix(matrix, nterms); 
     868      if ( status == MagickFalse ) { 
     869        coefficients = RelinquishMagickMemory(coefficients); 
     870        (void) ThrowMagickException(exception,GetMagickModule(),OptionError, 
     871                     "InvalidArgument","%s : '%s'","distort Polynomial", 
     872                     "Unsolvable Matrix"); 
     873        return((double *) NULL); 
     874      } 
     875      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) { 
     876        fprintf(stderr, "Polynomial (%ld terms) Reverse Map\n", (long)nterms); 
     877        fprintf(stderr, "  -fx 'xx="); 
     878        for (i=0; i<nterms; i++) { 
     879          if ( i != 0 && i%3 == 0 ) fprintf(stderr, "\n          "); 
     880          fprintf(stderr, " %lf%s", coefficients[i+1], poly_term_str(i)); 
     881        } 
     882        fprintf(stderr, "\n       yy="); 
     883        for (i=0; i<nterms; i++) { 
     884          if ( i != 0 && i%3 == 0 ) fprintf(stderr, "\n          "); 
     885          fprintf(stderr, " %lf%s", coefficients[i+1+nterms], poly_term_str(i)); 
     886        } 
     887        fprintf(stderr, "\n       p{xx,yy}'\n"); 
     888      } 
    593889      return(coefficients); 
    594890    } 
    … …  
    599895        x,y,sx,sy,a,nx,ny; 
    600896 
    601       /* 
     897      /* Scale, Rotate and Translate Distortion 
    602898         Argument options, by number of arguments given: 
    603899           7: x,y, sx,sy, a, nx,ny 
    … …  
    7111007          fprintf(stderr, "  %lf\n", inverse[i]); 
    7121008      } 
     1009      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) { 
     1010        double *inverse = AcquireQuantumMemory(8,sizeof(coefficients[i])); 
     1011        if (inverse == (double *) NULL) { 
     1012          coefficients = RelinquishMagickMemory(coefficients); 
     1013          (void) ThrowMagickException(exception,GetMagickModule(), 
     1014                  ResourceLimitError,"MemoryAllocationFailed", 
     1015                  "%s", "DistortCoefficients"); 
     1016          return((double *) NULL); 
     1017        } 
     1018        InvertAffineCoefficients(coefficients, inverse); 
     1019        fprintf(stderr, "Affine Reverse Map\n"); 
     1020        fprintf(stderr, "  -fx 'xx=%+lf*i %+lf*j %+lf;\n", 
     1021            coefficients[0], coefficients[2], coefficients[4]); 
     1022        fprintf(stderr, "       yy=%+lf*i %+lf*j %+lf; p{xx,yy}'\n", 
     1023            coefficients[1], coefficients[3], coefficients[5]); 
     1024        fprintf(stderr, "Affine Forward Map\n"); 
     1025        fprintf(stderr, "  -distort AffineProjection \\\n  '"); 
     1026        for (i=0; i<5; i++) 
     1027          fprintf(stderr, "%lg,", inverse[i]); 
     1028        fprintf(stderr, "%lg'\n", inverse[5]); 
     1029        inverse = RelinquishMagickMemory(inverse); 
     1030      } 
    7131031      return(coefficients); 
    7141032    } 
    … …  
    7771095      } 
    7781096      coefficients[4] = (1.0*image->columns-1.0)/2.0; 
     1097      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) { 
     1098        fprintf(stderr, "Arc Reverse Map\n"); 
     1099#if 0 
     1100  Not working yet 
     1101        fprintf(stderr, "c0=%-8lg  # angle for center line in source image\n", 
     1102             coefficients[0]); 
     1103        fprintf(stderr, "c1=%-8lg  # angle scale for mapping to source image\n", 
     1104             coefficients[1]); 
     1105        fprintf(stderr, "c2=%-8lg  # radius for top of source imagen\n", 
     1106             coefficients[2]); 
     1107        fprintf(stderr, "c3=%-8lg  # radius scale for mapping source image\n", 
     1108             coefficients[3]); 
     1109        fprintf(stderr, "c4=%-8lg  # centerline of arc within source image\n", 
     1110             coefficients[4]); 
     1111 
     1112        coefficients[1] = 2.0*MagickPI*image->columns/coefficients[1]; 
     1113        coefficients[3] = 1.0*image->rows/coefficients[3]; 
     1114          point.x = (atan2((double)y,(double)x) - coefficients[0])/(2*MagickPI); 
     1115          point.x -= MagickRound(point.x);     /* angle */ 
     1116          point.y = sqrt((double) x*x+y*y);    /* radius */ 
     1117          point.x = point.x*coefficients[1] + coefficients[4]; 
     1118          point.y = (coefficients[2] - point.y) * coefficients[3]; 
     1119 
     1120 
     1121        fprintf(stderr, "  -size %ldx%ld -page -%ld-%ld xc: +swap \\\n", 
     1122            (long)coefficients[2]*2+2, (long)coefficients[2]*2+2, 
     1123            (long)coefficients[2]+1, (long)coefficients[2]+1); 
     1124        fprintf(stderr, "  -fx 'xx=(atan2(j-h/2,i-w/2) %+lg)*2*pi;\n", 
     1125            coefficients[0]); 
     1126        fprintf(stderr, "       xx=(xx%%1+1.5)%%1-0.5;\n"); 
     1127        fprintf(stderr, "       yy=hypot(i-w/2,j-h/2);\n"); 
     1128        fprintf(stderr, "       xx=xx*%lg + %lg;\n", 
     1129            2.0*MagickPI*image->columns/coefficients[1], coefficients[4]); 
     1130        fprintf(stderr, "       yy=(%lg - yy) * %lg;\n", 
     1131            coefficients[2], ((double)image->rows)/coefficients[3] ); 
     1132        fprintf(stderr, "       v.p{xx,yy}'\n"); 
     1133#endif 
     1134      } 
    7791135      return(coefficients); 
    7801136    } 
    … …  
    8541210  point.y=0.0; 
    8551211 
     1212  /* These distotions must be bestfit, as display is warped toomuch */ 
     1213  if ( method == ArcDistortion ) 
     1214    bestfit = MagickTrue; 
     1215 
    8561216  /* Work out the 'best fit', (required for ArcDistortion) */ 
    857   if ( bestfit || method == ArcDistortion ) { 
     1217  if ( bestfit ) { 
    8581218    double 
    8591219      min_x,max_x,min_y,max_y; 
    … …  
    9701330      } 
    9711331      case BilinearDistortion: 
     1332      case PolynomialDistortion: 
    9721333      default: 
    9731334        /* no bestfit available for this distortion YET */ 
    … …  
    9851346  /* User provided a 'viewport' expert option which may 
    9861347     overrides some parts of the current output image geometry. 
    987      For ArcDistortion, this also overrides it default 'bestfit' setting. 
     1348     For ArcDistortion, this also overrides its default 'bestfit' setting. 
    9881349  */ 
    9891350  artifact = GetImageArtifact(image,"distort:viewport"); 
    9901351  if (artifact != (const char *) NULL) 
    9911352    (void) ParseAbsoluteGeometry(artifact,&geometry); 
    992  
    9931353 
    9941354  /* 
    … …  
    10931453          point.y=coefficients[1]*x+coefficients[3]*y+coefficients[5]; 
    10941454          /* Affine partial derivitives are constant -- set above */ 
    1095           break; 
    1096         } 
    1097         case BilinearDistortion: 
    1098         { 
    1099           point.x=coefficients[0]*x+coefficients[1]*y+coefficients[2]*x*y+ 
    1100             coefficients[3]; 
    1101           point.y=coefficients[4]*x+coefficients[5]*y+coefficients[6]*x*y+ 
    1102             coefficients[7]; 
    1103           /* Bilinear partial derivitives of scaling vectors */ 
    1104           ScaleResampleFilter( resample_filter[id], 
    1105               coefficients[0] + coefficients[2]*y, 
    1106               coefficients[1] + coefficients[2]*x, 
    1107               coefficients[4] + coefficients[6]*y, 
    1108               coefficients[5] + coefficients[6]*x ); 
    11091455          break; 
    11101456        } 
    … …  
    11451491          break; 
    11461492        } 
     1493        case BilinearDistortion: 
     1494        { 
     1495          point.x=coefficients[0]*x+coefficients[1]*y+coefficients[2]*x*y+ 
     1496            coefficients[3]; 
     1497          point.y=coefficients[4]*x+coefficients[5]*y+coefficients[6]*x*y+ 
     1498            coefficients[7]; 
     1499          /* Bilinear partial derivitives of scaling vectors */ 
     1500          ScaleResampleFilter( resample_filter[id], 
     1501              coefficients[0] + coefficients[2]*y, 
     1502              coefficients[1] + coefficients[2]*x, 
     1503              coefficients[4] + coefficients[6]*y, 
     1504              coefficients[5] + coefficients[6]*x ); 
     1505          break; 
     1506        } 
     1507        case PolynomialDistortion: 
     1508        { 
     1509          register long 
     1510            k; 
     1511 
     1512          double 
     1513            dudx,dudy,dvdx,dvdy; 
     1514 
     1515          point.x=point.y=dudx=dudy=dvdx=dvdy=0.0; 
     1516          for(k=0; k < coefficients[0]; k++) { 
     1517            point.x += poly_term(k,x,y)*coefficients[k+1]; 
     1518            dudx += poly_term_dx(k,x,y)*coefficients[k+1]; 
     1519            dudy += poly_term_dy(k,x,y)*coefficients[k+1]; 
     1520            point.y += poly_term(k,x,y)*coefficients[k+1+(long)coefficients[0]]; 
     1521            dvdx += poly_term_dx(k,x,y)*coefficients[k+1+(long)coefficients[0]]; 
     1522            dvdy += poly_term_dy(k,x,y)*coefficients[k+1+(long)coefficients[0]]; 
     1523          } 
     1524          ScaleResampleFilter( resample_filter[id], dudx,dudy,dvdx,dvdy ); 
     1525          break; 
     1526        } 
    11471527        case ArcDistortion: 
    11481528        { 
    1149           double radius = sqrt((double) x*x+y*y); 
     1529          /* what is the angle and radius in the destination image */ 
    11501530          point.x = (atan2((double)y,(double)x) - coefficients[0])/(2*MagickPI); 
    1151           point.x -= MagickRound(point.x); 
    1152           point.x = point.x*coefficients[1] + coefficients[4];