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