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 fuzz;
448
449 MagickBooleanType
450 status;
451
452 size_t
453 columns,
454 rows;
455
456 ssize_t
457 j,
458 y;
459
460 /*
461 Compute the absolute difference in pixels between two images.
462 */
463 status=MagickTrue;
464 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
465 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
466 image_view=AcquireVirtualCacheView(image,exception);
467 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
468#if defined(MAGICKCORE_OPENMP_SUPPORT)
469 #pragma omp parallel for schedule(static) shared(status) \
470 magick_number_threads(image,image,rows,1)
471#endif
472 for (y=0; y < (ssize_t) rows; y++)
473 {
474 const IndexPacket
475 *magick_restrict indexes,
476 *magick_restrict reconstruct_indexes;
477
478 const PixelPacket
479 *magick_restrict p,
480 *magick_restrict q;
481
482 double
483 channel_distortion[CompositeChannels+1];
484
485 ssize_t
486 i,
487 x;
488
489 if (status == MagickFalse)
490 continue;
491 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
492 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
493 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
494 {
495 status=MagickFalse;
496 continue;
497 }
498 indexes=GetCacheViewVirtualIndexQueue(image_view);
499 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
500 (void) memset(channel_distortion,0,sizeof(channel_distortion));
501 for (x=0; x < (ssize_t) columns; x++)
502 {
503 double
504 Da,
505 delta,
506 Sa;
507
508 size_t
509 count = 0;
510
511 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
512 ((double) QuantumRange-(double) OpaqueOpacity));
513 Da=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(q) :
514 ((double) QuantumRange-(double) OpaqueOpacity));
515 if ((channel & RedChannel) != 0)
516 {
517 delta=Sa*(double) GetPixelRed(p)-Da*(double)
518 GetPixelRed(q);
519 if ((delta*delta) > fuzz)
520 {
521 channel_distortion[RedChannel]++;
522 count++;
523 }
524 }
525 if ((channel & GreenChannel) != 0)
526 {
527 delta=Sa*(double) GetPixelGreen(p)-Da*(double)
528 GetPixelGreen(q);
529 if ((delta*delta) > fuzz)
530 {
531 channel_distortion[GreenChannel]++;
532 count++;
533 }
534 }
535 if ((channel & BlueChannel) != 0)
536 {
537 delta=Sa*(double) GetPixelBlue(p)-Da*(double)
538 GetPixelBlue(q);
539 if ((delta*delta) > fuzz)
540 {
541 channel_distortion[BlueChannel]++;
542 count++;
543 }
544 }
545 if (((channel & OpacityChannel) != 0) &&
546 (image->matte != MagickFalse))
547 {
548 delta=(double) GetPixelOpacity(p)-(double)
549 GetPixelOpacity(q);
550 if ((delta*delta) > fuzz)
551 {
552 channel_distortion[OpacityChannel]++;
553 count++;
554 }
555 }
556 if (((channel & IndexChannel) != 0) &&
557 (image->colorspace == CMYKColorspace))
558 {
559 delta=Sa*(double) indexes[x]-Da*(double)
560 reconstruct_indexes[x];
561 if ((delta*delta) > fuzz)
562 {
563 channel_distortion[IndexChannel]++;
564 count++;
565 }
566 }
567 if (count != 0)
568 channel_distortion[CompositeChannels]++;
569 p++;
570 q++;
571 }
572#if defined(MAGICKCORE_OPENMP_SUPPORT)
573 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
574#endif
575 for (i=0; i <= (ssize_t) CompositeChannels; i++)
576 distortion[i]+=channel_distortion[i];
577 }
578 reconstruct_view=DestroyCacheView(reconstruct_view);
579 image_view=DestroyCacheView(image_view);
580 area=PerceptibleReciprocal((double) columns*rows);
581 for (j=0; j <= CompositeChannels; j++)
582 distortion[j]*=area;
583 return(status);
584}
585
586static MagickBooleanType GetFuzzDistortion(const Image *image,
587 const Image *reconstruct_image,const ChannelType channel,
588 double *distortion,ExceptionInfo *exception)
589{
590 CacheView
591 *image_view,
592 *reconstruct_view;
593
594 MagickBooleanType
595 status;
596
597 ssize_t
598 i;
599
600 size_t
601 columns,
602 rows;
603
604 ssize_t
605 y;
606
607 status=MagickTrue;
608 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
609 image_view=AcquireVirtualCacheView(image,exception);
610 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
611#if defined(MAGICKCORE_OPENMP_SUPPORT)
612 #pragma omp parallel for schedule(static) shared(status) \
613 magick_number_threads(image,image,rows,1)
614#endif
615 for (y=0; y < (ssize_t) rows; y++)
616 {
617 double
618 channel_distortion[CompositeChannels+1];
619
620 const IndexPacket
621 *magick_restrict indexes,
622 *magick_restrict reconstruct_indexes;
623
624 const PixelPacket
625 *magick_restrict p,
626 *magick_restrict q;
627
628 ssize_t
629 i,
630 x;
631
632 if (status == MagickFalse)
633 continue;
634 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
635 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
636 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
637 {
638 status=MagickFalse;
639 continue;
640 }
641 indexes=GetCacheViewVirtualIndexQueue(image_view);
642 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
643 (void) memset(channel_distortion,0,sizeof(channel_distortion));
644 for (x=0; x < (ssize_t) columns; x++)
645 {
646 MagickRealType
647 distance,
648 Da,
649 Sa;
650
651 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
652 ((double) QuantumRange-(double) OpaqueOpacity));
653 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
654 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
655 OpaqueOpacity));
656 if ((channel & RedChannel) != 0)
657 {
658 distance=QuantumScale*(Sa*(double) GetPixelRed(p)-Da*(double)
659 GetPixelRed(q));
660 channel_distortion[RedChannel]+=distance*distance;
661 channel_distortion[CompositeChannels]+=distance*distance;
662 }
663 if ((channel & GreenChannel) != 0)
664 {
665 distance=QuantumScale*(Sa*(double) GetPixelGreen(p)-Da*(double)
666 GetPixelGreen(q));
667 channel_distortion[GreenChannel]+=distance*distance;
668 channel_distortion[CompositeChannels]+=distance*distance;
669 }
670 if ((channel & BlueChannel) != 0)
671 {
672 distance=QuantumScale*(Sa*(double) GetPixelBlue(p)-Da*(double)
673 GetPixelBlue(q));
674 channel_distortion[BlueChannel]+=distance*distance;
675 channel_distortion[CompositeChannels]+=distance*distance;
676 }
677 if (((channel & OpacityChannel) != 0) && ((image->matte != MagickFalse) ||
678 (reconstruct_image->matte != MagickFalse)))
679 {
680 distance=QuantumScale*((image->matte != MagickFalse ? (double)
681 GetPixelOpacity(p) : (double) OpaqueOpacity)-
682 (reconstruct_image->matte != MagickFalse ?
683 (double) GetPixelOpacity(q): (double) OpaqueOpacity));
684 channel_distortion[OpacityChannel]+=distance*distance;
685 channel_distortion[CompositeChannels]+=distance*distance;
686 }
687 if (((channel & IndexChannel) != 0) &&
688 (image->colorspace == CMYKColorspace) &&
689 (reconstruct_image->colorspace == CMYKColorspace))
690 {
691 distance=QuantumScale*(Sa*(double) GetPixelIndex(indexes+x)-
692 Da*(double) GetPixelIndex(reconstruct_indexes+x));
693 channel_distortion[BlackChannel]+=distance*distance;
694 channel_distortion[CompositeChannels]+=distance*distance;
695 }
696 p++;
697 q++;
698 }
699#if defined(MAGICKCORE_OPENMP_SUPPORT)
700 #pragma omp critical (MagickCore_GetFuzzDistortion)
701#endif
702 for (i=0; i <= (ssize_t) CompositeChannels; i++)
703 distortion[i]+=channel_distortion[i];
704 }
705 reconstruct_view=DestroyCacheView(reconstruct_view);
706 image_view=DestroyCacheView(image_view);
707 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
708 for (i=0; i <= (ssize_t) CompositeChannels; i++)
709 distortion[i]/=((double) columns*rows);
710 distortion[CompositeChannels]=sqrt(distortion[CompositeChannels]);
711 return(status);
712}
713
714static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
715 const Image *reconstruct_image,const ChannelType channel,
716 double *distortion,ExceptionInfo *exception)
717{
718 CacheView
719 *image_view,
720 *reconstruct_view;
721
722 MagickBooleanType
723 status;
724
725 size_t
726 columns,
727 rows;
728
729 ssize_t
730 i,
731 y;
732
733 status=MagickTrue;
734 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
735 image_view=AcquireVirtualCacheView(image,exception);
736 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
737#if defined(MAGICKCORE_OPENMP_SUPPORT)
738 #pragma omp parallel for schedule(static) shared(status) \
739 magick_number_threads(image,image,rows,1)
740#endif
741 for (y=0; y < (ssize_t) rows; y++)
742 {
743 double
744 channel_distortion[CompositeChannels+1];
745
746 const IndexPacket
747 *magick_restrict indexes,
748 *magick_restrict reconstruct_indexes;
749
750 const PixelPacket
751 *magick_restrict p,
752 *magick_restrict q;
753
754 ssize_t
755 i,
756 x;
757
758 if (status == MagickFalse)
759 continue;
760 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
761 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
762 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
763 {
764 status=MagickFalse;
765 continue;
766 }
767 indexes=GetCacheViewVirtualIndexQueue(image_view);
768 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
769 (void) memset(channel_distortion,0,sizeof(channel_distortion));
770 for (x=0; x < (ssize_t) columns; x++)
771 {
772 MagickRealType
773 distance,
774 Da,
775 Sa;
776
777 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
778 ((double) QuantumRange-(double) OpaqueOpacity));
779 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
780 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
781 OpaqueOpacity));
782 if ((channel & RedChannel) != 0)
783 {
784 distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
785 (double) GetPixelRed(q));
786 channel_distortion[RedChannel]+=distance;
787 channel_distortion[CompositeChannels]+=distance;
788 }
789 if ((channel & GreenChannel) != 0)
790 {
791 distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
792 (double) GetPixelGreen(q));
793 channel_distortion[GreenChannel]+=distance;
794 channel_distortion[CompositeChannels]+=distance;
795 }
796 if ((channel & BlueChannel) != 0)
797 {
798 distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
799 (double) GetPixelBlue(q));
800 channel_distortion[BlueChannel]+=distance;
801 channel_distortion[CompositeChannels]+=distance;
802 }
803 if (((channel & OpacityChannel) != 0) &&
804 (image->matte != MagickFalse))
805 {
806 distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
807 GetPixelOpacity(q));
808 channel_distortion[OpacityChannel]+=distance;
809 channel_distortion[CompositeChannels]+=distance;
810 }
811 if (((channel & IndexChannel) != 0) &&
812 (image->colorspace == CMYKColorspace))
813 {
814 distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
815 (double) GetPixelIndex(reconstruct_indexes+x));
816 channel_distortion[BlackChannel]+=distance;
817 channel_distortion[CompositeChannels]+=distance;
818 }
819 p++;
820 q++;
821 }
822#if defined(MAGICKCORE_OPENMP_SUPPORT)
823 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
824#endif
825 for (i=0; i <= (ssize_t) CompositeChannels; i++)
826 distortion[i]+=channel_distortion[i];
827 }
828 reconstruct_view=DestroyCacheView(reconstruct_view);
829 image_view=DestroyCacheView(image_view);
830 for (i=0; i <= (ssize_t) CompositeChannels; i++)
831 distortion[i]/=((double) columns*rows);
832 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
833 return(status);
834}
835
836static MagickBooleanType GetMeanErrorPerPixel(Image *image,
837 const Image *reconstruct_image,const ChannelType channel,double *distortion,
838 ExceptionInfo *exception)
839{
840 CacheView
841 *image_view,
842 *reconstruct_view;
843
844 double
845 maximum_error = MagickMinimumValue,
846 mean_error = 0.0;
847
848 MagickBooleanType
849 status;
850
851 size_t
852 columns,
853 rows;
854
855 ssize_t
856 i,
857 y;
858
859 status=MagickTrue;
860 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
861 image_view=AcquireVirtualCacheView(image,exception);
862 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
863#if defined(MAGICKCORE_OPENMP_SUPPORT)
864 #pragma omp parallel for schedule(static) shared(status) \
865 magick_number_threads(image,image,rows,1)
866#endif
867 for (y=0; y < (ssize_t) rows; y++)
868 {
869 double
870 channel_distortion[CompositeChannels+1],
871 local_maximum = MinimumValue,
872 local_mean_error = 0.0;
873
874 const IndexPacket
875 *magick_restrict indexes,
876 *magick_restrict reconstruct_indexes;
877
878 const PixelPacket
879 *magick_restrict p,
880 *magick_restrict q;
881
882 ssize_t
883 i,
884 x;
885
886 if (status == MagickFalse)
887 continue;
888 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
889 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
890 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
891 {
892 status=MagickFalse;
893 continue;
894 }
895 indexes=GetCacheViewVirtualIndexQueue(image_view);
896 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
897 (void) memset(channel_distortion,0,sizeof(channel_distortion));
898 for (x=0; x < (ssize_t) columns; x++)
899 {
900 MagickRealType
901 distance,
902 Da,
903 Sa;
904
905 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
906 ((double) QuantumRange-(double) OpaqueOpacity));
907 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
908 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
909 OpaqueOpacity));
910 if ((channel & RedChannel) != 0)
911 {
912 distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
913 (double) GetPixelRed(q));
914 channel_distortion[RedChannel]+=distance;
915 channel_distortion[CompositeChannels]+=distance;
916 local_mean_error+=distance*distance;
917 if (distance > local_maximum)
918 local_maximum=distance;
919 }
920 if ((channel & GreenChannel) != 0)
921 {
922 distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
923 (double) GetPixelGreen(q));
924 channel_distortion[GreenChannel]+=distance;
925 channel_distortion[CompositeChannels]+=distance;
926 local_mean_error+=distance*distance;
927 if (distance > local_maximum)
928 local_maximum=distance;
929 }
930 if ((channel & BlueChannel) != 0)
931 {
932 distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
933 (double) GetPixelBlue(q));
934 channel_distortion[BlueChannel]+=distance;
935 channel_distortion[CompositeChannels]+=distance;
936 local_mean_error+=distance*distance;
937 if (distance > local_maximum)
938 local_maximum=distance;
939 }
940 if (((channel & OpacityChannel) != 0) &&
941 (image->matte != MagickFalse))
942 {
943 distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
944 GetPixelOpacity(q));
945 channel_distortion[OpacityChannel]+=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 & IndexChannel) != 0) &&
952 (image->colorspace == CMYKColorspace))
953 {
954 distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
955 (double) GetPixelIndex(reconstruct_indexes+x));
956 channel_distortion[BlackChannel]+=distance;
957 channel_distortion[CompositeChannels]+=distance;
958 local_mean_error+=distance*distance;
959 if (distance > local_maximum)
960 local_maximum=distance;
961 }
962 p++;
963 q++;
964 }
965#if defined(MAGICKCORE_OPENMP_SUPPORT)
966 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
967#endif
968 {
969 for (i=0; i <= (ssize_t) CompositeChannels; i++)
970 distortion[i]+=channel_distortion[i];
971 mean_error+=local_mean_error;
972 if (local_maximum > maximum_error)
973 maximum_error=local_maximum;
974 }
975 }
976 reconstruct_view=DestroyCacheView(reconstruct_view);
977 image_view=DestroyCacheView(image_view);
978 for (i=0; i <= (ssize_t) CompositeChannels; i++)
979 distortion[i]/=((double) columns*rows);
980 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
981 image->error.mean_error_per_pixel=QuantumRange*distortion[CompositeChannels];
982 image->error.normalized_mean_error=mean_error;
983 image->error.normalized_maximum_error=maximum_error;
984 return(status);
985}
986
987static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
988 const Image *reconstruct_image,const ChannelType channel,
989 double *distortion,ExceptionInfo *exception)
990{
991 CacheView
992 *image_view,
993 *reconstruct_view;
994
995 double
996 area;
997
998 MagickBooleanType
999 status;
1000
1001 size_t
1002 columns,
1003 rows;
1004
1005 ssize_t
1006 i,
1007 y;
1008
1009 status=MagickTrue;
1010 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
1011 image_view=AcquireVirtualCacheView(image,exception);
1012 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1013#if defined(MAGICKCORE_OPENMP_SUPPORT)
1014 #pragma omp parallel for schedule(static) shared(status) \
1015 magick_number_threads(image,image,rows,1)
1016#endif
1017 for (y=0; y < (ssize_t) rows; y++)
1018 {
1019 double
1020 channel_distortion[CompositeChannels+1];
1021
1022 const IndexPacket
1023 *magick_restrict indexes,
1024 *magick_restrict reconstruct_indexes;
1025
1026 const PixelPacket
1027 *magick_restrict p,
1028 *magick_restrict q;
1029
1030 ssize_t
1031 i,
1032 x;
1033
1034 if (status == MagickFalse)
1035 continue;
1036 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1037 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1038 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1039 {
1040 status=MagickFalse;
1041 continue;
1042 }
1043 indexes=GetCacheViewVirtualIndexQueue(image_view);
1044 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1045 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1046 for (x=0; x < (ssize_t) columns; x++)
1047 {
1048 MagickRealType
1049 distance,
1050 Da,
1051 Sa;
1052
1053 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1054 ((double) QuantumRange-(double) OpaqueOpacity));
1055 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1056 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
1057 OpaqueOpacity));
1058 if ((channel & RedChannel) != 0)
1059 {
1060 distance=QuantumScale*(Sa*(double) GetPixelRed(p)-Da*(double)
1061 GetPixelRed(q));
1062 channel_distortion[RedChannel]+=distance*distance;
1063 channel_distortion[CompositeChannels]+=distance*distance;
1064 }
1065 if ((channel & GreenChannel) != 0)
1066 {
1067 distance=QuantumScale*(Sa*(double) GetPixelGreen(p)-Da*(double)
1068 GetPixelGreen(q));
1069 channel_distortion[GreenChannel]+=distance*distance;
1070 channel_distortion[CompositeChannels]+=distance*distance;
1071 }
1072 if ((channel & BlueChannel) != 0)
1073 {
1074 distance=QuantumScale*(Sa*(double) GetPixelBlue(p)-Da*(double)
1075 GetPixelBlue(q));
1076 channel_distortion[BlueChannel]+=distance*distance;
1077 channel_distortion[CompositeChannels]+=distance*distance;
1078 }
1079 if (((channel & OpacityChannel) != 0) &&
1080 (image->matte != MagickFalse))
1081 {
1082 distance=QuantumScale*((double) GetPixelOpacity(p)-(double)
1083 GetPixelOpacity(q));
1084 channel_distortion[OpacityChannel]+=distance*distance;
1085 channel_distortion[CompositeChannels]+=distance*distance;
1086 }
1087 if (((channel & IndexChannel) != 0) &&
1088 (image->colorspace == CMYKColorspace) &&
1089 (reconstruct_image->colorspace == CMYKColorspace))
1090 {
1091 distance=QuantumScale*(Sa*(double) GetPixelIndex(indexes+x)-Da*
1092 (double) GetPixelIndex(reconstruct_indexes+x));
1093 channel_distortion[BlackChannel]+=distance*distance;
1094 channel_distortion[CompositeChannels]+=distance*distance;
1095 }
1096 p++;
1097 q++;
1098 }
1099#if defined(MAGICKCORE_OPENMP_SUPPORT)
1100 #pragma omp critical (MagickCore_GetMeanSquaredError)
1101#endif
1102 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1103 distortion[i]+=channel_distortion[i];
1104 }
1105 reconstruct_view=DestroyCacheView(reconstruct_view);
1106 image_view=DestroyCacheView(image_view);
1107 area=PerceptibleReciprocal((double) columns*rows);
1108 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1109 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1110 distortion[i]*=area;
1111 return(status);
1112}
1113
1114static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
1115 const Image *image,const Image *reconstruct_image,const ChannelType channel,
1116 double *distortion,ExceptionInfo *exception)
1117{
1118#define SimilarityImageTag "Similarity/Image"
1119
1120 CacheView
1121 *image_view,
1122 *reconstruct_view;
1123
1124 ChannelStatistics
1125 *image_statistics,
1126 *reconstruct_statistics;
1127
1128 double
1129 alpha_variance[CompositeChannels+1] = { 0.0 },
1130 area,
1131 beta_variance[CompositeChannels+1] = { 0.0 };
1132
1133 MagickBooleanType
1134 status;
1135
1136 MagickOffsetType
1137 progress;
1138
1139 size_t
1140 columns,
1141 rows;
1142
1143 ssize_t
1144 i,
1145 y;
1146
1147 /*
1148 Normalize to account for variation due to lighting and exposure condition.
1149 */
1150 image_statistics=GetImageChannelStatistics(image,exception);
1151 reconstruct_statistics=GetImageChannelStatistics(reconstruct_image,exception);
1152 if ((image_statistics == (ChannelStatistics *) NULL) ||
1153 (reconstruct_statistics == (ChannelStatistics *) NULL))
1154 {
1155 if (image_statistics != (ChannelStatistics *) NULL)
1156 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1157 image_statistics);
1158 if (reconstruct_statistics != (ChannelStatistics *) NULL)
1159 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1160 reconstruct_statistics);
1161 return(MagickFalse);
1162 }
1163 (void) memset(distortion,0,(CompositeChannels+1)*sizeof(*distortion));
1164 status=MagickTrue;
1165 progress=0;
1166 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
1167 image_view=AcquireVirtualCacheView(image,exception);
1168 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1169#if defined(MAGICKCORE_OPENMP_SUPPORT)
1170 #pragma omp parallel for schedule(static) shared(status) \
1171 magick_number_threads(image,image,rows,1)
1172#endif
1173 for (y=0; y < (ssize_t) rows; y++)
1174 {
1175 const IndexPacket
1176 *magick_restrict indexes,
1177 *magick_restrict reconstruct_indexes;
1178
1179 const PixelPacket
1180 *magick_restrict p,
1181 *magick_restrict q;
1182
1183 double
1184 channel_alpha_variance[CompositeChannels+1] = { 0.0 },
1185 channel_beta_variance[CompositeChannels+1] = { 0.0 },
1186 channel_distortion[CompositeChannels+1] = { 0.0 };
1187
1188 ssize_t
1189 x;
1190
1191 if (status == MagickFalse)
1192 continue;
1193 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1194 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1195 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1196 {
1197 status=MagickFalse;
1198 continue;
1199 }
1200 indexes=GetCacheViewVirtualIndexQueue(image_view);
1201 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1202 for (x=0; x < (ssize_t) columns; x++)
1203 {
1204 MagickRealType
1205 alpha,
1206 beta,
1207 Da,
1208 Sa;
1209
1210 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1211 (double) QuantumRange);
1212 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1213 (double) GetPixelAlpha(q) : (double) QuantumRange);
1214 if ((channel & RedChannel) != 0)
1215 {
1216 alpha=QuantumScale*(Sa*(double) GetPixelRed(p)-
1217 image_statistics[RedChannel].mean);
1218 beta=QuantumScale*(Da*(double) GetPixelRed(q)-
1219 reconstruct_statistics[RedChannel].mean);
1220 channel_distortion[RedChannel]+=alpha*beta;
1221 channel_alpha_variance[RedChannel]+=alpha*alpha;
1222 channel_beta_variance[RedChannel]+=beta*beta;
1223 }
1224 if ((channel & GreenChannel) != 0)
1225 {
1226 alpha=QuantumScale*(Sa*(double) GetPixelGreen(p)-
1227 image_statistics[GreenChannel].mean);
1228 beta=QuantumScale*(Da*(double) GetPixelGreen(q)-
1229 reconstruct_statistics[GreenChannel].mean);
1230 channel_distortion[GreenChannel]+=alpha*beta;
1231 channel_alpha_variance[GreenChannel]+=alpha*alpha;
1232 channel_beta_variance[GreenChannel]+=beta*beta;
1233 }
1234 if ((channel & BlueChannel) != 0)
1235 {
1236 alpha=QuantumScale*(Sa*(double) GetPixelBlue(p)-
1237 image_statistics[BlueChannel].mean);
1238 beta=QuantumScale*(Da*(double) GetPixelBlue(q)-
1239 reconstruct_statistics[BlueChannel].mean);
1240 channel_distortion[BlueChannel]+=alpha*beta;
1241 channel_alpha_variance[BlueChannel]+=alpha*alpha;
1242 channel_beta_variance[BlueChannel]+=beta*beta;
1243 }
1244 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1245 {
1246 alpha=QuantumScale*((double) GetPixelAlpha(p)-
1247 image_statistics[AlphaChannel].mean);
1248 beta=QuantumScale*((double) GetPixelAlpha(q)-
1249 reconstruct_statistics[AlphaChannel].mean);
1250 channel_distortion[OpacityChannel]+=alpha*beta;
1251 channel_alpha_variance[OpacityChannel]+=alpha*alpha;
1252 channel_beta_variance[OpacityChannel]+=beta*beta;
1253 }
1254 if (((channel & IndexChannel) != 0) &&
1255 (image->colorspace == CMYKColorspace) &&
1256 (reconstruct_image->colorspace == CMYKColorspace))
1257 {
1258 alpha=QuantumScale*(Sa*(double) GetPixelIndex(indexes+x)-
1259 image_statistics[BlackChannel].mean);
1260 beta=QuantumScale*(Da*(double) GetPixelIndex(reconstruct_indexes+
1261 x)-reconstruct_statistics[BlackChannel].mean);
1262 channel_distortion[BlackChannel]+=alpha*beta;
1263 channel_alpha_variance[BlackChannel]+=alpha*alpha;
1264 channel_beta_variance[BlackChannel]+=beta*beta;
1265 }
1266 p++;
1267 q++;
1268 }
1269#if defined(MAGICKCORE_OPENMP_SUPPORT)
1270 #pragma omp critical (GetNormalizedCrossCorrelationDistortion)
1271#endif
1272 {
1273 ssize_t
1274 j;
1275
1276 for (j=0; j < (ssize_t) CompositeChannels; j++)
1277 {
1278 distortion[j]+=channel_distortion[j];
1279 alpha_variance[j]+=channel_alpha_variance[j];
1280 beta_variance[j]+=channel_beta_variance[j];
1281 }
1282 }
1283 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1284 {
1285 MagickBooleanType
1286 proceed;
1287
1288#if defined(MAGICKCORE_OPENMP_SUPPORT)
1289 #pragma omp atomic
1290#endif
1291 progress++;
1292 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1293 if (proceed == MagickFalse)
1294 status=MagickFalse;
1295 }
1296 }
1297 reconstruct_view=DestroyCacheView(reconstruct_view);
1298 image_view=DestroyCacheView(image_view);
1299 /*
1300 Divide by the standard deviation.
1301 */
1302 area=PerceptibleReciprocal((double) columns*rows);
1303 for (i=0; i < (ssize_t) CompositeChannels; i++)
1304 {
1305 distortion[i]*=area;
1306 alpha_variance[i]*=area;
1307 beta_variance[i]*=area;
1308 distortion[i]*=PerceptibleReciprocal(sqrt(alpha_variance[i]*
1309 beta_variance[i]));
1310 distortion[i]=1.0-distortion[i];
1311 if (fabs(distortion[i]) < MagickEpsilon)
1312 distortion[i]=0.0;
1313 distortion[CompositeChannels]+=distortion[i];
1314 }
1315 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1316 /*
1317 Free resources.
1318 */
1319 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1320 reconstruct_statistics);
1321 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1322 image_statistics);
1323 return(status);
1324}
1325
1326static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1327 const Image *reconstruct_image,const ChannelType channel,
1328 double *distortion,ExceptionInfo *exception)
1329{
1330 CacheView
1331 *image_view,
1332 *reconstruct_view;
1333
1334 MagickBooleanType
1335 status;
1336
1337 size_t
1338 columns,
1339 rows;
1340
1341 ssize_t
1342 y;
1343
1344 status=MagickTrue;
1345 (void) memset(distortion,0,(CompositeChannels+1)*sizeof(*distortion));
1346 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
1347 image_view=AcquireVirtualCacheView(image,exception);
1348 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1349#if defined(MAGICKCORE_OPENMP_SUPPORT)
1350 #pragma omp parallel for schedule(static) shared(status) \
1351 magick_number_threads(image,image,rows,1)
1352#endif
1353 for (y=0; y < (ssize_t) rows; y++)
1354 {
1355 double
1356 channel_distortion[CompositeChannels+1];
1357
1358 const IndexPacket
1359 *magick_restrict indexes,
1360 *magick_restrict reconstruct_indexes;
1361
1362 const PixelPacket
1363 *magick_restrict p,
1364 *magick_restrict q;
1365
1366 ssize_t
1367 i,
1368 x;
1369
1370 if (status == MagickFalse)
1371 continue;
1372 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1373 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1374 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1375 {
1376 status=MagickFalse;
1377 continue;
1378 }
1379 indexes=GetCacheViewVirtualIndexQueue(image_view);
1380 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1381 (void) memset(channel_distortion,0,(CompositeChannels+1)*
1382 sizeof(*channel_distortion));
1383 for (x=0; x < (ssize_t) columns; x++)
1384 {
1385 MagickRealType
1386 distance,
1387 Da,
1388 Sa;
1389
1390 Sa=QuantumScale*(image->matte != MagickFalse ? (double) GetPixelAlpha(p) :
1391 ((double) QuantumRange-(double) OpaqueOpacity));
1392 Da=QuantumScale*(reconstruct_image->matte != MagickFalse ?
1393 (double) GetPixelAlpha(q) : ((double) QuantumRange-(double)
1394 OpaqueOpacity));
1395 if ((channel & RedChannel) != 0)
1396 {
1397 distance=QuantumScale*fabs(Sa*(double) GetPixelRed(p)-Da*
1398 (double) GetPixelRed(q));
1399 if (distance > channel_distortion[RedChannel])
1400 channel_distortion[RedChannel]=distance;
1401 if (distance > channel_distortion[CompositeChannels])
1402 channel_distortion[CompositeChannels]=distance;
1403 }
1404 if ((channel & GreenChannel) != 0)
1405 {
1406 distance=QuantumScale*fabs(Sa*(double) GetPixelGreen(p)-Da*
1407 (double) GetPixelGreen(q));
1408 if (distance > channel_distortion[GreenChannel])
1409 channel_distortion[GreenChannel]=distance;
1410 if (distance > channel_distortion[CompositeChannels])
1411 channel_distortion[CompositeChannels]=distance;
1412 }
1413 if ((channel & BlueChannel) != 0)
1414 {
1415 distance=QuantumScale*fabs(Sa*(double) GetPixelBlue(p)-Da*
1416 (double) GetPixelBlue(q));
1417 if (distance > channel_distortion[BlueChannel])
1418 channel_distortion[BlueChannel]=distance;
1419 if (distance > channel_distortion[CompositeChannels])
1420 channel_distortion[CompositeChannels]=distance;
1421 }
1422 if (((channel & OpacityChannel) != 0) &&
1423 (image->matte != MagickFalse))
1424 {
1425 distance=QuantumScale*fabs((double) GetPixelOpacity(p)-(double)
1426 GetPixelOpacity(q));
1427 if (distance > channel_distortion[OpacityChannel])
1428 channel_distortion[OpacityChannel]=distance;
1429 if (distance > channel_distortion[CompositeChannels])
1430 channel_distortion[CompositeChannels]=distance;
1431 }
1432 if (((channel & IndexChannel) != 0) &&
1433 (image->colorspace == CMYKColorspace) &&
1434 (reconstruct_image->colorspace == CMYKColorspace))
1435 {
1436 distance=QuantumScale*fabs(Sa*(double) GetPixelIndex(indexes+x)-Da*
1437 (double) GetPixelIndex(reconstruct_indexes+x));
1438 if (distance > channel_distortion[BlackChannel])
1439 channel_distortion[BlackChannel]=distance;
1440 if (distance > channel_distortion[CompositeChannels])
1441 channel_distortion[CompositeChannels]=distance;
1442 }
1443 p++;
1444 q++;
1445 }
1446#if defined(MAGICKCORE_OPENMP_SUPPORT)
1447 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1448#endif
1449 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1450 if (channel_distortion[i] > distortion[i])
1451 distortion[i]=channel_distortion[i];
1452 }
1453 reconstruct_view=DestroyCacheView(reconstruct_view);
1454 image_view=DestroyCacheView(image_view);
1455 return(status);
1456}
1457
1458static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1459 const Image *reconstruct_image,const ChannelType channel,
1460 double *distortion,ExceptionInfo *exception)
1461{
1462 MagickBooleanType
1463 status;
1464
1465 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1466 exception);
1467 if ((channel & RedChannel) != 0)
1468 distortion[RedChannel]=10.0*PerceptibleLog10(PerceptibleReciprocal(
1469 distortion[RedChannel]))/MagickPSNRDistortion;
1470 if ((channel & GreenChannel) != 0)
1471 distortion[GreenChannel]=10.0*PerceptibleLog10(PerceptibleReciprocal(
1472 distortion[GreenChannel]))/MagickPSNRDistortion;
1473 if ((channel & BlueChannel) != 0)
1474 distortion[BlueChannel]=10.0*PerceptibleLog10(PerceptibleReciprocal(
1475 distortion[BlueChannel]))/MagickPSNRDistortion;
1476 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
1477 distortion[OpacityChannel]=10.0*PerceptibleLog10(PerceptibleReciprocal(
1478 distortion[OpacityChannel]))/MagickPSNRDistortion;
1479 if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace))
1480 distortion[BlackChannel]=10.0*PerceptibleLog10(PerceptibleReciprocal(
1481 distortion[BlackChannel]))/MagickPSNRDistortion;
1482 distortion[CompositeChannels]=10.0*PerceptibleLog10(PerceptibleReciprocal(
1483 distortion[CompositeChannels]))/MagickPSNRDistortion;
1484 return(status);
1485}
1486
1487static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1488 const Image *reconstruct_image,const ChannelType channel,double *distortion,
1489 ExceptionInfo *exception)
1490{
1491#define PHASHNormalizationFactor 389.373723242
1492
1493 ChannelPerceptualHash
1494 *image_phash,
1495 *reconstruct_phash;
1496
1497 double
1498 difference;
1499
1500 ssize_t
1501 i;
1502
1503 /*
1504 Compute perceptual hash in the sRGB colorspace.
1505 */
1506 image_phash=GetImageChannelPerceptualHash(image,exception);
1507 if (image_phash == (ChannelPerceptualHash *) NULL)
1508 return(MagickFalse);
1509 reconstruct_phash=GetImageChannelPerceptualHash(reconstruct_image,exception);
1510 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1511 {
1512 image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1513 return(MagickFalse);
1514 }
1515 for (i=0; i < MaximumNumberOfImageMoments; i++)
1516 {
1517 /*
1518 Compute sum of moment differences squared.
1519 */
1520 if ((channel & RedChannel) != 0)
1521 {
1522 difference=reconstruct_phash[RedChannel].P[i]-
1523 image_phash[RedChannel].P[i];
1524 if (IsNaN(difference) != 0)
1525 difference=0.0;
1526 distortion[RedChannel]+=difference*difference;
1527 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1528 }
1529 if ((channel & GreenChannel) != 0)
1530 {
1531 difference=reconstruct_phash[GreenChannel].P[i]-
1532 image_phash[GreenChannel].P[i];
1533 if (IsNaN(difference) != 0)
1534 difference=0.0;
1535 distortion[GreenChannel]+=difference*difference;
1536 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1537 }
1538 if ((channel & BlueChannel) != 0)
1539 {
1540 difference=reconstruct_phash[BlueChannel].P[i]-
1541 image_phash[BlueChannel].P[i];
1542 if (IsNaN(difference) != 0)
1543 difference=0.0;
1544 distortion[BlueChannel]+=difference*difference;
1545 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1546 }
1547 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1548 (reconstruct_image->matte != MagickFalse))
1549 {
1550 difference=reconstruct_phash[OpacityChannel].P[i]-
1551 image_phash[OpacityChannel].P[i];
1552 if (IsNaN(difference) != 0)
1553 difference=0.0;
1554 distortion[OpacityChannel]+=difference*difference;
1555 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1556 }
1557 if (((channel & IndexChannel) != 0) &&
1558 (image->colorspace == CMYKColorspace) &&
1559 (reconstruct_image->colorspace == CMYKColorspace))
1560 {
1561 difference=reconstruct_phash[IndexChannel].P[i]-
1562 image_phash[IndexChannel].P[i];
1563 if (IsNaN(difference) != 0)
1564 difference=0.0;
1565 distortion[IndexChannel]+=difference*difference;
1566 distortion[CompositeChannels]+=QuantumScale*difference*difference;
1567 }
1568 }
1569 /*
1570 Compute perceptual hash in the HCLP colorspace.
1571 */
1572 for (i=0; i < MaximumNumberOfImageMoments; i++)
1573 {
1574 /*
1575 Compute sum of moment differences squared.
1576 */
1577 if ((channel & RedChannel) != 0)
1578 {
1579 difference=reconstruct_phash[RedChannel].Q[i]-
1580 image_phash[RedChannel].Q[i];
1581 if (IsNaN(difference) != 0)
1582 difference=0.0;
1583 distortion[RedChannel]+=difference*difference;
1584 distortion[CompositeChannels]+=difference*difference;
1585 }
1586 if ((channel & GreenChannel) != 0)
1587 {
1588 difference=reconstruct_phash[GreenChannel].Q[i]-
1589 image_phash[GreenChannel].Q[i];
1590 if (IsNaN(difference) != 0)
1591 difference=0.0;
1592 distortion[GreenChannel]+=difference*difference;
1593 distortion[CompositeChannels]+=difference*difference;
1594 }
1595 if ((channel & BlueChannel) != 0)
1596 {
1597 difference=reconstruct_phash[BlueChannel].Q[i]-
1598 image_phash[BlueChannel].Q[i];
1599 distortion[BlueChannel]+=difference*difference;
1600 distortion[CompositeChannels]+=difference*difference;
1601 }
1602 if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse) &&
1603 (reconstruct_image->matte != MagickFalse))
1604 {
1605 difference=reconstruct_phash[OpacityChannel].Q[i]-
1606 image_phash[OpacityChannel].Q[i];
1607 if (IsNaN(difference) != 0)
1608 difference=0.0;
1609 distortion[OpacityChannel]+=difference*difference;
1610 distortion[CompositeChannels]+=difference*difference;
1611 }
1612 if (((channel & IndexChannel) != 0) &&
1613 (image->colorspace == CMYKColorspace) &&
1614 (reconstruct_image->colorspace == CMYKColorspace))
1615 {
1616 difference=reconstruct_phash[IndexChannel].Q[i]-
1617 image_phash[IndexChannel].Q[i];
1618 if (IsNaN(difference) != 0)
1619 difference=0.0;
1620 distortion[IndexChannel]+=difference*difference;
1621 distortion[CompositeChannels]+=difference*difference;
1622 }
1623 }
1624 distortion[CompositeChannels]/=(double) GetNumberChannels(image,channel);
1625 for (i=0; i <= (ssize_t) CompositeChannels; i++)
1626 distortion[i]/=PHASHNormalizationFactor;
1627 /*
1628 Free resources.
1629 */
1630 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1631 reconstruct_phash);
1632 image_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(image_phash);
1633 return(MagickTrue);
1634}
1635
1636static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1637 const Image *reconstruct_image,const ChannelType channel,double *distortion,
1638 ExceptionInfo *exception)
1639{
1640 MagickBooleanType
1641 status;
1642
1643 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,distortion,
1644 exception);
1645 if ((channel & RedChannel) != 0)
1646 distortion[RedChannel]=sqrt(distortion[RedChannel]);
1647 if ((channel & GreenChannel) != 0)
1648 distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
1649 if ((channel & BlueChannel) != 0)
1650 distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
1651 if (((channel & OpacityChannel) != 0) &&
1652 (image->matte != MagickFalse))
1653 distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
1654 if (((channel & IndexChannel) != 0) &&
1655 (image->colorspace == CMYKColorspace))
1656 distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
1657 distortion[CompositeChannels]=sqrt(distortion[CompositeChannels]);
1658 return(status);
1659}
1660
1661MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1662 const Image *reconstruct_image,const ChannelType channel,
1663 const MetricType metric,double *distortion,ExceptionInfo *exception)
1664{
1665 double
1666 *channel_distortion;
1667
1668 MagickBooleanType
1669 status;
1670
1671 size_t
1672 length;
1673
1674 assert(image != (Image *) NULL);
1675 assert(image->signature == MagickCoreSignature);
1676 assert(reconstruct_image != (const Image *) NULL);
1677 assert(reconstruct_image->signature == MagickCoreSignature);
1678 assert(distortion != (double *) NULL);
1679 if (IsEventLogging() != MagickFalse)
1680 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1681 *distortion=0.0;
1682 if (metric != PerceptualHashErrorMetric)
1683 if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1684 ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1685 /*
1686 Get image distortion.
1687 */
1688 length=CompositeChannels+1UL;
1689 channel_distortion=(double *) AcquireQuantumMemory(length,
1690 sizeof(*channel_distortion));
1691 if (channel_distortion == (double *) NULL)
1692 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1693 (void) memset(channel_distortion,0,length*sizeof(*channel_distortion));
1694 switch (metric)
1695 {
1696 case AbsoluteErrorMetric:
1697 {
1698 status=GetAbsoluteDistortion(image,reconstruct_image,channel,
1699 channel_distortion,exception);
1700 break;
1701 }
1702 case FuzzErrorMetric:
1703 {
1704 status=GetFuzzDistortion(image,reconstruct_image,channel,
1705 channel_distortion,exception);
1706 break;
1707 }
1708 case MeanAbsoluteErrorMetric:
1709 {
1710 status=GetMeanAbsoluteDistortion(image,reconstruct_image,channel,
1711 channel_distortion,exception);
1712 break;
1713 }
1714 case MeanErrorPerPixelMetric:
1715 {
1716 status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1717 channel_distortion,exception);
1718 break;
1719 }
1720 case MeanSquaredErrorMetric:
1721 {
1722 status=GetMeanSquaredDistortion(image,reconstruct_image,channel,
1723 channel_distortion,exception);
1724 break;
1725 }
1726 case NormalizedCrossCorrelationErrorMetric:
1727 {
1728 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1729 channel,channel_distortion,exception);
1730 break;
1731 }
1732 case PeakAbsoluteErrorMetric:
1733 {
1734 status=GetPeakAbsoluteDistortion(image,reconstruct_image,channel,
1735 channel_distortion,exception);
1736 break;
1737 }
1738 case PeakSignalToNoiseRatioMetric:
1739 {
1740 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1741 channel_distortion,exception);
1742 break;
1743 }
1744 case PerceptualHashErrorMetric:
1745 {
1746 status=GetPerceptualHashDistortion(image,reconstruct_image,channel,
1747 channel_distortion,exception);
1748 break;
1749 }
1750 case RootMeanSquaredErrorMetric:
1751 case UndefinedErrorMetric:
1752 default:
1753 {
1754 status=GetRootMeanSquaredDistortion(image,reconstruct_image,channel,
1755 channel_distortion,exception);
1756 break;
1757 }
1758 }
1759 *distortion=channel_distortion[CompositeChannels];
1760 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1761 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1762 *distortion);
1763 return(status);
1764}
1765
1766/*
1767%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1768% %
1769% %
1770% %
1771% 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 %
1772% %
1773% %
1774% %
1775%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1776%
1777% GetImageChannelDistortions() compares the image channels of an image to a
1778% reconstructed image and returns the specified distortion metric for each
1779% channel.
1780%
1781% The format of the GetImageChannelDistortions method is:
1782%
1783% double *GetImageChannelDistortions(const Image *image,
1784% const Image *reconstruct_image,const MetricType metric,
1785% ExceptionInfo *exception)
1786%
1787% A description of each parameter follows:
1788%
1789% o image: the image.
1790%
1791% o reconstruct_image: the reconstruct image.
1792%
1793% o metric: the metric.
1794%
1795% o exception: return any errors or warnings in this structure.
1796%
1797*/
1798MagickExport double *GetImageChannelDistortions(Image *image,
1799 const Image *reconstruct_image,const MetricType metric,
1800 ExceptionInfo *exception)
1801{
1802 double
1803 *channel_distortion;
1804
1805 MagickBooleanType
1806 status;
1807
1808 size_t
1809 length;
1810
1811 assert(image != (Image *) NULL);
1812 assert(image->signature == MagickCoreSignature);
1813 assert(reconstruct_image != (const Image *) NULL);
1814 assert(reconstruct_image->signature == MagickCoreSignature);
1815 if (IsEventLogging() != MagickFalse)
1816 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1817 if (metric != PerceptualHashErrorMetric)
1818 if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1819 {
1820 (void) ThrowMagickException(&image->exception,GetMagickModule(),
1821 ImageError,"ImageMorphologyDiffers","`%s'",image->filename);
1822 return((double *) NULL);
1823 }
1824 /*
1825 Get image distortion.
1826 */
1827 length=CompositeChannels+1UL;
1828 channel_distortion=(double *) AcquireQuantumMemory(length,
1829 sizeof(*channel_distortion));
1830 if (channel_distortion == (double *) NULL)
1831 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1832 (void) memset(channel_distortion,0,length*
1833 sizeof(*channel_distortion));
1834 status=MagickTrue;
1835 switch (metric)
1836 {
1837 case AbsoluteErrorMetric:
1838 {
1839 status=GetAbsoluteDistortion(image,reconstruct_image,CompositeChannels,
1840 channel_distortion,exception);
1841 break;
1842 }
1843 case FuzzErrorMetric:
1844 {
1845 status=GetFuzzDistortion(image,reconstruct_image,CompositeChannels,
1846 channel_distortion,exception);
1847 break;
1848 }
1849 case MeanAbsoluteErrorMetric:
1850 {
1851 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1852 CompositeChannels,channel_distortion,exception);
1853 break;
1854 }
1855 case MeanErrorPerPixelMetric:
1856 {
1857 status=GetMeanErrorPerPixel(image,reconstruct_image,CompositeChannels,
1858 channel_distortion,exception);
1859 break;
1860 }
1861 case MeanSquaredErrorMetric:
1862 {
1863 status=GetMeanSquaredDistortion(image,reconstruct_image,CompositeChannels,
1864 channel_distortion,exception);
1865 break;
1866 }
1867 case NormalizedCrossCorrelationErrorMetric:
1868 {
1869 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1870 CompositeChannels,channel_distortion,exception);
1871 break;
1872 }
1873 case PeakAbsoluteErrorMetric:
1874 {
1875 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1876 CompositeChannels,channel_distortion,exception);
1877 break;
1878 }
1879 case PeakSignalToNoiseRatioMetric:
1880 {
1881 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1882 CompositeChannels,channel_distortion,exception);
1883 break;
1884 }
1885 case PerceptualHashErrorMetric:
1886 {
1887 status=GetPerceptualHashDistortion(image,reconstruct_image,
1888 CompositeChannels,channel_distortion,exception);
1889 break;
1890 }
1891 case RootMeanSquaredErrorMetric:
1892 case UndefinedErrorMetric:
1893 default:
1894 {
1895 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1896 CompositeChannels,channel_distortion,exception);
1897 break;
1898 }
1899 }
1900 if (status == MagickFalse)
1901 {
1902 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1903 return((double *) NULL);
1904 }
1905 return(channel_distortion);
1906}
1907
1908/*
1909%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1910% %
1911% %
1912% %
1913% I s I m a g e s E q u a l %
1914% %
1915% %
1916% %
1917%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1918%
1919% IsImagesEqual() measures the difference between colors at each pixel
1920% location of two images. A value other than 0 means the colors match
1921% exactly. Otherwise an error measure is computed by summing over all
1922% pixels in an image the distance squared in RGB space between each image
1923% pixel and its corresponding pixel in the reconstruct image. The error
1924% measure is assigned to these image members:
1925%
1926% o mean_error_per_pixel: The mean error for any single pixel in
1927% the image.
1928%
1929% o normalized_mean_error: The normalized mean quantization error for
1930% any single pixel in the image. This distance measure is normalized to
1931% a range between 0 and 1. It is independent of the range of red, green,
1932% and blue values in the image.
1933%
1934% o normalized_maximum_error: The normalized maximum quantization
1935% error for any single pixel in the image. This distance measure is
1936% normalized to a range between 0 and 1. It is independent of the range
1937% of red, green, and blue values in your image.
1938%
1939% A small normalized mean square error, accessed as
1940% image->normalized_mean_error, suggests the images are very similar in
1941% spatial layout and color.
1942%
1943% The format of the IsImagesEqual method is:
1944%
1945% MagickBooleanType IsImagesEqual(Image *image,
1946% const Image *reconstruct_image)
1947%
1948% A description of each parameter follows.
1949%
1950% o image: the image.
1951%
1952% o reconstruct_image: the reconstruct image.
1953%
1954*/
1955MagickExport MagickBooleanType IsImagesEqual(Image *image,
1956 const Image *reconstruct_image)
1957{
1958 CacheView
1959 *image_view,
1960 *reconstruct_view;
1961
1962 ExceptionInfo
1963 *exception;
1964
1965 MagickBooleanType
1966 status;
1967
1968 MagickRealType
1969 area,
1970 gamma,
1971 maximum_error,
1972 mean_error,
1973 mean_error_per_pixel;
1974
1975 size_t
1976 columns,
1977 rows;
1978
1979 ssize_t
1980 y;
1981
1982 assert(image != (Image *) NULL);
1983 assert(image->signature == MagickCoreSignature);
1984 assert(reconstruct_image != (const Image *) NULL);
1985 assert(reconstruct_image->signature == MagickCoreSignature);
1986 exception=(&image->exception);
1987 if (ValidateImageMorphology(image,reconstruct_image) == MagickFalse)
1988 ThrowBinaryException(ImageError,"ImageMorphologyDiffers",image->filename);
1989 area=0.0;
1990 maximum_error=0.0;
1991 mean_error_per_pixel=0.0;
1992 mean_error=0.0;
1993 SetImageDistortionBounds(image,reconstruct_image,&columns,&rows);
1994 image_view=AcquireVirtualCacheView(image,exception);
1995 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1996 for (y=0; y < (ssize_t) rows; y++)
1997 {
1998 const IndexPacket
1999 *magick_restrict indexes,
2000 *magick_restrict reconstruct_indexes;
2001
2002 const PixelPacket
2003 *magick_restrict p,
2004 *magick_restrict q;
2005
2006 ssize_t
2007 x;
2008
2009 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2010 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2011 if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
2012 break;
2013 indexes=GetCacheViewVirtualIndexQueue(image_view);
2014 reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
2015 for (x=0; x < (ssize_t) columns; x++)
2016 {
2017 MagickRealType
2018 distance;
2019
2020 distance=fabs((double) GetPixelRed(p)-(double) GetPixelRed(q));
2021 mean_error_per_pixel+=distance;
2022 mean_error+=distance*distance;
2023 if (distance > maximum_error)
2024 maximum_error=distance;
2025 area++;
2026 distance=fabs((double) GetPixelGreen(p)-(double) GetPixelGreen(q));
2027 mean_error_per_pixel+=distance;
2028 mean_error+=distance*distance;
2029 if (distance > maximum_error)
2030 maximum_error=distance;
2031 area++;
2032 distance=fabs((double) GetPixelBlue(p)-(double) GetPixelBlue(q));
2033 mean_error_per_pixel+=distance;
2034 mean_error+=distance*distance;
2035 if (distance > maximum_error)
2036 maximum_error=distance;
2037 area++;
2038 if (image->matte != MagickFalse)
2039 {
2040 distance=fabs((double) GetPixelOpacity(p)-(double)
2041 GetPixelOpacity(q));
2042 mean_error_per_pixel+=distance;
2043 mean_error+=distance*distance;
2044 if (distance > maximum_error)
2045 maximum_error=distance;
2046 area++;
2047 }
2048 if ((image->colorspace == CMYKColorspace) &&
2049 (reconstruct_image->colorspace == CMYKColorspace))
2050 {
2051 distance=fabs((double) GetPixelIndex(indexes+x)-(double)
2052 GetPixelIndex(reconstruct_indexes+x));
2053 mean_error_per_pixel+=distance;
2054 mean_error+=distance*distance;
2055 if (distance > maximum_error)
2056 maximum_error=distance;
2057 area++;
2058 }
2059 p++;
2060 q++;
2061 }
2062 }
2063 reconstruct_view=DestroyCacheView(reconstruct_view);
2064 image_view=DestroyCacheView(image_view);
2065 gamma=PerceptibleReciprocal(area);
2066 image->error.mean_error_per_pixel=gamma*mean_error_per_pixel;
2067 image->error.normalized_mean_error=gamma*QuantumScale*QuantumScale*mean_error;
2068 image->error.normalized_maximum_error=QuantumScale*maximum_error;
2069 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2070 return(status);
2071}
2072
2073/*
2074%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2075% %
2076% %
2077% %
2078% S i m i l a r i t y I m a g e %
2079% %
2080% %
2081% %
2082%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2083%
2084% SimilarityImage() compares the reference image of the image and returns the
2085% best match offset. In addition, it returns a similarity image such that an
2086% exact match location is completely white and if none of the pixels match,
2087% black, otherwise some gray level in-between.
2088%
2089% The format of the SimilarityImageImage method is:
2090%
2091% Image *SimilarityImage(const Image *image,const Image *reference,
2092% RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2093%
2094% A description of each parameter follows:
2095%
2096% o image: the image.
2097%
2098% o reference: find an area of the image that closely resembles this image.
2099%
2100% o the best match offset of the reference image within the image.
2101%
2102% o similarity: the computed similarity between the images.
2103%
2104% o exception: return any errors or warnings in this structure.
2105%
2106*/
2107
2108static double GetSimilarityMetric(const Image *image,const Image *reference,
2109 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2110 ExceptionInfo *exception)
2111{
2112 double
2113 distortion;
2114
2115 Image
2116 *similarity_image;
2117
2118 MagickBooleanType
2119 status;
2120
2121 RectangleInfo
2122 geometry;
2123
2124 SetGeometry(reference,&geometry);
2125 geometry.x=x_offset;
2126 geometry.y=y_offset;
2127 similarity_image=CropImage(image,&geometry,exception);
2128 if (similarity_image == (Image *) NULL)
2129 return(0.0);
2130 distortion=0.0;
2131 status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2132 exception);
2133 (void) status;
2134 similarity_image=DestroyImage(similarity_image);
2135 return(distortion);
2136}
2137
2138MagickExport Image *SimilarityImage(Image *image,const Image *reference,
2139 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2140{
2141 Image
2142 *similarity_image;
2143
2144 similarity_image=SimilarityMetricImage(image,reference,
2145 RootMeanSquaredErrorMetric,offset,similarity_metric,exception);
2146 return(similarity_image);
2147}
2148
2149MagickExport Image *SimilarityMetricImage(Image *image,const Image *reconstruct,
2150 const MetricType metric,RectangleInfo *offset,double *similarity_metric,
2151 ExceptionInfo *exception)
2152{
2153#define SimilarityImageTag "Similarity/Image"
2154
2155 typedef struct
2156 {
2157 double
2158 similarity;
2159
2160 ssize_t
2161 x,
2162 y;
2163 } SimilarityInfo;
2164
2165 CacheView
2166 *similarity_view;
2167
2168 const char
2169 *artifact;
2170
2171 double
2172 similarity_threshold;
2173
2174 Image
2175 *similarity_image = (Image *) NULL;
2176
2177 MagickBooleanType
2178 status;
2179
2180 MagickOffsetType
2181 progress;
2182
2183 SimilarityInfo
2184 similarity_info = { 0 };
2185
2186 ssize_t
2187 y;
2188
2189 assert(image != (const Image *) NULL);
2190 assert(image->signature == MagickCoreSignature);
2191 assert(exception != (ExceptionInfo *) NULL);
2192 assert(exception->signature == MagickCoreSignature);
2193 assert(offset != (RectangleInfo *) NULL);
2194 if (IsEventLogging() != MagickFalse)
2195 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2196 SetGeometry(reconstruct,offset);
2197 *similarity_metric=0.0;
2198 offset->x=0;
2199 offset->y=0;
2200 if (ValidateImageMorphology(image,reconstruct) == MagickFalse)
2201 ThrowImageException(ImageError,"ImageMorphologyDiffers");
2202 if ((image->columns < reconstruct->columns) ||
2203 (image->rows < reconstruct->rows))
2204 {
2205 (void) ThrowMagickException(&image->exception,GetMagickModule(),
2206 OptionWarning,"GeometryDoesNotContainImage","`%s'",image->filename);
2207 return((Image *) NULL);
2208 }
2209 similarity_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2210 exception);
2211 if (similarity_image == (Image *) NULL)
2212 return((Image *) NULL);
2213 similarity_image->depth=32;
2214 similarity_image->colorspace=GRAYColorspace;
2215 similarity_image->matte=MagickFalse;
2216 status=SetImageStorageClass(similarity_image,DirectClass);
2217 if (status == MagickFalse)
2218 {
2219 InheritException(exception,&similarity_image->exception);
2220 return(DestroyImage(similarity_image));
2221 }
2222 /*
2223 Measure similarity of reconstruction image against image.
2224 */
2225 similarity_threshold=DefaultSimilarityThreshold;
2226 artifact=GetImageArtifact(image,"compare:similarity-threshold");
2227 if (artifact != (const char *) NULL)
2228 similarity_threshold=StringToDouble(artifact,(char **) NULL);
2229 status=MagickTrue;
2230 similarity_info.similarity=GetSimilarityMetric(image,reconstruct,metric,
2231 similarity_info.x,similarity_info.y,exception);
2232 progress=0;
2233 similarity_view=AcquireVirtualCacheView(similarity_image,exception);
2234#if defined(MAGICKCORE_OPENMP_SUPPORT)
2235 #pragma omp parallel for schedule(static) shared(status,similarity_info) \
2236 magick_number_threads(image,reconstruct,similarity_image->rows << 2,1)
2237#endif
2238 for (y=0; y < (ssize_t) similarity_image->rows; y++)
2239 {
2240 double
2241 similarity;
2242
2243 MagickBooleanType
2244 threshold_trigger = MagickFalse;
2245
2246 PixelPacket
2247 *magick_restrict q;
2248
2249 SimilarityInfo
2250 channel_info = similarity_info;
2251
2252 ssize_t
2253 x;
2254
2255 if (status == MagickFalse)
2256 continue;
2257 if (threshold_trigger != MagickFalse)
2258 continue;
2259 q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
2260 similarity_image->columns,1,exception);
2261 if (q == (PixelPacket *) NULL)
2262 {
2263 status=MagickFalse;
2264 continue;
2265 }
2266 for (x=0; x < (ssize_t) similarity_image->columns; x++)
2267 {
2268 MagickBooleanType
2269 update = MagickFalse;
2270
2271 similarity=GetSimilarityMetric(image,reconstruct,metric,x,y,exception);
2272 switch (metric)
2273 {
2274 case NormalizedCrossCorrelationErrorMetric:
2275 case PeakSignalToNoiseRatioMetric:
2276 {
2277 if (similarity > channel_info.similarity)
2278 update=MagickTrue;
2279 break;
2280 }
2281 default:
2282 {
2283 if (similarity < channel_info.similarity)
2284 update=MagickTrue;
2285 break;
2286 }
2287 }
2288 if (update != MagickFalse)
2289 {
2290 channel_info.similarity=similarity;
2291 channel_info.x=x;
2292 channel_info.y=y;
2293 }
2294 switch (metric)
2295 {
2296 case NormalizedCrossCorrelationErrorMetric:
2297 case PeakSignalToNoiseRatioMetric:
2298 {
2299 SetPixelRed(q,ClampToQuantum((double) QuantumRange*similarity));
2300 break;
2301 }
2302 default:
2303 {
2304 SetPixelRed(q,ClampToQuantum((double) QuantumRange*(1.0-similarity)));
2305 break;
2306 }
2307 }
2308 SetPixelGreen(q,GetPixelRed(q));
2309 SetPixelBlue(q,GetPixelRed(q));
2310 q++;
2311 }
2312#if defined(MAGICKCORE_OPENMP_SUPPORT)
2313 #pragma omp critical (MagickCore_SimilarityMetricImage)
2314#endif
2315 switch (metric)
2316 {
2317 case NormalizedCrossCorrelationErrorMetric:
2318 case PeakSignalToNoiseRatioMetric:
2319 {
2320 if (similarity_threshold != DefaultSimilarityThreshold)
2321 if (channel_info.similarity >= similarity_threshold)
2322 threshold_trigger=MagickTrue;
2323 if (channel_info.similarity >= similarity_info.similarity)
2324 similarity_info=channel_info;
2325 break;
2326 }
2327 default:
2328 {
2329 if (similarity_threshold != DefaultSimilarityThreshold)
2330 if (channel_info.similarity < similarity_threshold)
2331 threshold_trigger=MagickTrue;
2332 if (channel_info.similarity < similarity_info.similarity)
2333 similarity_info=channel_info;
2334 break;
2335 }
2336 }
2337 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2338 status=MagickFalse;
2339 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2340 {
2341 MagickBooleanType
2342 proceed;
2343
2344 progress++;
2345 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2346 if (proceed == MagickFalse)
2347 status=MagickFalse;
2348 }
2349 }
2350 similarity_view=DestroyCacheView(similarity_view);
2351 *similarity_metric=similarity_info.similarity;
2352 offset->x=similarity_info.x;
2353 offset->y=similarity_info.y;
2354 if (status == MagickFalse)
2355 similarity_image=DestroyImage(similarity_image);
2356 return(similarity_image);
2357}