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