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