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