granicus.if.org Git - imagemagick/blob - MagickCore/resample.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %           RRRR    EEEEE   SSSSS   AAA   M   M  PPPP   L      EEEEE          %
   7 %           R   R   E       SS     A   A  MM MM  P   P  L      E              %
   8 %           RRRR    EEE      SSS   AAAAA  M M M  PPPP   L      EEE            %
   9 %           R R     E          SS  A   A  M   M  P      L      E              %
  10 %           R  R    EEEEE   SSSSS  A   A  M   M  P      LLLLL  EEEEE          %
  11 %                                                                             %
  12 %                                                                             %
  13 %                      MagickCore Pixel Resampling Methods                    %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                John Cristy                                  %
  17 %                              Anthony Thyssen                                %
  18 %                                August 2007                                  %
  19 %                                                                             %
  20 %                                                                             %
  21 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  22 %  dedicated to making software imaging solutions freely available.           %
  23 %                                                                             %
  24 %  You may not use this file except in compliance with the License.  You may  %
  25 %  obtain a copy of the License at                                            %
  26 %                                                                             %
  27 %    http://www.imagemagick.org/script/license.php                            %
  28 %                                                                             %
  29 %  Unless required by applicable law or agreed to in writing, software        %
  30 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  31 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  32 %  See the License for the specific language governing permissions and        %
  33 %  limitations under the License.                                             %
  34 %                                                                             %
  35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/artifact.h"
  45 #include "MagickCore/color-private.h"
  46 #include "MagickCore/cache.h"
  47 #include "MagickCore/draw.h"
  48 #include "MagickCore/exception-private.h"
  49 #include "MagickCore/gem.h"
  50 #include "MagickCore/image.h"
  51 #include "MagickCore/image-private.h"
  52 #include "MagickCore/log.h"
  53 #include "MagickCore/magick.h"
  54 #include "MagickCore/memory_.h"
  55 #include "MagickCore/pixel.h"
  56 #include "MagickCore/pixel-accessor.h"
  57 #include "MagickCore/quantum.h"
  58 #include "MagickCore/random_.h"
  59 #include "MagickCore/resample.h"
  60 #include "MagickCore/resize.h"
  61 #include "MagickCore/resize-private.h"
  62 #include "MagickCore/transform.h"
  63 #include "MagickCore/signature-private.h"
  64 #include "MagickCore/utility.h"
  65 #include "MagickCore/utility-private.h"
  66 /*
  67   EWA Resampling Options
  68 */
  69
  70 /* select ONE resampling method */
  71 #define EWA 1                 /* Normal EWA handling - raw or clamped */
  72                               /* if 0 then use "High Quality EWA" */
  73 #define EWA_CLAMP 1           /* EWA Clamping from Nicolas Robidoux */
  74
  75 #define FILTER_LUT 1          /* Use a LUT rather then direct filter calls */
  76
  77 /* output debugging information */
  78 #define DEBUG_ELLIPSE 0       /* output ellipse info for debug */
  79 #define DEBUG_HIT_MISS 0      /* output hit/miss pixels (as gnuplot commands) */
  80 #define DEBUG_NO_PIXEL_HIT 0  /* Make pixels that fail to hit anything - RED */
  81
  82 #if ! FILTER_DIRECT
  83 #define WLUT_WIDTH 1024       /* size of the filter cache */
  84 #endif
  85
  86 /*
  87   Typedef declarations.
  88 */
  89 struct _ResampleFilter
  90 {
  91   CacheView
  92     *view;
  93
  94   Image
  95     *image;
  96
  97   ExceptionInfo
  98     *exception;
  99
 100   MagickBooleanType
 101     debug;
 102
 103   /* Information about image being resampled */
 104   ssize_t
 105     image_area;
 106
 107   PixelInterpolateMethod
 108     interpolate;
 109
 110   VirtualPixelMethod
 111     virtual_pixel;
 112
 113   FilterTypes
 114     filter;
 115
 116   /* processing settings needed */
 117   MagickBooleanType
 118     limit_reached,
 119     do_interpolate,
 120     average_defined;
 121
 122   PixelInfo
 123     average_pixel;
 124
 125   /* current ellipitical area being resampled around center point */
 126   double
 127     A, B, C,
 128     Vlimit, Ulimit, Uwidth, slope;
 129
 130 #if FILTER_LUT
 131   /* LUT of weights for filtered average in elliptical area */
 132   double
 133     filter_lut[WLUT_WIDTH];
 134 #else
 135   /* Use a Direct call to the filter functions */
 136   ResizeFilter
 137     *filter_def;
 138
 139   double
 140     F;
 141 #endif
 142
 143   /* the practical working support of the filter */
 144   double
 145     support;
 146
 147   size_t
 148     signature;
 149 };
 150 \f
 151 /*
 152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 153 %                                                                             %
 154 %                                                                             %
 155 %                                                                             %
 156 %   A c q u i r e R e s a m p l e I n f o                                     %
 157 %                                                                             %
 158 %                                                                             %
 159 %                                                                             %
 160 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 161 %
 162 %  AcquireResampleFilter() initializes the information resample needs do to a
 163 %  scaled lookup of a color from an image, using area sampling.
 164 %
 165 %  The algorithm is based on a Elliptical Weighted Average, where the pixels
 166 %  found in a large elliptical area is averaged together according to a
 167 %  weighting (filter) function.  For more details see "Fundamentals of Texture
 168 %  Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
 169 %  1989.  Available for free from, http://www.cs.cmu.edu/~ph/
 170 %
 171 %  As EWA resampling (or any sort of resampling) can require a lot of
 172 %  calculations to produce a distorted scaling of the source image for each
 173 %  output pixel, the ResampleFilter structure generated holds that information
 174 %  between individual image resampling.
 175 %
 176 %  This function will make the appropriate AcquireCacheView() calls
 177 %  to view the image, calling functions do not need to open a cache view.
 178 %
 179 %  Usage Example...
 180 %      resample_filter=AcquireResampleFilter(image,exception);
 181 %      SetResampleFilter(resample_filter, GaussianFilter, 1.0);
 182 %      for (y=0; y < (ssize_t) image->rows; y++) {
 183 %        for (x=0; x < (ssize_t) image->columns; x++) {
 184 %          u= ....;   v= ....;
 185 %          ScaleResampleFilter(resample_filter, ... scaling vectors ...);
 186 %          (void) ResamplePixelColor(resample_filter,u,v,&pixel);
 187 %          ... assign resampled pixel value ...
 188 %        }
 189 %      }
 190 %      DestroyResampleFilter(resample_filter);
 191 %
 192 %  The format of the AcquireResampleFilter method is:
 193 %
 194 %     ResampleFilter *AcquireResampleFilter(const Image *image,
 195 %       ExceptionInfo *exception)
 196 %
 197 %  A description of each parameter follows:
 198 %
 199 %    o image: the image.
 200 %
 201 %    o exception: return any errors or warnings in this structure.
 202 %
 203 */
 204 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
 205   ExceptionInfo *exception)
 206 {
 207   register ResampleFilter
 208     *resample_filter;
 209
 210   assert(image != (Image *) NULL);
 211   assert(image->signature == MagickSignature);
 212   if (image->debug != MagickFalse)
 213     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 214   assert(exception != (ExceptionInfo *) NULL);
 215   assert(exception->signature == MagickSignature);
 216
 217   resample_filter=(ResampleFilter *) AcquireMagickMemory(
 218     sizeof(*resample_filter));
 219   if (resample_filter == (ResampleFilter *) NULL)
 220     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 221   (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
 222
 223   resample_filter->exception=exception;
 224   resample_filter->image=ReferenceImage((Image *) image);
 225   resample_filter->view=AcquireCacheView(resample_filter->image);
 226
 227   resample_filter->debug=IsEventLogging();
 228   resample_filter->signature=MagickSignature;
 229
 230   resample_filter->image_area=(ssize_t) (image->columns*image->rows);
 231   resample_filter->average_defined = MagickFalse;
 232
 233   /* initialise the resampling filter settings */
 234   SetResampleFilter(resample_filter, image->filter, image->blur);
 235   (void) SetResampleFilterInterpolateMethod(resample_filter,
 236     image->interpolate);
 237   (void) SetResampleFilterVirtualPixelMethod(resample_filter,
 238     GetImageVirtualPixelMethod(image));
 239
 240   return(resample_filter);
 241 }
 242 \f
 243 /*
 244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 245 %                                                                             %
 246 %                                                                             %
 247 %                                                                             %
 248 %   D e s t r o y R e s a m p l e I n f o                                     %
 249 %                                                                             %
 250 %                                                                             %
 251 %                                                                             %
 252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 253 %
 254 %  DestroyResampleFilter() finalizes and cleans up the resampling
 255 %  resample_filter as returned by AcquireResampleFilter(), freeing any memory
 256 %  or other information as needed.
 257 %
 258 %  The format of the DestroyResampleFilter method is:
 259 %
 260 %      ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
 261 %
 262 %  A description of each parameter follows:
 263 %
 264 %    o resample_filter: resampling information structure
 265 %
 266 */
 267 MagickExport ResampleFilter *DestroyResampleFilter(
 268   ResampleFilter *resample_filter)
 269 {
 270   assert(resample_filter != (ResampleFilter *) NULL);
 271   assert(resample_filter->signature == MagickSignature);
 272   assert(resample_filter->image != (Image *) NULL);
 273   if (resample_filter->debug != MagickFalse)
 274     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
 275       resample_filter->image->filename);
 276   resample_filter->view=DestroyCacheView(resample_filter->view);
 277   resample_filter->image=DestroyImage(resample_filter->image);
 278 #if ! FILTER_LUT
 279   resample_filter->filter_def=DestroyResizeFilter(resample_filter->filter_def);
 280 #endif
 281   resample_filter->signature=(~MagickSignature);
 282   resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
 283   return(resample_filter);
 284 }
 285 \f
 286 /*
 287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 288 %                                                                             %
 289 %                                                                             %
 290 %                                                                             %
 291 %   R e s a m p l e P i x e l C o l o r                                       %
 292 %                                                                             %
 293 %                                                                             %
 294 %                                                                             %
 295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 296 %
 297 %  ResamplePixelColor() samples the pixel values surrounding the location
 298 %  given using an elliptical weighted average, at the scale previously
 299 %  calculated, and in the most efficent manner possible for the
 300 %  VirtualPixelMethod setting.
 301 %
 302 %  The format of the ResamplePixelColor method is:
 303 %
 304 %     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
 305 %       const double u0,const double v0,PixelInfo *pixel)
 306 %
 307 %  A description of each parameter follows:
 308 %
 309 %    o resample_filter: the resample filter.
 310 %
 311 %    o u0,v0: A double representing the center of the area to resample,
 312 %        The distortion transformed transformed x,y coordinate.
 313 %
 314 %    o pixel: the resampled pixel is returned here.
 315 %
 316 */
 317 MagickExport MagickBooleanType ResamplePixelColor(
 318   ResampleFilter *resample_filter,const double u0,const double v0,
 319   PixelInfo *pixel)
 320 {
 321   MagickBooleanType
 322     status;
 323
 324   ssize_t u,v, v1, v2, uw, hit;
 325   double u1;
 326   double U,V,Q,DQ,DDQ;
 327   double divisor_c,divisor_m;
 328   register double weight;
 329   register const Quantum *pixels;
 330   assert(resample_filter != (ResampleFilter *) NULL);
 331   assert(resample_filter->signature == MagickSignature);
 332
 333   status=MagickTrue;
 334   /* GetPixelInfo(resample_filter->image,pixel); */
 335   if ( resample_filter->do_interpolate ) {
 336     status=InterpolatePixelInfo(resample_filter->image,
 337       resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
 338       resample_filter->exception);
 339     return(status);
 340   }
 341
 342 #if DEBUG_ELLIPSE
 343   (void) FormatLocaleFile(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
 344 #endif
 345
 346   /*
 347     Does resample area Miss the image?
 348     And is that area a simple solid color - then return that color
 349   */
 350   hit = 0;
 351   switch ( resample_filter->virtual_pixel ) {
 352     case BackgroundVirtualPixelMethod:
 353     case TransparentVirtualPixelMethod:
 354     case BlackVirtualPixelMethod:
 355     case GrayVirtualPixelMethod:
 356     case WhiteVirtualPixelMethod:
 357     case MaskVirtualPixelMethod:
 358       if ( resample_filter->limit_reached
 359            || u0 + resample_filter->Ulimit < 0.0
 360            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 361            || v0 + resample_filter->Vlimit < 0.0
 362            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 363            )
 364         hit++;
 365       break;
 366
 367     case UndefinedVirtualPixelMethod:
 368     case EdgeVirtualPixelMethod:
 369       if (    ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
 370            || ( u0 + resample_filter->Ulimit < 0.0
 371                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 372            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 373                 && v0 + resample_filter->Vlimit < 0.0 )
 374            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 375                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 376            )
 377         hit++;
 378       break;
 379     case HorizontalTileVirtualPixelMethod:
 380       if (    v0 + resample_filter->Vlimit < 0.0
 381            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 382            )
 383         hit++;  /* outside the horizontally tiled images. */
 384       break;
 385     case VerticalTileVirtualPixelMethod:
 386       if (    u0 + resample_filter->Ulimit < 0.0
 387            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 388            )
 389         hit++;  /* outside the vertically tiled images. */
 390       break;
 391     case DitherVirtualPixelMethod:
 392       if (    ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
 393            || ( u0 + resample_filter->Ulimit < -32.0
 394                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 395            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 396                 && v0 + resample_filter->Vlimit < -32.0 )
 397            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 398                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 399            )
 400         hit++;
 401       break;
 402     case TileVirtualPixelMethod:
 403     case MirrorVirtualPixelMethod:
 404     case RandomVirtualPixelMethod:
 405     case HorizontalTileEdgeVirtualPixelMethod:
 406     case VerticalTileEdgeVirtualPixelMethod:
 407     case CheckerTileVirtualPixelMethod:
 408       /* resampling of area is always needed - no VP limits */
 409       break;
 410   }
 411   if ( hit ) {
 412     /* whole area is a solid color -- just return that color */
 413     status=InterpolatePixelInfo(resample_filter->image,
 414       resample_filter->view,IntegerInterpolatePixel,u0,v0,pixel,
 415       resample_filter->exception);
 416     return(status);
 417   }
 418
 419   /*
 420     Scaling limits reached, return an 'averaged' result.
 421   */
 422   if ( resample_filter->limit_reached ) {
 423     switch ( resample_filter->virtual_pixel ) {
 424       /*  This is always handled by the above, so no need.
 425         case BackgroundVirtualPixelMethod:
 426         case ConstantVirtualPixelMethod:
 427         case TransparentVirtualPixelMethod:
 428         case GrayVirtualPixelMethod,
 429         case WhiteVirtualPixelMethod
 430         case MaskVirtualPixelMethod:
 431       */
 432       case UndefinedVirtualPixelMethod:
 433       case EdgeVirtualPixelMethod:
 434       case DitherVirtualPixelMethod:
 435       case HorizontalTileEdgeVirtualPixelMethod:
 436       case VerticalTileEdgeVirtualPixelMethod:
 437         /* We need an average edge pixel, from the correct edge!
 438            How should I calculate an average edge color?
 439            Just returning an averaged neighbourhood,
 440            works well in general, but falls down for TileEdge methods.
 441            This needs to be done properly!!!!!!
 442         */
 443         status=InterpolatePixelInfo(resample_filter->image,
 444           resample_filter->view,AverageInterpolatePixel,u0,v0,pixel,
 445           resample_filter->exception);
 446         break;
 447       case HorizontalTileVirtualPixelMethod:
 448       case VerticalTileVirtualPixelMethod:
 449         /* just return the background pixel - Is there more direct way? */
 450         status=InterpolatePixelInfo(resample_filter->image,
 451           resample_filter->view,IntegerInterpolatePixel,-1.0,-1.0,pixel,
 452           resample_filter->exception);
 453         break;
 454       case TileVirtualPixelMethod:
 455       case MirrorVirtualPixelMethod:
 456       case RandomVirtualPixelMethod:
 457       case CheckerTileVirtualPixelMethod:
 458       default:
 459         /* generate a average color of the WHOLE image */
 460         if ( resample_filter->average_defined == MagickFalse ) {
 461           Image
 462             *average_image;
 463
 464           CacheView
 465             *average_view;
 466
 467           GetPixelInfo(resample_filter->image,(PixelInfo *)
 468             &resample_filter->average_pixel);
 469           resample_filter->average_defined=MagickTrue;
 470
 471           /* Try to get an averaged pixel color of whole image */
 472           average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
 473             resample_filter->exception);
 474           if (average_image == (Image *) NULL)
 475             {
 476               *pixel=resample_filter->average_pixel; /* FAILED */
 477               break;
 478             }
 479           average_view=AcquireCacheView(average_image);
 480           pixels=GetCacheViewVirtualPixels(average_view,0,0,1,1,
 481             resample_filter->exception);
 482           if (pixels == (const Quantum *) NULL) {
 483             average_view=DestroyCacheView(average_view);
 484             average_image=DestroyImage(average_image);
 485             *pixel=resample_filter->average_pixel; /* FAILED */
 486             break;
 487           }
 488           SetPixelInfo(resample_filter->image,pixels,
 489             &(resample_filter->average_pixel));
 490           average_view=DestroyCacheView(average_view);
 491           average_image=DestroyImage(average_image);
 492
 493           if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod )
 494             {
 495               /* CheckerTile is avergae of image average half background */
 496               /* FUTURE: replace with a 50% blend of both pixels */
 497
 498               weight = QuantumScale*((MagickRealType)
 499                 resample_filter->average_pixel.alpha);
 500               resample_filter->average_pixel.red *= weight;
 501               resample_filter->average_pixel.green *= weight;
 502               resample_filter->average_pixel.blue *= weight;
 503               divisor_c = weight;
 504
 505               weight = QuantumScale*((MagickRealType)
 506                 resample_filter->image->background_color.alpha);
 507               resample_filter->average_pixel.red +=
 508                       weight*resample_filter->image->background_color.red;
 509               resample_filter->average_pixel.green +=
 510                       weight*resample_filter->image->background_color.green;
 511               resample_filter->average_pixel.blue +=
 512                       weight*resample_filter->image->background_color.blue;
 513               resample_filter->average_pixel.alpha +=
 514                       resample_filter->image->background_color.alpha;
 515               divisor_c += weight;
 516
 517               resample_filter->average_pixel.red /= divisor_c;
 518               resample_filter->average_pixel.green /= divisor_c;
 519               resample_filter->average_pixel.blue /= divisor_c;
 520               resample_filter->average_pixel.alpha /= 2;
 521
 522             }
 523         }
 524         *pixel=resample_filter->average_pixel;
 525         break;
 526     }
 527     return(status);
 528   }
 529
 530   /*
 531     Initialize weighted average data collection
 532   */
 533   hit = 0;
 534   divisor_c = 0.0;
 535   divisor_m = 0.0;
 536   pixel->red = pixel->green = pixel->blue = 0.0;
 537   if (pixel->colorspace == CMYKColorspace)
 538     pixel->black = 0.0;
 539   if (pixel->matte != MagickFalse)
 540     pixel->alpha = 0.0;
 541
 542   /*
 543     Determine the parellelogram bounding box fitted to the ellipse
 544     centered at u0,v0.  This area is bounding by the lines...
 545   */
 546   v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit);  /* range of scan lines */
 547   v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
 548
 549   /* scan line start and width accross the parallelogram */
 550   u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
 551   uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
 552
 553 #if DEBUG_ELLIPSE
 554   (void) FormatLocaleFile(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
 555   (void) FormatLocaleFile(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
 556 #else
 557 # define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
 558 #endif
 559
 560   /*
 561     Do weighted resampling of all pixels,  within the scaled ellipse,
 562     bound by a Parellelogram fitted to the ellipse.
 563   */
 564   DDQ = 2*resample_filter->A;
 565   for( v=v1; v<=v2;  v++ ) {
 566 #if DEBUG_HIT_MISS
 567     long uu = ceil(u1);   /* actual pixel location (for debug only) */
 568     (void) FormatLocaleFile(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
 569 #endif
 570     u = (ssize_t)ceil(u1);        /* first pixel in scanline */
 571     u1 += resample_filter->slope; /* start of next scan line */
 572
 573
 574     /* location of this first pixel, relative to u0,v0 */
 575     U = (double)u-u0;
 576     V = (double)v-v0;
 577
 578     /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
 579     Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
 580     DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
 581
 582     /* get the scanline of pixels for this v */
 583     pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
 584       1,resample_filter->exception);
 585     if (pixels == (const Quantum *) NULL)
 586       return(MagickFalse);
 587
 588     /* count up the weighted pixel colors */
 589     for( u=0; u<uw; u++ ) {
 590 #if FILTER_LUT
 591       /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
 592       if ( Q < (double)WLUT_WIDTH ) {
 593         weight = resample_filter->filter_lut[(int)Q];
 594 #else
 595       /* Note that the ellipse has been pre-scaled so F = support^2 */
 596       if ( Q < (double)resample_filter->F ) {
 597         weight = GetResizeFilterWeight(resample_filter->filter_def,
 598              sqrt(Q));    /* a SquareRoot!  Arrggghhhhh... */
 599 #endif
 600
 601         pixel->alpha  += weight*GetPixelAlpha(resample_filter->image,pixels);
 602         divisor_m += weight;
 603
 604         if (pixel->matte != MagickFalse)
 605           weight *= QuantumScale*((MagickRealType) GetPixelAlpha(resample_filter->image,pixels));
 606         pixel->red   += weight*GetPixelRed(resample_filter->image,pixels);
 607         pixel->green += weight*GetPixelGreen(resample_filter->image,pixels);
 608         pixel->blue  += weight*GetPixelBlue(resample_filter->image,pixels);
 609         if (pixel->colorspace == CMYKColorspace)
 610           pixel->black += weight*GetPixelBlack(resample_filter->image,pixels);
 611         divisor_c += weight;
 612
 613         hit++;
 614 #if DEBUG_HIT_MISS
 615         /* mark the pixel according to hit/miss of the ellipse */
 616         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
 617                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
 618         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
 619                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
 620       } else {
 621         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
 622                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
 623         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
 624                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
 625       }
 626       uu++;
 627 #else
 628       }
 629 #endif
 630       pixels+=GetPixelChannels(resample_filter->image);
 631       Q += DQ;
 632       DQ += DDQ;
 633     }
 634   }
 635 #if DEBUG_ELLIPSE
 636   (void) FormatLocaleFile(stderr, "Hit=%ld;  Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
 637 #endif
 638
 639   /*
 640     Result sanity check -- this should NOT happen
 641   */
 642   if ( hit == 0 ) {
 643     /* not enough pixels in resampling, resort to direct interpolation */
 644 #if DEBUG_NO_PIXEL_HIT
 645     pixel->alpha = pixel->red = pixel->green = pixel->blue = 0;
 646     pixel->red = QuantumRange; /* show pixels for which EWA fails */
 647 #else
 648     status=InterpolatePixelInfo(resample_filter->image,
 649       resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
 650       resample_filter->exception);
 651 #endif
 652     return status;
 653   }
 654
 655   /*
 656     Finialize results of resampling
 657   */
 658   divisor_m = 1.0/divisor_m;
 659   pixel->alpha = (MagickRealType) ClampToQuantum(divisor_m*pixel->alpha);
 660   divisor_c = 1.0/divisor_c;
 661   pixel->red   = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
 662   pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
 663   pixel->blue  = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
 664   if (pixel->colorspace == CMYKColorspace)
 665     pixel->black = (MagickRealType) ClampToQuantum(divisor_c*pixel->black);
 666   return(MagickTrue);
 667 }
 668 \f
 669 #if EWA && EWA_CLAMP
 670 /*
 671 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 672 %                                                                             %
 673 %                                                                             %
 674 %                                                                             %
 675 -   C l a m p U p A x e s                                                     %
 676 %                                                                             %
 677 %                                                                             %
 678 %                                                                             %
 679 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 680 %
 681 % ClampUpAxes() function converts the input vectors into a major and
 682 % minor axis unit vectors, and their magnitude.  This allows us to
 683 % ensure that the ellipse generated is never smaller than the unit
 684 % circle and thus never too small for use in EWA resampling.
 685 %
 686 % This purely mathematical 'magic' was provided by Professor Nicolas
 687 % Robidoux and his Masters student Chantal Racette.
 688 %
 689 % Reference: "We Recommend Singular Value Decomposition", David Austin
 690 %   http://www.ams.org/samplings/feature-column/fcarc-svd
 691 %
 692 % By generating major and minor axis vectors, we can actually use the
 693 % ellipse in its "canonical form", by remapping the dx,dy of the
 694 % sampled point into distances along the major and minor axis unit
 695 % vectors.
 696 %
 697 % Reference: http://en.wikipedia.org/wiki/Ellipse#Canonical_form
 698 */
 699 static inline void ClampUpAxes(const double dux,
 700                                const double dvx,
 701                                const double duy,
 702                                const double dvy,
 703                                double *major_mag,
 704                                double *minor_mag,
 705                                double *major_unit_x,
 706                                double *major_unit_y,
 707                                double *minor_unit_x,
 708                                double *minor_unit_y)
 709 {
 710   /*
 711    * ClampUpAxes takes an input 2x2 matrix
 712    *
 713    * [ a b ] = [ dux duy ]
 714    * [ c d ] = [ dvx dvy ]
 715    *
 716    * and computes from it the major and minor axis vectors [major_x,
 717    * major_y] and [minor_x,minor_y] of the smallest ellipse containing
 718    * both the unit disk and the ellipse which is the image of the unit
 719    * disk by the linear transformation
 720    *
 721    * [ dux duy ] [S] = [s]
 722    * [ dvx dvy ] [T] = [t]
 723    *
 724    * (The vector [S,T] is the difference between a position in output
 725    * space and [X,Y]; the vector [s,t] is the difference between a
 726    * position in input space and [x,y].)
 727    */
 728   /*
 729    * Output:
 730    *
 731    * major_mag is the half-length of the major axis of the "new"
 732    * ellipse.
 733    *
 734    * minor_mag is the half-length of the minor axis of the "new"
 735    * ellipse.
 736    *
 737    * major_unit_x is the x-coordinate of the major axis direction vector
 738    * of both the "old" and "new" ellipses.
 739    *
 740    * major_unit_y is the y-coordinate of the major axis direction vector.
 741    *
 742    * minor_unit_x is the x-coordinate of the minor axis direction vector.
 743    *
 744    * minor_unit_y is the y-coordinate of the minor axis direction vector.
 745    *
 746    * Unit vectors are useful for computing projections, in particular,
 747    * to compute the distance between a point in output space and the
 748    * center of a unit disk in output space, using the position of the
 749    * corresponding point [s,t] in input space. Following the clamping,
 750    * the square of this distance is
 751    *
 752    * ( ( s * major_unit_x + t * major_unit_y ) / major_mag )^2
 753    * +
 754    * ( ( s * minor_unit_x + t * minor_unit_y ) / minor_mag )^2
 755    *
 756    * If such distances will be computed for many [s,t]'s, it makes
 757    * sense to actually compute the reciprocal of major_mag and
 758    * minor_mag and multiply them by the above unit lengths.
 759    *
 760    * Now, if you want to modify the input pair of tangent vectors so
 761    * that it defines the modified ellipse, all you have to do is set
 762    *
 763    * newdux = major_mag * major_unit_x
 764    * newdvx = major_mag * major_unit_y
 765    * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
 766    * newdvy = minor_mag * minor_unit_y = minor_mag *  major_unit_x
 767    *
 768    * and use these tangent vectors as if they were the original ones.
 769    * Usually, this is a drastic change in the tangent vectors even if
 770    * the singular values are not clamped; for example, the minor axis
 771    * vector always points in a direction which is 90 degrees
 772    * counterclockwise from the direction of the major axis vector.
 773    */
 774   /*
 775    * Discussion:
 776    *
 777    * GOAL: Fix things so that the pullback, in input space, of a disk
 778    * of radius r in output space is an ellipse which contains, at
 779    * least, a disc of radius r. (Make this hold for any r>0.)
 780    *
 781    * ESSENCE OF THE METHOD: Compute the product of the first two
 782    * factors of an SVD of the linear transformation defining the
 783    * ellipse and make sure that both its columns have norm at least 1.
 784    * Because rotations and reflexions map disks to themselves, it is
 785    * not necessary to compute the third (rightmost) factor of the SVD.
 786    *
 787    * DETAILS: Find the singular values and (unit) left singular
 788    * vectors of Jinv, clampling up the singular values to 1, and
 789    * multiply the unit left singular vectors by the new singular
 790    * values in order to get the minor and major ellipse axis vectors.
 791    *
 792    * Image resampling context:
 793    *
 794    * The Jacobian matrix of the transformation at the output point
 795    * under consideration is defined as follows:
 796    *
 797    * Consider the transformation (x,y) -> (X,Y) from input locations
 798    * to output locations. (Anthony Thyssen, elsewhere in resample.c,
 799    * uses the notation (u,v) -> (x,y).)
 800    *
 801    * The Jacobian matrix of the transformation at (x,y) is equal to
 802    *
 803    *   J = [ A, B ] = [ dX/dx, dX/dy ]
 804    *       [ C, D ]   [ dY/dx, dY/dy ]
 805    *
 806    * that is, the vector [A,C] is the tangent vector corresponding to
 807    * input changes in the horizontal direction, and the vector [B,D]
 808    * is the tangent vector corresponding to input changes in the
 809    * vertical direction.
 810    *
 811    * In the context of resampling, it is natural to use the inverse
 812    * Jacobian matrix Jinv because resampling is generally performed by
 813    * pulling pixel locations in the output image back to locations in
 814    * the input image. Jinv is
 815    *
 816    *   Jinv = [ a, b ] = [ dx/dX, dx/dY ]
 817    *          [ c, d ]   [ dy/dX, dy/dY ]
 818    *
 819    * Note: Jinv can be computed from J with the following matrix
 820    * formula:
 821    *
 822    *   Jinv = 1/(A*D-B*C) [  D, -B ]
 823    *                      [ -C,  A ]
 824    *
 825    * What we do is modify Jinv so that it generates an ellipse which
 826    * is as close as possible to the original but which contains the
 827    * unit disk. This can be accomplished as follows:
 828    *
 829    * Let
 830    *
 831    *   Jinv = U Sigma V^T
 832    *
 833    * be an SVD decomposition of Jinv. (The SVD is not unique, but the
 834    * final ellipse does not depend on the particular SVD.)
 835    *
 836    * We could clamp up the entries of the diagonal matrix Sigma so
 837    * that they are at least 1, and then set
 838    *
 839    *   Jinv = U newSigma V^T.
 840    *
 841    * However, we do not need to compute V for the following reason:
 842    * V^T is an orthogonal matrix (that is, it represents a combination
 843    * of rotations and reflexions) so that it maps the unit circle to
 844    * itself. For this reason, the exact value of V does not affect the
 845    * final ellipse, and we can choose V to be the identity
 846    * matrix. This gives
 847    *
 848    *   Jinv = U newSigma.
 849    *
 850    * In the end, we return the two diagonal entries of newSigma
 851    * together with the two columns of U.
 852    */
 853   /*
 854    * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
 855    * of Laurentian University with insightful suggestions from Anthony
 856    * Thyssen and funding from the National Science and Engineering
 857    * Research Council of Canada. It is distinguished from its
 858    * predecessors by its efficient handling of degenerate cases.
 859    *
 860    * The idea of clamping up the EWA ellipse's major and minor axes so
 861    * that the result contains the reconstruction kernel filter support
 862    * is taken from Andreas Gustaffson's Masters thesis "Interactive
 863    * Image Warping", Helsinki University of Technology, Faculty of
 864    * Information Technology, 59 pages, 1993 (see Section 3.6).
 865    *
 866    * The use of the SVD to clamp up the singular values of the
 867    * Jacobian matrix of the pullback transformation for EWA resampling
 868    * is taken from the astrophysicist Craig DeForest.  It is
 869    * implemented in his PDL::Transform code (PDL = Perl Data
 870    * Language).
 871    */
 872   const double a = dux;
 873   const double b = duy;
 874   const double c = dvx;
 875   const double d = dvy;
 876   /*
 877    * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
 878    * squares of the singular values of Jinv.
 879    */
 880   const double aa = a*a;
 881   const double bb = b*b;
 882   const double cc = c*c;
 883   const double dd = d*d;
 884   /*
 885    * Eigenvectors of n are left singular vectors of Jinv.
 886    */
 887   const double n11 = aa+bb;
 888   const double n12 = a*c+b*d;
 889   const double n21 = n12;
 890   const double n22 = cc+dd;
 891   const double det = a*d-b*c;
 892   const double twice_det = det+det;
 893   const double frobenius_squared = n11+n22;
 894   const double discriminant =
 895     (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
 896   const double sqrt_discriminant = sqrt(discriminant);
 897   /*
 898    * s1 is the largest singular value of the inverse Jacobian
 899    * matrix. In other words, its reciprocal is the smallest singular
 900    * value of the Jacobian matrix itself.
 901    * If s1 = 0, both singular values are 0, and any orthogonal pair of
 902    * left and right factors produces a singular decomposition of Jinv.
 903    */
 904   /*
 905    * Initially, we only compute the squares of the singular values.
 906    */
 907   const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
 908   /*
 909    * s2 the smallest singular value of the inverse Jacobian
 910    * matrix. Its reciprocal is the largest singular value of the
 911    * Jacobian matrix itself.
 912    */
 913   const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
 914   const double s1s1minusn11 = s1s1-n11;
 915   const double s1s1minusn22 = s1s1-n22;
 916   /*
 917    * u1, the first column of the U factor of a singular decomposition
 918    * of Jinv, is a (non-normalized) left singular vector corresponding
 919    * to s1. It has entries u11 and u21. We compute u1 from the fact
 920    * that it is an eigenvector of n corresponding to the eigenvalue
 921    * s1^2.
 922    */
 923   const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
 924   const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
 925   /*
 926    * The following selects the largest row of n-s1^2 I as the one
 927    * which is used to find the eigenvector. If both s1^2-n11 and
 928    * s1^2-n22 are zero, n-s1^2 I is the zero matrix.  In that case,
 929    * any vector is an eigenvector; in addition, norm below is equal to
 930    * zero, and, in exact arithmetic, this is the only case in which
 931    * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
 932    * if norm = 0 safely takes care of all cases.
 933    */
 934   const double temp_u11 =
 935     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
 936   const double temp_u21 =
 937     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
 938   const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
 939   /*
 940    * Finalize the entries of first left singular vector (associated
 941    * with the largest singular value).
 942    */
 943   const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
 944   const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
 945   /*
 946    * Clamp the singular values up to 1.
 947    */
 948   *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
 949   *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
 950   /*
 951    * Return the unit major and minor axis direction vectors.
 952    */
 953   *major_unit_x = u11;
 954   *major_unit_y = u21;
 955   *minor_unit_x = -u21;
 956   *minor_unit_y = u11;
 957 }
 958 \f
 959 #endif
 960 /*
 961 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 962 %                                                                             %
 963 %                                                                             %
 964 %                                                                             %
 965 %   S c a l e R e s a m p l e F i l t e r                                     %
 966 %                                                                             %
 967 %                                                                             %
 968 %                                                                             %
 969 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 970 %
 971 %  ScaleResampleFilter() does all the calculations needed to resample an image
 972 %  at a specific scale, defined by two scaling vectors.  This not using
 973 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
 974 %  generation of a angled ellipse.
 975 %
 976 %  As only two deritive scaling vectors are used the center of the ellipse
 977 %  must be the center of the lookup.  That is any curvature that the
 978 %  distortion may produce is discounted.
 979 %
 980 %  The input vectors are produced by either finding the derivitives of the
 981 %  distortion function, or the partial derivitives from a distortion mapping.
 982 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
 983 %  calculated from other derivatives.  For example you could use  dr,da/r
 984 %  polar coordinate vector scaling vectors
 985 %
 986 %  If   u,v =  DistortEquation(x,y)   OR   u = Fu(x,y); v = Fv(x,y)
 987 %  Then the scaling vectors are determined from the deritives...
 988 %      du/dx, dv/dx     and    du/dy, dv/dy
 989 %  If the resulting scaling vectors is othogonally aligned then...
 990 %      dv/dx = 0   and   du/dy  =  0
 991 %  Producing an othogonally alligned ellipse in source space for the area to
 992 %  be resampled.
 993 %
 994 %  Note that scaling vectors are different to argument order.  Argument order
 995 %  is the general order the deritives are extracted from the distortion
 996 %  equations, and not the scaling vectors. As such the middle two vaules
 997 %  may be swapped from what you expect.  Caution is advised.
 998 %
 999 %  WARNING: It is assumed that any SetResampleFilter() method call will
1000 %  always be performed before the ScaleResampleFilter() method, so that the
1001 %  size of the ellipse will match the support for the resampling filter being
1002 %  used.
1003 %
1004 %  The format of the ScaleResampleFilter method is:
1005 %
1006 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
1007 %       const double dux,const double duy,const double dvx,const double dvy)
1008 %
1009 %  A description of each parameter follows:
1010 %
1011 %    o resample_filter: the resampling resample_filterrmation defining the
1012 %      image being resampled
1013 %
1014 %    o dux,duy,dvx,dvy:
1015 %         The deritives or scaling vectors defining the EWA ellipse.
1016 %         NOTE: watch the order, which is based on the order deritives
1017 %         are usally determined from distortion equations (see above).
1018 %         The middle two values may need to be swapped if you are thinking
1019 %         in terms of scaling vectors.
1020 %
1021 */
1022 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1023   const double dux,const double duy,const double dvx,const double dvy)
1024 {
1025   double A,B,C,F;
1026
1027   assert(resample_filter != (ResampleFilter *) NULL);
1028   assert(resample_filter->signature == MagickSignature);
1029
1030   resample_filter->limit_reached = MagickFalse;
1031
1032   /* A 'point' filter forces use of interpolation instead of area sampling */
1033   if ( resample_filter->filter == PointFilter )
1034     return; /* EWA turned off - nothing to do */
1035
1036 #if DEBUG_ELLIPSE
1037   (void) FormatLocaleFile(stderr, "# -----\n" );
1038   (void) FormatLocaleFile(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy=%lf;\n",
1039        dux, dvx, duy, dvy);
1040 #endif
1041
1042   /* Find Ellipse Coefficents such that
1043         A*u^2 + B*u*v + C*v^2 = F
1044      With u,v relative to point around which we are resampling.
1045      And the given scaling dx,dy vectors in u,v space
1046          du/dx,dv/dx   and  du/dy,dv/dy
1047   */
1048 #if EWA
1049   /* Direct conversion of derivatives into elliptical coefficients
1050      However when magnifying images, the scaling vectors will be small
1051      resulting in a ellipse that is too small to sample properly.
1052      As such we need to clamp the major/minor axis to a minumum of 1.0
1053      to prevent it getting too small.
1054   */
1055 #if EWA_CLAMP
1056   { double major_mag,
1057            minor_mag,
1058            major_x,
1059            major_y,
1060            minor_x,
1061            minor_y;
1062
1063   ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1064                 &major_x, &major_y, &minor_x, &minor_y);
1065   major_x *= major_mag;  major_y *= major_mag;
1066   minor_x *= minor_mag;  minor_y *= minor_mag;
1067 #if DEBUG_ELLIPSE
1068   (void) FormatLocaleFile(stderr, "major_x=%lf; major_y=%lf;  minor_x=%lf; minor_y=%lf;\n",
1069         major_x, major_y, minor_x, minor_y);
1070 #endif
1071   A = major_y*major_y+minor_y*minor_y;
1072   B = -2.0*(major_x*major_y+minor_x*minor_y);
1073   C = major_x*major_x+minor_x*minor_x;
1074   F = major_mag*minor_mag;
1075   F *= F; /* square it */
1076   }
1077 #else /* raw unclamped EWA */
1078   A = dvx*dvx+dvy*dvy;
1079   B = -2.0*(dux*dvx+duy*dvy);
1080   C = dux*dux+duy*duy;
1081   F = dux*dvy-duy*dvx;
1082   F *= F; /* square it */
1083 #endif /* EWA_CLAMP */
1084
1085 #else /* HQ_EWA */
1086   /*
1087     This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
1088     thesis, which adds a unit circle to the elliptical area so as to do both
1089     Reconstruction and Prefiltering of the pixels in the resampling.  It also
1090     means it is always likely to have at least 4 pixels within the area of the
1091     ellipse, for weighted averaging.  No scaling will result with F == 4.0 and
1092     a circle of radius 2.0, and F smaller than this means magnification is
1093     being used.
1094
1095     NOTE: This method produces a very blury result at near unity scale while
1096     producing perfect results for strong minitification and magnifications.
1097
1098     However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
1099   */
1100   A = dvx*dvx+dvy*dvy+1;
1101   B = -2.0*(dux*dvx+duy*dvy);
1102   C = dux*dux+duy*duy+1;
1103   F = A*C - B*B/4;
1104 #endif
1105
1106 #if DEBUG_ELLIPSE
1107   (void) FormatLocaleFile(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1108
1109   /* Figure out the various information directly about the ellipse.
1110      This information currently not needed at this time, but may be
1111      needed later for better limit determination.
1112
1113      It is also good to have as a record for future debugging
1114   */
1115   { double alpha, beta, gamma, Major, Minor;
1116     double Eccentricity, Ellipse_Area, Ellipse_Angle;
1117
1118     alpha = A+C;
1119     beta  = A-C;
1120     gamma = sqrt(beta*beta + B*B );
1121
1122     if ( alpha - gamma <= MagickEpsilon )
1123       Major = MagickHuge;
1124     else
1125       Major = sqrt(2*F/(alpha - gamma));
1126     Minor = sqrt(2*F/(alpha + gamma));
1127
1128     (void) FormatLocaleFile(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
1129
1130     /* other information about ellipse include... */
1131     Eccentricity = Major/Minor;
1132     Ellipse_Area = MagickPI*Major*Minor;
1133     Ellipse_Angle = atan2(B, A-C);
1134
1135     (void) FormatLocaleFile(stderr, "# Angle=%lf   Area=%lf\n",
1136          RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1137   }
1138 #endif
1139
1140   /* If one or both of the scaling vectors is impossibly large
1141      (producing a very large raw F value), we may as well not bother
1142      doing any form of resampling since resampled area is very large.
1143      In this case some alternative means of pixel sampling, such as
1144      the average of the whole image is needed to get a reasonable
1145      result. Calculate only as needed.
1146   */
1147   if ( (4*A*C - B*B) > MagickHuge ) {
1148     resample_filter->limit_reached = MagickTrue;
1149     return;
1150   }
1151
1152   /* Scale ellipse to match the filters support
1153      (that is, multiply F by the square of the support).
1154   */
1155   F *= resample_filter->support;
1156   F *= resample_filter->support;
1157
1158   /* Orthogonal bounds of the ellipse */
1159   resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
1160   resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
1161
1162   /* Horizontally aligned parallelogram fitted to Ellipse */
1163   resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
1164   resample_filter->slope = -B/(2.0*A); /* Reciprocal slope of the parallelogram */
1165
1166 #if DEBUG_ELLIPSE
1167   (void) FormatLocaleFile(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1168            resample_filter->Ulimit, resample_filter->Vlimit,
1169            resample_filter->Uwidth, resample_filter->slope );
1170 #endif
1171
1172   /* Check the absolute area of the parallelogram involved.
1173    * This limit needs more work, as it is too slow for larger images
1174    * with tiled views of the horizon.
1175   */
1176   if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1177          > (4.0*resample_filter->image_area)) {
1178     resample_filter->limit_reached = MagickTrue;
1179     return;
1180   }
1181
1182   /* Scale ellipse formula to directly index the Filter Lookup Table */
1183   { register double scale;
1184 #if FILTER_LUT
1185     /* scale so that F = WLUT_WIDTH; -- hardcoded */
1186     scale = (double)WLUT_WIDTH/F;
1187 #else
1188     /* scale so that F = resample_filter->F (support^2) */
1189     scale = resample_filter->F/F;
1190 #endif
1191     resample_filter->A = A*scale;
1192     resample_filter->B = B*scale;
1193     resample_filter->C = C*scale;
1194   }
1195 }
1196 \f
1197 /*
1198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1199 %                                                                             %
1200 %                                                                             %
1201 %                                                                             %
1202 %   S e t R e s a m p l e F i l t e r                                         %
1203 %                                                                             %
1204 %                                                                             %
1205 %                                                                             %
1206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1207 %
1208 %  SetResampleFilter() set the resampling filter lookup table based on a
1209 %  specific filter.  Note that the filter is used as a radial filter not as a
1210 %  two pass othogonally aligned resampling filter.
1211 %
1212 %  The default Filter, is Gaussian, which is the standard filter used by the
1213 %  original paper on the Elliptical Weighted Everage Algorithm. However other
1214 %  filters can also be used.
1215 %
1216 %  The format of the SetResampleFilter method is:
1217 %
1218 %    void SetResampleFilter(ResampleFilter *resample_filter,
1219 %      const FilterTypes filter,const double blur)
1220 %
1221 %  A description of each parameter follows:
1222 %
1223 %    o resample_filter: resampling resample_filterrmation structure
1224 %
1225 %    o filter: the resize filter for elliptical weighting LUT
1226 %
1227 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1228 %
1229 */
1230 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1231   const FilterTypes filter,const double blur)
1232 {
1233   ResizeFilter
1234      *resize_filter;
1235
1236   assert(resample_filter != (ResampleFilter *) NULL);
1237   assert(resample_filter->signature == MagickSignature);
1238
1239   resample_filter->do_interpolate = MagickFalse;
1240   resample_filter->filter = filter;
1241
1242   if ( filter == PointFilter )
1243     {
1244       resample_filter->do_interpolate = MagickTrue;
1245       return;  /* EWA turned off - nothing more to do */
1246     }
1247
1248   /* Set a default cylindrical filter of a 'low blur' Jinc windowed Jinc */
1249   if ( filter == UndefinedFilter )
1250     resample_filter->filter = RobidouxFilter;
1251
1252   resize_filter = AcquireResizeFilter(resample_filter->image,
1253        resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1254   if (resize_filter == (ResizeFilter *) NULL)
1255     {
1256       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1257            ModuleError, "UnableToSetFilteringValue",
1258            "Fall back to default EWA gaussian filter");
1259       resample_filter->filter = PointFilter;
1260     }
1261
1262   /* Get the practical working support for the filter,
1263    * after any API call blur factors have been accoded for.
1264    */
1265 #if EWA
1266   resample_filter->support = GetResizeFilterSupport(resize_filter);
1267 #else
1268   resample_filter->support = 2.0;  /* fixed support size for HQ-EWA */
1269 #endif
1270
1271 #if FILTER_LUT
1272   /* Fill the LUT with the weights from the selected filter function */
1273   { register int
1274        Q;
1275     double
1276        r_scale;
1277     /* Scale radius so the filter LUT covers the full support range */
1278     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1279     for(Q=0; Q<WLUT_WIDTH; Q++)
1280       resample_filter->filter_lut[Q] = (double)
1281            GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1282
1283     /* finished with the resize filter */
1284     resize_filter = DestroyResizeFilter(resize_filter);
1285   }
1286 #else
1287   /* save the filter and the scaled ellipse bounds needed for filter */
1288   resample_filter->filter_def = resize_filter;
1289   resample_filter->F = resample_filter->support*resample_filter->support;
1290 #endif
1291
1292   /*
1293     Adjust the scaling of the default unit circle
1294     This assumes that any real scaling changes will always
1295     take place AFTER the filter method has been initialized.
1296   */
1297   ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1298
1299 #if 0
1300   /* This is old code kept as a reference only.  It is very wrong,
1301      and I don't understand exactly what it was attempting to do.
1302   */
1303   /*
1304     Create Normal Gaussian 2D Filter Weighted Lookup Table.
1305     A normal EWA guassual lookup would use   exp(Q*ALPHA)
1306     where  Q = distance squared from 0.0 (center) to 1.0 (edge)
1307     and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
1308     The table is of length 1024, and equates to support radius of 2.0
1309     thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
1310
1311     The above came from some reference code provided by Fred Weinhaus
1312     and seems to have been a guess that was appropriate for its use
1313     in a 3d perspective landscape mapping program.
1314   */
1315   r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1316   for(Q=0; Q<WLUT_WIDTH; Q++)
1317     resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1318   resample_filter->support = WLUT_WIDTH;
1319   break;
1320 #endif
1321
1322 #if FILTER_LUT
1323 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1324   #pragma omp single
1325 #endif
1326   { register int
1327        Q;
1328     double
1329        r_scale;
1330
1331     /* Scale radius so the filter LUT covers the full support range */
1332     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1333     if (IsMagickTrue(GetImageArtifact(resample_filter->image,"resample:verbose")) )
1334       {
1335         /* Debug output of the filter weighting LUT
1336           Gnuplot the LUT with hoizontal adjusted to 'r' using...
1337             plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1338           The filter values is normalized for comparision
1339         */
1340         printf("#\n");
1341         printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH);
1342         printf("#\n");
1343         printf("# Note: values in table are using a squared radius lookup.\n");
1344         printf("# And the whole table represents the filters support.\n");
1345         printf("\n"); /* generates a 'break' in gnuplot if multiple outputs */
1346         for(Q=0; Q<WLUT_WIDTH; Q++)
1347           printf("%8.*g %.*g\n",
1348                GetMagickPrecision(),sqrt((double)Q)*r_scale,
1349                GetMagickPrecision(),resample_filter->filter_lut[Q] );
1350       }
1351     /* output the above once only for each image, and each setting */
1352     (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
1353   }
1354 #endif /* FILTER_LUT */
1355   return;
1356 }
1357 \f
1358 /*
1359 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1360 %                                                                             %
1361 %                                                                             %
1362 %                                                                             %
1363 %   S e t R e s a m p l e F i l t e r I n t e r p o l a t e M e t h o d       %
1364 %                                                                             %
1365 %                                                                             %
1366 %                                                                             %
1367 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1368 %
1369 %  SetResampleFilterInterpolateMethod() sets the resample filter interpolation
1370 %  method.
1371 %
1372 %  The format of the SetResampleFilterInterpolateMethod method is:
1373 %
1374 %      MagickBooleanType SetResampleFilterInterpolateMethod(
1375 %        ResampleFilter *resample_filter,const InterpolateMethod method)
1376 %
1377 %  A description of each parameter follows:
1378 %
1379 %    o resample_filter: the resample filter.
1380 %
1381 %    o method: the interpolation method.
1382 %
1383 */
1384 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1385   ResampleFilter *resample_filter,const PixelInterpolateMethod method)
1386 {
1387   assert(resample_filter != (ResampleFilter *) NULL);
1388   assert(resample_filter->signature == MagickSignature);
1389   assert(resample_filter->image != (Image *) NULL);
1390   if (resample_filter->debug != MagickFalse)
1391     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1392       resample_filter->image->filename);
1393   resample_filter->interpolate=method;
1394   return(MagickTrue);
1395 }
1396 \f
1397 /*
1398 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1399 %                                                                             %
1400 %                                                                             %
1401 %                                                                             %
1402 %   S e t R e s a m p l e F i l t e r V i r t u a l P i x e l M e t h o d     %
1403 %                                                                             %
1404 %                                                                             %
1405 %                                                                             %
1406 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1407 %
1408 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1409 %  associated with the specified resample filter.
1410 %
1411 %  The format of the SetResampleFilterVirtualPixelMethod method is:
1412 %
1413 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
1414 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
1415 %
1416 %  A description of each parameter follows:
1417 %
1418 %    o resample_filter: the resample filter.
1419 %
1420 %    o method: the virtual pixel method.
1421 %
1422 */
1423 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1424   ResampleFilter *resample_filter,const VirtualPixelMethod method)
1425 {
1426   assert(resample_filter != (ResampleFilter *) NULL);
1427   assert(resample_filter->signature == MagickSignature);
1428   assert(resample_filter->image != (Image *) NULL);
1429   if (resample_filter->debug != MagickFalse)
1430     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1431       resample_filter->image->filename);
1432   resample_filter->virtual_pixel=method;
1433   if (method != UndefinedVirtualPixelMethod)
1434     (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1435   return(MagickTrue);
1436 }