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