43 #include "magick/studio.h"
44 #include "magick/artifact.h"
45 #include "magick/color-private.h"
46 #include "magick/cache.h"
47 #include "magick/draw.h"
48 #include "magick/exception-private.h"
49 #include "magick/gem.h"
50 #include "magick/image.h"
51 #include "magick/image-private.h"
52 #include "magick/log.h"
53 #include "magick/magick.h"
54 #include "magick/memory_.h"
55 #include "magick/pixel.h"
56 #include "magick/pixel-private.h"
57 #include "magick/quantum.h"
58 #include "magick/random_.h"
59 #include "magick/resample.h"
60 #include "magick/resize.h"
61 #include "magick/resize-private.h"
62 #include "magick/resource_.h"
63 #include "magick/transform.h"
64 #include "magick/signature-private.h"
65 #include "magick/token.h"
66 #include "magick/utility.h"
67 #include "magick/option.h"
80 #define DEBUG_ELLIPSE 0
81 #define DEBUG_HIT_MISS 0
82 #define DEBUG_NO_PIXEL_HIT 0
85 #define WLUT_WIDTH 1024
109 InterpolatePixelMethod
130 Vlimit, Ulimit, Uwidth, slope;
135 filter_lut[WLUT_WIDTH];
212 assert(image != (
Image *) NULL);
213 assert(image->signature == MagickCoreSignature);
215 assert(exception->signature == MagickCoreSignature);
216 if (IsEventLogging() != MagickFalse)
217 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
220 sizeof(*resample_filter));
222 ThrowFatalException(ResourceLimitFatalError,
"MemoryAllocationFailed");
223 (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,exception);
229 resample_filter->debug=IsEventLogging();
230 resample_filter->signature=MagickCoreSignature;
232 resample_filter->image_area=(ssize_t) (image->columns*image->rows);
233 resample_filter->average_defined = MagickFalse;
236 SetResampleFilter(resample_filter, image->filter, image->blur);
237 (void) SetResampleFilterInterpolateMethod(resample_filter,
239 (void) SetResampleFilterVirtualPixelMethod(resample_filter,
240 GetImageVirtualPixelMethod(image));
242 return(resample_filter);
273 assert(resample_filter->signature == MagickCoreSignature);
274 assert(resample_filter->image != (
Image *) NULL);
275 if (IsEventLogging() != MagickFalse)
276 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
277 resample_filter->image->filename);
278 resample_filter->view=DestroyCacheView(resample_filter->view);
279 resample_filter->image=DestroyImage(resample_filter->image);
281 resample_filter->filter_def=DestroyResizeFilter(resample_filter->filter_def);
283 resample_filter->signature=(~MagickCoreSignature);
284 resample_filter=(
ResampleFilter *) RelinquishMagickMemory(resample_filter);
285 return(resample_filter);
319 MagickExport MagickBooleanType ResamplePixelColor(
326 ssize_t u,v, v1, v2, uw, hit;
329 double divisor_c,divisor_m;
332 const IndexPacket *indexes;
334 assert(resample_filter->signature == MagickCoreSignature);
338 if ( resample_filter->do_interpolate ) {
339 status=InterpolateMagickPixelPacket(resample_filter->image,
340 resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
341 resample_filter->exception);
346 (void) FormatLocaleFile(stderr,
"u0=%lf; v0=%lf;\n", u0, v0);
359 switch ( resample_filter->virtual_pixel ) {
360 case BackgroundVirtualPixelMethod:
361 case ConstantVirtualPixelMethod:
362 case TransparentVirtualPixelMethod:
363 case BlackVirtualPixelMethod:
364 case GrayVirtualPixelMethod:
365 case WhiteVirtualPixelMethod:
366 case MaskVirtualPixelMethod:
367 if ( resample_filter->limit_reached
368 || u0 + resample_filter->Ulimit < 0.0
369 || u0 - resample_filter->Ulimit > (
double) resample_filter->image->columns-1.0
370 || v0 + resample_filter->Vlimit < 0.0
371 || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows-1.0
376 case UndefinedVirtualPixelMethod:
377 case EdgeVirtualPixelMethod:
378 if ( ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
379 || ( u0 + resample_filter->Ulimit < 0.0
380 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows-1.0 )
381 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns-1.0
382 && v0 + resample_filter->Vlimit < 0.0 )
383 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns-1.0
384 && v0 - resample_filter->Vlimit > (
double) resample_filter->image->rows-1.0 )
388 case HorizontalTileVirtualPixelMethod:
389 if ( v0 + resample_filter->Vlimit < 0.0
390 || v0 - resample_filter->Vlimit > (
double) resample_filter->image->rows-1.0
394 case VerticalTileVirtualPixelMethod:
395 if ( u0 + resample_filter->Ulimit < 0.0
396 || u0 - resample_filter->Ulimit > (
double) resample_filter->image->columns-1.0
400 case DitherVirtualPixelMethod:
401 if ( ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
402 || ( u0 + resample_filter->Ulimit < -32.0
403 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+31.0 )
404 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+31.0
405 && v0 + resample_filter->Vlimit < -32.0 )
406 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+31.0
407 && v0 - resample_filter->Vlimit > (
double) resample_filter->image->rows+31.0 )
411 case TileVirtualPixelMethod:
412 case MirrorVirtualPixelMethod:
413 case RandomVirtualPixelMethod:
414 case HorizontalTileEdgeVirtualPixelMethod:
415 case VerticalTileEdgeVirtualPixelMethod:
416 case CheckerTileVirtualPixelMethod:
426 status=InterpolateMagickPixelPacket(resample_filter->image,
427 resample_filter->view,IntegerInterpolatePixel,u0,v0,pixel,
428 resample_filter->exception);
435 if ( resample_filter->limit_reached ) {
436 switch ( resample_filter->virtual_pixel ) {
445 case UndefinedVirtualPixelMethod:
446 case EdgeVirtualPixelMethod:
447 case DitherVirtualPixelMethod:
448 case HorizontalTileEdgeVirtualPixelMethod:
449 case VerticalTileEdgeVirtualPixelMethod:
456 status=InterpolateMagickPixelPacket(resample_filter->image,
457 resample_filter->view,AverageInterpolatePixel,u0,v0,pixel,
458 resample_filter->exception);
460 case HorizontalTileVirtualPixelMethod:
461 case VerticalTileVirtualPixelMethod:
463 status=InterpolateMagickPixelPacket(resample_filter->image,
464 resample_filter->view,IntegerInterpolatePixel,-1.0,-1.0,pixel,
465 resample_filter->exception);
467 case TileVirtualPixelMethod:
468 case MirrorVirtualPixelMethod:
469 case RandomVirtualPixelMethod:
470 case CheckerTileVirtualPixelMethod:
473 if ( resample_filter->average_defined == MagickFalse ) {
481 &resample_filter->average_pixel);
482 resample_filter->average_defined=MagickTrue;
485 average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
486 resample_filter->exception);
487 if (average_image == (
Image *) NULL)
489 *pixel=resample_filter->average_pixel;
492 average_view=AcquireVirtualCacheView(average_image,
493 &average_image->exception);
494 pixels=(
PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
495 resample_filter->exception);
497 average_view=DestroyCacheView(average_view);
498 average_image=DestroyImage(average_image);
499 *pixel=resample_filter->average_pixel;
502 indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
503 SetMagickPixelPacket(resample_filter->image,pixels,indexes,
504 &(resample_filter->average_pixel));
505 average_view=DestroyCacheView(average_view);
506 average_image=DestroyImage(average_image);
508 if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod )
514 weight = QuantumScale*((MagickRealType)(QuantumRange-
515 resample_filter->average_pixel.opacity));
516 resample_filter->average_pixel.red *= weight;
517 resample_filter->average_pixel.green *= weight;
518 resample_filter->average_pixel.blue *= weight;
522 weight = QuantumScale*((MagickRealType)(QuantumRange-
523 resample_filter->image->background_color.opacity));
524 resample_filter->average_pixel.red +=
525 weight*resample_filter->image->background_color.red;
526 resample_filter->average_pixel.green +=
527 weight*resample_filter->image->background_color.green;
528 resample_filter->average_pixel.blue +=
529 weight*resample_filter->image->background_color.blue;
530 resample_filter->average_pixel.opacity +=
531 resample_filter->image->background_color.opacity;
535 resample_filter->average_pixel.red /= divisor_c;
536 resample_filter->average_pixel.green /= divisor_c;
537 resample_filter->average_pixel.blue /= divisor_c;
538 resample_filter->average_pixel.opacity /= 2;
542 *pixel=resample_filter->average_pixel;
554 pixel->red = pixel->green = pixel->blue = 0.0;
555 if (pixel->matte != MagickFalse) pixel->opacity = 0.0;
556 if (pixel->colorspace == CMYKColorspace) pixel->index = 0.0;
562 v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit);
563 v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
566 u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
567 uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
570 (void) FormatLocaleFile(stderr,
"v1=%ld; v2=%ld\n", (
long)v1, (long)v2);
571 (void) FormatLocaleFile(stderr,
"u1=%ld; uw=%ld\n", (
long)u1, (long)uw);
573 # define DEBUG_HIT_MISS 0
580 DDQ = 2*resample_filter->A;
581 for( v=v1; v<=v2; v++ ) {
584 (void) FormatLocaleFile(stderr,
"# scan line from pixel %ld, %ld\n", (
long)uu, (long)v);
586 u = (ssize_t)ceil(u1);
587 u1 += resample_filter->slope;
595 Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
596 DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
599 pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(
size_t) uw,
600 1,resample_filter->exception);
603 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
606 for( u=0; u<uw; u++ ) {
610 if ( Q < (
double)WLUT_WIDTH ) {
611 weight = resample_filter->filter_lut[(int)Q];
614 if ( Q < (
double)resample_filter->F ) {
615 weight = GetResizeFilterWeight(resample_filter->filter_def,
619 if (pixel->matte != MagickFalse)
620 pixel->opacity += weight*pixels->opacity;
623 if (pixel->matte != MagickFalse)
624 weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
625 pixel->red += weight*pixels->red;
626 pixel->green += weight*pixels->green;
627 pixel->blue += weight*pixels->blue;
628 if (pixel->colorspace == CMYKColorspace)
629 pixel->index += weight*(*indexes);
635 (void) FormatLocaleFile(stderr,
"set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
636 (
long)uu-.1,(double)v-.1,(
long)uu+.1,(long)v+.1);
637 (void) FormatLocaleFile(stderr,
"set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
638 (
long)uu+.1,(double)v-.1,(
long)uu-.1,(long)v+.1);
640 (void) FormatLocaleFile(stderr,
"set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
641 (
long)uu-.1,(double)v-.1,(
long)uu+.1,(long)v+.1);
642 (void) FormatLocaleFile(stderr,
"set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
643 (
long)uu+.1,(double)v-.1,(
long)uu-.1,(long)v+.1);
656 (void) FormatLocaleFile(stderr,
"Hit=%ld; Total=%ld;\n", (
long)hit, (long)uw*(v2-v1) );
662 if ( hit == 0 || divisor_m <= MagickEpsilon || divisor_c <= MagickEpsilon ) {
665 #if DEBUG_NO_PIXEL_HIT
666 pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
667 pixel->red = QuantumRange;
669 status=InterpolateMagickPixelPacket(resample_filter->image,
670 resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
671 resample_filter->exception);
679 divisor_m = 1.0/divisor_m;
680 if (pixel->matte != MagickFalse)
681 pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
682 divisor_c = 1.0/divisor_c;
683 pixel->red = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
684 pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
685 pixel->blue = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
686 if (pixel->colorspace == CMYKColorspace)
687 pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
721 static inline void ClampUpAxes(
const double dux,
727 double *major_unit_x,
728 double *major_unit_y,
729 double *minor_unit_x,
730 double *minor_unit_y)
894 const double a = dux;
895 const double b = duy;
896 const double c = dvx;
897 const double d = dvy;
902 const double aa = a*a;
903 const double bb = b*b;
904 const double cc = c*c;
905 const double dd = d*d;
909 const double n11 = aa+bb;
910 const double n12 = a*c+b*d;
911 const double n21 = n12;
912 const double n22 = cc+dd;
913 const double det = a*d-b*c;
914 const double twice_det = det+det;
915 const double frobenius_squared = n11+n22;
916 const double discriminant =
917 (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
923 const double sqrt_discriminant =
924 sqrt(discriminant > 0.0 ? discriminant : 0.0);
935 const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
941 const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
942 const double s1s1minusn11 = s1s1-n11;
943 const double s1s1minusn22 = s1s1-n22;
951 const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
952 const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
962 const double temp_u11 =
963 ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
964 const double temp_u21 =
965 ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
966 const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
971 const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
972 const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
976 *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
977 *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
983 *minor_unit_x = -u21;
1050 MagickExport
void ScaleResampleFilter(
ResampleFilter *resample_filter,
1051 const double dux,
const double duy,
const double dvx,
const double dvy)
1056 assert(resample_filter->signature == MagickCoreSignature);
1058 resample_filter->limit_reached = MagickFalse;
1061 if ( resample_filter->filter == PointFilter )
1065 (void) FormatLocaleFile(stderr,
"# -----\n" );
1066 (void) FormatLocaleFile(stderr,
"dux=%lf; dvx=%lf; duy=%lf; dvy=%lf;\n",
1067 dux, dvx, duy, dvy);
1091 ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1092 &major_x, &major_y, &minor_x, &minor_y);
1093 major_x *= major_mag; major_y *= major_mag;
1094 minor_x *= minor_mag; minor_y *= minor_mag;
1096 (void) FormatLocaleFile(stderr,
"major_x=%lf; major_y=%lf; minor_x=%lf; minor_y=%lf;\n",
1097 major_x, major_y, minor_x, minor_y);
1099 A = major_y*major_y+minor_y*minor_y;
1100 B = -2.0*(major_x*major_y+minor_x*minor_y);
1101 C = major_x*major_x+minor_x*minor_x;
1102 F = major_mag*minor_mag;
1106 A = dvx*dvx+dvy*dvy;
1107 B = -2.0*(dux*dvx+duy*dvy);
1108 C = dux*dux+duy*duy;
1109 F = dux*dvy-duy*dvx;
1128 A = dvx*dvx+dvy*dvy+1;
1129 B = -2.0*(dux*dvx+duy*dvy);
1130 C = dux*dux+duy*duy+1;
1135 (void) FormatLocaleFile(stderr,
"A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1143 {
double alpha, beta, gamma, Major, Minor;
1144 double Eccentricity, Ellipse_Area, Ellipse_Angle;
1148 gamma = sqrt(beta*beta + B*B );
1150 if ( alpha - gamma <= MagickEpsilon )
1151 Major= MagickMaximumValue;
1153 Major= sqrt(2*F/(alpha - gamma));
1154 Minor = sqrt(2*F/(alpha + gamma));
1156 (void) FormatLocaleFile(stderr,
"# Major=%lf; Minor=%lf\n", Major, Minor );
1159 Eccentricity = Major/Minor;
1160 Ellipse_Area = MagickPI*Major*Minor;
1161 Ellipse_Angle = atan2(B, A-C);
1163 (void) FormatLocaleFile(stderr,
"# Angle=%lf Area=%lf\n",
1164 (
double) RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1175 if ( (4*A*C - B*B) > MagickMaximumValue ) {
1176 resample_filter->limit_reached = MagickTrue;
1184 F *= resample_filter->support;
1185 F *= resample_filter->support;
1188 resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
1189 resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
1192 resample_filter->Uwidth = sqrt(F/A);
1193 resample_filter->slope = -B/(2.0*A);
1196 (void) FormatLocaleFile(stderr,
"Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1197 resample_filter->Ulimit, resample_filter->Vlimit,
1198 resample_filter->Uwidth, resample_filter->slope );
1205 if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1206 > (4.0*resample_filter->image_area)) {
1207 resample_filter->limit_reached = MagickTrue;
1215 scale=(double) WLUT_WIDTH*PerceptibleReciprocal(F);
1218 scale=resample_filter->F*PerceptibleReciprocal(F);
1220 resample_filter->A = A*scale;
1221 resample_filter->B = B*scale;
1222 resample_filter->C = C*scale;
1255 MagickExport
void SetResampleFilter(
ResampleFilter *resample_filter,
1256 const FilterTypes filter,
const double blur)
1262 assert(resample_filter->signature == MagickCoreSignature);
1264 resample_filter->do_interpolate = MagickFalse;
1265 resample_filter->filter = filter;
1268 if ( filter == UndefinedFilter )
1269 resample_filter->filter = RobidouxFilter;
1271 if ( resample_filter->filter == PointFilter ) {
1272 resample_filter->do_interpolate = MagickTrue;
1276 resize_filter = AcquireResizeFilter(resample_filter->image,
1277 resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1279 (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1280 ModuleError,
"UnableToSetFilteringValue",
1281 "Fall back to Interpolated 'Point' filter");
1282 resample_filter->filter = PointFilter;
1283 resample_filter->do_interpolate = MagickTrue;
1291 resample_filter->support = GetResizeFilterSupport(resize_filter);
1293 resample_filter->support = 2.0;
1304 r_scale = resample_filter->support*sqrt(1.0/(
double)WLUT_WIDTH);
1305 for(Q=0; Q<WLUT_WIDTH; Q++)
1306 resample_filter->filter_lut[Q] = (
double)
1307 GetResizeFilterWeight(resize_filter,sqrt((
double)Q)*r_scale);
1310 resize_filter = DestroyResizeFilter(resize_filter);
1314 resample_filter->filter_def = resize_filter;
1315 resample_filter->F = resample_filter->support*resample_filter->support;
1323 ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1339 r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1340 for(Q=0; Q<WLUT_WIDTH; Q++)
1341 resample_filter->filter_lut[Q] = exp((
double)Q*r_scale);
1342 resample_filter->support = WLUT_WIDTH;
1346 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1350 if (IsMagickTrue(GetImageArtifact(resample_filter->image,
1351 "resample:verbose")) )
1364 printf(
"# Resampling Filter LUT (%d values) for '%s' filter\n",
1365 WLUT_WIDTH, CommandOptionToMnemonic(MagickFilterOptions,
1366 resample_filter->filter) );
1368 printf(
"# Note: values in table are using a squared radius lookup.\n");
1369 printf(
"# As such its distribution is not uniform.\n");
1371 printf(
"# The X value is the support distance for the Y weight\n");
1372 printf(
"# so you can use gnuplot to plot this cylindrical filter\n");
1373 printf(
"# plot [0:2][-.2:1] \"lut.dat\" with lines\n");
1377 r_scale = resample_filter->support*sqrt(1.0/(
double)WLUT_WIDTH);
1378 for(Q=0; Q<WLUT_WIDTH; Q++)
1379 printf(
"%8.*g %.*g\n",
1380 GetMagickPrecision(),sqrt((
double)Q)*r_scale,
1381 GetMagickPrecision(),resample_filter->filter_lut[Q] );
1418 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1419 ResampleFilter *resample_filter,
const InterpolatePixelMethod method)
1422 assert(resample_filter->signature == MagickCoreSignature);
1423 assert(resample_filter->image != (
Image *) NULL);
1424 if (IsEventLogging() != MagickFalse)
1425 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
1426 resample_filter->image->filename);
1427 resample_filter->interpolate=method;
1457 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1461 assert(resample_filter->signature == MagickCoreSignature);
1462 assert(resample_filter->image != (
Image *) NULL);
1463 if (IsEventLogging() != MagickFalse)
1464 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",
1465 resample_filter->image->filename);
1466 resample_filter->virtual_pixel=method;
1467 if (method != UndefinedVirtualPixelMethod)
1468 (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);