MagickCore 7.1.1
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 %
21% dedicated 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 "MagickCore/studio.h"
44#include "MagickCore/artifact.h"
45#include "MagickCore/attribute.h"
46#include "MagickCore/cache-view.h"
47#include "MagickCore/channel.h"
48#include "MagickCore/client.h"
49#include "MagickCore/color.h"
50#include "MagickCore/color-private.h"
51#include "MagickCore/colorspace.h"
52#include "MagickCore/colorspace-private.h"
53#include "MagickCore/compare.h"
54#include "MagickCore/composite-private.h"
55#include "MagickCore/constitute.h"
56#include "MagickCore/exception-private.h"
57#include "MagickCore/enhance.h"
58#include "MagickCore/fourier.h"
59#include "MagickCore/geometry.h"
60#include "MagickCore/image-private.h"
61#include "MagickCore/list.h"
62#include "MagickCore/log.h"
63#include "MagickCore/memory_.h"
64#include "MagickCore/monitor.h"
65#include "MagickCore/monitor-private.h"
66#include "MagickCore/option.h"
67#include "MagickCore/pixel-accessor.h"
68#include "MagickCore/property.h"
69#include "MagickCore/registry.h"
70#include "MagickCore/resource_.h"
71#include "MagickCore/string_.h"
72#include "MagickCore/statistic.h"
73#include "MagickCore/statistic-private.h"
74#include "MagickCore/string-private.h"
75#include "MagickCore/thread-private.h"
76#include "MagickCore/transform.h"
77#include "MagickCore/utility.h"
78#include "MagickCore/version.h"
79␌
80/*
81%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82% %
83% %
84% %
85% C o m p a r e I m a g e s %
86% %
87% %
88% %
89%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90%
91% CompareImages() compares one or more pixel channels of an image to a
92% reconstructed image and returns the difference image.
93%
94% The format of the CompareImages method is:
95%
96% Image *CompareImages(const Image *image,const Image *reconstruct_image,
97% const MetricType metric,double *distortion,ExceptionInfo *exception)
98%
99% A description of each parameter follows:
100%
101% o image: the image.
102%
103% o reconstruct_image: the reconstruct image.
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 CacheView
117 *highlight_view,
118 *image_view,
119 *reconstruct_view;
120
121 const char
122 *artifact;
123
124 double
125 fuzz;
126
127 Image
128 *clone_image,
129 *difference_image,
130 *highlight_image;
131
132 MagickBooleanType
133 status;
134
135 PixelInfo
136 highlight,
137 lowlight,
138 masklight;
139
140 RectangleInfo
141 geometry;
142
143 size_t
144 columns,
145 rows;
146
147 ssize_t
148 y;
149
150 assert(image != (Image *) NULL);
151 assert(image->signature == MagickCoreSignature);
152 assert(reconstruct_image != (const Image *) NULL);
153 assert(reconstruct_image->signature == MagickCoreSignature);
154 assert(distortion != (double *) NULL);
155 if (IsEventLogging() != MagickFalse)
156 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
157 *distortion=0.0;
158 status=GetImageDistortion(image,reconstruct_image,metric,distortion,
159 exception);
160 if (status == MagickFalse)
161 return((Image *) NULL);
162 columns=MagickMax(image->columns,reconstruct_image->columns);
163 rows=MagickMax(image->rows,reconstruct_image->rows);
164 SetGeometry(image,&geometry);
165 geometry.width=columns;
166 geometry.height=rows;
167 clone_image=CloneImage(image,0,0,MagickTrue,exception);
168 if (clone_image == (Image *) NULL)
169 return((Image *) NULL);
170 (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
171 difference_image=ExtentImage(clone_image,&geometry,exception);
172 clone_image=DestroyImage(clone_image);
173 if (difference_image == (Image *) NULL)
174 return((Image *) NULL);
175 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
176 highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
177 if (highlight_image == (Image *) NULL)
178 {
179 difference_image=DestroyImage(difference_image);
180 return((Image *) NULL);
181 }
182 status=SetImageStorageClass(highlight_image,DirectClass,exception);
183 if (status == MagickFalse)
184 {
185 difference_image=DestroyImage(difference_image);
186 highlight_image=DestroyImage(highlight_image);
187 return((Image *) NULL);
188 }
189 (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
190 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
191 (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
192 artifact=GetImageArtifact(image,"compare:highlight-color");
193 if (artifact != (const char *) NULL)
194 (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
195 (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
196 artifact=GetImageArtifact(image,"compare:lowlight-color");
197 if (artifact != (const char *) NULL)
198 (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
199 (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
200 artifact=GetImageArtifact(image,"compare:masklight-color");
201 if (artifact != (const char *) NULL)
202 (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
203 /*
204 Generate difference image.
205 */
206 status=MagickTrue;
207 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
208 image_view=AcquireVirtualCacheView(image,exception);
209 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
210 highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
211#if defined(MAGICKCORE_OPENMP_SUPPORT)
212 #pragma omp parallel for schedule(static) shared(status) \
213 magick_number_threads(image,highlight_image,rows,1)
214#endif
215 for (y=0; y < (ssize_t) rows; y++)
216 {
217 MagickBooleanType
218 sync;
219
220 const Quantum
221 *magick_restrict p,
222 *magick_restrict q;
223
224 Quantum
225 *magick_restrict r;
226
227 ssize_t
228 x;
229
230 if (status == MagickFalse)
231 continue;
232 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
233 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
234 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
235 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
236 (r == (Quantum *) NULL))
237 {
238 status=MagickFalse;
239 continue;
240 }
241 for (x=0; x < (ssize_t) columns; x++)
242 {
243 double
244 Da,
245 Sa;
246
247 MagickStatusType
248 difference;
249
250 ssize_t
251 i;
252
253 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
254 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
255 {
256 SetPixelViaPixelInfo(highlight_image,&masklight,r);
257 p+=(ptrdiff_t) GetPixelChannels(image);
258 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
259 r+=(ptrdiff_t) GetPixelChannels(highlight_image);
260 continue;
261 }
262 difference=MagickFalse;
263 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
264 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
265 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
266 {
267 double
268 distance,
269 pixel;
270
271 PixelChannel channel = GetPixelChannelChannel(image,i);
272 PixelTrait traits = GetPixelChannelTraits(image,channel);
273 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
274 channel);
275 if ((traits == UndefinedPixelTrait) ||
276 (reconstruct_traits == UndefinedPixelTrait) ||
277 ((reconstruct_traits & UpdatePixelTrait) == 0))
278 continue;
279 if (channel == AlphaPixelChannel)
280 pixel=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
281 channel,q);
282 else
283 pixel=Sa*(double) p[i]-Da*(double)
284 GetPixelChannel(reconstruct_image,channel,q);
285 distance=pixel*pixel;
286 if (distance >= fuzz)
287 {
288 difference=MagickTrue;
289 break;
290 }
291 }
292 if (difference == MagickFalse)
293 SetPixelViaPixelInfo(highlight_image,&lowlight,r);
294 else
295 SetPixelViaPixelInfo(highlight_image,&highlight,r);
296 p+=(ptrdiff_t) GetPixelChannels(image);
297 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
298 r+=(ptrdiff_t) GetPixelChannels(highlight_image);
299 }
300 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
301 if (sync == MagickFalse)
302 status=MagickFalse;
303 }
304 highlight_view=DestroyCacheView(highlight_view);
305 reconstruct_view=DestroyCacheView(reconstruct_view);
306 image_view=DestroyCacheView(image_view);
307 (void) CompositeImage(difference_image,highlight_image,image->compose,
308 MagickTrue,0,0,exception);
309 highlight_image=DestroyImage(highlight_image);
310 if (status == MagickFalse)
311 difference_image=DestroyImage(difference_image);
312 return(difference_image);
313}
314␌
315/*
316%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
317% %
318% %
319% %
320% G e t I m a g e D i s t o r t i o n %
321% %
322% %
323% %
324%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
325%
326% GetImageDistortion() compares one or more pixel channels of an image to a
327% reconstructed image and returns the specified distortion metric.
328%
329% The format of the GetImageDistortion method is:
330%
331% MagickBooleanType GetImageDistortion(const Image *image,
332% const Image *reconstruct_image,const MetricType metric,
333% double *distortion,ExceptionInfo *exception)
334%
335% A description of each parameter follows:
336%
337% o image: the image.
338%
339% o reconstruct_image: the reconstruct image.
340%
341% o metric: the metric.
342%
343% o distortion: the computed distortion between the images.
344%
345% o exception: return any errors or warnings in this structure.
346%
347*/
348
349static MagickBooleanType GetAbsoluteDistortion(const Image *image,
350 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
351{
352 CacheView
353 *image_view,
354 *reconstruct_view;
355
356 double
357 fuzz;
358
359 MagickBooleanType
360 status;
361
362 size_t
363 columns,
364 rows;
365
366 ssize_t
367 y;
368
369 /*
370 Compute the absolute difference in pixels between two images.
371 */
372 status=MagickTrue;
373 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
374 rows=MagickMax(image->rows,reconstruct_image->rows);
375 columns=MagickMax(image->columns,reconstruct_image->columns);
376 image_view=AcquireVirtualCacheView(image,exception);
377 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
378#if defined(MAGICKCORE_OPENMP_SUPPORT)
379 #pragma omp parallel for schedule(static) shared(status) \
380 magick_number_threads(image,image,rows,1)
381#endif
382 for (y=0; y < (ssize_t) rows; y++)
383 {
384 double
385 channel_distortion[MaxPixelChannels+1];
386
387 const Quantum
388 *magick_restrict p,
389 *magick_restrict q;
390
391 ssize_t
392 j,
393 x;
394
395 if (status == MagickFalse)
396 continue;
397 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
398 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
399 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
400 {
401 status=MagickFalse;
402 continue;
403 }
404 (void) memset(channel_distortion,0,sizeof(channel_distortion));
405 for (x=0; x < (ssize_t) columns; x++)
406 {
407 double
408 Da,
409 Sa;
410
411 MagickBooleanType
412 difference;
413
414 ssize_t
415 i;
416
417 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
418 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
419 {
420 p+=(ptrdiff_t) GetPixelChannels(image);
421 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
422 continue;
423 }
424 difference=MagickFalse;
425 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
426 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
427 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
428 {
429 double
430 distance,
431 pixel;
432
433 PixelChannel channel = GetPixelChannelChannel(image,i);
434 PixelTrait traits = GetPixelChannelTraits(image,channel);
435 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
436 channel);
437 if ((traits == UndefinedPixelTrait) ||
438 (reconstruct_traits == UndefinedPixelTrait) ||
439 ((reconstruct_traits & UpdatePixelTrait) == 0))
440 continue;
441 if (channel == AlphaPixelChannel)
442 pixel=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
443 channel,q);
444 else
445 pixel=Sa*(double) p[i]-Da*(double) GetPixelChannel(reconstruct_image,
446 channel,q);
447 distance=pixel*pixel;
448 if (distance >= fuzz)
449 {
450 channel_distortion[i]++;
451 difference=MagickTrue;
452 }
453 }
454 if (difference != MagickFalse)
455 channel_distortion[CompositePixelChannel]++;
456 p+=(ptrdiff_t) GetPixelChannels(image);
457 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
458 }
459#if defined(MAGICKCORE_OPENMP_SUPPORT)
460 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
461#endif
462 for (j=0; j <= MaxPixelChannels; j++)
463 distortion[j]+=channel_distortion[j];
464 }
465 reconstruct_view=DestroyCacheView(reconstruct_view);
466 image_view=DestroyCacheView(image_view);
467 return(status);
468}
469
470static MagickBooleanType GetFuzzDistortion(const Image *image,
471 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
472{
473 CacheView
474 *image_view,
475 *reconstruct_view;
476
477 double
478 area;
479
480 MagickBooleanType
481 status;
482
483 ssize_t
484 j;
485
486 size_t
487 columns,
488 rows;
489
490 ssize_t
491 y;
492
493 status=MagickTrue;
494 rows=MagickMax(image->rows,reconstruct_image->rows);
495 columns=MagickMax(image->columns,reconstruct_image->columns);
496 area=0.0;
497 image_view=AcquireVirtualCacheView(image,exception);
498 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
499#if defined(MAGICKCORE_OPENMP_SUPPORT)
500 #pragma omp parallel for schedule(static) shared(status) \
501 magick_number_threads(image,image,rows,1)
502#endif
503 for (y=0; y < (ssize_t) rows; y++)
504 {
505 double
506 channel_distortion[MaxPixelChannels+1];
507
508 const Quantum
509 *magick_restrict p,
510 *magick_restrict q;
511
512 size_t
513 local_area = 0;
514
515 ssize_t
516 x;
517
518 if (status == MagickFalse)
519 continue;
520 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
521 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
522 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
523 {
524 status=MagickFalse;
525 continue;
526 }
527 (void) memset(channel_distortion,0,sizeof(channel_distortion));
528 for (x=0; x < (ssize_t) columns; x++)
529 {
530 double
531 Da,
532 Sa;
533
534 ssize_t
535 i;
536
537 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
538 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
539 {
540 p+=(ptrdiff_t) GetPixelChannels(image);
541 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
542 continue;
543 }
544 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
545 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
546 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
547 {
548 double
549 distance;
550
551 PixelChannel channel = GetPixelChannelChannel(image,i);
552 PixelTrait traits = GetPixelChannelTraits(image,channel);
553 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
554 channel);
555 if ((traits == UndefinedPixelTrait) ||
556 (reconstruct_traits == UndefinedPixelTrait) ||
557 ((reconstruct_traits & UpdatePixelTrait) == 0))
558 continue;
559 if (channel == AlphaPixelChannel)
560 distance=QuantumScale*((double) p[i]-(double) GetPixelChannel(
561 reconstruct_image,channel,q));
562 else
563 distance=QuantumScale*(Sa*(double) p[i]-Da*(double) GetPixelChannel(
564 reconstruct_image,channel,q));
565 channel_distortion[i]+=distance*distance;
566 channel_distortion[CompositePixelChannel]+=distance*distance;
567 }
568 local_area++;
569 p+=(ptrdiff_t) GetPixelChannels(image);
570 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
571 }
572#if defined(MAGICKCORE_OPENMP_SUPPORT)
573 #pragma omp critical (MagickCore_GetFuzzDistortion)
574#endif
575 {
576 area+=local_area;
577 for (j=0; j <= MaxPixelChannels; j++)
578 distortion[j]+=channel_distortion[j];
579 }
580 }
581 reconstruct_view=DestroyCacheView(reconstruct_view);
582 image_view=DestroyCacheView(image_view);
583 area=PerceptibleReciprocal(area);
584 for (j=0; j <= MaxPixelChannels; j++)
585 distortion[j]*=area;
586 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
587 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
588 return(status);
589}
590
591static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
592 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
593{
594 CacheView
595 *image_view,
596 *reconstruct_view;
597
598 double
599 area;
600
601 MagickBooleanType
602 status;
603
604 ssize_t
605 j;
606
607 size_t
608 columns,
609 rows;
610
611 ssize_t
612 y;
613
614 status=MagickTrue;
615 rows=MagickMax(image->rows,reconstruct_image->rows);
616 columns=MagickMax(image->columns,reconstruct_image->columns);
617 area=0.0;
618 image_view=AcquireVirtualCacheView(image,exception);
619 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
620#if defined(MAGICKCORE_OPENMP_SUPPORT)
621 #pragma omp parallel for schedule(static) shared(status) \
622 magick_number_threads(image,image,rows,1)
623#endif
624 for (y=0; y < (ssize_t) rows; y++)
625 {
626 double
627 channel_distortion[MaxPixelChannels+1];
628
629 const Quantum
630 *magick_restrict p,
631 *magick_restrict q;
632
633 size_t
634 local_area = 0;
635
636 ssize_t
637 x;
638
639 if (status == MagickFalse)
640 continue;
641 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
642 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
643 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
644 {
645 status=MagickFalse;
646 continue;
647 }
648 (void) memset(channel_distortion,0,sizeof(channel_distortion));
649 for (x=0; x < (ssize_t) columns; x++)
650 {
651 double
652 Da,
653 Sa;
654
655 ssize_t
656 i;
657
658 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
659 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
660 {
661 p+=(ptrdiff_t) GetPixelChannels(image);
662 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
663 continue;
664 }
665 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
666 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
667 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
668 {
669 double
670 distance;
671
672 PixelChannel channel = GetPixelChannelChannel(image,i);
673 PixelTrait traits = GetPixelChannelTraits(image,channel);
674 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
675 channel);
676 if ((traits == UndefinedPixelTrait) ||
677 (reconstruct_traits == UndefinedPixelTrait) ||
678 ((reconstruct_traits & UpdatePixelTrait) == 0))
679 continue;
680 if (channel == AlphaPixelChannel)
681 distance=QuantumScale*fabs((double) p[i]-(double)
682 GetPixelChannel(reconstruct_image,channel,q));
683 else
684 distance=QuantumScale*fabs(Sa*(double) p[i]-Da*(double)
685 GetPixelChannel(reconstruct_image,channel,q));
686 channel_distortion[i]+=distance;
687 channel_distortion[CompositePixelChannel]+=distance;
688 }
689 local_area++;
690 p+=(ptrdiff_t) GetPixelChannels(image);
691 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
692 }
693#if defined(MAGICKCORE_OPENMP_SUPPORT)
694 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
695#endif
696 {
697 area+=local_area;
698 for (j=0; j <= MaxPixelChannels; j++)
699 distortion[j]+=channel_distortion[j];
700 }
701 }
702 reconstruct_view=DestroyCacheView(reconstruct_view);
703 image_view=DestroyCacheView(image_view);
704 area=PerceptibleReciprocal(area);
705 for (j=0; j <= MaxPixelChannels; j++)
706 distortion[j]*=area;
707 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
708 return(status);
709}
710
711static MagickBooleanType GetMeanErrorPerPixel(Image *image,
712 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
713{
714 CacheView
715 *image_view,
716 *reconstruct_view;
717
718 MagickBooleanType
719 status;
720
721 double
722 area,
723 maximum_error,
724 mean_error;
725
726 size_t
727 columns,
728 rows;
729
730 ssize_t
731 y;
732
733 status=MagickTrue;
734 area=0.0;
735 maximum_error=0.0;
736 mean_error=0.0;
737 rows=MagickMax(image->rows,reconstruct_image->rows);
738 columns=MagickMax(image->columns,reconstruct_image->columns);
739 image_view=AcquireVirtualCacheView(image,exception);
740 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
741 for (y=0; y < (ssize_t) rows; y++)
742 {
743 const Quantum
744 *magick_restrict p,
745 *magick_restrict q;
746
747 ssize_t
748 x;
749
750 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
751 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
752 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
753 {
754 status=MagickFalse;
755 break;
756 }
757 for (x=0; x < (ssize_t) columns; x++)
758 {
759 double
760 Da,
761 Sa;
762
763 ssize_t
764 i;
765
766 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
767 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
768 {
769 p+=(ptrdiff_t) GetPixelChannels(image);
770 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
771 continue;
772 }
773 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
774 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
775 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
776 {
777 double
778 distance;
779
780 PixelChannel channel = GetPixelChannelChannel(image,i);
781 PixelTrait traits = GetPixelChannelTraits(image,channel);
782 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
783 channel);
784 if ((traits == UndefinedPixelTrait) ||
785 (reconstruct_traits == UndefinedPixelTrait) ||
786 ((reconstruct_traits & UpdatePixelTrait) == 0))
787 continue;
788 if (channel == AlphaPixelChannel)
789 distance=fabs((double) p[i]-(double)
790 GetPixelChannel(reconstruct_image,channel,q));
791 else
792 distance=fabs(Sa*(double) p[i]-Da*(double)
793 GetPixelChannel(reconstruct_image,channel,q));
794 distortion[i]+=distance;
795 distortion[CompositePixelChannel]+=distance;
796 mean_error+=distance*distance;
797 if (distance > maximum_error)
798 maximum_error=distance;
799 area++;
800 }
801 p+=(ptrdiff_t) GetPixelChannels(image);
802 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
803 }
804 }
805 reconstruct_view=DestroyCacheView(reconstruct_view);
806 image_view=DestroyCacheView(image_view);
807 area=PerceptibleReciprocal(area);
808 image->error.mean_error_per_pixel=area*distortion[CompositePixelChannel];
809 image->error.normalized_mean_error=area*QuantumScale*QuantumScale*mean_error;
810 image->error.normalized_maximum_error=QuantumScale*maximum_error;
811 return(status);
812}
813
814static MagickBooleanType GetMeanSquaredDistortion(const Image *image,
815 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
816{
817 CacheView
818 *image_view,
819 *reconstruct_view;
820
821 double
822 area;
823
824 MagickBooleanType
825 status;
826
827 size_t
828 columns,
829 rows;
830
831 ssize_t
832 j,
833 y;
834
835 status=MagickTrue;
836 rows=MagickMax(image->rows,reconstruct_image->rows);
837 columns=MagickMax(image->columns,reconstruct_image->columns);
838 area=0.0;
839 image_view=AcquireVirtualCacheView(image,exception);
840 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
841#if defined(MAGICKCORE_OPENMP_SUPPORT)
842 #pragma omp parallel for schedule(static) shared(status) \
843 magick_number_threads(image,image,rows,1)
844#endif
845 for (y=0; y < (ssize_t) rows; y++)
846 {
847 const Quantum
848 *magick_restrict p,
849 *magick_restrict q;
850
851 double
852 channel_distortion[MaxPixelChannels+1];
853
854 size_t
855 local_area = 0;
856
857 ssize_t
858 x;
859
860 if (status == MagickFalse)
861 continue;
862 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
863 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
864 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
865 {
866 status=MagickFalse;
867 continue;
868 }
869 (void) memset(channel_distortion,0,sizeof(channel_distortion));
870 for (x=0; x < (ssize_t) columns; x++)
871 {
872 double
873 Da,
874 Sa;
875
876 ssize_t
877 i;
878
879 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
880 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
881 {
882 p+=(ptrdiff_t) GetPixelChannels(image);
883 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
884 continue;
885 }
886 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
887 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
888 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
889 {
890 double
891 distance;
892
893 PixelChannel channel = GetPixelChannelChannel(image,i);
894 PixelTrait traits = GetPixelChannelTraits(image,channel);
895 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
896 channel);
897 if ((traits == UndefinedPixelTrait) ||
898 (reconstruct_traits == UndefinedPixelTrait) ||
899 ((reconstruct_traits & UpdatePixelTrait) == 0))
900 continue;
901 if (channel == AlphaPixelChannel)
902 distance=QuantumScale*((double) p[i]-(double) GetPixelChannel(
903 reconstruct_image,channel,q));
904 else
905 distance=QuantumScale*(Sa*(double) p[i]-Da*(double) GetPixelChannel(
906 reconstruct_image,channel,q));
907 channel_distortion[i]+=distance*distance;
908 channel_distortion[CompositePixelChannel]+=distance*distance;
909 }
910 local_area++;
911 p+=(ptrdiff_t) GetPixelChannels(image);
912 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
913 }
914#if defined(MAGICKCORE_OPENMP_SUPPORT)
915 #pragma omp critical (MagickCore_GetMeanSquaredError)
916#endif
917 {
918 area+=local_area;
919 for (j=0; j <= MaxPixelChannels; j++)
920 distortion[j]+=channel_distortion[j];
921 }
922 }
923 reconstruct_view=DestroyCacheView(reconstruct_view);
924 image_view=DestroyCacheView(image_view);
925 area=PerceptibleReciprocal(area);
926 for (j=0; j <= MaxPixelChannels; j++)
927 distortion[j]*=area;
928 distortion[CompositePixelChannel]/=GetImageChannels(image);
929 return(status);
930}
931
932static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
933 const Image *image,const Image *reconstruct_image,double *distortion,
934 ExceptionInfo *exception)
935{
936#define SimilarityImageTag "Similarity/Image"
937
938 CacheView
939 *image_view,
940 *reconstruct_view;
941
942 ChannelStatistics
943 *image_statistics,
944 *reconstruct_statistics;
945
946 double
947 alpha_variance[MaxPixelChannels+1],
948 beta_variance[MaxPixelChannels+1];
949
950 MagickBooleanType
951 status;
952
953 MagickOffsetType
954 progress;
955
956 size_t
957 columns,
958 rows;
959
960 ssize_t
961 i,
962 y;
963
964 /*
965 Normalized cross correlation distortion.
966 */
967 image_statistics=GetImageStatistics(image,exception);
968 reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
969 if ((image_statistics == (ChannelStatistics *) NULL) ||
970 (reconstruct_statistics == (ChannelStatistics *) NULL))
971 {
972 if (image_statistics != (ChannelStatistics *) NULL)
973 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
974 image_statistics);
975 if (reconstruct_statistics != (ChannelStatistics *) NULL)
976 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
977 reconstruct_statistics);
978 return(MagickFalse);
979 }
980 (void) memset(distortion,0,(MaxPixelChannels+1)*sizeof(*distortion));
981 (void) memset(alpha_variance,0,(MaxPixelChannels+1)*sizeof(*alpha_variance));
982 (void) memset(beta_variance,0,(MaxPixelChannels+1)*sizeof(*beta_variance));
983 status=MagickTrue;
984 progress=0;
985 columns=MagickMax(image->columns,reconstruct_image->columns);
986 rows=MagickMax(image->rows,reconstruct_image->rows);
987 image_view=AcquireVirtualCacheView(image,exception);
988 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
989 for (y=0; y < (ssize_t) rows; y++)
990 {
991 const Quantum
992 *magick_restrict p,
993 *magick_restrict q;
994
995 ssize_t
996 x;
997
998 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
999 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1000 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1001 {
1002 status=MagickFalse;
1003 break;
1004 }
1005 for (x=0; x < (ssize_t) columns; x++)
1006 {
1007 double
1008 Da,
1009 Sa;
1010
1011 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1012 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1013 {
1014 p+=(ptrdiff_t) GetPixelChannels(image);
1015 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1016 continue;
1017 }
1018 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1019 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1020 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1021 {
1022 double
1023 alpha,
1024 beta;
1025
1026 PixelChannel channel = GetPixelChannelChannel(image,i);
1027 PixelTrait traits = GetPixelChannelTraits(image,channel);
1028 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1029 channel);
1030 if ((traits == UndefinedPixelTrait) ||
1031 (reconstruct_traits == UndefinedPixelTrait) ||
1032 ((reconstruct_traits & UpdatePixelTrait) == 0))
1033 continue;
1034 if (channel == AlphaPixelChannel)
1035 {
1036 alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
1037 beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
1038 channel,q)-reconstruct_statistics[channel].mean);
1039 }
1040 else
1041 {
1042 alpha=QuantumScale*(Sa*(double) p[i]-
1043 image_statistics[channel].mean);
1044 beta=QuantumScale*(Da*(double) GetPixelChannel(reconstruct_image,
1045 channel,q)-reconstruct_statistics[channel].mean);
1046 }
1047 distortion[i]+=alpha*beta;
1048 alpha_variance[i]+=alpha*alpha;
1049 beta_variance[i]+=beta*beta;
1050 }
1051 p+=(ptrdiff_t) GetPixelChannels(image);
1052 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1053 }
1054 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1055 {
1056 MagickBooleanType
1057 proceed;
1058
1059#if defined(MAGICKCORE_OPENMP_SUPPORT)
1060 #pragma omp atomic
1061#endif
1062 progress++;
1063 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1064 if (proceed == MagickFalse)
1065 {
1066 status=MagickFalse;
1067 break;
1068 }
1069 }
1070 }
1071 reconstruct_view=DestroyCacheView(reconstruct_view);
1072 image_view=DestroyCacheView(image_view);
1073 /*
1074 Compute normalized cross correlation: divide by standard deviation.
1075 */
1076 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1077 {
1078 distortion[i]/=sqrt(alpha_variance[i]*beta_variance[i]);
1079 if (fabs(distortion[i]) > MagickEpsilon)
1080 distortion[CompositePixelChannel]+=distortion[i];
1081 }
1082 distortion[CompositePixelChannel]=distortion[CompositePixelChannel]/
1083 GetImageChannels(image);
1084 /*
1085 Free resources.
1086 */
1087 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1088 reconstruct_statistics);
1089 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1090 image_statistics);
1091 return(status);
1092}
1093
1094static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1095 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1096{
1097 CacheView
1098 *image_view,
1099 *reconstruct_view;
1100
1101 MagickBooleanType
1102 status;
1103
1104 size_t
1105 columns,
1106 rows;
1107
1108 ssize_t
1109 y;
1110
1111 status=MagickTrue;
1112 rows=MagickMax(image->rows,reconstruct_image->rows);
1113 columns=MagickMax(image->columns,reconstruct_image->columns);
1114 image_view=AcquireVirtualCacheView(image,exception);
1115 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1116#if defined(MAGICKCORE_OPENMP_SUPPORT)
1117 #pragma omp parallel for schedule(static) shared(status) \
1118 magick_number_threads(image,image,rows,1)
1119#endif
1120 for (y=0; y < (ssize_t) rows; y++)
1121 {
1122 double
1123 channel_distortion[MaxPixelChannels+1];
1124
1125 const Quantum
1126 *magick_restrict p,
1127 *magick_restrict q;
1128
1129 ssize_t
1130 j,
1131 x;
1132
1133 if (status == MagickFalse)
1134 continue;
1135 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1136 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1137 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1138 {
1139 status=MagickFalse;
1140 continue;
1141 }
1142 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1143 for (x=0; x < (ssize_t) columns; x++)
1144 {
1145 double
1146 Da,
1147 Sa;
1148
1149 ssize_t
1150 i;
1151
1152 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1153 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1154 {
1155 p+=(ptrdiff_t) GetPixelChannels(image);
1156 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1157 continue;
1158 }
1159 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1160 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1161 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1162 {
1163 double
1164 distance;
1165
1166 PixelChannel channel = GetPixelChannelChannel(image,i);
1167 PixelTrait traits = GetPixelChannelTraits(image,channel);
1168 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1169 channel);
1170 if ((traits == UndefinedPixelTrait) ||
1171 (reconstruct_traits == UndefinedPixelTrait) ||
1172 ((reconstruct_traits & UpdatePixelTrait) == 0))
1173 continue;
1174 if (channel == AlphaPixelChannel)
1175 distance=QuantumScale*fabs((double) p[i]-(double)
1176 GetPixelChannel(reconstruct_image,channel,q));
1177 else
1178 distance=QuantumScale*fabs(Sa*(double) p[i]-Da*(double)
1179 GetPixelChannel(reconstruct_image,channel,q));
1180 if (distance > channel_distortion[i])
1181 channel_distortion[i]=distance;
1182 if (distance > channel_distortion[CompositePixelChannel])
1183 channel_distortion[CompositePixelChannel]=distance;
1184 }
1185 p+=(ptrdiff_t) GetPixelChannels(image);
1186 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1187 }
1188#if defined(MAGICKCORE_OPENMP_SUPPORT)
1189 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1190#endif
1191 for (j=0; j <= MaxPixelChannels; j++)
1192 if (channel_distortion[j] > distortion[j])
1193 distortion[j]=channel_distortion[j];
1194 }
1195 reconstruct_view=DestroyCacheView(reconstruct_view);
1196 image_view=DestroyCacheView(image_view);
1197 return(status);
1198}
1199
1200static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1201 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1202{
1203 MagickBooleanType
1204 status;
1205
1206 ssize_t
1207 i;
1208
1209 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1210 for (i=0; i <= MaxPixelChannels; i++)
1211 if (fabs(distortion[i]) >= MagickEpsilon)
1212 distortion[i]=(-10.0*MagickLog10(distortion[i]));
1213 return(status);
1214}
1215
1216static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1217 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1218{
1219 ChannelPerceptualHash
1220 *channel_phash,
1221 *reconstruct_phash;
1222
1223 const char
1224 *artifact;
1225
1226 ssize_t
1227 channel;
1228
1229 /*
1230 Compute perceptual hash in the sRGB colorspace.
1231 */
1232 channel_phash=GetImagePerceptualHash(image,exception);
1233 if (channel_phash == (ChannelPerceptualHash *) NULL)
1234 return(MagickFalse);
1235 reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1236 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1237 {
1238 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1239 channel_phash);
1240 return(MagickFalse);
1241 }
1242#if defined(MAGICKCORE_OPENMP_SUPPORT)
1243 #pragma omp parallel for schedule(static)
1244#endif
1245 for (channel=0; channel < MaxPixelChannels; channel++)
1246 {
1247 double
1248 difference;
1249
1250 ssize_t
1251 i;
1252
1253 difference=0.0;
1254 for (i=0; i < MaximumNumberOfImageMoments; i++)
1255 {
1256 double
1257 alpha,
1258 beta;
1259
1260 ssize_t
1261 j;
1262
1263 for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1264 {
1265 alpha=channel_phash[channel].phash[j][i];
1266 beta=reconstruct_phash[channel].phash[j][i];
1267 difference+=(beta-alpha)*(beta-alpha);
1268 }
1269 }
1270 distortion[channel]+=difference;
1271#if defined(MAGICKCORE_OPENMP_SUPPORT)
1272 #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1273#endif
1274 distortion[CompositePixelChannel]+=difference;
1275 }
1276 artifact=GetImageArtifact(image,"phash:normalize");
1277 if ((artifact != (const char *) NULL) &&
1278 (IsStringTrue(artifact) != MagickFalse))
1279 {
1280 ssize_t
1281 j;
1282
1283 for (j=0; j <= MaxPixelChannels; j++)
1284 distortion[j]=sqrt(distortion[j]/channel_phash[0].number_channels);
1285 }
1286 /*
1287 Free resources.
1288 */
1289 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1290 reconstruct_phash);
1291 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1292 return(MagickTrue);
1293}
1294
1295static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1296 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1297{
1298 MagickBooleanType
1299 status;
1300
1301 ssize_t
1302 i;
1303
1304 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1305 for (i=0; i <= MaxPixelChannels; i++)
1306 distortion[i]=sqrt(distortion[i]);
1307 return(status);
1308}
1309
1310static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image,
1311 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1312{
1313#define SSIMRadius 5.0
1314#define SSIMSigma 1.5
1315#define SSIMBlocksize 8
1316#define SSIMK1 0.01
1317#define SSIMK2 0.03
1318#define SSIML 1.0
1319
1320 CacheView
1321 *image_view,
1322 *reconstruct_view;
1323
1324 char
1325 geometry[MagickPathExtent];
1326
1327 const char
1328 *artifact;
1329
1330 double
1331 area,
1332 c1,
1333 c2,
1334 radius,
1335 sigma;
1336
1337 KernelInfo
1338 *kernel_info;
1339
1340 MagickBooleanType
1341 status;
1342
1343 ssize_t
1344 j;
1345
1346 size_t
1347 columns,
1348 rows;
1349
1350 ssize_t
1351 y;
1352
1353 /*
1354 Compute structural similarity index @
1355 https://en.wikipedia.org/wiki/Structural_similarity.
1356 */
1357 radius=SSIMRadius;
1358 artifact=GetImageArtifact(image,"compare:ssim-radius");
1359 if (artifact != (const char *) NULL)
1360 radius=StringToDouble(artifact,(char **) NULL);
1361 sigma=SSIMSigma;
1362 artifact=GetImageArtifact(image,"compare:ssim-sigma");
1363 if (artifact != (const char *) NULL)
1364 sigma=StringToDouble(artifact,(char **) NULL);
1365 (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1366 radius,sigma);
1367 kernel_info=AcquireKernelInfo(geometry,exception);
1368 if (kernel_info == (KernelInfo *) NULL)
1369 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1370 image->filename);
1371 c1=pow(SSIMK1*SSIML,2.0);
1372 artifact=GetImageArtifact(image,"compare:ssim-k1");
1373 if (artifact != (const char *) NULL)
1374 c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1375 c2=pow(SSIMK2*SSIML,2.0);
1376 artifact=GetImageArtifact(image,"compare:ssim-k2");
1377 if (artifact != (const char *) NULL)
1378 c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1379 status=MagickTrue;
1380 area=0.0;
1381 rows=MagickMax(image->rows,reconstruct_image->rows);
1382 columns=MagickMax(image->columns,reconstruct_image->columns);
1383 image_view=AcquireVirtualCacheView(image,exception);
1384 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1385#if defined(MAGICKCORE_OPENMP_SUPPORT)
1386 #pragma omp parallel for schedule(static) shared(status) \
1387 magick_number_threads(image,reconstruct_image,rows,1)
1388#endif
1389 for (y=0; y < (ssize_t) rows; y++)
1390 {
1391 double
1392 channel_distortion[MaxPixelChannels+1];
1393
1394 const Quantum
1395 *magick_restrict p,
1396 *magick_restrict q;
1397
1398 size_t
1399 local_area = 0;
1400
1401 ssize_t
1402 i,
1403 x;
1404
1405 if (status == MagickFalse)
1406 continue;
1407 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1408 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1409 kernel_info->height,exception);
1410 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1411 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1412 kernel_info->height,exception);
1413 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1414 {
1415 status=MagickFalse;
1416 continue;
1417 }
1418 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1419 for (x=0; x < (ssize_t) columns; x++)
1420 {
1421 double
1422 x_pixel_mu[MaxPixelChannels+1],
1423 x_pixel_sigma_squared[MaxPixelChannels+1],
1424 xy_sigma[MaxPixelChannels+1],
1425 y_pixel_mu[MaxPixelChannels+1],
1426 y_pixel_sigma_squared[MaxPixelChannels+1];
1427
1428 const Quantum
1429 *magick_restrict reference,
1430 *magick_restrict target;
1431
1432 MagickRealType
1433 *k;
1434
1435 ssize_t
1436 v;
1437
1438 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1439 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1440 {
1441 p+=(ptrdiff_t) GetPixelChannels(image);
1442 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1443 continue;
1444 }
1445 (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1446 (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1447 (void) memset(xy_sigma,0,sizeof(xy_sigma));
1448 (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1449 (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1450 (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1451 k=kernel_info->values;
1452 reference=p;
1453 target=q;
1454 for (v=0; v < (ssize_t) kernel_info->height; v++)
1455 {
1456 ssize_t
1457 u;
1458
1459 for (u=0; u < (ssize_t) kernel_info->width; u++)
1460 {
1461 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1462 {
1463 double
1464 x_pixel,
1465 y_pixel;
1466
1467 PixelChannel channel = GetPixelChannelChannel(image,i);
1468 PixelTrait traits = GetPixelChannelTraits(image,channel);
1469 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1470 reconstruct_image,channel);
1471 if ((traits == UndefinedPixelTrait) ||
1472 (reconstruct_traits == UndefinedPixelTrait) ||
1473 ((reconstruct_traits & UpdatePixelTrait) == 0))
1474 continue;
1475 x_pixel=QuantumScale*(double) reference[i];
1476 x_pixel_mu[i]+=(*k)*x_pixel;
1477 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1478 y_pixel=QuantumScale*(double)
1479 GetPixelChannel(reconstruct_image,channel,target);
1480 y_pixel_mu[i]+=(*k)*y_pixel;
1481 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1482 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1483 }
1484 k++;
1485 reference+=(ptrdiff_t) GetPixelChannels(image);
1486 target+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1487 }
1488 reference+=(ptrdiff_t) GetPixelChannels(image)*columns;
1489 target+=(ptrdiff_t) GetPixelChannels(reconstruct_image)*columns;
1490 }
1491 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1492 {
1493 double
1494 ssim,
1495 x_pixel_mu_squared,
1496 x_pixel_sigmas_squared,
1497 xy_mu,
1498 xy_sigmas,
1499 y_pixel_mu_squared,
1500 y_pixel_sigmas_squared;
1501
1502 PixelChannel channel = GetPixelChannelChannel(image,i);
1503 PixelTrait traits = GetPixelChannelTraits(image,channel);
1504 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1505 reconstruct_image,channel);
1506 if ((traits == UndefinedPixelTrait) ||
1507 (reconstruct_traits == UndefinedPixelTrait) ||
1508 ((reconstruct_traits & UpdatePixelTrait) == 0))
1509 continue;
1510 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1511 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1512 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1513 xy_sigmas=xy_sigma[i]-xy_mu;
1514 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1515 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1516 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1517 ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1518 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1519 channel_distortion[i]+=ssim;
1520 channel_distortion[CompositePixelChannel]+=ssim;
1521 }
1522 p+=(ptrdiff_t) GetPixelChannels(image);
1523 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1524 local_area++;
1525 }
1526#if defined(MAGICKCORE_OPENMP_SUPPORT)
1527 #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1528#endif
1529 {
1530 area+=local_area;
1531 for (i=0; i <= MaxPixelChannels; i++)
1532 distortion[i]+=channel_distortion[i];
1533 }
1534 }
1535 image_view=DestroyCacheView(image_view);
1536 reconstruct_view=DestroyCacheView(reconstruct_view);
1537 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1538 {
1539 PixelChannel channel = GetPixelChannelChannel(image,j);
1540 PixelTrait traits = GetPixelChannelTraits(image,channel);
1541 if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1542 continue;
1543 distortion[j]/=area;
1544 }
1545 distortion[CompositePixelChannel]/=area;
1546 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1547 kernel_info=DestroyKernelInfo(kernel_info);
1548 return(status);
1549}
1550
1551static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image,
1552 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1553{
1554 MagickBooleanType
1555 status;
1556
1557 ssize_t
1558 i;
1559
1560 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1561 distortion,exception);
1562 for (i=0; i <= MaxPixelChannels; i++)
1563 distortion[i]=(1.0-(distortion[i]))/2.0;
1564 return(status);
1565}
1566
1567MagickExport MagickBooleanType GetImageDistortion(Image *image,
1568 const Image *reconstruct_image,const MetricType metric,double *distortion,
1569 ExceptionInfo *exception)
1570{
1571 double
1572 *channel_distortion;
1573
1574 MagickBooleanType
1575 status;
1576
1577 size_t
1578 length;
1579
1580 assert(image != (Image *) NULL);
1581 assert(image->signature == MagickCoreSignature);
1582 assert(reconstruct_image != (const Image *) NULL);
1583 assert(reconstruct_image->signature == MagickCoreSignature);
1584 assert(distortion != (double *) NULL);
1585 if (IsEventLogging() != MagickFalse)
1586 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1587 *distortion=0.0;
1588 /*
1589 Get image distortion.
1590 */
1591 length=MaxPixelChannels+1UL;
1592 channel_distortion=(double *) AcquireQuantumMemory(length,
1593 sizeof(*channel_distortion));
1594 if (channel_distortion == (double *) NULL)
1595 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1596 (void) memset(channel_distortion,0,length*
1597 sizeof(*channel_distortion));
1598 switch (metric)
1599 {
1600 case AbsoluteErrorMetric:
1601 {
1602 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1603 exception);
1604 break;
1605 }
1606 case FuzzErrorMetric:
1607 {
1608 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1609 exception);
1610 break;
1611 }
1612 case MeanAbsoluteErrorMetric:
1613 {
1614 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1615 channel_distortion,exception);
1616 break;
1617 }
1618 case MeanErrorPerPixelErrorMetric:
1619 {
1620 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1621 exception);
1622 break;
1623 }
1624 case MeanSquaredErrorMetric:
1625 {
1626 status=GetMeanSquaredDistortion(image,reconstruct_image,
1627 channel_distortion,exception);
1628 break;
1629 }
1630 case NormalizedCrossCorrelationErrorMetric:
1631 default:
1632 {
1633 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1634 channel_distortion,exception);
1635 break;
1636 }
1637 case PeakAbsoluteErrorMetric:
1638 {
1639 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1640 channel_distortion,exception);
1641 break;
1642 }
1643 case PeakSignalToNoiseRatioErrorMetric:
1644 {
1645 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1646 channel_distortion,exception);
1647 break;
1648 }
1649 case PerceptualHashErrorMetric:
1650 {
1651 status=GetPerceptualHashDistortion(image,reconstruct_image,
1652 channel_distortion,exception);
1653 break;
1654 }
1655 case RootMeanSquaredErrorMetric:
1656 {
1657 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1658 channel_distortion,exception);
1659 break;
1660 }
1661 case StructuralSimilarityErrorMetric:
1662 {
1663 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1664 channel_distortion,exception);
1665 break;
1666 }
1667 case StructuralDissimilarityErrorMetric:
1668 {
1669 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1670 channel_distortion,exception);
1671 break;
1672 }
1673 }
1674 *distortion=channel_distortion[CompositePixelChannel];
1675 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1676 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1677 *distortion);
1678 return(status);
1679}
1680␌
1681/*
1682%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1683% %
1684% %
1685% %
1686% G e t I m a g e D i s t o r t i o n s %
1687% %
1688% %
1689% %
1690%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1691%
1692% GetImageDistortions() compares the pixel channels of an image to a
1693% reconstructed image and returns the specified distortion metric for each
1694% channel.
1695%
1696% The format of the GetImageDistortions method is:
1697%
1698% double *GetImageDistortions(const Image *image,
1699% const Image *reconstruct_image,const MetricType metric,
1700% ExceptionInfo *exception)
1701%
1702% A description of each parameter follows:
1703%
1704% o image: the image.
1705%
1706% o reconstruct_image: the reconstruct image.
1707%
1708% o metric: the metric.
1709%
1710% o exception: return any errors or warnings in this structure.
1711%
1712*/
1713MagickExport double *GetImageDistortions(Image *image,
1714 const Image *reconstruct_image,const MetricType metric,
1715 ExceptionInfo *exception)
1716{
1717 double
1718 *channel_distortion;
1719
1720 MagickBooleanType
1721 status;
1722
1723 size_t
1724 length;
1725
1726 assert(image != (Image *) NULL);
1727 assert(image->signature == MagickCoreSignature);
1728 assert(reconstruct_image != (const Image *) NULL);
1729 assert(reconstruct_image->signature == MagickCoreSignature);
1730 if (IsEventLogging() != MagickFalse)
1731 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1732 /*
1733 Get image distortion.
1734 */
1735 length=MaxPixelChannels+1UL;
1736 channel_distortion=(double *) AcquireQuantumMemory(length,
1737 sizeof(*channel_distortion));
1738 if (channel_distortion == (double *) NULL)
1739 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1740 (void) memset(channel_distortion,0,length*
1741 sizeof(*channel_distortion));
1742 status=MagickTrue;
1743 switch (metric)
1744 {
1745 case AbsoluteErrorMetric:
1746 {
1747 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1748 exception);
1749 break;
1750 }
1751 case FuzzErrorMetric:
1752 {
1753 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1754 exception);
1755 break;
1756 }
1757 case MeanAbsoluteErrorMetric:
1758 {
1759 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1760 channel_distortion,exception);
1761 break;
1762 }
1763 case MeanErrorPerPixelErrorMetric:
1764 {
1765 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1766 exception);
1767 break;
1768 }
1769 case MeanSquaredErrorMetric:
1770 {
1771 status=GetMeanSquaredDistortion(image,reconstruct_image,
1772 channel_distortion,exception);
1773 break;
1774 }
1775 case NormalizedCrossCorrelationErrorMetric:
1776 default:
1777 {
1778 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1779 channel_distortion,exception);
1780 break;
1781 }
1782 case PeakAbsoluteErrorMetric:
1783 {
1784 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1785 channel_distortion,exception);
1786 break;
1787 }
1788 case PeakSignalToNoiseRatioErrorMetric:
1789 {
1790 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1791 channel_distortion,exception);
1792 break;
1793 }
1794 case PerceptualHashErrorMetric:
1795 {
1796 status=GetPerceptualHashDistortion(image,reconstruct_image,
1797 channel_distortion,exception);
1798 break;
1799 }
1800 case RootMeanSquaredErrorMetric:
1801 {
1802 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1803 channel_distortion,exception);
1804 break;
1805 }
1806 case StructuralSimilarityErrorMetric:
1807 {
1808 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1809 channel_distortion,exception);
1810 break;
1811 }
1812 case StructuralDissimilarityErrorMetric:
1813 {
1814 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1815 channel_distortion,exception);
1816 break;
1817 }
1818 }
1819 if (status == MagickFalse)
1820 {
1821 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1822 return((double *) NULL);
1823 }
1824 return(channel_distortion);
1825}
1826␌
1827/*
1828%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1829% %
1830% %
1831% %
1832% I s I m a g e s E q u a l %
1833% %
1834% %
1835% %
1836%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1837%
1838% IsImagesEqual() compare the pixels of two images and returns immediately
1839% if any pixel is not identical.
1840%
1841% The format of the IsImagesEqual method is:
1842%
1843% MagickBooleanType IsImagesEqual(const Image *image,
1844% const Image *reconstruct_image,ExceptionInfo *exception)
1845%
1846% A description of each parameter follows.
1847%
1848% o image: the image.
1849%
1850% o reconstruct_image: the reconstruct image.
1851%
1852% o exception: return any errors or warnings in this structure.
1853%
1854*/
1855MagickExport MagickBooleanType IsImagesEqual(const Image *image,
1856 const Image *reconstruct_image,ExceptionInfo *exception)
1857{
1858 CacheView
1859 *image_view,
1860 *reconstruct_view;
1861
1862 size_t
1863 columns,
1864 rows;
1865
1866 ssize_t
1867 y;
1868
1869 assert(image != (Image *) NULL);
1870 assert(image->signature == MagickCoreSignature);
1871 assert(reconstruct_image != (const Image *) NULL);
1872 assert(reconstruct_image->signature == MagickCoreSignature);
1873 rows=MagickMax(image->rows,reconstruct_image->rows);
1874 columns=MagickMax(image->columns,reconstruct_image->columns);
1875 image_view=AcquireVirtualCacheView(image,exception);
1876 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1877 for (y=0; y < (ssize_t) rows; y++)
1878 {
1879 const Quantum
1880 *magick_restrict p,
1881 *magick_restrict q;
1882
1883 ssize_t
1884 x;
1885
1886 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1887 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1888 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1889 break;
1890 for (x=0; x < (ssize_t) columns; x++)
1891 {
1892 ssize_t
1893 i;
1894
1895 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1896 {
1897 double
1898 distance;
1899
1900 PixelChannel channel = GetPixelChannelChannel(image,i);
1901 PixelTrait traits = GetPixelChannelTraits(image,channel);
1902 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1903 channel);
1904 if ((traits == UndefinedPixelTrait) ||
1905 (reconstruct_traits == UndefinedPixelTrait) ||
1906 ((reconstruct_traits & UpdatePixelTrait) == 0))
1907 continue;
1908 distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
1909 channel,q));
1910 if (distance >= MagickEpsilon)
1911 break;
1912 }
1913 if (i < (ssize_t) GetPixelChannels(image))
1914 break;
1915 p+=(ptrdiff_t) GetPixelChannels(image);
1916 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1917 }
1918 if (x < (ssize_t) columns)
1919 break;
1920 }
1921 reconstruct_view=DestroyCacheView(reconstruct_view);
1922 image_view=DestroyCacheView(image_view);
1923 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1924}
1925␌
1926/*
1927%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1928% %
1929% %
1930% %
1931% S e t I m a g e C o l o r M e t r i c %
1932% %
1933% %
1934% %
1935%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1936%
1937% SetImageColorMetric() measures the difference between colors at each pixel
1938% location of two images. A value other than 0 means the colors match
1939% exactly. Otherwise an error measure is computed by summing over all
1940% pixels in an image the distance squared in RGB space between each image
1941% pixel and its corresponding pixel in the reconstruct image. The error
1942% measure is assigned to these image members:
1943%
1944% o mean_error_per_pixel: The mean error for any single pixel in
1945% the image.
1946%
1947% o normalized_mean_error: The normalized mean quantization error for
1948% any single pixel in the image. This distance measure is normalized to
1949% a range between 0 and 1. It is independent of the range of red, green,
1950% and blue values in the image.
1951%
1952% o normalized_maximum_error: The normalized maximum quantization
1953% error for any single pixel in the image. This distance measure is
1954% normalized to a range between 0 and 1. It is independent of the range
1955% of red, green, and blue values in your image.
1956%
1957% A small normalized mean square error, accessed as
1958% image->normalized_mean_error, suggests the images are very similar in
1959% spatial layout and color.
1960%
1961% The format of the SetImageColorMetric method is:
1962%
1963% MagickBooleanType SetImageColorMetric(Image *image,
1964% const Image *reconstruct_image,ExceptionInfo *exception)
1965%
1966% A description of each parameter follows.
1967%
1968% o image: the image.
1969%
1970% o reconstruct_image: the reconstruct image.
1971%
1972% o exception: return any errors or warnings in this structure.
1973%
1974*/
1975MagickExport MagickBooleanType SetImageColorMetric(Image *image,
1976 const Image *reconstruct_image,ExceptionInfo *exception)
1977{
1978 CacheView
1979 *image_view,
1980 *reconstruct_view;
1981
1982 double
1983 area,
1984 maximum_error,
1985 mean_error,
1986 mean_error_per_pixel;
1987
1988 MagickBooleanType
1989 status;
1990
1991 size_t
1992 columns,
1993 rows;
1994
1995 ssize_t
1996 y;
1997
1998 assert(image != (Image *) NULL);
1999 assert(image->signature == MagickCoreSignature);
2000 assert(reconstruct_image != (const Image *) NULL);
2001 assert(reconstruct_image->signature == MagickCoreSignature);
2002 area=0.0;
2003 maximum_error=0.0;
2004 mean_error_per_pixel=0.0;
2005 mean_error=0.0;
2006 rows=MagickMax(image->rows,reconstruct_image->rows);
2007 columns=MagickMax(image->columns,reconstruct_image->columns);
2008 image_view=AcquireVirtualCacheView(image,exception);
2009 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2010 for (y=0; y < (ssize_t) rows; y++)
2011 {
2012 const Quantum
2013 *magick_restrict p,
2014 *magick_restrict q;
2015
2016 ssize_t
2017 x;
2018
2019 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2020 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2021 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2022 break;
2023 for (x=0; x < (ssize_t) columns; x++)
2024 {
2025 ssize_t
2026 i;
2027
2028 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2029 {
2030 double
2031 distance;
2032
2033 PixelChannel channel = GetPixelChannelChannel(image,i);
2034 PixelTrait traits = GetPixelChannelTraits(image,channel);
2035 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2036 channel);
2037 if ((traits == UndefinedPixelTrait) ||
2038 (reconstruct_traits == UndefinedPixelTrait) ||
2039 ((reconstruct_traits & UpdatePixelTrait) == 0))
2040 continue;
2041 distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
2042 channel,q));
2043 if (distance >= MagickEpsilon)
2044 {
2045 mean_error_per_pixel+=distance;
2046 mean_error+=distance*distance;
2047 if (distance > maximum_error)
2048 maximum_error=distance;
2049 }
2050 area++;
2051 }
2052 p+=(ptrdiff_t) GetPixelChannels(image);
2053 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2054 }
2055 }
2056 reconstruct_view=DestroyCacheView(reconstruct_view);
2057 image_view=DestroyCacheView(image_view);
2058 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2059 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
2060 mean_error/area);
2061 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2062 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2063 return(status);
2064}
2065␌
2066/*
2067%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2068% %
2069% %
2070% %
2071% S i m i l a r i t y I m a g e %
2072% %
2073% %
2074% %
2075%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2076%
2077% SimilarityImage() compares the reference image of the image and returns the
2078% best match offset. In addition, it returns a similarity image such that an
2079% exact match location is completely white and if none of the pixels match,
2080% black, otherwise some gray level in-between.
2081%
2082% The format of the SimilarityImageImage method is:
2083%
2084% Image *SimilarityImage(const Image *image,const Image *reference,
2085% const MetricType metric,const double similarity_threshold,
2086% RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2087%
2088% A description of each parameter follows:
2089%
2090% o image: the image.
2091%
2092% o reference: find an area of the image that closely resembles this image.
2093%
2094% o metric: the metric.
2095%
2096% o similarity_threshold: minimum distortion for (sub)image match.
2097%
2098% o offset: the best match offset of the reference image within the image.
2099%
2100% o similarity: the computed similarity between the images.
2101%
2102% o exception: return any errors or warnings in this structure.
2103%
2104*/
2105
2106#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2107static Image *CrossCorrelationImage(const Image *alpha_image,
2108 const Image *beta_image,ExceptionInfo *exception)
2109{
2110 Image
2111 *clone_image,
2112 *complex_conjugate,
2113 *complex_multiplication,
2114 *cross_correlation,
2115 *fft_images;
2116
2117 /*
2118 Take the FFT of beta image.
2119 */
2120 clone_image=CloneImage(beta_image,0,0,MagickTrue,exception);
2121 if (clone_image == (Image *) NULL)
2122 return(clone_image);
2123 (void) SetImageArtifact(clone_image,"fourier:normalize","inverse");
2124 fft_images=ForwardFourierTransformImage(clone_image,MagickFalse,
2125 exception);
2126 clone_image=DestroyImageList(clone_image);
2127 if (fft_images == (Image *) NULL)
2128 return(fft_images);
2129 /*
2130 Take the complex conjugate of beta image.
2131 */
2132 complex_conjugate=ComplexImages(fft_images,ConjugateComplexOperator,
2133 exception);
2134 fft_images=DestroyImageList(fft_images);
2135 if (complex_conjugate == (Image *) NULL)
2136 return(complex_conjugate);
2137 /*
2138 Take the FFT of the alpha image.
2139 */
2140 clone_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2141 if (clone_image == (Image *) NULL)
2142 {
2143 complex_conjugate=DestroyImageList(complex_conjugate);
2144 return(clone_image);
2145 }
2146 (void) SetImageArtifact(clone_image,"fourier:normalize","inverse");
2147 fft_images=ForwardFourierTransformImage(clone_image,MagickFalse,exception);
2148 clone_image=DestroyImageList(clone_image);
2149 if (fft_images == (Image *) NULL)
2150 {
2151 complex_conjugate=DestroyImageList(complex_conjugate);
2152 return(fft_images);
2153 }
2154 complex_conjugate->next->next=fft_images;
2155 /*
2156 Do complex multiplication.
2157 */
2158 DisableCompositeClampUnlessSpecified(complex_conjugate);
2159 complex_multiplication=ComplexImages(complex_conjugate,
2160 MultiplyComplexOperator,exception);
2161 complex_conjugate=DestroyImageList(complex_conjugate);
2162 if (fft_images == (Image *) NULL)
2163 return(fft_images);
2164 /*
2165 Do the IFT and return the cross-correlation result.
2166 */
2167 cross_correlation=InverseFourierTransformImage(complex_multiplication,
2168 complex_multiplication->next,MagickFalse,exception);
2169 complex_multiplication=DestroyImageList(complex_multiplication);
2170 return(cross_correlation);
2171}
2172
2173static Image *NCCDivideImage(const Image *alpha_image,const Image *beta_image,
2174 ExceptionInfo *exception)
2175{
2176 CacheView
2177 *alpha_view,
2178 *beta_view;
2179
2180 Image
2181 *divide_image;
2182
2183 MagickBooleanType
2184 status;
2185
2186 ssize_t
2187 y;
2188
2189 /*
2190 Divide one image into another.
2191 */
2192 divide_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2193 if (divide_image == (Image *) NULL)
2194 return(divide_image);
2195 status=MagickTrue;
2196 alpha_view=AcquireAuthenticCacheView(divide_image,exception);
2197 beta_view=AcquireVirtualCacheView(beta_image,exception);
2198#if defined(MAGICKCORE_OPENMP_SUPPORT)
2199 #pragma omp parallel for schedule(static) shared(status) \
2200 magick_number_threads(beta_image,divide_image,divide_image->rows,1)
2201#endif
2202 for (y=0; y < (ssize_t) divide_image->rows; y++)
2203 {
2204 const Quantum
2205 *magick_restrict p;
2206
2207 Quantum
2208 *magick_restrict q;
2209
2210 ssize_t
2211 x;
2212
2213 if (status == MagickFalse)
2214 continue;
2215 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
2216 exception);
2217 q=GetCacheViewAuthenticPixels(alpha_view,0,y,divide_image->columns,1,
2218 exception);
2219 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2220 {
2221 status=MagickFalse;
2222 continue;
2223 }
2224 for (x=0; x < (ssize_t) divide_image->columns; x++)
2225 {
2226 ssize_t
2227 i;
2228
2229 for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
2230 {
2231 PixelChannel channel = GetPixelChannelChannel(divide_image,i);
2232 PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
2233 if ((traits & UpdatePixelTrait) == 0)
2234 continue;
2235 if (fabs(p[i]) >= MagickEpsilon)
2236 q[i]=(Quantum) ((double) q[i]*PerceptibleReciprocal(QuantumScale*
2237 (double) p[i]));
2238 }
2239 p+=(ptrdiff_t) GetPixelChannels(beta_image);
2240 q+=(ptrdiff_t) GetPixelChannels(divide_image);
2241 }
2242 if (SyncCacheViewAuthenticPixels(alpha_view,exception) == MagickFalse)
2243 status=MagickFalse;
2244 }
2245 beta_view=DestroyCacheView(beta_view);
2246 alpha_view=DestroyCacheView(alpha_view);
2247 if (status == MagickFalse)
2248 divide_image=DestroyImage(divide_image);
2249 return(divide_image);
2250}
2251
2252static MagickBooleanType NCCMaximaImage(const Image *image,double *maxima,
2253 RectangleInfo *offset,ExceptionInfo *exception)
2254{
2255 CacheView
2256 *image_view;
2257
2258 MagickBooleanType
2259 status;
2260
2261 ssize_t
2262 y;
2263
2264 /*
2265 Identify the maxima value in the image and its location.
2266 */
2267 status=MagickTrue;
2268 *maxima=0.0;
2269 offset->x=0;
2270 offset->y=0;
2271 image_view=AcquireVirtualCacheView(image,exception);
2272 for (y=0; y < (ssize_t) image->rows; y++)
2273 {
2274 const Quantum
2275 *magick_restrict p;
2276
2277 ssize_t
2278 x;
2279
2280 if (status == MagickFalse)
2281 continue;
2282 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2283 if (p == (const Quantum *) NULL)
2284 {
2285 status=MagickFalse;
2286 continue;
2287 }
2288 for (x=0; x < (ssize_t) image->columns; x++)
2289 {
2290 double
2291 sum = 0.0;
2292
2293 ssize_t
2294 channels = 0,
2295 i;
2296
2297 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2298 {
2299 PixelChannel channel = GetPixelChannelChannel(image,i);
2300 PixelTrait traits = GetPixelChannelTraits(image,channel);
2301 if ((traits & UpdatePixelTrait) == 0)
2302 continue;
2303 sum+=(double) p[i];
2304 channels++;
2305 }
2306 if ((channels != 0) && ((sum/channels) > *maxima))
2307 {
2308 *maxima=sum/channels;
2309 offset->x=x;
2310 offset->y=y;
2311 }
2312 p+=(ptrdiff_t) GetPixelChannels(image);
2313 }
2314 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2315 status=MagickFalse;
2316 }
2317 image_view=DestroyCacheView(image_view);
2318 return(status);
2319}
2320
2321static MagickBooleanType NCCMultiplyImage(Image *image,const double factor,
2322 const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
2323{
2324 CacheView
2325 *image_view;
2326
2327 MagickBooleanType
2328 status;
2329
2330 ssize_t
2331 y;
2332
2333 /*
2334 Multiply each pixel by a factor.
2335 */
2336 status=MagickTrue;
2337 image_view=AcquireAuthenticCacheView(image,exception);
2338#if defined(MAGICKCORE_OPENMP_SUPPORT)
2339 #pragma omp parallel for schedule(static) shared(status) \
2340 magick_number_threads(image,image,image->rows,1)
2341#endif
2342 for (y=0; y < (ssize_t) image->rows; y++)
2343 {
2344 Quantum
2345 *magick_restrict q;
2346
2347 ssize_t
2348 x;
2349
2350 if (status == MagickFalse)
2351 continue;
2352 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2353 if (q == (Quantum *) NULL)
2354 {
2355 status=MagickFalse;
2356 continue;
2357 }
2358 for (x=0; x < (ssize_t) image->columns; x++)
2359 {
2360 ssize_t
2361 i;
2362
2363 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2364 {
2365 PixelChannel channel = GetPixelChannelChannel(image,i);
2366 PixelTrait traits = GetPixelChannelTraits(image,channel);
2367 if ((traits & UpdatePixelTrait) == 0)
2368 continue;
2369 if (channel_statistics != (const ChannelStatistics *) NULL)
2370 q[i]=(Quantum) ((double) q[i]*QuantumScale*
2371 channel_statistics[channel].standard_deviation);
2372 q[i]=(Quantum) ((double) q[i]*factor);
2373 }
2374 q+=(ptrdiff_t) GetPixelChannels(image);
2375 }
2376 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2377 status=MagickFalse;
2378 }
2379 image_view=DestroyCacheView(image_view);
2380 return(status);
2381}
2382
2383static Image *NCCSquareImage(const Image *image,ExceptionInfo *exception)
2384{
2385 CacheView
2386 *image_view;
2387
2388 Image
2389 *square_image;
2390
2391 MagickBooleanType
2392 status;
2393
2394 ssize_t
2395 y;
2396
2397 /*
2398 Square each pixel in the image.
2399 */
2400 square_image=CloneImage(image,0,0,MagickTrue,exception);
2401 if (square_image == (Image *) NULL)
2402 return(square_image);
2403 status=MagickTrue;
2404 image_view=AcquireAuthenticCacheView(square_image,exception);
2405#if defined(MAGICKCORE_OPENMP_SUPPORT)
2406 #pragma omp parallel for schedule(static) shared(status) \
2407 magick_number_threads(square_image,square_image,square_image->rows,1)
2408#endif
2409 for (y=0; y < (ssize_t) square_image->rows; y++)
2410 {
2411 Quantum
2412 *magick_restrict q;
2413
2414 ssize_t
2415 x;
2416
2417 if (status == MagickFalse)
2418 continue;
2419 q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
2420 exception);
2421 if (q == (Quantum *) NULL)
2422 {
2423 status=MagickFalse;
2424 continue;
2425 }
2426 for (x=0; x < (ssize_t) square_image->columns; x++)
2427 {
2428 ssize_t
2429 i;
2430
2431 for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
2432 {
2433 PixelChannel channel = GetPixelChannelChannel(square_image,i);
2434 PixelTrait traits = GetPixelChannelTraits(square_image,channel);
2435 if ((traits & UpdatePixelTrait) == 0)
2436 continue;
2437 q[i]=(Quantum) ((double) q[i]*QuantumScale*(double) q[i]);
2438 }
2439 q+=(ptrdiff_t) GetPixelChannels(square_image);
2440 }
2441 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2442 status=MagickFalse;
2443 }
2444 image_view=DestroyCacheView(image_view);
2445 if (status == MagickFalse)
2446 square_image=DestroyImage(square_image);
2447 return(square_image);
2448}
2449
2450static Image *NCCSubtractImageMean(const Image *alpha_image,
2451 const Image *beta_image,const ChannelStatistics *channel_statistics,
2452 ExceptionInfo *exception)
2453{
2454 CacheView
2455 *beta_view,
2456 *image_view;
2457
2458 Image
2459 *gamma_image;
2460
2461 MagickBooleanType
2462 status;
2463
2464 ssize_t
2465 y;
2466
2467 /*
2468 Subtract the image mean and pad.
2469 */
2470 gamma_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
2471 MagickTrue,exception);
2472 if (gamma_image == (Image *) NULL)
2473 return(gamma_image);
2474 status=MagickTrue;
2475 image_view=AcquireAuthenticCacheView(gamma_image,exception);
2476 beta_view=AcquireVirtualCacheView(beta_image,exception);
2477#if defined(MAGICKCORE_OPENMP_SUPPORT)
2478 #pragma omp parallel for schedule(static) shared(status) \
2479 magick_number_threads(beta_image,gamma_image,gamma_image->rows,1)
2480#endif
2481 for (y=0; y < (ssize_t) gamma_image->rows; y++)
2482 {
2483 const Quantum
2484 *magick_restrict p;
2485
2486 Quantum
2487 *magick_restrict q;
2488
2489 ssize_t
2490 x;
2491
2492 if (status == MagickFalse)
2493 continue;
2494 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
2495 exception);
2496 q=GetCacheViewAuthenticPixels(image_view,0,y,gamma_image->columns,1,
2497 exception);
2498 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2499 {
2500 status=MagickFalse;
2501 continue;
2502 }
2503 for (x=0; x < (ssize_t) gamma_image->columns; x++)
2504 {
2505 ssize_t
2506 i;
2507
2508 for (i=0; i < (ssize_t) GetPixelChannels(gamma_image); i++)
2509 {
2510 PixelChannel channel = GetPixelChannelChannel(gamma_image,i);
2511 PixelTrait traits = GetPixelChannelTraits(gamma_image,channel);
2512 if ((traits & UpdatePixelTrait) == 0)
2513 continue;
2514 if ((x >= (ssize_t) beta_image->columns) ||
2515 (y >= (ssize_t) beta_image->rows))
2516 q[i]=(Quantum) 0;
2517 else
2518 q[i]=(Quantum) ((double) p[i]-channel_statistics[channel].mean);
2519 }
2520 p+=(ptrdiff_t) GetPixelChannels(beta_image);
2521 q+=(ptrdiff_t) GetPixelChannels(gamma_image);
2522 }
2523 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2524 status=MagickFalse;
2525 }
2526 beta_view=DestroyCacheView(beta_view);
2527 image_view=DestroyCacheView(image_view);
2528 if (status == MagickFalse)
2529 gamma_image=DestroyImage(gamma_image);
2530 return(gamma_image);
2531}
2532
2533static Image *NCCUnityImage(const Image *alpha_image,const Image *beta_image,
2534 ExceptionInfo *exception)
2535{
2536 CacheView
2537 *image_view;
2538
2539 Image
2540 *unity_image;
2541
2542 MagickBooleanType
2543 status;
2544
2545 ssize_t
2546 y;
2547
2548 /*
2549 Create a padded unity image.
2550 */
2551 unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
2552 MagickTrue,exception);
2553 if (unity_image == (Image *) NULL)
2554 return(unity_image);
2555 if (SetImageStorageClass(unity_image,DirectClass,exception) == MagickFalse)
2556 return(DestroyImage(unity_image));
2557 status=MagickTrue;
2558 image_view=AcquireAuthenticCacheView(unity_image,exception);
2559#if defined(MAGICKCORE_OPENMP_SUPPORT)
2560 #pragma omp parallel for schedule(static) shared(status) \
2561 magick_number_threads(unity_image,unity_image,unity_image->rows,1)
2562#endif
2563 for (y=0; y < (ssize_t) unity_image->rows; y++)
2564 {
2565 Quantum
2566 *magick_restrict q;
2567
2568 ssize_t
2569 x;
2570
2571 if (status == MagickFalse)
2572 continue;
2573 q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
2574 exception);
2575 if (q == (Quantum *) NULL)
2576 {
2577 status=MagickFalse;
2578 continue;
2579 }
2580 for (x=0; x < (ssize_t) unity_image->columns; x++)
2581 {
2582 ssize_t
2583 i;
2584
2585 for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
2586 {
2587 PixelChannel channel = GetPixelChannelChannel(unity_image,i);
2588 PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
2589 if ((traits & UpdatePixelTrait) == 0)
2590 continue;
2591 q[i]=QuantumRange;
2592 if ((x >= (ssize_t) beta_image->columns) ||
2593 (y >= (ssize_t) beta_image->rows))
2594 q[i]=(Quantum) 0;
2595 }
2596 q+=(ptrdiff_t) GetPixelChannels(unity_image);
2597 }
2598 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2599 status=MagickFalse;
2600 }
2601 image_view=DestroyCacheView(image_view);
2602 if (status == MagickFalse)
2603 unity_image=DestroyImage(unity_image);
2604 return(unity_image);
2605}
2606
2607static Image *NCCVarianceImage(Image *alpha_image,const Image *beta_image,
2608 ExceptionInfo *exception)
2609{
2610 CacheView
2611 *beta_view,
2612 *image_view;
2613
2614 Image
2615 *variance_image;
2616
2617 MagickBooleanType
2618 status;
2619
2620 ssize_t
2621 y;
2622
2623 /*
2624 Compute the variance of the two images.
2625 */
2626 variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2627 if (variance_image == (Image *) NULL)
2628 return(variance_image);
2629 status=MagickTrue;
2630 image_view=AcquireAuthenticCacheView(variance_image,exception);
2631 beta_view=AcquireVirtualCacheView(beta_image,exception);
2632#if defined(MAGICKCORE_OPENMP_SUPPORT)
2633 #pragma omp parallel for schedule(static) shared(status) \
2634 magick_number_threads(beta_image,variance_image,variance_image->rows,1)
2635#endif
2636 for (y=0; y < (ssize_t) variance_image->rows; y++)
2637 {
2638 const Quantum
2639 *magick_restrict p;
2640
2641 Quantum
2642 *magick_restrict q;
2643
2644 ssize_t
2645 x;
2646
2647 if (status == MagickFalse)
2648 continue;
2649 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
2650 exception);
2651 q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
2652 exception);
2653 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2654 {
2655 status=MagickFalse;
2656 continue;
2657 }
2658 for (x=0; x < (ssize_t) variance_image->columns; x++)
2659 {
2660 ssize_t
2661 i;
2662
2663 for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
2664 {
2665 PixelChannel channel = GetPixelChannelChannel(variance_image,i);
2666 PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
2667 if ((traits & UpdatePixelTrait) == 0)
2668 continue;
2669 q[i]=(Quantum) ((double) ClampToQuantum((double) QuantumRange*sqrt(fabs(
2670 QuantumScale*((double) q[i]-(double) p[i]))))/
2671 sqrt((double) QuantumRange));
2672 }
2673 p+=(ptrdiff_t) GetPixelChannels(beta_image);
2674 q+=(ptrdiff_t) GetPixelChannels(variance_image);
2675 }
2676 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2677 status=MagickFalse;
2678 }
2679 beta_view=DestroyCacheView(beta_view);
2680 image_view=DestroyCacheView(image_view);
2681 if (status == MagickFalse)
2682 variance_image=DestroyImage(variance_image);
2683 return(variance_image);
2684}
2685
2686static Image *NCCSimilarityImage(const Image *image,const Image *reference,
2687 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2688{
2689#define DestroySimilarityResources() \
2690{ \
2691 if (channel_statistics != (ChannelStatistics *) NULL) \
2692 channel_statistics=(ChannelStatistics *) \
2693 RelinquishMagickMemory(channel_statistics); \
2694 if (beta_image != (Image *) NULL) \
2695 beta_image=DestroyImage(beta_image); \
2696 if (gamma_image != (Image *) NULL) \
2697 gamma_image=DestroyImage(gamma_image); \
2698 if (ncc_image != (Image *) NULL) \
2699 ncc_image=DestroyImage(ncc_image); \
2700 if (normalize_image != (Image *) NULL) \
2701 normalize_image=DestroyImage(normalize_image); \
2702 if (square_image != (Image *) NULL) \
2703 square_image=DestroyImage(square_image); \
2704 if (unity_image != (Image *) NULL) \
2705 unity_image=DestroyImage(unity_image); \
2706}
2707#define ThrowSimilarityException() \
2708{ \
2709 DestroySimilarityResources() \
2710 return((Image *) NULL); \
2711}
2712
2713 ChannelStatistics
2714 *channel_statistics = (ChannelStatistics *) NULL;
2715
2716 double
2717 maxima = 0.0;
2718
2719 Image
2720 *beta_image = (Image *) NULL,
2721 *correlation_image = (Image *) NULL,
2722 *gamma_image = (Image *) NULL,
2723 *ncc_image = (Image *) NULL,
2724 *normalize_image = (Image *) NULL,
2725 *square_image = (Image *) NULL,
2726 *unity_image = (Image *) NULL;
2727
2728 MagickBooleanType
2729 status;
2730
2731 RectangleInfo
2732 geometry;
2733
2734 /*
2735 Accelerated correlation-based image similary using FFT local statistics.
2736 Contributed by Fred Weinhaus.
2737 */
2738 square_image=NCCSquareImage(image,exception);
2739 if (square_image == (Image *) NULL)
2740 ThrowSimilarityException();
2741 unity_image=NCCUnityImage(image,reference,exception);
2742 if (unity_image == (Image *) NULL)
2743 ThrowSimilarityException();
2744 /*
2745 Compute the cross correlation of the square and unity images.
2746 */
2747 ncc_image=CrossCorrelationImage(square_image,unity_image,exception);
2748 square_image=DestroyImage(square_image);
2749 if (ncc_image == (Image *) NULL)
2750 ThrowSimilarityException();
2751 status=NCCMultiplyImage(ncc_image,(double) QuantumRange*reference->columns*
2752 reference->rows,(const ChannelStatistics *) NULL,exception);
2753 if (status == MagickFalse)
2754 ThrowSimilarityException();
2755 /*
2756 Compute the cross correlation of the source and unity images.
2757 */
2758 gamma_image=CrossCorrelationImage(image,unity_image,exception);
2759 unity_image=DestroyImage(unity_image);
2760 if (gamma_image == (Image *) NULL)
2761 ThrowSimilarityException();
2762 square_image=NCCSquareImage(gamma_image,exception);
2763 gamma_image=DestroyImage(gamma_image);
2764 status=NCCMultiplyImage(square_image,(double) QuantumRange,
2765 (const ChannelStatistics *) NULL,exception);
2766 if (status == MagickFalse)
2767 ThrowSimilarityException();
2768 /*
2769 Compute the variance of the two images.
2770 */
2771 gamma_image=NCCVarianceImage(ncc_image,square_image,exception);
2772 square_image=DestroyImage(square_image);
2773 ncc_image=DestroyImage(ncc_image);
2774 if (gamma_image == (Image *) NULL)
2775 ThrowSimilarityException();
2776 channel_statistics=GetImageStatistics(reference,exception);
2777 if (channel_statistics == (ChannelStatistics *) NULL)
2778 ThrowSimilarityException();
2779 /*
2780 Subtract the image mean.
2781 */
2782 status=NCCMultiplyImage(gamma_image,1.0,channel_statistics,exception);
2783 if (status == MagickFalse)
2784 ThrowSimilarityException();
2785 normalize_image=NCCSubtractImageMean(image,reference,channel_statistics,
2786 exception);
2787 if (normalize_image == (Image *) NULL)
2788 ThrowSimilarityException();
2789 ncc_image=CrossCorrelationImage(image,normalize_image,exception);
2790 normalize_image=DestroyImage(normalize_image);
2791 if (ncc_image == (Image *) NULL)
2792 ThrowSimilarityException();
2793 /*
2794 Divide the two images.
2795 */
2796 beta_image=NCCDivideImage(ncc_image,gamma_image,exception);
2797 ncc_image=DestroyImage(ncc_image);
2798 gamma_image=DestroyImage(gamma_image);
2799 if (beta_image == (Image *) NULL)
2800 ThrowSimilarityException();
2801 (void) ResetImagePage(beta_image,"0x0+0+0");
2802 SetGeometry(image,&geometry);
2803 geometry.width=image->columns-reference->columns;
2804 geometry.height=image->rows-reference->rows;
2805 /*
2806 Crop padding.
2807 */
2808 correlation_image=CropImage(beta_image,&geometry,exception);
2809 beta_image=DestroyImage(beta_image);
2810 if (correlation_image == (Image *) NULL)
2811 ThrowSimilarityException();
2812 (void) ResetImagePage(correlation_image,"0x0+0+0");
2813 /*
2814 Identify the maxima value in the image and its location.
2815 */
2816 status=GrayscaleImage(correlation_image,AveragePixelIntensityMethod,
2817 exception);
2818 if (status == MagickFalse)
2819 ThrowSimilarityException();
2820 status=NCCMaximaImage(correlation_image,&maxima,offset,exception);
2821 if (status == MagickFalse)
2822 {
2823 correlation_image=DestroyImage(correlation_image);
2824 ThrowSimilarityException();
2825 }
2826 *similarity_metric=1.0-QuantumScale*maxima;
2827 DestroySimilarityResources();
2828 return(correlation_image);
2829}
2830#endif
2831
2832static double GetSimilarityMetric(const Image *image,const Image *reference,
2833 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2834 ExceptionInfo *exception)
2835{
2836 double
2837 distortion;
2838
2839 Image
2840 *similarity_image;
2841
2842 MagickBooleanType
2843 status;
2844
2845 RectangleInfo
2846 geometry;
2847
2848 SetGeometry(reference,&geometry);
2849 geometry.x=x_offset;
2850 geometry.y=y_offset;
2851 similarity_image=CropImage(image,&geometry,exception);
2852 if (similarity_image == (Image *) NULL)
2853 return(0.0);
2854 distortion=0.0;
2855 status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2856 exception);
2857 similarity_image=DestroyImage(similarity_image);
2858 if (status == MagickFalse)
2859 return(0.0);
2860 return(distortion);
2861}
2862
2863MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2864 const MetricType metric,const double similarity_threshold,
2865 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2866{
2867#define SimilarityImageTag "Similarity/Image"
2868
2869 CacheView
2870 *similarity_view;
2871
2872 Image
2873 *similarity_image = (Image *) NULL;
2874
2875 MagickBooleanType
2876 status;
2877
2878 MagickOffsetType
2879 progress;
2880
2881 ssize_t
2882 y;
2883
2884 assert(image != (const Image *) NULL);
2885 assert(image->signature == MagickCoreSignature);
2886 assert(exception != (ExceptionInfo *) NULL);
2887 assert(exception->signature == MagickCoreSignature);
2888 assert(offset != (RectangleInfo *) NULL);
2889 if (IsEventLogging() != MagickFalse)
2890 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2891 SetGeometry(reference,offset);
2892 *similarity_metric=MagickMaximumValue;
2893#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2894 if ((image->channels & ReadMaskChannel) == 0)
2895 {
2896 const char *artifact = GetImageArtifact(image,"compare:accelerate-ncc");
2897 MagickBooleanType accelerate = (artifact != (const char *) NULL) &&
2898 (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
2899 if ((accelerate != MagickFalse) &&
2900 (metric == NormalizedCrossCorrelationErrorMetric))
2901 {
2902 similarity_image=NCCSimilarityImage(image,reference,offset,
2903 similarity_metric,exception);
2904 return(similarity_image);
2905 }
2906 }
2907#endif
2908 if ((image->columns >= reference->columns) &&
2909 (image->rows >= reference->rows))
2910 similarity_image=CloneImage(image,image->columns-reference->columns+1,
2911 image->rows-reference->rows+1,MagickTrue,exception);
2912 if (similarity_image == (Image *) NULL)
2913 return((Image *) NULL);
2914 status=SetImageStorageClass(similarity_image,DirectClass,exception);
2915 if (status == MagickFalse)
2916 {
2917 similarity_image=DestroyImage(similarity_image);
2918 return((Image *) NULL);
2919 }
2920 (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2921 exception);
2922 /*
2923 Measure similarity of reference image against image.
2924 */
2925 status=MagickTrue;
2926 progress=0;
2927 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2928#if defined(MAGICKCORE_OPENMP_SUPPORT)
2929 #pragma omp parallel for schedule(static) \
2930 shared(progress,status,similarity_metric) \
2931 magick_number_threads(image,image,image->rows-reference->rows+1,1)
2932#endif
2933 for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2934 {
2935 double
2936 similarity;
2937
2938 Quantum
2939 *magick_restrict q;
2940
2941 ssize_t
2942 x;
2943
2944 if (status == MagickFalse)
2945 continue;
2946#if defined(MAGICKCORE_OPENMP_SUPPORT)
2947 #pragma omp flush(similarity_metric)
2948#endif
2949 if (*similarity_metric <= similarity_threshold)
2950 continue;
2951 q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2952 1,exception);
2953 if (q == (Quantum *) NULL)
2954 {
2955 status=MagickFalse;
2956 continue;
2957 }
2958 for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2959 {
2960 ssize_t
2961 i;
2962
2963#if defined(MAGICKCORE_OPENMP_SUPPORT)
2964 #pragma omp flush(similarity_metric)
2965#endif
2966 if (*similarity_metric <= similarity_threshold)
2967 break;
2968 similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2969 if (metric == PeakSignalToNoiseRatioErrorMetric)
2970 similarity*=0.01;
2971 if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2972 (metric == UndefinedErrorMetric))
2973 similarity=1.0-similarity;
2974#if defined(MAGICKCORE_OPENMP_SUPPORT)
2975 #pragma omp critical (MagickCore_SimilarityImage)
2976#endif
2977 if (similarity < *similarity_metric)
2978 {
2979 offset->x=x;
2980 offset->y=y;
2981 *similarity_metric=similarity;
2982 }
2983 if (metric == PerceptualHashErrorMetric)
2984 similarity=MagickMin(0.01*similarity,1.0);
2985 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2986 {
2987 PixelChannel channel = GetPixelChannelChannel(image,i);
2988 PixelTrait traits = GetPixelChannelTraits(image,channel);
2989 PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2990 channel);
2991 if ((traits == UndefinedPixelTrait) ||
2992 (similarity_traits == UndefinedPixelTrait) ||
2993 ((similarity_traits & UpdatePixelTrait) == 0))
2994 continue;
2995 SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
2996 QuantumRange-(double) QuantumRange*similarity),q);
2997 }
2998 q+=(ptrdiff_t) GetPixelChannels(similarity_image);
2999 }
3000 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
3001 status=MagickFalse;
3002 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3003 {
3004 MagickBooleanType
3005 proceed;
3006
3007#if defined(MAGICKCORE_OPENMP_SUPPORT)
3008 #pragma omp atomic
3009#endif
3010 progress++;
3011 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
3012 if (proceed == MagickFalse)
3013 status=MagickFalse;
3014 }
3015 }
3016 similarity_view=DestroyCacheView(similarity_view);
3017 if (status == MagickFalse)
3018 similarity_image=DestroyImage(similarity_image);
3019 return(similarity_image);
3020}
_CacheView
Definition cache-view.c:66
_ChannelPerceptualHash
Definition statistic.h:74
_ChannelStatistics
Definition statistic.h:30
_ExceptionInfo
Definition exception.h:102
_Image
Definition image.h:132
_KernelInfo
Definition morphology.h:103
_PixelInfo
Definition pixel.h:182
_RectangleInfo
Definition geometry.h:130