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