MagickCore 6.9.13-26
Convert, Edit, Or Compose Bitmap Images
Loading...
Searching...
No Matches
compare.c
1/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3% %
4% %
5% %
6% CCCC OOO M M PPPP AAA RRRR EEEEE %
7% C O O MM MM P P A A R R E %
8% C O O M M M PPPP AAAAA RRRR EEE %
9% C O O M M P A A R R E %
10% CCCC OOO M M P A A R R EEEEE %
11% %
12% %
13% MagickCore Image Comparison Methods %
14% %
15% Software Design %
16% Cristy %
17% December 2003 %
18% %
19% %
20% Copyright 1999 ImageMagick Studio LLC, a non-profit organization dedicated %
21% to making software imaging solutions freely available. %
22% %
23% You may not use this file except in compliance with the License. You may %
24% obtain a copy of the License at %
25% %
26% https://imagemagick.org/script/license.php %
27% %
28% Unless required by applicable law or agreed to in writing, software %
29% distributed under the License is distributed on an "AS IS" BASIS, %
30% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31% See the License for the specific language governing permissions and %
32% limitations under the License. %
33% %
34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35%
36%
37%
38*/
39
40/*
41 Include declarations.
42*/
43#include "magick/studio.h"
44#include "magick/artifact.h"
45#include "magick/attribute.h"
46#include "magick/cache-view.h"
47#include "magick/channel.h"
48#include "magick/client.h"
49#include "magick/color.h"
50#include "magick/color-private.h"
51#include "magick/colorspace.h"
52#include "magick/colorspace-private.h"
53#include "magick/compare.h"
54#include "magick/compare-private.h"
55#include "magick/composite-private.h"
56#include "magick/constitute.h"
57#include "magick/exception-private.h"
58#include "magick/geometry.h"
59#include "magick/image-private.h"
60#include "magick/list.h"
61#include "magick/log.h"
62#include "magick/memory_.h"
63#include "magick/monitor.h"
64#include "magick/monitor-private.h"
65#include "magick/option.h"
66#include "magick/pixel-private.h"
67#include "magick/property.h"
68#include "magick/resource_.h"
69#include "magick/statistic-private.h"
70#include "magick/string_.h"
71#include "magick/string-private.h"
72#include "magick/statistic.h"
73#include "magick/thread-private.h"
74#include "magick/transform.h"
75#include "magick/utility.h"
76#include "magick/version.h"
77
78/*
79%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80% %
81% %
82% %
83% C o m p a r e I m a g e C h a n n e l s %
84% %
85% %
86% %
87%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
88%
89% CompareImageChannels() compares one or more image channels of an image
90% to a reconstructed image and returns the difference image.
91%
92% The format of the CompareImageChannels method is:
93%
94% Image *CompareImageChannels(const Image *image,
95% const Image *reconstruct_image,const ChannelType channel,
96% const MetricType metric,double *distortion,ExceptionInfo *exception)
97%
98% A description of each parameter follows:
99%
100% o image: the image.
101%
102% o reconstruct_image: the reconstruct image.
103%
104% o channel: the channel.
105%
106% o metric: the metric.
107%
108% o distortion: the computed distortion between the images.
109%
110% o exception: return any errors or warnings in this structure.
111%
112*/
113
114MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
115 const MetricType metric,double *distortion,ExceptionInfo *exception)
116{
117 Image
118 *highlight_image;
119
120 highlight_image=CompareImageChannels(image,reconstruct_image,
121 CompositeChannels,metric,distortion,exception);
122 return(highlight_image);
123}
124
125static size_t GetNumberChannels(const Image *image,const ChannelType channel)
126{
127 size_t
128 channels;
129
130 channels=0;
131 if ((channel & RedChannel) != 0)
132 channels++;
133 if ((channel & GreenChannel) != 0)
134 channels++;
135 if ((channel & BlueChannel) != 0)
136 channels++;
137 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
138 channels++;
139 if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
140 channels++;
141 return(channels == 0 ? 1UL : channels);
142}
143
144static inline MagickBooleanType ValidateImageMorphology(
145 const Image *magick_restrict image,
146 const Image *magick_restrict reconstruct_image)
147{
148 /*
149 Does the image match the reconstructed image morphology?
150 */
151 if (GetNumberChannels(image,DefaultChannels) !=
152 GetNumberChannels(reconstruct_image,DefaultChannels))
153 return(MagickFalse);
154 return(MagickTrue);
155}
156
157MagickExport Image *CompareImageChannels(Image *image,
158 const Image *reconstruct_image,const ChannelType channel,
159 const MetricType metric,double *distortion,ExceptionInfo *exception)
160{
161 CacheView
162 *highlight_view,
163 *image_view,
164 *reconstruct_view;
165
166 const char
167 *artifact;
168
169 double
170 fuzz;
171
172 Image
173 *clone_image,
174 *difference_image,
175 *highlight_image;
176
177 MagickBooleanType
178 status;
179
180 MagickPixelPacket
181 highlight,
182 lowlight,
183 zero;
184
185 size_t
186 columns,
187 rows;
188
189 ssize_t
190 y;
191
192 assert(image != (Image *) NULL);
193 assert(image->signature == MagickCoreSignature);
194 assert(reconstruct_image != (const Image *) NULL);
195 assert(reconstruct_image->signature == MagickCoreSignature);
196 assert(distortion != (double *) NULL);
197 if (IsEventLogging() != MagickFalse)
198 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
199 *distortion=0.0;
200 if (metric != PerceptualHashErrorMetric)
201 if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
202 ThrowImageException(ImageError,"ImageMorphologyDiffers");
203 status=GetImageChannelDistortion(image,reconstruct_image,channel,metric,
204 distortion,exception);
205 if (status == MagickFalse)
206 return((Image *) NULL);
207 clone_image=CloneImage(image,0,0,MagickTrue,exception);
208 if (clone_image == (Image *) NULL)
209 return((Image *) NULL);
210 (void) SetImageMask(clone_image,(Image *) NULL);
211 difference_image=CloneImage(clone_image,0,0,MagickTrue,exception);
212 clone_image=DestroyImage(clone_image);
213 if (difference_image == (Image *) NULL)
214 return((Image *) NULL);
215 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel);
216 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
217 highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
218 if (highlight_image == (Image *) NULL)
219 {
220 difference_image=DestroyImage(difference_image);
221 return((Image *) NULL);
222 }
223 if (SetImageStorageClass(highlight_image,DirectClass) == MagickFalse)
224 {
225 InheritException(exception,&highlight_image->exception);
226 difference_image=DestroyImage(difference_image);
227 highlight_image=DestroyImage(highlight_image);
228 return((Image *) NULL);
229 }
230 (void) SetImageMask(highlight_image,(Image *) NULL);
231 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel);
232 (void) QueryMagickColor("#f1001ecc",&highlight,exception);
233 artifact=GetImageArtifact(image,"compare:highlight-color");
234 if (artifact != (const char *) NULL)
235 (void) QueryMagickColor(artifact,&highlight,exception);
236 (void) QueryMagickColor("#ffffffcc",&lowlight,exception);
237 artifact=GetImageArtifact(image,"compare:lowlight-color");
238 if (artifact != (const char *) NULL)
239 (void) QueryMagickColor(artifact,&lowlight,exception);
240 if (highlight_image->colorspace == CMYKColorspace)
241 {
242 ConvertRGBToCMYK(&highlight);
243 ConvertRGBToCMYK(&lowlight);
244 }
245 /*
246 Generate difference image.
247 */
248 status=MagickTrue;
249 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
250 GetMagickPixelPacket(image,&zero);
251 image_view=AcquireVirtualCacheView(image,exception);
252 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
253 highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
254#if defined(MAGICKCORE_OPENMP_SUPPORT)
255 #pragma omp parallel for schedule(static) shared(status) \
256 magick_number_threads(image,highlight_image,rows,1)
257#endif
258 for (y=0; y < (ssize_t) rows; y++)
259 {
260 MagickBooleanType
261 sync;
262
263 MagickPixelPacket
264 pixel,
265 reconstruct_pixel;
266
267 const IndexPacket
268 *magick_restrict indexes,
269 *magick_restrict reconstruct_indexes;
270
271 const PixelPacket
272 *magick_restrict p,
273 *magick_restrict q;
274
275 IndexPacket
276 *magick_restrict highlight_indexes;
277
278 PixelPacket
279 *magick_restrict r;
280
281 ssize_t
282 x;
283
284 if (status == MagickFalse)
285 continue;
286 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
287 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
288 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
289 if ((p == (const PixelPacket *) NULL) ||
290 (q == (const PixelPacket *) NULL) || (r == (PixelPacket *) NULL))
291 {
292 status=MagickFalse;
293 continue;
294 }
295 indexes=GetCacheViewVirtualIndexQueue(image_view);
296 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
297 highlight_indexes=GetCacheViewAuthenticIndexQueue(highlight_view);
298 pixel=zero;
299 reconstruct_pixel=zero;
300 for (x=0; x < (ssize_t) columns; x++)
301 {
302 MagickStatusType
303 difference;
304
305 SetMagickPixelPacket(image,p,indexes == (IndexPacket *) NULL ? NULL :
306 indexes+x,&pixel);
307 SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes ==
308 (IndexPacket *) NULL ? NULL : reconstruct_indexes+x,&reconstruct_pixel);
309 difference=MagickFalse;
310 if (channel == CompositeChannels)
311 {
312 if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
313 difference=MagickTrue;
314 }
315 else
316 {
317 double
318 Da,
319 distance,
320 pixel,
321 Sa;
322
323 Sa=QuantumScale*(image->matte != MagickFalse ? (double)
324 GetPixelAlpha(p) : ((double) QuantumRange-(double) OpaqueOpacity));
325 Da=QuantumScale*(image->matte != MagickFalse ? (double)
326 GetPixelAlpha(q) : ((double) QuantumRange-(double) OpaqueOpacity));
327 if ((channel & RedChannel) != 0)
328 {
329 pixel=Sa*(double) GetPixelRed(p)-Da*(double) GetPixelRed(q);
330 distance=pixel*pixel;
331 if (distance > fuzz)
332 difference=MagickTrue;
333 }
334 if ((channel & GreenChannel) != 0)
335 {
336 pixel=Sa*(double) GetPixelGreen(p)-Da*(double) GetPixelGreen(q);
337 distance=pixel*pixel;
338 if (distance > fuzz)
339 difference=MagickTrue;
340 }
341 if ((channel & BlueChannel) != 0)
342 {
343 pixel=Sa*(double) GetPixelBlue(p)-Da*(double) GetPixelBlue(q);
344 distance=pixel*pixel;
345 if (distance > fuzz)
346 difference=MagickTrue;
347 }
348 if (((channel & OpacityChannel) != 0) &&
349 (image->matte != MagickFalse))
350 {
351 pixel=(double) GetPixelOpacity(p)-(double) GetPixelOpacity(q);
352 distance=pixel*pixel;
353 if (distance > fuzz)
354 difference=MagickTrue;
355 }
356 if (((channel & IndexChannel) != 0) &&
357 (image->colorspace == CMYKColorspace))
358 {
359 pixel=Sa*(double) indexes[x]-Da*(double) reconstruct_indexes[x];
360 distance=pixel*pixel;
361 if (distance > fuzz)
362 difference=MagickTrue;
363 }
364 }
365 if (difference != MagickFalse)
366 SetPixelPacket(highlight_image,&highlight,r,highlight_indexes ==
367 (IndexPacket *) NULL ? NULL : highlight_indexes+x);
368 else
369 SetPixelPacket(highlight_image,&lowlight,r,highlight_indexes ==
370 (IndexPacket *) NULL ? NULL : highlight_indexes+x);
371 p++;
372 q++;
373 r++;
374 }
375 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
376 if (sync == MagickFalse)
377 status=MagickFalse;
378 }
379 highlight_view=DestroyCacheView(highlight_view);
380 reconstruct_view=DestroyCacheView(reconstruct_view);
381 image_view=DestroyCacheView(image_view);
382 (void) CompositeImage(difference_image,image->compose,highlight_image,0,0);
383 highlight_image=DestroyImage(highlight_image);
384 if (status == MagickFalse)
385 difference_image=DestroyImage(difference_image);
386 return(difference_image);
387}
388
389/*
390%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
391% %
392% %
393% %
394% G e t I m a g e C h a n n e l D i s t o r t i o n %
395% %
396% %
397% %
398%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
399%
400% GetImageChannelDistortion() compares one or more image channels of an image
401% to a reconstructed image and returns the specified distortion metric.
402%
403% The format of the GetImageChannelDistortion method is:
404%
405% MagickBooleanType GetImageChannelDistortion(const Image *image,
406% const Image *reconstruct_image,const ChannelType channel,
407% const MetricType metric,double *distortion,ExceptionInfo *exception)
408%
409% A description of each parameter follows:
410%
411% o image: the image.
412%
413% o reconstruct_image: the reconstruct image.
414%
415% o channel: the channel.
416%
417% o metric: the metric.
418%
419% o distortion: the computed distortion between the images.
420%
421% o exception: return any errors or warnings in this structure.
422%
423*/
424
425MagickExport MagickBooleanType GetImageDistortion(Image *image,
426 const Image *reconstruct_image,const MetricType metric,double *distortion,
427 ExceptionInfo *exception)
428{
429 MagickBooleanType
430 status;
431
432 status=GetImageChannelDistortion(image,reconstruct_image,CompositeChannels,
433 metric,distortion,exception);
434 return(status);
435}
436
437static MagickBooleanType GetAbsoluteDistortion(const Image *image,
438 const Image *reconstruct_image,const ChannelType channel,double *distortion,
439 ExceptionInfo *exception)
440{
441 CacheView
442 *image_view,
443 *reconstruct_view;
444
445 double
446 area;
447
448 MagickBooleanType
449 status;
450
451 size_t
452 columns,
453 rows;
454
455 ssize_t
456 j,
457 y;
458
459 /*
460 Compute the absolute difference in pixels between two images.
461 */
462 status=MagickTrue;
463 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
464 image_view=AcquireVirtualCacheView(image,exception);
465 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
466#if defined(MAGICKCORE_OPENMP_SUPPORT)
467 #pragma omp parallel for schedule(static) shared(distortion,status) \
468 magick_number_threads(image,image,rows,1)
469#endif
470 for (y=0; y < (ssize_t) rows; y++)
471 {
472 const IndexPacket
473 *magick_restrict indexes,
474 *magick_restrict reconstruct_indexes;
475
476 const PixelPacket
477 *magick_restrict p,
478 *magick_restrict q;
479
480 double
481 channel_distortion[CompositeChannels+1] = { 0.0 };
482
483 ssize_t
484 i,
485 x;
486
487 if (status == MagickFalse)
488 continue;
489 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
490 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
491 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
492 {
493 status=MagickFalse;
494 continue;
495 }
496 indexes=GetCacheViewVirtualIndexQueue(image_view);
497 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
498 (void) memset(channel_distortion,0,sizeof(channel_distortion));
499 for (x=0; x < (ssize_t) columns; x++)
500 {
501 double
502 Da,
503 delta,
504 Sa;
505
506 size_t
507 count = 0;
508
509 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
510 ((double) QuantumRange-(double) OpaqueOpacity));
511 Da=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(q) :
512 ((double) QuantumRange-(double) OpaqueOpacity));
513 if ((channel & RedChannel) != 0)
514 {
515 delta=Sa*(double) GetPixelRed(p)-Da*(double)
516 GetPixelRed(q);
517 if ((delta*delta) >= MagickEpsilon)
518 {
519 channel_distortion[RedChannel]++;
520 count++;
521 }
522 }
523 if ((channel & GreenChannel) != 0)
524 {
525 delta=Sa*(double) GetPixelGreen(p)-Da*(double)
526 GetPixelGreen(q);
527 if ((delta*delta) >= MagickEpsilon)
528 {
529 channel_distortion[GreenChannel]++;
530 count++;
531 }
532 }
533 if ((channel & BlueChannel) != 0)
534 {
535 delta=Sa*(double) GetPixelBlue(p)-Da*(double)
536 GetPixelBlue(q);
537 if ((delta*delta) >= MagickEpsilon)
538 {
539 channel_distortion[BlueChannel]++;
540 count++;
541 }
542 }
543 if (((channel & OpacityChannel) != 0) &&
544 (image->matte != MagickFalse))
545 {
546 delta=(double) GetPixelOpacity(p)-(double)
547 GetPixelOpacity(q);
548 if ((delta*delta) >= MagickEpsilon)
549 {
550 channel_distortion[OpacityChannel]++;
551 count++;
552 }
553 }
554 if (((channel & IndexChannel) != 0) &&
555 (image->colorspace == CMYKColorspace))
556 {
557 delta=Sa*(double) indexes[x]-Da*(double)
558 reconstruct_indexes[x];
559 if ((delta*delta) >= MagickEpsilon)
560 {
561 channel_distortion[IndexChannel]++;
562 count++;
563 }
564 }
565 if (count != 0)
566 channel_distortion[CompositeChannels]++;
567 p++;
568 q++;
569 }
570#if defined(MAGICKCORE_OPENMP_SUPPORT)
571 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
572#endif
573 for (i=0; i <= (ssize_t) CompositeChannels; i++)
574 distortion[i]+=channel_distortion[i];
575 }
576 reconstruct_view=DestroyCacheView(reconstruct_view);
577 image_view=DestroyCacheView(image_view);
578 area=MagickSafeReciprocal((double) columns*rows);
579 for (j=0; j <= CompositeChannels; j++)
580 distortion[j]*=area;
581 return(status);
582}
583
584static MagickBooleanType GetFuzzDistortion(const Image *image,
585 const Image *reconstruct_image,const ChannelType channel,double *distortion,
586 ExceptionInfo *exception)
587{
588 CacheView
589 *image_view,
590 *reconstruct_view;
591
592 double
593 area,
594 fuzz;
595
596 MagickBooleanType
597 status;
598
599 size_t
600 columns,
601 rows;
602
603 ssize_t
604 j,
605 y;
606
607 /*
608 Compute the absolute difference in pixels between two images.
609 */
610 status=MagickTrue;
611 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
612 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
613 image_view=AcquireVirtualCacheView(image,exception);
614 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
615#if defined(MAGICKCORE_OPENMP_SUPPORT)
616 #pragma omp parallel for schedule(static) shared(distortion,status) \
617 magick_number_threads(image,image,rows,1)
618#endif
619 for (y=0; y < (ssize_t) rows; y++)
620 {
621 const IndexPacket
622 *magick_restrict indexes,
623 *magick_restrict reconstruct_indexes;
624
625 const PixelPacket
626 *magick_restrict p,
627 *magick_restrict q;
628
629 double
630 channel_distortion[CompositeChannels+1] = { 0.0 };
631
632 ssize_t
633 i,
634 x;
635
636 if (status == MagickFalse)
637 continue;
638 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
639 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
640 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
641 {
642 status=MagickFalse;
643 continue;
644 }
645 indexes=GetCacheViewVirtualIndexQueue(image_view);
646 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
647 (void) memset(channel_distortion,0,sizeof(channel_distortion));
648 for (x=0; x < (ssize_t) columns; x++)
649 {
650 double
651 Da,
652 delta,
653 Sa;
654
655 size_t
656 count = 0;
657
658 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
659 ((double) QuantumRange-(double) OpaqueOpacity));
660 Da=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(q) :
661 ((double) QuantumRange-(double) OpaqueOpacity));
662 if ((channel & RedChannel) != 0)
663 {
664 delta=Sa*(double) GetPixelRed(p)-Da*(double)
665 GetPixelRed(q);
666 if ((delta*delta) > fuzz)
667 {
668 channel_distortion[RedChannel]++;
669 count++;
670 }
671 }
672 if ((channel & GreenChannel) != 0)
673 {
674 delta=Sa*(double) GetPixelGreen(p)-Da*(double)
675 GetPixelGreen(q);
676 if ((delta*delta) > fuzz)
677 {
678 channel_distortion[GreenChannel]++;
679 count++;
680 }
681 }
682 if ((channel & BlueChannel) != 0)
683 {
684 delta=Sa*(double) GetPixelBlue(p)-Da*(double)
685 GetPixelBlue(q);
686 if ((delta*delta) > fuzz)
687 {
688 channel_distortion[BlueChannel]++;
689 count++;
690 }
691 }
692 if (((channel & OpacityChannel) != 0) &&
693 (image->matte != MagickFalse))
694 {
695 delta=(double) GetPixelOpacity(p)-(double)
696 GetPixelOpacity(q);
697 if ((delta*delta) > fuzz)
698 {
699 channel_distortion[OpacityChannel]++;
700 count++;
701 }
702 }
703 if (((channel & IndexChannel) != 0) &&
704 (image->colorspace == CMYKColorspace))
705 {
706 delta=Sa*(double) indexes[x]-Da*(double)
707 reconstruct_indexes[x];
708 if ((delta*delta) > fuzz)
709 {
710 channel_distortion[IndexChannel]++;
711 count++;
712 }
713 }
714 if (count != 0)
715 channel_distortion[CompositeChannels]++;
716 p++;
717 q++;
718 }
719#if defined(MAGICKCORE_OPENMP_SUPPORT)
720 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
721#endif
722 for (i=0; i <= (ssize_t) CompositeChannels; i++)
723 distortion[i]+=channel_distortion[i];
724 }
725 reconstruct_view=DestroyCacheView(reconstruct_view);
726 image_view=DestroyCacheView(image_view);
727 area=MagickSafeReciprocal((double) columns*rows);
728 for (j=0; j <= CompositeChannels; j++)
729 distortion[j]*=area;
730 return(status);
731}
732
733
734
735static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
736 const Image *reconstruct_image,const ChannelType channel,
737 double *distortion,ExceptionInfo *exception)
738{
739 CacheView
740 *image_view,
741 *reconstruct_view;
742
743 MagickBooleanType
744 status;
745
746 size_t
747 columns,
748 rows;
749
750 ssize_t
751 i,
752 y;
753
754 status=MagickTrue;
755 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
756 image_view=AcquireVirtualCacheView(image,exception);
757 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
758#if defined(MAGICKCORE_OPENMP_SUPPORT)
759 #pragma omp parallel for schedule(static) shared(status) \
760 magick_number_threads(image,image,rows,1)
761#endif
762 for (y=0; y < (ssize_t) rows; y++)
763 {
764 double
765 channel_distortion[CompositeChannels+1];
766
767 const IndexPacket
768 *magick_restrict indexes,
769 *magick_restrict reconstruct_indexes;
770
771 const PixelPacket
772 *magick_restrict p,
773 *magick_restrict q;
774
775 ssize_t
776 i,
777 x;
778
779 if (status == MagickFalse)
780 continue;
781 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
782 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
783 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
784 {
785 status=MagickFalse;
786 continue;
787 }
788 indexes=GetCacheViewVirtualIndexQueue(image_view);
789 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
790 (void) memset(channel_distortion,0,sizeof(channel_distortion));
791 for (x=0; x < (ssize_t) columns; x++)
792 {
793 MagickRealType
794 distance,
795 Da,
796 Sa;
797
798 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
799 ((double) QuantumRange-(double) OpaqueOpacity));
800 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
801 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
802 OpaqueOpacity));
803 if ((channel & RedChannel) != 0)
804 {
805 distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
806 (double) GetPixelRed(q));
807 channel_distortion[RedChannel]+=distance;
808 channel_distortion[CompositeChannels]+=distance;
809 }
810 if ((channel & GreenChannel) != 0)
811 {
812 distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
813 (double) GetPixelGreen(q));
814 channel_distortion[GreenChannel]+=distance;
815 channel_distortion[CompositeChannels]+=distance;
816 }
817 if ((channel & BlueChannel) != 0)
818 {
819 distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
820 (double) GetPixelBlue(q));
821 channel_distortion[BlueChannel]+=distance;
822 channel_distortion[CompositeChannels]+=distance;
823 }
824 if (((channel & OpacityChannel) != 0) &&
825 (image->matte != MagickFalse))
826 {
827 distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
828 GetPixelOpacity(q));
829 channel_distortion[OpacityChannel]+=distance;
830 channel_distortion[CompositeChannels]+=distance;
831 }
832 if (((channel & IndexChannel) != 0) &&
833 (image->colorspace == CMYKColorspace))
834 {
835 distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
836 (double) GetPixelIndex(reconstruct_indexes+x));
837 channel_distortion[BlackChannel]+=distance;
838 channel_distortion[CompositeChannels]+=distance;
839 }
840 p++;
841 q++;
842 }
843#if defined(MAGICKCORE_OPENMP_SUPPORT)
844 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
845#endif
846 for (i=0; i <= (ssize_t) CompositeChannels; i++)
847 distortion[i]+=channel_distortion[i];
848 }
849 reconstruct_view=DestroyCacheView(reconstruct_view);
850 image_view=DestroyCacheView(image_view);
851 for (i=0; i <= (ssize_t) CompositeChannels; i++)
852 distortion[i]/=((double) columns*rows);
853 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
854 return(status);
855}
856
857static MagickBooleanType GetMeanErrorPerPixel(Image *image,
858 const Image *reconstruct_image,const ChannelType channel,double *distortion,
859 ExceptionInfo *exception)
860{
861 CacheView
862 *image_view,
863 *reconstruct_view;
864
865 double
866 maximum_error = -MagickMaximumValue,
867 mean_error = 0.0;
868
869 MagickBooleanType
870 status;
871
872 size_t
873 columns,
874 rows;
875
876 ssize_t
877 i,
878 y;
879
880 status=MagickTrue;
881 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
882 image_view=AcquireVirtualCacheView(image,exception);
883 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
884#if defined(MAGICKCORE_OPENMP_SUPPORT)
885 #pragma omp parallel for schedule(static) shared(maximum_error,status) \
886 magick_number_threads(image,image,rows,1)
887#endif
888 for (y=0; y < (ssize_t) rows; y++)
889 {
890 double
891 channel_distortion[CompositeChannels+1] = { 0.0 },
892 local_maximum = maximum_error,
893 local_mean_error = 0.0;
894
895 const IndexPacket
896 *magick_restrict indexes,
897 *magick_restrict reconstruct_indexes;
898
899 const PixelPacket
900 *magick_restrict p,
901 *magick_restrict q;
902
903 ssize_t
904 i,
905 x;
906
907 if (status == MagickFalse)
908 continue;
909 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
910 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
911 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
912 {
913 status=MagickFalse;
914 continue;
915 }
916 indexes=GetCacheViewVirtualIndexQueue(image_view);
917 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
918 (void) memset(channel_distortion,0,sizeof(channel_distortion));
919 for (x=0; x < (ssize_t) columns; x++)
920 {
921 MagickRealType
922 distance,
923 Da,
924 Sa;
925
926 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
927 ((double) QuantumRange-(double) OpaqueOpacity));
928 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
929 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
930 OpaqueOpacity));
931 if ((channel & RedChannel) != 0)
932 {
933 distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
934 (double) GetPixelRed(q));
935 channel_distortion[RedChannel]+=distance;
936 channel_distortion[CompositeChannels]+=distance;
937 local_mean_error+=distance*distance;
938 if (distance > local_maximum)
939 local_maximum=distance;
940 }
941 if ((channel & GreenChannel) != 0)
942 {
943 distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
944 (double) GetPixelGreen(q));
945 channel_distortion[GreenChannel]+=distance;
946 channel_distortion[CompositeChannels]+=distance;
947 local_mean_error+=distance*distance;
948 if (distance > local_maximum)
949 local_maximum=distance;
950 }
951 if ((channel & BlueChannel) != 0)
952 {
953 distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
954 (double) GetPixelBlue(q));
955 channel_distortion[BlueChannel]+=distance;
956 channel_distortion[CompositeChannels]+=distance;
957 local_mean_error+=distance*distance;
958 if (distance > local_maximum)
959 local_maximum=distance;
960 }
961 if (((channel & OpacityChannel) != 0) &&
962 (image->matte != MagickFalse))
963 {
964 distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
965 GetPixelOpacity(q));
966 channel_distortion[OpacityChannel]+=distance;
967 channel_distortion[CompositeChannels]+=distance;
968 local_mean_error+=distance*distance;
969 if (distance > local_maximum)
970 local_maximum=distance;
971 }
972 if (((channel & IndexChannel) != 0) &&
973 (image->colorspace == CMYKColorspace))
974 {
975 distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
976 (double) GetPixelIndex(reconstruct_indexes+x));
977 channel_distortion[BlackChannel]+=distance;
978 channel_distortion[CompositeChannels]+=distance;
979 local_mean_error+=distance*distance;
980 if (distance > local_maximum)
981 local_maximum=distance;
982 }
983 p++;
984 q++;
985 }
986#if defined(MAGICKCORE_OPENMP_SUPPORT)
987 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
988#endif
989 {
990 for (i=0; i <= (ssize_t) CompositeChannels; i++)
991 distortion[i]+=channel_distortion[i];
992 mean_error+=local_mean_error;
993 if (local_maximum > maximum_error)
994 maximum_error=local_maximum;
995 }
996 }
997 reconstruct_view=DestroyCacheView(reconstruct_view);
998 image_view=DestroyCacheView(image_view);
999 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1000 distortion[i]/=((double) columns*rows);
1001 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1002 image->error.mean_error_per_pixel=QuantumRange*distortion[CompositeChannels];
1003 image->error.normalized_mean_error=mean_error/((double) columns*rows);
1004 image->error.normalized_maximum_error=maximum_error;
1005 return(status);
1006}
1007
1008static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
1009 const Image *reconstruct_image,const ChannelType channel,
1010 double *distortion,ExceptionInfo *exception)
1011{
1012 CacheView
1013 *image_view,
1014 *reconstruct_view;
1015
1016 double
1017 area;
1018
1019 MagickBooleanType
1020 status;
1021
1022 size_t
1023 columns,
1024 rows;
1025
1026 ssize_t
1027 i,
1028 y;
1029
1030 status=MagickTrue;
1031 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
1032 image_view=AcquireVirtualCacheView(image,exception);
1033 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1034#if defined(MAGICKCORE_OPENMP_SUPPORT)
1035 #pragma omp parallel for schedule(static) shared(distortion,status) \
1036 magick_number_threads(image,image,rows,1)
1037#endif
1038 for (y=0; y < (ssize_t) rows; y++)
1039 {
1040 double
1041 channel_distortion[CompositeChannels+1] = { 0.0 };
1042
1043 const IndexPacket
1044 *magick_restrict indexes,
1045 *magick_restrict reconstruct_indexes;
1046
1047 const PixelPacket
1048 *magick_restrict p,
1049 *magick_restrict q;
1050
1051 ssize_t
1052 i,
1053 x;
1054
1055 if (status == MagickFalse)
1056 continue;
1057 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1058 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1059 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1060 {
1061 status=MagickFalse;
1062 continue;
1063 }
1064 indexes=GetCacheViewVirtualIndexQueue(image_view);
1065 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1066 for (x=0; x < (ssize_t) columns; x++)
1067 {
1068 double
1069 distance,
1070 Da,
1071 Sa;
1072
1073 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1074 ((double) QuantumRange-(double) OpaqueOpacity));
1075 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1076 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
1077 OpaqueOpacity));
1078 if ((channel & RedChannel) != 0)
1079 {
1080 distance=QuantumScale*(Sa*(double) GetPixelRed(p)-Da*(double)
1081 GetPixelRed(q));
1082 channel_distortion[RedChannel]+=distance*distance;
1083 channel_distortion[CompositeChannels]+=distance*distance;
1084 }
1085 if ((channel & GreenChannel) != 0)
1086 {
1087 distance=QuantumScale*(Sa*(double) GetPixelGreen(p)-Da*(double)
1088 GetPixelGreen(q));
1089 channel_distortion[GreenChannel]+=distance*distance;
1090 channel_distortion[CompositeChannels]+=distance*distance;
1091 }
1092 if ((channel & BlueChannel) != 0)
1093 {
1094 distance=QuantumScale*(Sa*(double) GetPixelBlue(p)-Da*(double)
1095 GetPixelBlue(q));
1096 channel_distortion[BlueChannel]+=distance*distance;
1097 channel_distortion[CompositeChannels]+=distance*distance;
1098 }
1099 if (((channel & OpacityChannel) != 0) &&
1100 (image->matte != MagickFalse))
1101 {
1102 distance=QuantumScale*((double) GetPixelOpacity(p)-(double)
1103 GetPixelOpacity(q));
1104 channel_distortion[OpacityChannel]+=distance*distance;
1105 channel_distortion[CompositeChannels]+=distance*distance;
1106 }
1107 if (((channel & IndexChannel) != 0) &&
1108 (image->colorspace == CMYKColorspace) &&
1109 (reconstruct_image->colorspace == CMYKColorspace))
1110 {
1111 distance=QuantumScale*(Sa*(double) GetPixelIndex(indexes+x)-Da*
1112 (double) GetPixelIndex(reconstruct_indexes+x));
1113 channel_distortion[BlackChannel]+=distance*distance;
1114 channel_distortion[CompositeChannels]+=distance*distance;
1115 }
1116 p++;
1117 q++;
1118 }
1119#if defined(MAGICKCORE_OPENMP_SUPPORT)
1120 #pragma omp critical (MagickCore_GetMeanSquaredError)
1121#endif
1122 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1123 distortion[i]+=channel_distortion[i];
1124 }
1125 reconstruct_view=DestroyCacheView(reconstruct_view);
1126 image_view=DestroyCacheView(image_view);
1127 area=MagickSafeReciprocal((double) columns*rows);
1128 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1129 distortion[i]*=area;
1130 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1131 return(status);
1132}
1133
1134static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
1135 const Image *image,const Image *reconstruct_image,const ChannelType channel,
1136 double *distortion,ExceptionInfo *exception)
1137{
1138#define SimilarityImageTag "Similarity/Image"
1139
1140 CacheView
1141 *image_view,
1142 *reconstruct_view;
1143
1144 ChannelStatistics
1145 *image_statistics,
1146 *reconstruct_statistics;
1147
1148 double
1149 alpha_variance[CompositeChannels+1] = { 0.0 },
1150 area,
1151 beta_variance[CompositeChannels+1] = { 0.0 };
1152
1153 MagickBooleanType
1154 status;
1155
1156 MagickOffsetType
1157 progress;
1158
1159 size_t
1160 columns,
1161 rows;
1162
1163 ssize_t
1164 i,
1165 y;
1166
1167 /*
1168 Normalize to account for variation due to lighting and exposure condition.
1169 */
1170 image_statistics=GetImageChannelStatistics(image,exception);
1171 reconstruct_statistics=GetImageChannelStatistics(reconstruct_image,exception);
1172 if ((image_statistics == (ChannelStatistics *) NULL) ||
1173 (reconstruct_statistics == (ChannelStatistics *) NULL))
1174 {
1175 if (image_statistics != (ChannelStatistics *) NULL)
1176 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1177 image_statistics);
1178 if (reconstruct_statistics != (ChannelStatistics *) NULL)
1179 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1180 reconstruct_statistics);
1181 return(MagickFalse);
1182 }
1183 (void) memset(distortion,0,(CompositeChannels+1)*sizeof(*distortion));
1184 status=MagickTrue;
1185 progress=0;
1186 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
1187 image_view=AcquireVirtualCacheView(image,exception);
1188 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1189#if defined(MAGICKCORE_OPENMP_SUPPORT)
1190 #pragma omp parallel for schedule(static) shared(status) \
1191 magick_number_threads(image,image,rows,1)
1192#endif
1193 for (y=0; y < (ssize_t) rows; y++)
1194 {
1195 const IndexPacket
1196 *magick_restrict indexes,
1197 *magick_restrict reconstruct_indexes;
1198
1199 const PixelPacket
1200 *magick_restrict p,
1201 *magick_restrict q;
1202
1203 double
1204 channel_alpha_variance[CompositeChannels+1] = { 0.0 },
1205 channel_beta_variance[CompositeChannels+1] = { 0.0 },
1206 channel_distortion[CompositeChannels+1] = { 0.0 };
1207
1208 ssize_t
1209 x;
1210
1211 if (status == MagickFalse)
1212 continue;
1213 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1214 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1215 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1216 {
1217 status=MagickFalse;
1218 continue;
1219 }
1220 indexes=GetCacheViewVirtualIndexQueue(image_view);
1221 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1222 for (x=0; x < (ssize_t) columns; x++)
1223 {
1224 MagickRealType
1225 alpha,
1226 beta,
1227 Da,
1228 Sa;
1229
1230 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1231 (double) QuantumRange);
1232 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1233 (double) GetPixelAlpha(q) : (double) QuantumRange);
1234 if ((channel & RedChannel) != 0)
1235 {
1236 alpha=QuantumScale*(Sa*(double) GetPixelRed(p)-
1237 image_statistics[RedChannel].mean);
1238 beta=QuantumScale*(Da*(double) GetPixelRed(q)-
1239 reconstruct_statistics[RedChannel].mean);
1240 channel_distortion[RedChannel]+=alpha*beta;
1241 channel_alpha_variance[RedChannel]+=alpha*alpha;
1242 channel_beta_variance[RedChannel]+=beta*beta;
1243 }
1244 if ((channel & GreenChannel) != 0)
1245 {
1246 alpha=QuantumScale*(Sa*(double) GetPixelGreen(p)-
1247 image_statistics[GreenChannel].mean);
1248 beta=QuantumScale*(Da*(double) GetPixelGreen(q)-
1249 reconstruct_statistics[GreenChannel].mean);
1250 channel_distortion[GreenChannel]+=alpha*beta;
1251 channel_alpha_variance[GreenChannel]+=alpha*alpha;
1252 channel_beta_variance[GreenChannel]+=beta*beta;
1253 }
1254 if ((channel & BlueChannel) != 0)
1255 {
1256 alpha=QuantumScale*(Sa*(double) GetPixelBlue(p)-
1257 image_statistics[BlueChannel].mean);
1258 beta=QuantumScale*(Da*(double) GetPixelBlue(q)-
1259 reconstruct_statistics[BlueChannel].mean);
1260 channel_distortion[BlueChannel]+=alpha*beta;
1261 channel_alpha_variance[BlueChannel]+=alpha*alpha;
1262 channel_beta_variance[BlueChannel]+=beta*beta;
1263 }
1264 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1265 {
1266 alpha=QuantumScale*((double) GetPixelAlpha(p)-
1267 image_statistics[AlphaChannel].mean);
1268 beta=QuantumScale*((double) GetPixelAlpha(q)-
1269 reconstruct_statistics[AlphaChannel].mean);
1270 channel_distortion[OpacityChannel]+=alpha*beta;
1271 channel_alpha_variance[OpacityChannel]+=alpha*alpha;
1272 channel_beta_variance[OpacityChannel]+=beta*beta;
1273 }
1274 if (((channel & IndexChannel) != 0) &&
1275 (image->colorspace == CMYKColorspace) &&
1276 (reconstruct_image->colorspace == CMYKColorspace))
1277 {
1278 alpha=QuantumScale*(Sa*(double) GetPixelIndex(indexes+x)-
1279 image_statistics[BlackChannel].mean);
1280 beta=QuantumScale*(Da*(double) GetPixelIndex(reconstruct_indexes+
1281 x)-reconstruct_statistics[BlackChannel].mean);
1282 channel_distortion[BlackChannel]+=alpha*beta;
1283 channel_alpha_variance[BlackChannel]+=alpha*alpha;
1284 channel_beta_variance[BlackChannel]+=beta*beta;
1285 }
1286 p++;
1287 q++;
1288 }
1289#if defined(MAGICKCORE_OPENMP_SUPPORT)
1290 #pragma omp critical (GetNormalizedCrossCorrelationDistortion)
1291#endif
1292 {
1293 ssize_t
1294 j;
1295
1296 for (j=0; j < (ssize_t) CompositeChannels; j++)
1297 {
1298 distortion[j]+=channel_distortion[j];
1299 alpha_variance[j]+=channel_alpha_variance[j];
1300 beta_variance[j]+=channel_beta_variance[j];
1301 }
1302 }
1303 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1304 {
1305 MagickBooleanType
1306 proceed;
1307
1308#if defined(MAGICKCORE_OPENMP_SUPPORT)
1309 #pragma omp atomic
1310#endif
1311 progress++;
1312 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1313 if (proceed == MagickFalse)
1314 status=MagickFalse;
1315 }
1316 }
1317 reconstruct_view=DestroyCacheView(reconstruct_view);
1318 image_view=DestroyCacheView(image_view);
1319 /*
1320 Divide by the standard deviation.
1321 */
1322 area=MagickSafeReciprocal((double) columns*rows);
1323 for (i=0; i < (ssize_t) CompositeChannels; i++)
1324 {
1325 distortion[i]*=area;
1326 alpha_variance[i]*=area;
1327 beta_variance[i]*=area;
1328 distortion[i]*=MagickSafeReciprocal(sqrt(alpha_variance[i]*
1329 beta_variance[i]));
1330 distortion[CompositeChannels]+=distortion[i];
1331 }
1332 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1333 /*
1334 Free resources.
1335 */
1336 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1337 reconstruct_statistics);
1338 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1339 image_statistics);
1340 return(status);
1341}
1342
1343static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1344 const Image *reconstruct_image,const ChannelType channel,
1345 double *distortion,ExceptionInfo *exception)
1346{
1347 CacheView
1348 *image_view,
1349 *reconstruct_view;
1350
1351 MagickBooleanType
1352 status;
1353
1354 size_t
1355 columns,
1356 rows;
1357
1358 ssize_t
1359 y;
1360
1361 status=MagickTrue;
1362 (void) memset(distortion,0,(CompositeChannels+1)*sizeof(*distortion));
1363 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
1364 image_view=AcquireVirtualCacheView(image,exception);
1365 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1366#if defined(MAGICKCORE_OPENMP_SUPPORT)
1367 #pragma omp parallel for schedule(static) shared(status) \
1368 magick_number_threads(image,image,rows,1)
1369#endif
1370 for (y=0; y < (ssize_t) rows; y++)
1371 {
1372 double
1373 channel_distortion[CompositeChannels+1];
1374
1375 const IndexPacket
1376 *magick_restrict indexes,
1377 *magick_restrict reconstruct_indexes;
1378
1379 const PixelPacket
1380 *magick_restrict p,
1381 *magick_restrict q;
1382
1383 ssize_t
1384 i,
1385 x;
1386
1387 if (status == MagickFalse)
1388 continue;
1389 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1390 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1391 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1392 {
1393 status=MagickFalse;
1394 continue;
1395 }
1396 indexes=GetCacheViewVirtualIndexQueue(image_view);
1397 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1398 (void) memset(channel_distortion,0,(CompositeChannels+1)*
1399 sizeof(*channel_distortion));
1400 for (x=0; x < (ssize_t) columns; x++)
1401 {
1402 MagickRealType
1403 distance,
1404 Da,
1405 Sa;
1406
1407 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1408 ((double) QuantumRange-(double) OpaqueOpacity));
1409 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1410 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
1411 OpaqueOpacity));
1412 if ((channel & RedChannel) != 0)
1413 {
1414 distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
1415 (double) GetPixelRed(q));
1416 if (distance > channel_distortion[RedChannel])
1417 channel_distortion[RedChannel]=distance;
1418 if (distance > channel_distortion[CompositeChannels])
1419 channel_distortion[CompositeChannels]=distance;
1420 }
1421 if ((channel & GreenChannel) != 0)
1422 {
1423 distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
1424 (double) GetPixelGreen(q));
1425 if (distance > channel_distortion[GreenChannel])
1426 channel_distortion[GreenChannel]=distance;
1427 if (distance > channel_distortion[CompositeChannels])
1428 channel_distortion[CompositeChannels]=distance;
1429 }
1430 if ((channel & BlueChannel) != 0)
1431 {
1432 distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
1433 (double) GetPixelBlue(q));
1434 if (distance > channel_distortion[BlueChannel])
1435 channel_distortion[BlueChannel]=distance;
1436 if (distance > channel_distortion[CompositeChannels])
1437 channel_distortion[CompositeChannels]=distance;
1438 }
1439 if (((channel & OpacityChannel) != 0) &&
1440 (image->matte != MagickFalse))
1441 {
1442 distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
1443 GetPixelOpacity(q));
1444 if (distance > channel_distortion[OpacityChannel])
1445 channel_distortion[OpacityChannel]=distance;
1446 if (distance > channel_distortion[CompositeChannels])
1447 channel_distortion[CompositeChannels]=distance;
1448 }
1449 if (((channel & IndexChannel) != 0) &&
1450 (image->colorspace == CMYKColorspace) &&
1451 (reconstruct_image->colorspace == CMYKColorspace))
1452 {
1453 distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
1454 (double) GetPixelIndex(reconstruct_indexes+x));
1455 if (distance > channel_distortion[BlackChannel])
1456 channel_distortion[BlackChannel]=distance;
1457 if (distance > channel_distortion[CompositeChannels])
1458 channel_distortion[CompositeChannels]=distance;
1459 }
1460 p++;
1461 q++;
1462 }
1463#if defined(MAGICKCORE_OPENMP_SUPPORT)
1464 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1465#endif
1466 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1467 if (channel_distortion[i] > distortion[i])
1468 distortion[i]=channel_distortion[i];
1469 }
1470 reconstruct_view=DestroyCacheView(reconstruct_view);
1471 image_view=DestroyCacheView(image_view);
1472 return(status);
1473}
1474
1475static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1476 const Image *reconstruct_image,const ChannelType channel,
1477 double *distortion,ExceptionInfo *exception)
1478{
1479 MagickBooleanType
1480 status;
1481
1482 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1483 exception);
1484 if ((channel & RedChannel) != 0)
1485 distortion[RedChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1486 distortion[RedChannel]))/MagickPSNRDistortion;
1487 if ((channel & GreenChannel) != 0)
1488 distortion[GreenChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1489 distortion[GreenChannel]))/MagickPSNRDistortion;
1490 if ((channel & BlueChannel) != 0)
1491 distortion[BlueChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1492 distortion[BlueChannel]))/MagickPSNRDistortion;
1493 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1494 distortion[OpacityChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1495 distortion[OpacityChannel]))/MagickPSNRDistortion;
1496 if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1497 distortion[BlackChannel]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1498 distortion[BlackChannel]))/MagickPSNRDistortion;
1499 distortion[CompositeChannels]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1500 distortion[CompositeChannels]))/MagickPSNRDistortion;
1501 return(status);
1502}
1503
1504static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1505 const Image *reconstruct_image,const ChannelType channel,double *distortion,
1506 ExceptionInfo *exception)
1507{
1508#define PHASHNormalizationFactor 389.373723242
1509
1510 ChannelPerceptualHash
1511 *image_phash,
1512 *reconstruct_phash;
1513
1514 double
1515 difference;
1516
1517 ssize_t
1518 i;
1519
1520 /*
1521 Compute perceptual hash in the sRGB colorspace.
1522 */
1523 image_phash=GetImageChannelPerceptualHash(image,exception);
1524 if (image_phash == (ChannelPerceptualHash *) NULL)
1525 return(MagickFalse);
1526 reconstruct_phash=GetImageChannelPerceptualHash(reconstruct_image,exception);
1527 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1528 {
1529 image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1530 return(MagickFalse);
1531 }
1532 for (i=0; i < MaximumNumberOfImageMoments; i++)
1533 {
1534 /*
1535 Compute sum of moment differences squared.
1536 */
1537 if ((channel & RedChannel) != 0)
1538 {
1539 difference=reconstruct_phash[RedChannel].P[i]-
1540 image_phash[RedChannel].P[i];
1541 if (IsNaN(difference) != 0)
1542 difference=0.0;
1543 distortion[RedChannel]+=difference*difference;
1544 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1545 }
1546 if ((channel & GreenChannel) != 0)
1547 {
1548 difference=reconstruct_phash[GreenChannel].P[i]-
1549 image_phash[GreenChannel].P[i];
1550 if (IsNaN(difference) != 0)
1551 difference=0.0;
1552 distortion[GreenChannel]+=difference*difference;
1553 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1554 }
1555 if ((channel & BlueChannel) != 0)
1556 {
1557 difference=reconstruct_phash[BlueChannel].P[i]-
1558 image_phash[BlueChannel].P[i];
1559 if (IsNaN(difference) != 0)
1560 difference=0.0;
1561 distortion[BlueChannel]+=difference*difference;
1562 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1563 }
1564 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1565 (reconstruct_image->matte != MagickFalse))
1566 {
1567 difference=reconstruct_phash[OpacityChannel].P[i]-
1568 image_phash[OpacityChannel].P[i];
1569 if (IsNaN(difference) != 0)
1570 difference=0.0;
1571 distortion[OpacityChannel]+=difference*difference;
1572 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1573 }
1574 if (((channel & IndexChannel) != 0) &&
1575 (image->colorspace == CMYKColorspace) &&
1576 (reconstruct_image->colorspace == CMYKColorspace))
1577 {
1578 difference=reconstruct_phash[IndexChannel].P[i]-
1579 image_phash[IndexChannel].P[i];
1580 if (IsNaN(difference) != 0)
1581 difference=0.0;
1582 distortion[IndexChannel]+=difference*difference;
1583 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1584 }
1585 }
1586 /*
1587 Compute perceptual hash in the HCLP colorspace.
1588 */
1589 for (i=0; i < MaximumNumberOfImageMoments; i++)
1590 {
1591 /*
1592 Compute sum of moment differences squared.
1593 */
1594 if ((channel & RedChannel) != 0)
1595 {
1596 difference=reconstruct_phash[RedChannel].Q[i]-
1597 image_phash[RedChannel].Q[i];
1598 if (IsNaN(difference) != 0)
1599 difference=0.0;
1600 distortion[RedChannel]+=difference*difference;
1601 distortion[CompositeChannels]+=difference*difference;
1602 }
1603 if ((channel & GreenChannel) != 0)
1604 {
1605 difference=reconstruct_phash[GreenChannel].Q[i]-
1606 image_phash[GreenChannel].Q[i];
1607 if (IsNaN(difference) != 0)
1608 difference=0.0;
1609 distortion[GreenChannel]+=difference*difference;
1610 distortion[CompositeChannels]+=difference*difference;
1611 }
1612 if ((channel & BlueChannel) != 0)
1613 {
1614 difference=reconstruct_phash[BlueChannel].Q[i]-
1615 image_phash[BlueChannel].Q[i];
1616 distortion[BlueChannel]+=difference*difference;
1617 distortion[CompositeChannels]+=difference*difference;
1618 }
1619 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1620 (reconstruct_image->matte != MagickFalse))
1621 {
1622 difference=reconstruct_phash[OpacityChannel].Q[i]-
1623 image_phash[OpacityChannel].Q[i];
1624 if (IsNaN(difference) != 0)
1625 difference=0.0;
1626 distortion[OpacityChannel]+=difference*difference;
1627 distortion[CompositeChannels]+=difference*difference;
1628 }
1629 if (((channel & IndexChannel) != 0) &&
1630 (image->colorspace == CMYKColorspace) &&
1631 (reconstruct_image->colorspace == CMYKColorspace))
1632 {
1633 difference=reconstruct_phash[IndexChannel].Q[i]-
1634 image_phash[IndexChannel].Q[i];
1635 if (IsNaN(difference) != 0)
1636 difference=0.0;
1637 distortion[IndexChannel]+=difference*difference;
1638 distortion[CompositeChannels]+=difference*difference;
1639 }
1640 }
1641 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1642 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1643 distortion[i]/=PHASHNormalizationFactor;
1644 /*
1645 Free resources.
1646 */
1647 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1648 reconstruct_phash);
1649 image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1650 return(MagickTrue);
1651}
1652
1653static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1654 const Image *reconstruct_image,const ChannelType channel,double *distortion,
1655 ExceptionInfo *exception)
1656{
1657#define RMSESquareRoot(x) sqrt((x) < 0.0 ? 0.0 : (x))
1658
1659 MagickBooleanType
1660 status;
1661
1662 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1663 exception);
1664 if ((channel & RedChannel) != 0)
1665 distortion[RedChannel]=RMSESquareRoot(distortion[RedChannel]);
1666 if ((channel & GreenChannel) != 0)
1667 distortion[GreenChannel]=RMSESquareRoot(distortion[GreenChannel]);
1668 if ((channel & BlueChannel) != 0)
1669 distortion[BlueChannel]=RMSESquareRoot(distortion[BlueChannel]);
1670 if (((channel & OpacityChannel) != 0) &&
1671 (image->matte != MagickFalse))
1672 distortion[OpacityChannel]=RMSESquareRoot(distortion[OpacityChannel]);
1673 if (((channel & IndexChannel) != 0) &&
1674 (image->colorspace == CMYKColorspace))
1675 distortion[BlackChannel]=RMSESquareRoot(distortion[BlackChannel]);
1676 distortion[CompositeChannels]=RMSESquareRoot(distortion[CompositeChannels]);
1677 return(status);
1678}
1679
1680MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1681 const Image *reconstruct_image,const ChannelType channel,
1682 const MetricType metric,double *distortion,ExceptionInfo *exception)
1683{
1684 double
1685 *channel_distortion;
1686
1687 MagickBooleanType
1688 status;
1689
1690 size_t
1691 length;
1692
1693 ssize_t
1694 i;
1695
1696 assert(image != (Image *) NULL);
1697 assert(image->signature == MagickCoreSignature);
1698 assert(reconstruct_image != (const Image *) NULL);
1699 assert(reconstruct_image->signature == MagickCoreSignature);
1700 assert(distortion != (double *) NULL);
1701 if (IsEventLogging() != MagickFalse)
1702 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1703 *distortion=0.0;
1704 if (metric != PerceptualHashErrorMetric)
1705 if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1706 ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1707 /*
1708 Get image distortion.
1709 */
1710 length=CompositeChannels+1UL;
1711 channel_distortion=(double *) AcquireQuantumMemory(length,
1712 sizeof(*channel_distortion));
1713 if (channel_distortion == (double *) NULL)
1714 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1715 (void) memset(channel_distortion,0,length*sizeof(*channel_distortion));
1716 switch (metric)
1717 {
1718 case AbsoluteErrorMetric:
1719 {
1720 status=GetAbsoluteDistortion(image,reconstruct_image,channel,
1721 channel_distortion,exception);
1722 break;
1723 }
1724 case FuzzErrorMetric:
1725 {
1726 status=GetFuzzDistortion(image,reconstruct_image,channel,
1727 channel_distortion,exception);
1728 break;
1729 }
1730 case MeanAbsoluteErrorMetric:
1731 {
1732 status=GetMeanAbsoluteDistortion(image,reconstruct_image,channel,
1733 channel_distortion,exception);
1734 break;
1735 }
1736 case MeanErrorPerPixelMetric:
1737 {
1738 status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1739 channel_distortion,exception);
1740 break;
1741 }
1742 case MeanSquaredErrorMetric:
1743 {
1744 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,
1745 channel_distortion,exception);
1746 break;
1747 }
1748 case NormalizedCrossCorrelationErrorMetric:
1749 {
1750 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1751 channel,channel_distortion,exception);
1752 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1753 channel_distortion[i]=1.0-channel_distortion[i];
1754 break;
1755 }
1756 case PeakAbsoluteErrorMetric:
1757 {
1758 status=GetPeakAbsoluteDistortion(image,reconstruct_image,channel,
1759 channel_distortion,exception);
1760 break;
1761 }
1762 case PeakSignalToNoiseRatioMetric:
1763 {
1764 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1765 channel_distortion,exception);
1766 break;
1767 }
1768 case PerceptualHashErrorMetric:
1769 {
1770 status=GetPerceptualHashDistortion(image,reconstruct_image,channel,
1771 channel_distortion,exception);
1772 break;
1773 }
1774 case RootMeanSquaredErrorMetric:
1775 case UndefinedErrorMetric:
1776 default:
1777 {
1778 status=GetRootMeanSquaredDistortion(image,reconstruct_image,channel,
1779 channel_distortion,exception);
1780 break;
1781 }
1782 }
1783 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1784 channel_distortion[i]=MagickMin(MagickMax(channel_distortion[i],0.0),1.0);
1785 *distortion=channel_distortion[CompositeChannels];
1786 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1787 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1788 *distortion);
1789 return(status);
1790}
1791
1792/*
1793%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1794% %
1795% %
1796% %
1797% G e t I m a g e C h a n n e l D i s t o r t i o n s %
1798% %
1799% %
1800% %
1801%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1802%
1803% GetImageChannelDistortions() compares the image channels of an image to a
1804% reconstructed image and returns the specified distortion metric for each
1805% channel.
1806%
1807% The format of the GetImageChannelDistortions method is:
1808%
1809% double *GetImageChannelDistortions(const Image *image,
1810% const Image *reconstruct_image,const MetricType metric,
1811% ExceptionInfo *exception)
1812%
1813% A description of each parameter follows:
1814%
1815% o image: the image.
1816%
1817% o reconstruct_image: the reconstruct image.
1818%
1819% o metric: the metric.
1820%
1821% o exception: return any errors or warnings in this structure.
1822%
1823*/
1824MagickExport double *GetImageChannelDistortions(Image *image,
1825 const Image *reconstruct_image,const MetricType metric,
1826 ExceptionInfo *exception)
1827{
1828 double
1829 *channel_distortion;
1830
1831 MagickBooleanType
1832 status;
1833
1834 size_t
1835 length;
1836
1837 ssize_t
1838 i;
1839
1840 assert(image != (Image *) NULL);
1841 assert(image->signature == MagickCoreSignature);
1842 assert(reconstruct_image != (const Image *) NULL);
1843 assert(reconstruct_image->signature == MagickCoreSignature);
1844 if (IsEventLogging() != MagickFalse)
1845 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1846 if (metric != PerceptualHashErrorMetric)
1847 if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1848 {
1849 (void) ThrowMagickException(&image->exception,GetMagickModule(),
1850 ImageError,"ImageMorphologyDiffers","`%s'",image->filename);
1851 return((double *) NULL);
1852 }
1853 /*
1854 Get image distortion.
1855 */
1856 length=CompositeChannels+1UL;
1857 channel_distortion=(double *) AcquireQuantumMemory(length,
1858 sizeof(*channel_distortion));
1859 if (channel_distortion == (double *) NULL)
1860 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1861 (void) memset(channel_distortion,0,length*
1862 sizeof(*channel_distortion));
1863 status=MagickTrue;
1864 switch (metric)
1865 {
1866 case AbsoluteErrorMetric:
1867 {
1868 status=GetAbsoluteDistortion(image,reconstruct_image,CompositeChannels,
1869 channel_distortion,exception);
1870 break;
1871 }
1872 case FuzzErrorMetric:
1873 {
1874 status=GetFuzzDistortion(image,reconstruct_image,CompositeChannels,
1875 channel_distortion,exception);
1876 break;
1877 }
1878 case MeanAbsoluteErrorMetric:
1879 {
1880 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1881 CompositeChannels,channel_distortion,exception);
1882 break;
1883 }
1884 case MeanErrorPerPixelMetric:
1885 {
1886 status=GetMeanErrorPerPixel(image,reconstruct_image,CompositeChannels,
1887 channel_distortion,exception);
1888 break;
1889 }
1890 case MeanSquaredErrorMetric:
1891 {
1892 status=GetMeanSquaredDistortion(image,reconstruct_image,CompositeChannels,
1893 channel_distortion,exception);
1894 break;
1895 }
1896 case NormalizedCrossCorrelationErrorMetric:
1897 {
1898 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1899 CompositeChannels,channel_distortion,exception);
1900 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1901 channel_distortion[i]=1.0-channel_distortion[i];
1902 break;
1903 }
1904 case PeakAbsoluteErrorMetric:
1905 {
1906 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1907 CompositeChannels,channel_distortion,exception);
1908 break;
1909 }
1910 case PeakSignalToNoiseRatioMetric:
1911 {
1912 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1913 CompositeChannels,channel_distortion,exception);
1914 break;
1915 }
1916 case PerceptualHashErrorMetric:
1917 {
1918 status=GetPerceptualHashDistortion(image,reconstruct_image,
1919 CompositeChannels,channel_distortion,exception);
1920 break;
1921 }
1922 case RootMeanSquaredErrorMetric:
1923 case UndefinedErrorMetric:
1924 default:
1925 {
1926 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1927 CompositeChannels,channel_distortion,exception);
1928 break;
1929 }
1930 }
1931 if (status == MagickFalse)
1932 {
1933 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1934 return((double *) NULL);
1935 }
1936 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1937 channel_distortion[i]=MagickMin(MagickMax(channel_distortion[i],0.0),1.0);
1938 return(channel_distortion);
1939}
1940
1941/*
1942%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1943% %
1944% %
1945% %
1946% I s I m a g e s E q u a l %
1947% %
1948% %
1949% %
1950%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1951%
1952% IsImagesEqual() measures the difference between colors at each pixel
1953% location of two images. A value other than 0 means the colors match
1954% exactly. Otherwise an error measure is computed by summing over all
1955% pixels in an image the distance squared in RGB space between each image
1956% pixel and its corresponding pixel in the reconstruct image. The error
1957% measure is assigned to these image members:
1958%
1959% o mean_error_per_pixel: The mean error for any single pixel in
1960% the image.
1961%
1962% o normalized_mean_error: The normalized mean quantization error for
1963% any single pixel in the image. This distance measure is normalized to
1964% a range between 0 and 1. It is independent of the range of red, green,
1965% and blue values in the image.
1966%
1967% o normalized_maximum_error: The normalized maximum quantization
1968% error for any single pixel in the image. This distance measure is
1969% normalized to a range between 0 and 1. It is independent of the range
1970% of red, green, and blue values in your image.
1971%
1972% A small normalized mean square error, accessed as
1973% image->normalized_mean_error, suggests the images are very similar in
1974% spatial layout and color.
1975%
1976% The format of the IsImagesEqual method is:
1977%
1978% MagickBooleanType IsImagesEqual(Image *image,
1979% const Image *reconstruct_image)
1980%
1981% A description of each parameter follows.
1982%
1983% o image: the image.
1984%
1985% o reconstruct_image: the reconstruct image.
1986%
1987*/
1988MagickExport MagickBooleanType IsImagesEqual(Image *image,
1989 const Image *reconstruct_image)
1990{
1991 CacheView
1992 *image_view,
1993 *reconstruct_view;
1994
1995 ExceptionInfo
1996 *exception;
1997
1998 MagickBooleanType
1999 status;
2000
2001 MagickRealType
2002 area,
2003 gamma,
2004 maximum_error,
2005 mean_error,
2006 mean_error_per_pixel;
2007
2008 size_t
2009 columns,
2010 rows;
2011
2012 ssize_t
2013 y;
2014
2015 assert(image != (Image *) NULL);
2016 assert(image->signature == MagickCoreSignature);
2017 assert(reconstruct_image != (const Image *) NULL);
2018 assert(reconstruct_image->signature == MagickCoreSignature);
2019 exception=(&image->exception);
2020 if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
2021 ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
2022 area=0.0;
2023 maximum_error=0.0;
2024 mean_error_per_pixel=0.0;
2025 mean_error=0.0;
2026 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
2027 image_view=AcquireVirtualCacheView(image,exception);
2028 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2029 for (y=0; y < (ssize_t) rows; y++)
2030 {
2031 const IndexPacket
2032 *magick_restrict indexes,
2033 *magick_restrict reconstruct_indexes;
2034
2035 const PixelPacket
2036 *magick_restrict p,
2037 *magick_restrict q;
2038
2039 ssize_t
2040 x;
2041
2042 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2043 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2044 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
2045 break;
2046 indexes=GetCacheViewVirtualIndexQueue(image_view);
2047 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
2048 for (x=0; x < (ssize_t) columns; x++)
2049 {
2050 MagickRealType
2051 distance;
2052
2053 distance=fabs((double) GetPixelRed(p)-(double) GetPixelRed(q));
2054 mean_error_per_pixel+=distance;
2055 mean_error+=distance*distance;
2056 if (distance > maximum_error)
2057 maximum_error=distance;
2058 area++;
2059 distance=fabs((double) GetPixelGreen(p)-(double) GetPixelGreen(q));
2060 mean_error_per_pixel+=distance;
2061 mean_error+=distance*distance;
2062 if (distance > maximum_error)
2063 maximum_error=distance;
2064 area++;
2065 distance=fabs((double) GetPixelBlue(p)-(double) GetPixelBlue(q));
2066 mean_error_per_pixel+=distance;
2067 mean_error+=distance*distance;
2068 if (distance > maximum_error)
2069 maximum_error=distance;
2070 area++;
2071 if (image->matte != MagickFalse)
2072 {
2073 distance=fabs((double) GetPixelOpacity(p)-(double)
2074 GetPixelOpacity(q));
2075 mean_error_per_pixel+=distance;
2076 mean_error+=distance*distance;
2077 if (distance > maximum_error)
2078 maximum_error=distance;
2079 area++;
2080 }
2081 if ((image->colorspace == CMYKColorspace) &&
2082 (reconstruct_image->colorspace == CMYKColorspace))
2083 {
2084 distance=fabs((double) GetPixelIndex(indexes+x)-(double)
2085 GetPixelIndex(reconstruct_indexes+x));
2086 mean_error_per_pixel+=distance;
2087 mean_error+=distance*distance;
2088 if (distance > maximum_error)
2089 maximum_error=distance;
2090 area++;
2091 }
2092 p++;
2093 q++;
2094 }
2095 }
2096 reconstruct_view=DestroyCacheView(reconstruct_view);
2097 image_view=DestroyCacheView(image_view);
2098 gamma=MagickSafeReciprocal(area);
2099 image->error.mean_error_per_pixel=gamma*mean_error_per_pixel;
2100 image->error.normalized_mean_error=gamma*QuantumScale*QuantumScale*mean_error;
2101 image->error.normalized_maximum_error=QuantumScale*maximum_error;
2102 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2103 return(status);
2104}
2105
2106/*
2107%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2108% %
2109% %
2110% %
2111% S i m i l a r i t y I m a g e %
2112% %
2113% %
2114% %
2115%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2116%
2117% SimilarityImage() compares the reference image of the image and returns the
2118% best match offset. In addition, it returns a similarity image such that an
2119% exact match location is completely white and if none of the pixels match,
2120% black, otherwise some gray level in-between.
2121%
2122% The format of the SimilarityImageImage method is:
2123%
2124% Image *SimilarityImage(const Image *image,const Image *reference,
2125% RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2126%
2127% A description of each parameter follows:
2128%
2129% o image: the image.
2130%
2131% o reference: find an area of the image that closely resembles this image.
2132%
2133% o the best match offset of the reference image within the image.
2134%
2135% o similarity: the computed similarity between the images.
2136%
2137% o exception: return any errors or warnings in this structure.
2138%
2139*/
2140
2141static double GetSimilarityMetric(const Image *image,const Image *reference,
2142 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2143 ExceptionInfo *exception)
2144{
2145 double
2146 distortion;
2147
2148 Image
2149 *similarity_image;
2150
2151 MagickBooleanType
2152 status;
2153
2154 RectangleInfo
2155 geometry;
2156
2157 SetGeometry(reference,&geometry);
2158 geometry.x=x_offset;
2159 geometry.y=y_offset;
2160 similarity_image=CropImage(image,&geometry,exception);
2161 if (similarity_image == (Image *) NULL)
2162 return(0.0);
2163 distortion=0.0;
2164 status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2165 exception);
2166 (void) status;
2167 similarity_image=DestroyImage(similarity_image);
2168 return(distortion);
2169}
2170
2171MagickExport Image *SimilarityImage(Image *image,const Image *reference,
2172 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2173{
2174 Image
2175 *similarity_image;
2176
2177 similarity_image=SimilarityMetricImage(image,reference,
2178 RootMeanSquaredErrorMetric,offset,similarity_metric,exception);
2179 return(similarity_image);
2180}
2181
2182MagickExport Image *SimilarityMetricImage(Image *image,const Image *reconstruct,
2183 const MetricType metric,RectangleInfo *offset,double *similarity_metric,
2184 ExceptionInfo *exception)
2185{
2186#define SimilarityImageTag "Similarity/Image"
2187
2188 typedef struct
2189 {
2190 double
2191 similarity;
2192
2193 ssize_t
2194 x,
2195 y;
2196 } SimilarityInfo;
2197
2198 CacheView
2199 *similarity_view;
2200
2201 const char
2202 *artifact;
2203
2204 double
2205 similarity_threshold;
2206
2207 Image
2208 *similarity_image = (Image *) NULL;
2209
2210 MagickBooleanType
2211 status;
2212
2213 MagickOffsetType
2214 progress;
2215
2216 SimilarityInfo
2217 similarity_info = { 0 };
2218
2219 ssize_t
2220 y;
2221
2222 assert(image != (const Image *) NULL);
2223 assert(image->signature == MagickCoreSignature);
2224 assert(exception != (ExceptionInfo *) NULL);
2225 assert(exception->signature == MagickCoreSignature);
2226 assert(offset != (RectangleInfo *) NULL);
2227 if (IsEventLogging() != MagickFalse)
2228 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2229 SetGeometry(reconstruct,offset);
2230 *similarity_metric=0.0;
2231 offset->x=0;
2232 offset->y=0;
2233 if (ValidateImageMorphology(image,reconstruct) == MagickFalse)
2234 ThrowImageException(ImageError,"ImageMorphologyDiffers");
2235 if ((image->columns < reconstruct->columns) ||
2236 (image->rows < reconstruct->rows))
2237 {
2238 (void) ThrowMagickException(&image->exception,GetMagickModule(),
2239 OptionWarning,"GeometryDoesNotContainImage","`%s'",image->filename);
2240 return((Image *) NULL);
2241 }
2242 similarity_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2243 exception);
2244 if (similarity_image == (Image *) NULL)
2245 return((Image *) NULL);
2246 similarity_image->depth=32;
2247 similarity_image->colorspace=GRAYColorspace;
2248 similarity_image->matte=MagickFalse;
2249 status=SetImageStorageClass(similarity_image,DirectClass);
2250 if (status == MagickFalse)
2251 {
2252 InheritException(exception,&similarity_image->exception);
2253 return(DestroyImage(similarity_image));
2254 }
2255 /*
2256 Measure similarity of reconstruction image against image.
2257 */
2258 similarity_threshold=DefaultSimilarityThreshold;
2259 artifact=GetImageArtifact(image,"compare:similarity-threshold");
2260 if (artifact != (const char *) NULL)
2261 similarity_threshold=StringToDouble(artifact,(char **) NULL);
2262 status=MagickTrue;
2263 similarity_info.similarity=GetSimilarityMetric(image,reconstruct,metric,
2264 similarity_info.x,similarity_info.y,exception);
2265 progress=0;
2266 similarity_view=AcquireVirtualCacheView(similarity_image,exception);
2267#if defined(MAGICKCORE_OPENMP_SUPPORT)
2268 #pragma omp parallel for schedule(static) shared(status,similarity_info) \
2269 magick_number_threads(image,reconstruct,similarity_image->rows << 2,1)
2270#endif
2271 for (y=0; y < (ssize_t) similarity_image->rows; y++)
2272 {
2273 double
2274 similarity;
2275
2276 MagickBooleanType
2277 threshold_trigger = MagickFalse;
2278
2279 PixelPacket
2280 *magick_restrict q;
2281
2282 SimilarityInfo
2283 channel_info = similarity_info;
2284
2285 ssize_t
2286 x;
2287
2288 if (status == MagickFalse)
2289 continue;
2290 if (threshold_trigger != MagickFalse)
2291 continue;
2292 q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
2293 similarity_image->columns,1,exception);
2294 if (q == (PixelPacket *) NULL)
2295 {
2296 status=MagickFalse;
2297 continue;
2298 }
2299 for (x=0; x < (ssize_t) similarity_image->columns; x++)
2300 {
2301 MagickBooleanType
2302 update = MagickFalse;
2303
2304 similarity=GetSimilarityMetric(image,reconstruct,metric,x,y,exception);
2305 switch (metric)
2306 {
2307 case PeakSignalToNoiseRatioMetric:
2308 {
2309 if (similarity > channel_info.similarity)
2310 update=MagickTrue;
2311 break;
2312 }
2313 default:
2314 {
2315 if (similarity < channel_info.similarity)
2316 update=MagickTrue;
2317 break;
2318 }
2319 }
2320 if (update != MagickFalse)
2321 {
2322 channel_info.similarity=similarity;
2323 channel_info.x=x;
2324 channel_info.y=y;
2325 }
2326 switch (metric)
2327 {
2328 case PeakSignalToNoiseRatioMetric:
2329 {
2330 SetPixelRed(q,ClampToQuantum((double) QuantumRange*similarity));
2331 break;
2332 }
2333 default:
2334 {
2335 SetPixelRed(q,ClampToQuantum((double) QuantumRange*(1.0-similarity)));
2336 break;
2337 }
2338 }
2339 SetPixelGreen(q,GetPixelRed(q));
2340 SetPixelBlue(q,GetPixelRed(q));
2341 q++;
2342 }
2343#if defined(MAGICKCORE_OPENMP_SUPPORT)
2344 #pragma omp critical (MagickCore_SimilarityMetricImage)
2345#endif
2346 switch (metric)
2347 {
2348 case PeakSignalToNoiseRatioMetric:
2349 {
2350 if (similarity_threshold != DefaultSimilarityThreshold)
2351 if (channel_info.similarity >= similarity_threshold)
2352 threshold_trigger=MagickTrue;
2353 if (channel_info.similarity >= similarity_info.similarity)
2354 similarity_info=channel_info;
2355 break;
2356 }
2357 default:
2358 {
2359 if (similarity_threshold != DefaultSimilarityThreshold)
2360 if (channel_info.similarity < similarity_threshold)
2361 threshold_trigger=MagickTrue;
2362 if (channel_info.similarity < similarity_info.similarity)
2363 similarity_info=channel_info;
2364 break;
2365 }
2366 }
2367 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2368 status=MagickFalse;
2369 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2370 {
2371 MagickBooleanType
2372 proceed;
2373
2374 progress++;
2375 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2376 if (proceed == MagickFalse)
2377 status=MagickFalse;
2378 }
2379 }
2380 similarity_view=DestroyCacheView(similarity_view);
2381 switch (metric)
2382 {
2383 case NormalizedCrossCorrelationErrorMetric:
2384 {
2385 similarity_info.similarity=1.0-similarity_info.similarity;
2386 break;
2387 }
2388 default:
2389 break;
2390 }
2391 if (status == MagickFalse)
2392 similarity_image=DestroyImage(similarity_image);
2393 *similarity_metric=similarity_info.similarity;
2394 offset->x=similarity_info.x;
2395 offset->y=similarity_info.y;
2396 return(similarity_image);
2397}