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