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