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