MagickCore 7.1.2-25
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/license/ %
27% %
28% Unless required by applicable law or agreed to in writing, software %
29% distributed under the License is distributed on an "AS IS" BASIS, %
30% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31% See the License for the specific language governing permissions and %
32% limitations under the License. %
33% %
34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35%
36%
37%
38*/
39␌
40/*
41 Include declarations.
42*/
43#include "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/compare-private.h"
55#include "MagickCore/composite-private.h"
56#include "MagickCore/constitute.h"
57#include "MagickCore/distort.h"
58#include "MagickCore/exception-private.h"
59#include "MagickCore/enhance.h"
60#include "MagickCore/fourier.h"
61#include "MagickCore/geometry.h"
62#include "MagickCore/image-private.h"
63#include "MagickCore/list.h"
64#include "MagickCore/log.h"
65#include "MagickCore/memory_.h"
66#include "MagickCore/monitor.h"
67#include "MagickCore/monitor-private.h"
68#include "MagickCore/option.h"
69#include "MagickCore/pixel-accessor.h"
70#include "MagickCore/property.h"
71#include "MagickCore/registry.h"
72#include "MagickCore/resource_.h"
73#include "MagickCore/string_.h"
74#include "MagickCore/statistic.h"
75#include "MagickCore/statistic-private.h"
76#include "MagickCore/string-private.h"
77#include "MagickCore/thread-private.h"
78#include "MagickCore/threshold.h"
79#include "MagickCore/transform.h"
80#include "MagickCore/utility.h"
81#include "MagickCore/version.h"
82␌
83/*
84%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85% %
86% %
87% %
88% C o m p a r e I m a g e s %
89% %
90% %
91% %
92%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93%
94% CompareImages() compares one or more pixel channels of an image to a
95% reconstructed image and returns the difference image.
96%
97% The format of the CompareImages method is:
98%
99% Image *CompareImages(const Image *image,const Image *reconstruct_image,
100% const MetricType metric,double *distortion,ExceptionInfo *exception)
101%
102% A description of each parameter follows:
103%
104% o image: the image.
105%
106% o reconstruct_image: the reconstruction image.
107%
108% o metric: the metric.
109%
110% o distortion: the computed distortion between the images.
111%
112% o exception: return any errors or warnings in this structure.
113%
114*/
115MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
116 const MetricType metric,double *distortion,ExceptionInfo *exception)
117{
118 CacheView
119 *highlight_view,
120 *image_view,
121 *reconstruct_view;
122
123 const char
124 *artifact;
125
126 Image
127 *clone_image,
128 *difference_image,
129 *highlight_image;
130
131 MagickBooleanType
132 status = MagickTrue;
133
134 PixelInfo
135 highlight,
136 lowlight,
137 masklight;
138
139 RectangleInfo
140 geometry;
141
142 size_t
143 columns,
144 rows;
145
146 ssize_t
147 y;
148
149 assert(image != (Image *) NULL);
150 assert(image->signature == MagickCoreSignature);
151 assert(reconstruct_image != (const Image *) NULL);
152 assert(reconstruct_image->signature == MagickCoreSignature);
153 assert(distortion != (double *) NULL);
154 if (IsEventLogging() != MagickFalse)
155 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
156 *distortion=0.0;
157 status=GetImageDistortion(image,reconstruct_image,metric,distortion,
158 exception);
159 if (status == MagickFalse)
160 return((Image *) NULL);
161 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
162 SetGeometry(image,&geometry);
163 geometry.width=columns;
164 geometry.height=rows;
165 clone_image=CloneImage(image,0,0,MagickTrue,exception);
166 if (clone_image == (Image *) NULL)
167 return((Image *) NULL);
168 (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
169 difference_image=ExtentImage(clone_image,&geometry,exception);
170 clone_image=DestroyImage(clone_image);
171 if (difference_image == (Image *) NULL)
172 return((Image *) NULL);
173 (void) ResetImagePage(difference_image,"0x0+0+0");
174 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
175 highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
176 if (highlight_image == (Image *) NULL)
177 {
178 difference_image=DestroyImage(difference_image);
179 return((Image *) NULL);
180 }
181 status=SetImageStorageClass(highlight_image,DirectClass,exception);
182 if (status == MagickFalse)
183 {
184 difference_image=DestroyImage(difference_image);
185 highlight_image=DestroyImage(highlight_image);
186 return((Image *) NULL);
187 }
188 (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
189 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
190 (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
191 artifact=GetImageArtifact(image,"compare:highlight-color");
192 if (artifact != (const char *) NULL)
193 (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
194 (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
195 artifact=GetImageArtifact(image,"compare:lowlight-color");
196 if (artifact != (const char *) NULL)
197 (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
198 (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
199 artifact=GetImageArtifact(image,"compare:masklight-color");
200 if (artifact != (const char *) NULL)
201 (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
202 /*
203 Generate difference image.
204 */
205 image_view=AcquireVirtualCacheView(image,exception);
206 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
207 highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
208#if defined(MAGICKCORE_OPENMP_SUPPORT)
209 #pragma omp parallel for schedule(static) shared(status) \
210 magick_number_threads(image,highlight_image,rows,1)
211#endif
212 for (y=0; y < (ssize_t) rows; y++)
213 {
214 const Quantum
215 *magick_restrict p,
216 *magick_restrict q;
217
218 MagickBooleanType
219 sync;
220
221 Quantum
222 *magick_restrict r;
223
224 ssize_t
225 x;
226
227 if (status == MagickFalse)
228 continue;
229 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
230 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
231 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
232 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
233 (r == (Quantum *) NULL))
234 {
235 status=MagickFalse;
236 continue;
237 }
238 for (x=0; x < (ssize_t) columns; x++)
239 {
240 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
241 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
242 {
243 SetPixelViaPixelInfo(highlight_image,&masklight,r);
244 p+=(ptrdiff_t) GetPixelChannels(image);
245 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
246 r+=(ptrdiff_t) GetPixelChannels(highlight_image);
247 continue;
248 }
249 if (IsFuzzyEquivalencePixel(image,p,reconstruct_image,q) == MagickFalse)
250 SetPixelViaPixelInfo(highlight_image,&highlight,r);
251 else
252 SetPixelViaPixelInfo(highlight_image,&lowlight,r);
253 p+=(ptrdiff_t) GetPixelChannels(image);
254 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
255 r+=(ptrdiff_t) GetPixelChannels(highlight_image);
256 }
257 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
258 if (sync == MagickFalse)
259 status=MagickFalse;
260 }
261 highlight_view=DestroyCacheView(highlight_view);
262 reconstruct_view=DestroyCacheView(reconstruct_view);
263 image_view=DestroyCacheView(image_view);
264 if ((status != MagickFalse) && (difference_image != (Image *) NULL))
265 status=CompositeImage(difference_image,highlight_image,image->compose,
266 MagickTrue,0,0,exception);
267 highlight_image=DestroyImage(highlight_image);
268 if (status == MagickFalse)
269 difference_image=DestroyImage(difference_image);
270 return(difference_image);
271}
272␌
273/*
274%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
275% %
276% %
277% %
278% G e t I m a g e D i s t o r t i o n %
279% %
280% %
281% %
282%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283%
284% GetImageDistortion() compares one or more pixel channels of an image to a
285% reconstructed image and returns the specified distortion metric.
286%
287% The format of the GetImageDistortion method is:
288%
289% MagickBooleanType GetImageDistortion(const Image *image,
290% const Image *reconstruct_image,const MetricType metric,
291% double *distortion,ExceptionInfo *exception)
292%
293% A description of each parameter follows:
294%
295% o image: the image.
296%
297% o reconstruct_image: the reconstruction image.
298%
299% o metric: the metric.
300%
301% o distortion: the computed distortion between the images.
302%
303% o exception: return any errors or warnings in this structure.
304%
305*/
306
307static MagickBooleanType GetAESimilarity(const Image *image,
308 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
309{
310 CacheView
311 *image_view,
312 *reconstruct_view;
313
314 double
315 area,
316 fuzz;
317
318 MagickBooleanType
319 status = MagickTrue;
320
321 size_t
322 columns,
323 rows;
324
325 ssize_t
326 k,
327 y;
328
329 /*
330 Compute the absolute error similarity.
331 */
332 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
333 (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
334 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
335 image_view=AcquireVirtualCacheView(image,exception);
336 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
337#if defined(MAGICKCORE_OPENMP_SUPPORT)
338 #pragma omp parallel for schedule(static) shared(similarity,status) \
339 magick_number_threads(image,image,rows,1)
340#endif
341 for (y=0; y < (ssize_t) rows; y++)
342 {
343 const Quantum
344 *magick_restrict p,
345 *magick_restrict q;
346
347 double
348 channel_similarity[MaxPixelChannels+1] = { 0.0 };
349
350 ssize_t
351 x;
352
353 if (status == MagickFalse)
354 continue;
355 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
356 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
357 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
358 {
359 status=MagickFalse;
360 continue;
361 }
362 for (x=0; x < (ssize_t) columns; x++)
363 {
364 double
365 Da,
366 Sa;
367
368 size_t
369 count = 0;
370
371 ssize_t
372 i;
373
374 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
375 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
376 {
377 p+=(ptrdiff_t) GetPixelChannels(image);
378 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
379 continue;
380 }
381 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
382 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
383 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
384 {
385 double
386 error;
387
388 PixelChannel channel = GetPixelChannelChannel(image,i);
389 PixelTrait traits = GetPixelChannelTraits(image,channel);
390 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
391 channel);
392 if (((traits & UpdatePixelTrait) == 0) ||
393 ((reconstruct_traits & UpdatePixelTrait) == 0))
394 continue;
395 if (channel == AlphaPixelChannel)
396 error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
397 channel,q);
398 else
399 error=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
400 if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
401 {
402 channel_similarity[i]++;
403 count++;
404 }
405 }
406 if (count != 0)
407 channel_similarity[CompositePixelChannel]++;
408 p+=(ptrdiff_t) GetPixelChannels(image);
409 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
410 }
411#if defined(MAGICKCORE_OPENMP_SUPPORT)
412 #pragma omp critical (MagickCore_GetAESimilarity)
413#endif
414 {
415 ssize_t
416 j;
417
418 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
419 {
420 PixelChannel channel = GetPixelChannelChannel(image,j);
421 PixelTrait traits = GetPixelChannelTraits(image,channel);
422 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
423 channel);
424 if (((traits & UpdatePixelTrait) == 0) ||
425 ((reconstruct_traits & UpdatePixelTrait) == 0))
426 continue;
427 similarity[j]+=channel_similarity[j];
428 }
429 similarity[CompositePixelChannel]+=
430 channel_similarity[CompositePixelChannel];
431 }
432 }
433 reconstruct_view=DestroyCacheView(reconstruct_view);
434 image_view=DestroyCacheView(image_view);
435 area=MagickSafeReciprocal((double) columns*rows);
436 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
437 similarity[k]*=area;
438 similarity[CompositePixelChannel]*=area;
439 return(status);
440}
441
442static MagickBooleanType GetDPCSimilarity(const Image *image,
443 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
444{
445#define SimilarityImageTag "Similarity/Image"
446
447 CacheView
448 *image_view,
449 *reconstruct_view;
450
451 ChannelStatistics
452 *image_statistics,
453 *reconstruct_statistics;
454
455 double
456 norm[MaxPixelChannels+1] = { 0.0 },
457 reconstruct_norm[MaxPixelChannels+1] = { 0.0 };
458
459 MagickBooleanType
460 status = MagickTrue;
461
462 MagickOffsetType
463 progress = 0;
464
465 size_t
466 columns,
467 rows;
468
469 ssize_t
470 k,
471 y;
472
473 /*
474 Compute the dot product correlation similarity.
475 */
476 image_statistics=GetImageStatistics(image,exception);
477 reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
478 if ((image_statistics == (ChannelStatistics *) NULL) ||
479 (reconstruct_statistics == (ChannelStatistics *) NULL))
480 {
481 if (image_statistics != (ChannelStatistics *) NULL)
482 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
483 image_statistics);
484 if (reconstruct_statistics != (ChannelStatistics *) NULL)
485 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
486 reconstruct_statistics);
487 return(MagickFalse);
488 }
489 (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
490 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
491 image_view=AcquireVirtualCacheView(image,exception);
492 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
493#if defined(MAGICKCORE_OPENMP_SUPPORT)
494 #pragma omp parallel for schedule(static) shared(norm,reconstruct_norm,similarity,status) \
495 magick_number_threads(image,image,rows,1)
496#endif
497 for (y=0; y < (ssize_t) rows; y++)
498 {
499 const Quantum
500 *magick_restrict p,
501 *magick_restrict q;
502
503 double
504 channel_norm[MaxPixelChannels+1] = { 0.0 },
505 channel_reconstruct_norm[MaxPixelChannels+1] = { 0.0 },
506 channel_similarity[MaxPixelChannels+1] = { 0.0 };
507
508 ssize_t
509 x;
510
511 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
512 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
513 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
514 {
515 status=MagickFalse;
516 continue;
517 }
518 for (x=0; x < (ssize_t) columns; x++)
519 {
520 double
521 Da,
522 Sa;
523
524 ssize_t
525 i;
526
527 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
528 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
529 {
530 p+=(ptrdiff_t) GetPixelChannels(image);
531 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
532 continue;
533 }
534 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
535 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
536 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
537 {
538 double
539 alpha,
540 beta;
541
542 PixelChannel channel = GetPixelChannelChannel(image,i);
543 PixelTrait traits = GetPixelChannelTraits(image,channel);
544 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
545 channel);
546 if (((traits & UpdatePixelTrait) == 0) ||
547 ((reconstruct_traits & UpdatePixelTrait) == 0))
548 continue;
549 if (channel == AlphaPixelChannel)
550 {
551 alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
552 beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
553 channel,q)-reconstruct_statistics[channel].mean);
554 }
555 else
556 {
557 alpha=QuantumScale*(Sa*p[i]-image_statistics[channel].mean);
558 beta=QuantumScale*(Da*GetPixelChannel(reconstruct_image,channel,
559 q)-reconstruct_statistics[channel].mean);
560 }
561 channel_similarity[i]+=alpha*beta;
562 channel_norm[i]+=alpha*alpha;
563 channel_reconstruct_norm[i]+=beta*beta;
564 }
565 p+=(ptrdiff_t) GetPixelChannels(image);
566 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
567 }
568#if defined(MAGICKCORE_OPENMP_SUPPORT)
569 #pragma omp critical (MagickCore_GetDPCSimilarity)
570#endif
571 {
572 ssize_t
573 j;
574
575 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
576 {
577 PixelChannel channel = GetPixelChannelChannel(image,j);
578 PixelTrait traits = GetPixelChannelTraits(image,channel);
579 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
580 channel);
581 if (((traits & UpdatePixelTrait) == 0) ||
582 ((reconstruct_traits & UpdatePixelTrait) == 0))
583 continue;
584 similarity[j]+=channel_similarity[j];
585 similarity[CompositePixelChannel]+=channel_similarity[j];
586 norm[j]+=channel_norm[j];
587 norm[CompositePixelChannel]+=channel_norm[j];
588 reconstruct_norm[j]+=channel_reconstruct_norm[j];
589 reconstruct_norm[CompositePixelChannel]+=channel_reconstruct_norm[j];
590 }
591 }
592 if (image->progress_monitor != (MagickProgressMonitor) NULL)
593 {
594 MagickBooleanType
595 proceed;
596
597#if defined(MAGICKCORE_OPENMP_SUPPORT)
598 #pragma omp atomic
599#endif
600 progress++;
601 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
602 if (proceed == MagickFalse)
603 {
604 status=MagickFalse;
605 continue;
606 }
607 }
608 }
609 reconstruct_view=DestroyCacheView(reconstruct_view);
610 image_view=DestroyCacheView(image_view);
611 /*
612 Compute dot product correlation: divide by mean.
613 */
614 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
615 {
616 PixelChannel channel = GetPixelChannelChannel(image,k);
617 PixelTrait traits = GetPixelChannelTraits(image,channel);
618 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
619 channel);
620 if (((traits & UpdatePixelTrait) == 0) ||
621 ((reconstruct_traits & UpdatePixelTrait) == 0))
622 continue;
623 similarity[k]*=MagickSafeReciprocal(sqrt(norm[k]*reconstruct_norm[k]));
624 }
625 similarity[CompositePixelChannel]*=MagickSafeReciprocal(sqrt(
626 norm[CompositePixelChannel]*reconstruct_norm[CompositePixelChannel]));
627 /*
628 Free resources.
629 */
630 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
631 reconstruct_statistics);
632 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
633 image_statistics);
634 return(status);
635}
636
637static MagickBooleanType GetFUZZSimilarity(const Image *image,
638 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
639{
640 CacheView
641 *image_view,
642 *reconstruct_view;
643
644 double
645 area = 0.0,
646 fuzz = 0.0;
647
648 MagickBooleanType
649 status = MagickTrue;
650
651 size_t
652 columns,
653 rows;
654
655 ssize_t
656 k,
657 y;
658
659 /*
660 Compute the MSE similarity within tolerance (fuzz).
661 */
662 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
663 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
664 image_view=AcquireVirtualCacheView(image,exception);
665 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
666#if defined(MAGICKCORE_OPENMP_SUPPORT)
667 #pragma omp parallel for schedule(static) shared(area,similarity,status) \
668 magick_number_threads(image,image,rows,1)
669#endif
670 for (y=0; y < (ssize_t) rows; y++)
671 {
672 const Quantum
673 *magick_restrict p,
674 *magick_restrict q;
675
676 double
677 channel_area = 0.0,
678 channel_similarity[MaxPixelChannels+1] = { 0.0 };
679
680 ssize_t
681 x;
682
683 if (status == MagickFalse)
684 continue;
685 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
686 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
687 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
688 {
689 status=MagickFalse;
690 continue;
691 }
692 for (x=0; x < (ssize_t) columns; x++)
693 {
694 double
695 Da,
696 Sa;
697
698 ssize_t
699 i;
700
701 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
702 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
703 {
704 p+=(ptrdiff_t) GetPixelChannels(image);
705 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
706 continue;
707 }
708 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
709 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
710 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
711 {
712 double
713 error;
714
715 PixelChannel channel = GetPixelChannelChannel(image,i);
716 PixelTrait traits = GetPixelChannelTraits(image,channel);
717 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
718 channel);
719 if (((traits & UpdatePixelTrait) == 0) ||
720 ((reconstruct_traits & UpdatePixelTrait) == 0))
721 continue;
722 if (channel == AlphaPixelChannel)
723 error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
724 channel,q);
725 else
726 error=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
727 if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
728 {
729 channel_similarity[i]+=QuantumScale*error*QuantumScale*error;
730 channel_similarity[CompositePixelChannel]+=QuantumScale*error*
731 QuantumScale*error;
732 channel_area++;
733 }
734 }
735 p+=(ptrdiff_t) GetPixelChannels(image);
736 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
737 }
738#if defined(MAGICKCORE_OPENMP_SUPPORT)
739 #pragma omp critical (MagickCore_GetFUZZSimilarity)
740#endif
741 {
742 ssize_t
743 j;
744
745 area+=channel_area;
746 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
747 {
748 PixelChannel channel = GetPixelChannelChannel(image,j);
749 PixelTrait traits = GetPixelChannelTraits(image,channel);
750 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
751 channel);
752 if (((traits & UpdatePixelTrait) == 0) ||
753 ((reconstruct_traits & UpdatePixelTrait) == 0))
754 continue;
755 similarity[j]+=channel_similarity[j];
756 }
757 similarity[CompositePixelChannel]+=
758 channel_similarity[CompositePixelChannel];
759 }
760 }
761 reconstruct_view=DestroyCacheView(reconstruct_view);
762 image_view=DestroyCacheView(image_view);
763 area=MagickSafeReciprocal(area);
764 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
765 {
766 PixelChannel channel = GetPixelChannelChannel(image,k);
767 PixelTrait traits = GetPixelChannelTraits(image,channel);
768 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
769 channel);
770 if (((traits & UpdatePixelTrait) == 0) ||
771 ((reconstruct_traits & UpdatePixelTrait) == 0))
772 continue;
773 similarity[k]*=area;
774 }
775 similarity[CompositePixelChannel]*=area;
776 return(status);
777}
778
779static MagickBooleanType GetMAESimilarity(const Image *image,
780 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
781{
782 CacheView
783 *image_view,
784 *reconstruct_view;
785
786 double
787 area = 0.0;
788
789 MagickBooleanType
790 status = MagickTrue;
791
792 size_t
793 columns,
794 rows;
795
796 ssize_t
797 k,
798 y;
799
800 /*
801 Compute the mean absolute error similarity.
802 */
803 (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
804 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
805 image_view=AcquireVirtualCacheView(image,exception);
806 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
807#if defined(MAGICKCORE_OPENMP_SUPPORT)
808 #pragma omp parallel for schedule(static) shared(area,similarity,status) \
809 magick_number_threads(image,image,rows,1)
810#endif
811 for (y=0; y < (ssize_t) rows; y++)
812 {
813 const Quantum
814 *magick_restrict p,
815 *magick_restrict q;
816
817 double
818 channel_area = 0.0,
819 channel_similarity[MaxPixelChannels+1] = { 0.0 };
820
821 ssize_t
822 x;
823
824 if (status == MagickFalse)
825 continue;
826 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
827 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
828 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
829 {
830 status=MagickFalse;
831 continue;
832 }
833 for (x=0; x < (ssize_t) columns; x++)
834 {
835 double
836 Da,
837 Sa;
838
839 ssize_t
840 i;
841
842 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
843 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
844 {
845 p+=(ptrdiff_t) GetPixelChannels(image);
846 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
847 continue;
848 }
849 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
850 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
851 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
852 {
853 double
854 error;
855
856 PixelChannel channel = GetPixelChannelChannel(image,i);
857 PixelTrait traits = GetPixelChannelTraits(image,channel);
858 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
859 channel);
860 if (((traits & UpdatePixelTrait) == 0) ||
861 ((reconstruct_traits & UpdatePixelTrait) == 0))
862 continue;
863 if (channel == AlphaPixelChannel)
864 error=QuantumScale*fabs((double) p[i]-(double) GetPixelChannel(
865 reconstruct_image,channel,q));
866 else
867 error=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
868 channel,q));
869 channel_similarity[i]+=error;
870 channel_similarity[CompositePixelChannel]+=error;
871 }
872 channel_area++;
873 p+=(ptrdiff_t) GetPixelChannels(image);
874 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
875 }
876#if defined(MAGICKCORE_OPENMP_SUPPORT)
877 #pragma omp critical (MagickCore_GetMAESimilarity)
878#endif
879 {
880 ssize_t
881 j;
882
883 area+=channel_area;
884 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
885 {
886 PixelChannel channel = GetPixelChannelChannel(image,j);
887 PixelTrait traits = GetPixelChannelTraits(image,channel);
888 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
889 channel);
890 if (((traits & UpdatePixelTrait) == 0) ||
891 ((reconstruct_traits & UpdatePixelTrait) == 0))
892 continue;
893 similarity[j]+=channel_similarity[j];
894 }
895 similarity[CompositePixelChannel]+=
896 channel_similarity[CompositePixelChannel];
897 }
898 }
899 reconstruct_view=DestroyCacheView(reconstruct_view);
900 image_view=DestroyCacheView(image_view);
901 area=MagickSafeReciprocal(area);
902 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
903 {
904 PixelChannel channel = GetPixelChannelChannel(image,k);
905 PixelTrait traits = GetPixelChannelTraits(image,channel);
906 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
907 channel);
908 if (((traits & UpdatePixelTrait) == 0) ||
909 ((reconstruct_traits & UpdatePixelTrait) == 0))
910 continue;
911 similarity[k]*=area;
912 }
913 similarity[CompositePixelChannel]*=area;
914 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
915 return(status);
916}
917
918static MagickBooleanType GetMEPPSimilarity(Image *image,
919 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
920{
921 CacheView
922 *image_view,
923 *reconstruct_view;
924
925 double
926 area = 0.0,
927 maximum_error = -MagickMaximumValue,
928 mean_error = 0.0;
929
930 MagickBooleanType
931 status = MagickTrue;
932
933 size_t
934 columns,
935 rows;
936
937 ssize_t
938 k,
939 y;
940
941 /*
942 Compute the mean error per pixel similarity.
943 */
944 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
945 image_view=AcquireVirtualCacheView(image,exception);
946 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
947#if defined(MAGICKCORE_OPENMP_SUPPORT)
948 #pragma omp parallel for schedule(static) shared(area,similarity,maximum_error,mean_error,status) \
949 magick_number_threads(image,image,rows,1)
950#endif
951 for (y=0; y < (ssize_t) rows; y++)
952 {
953 const Quantum
954 *magick_restrict p,
955 *magick_restrict q;
956
957 double
958 channel_area = 0.0,
959 channel_similarity[MaxPixelChannels+1] = { 0.0 },
960 channel_maximum_error = maximum_error,
961 channel_mean_error = 0.0;
962
963 ssize_t
964 x;
965
966 if (status == MagickFalse)
967 continue;
968 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
969 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
970 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
971 {
972 status=MagickFalse;
973 continue;
974 }
975 for (x=0; x < (ssize_t) columns; x++)
976 {
977 double
978 Da,
979 Sa;
980
981 ssize_t
982 i;
983
984 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
985 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
986 {
987 p+=(ptrdiff_t) GetPixelChannels(image);
988 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
989 continue;
990 }
991 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
992 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
993 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
994 {
995 double
996 error;
997
998 PixelChannel channel = GetPixelChannelChannel(image,i);
999 PixelTrait traits = GetPixelChannelTraits(image,channel);
1000 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1001 channel);
1002 if (((traits & UpdatePixelTrait) == 0) ||
1003 ((reconstruct_traits & UpdatePixelTrait) == 0))
1004 continue;
1005 if (channel == AlphaPixelChannel)
1006 error=QuantumScale*fabs((double) p[i]-(double) GetPixelChannel(
1007 reconstruct_image,channel,q));
1008 else
1009 error=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
1010 channel,q));
1011 channel_similarity[i]+=error;
1012 channel_similarity[CompositePixelChannel]+=error;
1013 channel_mean_error+=error*error;
1014 if (error > channel_maximum_error)
1015 channel_maximum_error=error;
1016 }
1017 channel_area++;
1018 p+=(ptrdiff_t) GetPixelChannels(image);
1019 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1020 }
1021#if defined(MAGICKCORE_OPENMP_SUPPORT)
1022 #pragma omp critical (MagickCore_GetMEPPSimilarity)
1023#endif
1024 {
1025 ssize_t
1026 j;
1027
1028 area+=channel_area;
1029 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1030 {
1031 PixelChannel channel = GetPixelChannelChannel(image,j);
1032 PixelTrait traits = GetPixelChannelTraits(image,channel);
1033 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1034 channel);
1035 if (((traits & UpdatePixelTrait) == 0) ||
1036 ((reconstruct_traits & UpdatePixelTrait) == 0))
1037 continue;
1038 similarity[j]+=channel_similarity[j];
1039 }
1040 similarity[CompositePixelChannel]+=
1041 channel_similarity[CompositePixelChannel];
1042 mean_error+=channel_mean_error;
1043 if (channel_maximum_error > maximum_error)
1044 maximum_error=channel_maximum_error;
1045 }
1046 }
1047 reconstruct_view=DestroyCacheView(reconstruct_view);
1048 image_view=DestroyCacheView(image_view);
1049 area=MagickSafeReciprocal(area);
1050 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1051 {
1052 PixelChannel channel = GetPixelChannelChannel(image,k);
1053 PixelTrait traits = GetPixelChannelTraits(image,channel);
1054 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1055 channel);
1056 if (((traits & UpdatePixelTrait) == 0) ||
1057 ((reconstruct_traits & UpdatePixelTrait) == 0))
1058 continue;
1059 similarity[k]*=area;
1060 }
1061 similarity[CompositePixelChannel]*=area;
1062 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
1063 image->error.mean_error_per_pixel=QuantumRange*
1064 similarity[CompositePixelChannel];
1065 image->error.normalized_mean_error=mean_error*area;
1066 image->error.normalized_maximum_error=maximum_error;
1067 return(status);
1068}
1069
1070static MagickBooleanType GetMSESimilarity(const Image *image,
1071 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1072{
1073 CacheView
1074 *image_view,
1075 *reconstruct_view;
1076
1077 double
1078 area = 0.0;
1079
1080 MagickBooleanType
1081 status = MagickTrue;
1082
1083 size_t
1084 columns,
1085 rows;
1086
1087 ssize_t
1088 k,
1089 y;
1090
1091 /*
1092 Compute the mean sequared error similarity.
1093 */
1094 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1095 image_view=AcquireVirtualCacheView(image,exception);
1096 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1097#if defined(MAGICKCORE_OPENMP_SUPPORT)
1098 #pragma omp parallel for schedule(static) shared(area,similarity,status) \
1099 magick_number_threads(image,image,rows,1)
1100#endif
1101 for (y=0; y < (ssize_t) rows; y++)
1102 {
1103 const Quantum
1104 *magick_restrict p,
1105 *magick_restrict q;
1106
1107 double
1108 channel_area = 0.0,
1109 channel_similarity[MaxPixelChannels+1] = { 0.0 };
1110
1111 ssize_t
1112 x;
1113
1114 if (status == MagickFalse)
1115 continue;
1116 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1117 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1118 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1119 {
1120 status=MagickFalse;
1121 continue;
1122 }
1123 for (x=0; x < (ssize_t) columns; x++)
1124 {
1125 double
1126 Da,
1127 Sa;
1128
1129 ssize_t
1130 i;
1131
1132 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1133 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1134 {
1135 p+=(ptrdiff_t) GetPixelChannels(image);
1136 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1137 continue;
1138 }
1139 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1140 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1141 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1142 {
1143 double
1144 error;
1145
1146 PixelChannel channel = GetPixelChannelChannel(image,i);
1147 PixelTrait traits = GetPixelChannelTraits(image,channel);
1148 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1149 channel);
1150 if (((traits & UpdatePixelTrait) == 0) ||
1151 ((reconstruct_traits & UpdatePixelTrait) == 0))
1152 continue;
1153 if (channel == AlphaPixelChannel)
1154 error=QuantumScale*((double) p[i]-(double) GetPixelChannel(
1155 reconstruct_image,channel,q));
1156 else
1157 error=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
1158 channel,q));
1159 channel_similarity[i]+=error*error;
1160 channel_similarity[CompositePixelChannel]+=error*error;
1161 }
1162 channel_area++;
1163 p+=(ptrdiff_t) GetPixelChannels(image);
1164 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1165 }
1166#if defined(MAGICKCORE_OPENMP_SUPPORT)
1167 #pragma omp critical (MagickCore_GetMSESimilarity)
1168#endif
1169 {
1170 ssize_t
1171 j;
1172
1173 area+=channel_area;
1174 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1175 {
1176 PixelChannel channel = GetPixelChannelChannel(image,j);
1177 PixelTrait traits = GetPixelChannelTraits(image,channel);
1178 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1179 channel);
1180 if (((traits & UpdatePixelTrait) == 0) ||
1181 ((reconstruct_traits & UpdatePixelTrait) == 0))
1182 continue;
1183 similarity[j]+=channel_similarity[j];
1184 }
1185 similarity[CompositePixelChannel]+=
1186 channel_similarity[CompositePixelChannel];
1187 }
1188 }
1189 reconstruct_view=DestroyCacheView(reconstruct_view);
1190 image_view=DestroyCacheView(image_view);
1191 area=MagickSafeReciprocal(area);
1192 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1193 {
1194 PixelChannel channel = GetPixelChannelChannel(image,k);
1195 PixelTrait traits = GetPixelChannelTraits(image,channel);
1196 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1197 channel);
1198 if (((traits & UpdatePixelTrait) == 0) ||
1199 ((reconstruct_traits & UpdatePixelTrait) == 0))
1200 continue;
1201 similarity[k]*=area;
1202 }
1203 similarity[CompositePixelChannel]*=area;
1204 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
1205 return(status);
1206}
1207
1208static MagickBooleanType GetNCCSimilarity(const Image *image,
1209 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1210{
1211 CacheView
1212 *image_view,
1213 *reconstruct_view;
1214
1215 ChannelStatistics
1216 *image_statistics,
1217 *reconstruct_statistics;
1218
1219 double
1220 reconstruct_variance[MaxPixelChannels+1] = { 0.0 },
1221 variance[MaxPixelChannels+1] = { 0.0 };
1222
1223 MagickBooleanType
1224 status = MagickTrue;
1225
1226 MagickOffsetType
1227 progress = 0;
1228
1229 size_t
1230 columns,
1231 rows;
1232
1233 ssize_t
1234 k,
1235 y;
1236
1237 /*
1238 Compute the normalized criss-correlation similarity.
1239 */
1240 image_statistics=GetImageStatistics(image,exception);
1241 reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
1242 if ((image_statistics == (ChannelStatistics *) NULL) ||
1243 (reconstruct_statistics == (ChannelStatistics *) NULL))
1244 {
1245 if (image_statistics != (ChannelStatistics *) NULL)
1246 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1247 image_statistics);
1248 if (reconstruct_statistics != (ChannelStatistics *) NULL)
1249 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1250 reconstruct_statistics);
1251 return(MagickFalse);
1252 }
1253 (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
1254 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1255 image_view=AcquireVirtualCacheView(image,exception);
1256 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1257#if defined(MAGICKCORE_OPENMP_SUPPORT)
1258 #pragma omp parallel for schedule(static) shared(variance,reconstruct_variance,similarity,status) \
1259 magick_number_threads(image,image,rows,1)
1260#endif
1261 for (y=0; y < (ssize_t) rows; y++)
1262 {
1263 const Quantum
1264 *magick_restrict p,
1265 *magick_restrict q;
1266
1267 double
1268 channel_reconstruct_variance[MaxPixelChannels+1] = { 0.0 },
1269 channel_similarity[MaxPixelChannels+1] = { 0.0 },
1270 channel_variance[MaxPixelChannels+1] = { 0.0 };
1271
1272 ssize_t
1273 x;
1274
1275 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1276 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1277 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1278 {
1279 status=MagickFalse;
1280 continue;
1281 }
1282 for (x=0; x < (ssize_t) columns; x++)
1283 {
1284 double
1285 Da,
1286 Sa;
1287
1288 ssize_t
1289 i;
1290
1291 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1292 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1293 {
1294 p+=(ptrdiff_t) GetPixelChannels(image);
1295 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1296 continue;
1297 }
1298 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1299 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1300 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1301 {
1302 double
1303 alpha,
1304 beta;
1305
1306 PixelChannel channel = GetPixelChannelChannel(image,i);
1307 PixelTrait traits = GetPixelChannelTraits(image,channel);
1308 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1309 channel);
1310 if (((traits & UpdatePixelTrait) == 0) ||
1311 ((reconstruct_traits & UpdatePixelTrait) == 0))
1312 continue;
1313 if (channel == AlphaPixelChannel)
1314 {
1315 alpha=QuantumScale*((double) p[i]-image_statistics[channel].mean);
1316 beta=QuantumScale*((double) GetPixelChannel(reconstruct_image,
1317 channel,q)-reconstruct_statistics[channel].mean);
1318 }
1319 else
1320 {
1321 alpha=QuantumScale*(Sa*p[i]-image_statistics[channel].mean);
1322 beta=QuantumScale*(Da*GetPixelChannel(reconstruct_image,channel,
1323 q)-reconstruct_statistics[channel].mean);
1324 }
1325 channel_similarity[i]+=alpha*beta;
1326 channel_variance[i]+=alpha*alpha;
1327 channel_reconstruct_variance[i]+=beta*beta;
1328 }
1329 p+=(ptrdiff_t) GetPixelChannels(image);
1330 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1331 }
1332#if defined(MAGICKCORE_OPENMP_SUPPORT)
1333 #pragma omp critical (MagickCore_GetNCCSimilarity)
1334#endif
1335 {
1336 ssize_t
1337 j;
1338
1339 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1340 {
1341 PixelChannel channel = GetPixelChannelChannel(image,j);
1342 PixelTrait traits = GetPixelChannelTraits(image,channel);
1343 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1344 channel);
1345 if (((traits & UpdatePixelTrait) == 0) ||
1346 ((reconstruct_traits & UpdatePixelTrait) == 0))
1347 continue;
1348 similarity[j]+=channel_similarity[j];
1349 similarity[CompositePixelChannel]+=channel_similarity[j];
1350 variance[j]+=channel_variance[j];
1351 variance[CompositePixelChannel]+=channel_variance[j];
1352 reconstruct_variance[j]+=channel_reconstruct_variance[j];
1353 reconstruct_variance[CompositePixelChannel]+=
1354 channel_reconstruct_variance[j];
1355 }
1356 }
1357 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1358 {
1359 MagickBooleanType
1360 proceed;
1361
1362#if defined(MAGICKCORE_OPENMP_SUPPORT)
1363 #pragma omp atomic
1364#endif
1365 progress++;
1366 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1367 if (proceed == MagickFalse)
1368 {
1369 status=MagickFalse;
1370 continue;
1371 }
1372 }
1373 }
1374 reconstruct_view=DestroyCacheView(reconstruct_view);
1375 image_view=DestroyCacheView(image_view);
1376 /*
1377 Compute normalized cross correlation: divide by standard deviation.
1378 */
1379 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1380 {
1381 PixelChannel channel = GetPixelChannelChannel(image,k);
1382 PixelTrait traits = GetPixelChannelTraits(image,channel);
1383 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1384 channel);
1385 if (((traits & UpdatePixelTrait) == 0) ||
1386 ((reconstruct_traits & UpdatePixelTrait) == 0))
1387 continue;
1388 similarity[k]*=MagickSafeReciprocal(sqrt(variance[k])*
1389 sqrt(reconstruct_variance[k]));
1390 }
1391 similarity[CompositePixelChannel]*=MagickSafeReciprocal(sqrt(
1392 variance[CompositePixelChannel])*sqrt(
1393 reconstruct_variance[CompositePixelChannel]));
1394 /*
1395 Free resources.
1396 */
1397 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1398 reconstruct_statistics);
1399 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1400 image_statistics);
1401 return(status);
1402}
1403
1404static MagickBooleanType GetPASimilarity(const Image *image,
1405 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1406{
1407 CacheView
1408 *image_view,
1409 *reconstruct_view;
1410
1411 MagickBooleanType
1412 status = MagickTrue;
1413
1414 size_t
1415 columns,
1416 rows;
1417
1418 ssize_t
1419 y;
1420
1421 /*
1422 Compute the peak absolute similarity.
1423 */
1424 (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
1425 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1426 image_view=AcquireVirtualCacheView(image,exception);
1427 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1428#if defined(MAGICKCORE_OPENMP_SUPPORT)
1429 #pragma omp parallel for schedule(static) shared(similarity,status) \
1430 magick_number_threads(image,image,rows,1)
1431#endif
1432 for (y=0; y < (ssize_t) rows; y++)
1433 {
1434 const Quantum
1435 *magick_restrict p,
1436 *magick_restrict q;
1437
1438 double
1439 channel_similarity[MaxPixelChannels+1] = { 0.0 };
1440
1441 ssize_t
1442 x;
1443
1444 if (status == MagickFalse)
1445 continue;
1446 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1447 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1448 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1449 {
1450 status=MagickFalse;
1451 continue;
1452 }
1453 for (x=0; x < (ssize_t) columns; x++)
1454 {
1455 double
1456 Da,
1457 Sa;
1458
1459 ssize_t
1460 i;
1461
1462 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1463 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1464 {
1465 p+=(ptrdiff_t) GetPixelChannels(image);
1466 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1467 continue;
1468 }
1469 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1470 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1471 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1472 {
1473 double
1474 distance;
1475
1476 PixelChannel channel = GetPixelChannelChannel(image,i);
1477 PixelTrait traits = GetPixelChannelTraits(image,channel);
1478 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1479 channel);
1480 if (((traits & UpdatePixelTrait) == 0) ||
1481 ((reconstruct_traits & UpdatePixelTrait) == 0))
1482 continue;
1483 if (channel == AlphaPixelChannel)
1484 distance=QuantumScale*fabs((double) p[i]-(double)
1485 GetPixelChannel(reconstruct_image,channel,q));
1486 else
1487 distance=QuantumScale*fabs(Sa*p[i]-Da*GetPixelChannel(
1488 reconstruct_image,channel,q));
1489 if (distance > channel_similarity[i])
1490 channel_similarity[i]=distance;
1491 if (distance > channel_similarity[CompositePixelChannel])
1492 channel_similarity[CompositePixelChannel]=distance;
1493 }
1494 p+=(ptrdiff_t) GetPixelChannels(image);
1495 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1496 }
1497#if defined(MAGICKCORE_OPENMP_SUPPORT)
1498 #pragma omp critical (MagickCore_GetPASimilarity)
1499#endif
1500 {
1501 ssize_t
1502 j;
1503
1504 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1505 {
1506 PixelChannel channel = GetPixelChannelChannel(image,j);
1507 PixelTrait traits = GetPixelChannelTraits(image,channel);
1508 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1509 channel);
1510 if (((traits & UpdatePixelTrait) == 0) ||
1511 ((reconstruct_traits & UpdatePixelTrait) == 0))
1512 continue;
1513 if (channel_similarity[j] > similarity[j])
1514 similarity[j]=channel_similarity[j];
1515 }
1516 if (channel_similarity[CompositePixelChannel] > similarity[CompositePixelChannel])
1517 similarity[CompositePixelChannel]=
1518 channel_similarity[CompositePixelChannel];
1519 }
1520 }
1521 reconstruct_view=DestroyCacheView(reconstruct_view);
1522 image_view=DestroyCacheView(image_view);
1523 return(status);
1524}
1525
1526static MagickBooleanType GetPDCSimilarity(const Image *image,
1527 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1528{
1529 CacheView
1530 *image_view,
1531 *reconstruct_view;
1532
1533 double
1534 area,
1535 fuzz;
1536
1537 MagickBooleanType
1538 status = MagickTrue;
1539
1540 size_t
1541 columns,
1542 rows;
1543
1544 ssize_t
1545 k,
1546 y;
1547
1548 /*
1549 Compute the pixel difference count similarity.
1550 */
1551 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
1552 (void) memset(similarity,0,(MaxPixelChannels+1)*sizeof(*similarity));
1553 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1554 image_view=AcquireVirtualCacheView(image,exception);
1555 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1556#if defined(MMAGICKCORE_OPENMP_SUPPORT)
1557 #pragma omp parallel for schedule(static) shared(similarity,status) \
1558 magick_number_threads(image,image,rows,1)
1559#endif
1560 for (y=0; y < (ssize_t) rows; y++)
1561 {
1562 const Quantum
1563 *magick_restrict p,
1564 *magick_restrict q;
1565
1566 double
1567 channel_similarity[MaxPixelChannels+1] = { 0.0 };
1568
1569 ssize_t
1570 x;
1571
1572 if (status == MagickFalse)
1573 continue;
1574 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1575 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1576 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1577 {
1578 status=MagickFalse;
1579 continue;
1580 }
1581 for (x=0; x < (ssize_t) columns; x++)
1582 {
1583 double
1584 Da,
1585 Sa;
1586
1587 size_t
1588 count = 0;
1589
1590 ssize_t
1591 i;
1592
1593 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1594 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1595 {
1596 p+=(ptrdiff_t) GetPixelChannels(image);
1597 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1598 continue;
1599 }
1600 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1601 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1602 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1603 {
1604 double
1605 error;
1606
1607 PixelChannel channel = GetPixelChannelChannel(image,i);
1608 PixelTrait traits = GetPixelChannelTraits(image,channel);
1609 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1610 channel);
1611 if (((traits & UpdatePixelTrait) == 0) ||
1612 ((reconstruct_traits & UpdatePixelTrait) == 0))
1613 continue;
1614 if (channel == AlphaPixelChannel)
1615 error=(double) p[i]-(double) GetPixelChannel(reconstruct_image,
1616 channel,q);
1617 else
1618 error=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
1619 if (MagickSafeSignificantError(error*error,fuzz) != MagickFalse)
1620 {
1621 channel_similarity[i]++;
1622 count++;
1623 }
1624 }
1625 if (count != 0)
1626 channel_similarity[CompositePixelChannel]++;
1627 p+=(ptrdiff_t) GetPixelChannels(image);
1628 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1629 }
1630#if defined(MAGICKCORE_OPENMP_SUPPORT)
1631 #pragma omp critical (MagickCore_GetPDCSimilarity)
1632#endif
1633 {
1634 ssize_t
1635 j;
1636
1637 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1638 {
1639 PixelChannel channel = GetPixelChannelChannel(image,j);
1640 PixelTrait traits = GetPixelChannelTraits(image,channel);
1641 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1642 channel);
1643 if (((traits & UpdatePixelTrait) == 0) ||
1644 ((reconstruct_traits & UpdatePixelTrait) == 0))
1645 continue;
1646 similarity[j]+=channel_similarity[j];
1647 }
1648 similarity[CompositePixelChannel]+=
1649 channel_similarity[CompositePixelChannel];
1650 }
1651 }
1652 reconstruct_view=DestroyCacheView(reconstruct_view);
1653 image_view=DestroyCacheView(image_view);
1654 area=MagickSafeReciprocal((double) columns*rows);
1655 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1656 similarity[k]*=area;
1657 similarity[CompositePixelChannel]*=area;
1658 return(status);
1659}
1660
1661static MagickBooleanType DFTPhaseSpectrum(const Image *image,const ssize_t u,
1662 const ssize_t v,double *phase,ExceptionInfo *exception)
1663{
1664#define PhaseImageTag "Phase/Image"
1665
1666 CacheView
1667 *image_view;
1668
1669 double
1670 channel_imag[MaxPixelChannels+1] = { 0.0 },
1671 channel_real[MaxPixelChannels+1] = { 0.0 };
1672
1673 MagickBooleanType
1674 status;
1675
1676 ssize_t
1677 k,
1678 y;
1679
1680 /*
1681 Compute DFT phase spectrum of an image.
1682 */
1683 status=MagickTrue;
1684 image_view=AcquireVirtualCacheView(image,exception);
1685 for (y=0; y < (ssize_t) image->rows; y++)
1686 {
1687 const Quantum
1688 *magick_restrict p;
1689
1690 ssize_t
1691 x;
1692
1693 if (status == MagickFalse)
1694 continue;
1695 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1696 if (p == (const Quantum *) NULL)
1697 {
1698 status=MagickFalse;
1699 continue;
1700 }
1701 for (x=0; x < (ssize_t) image->columns; x++)
1702 {
1703 double
1704 angle,
1705 Sa;
1706
1707 ssize_t
1708 i;
1709
1710 angle=MagickPI*((u*x/(double) image->rows)+(v*y/(double) image->columns));
1711 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1712 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1713 {
1714 PixelChannel channel = GetPixelChannelChannel(image,i);
1715 PixelTrait traits = GetPixelChannelTraits(image,channel);
1716 if (traits == UndefinedPixelTrait)
1717 continue;
1718 if (channel == AlphaPixelChannel)
1719 {
1720 channel_real[i]+=(QuantumScale*p[i])*cos(angle);
1721 channel_imag[i]-=(QuantumScale*p[i])*sin(angle);
1722 }
1723 else
1724 {
1725 channel_real[i]+=(QuantumScale*Sa*p[i])*cos(angle);
1726 channel_imag[i]-=(QuantumScale*Sa*p[i])*sin(angle);
1727 }
1728 }
1729 p+=(ptrdiff_t) GetPixelChannels(image);
1730 }
1731 }
1732 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1733 phase[k]=atan2(channel_imag[k],channel_real[k]);
1734 phase[CompositePixelChannel]=atan2(channel_imag[CompositePixelChannel],
1735 channel_real[CompositePixelChannel]);
1736 image_view=DestroyCacheView(image_view);
1737 return(status);
1738}
1739
1740static MagickBooleanType GetPHASESimilarity(const Image *image,
1741 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1742{
1743 CacheView
1744 *image_view,
1745 *reconstruct_view;
1746
1747 double
1748 area = 0.0;
1749
1750 MagickBooleanType
1751 status = MagickTrue;
1752
1753 size_t
1754 columns,
1755 rows;
1756
1757 ssize_t
1758 k,
1759 y;
1760
1761 /*
1762 Compute the phase congruency similarity.
1763 */
1764 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
1765 image_view=AcquireVirtualCacheView(image,exception);
1766 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1767#if defined(MAGICKCORE_OPENMP_SUPPORT)
1768 #pragma omp parallel for schedule(static) shared(area,similarity,status) \
1769 magick_number_threads(image,image,rows,1)
1770#endif
1771 for (y=0; y < (ssize_t) rows; y++)
1772 {
1773 const Quantum
1774 *magick_restrict p,
1775 *magick_restrict q;
1776
1777 double
1778 channel_area = 0.0,
1779 channel_similarity[MaxPixelChannels+1] = { 0.0 };
1780
1781 ssize_t
1782 x;
1783
1784 if (status == MagickFalse)
1785 continue;
1786 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1787 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1788 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1789 {
1790 status=MagickFalse;
1791 continue;
1792 }
1793 for (x=0; x < (ssize_t) columns; x++)
1794 {
1795 double
1796 phase[MaxPixelChannels+1] = { 0.0 },
1797 reconstruct_phase[MaxPixelChannels+1] = { 0.0 };
1798
1799 ssize_t
1800 i;
1801
1802 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1803 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1804 {
1805 p+=(ptrdiff_t) GetPixelChannels(image);
1806 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1807 continue;
1808 }
1809 status=DFTPhaseSpectrum(image,x,y,phase,exception);
1810 if (status == MagickFalse)
1811 break;
1812 status=DFTPhaseSpectrum(reconstruct_image,x,y,reconstruct_phase,
1813 exception);
1814 if (status == MagickFalse)
1815 break;
1816 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1817 {
1818 double
1819 delta;
1820
1821 PixelChannel channel = GetPixelChannelChannel(image,i);
1822 PixelTrait traits = GetPixelChannelTraits(image,channel);
1823 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1824 channel);
1825 if (((traits & UpdatePixelTrait) == 0) ||
1826 ((reconstruct_traits & UpdatePixelTrait) == 0))
1827 continue;
1828 delta=phase[i]-reconstruct_phase[i];
1829 channel_similarity[i]+=cos(delta);
1830 channel_similarity[CompositePixelChannel]+=cos(delta);
1831 }
1832 channel_area++;
1833 p+=(ptrdiff_t) GetPixelChannels(image);
1834 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1835 }
1836#if defined(MAGICKCORE_OPENMP_SUPPORT)
1837 #pragma omp critical (MagickCore_GetPHASESimilarity)
1838#endif
1839 {
1840 ssize_t
1841 j;
1842
1843 area+=channel_area;
1844 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1845 {
1846 PixelChannel channel = GetPixelChannelChannel(image,j);
1847 PixelTrait traits = GetPixelChannelTraits(image,channel);
1848 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1849 channel);
1850 if (((traits & UpdatePixelTrait) == 0) ||
1851 ((reconstruct_traits & UpdatePixelTrait) == 0))
1852 continue;
1853 similarity[j]+=channel_similarity[j];
1854 }
1855 similarity[CompositePixelChannel]+=
1856 channel_similarity[CompositePixelChannel];
1857 }
1858 }
1859 reconstruct_view=DestroyCacheView(reconstruct_view);
1860 image_view=DestroyCacheView(image_view);
1861 area=MagickSafeReciprocal(area);
1862 for (k=0; k < (ssize_t) GetPixelChannels(image); k++)
1863 {
1864 PixelChannel channel = GetPixelChannelChannel(image,k);
1865 PixelTrait traits = GetPixelChannelTraits(image,channel);
1866 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1867 channel);
1868 if (((traits & UpdatePixelTrait) == 0) ||
1869 ((reconstruct_traits & UpdatePixelTrait) == 0))
1870 continue;
1871 similarity[k]*=area;
1872 }
1873 similarity[CompositePixelChannel]*=area;
1874 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
1875 return(status);
1876}
1877
1878static MagickBooleanType GetPSNRSimilarity(const Image *image,
1879 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1880{
1881 MagickBooleanType
1882 status = MagickTrue;
1883
1884 ssize_t
1885 i;
1886
1887 /*
1888 Compute the peak signal-to-noise ratio similarity.
1889 */
1890 status=GetMSESimilarity(image,reconstruct_image,similarity,exception);
1891 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1892 {
1893 PixelChannel channel = GetPixelChannelChannel(image,i);
1894 PixelTrait traits = GetPixelChannelTraits(image,channel);
1895 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1896 channel);
1897 if (((traits & UpdatePixelTrait) == 0) ||
1898 ((reconstruct_traits & UpdatePixelTrait) == 0))
1899 continue;
1900 similarity[i]=10.0*MagickSafeLog10(MagickSafeReciprocal(
1901 similarity[i]))/MagickSafePSNRRecipicol(10.0);
1902 }
1903 similarity[CompositePixelChannel]=10.0*MagickSafeLog10(
1904 MagickSafeReciprocal(similarity[CompositePixelChannel]))/
1905 MagickSafePSNRRecipicol(10.0);
1906 return(status);
1907}
1908
1909static MagickBooleanType GetPHASHSimilarity(const Image *image,
1910 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
1911{
1912 ChannelPerceptualHash
1913 *channel_phash,
1914 *reconstruct_phash;
1915
1916 const char
1917 *artifact;
1918
1919 ssize_t
1920 i;
1921
1922 /*
1923 Compute the perceptual hash similarity.
1924 */
1925 channel_phash=GetImagePerceptualHash(image,exception);
1926 if (channel_phash == (ChannelPerceptualHash *) NULL)
1927 return(MagickFalse);
1928 reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1929 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1930 {
1931 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1932 channel_phash);
1933 return(MagickFalse);
1934 }
1935 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1936 {
1937 double
1938 difference = 0.0;
1939
1940 ssize_t
1941 j;
1942
1943 PixelChannel channel = GetPixelChannelChannel(image,i);
1944 PixelTrait traits = GetPixelChannelTraits(image,channel);
1945 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1946 channel);
1947 if (((traits & UpdatePixelTrait) == 0) ||
1948 ((reconstruct_traits & UpdatePixelTrait) == 0))
1949 continue;
1950 for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1951 {
1952 double
1953 alpha,
1954 beta;
1955
1956 ssize_t
1957 k;
1958
1959 for (k=0; k < MaximumNumberOfPerceptualHashes; k++)
1960 {
1961 double
1962 error;
1963
1964 alpha=channel_phash[i].phash[j][k];
1965 beta=reconstruct_phash[i].phash[j][k];
1966 error=beta-alpha;
1967 if (IsNaN(error) != 0)
1968 error=0.0;
1969 difference+=error*error;
1970 }
1971 }
1972 similarity[i]+=difference;
1973 similarity[CompositePixelChannel]+=difference;
1974 }
1975 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
1976 artifact=GetImageArtifact(image,"phash:normalize");
1977 if (IsStringTrue(artifact) != MagickFalse)
1978 {
1979 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1980 {
1981 PixelChannel channel = GetPixelChannelChannel(image,i);
1982 PixelTrait traits = GetPixelChannelTraits(image,channel);
1983 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1984 channel);
1985 if (((traits & UpdatePixelTrait) == 0) ||
1986 ((reconstruct_traits & UpdatePixelTrait) == 0))
1987 continue;
1988 similarity[i]=sqrt(similarity[i]/channel_phash[0].number_colorspaces);
1989 }
1990 similarity[CompositePixelChannel]=sqrt(similarity[CompositePixelChannel]/
1991 channel_phash[0].number_colorspaces);
1992 }
1993 /*
1994 Free resources.
1995 */
1996 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1997 reconstruct_phash);
1998 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1999 return(MagickTrue);
2000}
2001
2002static MagickBooleanType GetRMSESimilarity(const Image *image,
2003 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2004{
2005#define RMSESquareRoot(x) sqrt((x) < 0.0 ? 0.0 : (x))
2006
2007 MagickBooleanType
2008 status = MagickTrue;
2009
2010 ssize_t
2011 i;
2012
2013 /*
2014 Compute the root mean-squared error similarity.
2015 */
2016 status=GetMSESimilarity(image,reconstruct_image,similarity,exception);
2017 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2018 {
2019 PixelChannel channel = GetPixelChannelChannel(image,i);
2020 PixelTrait traits = GetPixelChannelTraits(image,channel);
2021 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2022 channel);
2023 if (((traits & UpdatePixelTrait) == 0) ||
2024 ((reconstruct_traits & UpdatePixelTrait) == 0))
2025 continue;
2026 similarity[i]=RMSESquareRoot(similarity[i]);
2027 }
2028 similarity[CompositePixelChannel]=RMSESquareRoot(
2029 similarity[CompositePixelChannel]);
2030 return(status);
2031}
2032
2033static MagickBooleanType GetSSIMSimularity(const Image *image,
2034 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2035{
2036#define SSIMRadius 5.0
2037#define SSIMSigma 1.5
2038#define SSIMK1 0.01
2039#define SSIMK2 0.03
2040#define SSIML 1.0
2041
2042 CacheView
2043 *image_view,
2044 *reconstruct_view;
2045
2046 char
2047 geometry[MagickPathExtent];
2048
2049 const char
2050 *artifact;
2051
2052 double
2053 area = 0.0,
2054 c1,
2055 c2,
2056 radius,
2057 sigma;
2058
2059 KernelInfo
2060 *kernel_info;
2061
2062 MagickBooleanType
2063 status = MagickTrue;
2064
2065 size_t
2066 columns,
2067 rows;
2068
2069 ssize_t
2070 l,
2071 y;
2072
2073 /*
2074 Compute the structual similarity index similarity.
2075 */
2076 radius=SSIMRadius;
2077 artifact=GetImageArtifact(image,"compare:ssim-radius");
2078 if (artifact != (const char *) NULL)
2079 radius=StringToDouble(artifact,(char **) NULL);
2080 sigma=SSIMSigma;
2081 artifact=GetImageArtifact(image,"compare:ssim-sigma");
2082 if (artifact != (const char *) NULL)
2083 sigma=StringToDouble(artifact,(char **) NULL);
2084 (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
2085 radius,sigma);
2086 kernel_info=AcquireKernelInfo(geometry,exception);
2087 if (kernel_info == (KernelInfo *) NULL)
2088 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
2089 image->filename);
2090 c1=pow(SSIMK1*SSIML,2.0);
2091 artifact=GetImageArtifact(image,"compare:ssim-k1");
2092 if (artifact != (const char *) NULL)
2093 c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
2094 c2=pow(SSIMK2*SSIML,2.0);
2095 artifact=GetImageArtifact(image,"compare:ssim-k2");
2096 if (artifact != (const char *) NULL)
2097 c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
2098 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
2099 image_view=AcquireVirtualCacheView(image,exception);
2100 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2101#if defined(MAGICKCORE_OPENMP_SUPPORT)
2102 #pragma omp parallel for schedule(static) shared(area,similarity,status) \
2103 magick_number_threads(image,reconstruct_image,rows,1)
2104#endif
2105 for (y=0; y < (ssize_t) rows; y++)
2106 {
2107 const Quantum
2108 *magick_restrict p,
2109 *magick_restrict q;
2110
2111 double
2112 channel_area = 0.0,
2113 channel_similarity[MaxPixelChannels+1] = { 0.0 };
2114
2115 ssize_t
2116 i,
2117 x;
2118
2119 if (status == MagickFalse)
2120 continue;
2121 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
2122 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
2123 kernel_info->height,exception);
2124 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
2125 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
2126 kernel_info->height,exception);
2127 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
2128 {
2129 status=MagickFalse;
2130 continue;
2131 }
2132 for (x=0; x < (ssize_t) columns; x++)
2133 {
2134 const Quantum
2135 *magick_restrict reconstruct,
2136 *magick_restrict test;
2137
2138 double
2139 x_pixel_mu[MaxPixelChannels+1] = { 0.0 },
2140 x_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 },
2141 xy_sigma[MaxPixelChannels+1] = { 0.0 },
2142 y_pixel_mu[MaxPixelChannels+1] = { 0.0 },
2143 y_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 };
2144
2145 MagickRealType
2146 *k;
2147
2148 ssize_t
2149 v;
2150
2151 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
2152 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
2153 {
2154 p+=(ptrdiff_t) GetPixelChannels(image);
2155 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2156 continue;
2157 }
2158 k=kernel_info->values;
2159 test=p;
2160 reconstruct=q;
2161 for (v=0; v < (ssize_t) kernel_info->height; v++)
2162 {
2163 ssize_t
2164 u;
2165
2166 for (u=0; u < (ssize_t) kernel_info->width; u++)
2167 {
2168 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2169 {
2170 double
2171 x_pixel,
2172 y_pixel;
2173
2174 PixelChannel channel = GetPixelChannelChannel(image,i);
2175 PixelTrait traits = GetPixelChannelTraits(image,channel);
2176 PixelTrait reconstruct_traits = GetPixelChannelTraits(
2177 reconstruct_image,channel);
2178 if (((traits & UpdatePixelTrait) == 0) ||
2179 ((reconstruct_traits & UpdatePixelTrait) == 0))
2180 continue;
2181 x_pixel=QuantumScale*(double) test[i];
2182 x_pixel_mu[i]+=(*k)*x_pixel;
2183 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
2184 y_pixel=QuantumScale*(double)
2185 GetPixelChannel(reconstruct_image,channel,reconstruct);
2186 y_pixel_mu[i]+=(*k)*y_pixel;
2187 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
2188 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
2189 }
2190 k++;
2191 test+=(ptrdiff_t) GetPixelChannels(image);
2192 reconstruct+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2193 }
2194 test+=(ptrdiff_t) GetPixelChannels(image)*columns;
2195 reconstruct+=(ptrdiff_t) GetPixelChannels(reconstruct_image)*columns;
2196 }
2197 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2198 {
2199 double
2200 ssim,
2201 x_pixel_mu_squared,
2202 x_pixel_sigmas_squared,
2203 xy_mu,
2204 xy_sigmas,
2205 y_pixel_mu_squared,
2206 y_pixel_sigmas_squared;
2207
2208 PixelChannel channel = GetPixelChannelChannel(image,i);
2209 PixelTrait traits = GetPixelChannelTraits(image,channel);
2210 PixelTrait reconstruct_traits = GetPixelChannelTraits(
2211 reconstruct_image,channel);
2212 if (((traits & UpdatePixelTrait) == 0) ||
2213 ((reconstruct_traits & UpdatePixelTrait) == 0))
2214 continue;
2215 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
2216 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
2217 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
2218 xy_sigmas=xy_sigma[i]-xy_mu;
2219 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
2220 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
2221 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))*
2222 MagickSafeReciprocal((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
2223 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
2224 channel_similarity[i]+=ssim;
2225 channel_similarity[CompositePixelChannel]+=ssim;
2226 }
2227 p+=(ptrdiff_t) GetPixelChannels(image);
2228 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2229 channel_area++;
2230 }
2231#if defined(MAGICKCORE_OPENMP_SUPPORT)
2232 #pragma omp critical (MagickCore_GetSSIMSimularity)
2233#endif
2234 {
2235 ssize_t
2236 j;
2237
2238 area+=channel_area;
2239 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
2240 {
2241 PixelChannel channel = GetPixelChannelChannel(image,j);
2242 PixelTrait traits = GetPixelChannelTraits(image,channel);
2243 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2244 channel);
2245 if (((traits & UpdatePixelTrait) == 0) ||
2246 ((reconstruct_traits & UpdatePixelTrait) == 0))
2247 continue;
2248 similarity[j]+=channel_similarity[j];
2249 }
2250 similarity[CompositePixelChannel]+=
2251 channel_similarity[CompositePixelChannel];
2252 }
2253 }
2254 image_view=DestroyCacheView(image_view);
2255 reconstruct_view=DestroyCacheView(reconstruct_view);
2256 area=MagickSafeReciprocal(area);
2257 for (l=0; l < (ssize_t) GetPixelChannels(image); l++)
2258 {
2259 PixelChannel channel = GetPixelChannelChannel(image,l);
2260 PixelTrait traits = GetPixelChannelTraits(image,channel);
2261 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2262 channel);
2263 if (((traits & UpdatePixelTrait) == 0) ||
2264 ((reconstruct_traits & UpdatePixelTrait) == 0))
2265 continue;
2266 similarity[l]*=area;
2267 }
2268 similarity[CompositePixelChannel]*=area;
2269 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
2270 kernel_info=DestroyKernelInfo(kernel_info);
2271 return(status);
2272}
2273
2274static MagickBooleanType GetDSSIMSimilarity(const Image *image,
2275 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2276{
2277 MagickBooleanType
2278 status = MagickTrue;
2279
2280 ssize_t
2281 i;
2282
2283 /*
2284 Compute the structual dissimilarity index similarity.
2285 */
2286 status=GetSSIMSimularity(image,reconstruct_image,similarity,exception);
2287 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2288 {
2289 PixelChannel channel = GetPixelChannelChannel(image,i);
2290 PixelTrait traits = GetPixelChannelTraits(image,channel);
2291 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2292 channel);
2293 if (((traits & UpdatePixelTrait) == 0) ||
2294 ((reconstruct_traits & UpdatePixelTrait) == 0))
2295 continue;
2296 similarity[i]=(1.0-similarity[i])/2.0;
2297 }
2298 similarity[CompositePixelChannel]=(1.0-similarity[CompositePixelChannel])/2.0;
2299 return(status);
2300}
2301
2302MagickExport MagickBooleanType GetImageDistortion(Image *image,
2303 const Image *reconstruct_image,const MetricType metric,double *distortion,
2304 ExceptionInfo *exception)
2305{
2306#define CompareMetricNotSupportedException "metric not supported"
2307
2308 double
2309 *channel_similarity;
2310
2311 MagickBooleanType
2312 status = MagickTrue;
2313
2314 size_t
2315 length;
2316
2317 assert(image != (Image *) NULL);
2318 assert(image->signature == MagickCoreSignature);
2319 assert(reconstruct_image != (const Image *) NULL);
2320 assert(reconstruct_image->signature == MagickCoreSignature);
2321 assert(distortion != (double *) NULL);
2322 if (IsEventLogging() != MagickFalse)
2323 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2324 /*
2325 Get image distortion.
2326 */
2327 *distortion=0.0;
2328 length=MaxPixelChannels+1UL;
2329 channel_similarity=(double *) AcquireQuantumMemory(length,
2330 sizeof(*channel_similarity));
2331 if (channel_similarity == (double *) NULL)
2332 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2333 (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
2334 switch (metric)
2335 {
2336 case AbsoluteErrorMetric:
2337 {
2338 status=GetAESimilarity(image,reconstruct_image,channel_similarity,
2339 exception);
2340 break;
2341 }
2342 case DotProductCorrelationErrorMetric:
2343 {
2344 status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
2345 exception);
2346 break;
2347 }
2348 case FuzzErrorMetric:
2349 {
2350 status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
2351 exception);
2352 break;
2353 }
2354 case MeanAbsoluteErrorMetric:
2355 {
2356 status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
2357 exception);
2358 break;
2359 }
2360 case MeanErrorPerPixelErrorMetric:
2361 {
2362 status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2363 exception);
2364 break;
2365 }
2366 case MeanSquaredErrorMetric:
2367 {
2368 status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
2369 exception);
2370 break;
2371 }
2372 case NormalizedCrossCorrelationErrorMetric:
2373 {
2374 status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
2375 exception);
2376 break;
2377 }
2378 case PeakAbsoluteErrorMetric:
2379 {
2380 status=GetPASimilarity(image,reconstruct_image,channel_similarity,
2381 exception);
2382 break;
2383 }
2384 case PeakSignalToNoiseRatioErrorMetric:
2385 {
2386 status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
2387 exception);
2388 break;
2389 }
2390 case PerceptualHashErrorMetric:
2391 {
2392 status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
2393 exception);
2394 break;
2395 }
2396 case PhaseCorrelationErrorMetric:
2397 {
2398 status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
2399 exception);
2400 break;
2401 }
2402 case PixelDifferenceCountErrorMetric:
2403 {
2404 status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
2405 exception);
2406 break;
2407 }
2408 case RootMeanSquaredErrorMetric:
2409 case UndefinedErrorMetric:
2410 default:
2411 {
2412 status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
2413 exception);
2414 break;
2415 }
2416 case StructuralDissimilarityErrorMetric:
2417 {
2418 status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
2419 exception);
2420 break;
2421 }
2422 case StructuralSimilarityErrorMetric:
2423 {
2424 status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
2425 exception);
2426 break;
2427 }
2428 }
2429 *distortion=channel_similarity[CompositePixelChannel];
2430 switch (metric)
2431 {
2432 case DotProductCorrelationErrorMetric:
2433 case NormalizedCrossCorrelationErrorMetric:
2434 case PhaseCorrelationErrorMetric:
2435 case StructuralSimilarityErrorMetric:
2436 {
2437 *distortion=(1.0-(*distortion))/2.0;
2438 break;
2439 }
2440 default: break;
2441 }
2442 channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
2443 if (fabs(*distortion) < MagickEpsilon)
2444 *distortion=0.0;
2445 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
2446 *distortion);
2447 return(status);
2448}
2449␌
2450/*
2451%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2452% %
2453% %
2454% %
2455% G e t I m a g e D i s t o r t i o n s %
2456% %
2457% %
2458% %
2459%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2460%
2461% GetImageDistortions() compares the pixel channels of an image to a
2462% reconstructed image and returns the specified metric for each channel.
2463%
2464% The format of the GetImageDistortions method is:
2465%
2466% double *GetImageDistortions(const Image *image,
2467% const Image *reconstruct_image,const MetricType metric,
2468% ExceptionInfo *exception)
2469%
2470% A description of each parameter follows:
2471%
2472% o image: the image.
2473%
2474% o reconstruct_image: the reconstruction image.
2475%
2476% o metric: the metric.
2477%
2478% o exception: return any errors or warnings in this structure.
2479%
2480*/
2481MagickExport double *GetImageDistortions(Image *image,
2482 const Image *reconstruct_image,const MetricType metric,
2483 ExceptionInfo *exception)
2484{
2485 double
2486 *distortion,
2487 *channel_similarity;
2488
2489 MagickBooleanType
2490 status = MagickTrue;
2491
2492 size_t
2493 length;
2494
2495 ssize_t
2496 i;
2497
2498 assert(image != (Image *) NULL);
2499 assert(image->signature == MagickCoreSignature);
2500 assert(reconstruct_image != (const Image *) NULL);
2501 assert(reconstruct_image->signature == MagickCoreSignature);
2502 if (IsEventLogging() != MagickFalse)
2503 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2504 /*
2505 Get image distortion.
2506 */
2507 length=MaxPixelChannels+1UL;
2508 channel_similarity=(double *) AcquireQuantumMemory(length,
2509 sizeof(*channel_similarity));
2510 if (channel_similarity == (double *) NULL)
2511 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2512 (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
2513 switch (metric)
2514 {
2515 case AbsoluteErrorMetric:
2516 {
2517 status=GetAESimilarity(image,reconstruct_image,channel_similarity,
2518 exception);
2519 break;
2520 }
2521 case DotProductCorrelationErrorMetric:
2522 {
2523 status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
2524 exception);
2525 break;
2526 }
2527 case FuzzErrorMetric:
2528 {
2529 status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
2530 exception);
2531 break;
2532 }
2533 case MeanAbsoluteErrorMetric:
2534 {
2535 status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
2536 exception);
2537 break;
2538 }
2539 case MeanErrorPerPixelErrorMetric:
2540 {
2541 status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2542 exception);
2543 break;
2544 }
2545 case MeanSquaredErrorMetric:
2546 {
2547 status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
2548 exception);
2549 break;
2550 }
2551 case NormalizedCrossCorrelationErrorMetric:
2552 {
2553 status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
2554 exception);
2555 break;
2556 }
2557 case PeakAbsoluteErrorMetric:
2558 {
2559 status=GetPASimilarity(image,reconstruct_image,channel_similarity,
2560 exception);
2561 break;
2562 }
2563 case PeakSignalToNoiseRatioErrorMetric:
2564 {
2565 status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
2566 exception);
2567 break;
2568 }
2569 case PerceptualHashErrorMetric:
2570 {
2571 status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
2572 exception);
2573 break;
2574 }
2575 case PhaseCorrelationErrorMetric:
2576 {
2577 status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
2578 exception);
2579 break;
2580 }
2581 case PixelDifferenceCountErrorMetric:
2582 {
2583 status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
2584 exception);
2585 break;
2586 }
2587 case RootMeanSquaredErrorMetric:
2588 case UndefinedErrorMetric:
2589 default:
2590 {
2591 status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
2592 exception);
2593 break;
2594 }
2595 case StructuralDissimilarityErrorMetric:
2596 {
2597 status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
2598 exception);
2599 break;
2600 }
2601 case StructuralSimilarityErrorMetric:
2602 {
2603 status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
2604 exception);
2605 break;
2606 }
2607 }
2608 if (status == MagickFalse)
2609 {
2610 channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
2611 return((double *) NULL);
2612 }
2613 distortion=channel_similarity;
2614 switch (metric)
2615 {
2616 case DotProductCorrelationErrorMetric:
2617 case NormalizedCrossCorrelationErrorMetric:
2618 case PhaseCorrelationErrorMetric:
2619 case StructuralSimilarityErrorMetric:
2620 {
2621 for (i=0; i <= MaxPixelChannels; i++)
2622 distortion[i]=(1.0-distortion[i])/2.0;
2623 break;
2624 }
2625 default: break;
2626 }
2627 for (i=0; i <= MaxPixelChannels; i++)
2628 if (fabs(distortion[i]) < MagickEpsilon)
2629 distortion[i]=0.0;
2630 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
2631 distortion[CompositePixelChannel]);
2632 return(distortion);
2633}
2634␌
2635/*
2636%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2637% %
2638% %
2639% %
2640% I s I m a g e s E q u a l %
2641% %
2642% %
2643% %
2644%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2645%
2646% IsImagesEqual() compare the pixels of two images and returns immediately
2647% if any pixel is not identical.
2648%
2649% The format of the IsImagesEqual method is:
2650%
2651% MagickBooleanType IsImagesEqual(const Image *image,
2652% const Image *reconstruct_image,ExceptionInfo *exception)
2653%
2654% A description of each parameter follows.
2655%
2656% o image: the image.
2657%
2658% o reconstruct_image: the reconstruction image.
2659%
2660% o exception: return any errors or warnings in this structure.
2661%
2662*/
2663MagickExport MagickBooleanType IsImagesEqual(const Image *image,
2664 const Image *reconstruct_image,ExceptionInfo *exception)
2665{
2666 CacheView
2667 *image_view,
2668 *reconstruct_view;
2669
2670 size_t
2671 columns,
2672 rows;
2673
2674 ssize_t
2675 y;
2676
2677 assert(image != (Image *) NULL);
2678 assert(image->signature == MagickCoreSignature);
2679 assert(reconstruct_image != (const Image *) NULL);
2680 assert(reconstruct_image->signature == MagickCoreSignature);
2681 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
2682 image_view=AcquireVirtualCacheView(image,exception);
2683 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2684 for (y=0; y < (ssize_t) rows; y++)
2685 {
2686 const Quantum
2687 *magick_restrict p,
2688 *magick_restrict q;
2689
2690 ssize_t
2691 x;
2692
2693 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2694 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2695 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
2696 break;
2697 for (x=0; x < (ssize_t) columns; x++)
2698 {
2699 ssize_t
2700 i;
2701
2702 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2703 {
2704 double
2705 distance;
2706
2707 PixelChannel channel = GetPixelChannelChannel(image,i);
2708 PixelTrait traits = GetPixelChannelTraits(image,channel);
2709 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2710 channel);
2711 if (((traits & UpdatePixelTrait) == 0) ||
2712 ((reconstruct_traits & UpdatePixelTrait) == 0))
2713 continue;
2714 distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
2715 channel,q));
2716 if (distance >= MagickEpsilon)
2717 break;
2718 }
2719 if (i < (ssize_t) GetPixelChannels(image))
2720 break;
2721 p+=(ptrdiff_t) GetPixelChannels(image);
2722 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2723 }
2724 if (x < (ssize_t) columns)
2725 break;
2726 }
2727 reconstruct_view=DestroyCacheView(reconstruct_view);
2728 image_view=DestroyCacheView(image_view);
2729 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
2730}
2731␌
2732/*
2733%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2734% %
2735% %
2736% %
2737% S e t I m a g e C o l o r M e t r i c %
2738% %
2739% %
2740% %
2741%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2742%
2743% SetImageColorMetric() measures the difference between colors at each pixel
2744% location of two images. A value other than 0 means the colors match
2745% exactly. Otherwise an error measure is computed by summing over all
2746% pixels in an image the distance squared in RGB space between each image
2747% pixel and its corresponding pixel in the reconstruction image. The error
2748% measure is assigned to these image members:
2749%
2750% o mean_error_per_pixel: The mean error for any single pixel in
2751% the image.
2752%
2753% o normalized_mean_error: The normalized mean quantization error for
2754% any single pixel in the image. This distance measure is normalized to
2755% a range between 0 and 1. It is independent of the range of red, green,
2756% and blue values in the image.
2757%
2758% o normalized_maximum_error: The normalized maximum quantization
2759% error for any single pixel in the image. This distance measure is
2760% normalized to a range between 0 and 1. It is independent of the range
2761% of red, green, and blue values in your image.
2762%
2763% A small normalized mean square error, accessed as
2764% image->normalized_mean_error, suggests the images are very similar in
2765% spatial layout and color.
2766%
2767% The format of the SetImageColorMetric method is:
2768%
2769% MagickBooleanType SetImageColorMetric(Image *image,
2770% const Image *reconstruct_image,ExceptionInfo *exception)
2771%
2772% A description of each parameter follows.
2773%
2774% o image: the image.
2775%
2776% o reconstruct_image: the reconstruction image.
2777%
2778% o exception: return any errors or warnings in this structure.
2779%
2780*/
2781MagickExport MagickBooleanType SetImageColorMetric(Image *image,
2782 const Image *reconstruct_image,ExceptionInfo *exception)
2783{
2784 double
2785 channel_similarity[MaxPixelChannels+1] = { 0.0 };
2786
2787 MagickBooleanType
2788 status;
2789
2790 status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2791 exception);
2792 if (status == MagickFalse)
2793 return(MagickFalse);
2794 status=fabs(image->error.mean_error_per_pixel) < MagickEpsilon ?
2795 MagickTrue : MagickFalse;
2796 return(status);
2797}
2798␌
2799/*
2800%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2801% %
2802% %
2803% %
2804% S i m i l a r i t y I m a g e %
2805% %
2806% %
2807% %
2808%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2809%
2810% SimilarityImage() compares the reconstruction of the image and returns the
2811% best match offset. In addition, it returns a similarity image such that an
2812% exact match location is completely white and if none of the pixels match,
2813% black, otherwise some gray level in-between.
2814%
2815% Contributed by Fred Weinhaus.
2816%
2817% The format of the SimilarityImageImage method is:
2818%
2819% Image *SimilarityImage(const Image *image,const Image *reconstruct,
2820% const MetricType metric,const double similarity_threshold,
2821% RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2822%
2823% A description of each parameter follows:
2824%
2825% o image: the image.
2826%
2827% o reconstruct: find an area of the image that closely resembles this image.
2828%
2829% o metric: the metric.
2830%
2831% o similarity_threshold: minimum similarity for (sub)image match.
2832%
2833% o offset: the best match offset of the reconstruction image within the
2834% image.
2835%
2836% o similarity: the computed similarity between the images.
2837%
2838% o exception: return any errors or warnings in this structure.
2839%
2840*/
2841
2842#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2843static Image *SIMCrossCorrelationImage(const Image *alpha_image,
2844 const Image *beta_image,ExceptionInfo *exception)
2845{
2846 Image
2847 *alpha_fft = (Image *) NULL,
2848 *beta_fft = (Image *) NULL,
2849 *complex_conjugate = (Image *) NULL,
2850 *complex_multiplication = (Image *) NULL,
2851 *cross_correlation = (Image *) NULL,
2852 *temp_image = (Image *) NULL;
2853
2854 /*
2855 Take the FFT of beta (reconstruction) image.
2856 */
2857 temp_image=CloneImage(beta_image,0,0,MagickTrue,exception);
2858 if (temp_image == (Image *) NULL)
2859 return((Image *) NULL);
2860 (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2861 beta_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2862 temp_image=DestroyImageList(temp_image);
2863 if (beta_fft == (Image *) NULL)
2864 return((Image *) NULL);
2865 /*
2866 Take the complex conjugate of beta_fft.
2867 */
2868 complex_conjugate=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
2869 beta_fft=DestroyImageList(beta_fft);
2870 if (complex_conjugate == (Image *) NULL)
2871 return((Image *) NULL);
2872 /*
2873 Take the FFT of the alpha (test) image.
2874 */
2875 temp_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2876 if (temp_image == (Image *) NULL)
2877 {
2878 complex_conjugate=DestroyImageList(complex_conjugate);
2879 return((Image *) NULL);
2880 }
2881 (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2882 alpha_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2883 temp_image=DestroyImageList(temp_image);
2884 if (alpha_fft == (Image *) NULL)
2885 {
2886 complex_conjugate=DestroyImageList(complex_conjugate);
2887 return((Image *) NULL);
2888 }
2889 /*
2890 Do complex multiplication.
2891 */
2892 DisableCompositeClampUnlessSpecified(complex_conjugate);
2893 DisableCompositeClampUnlessSpecified(complex_conjugate->next);
2894 AppendImageToList(&complex_conjugate,alpha_fft);
2895 complex_multiplication=ComplexImages(complex_conjugate,
2896 MultiplyComplexOperator,exception);
2897 complex_conjugate=DestroyImageList(complex_conjugate);
2898 if (complex_multiplication == (Image *) NULL)
2899 return((Image *) NULL);
2900 /*
2901 Do the IFT and return the cross-correlation result.
2902 */
2903 cross_correlation=InverseFourierTransformImage(complex_multiplication,
2904 complex_multiplication->next,MagickFalse,exception);
2905 complex_multiplication=DestroyImageList(complex_multiplication);
2906 return(cross_correlation);
2907}
2908
2909static Image *SIMDerivativeImage(const Image *image,const char *kernel,
2910 ExceptionInfo *exception)
2911{
2912 Image
2913 *derivative_image;
2914
2915 KernelInfo
2916 *kernel_info;
2917
2918 kernel_info=AcquireKernelInfo(kernel,exception);
2919 if (kernel_info == (KernelInfo *) NULL)
2920 return((Image *) NULL);
2921 derivative_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,
2922 exception);
2923 kernel_info=DestroyKernelInfo(kernel_info);
2924 return(derivative_image);
2925}
2926
2927static Image *SIMDivideImage(const Image *numerator_image,
2928 const Image *denominator_image,ExceptionInfo *exception)
2929{
2930 CacheView
2931 *denominator_view,
2932 *numerator_view;
2933
2934 Image
2935 *divide_image;
2936
2937 MagickBooleanType
2938 status = MagickTrue;
2939
2940 ssize_t
2941 y;
2942
2943 /*
2944 Divide one image into another.
2945 */
2946 divide_image=CloneImage(numerator_image,0,0,MagickTrue,exception);
2947 if (divide_image == (Image *) NULL)
2948 return(divide_image);
2949 numerator_view=AcquireAuthenticCacheView(divide_image,exception);
2950 denominator_view=AcquireVirtualCacheView(denominator_image,exception);
2951#if defined(MAGICKCORE_OPENMP_SUPPORT)
2952 #pragma omp parallel for schedule(static) shared(status) \
2953 magick_number_threads(denominator_image,divide_image,divide_image->rows,1)
2954#endif
2955 for (y=0; y < (ssize_t) divide_image->rows; y++)
2956 {
2957 const Quantum
2958 *magick_restrict p;
2959
2960 Quantum
2961 *magick_restrict q;
2962
2963 ssize_t
2964 x;
2965
2966 if (status == MagickFalse)
2967 continue;
2968 p=GetCacheViewVirtualPixels(denominator_view,0,y,
2969 denominator_image->columns,1,exception);
2970 q=GetCacheViewAuthenticPixels(numerator_view,0,y,divide_image->columns,1,
2971 exception);
2972 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2973 {
2974 status=MagickFalse;
2975 continue;
2976 }
2977 for (x=0; x < (ssize_t) divide_image->columns; x++)
2978 {
2979 ssize_t
2980 i;
2981
2982 for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
2983 {
2984 PixelChannel channel = GetPixelChannelChannel(divide_image,i);
2985 PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
2986 PixelTrait denominator_traits = GetPixelChannelTraits(denominator_image,
2987 channel);
2988 if (((traits & UpdatePixelTrait) == 0) ||
2989 ((denominator_traits & UpdatePixelTrait) == 0))
2990 continue;
2991 q[i]=(Quantum) ((double) q[i]*MagickSafeReciprocal(QuantumScale*
2992 (double) GetPixelChannel(denominator_image,channel,p)));
2993 }
2994 p+=(ptrdiff_t) GetPixelChannels(denominator_image);
2995 q+=(ptrdiff_t) GetPixelChannels(divide_image);
2996 }
2997 if (SyncCacheViewAuthenticPixels(numerator_view,exception) == MagickFalse)
2998 status=MagickFalse;
2999 }
3000 denominator_view=DestroyCacheView(denominator_view);
3001 numerator_view=DestroyCacheView(numerator_view);
3002 if (status == MagickFalse)
3003 divide_image=DestroyImage(divide_image);
3004 return(divide_image);
3005}
3006
3007static Image *SIMDivideByMagnitude(Image *image,Image *magnitude_image,
3008 const Image *source_image,ExceptionInfo *exception)
3009{
3010 Image
3011 *divide_image,
3012 *result_image;
3013
3014 RectangleInfo
3015 geometry;
3016
3017 divide_image=SIMDivideImage(image,magnitude_image,exception);
3018 if (divide_image == (Image *) NULL)
3019 return((Image *) NULL);
3020 GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
3021 &divide_image->background_color);
3022 SetGeometry(source_image,&geometry);
3023 geometry.width=MagickMax(source_image->columns,divide_image->columns);
3024 geometry.height=MagickMax(source_image->rows,divide_image->rows);
3025 result_image=ExtentImage(divide_image,&geometry,exception);
3026 divide_image=DestroyImage(divide_image);
3027 return(result_image);
3028}
3029
3030static MagickBooleanType SIMFilterImageNaNs(Image *image,
3031 ExceptionInfo *exception)
3032{
3033 CacheView
3034 *image_view;
3035
3036 MagickBooleanType
3037 status = MagickTrue;
3038
3039 ssize_t
3040 y;
3041
3042 /*
3043 Square each pixel in the image.
3044 */
3045 image_view=AcquireAuthenticCacheView(image,exception);
3046#if defined(MAGICKCORE_OPENMP_SUPPORT)
3047 #pragma omp parallel for schedule(static) shared(status) \
3048 magick_number_threads(image,image,image->rows,1)
3049#endif
3050 for (y=0; y < (ssize_t) image->rows; y++)
3051 {
3052 Quantum
3053 *magick_restrict q;
3054
3055 ssize_t
3056 x;
3057
3058 if (status == MagickFalse)
3059 continue;
3060 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3061 if (q == (Quantum *) NULL)
3062 {
3063 status=MagickFalse;
3064 continue;
3065 }
3066 for (x=0; x < (ssize_t) image->columns; x++)
3067 {
3068 ssize_t
3069 i;
3070
3071 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3072 {
3073 PixelChannel channel = GetPixelChannelChannel(image,i);
3074 PixelTrait traits = GetPixelChannelTraits(image,channel);
3075 if ((traits & UpdatePixelTrait) == 0)
3076 continue;
3077 if (IsNaN((double) q[i]) != 0)
3078 q[i]=(Quantum) 0;
3079 }
3080 q+=(ptrdiff_t) GetPixelChannels(image);
3081 }
3082 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3083 status=MagickFalse;
3084 }
3085 image_view=DestroyCacheView(image_view);
3086 return(status);
3087}
3088
3089static Image *SIMSquareImage(const Image *image,ExceptionInfo *exception)
3090{
3091 CacheView
3092 *image_view;
3093
3094 Image
3095 *square_image;
3096
3097 MagickBooleanType
3098 status = MagickTrue;
3099
3100 ssize_t
3101 y;
3102
3103 /*
3104 Square each pixel in the image.
3105 */
3106 square_image=CloneImage(image,0,0,MagickTrue,exception);
3107 if (square_image == (Image *) NULL)
3108 return(square_image);
3109 image_view=AcquireAuthenticCacheView(square_image,exception);
3110#if defined(MAGICKCORE_OPENMP_SUPPORT)
3111 #pragma omp parallel for schedule(static) shared(status) \
3112 magick_number_threads(square_image,square_image,square_image->rows,1)
3113#endif
3114 for (y=0; y < (ssize_t) square_image->rows; y++)
3115 {
3116 Quantum
3117 *magick_restrict q;
3118
3119 ssize_t
3120 x;
3121
3122 if (status == MagickFalse)
3123 continue;
3124 q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
3125 exception);
3126 if (q == (Quantum *) NULL)
3127 {
3128 status=MagickFalse;
3129 continue;
3130 }
3131 for (x=0; x < (ssize_t) square_image->columns; x++)
3132 {
3133 ssize_t
3134 i;
3135
3136 for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
3137 {
3138 PixelChannel channel = GetPixelChannelChannel(square_image,i);
3139 PixelTrait traits = GetPixelChannelTraits(square_image,channel);
3140 if ((traits & UpdatePixelTrait) == 0)
3141 continue;
3142 q[i]=(Quantum) (QuantumScale*q[i]*q[i]);
3143 }
3144 q+=(ptrdiff_t) GetPixelChannels(square_image);
3145 }
3146 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3147 status=MagickFalse;
3148 }
3149 image_view=DestroyCacheView(image_view);
3150 if (status == MagickFalse)
3151 square_image=DestroyImage(square_image);
3152 return(square_image);
3153}
3154
3155static Image *SIMMagnitudeImage(Image *alpha_image,Image *beta_image,
3156 ExceptionInfo *exception)
3157{
3158 Image
3159 *magnitude_image,
3160 *xsq_image,
3161 *ysq_image;
3162
3163 MagickBooleanType
3164 status = MagickTrue;
3165
3166 (void) SetImageArtifact(alpha_image,"compose:clamp","False");
3167 xsq_image=SIMSquareImage(alpha_image,exception);
3168 if (xsq_image == (Image *) NULL)
3169 return((Image *) NULL);
3170 (void) SetImageArtifact(beta_image,"compose:clamp","False");
3171 ysq_image=SIMSquareImage(beta_image,exception);
3172 if (ysq_image == (Image *) NULL)
3173 {
3174 xsq_image=DestroyImage(xsq_image);
3175 return((Image *) NULL);
3176 }
3177 status=CompositeImage(xsq_image,ysq_image,PlusCompositeOp,MagickTrue,0,0,
3178 exception);
3179 magnitude_image=xsq_image;
3180 ysq_image=DestroyImage(ysq_image);
3181 if (status == MagickFalse)
3182 {
3183 magnitude_image=DestroyImage(magnitude_image);
3184 return((Image *) NULL);
3185 }
3186 status=EvaluateImage(magnitude_image,PowEvaluateOperator,0.5,exception);
3187 if (status == MagickFalse)
3188 {
3189 magnitude_image=DestroyImage(magnitude_image);
3190 return (Image *) NULL;
3191 }
3192 return(magnitude_image);
3193}
3194
3195static MagickBooleanType SIMMaximaImage(const Image *image,double *maxima,
3196 RectangleInfo *offset,ExceptionInfo *exception)
3197{
3198 typedef struct
3199 {
3200 double
3201 maxima;
3202
3203 ssize_t
3204 x,
3205 y;
3206 } MaximaInfo;
3207
3208 CacheView
3209 *image_view;
3210
3211 const Quantum
3212 *magick_restrict q;
3213
3214 MagickBooleanType
3215 status = MagickTrue;
3216
3217 MaximaInfo
3218 maxima_info = { -MagickMaximumValue, 0, 0 };
3219
3220 ssize_t
3221 y;
3222
3223 /*
3224 Identify the maxima value in the image and its location.
3225 */
3226 image_view=AcquireVirtualCacheView(image,exception);
3227 q=GetCacheViewVirtualPixels(image_view,maxima_info.x,maxima_info.y,1,1,
3228 exception);
3229 if (q != (const Quantum *) NULL)
3230 maxima_info.maxima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
3231#if defined(MAGICKCORE_OPENMP_SUPPORT)
3232 #pragma omp parallel for schedule(static) shared(maxima_info,status) \
3233 magick_number_threads(image,image,image->rows,1)
3234#endif
3235 for (y=0; y < (ssize_t) image->rows; y++)
3236 {
3237 const Quantum
3238 *magick_restrict p;
3239
3240 MaximaInfo
3241 channel_maxima = { -MagickMaximumValue, 0, 0 };
3242
3243 ssize_t
3244 x;
3245
3246 if (status == MagickFalse)
3247 continue;
3248 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3249 if (p == (const Quantum *) NULL)
3250 {
3251 status=MagickFalse;
3252 continue;
3253 }
3254 channel_maxima=maxima_info;
3255 for (x=0; x < (ssize_t) image->columns; x++)
3256 {
3257 ssize_t
3258 i;
3259
3260 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3261 {
3262 double
3263 pixel;
3264
3265 PixelChannel channel = GetPixelChannelChannel(image,i);
3266 PixelTrait traits = GetPixelChannelTraits(image,channel);
3267 if ((traits & UpdatePixelTrait) == 0)
3268 continue;
3269 pixel=(double) p[i];
3270 if (IsNaN(pixel) != 0)
3271 pixel=0.0;
3272 if (pixel > channel_maxima.maxima)
3273 {
3274 channel_maxima.maxima=(double) p[i];
3275 channel_maxima.x=x;
3276 channel_maxima.y=y;
3277 }
3278 }
3279 p+=(ptrdiff_t) GetPixelChannels(image);
3280 }
3281#if defined(MAGICKCORE_OPENMP_SUPPORT)
3282 #pragma omp critical (MagickCore_SIMMaximaImage)
3283#endif
3284 if (channel_maxima.maxima > maxima_info.maxima)
3285 maxima_info=channel_maxima;
3286 }
3287 image_view=DestroyCacheView(image_view);
3288 *maxima=maxima_info.maxima;
3289 offset->x=maxima_info.x;
3290 offset->y=maxima_info.y;
3291 return(status);
3292}
3293
3294static MagickBooleanType SIMMinimaImage(const Image *image,double *minima,
3295 RectangleInfo *offset,ExceptionInfo *exception)
3296{
3297 typedef struct
3298 {
3299 double
3300 minima;
3301
3302 ssize_t
3303 x,
3304 y;
3305 } MinimaInfo;
3306
3307 CacheView
3308 *image_view;
3309
3310 const Quantum
3311 *magick_restrict q;
3312
3313 MagickBooleanType
3314 status = MagickTrue;
3315
3316 MinimaInfo
3317 minima_info = { MagickMaximumValue, 0, 0 };
3318
3319 ssize_t
3320 y;
3321
3322 /*
3323 Identify the minima value in the image and its location.
3324 */
3325 image_view=AcquireVirtualCacheView(image,exception);
3326 q=GetCacheViewVirtualPixels(image_view,minima_info.x,minima_info.y,1,1,
3327 exception);
3328 if (q != (const Quantum *) NULL)
3329 minima_info.minima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
3330#if defined(MAGICKCORE_OPENMP_SUPPORT)
3331 #pragma omp parallel for schedule(static) shared(minima_info,status) \
3332 magick_number_threads(image,image,image->rows,1)
3333#endif
3334 for (y=0; y < (ssize_t) image->rows; y++)
3335 {
3336 const Quantum
3337 *magick_restrict p;
3338
3339 MinimaInfo
3340 channel_minima = { MagickMaximumValue, 0, 0 };
3341
3342 ssize_t
3343 x;
3344
3345 if (status == MagickFalse)
3346 continue;
3347 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3348 if (p == (const Quantum *) NULL)
3349 {
3350 status=MagickFalse;
3351 continue;
3352 }
3353 channel_minima=minima_info;
3354 for (x=0; x < (ssize_t) image->columns; x++)
3355 {
3356 ssize_t
3357 i;
3358
3359 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3360 {
3361 double
3362 pixel;
3363
3364 PixelChannel channel = GetPixelChannelChannel(image,i);
3365 PixelTrait traits = GetPixelChannelTraits(image,channel);
3366 if ((traits & UpdatePixelTrait) == 0)
3367 continue;
3368 pixel=(double) p[i];
3369 if (IsNaN(pixel) != 0)
3370 pixel=0.0;
3371 if (pixel < channel_minima.minima)
3372 {
3373 channel_minima.minima=pixel;
3374 channel_minima.x=x;
3375 channel_minima.y=y;
3376 }
3377 }
3378 p+=(ptrdiff_t) GetPixelChannels(image);
3379 }
3380#if defined(MAGICKCORE_OPENMP_SUPPORT)
3381 #pragma omp critical (MagickCore_SIMMinimaImage)
3382#endif
3383 if (channel_minima.minima < minima_info.minima)
3384 minima_info=channel_minima;
3385 }
3386 image_view=DestroyCacheView(image_view);
3387 *minima=minima_info.minima;
3388 offset->x=minima_info.x;
3389 offset->y=minima_info.y;
3390 return(status);
3391}
3392
3393static MagickBooleanType SIMMultiplyImage(Image *image,const double factor,
3394 const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
3395{
3396 CacheView
3397 *image_view;
3398
3399 MagickBooleanType
3400 status = MagickTrue;
3401
3402 ssize_t
3403 y;
3404
3405 /*
3406 Multiply each pixel by a factor.
3407 */
3408 image_view=AcquireAuthenticCacheView(image,exception);
3409#if defined(MAGICKCORE_OPENMP_SUPPORT)
3410 #pragma omp parallel for schedule(static) shared(status) \
3411 magick_number_threads(image,image,image->rows,1)
3412#endif
3413 for (y=0; y < (ssize_t) image->rows; y++)
3414 {
3415 Quantum
3416 *magick_restrict q;
3417
3418 ssize_t
3419 x;
3420
3421 if (status == MagickFalse)
3422 continue;
3423 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3424 if (q == (Quantum *) NULL)
3425 {
3426 status=MagickFalse;
3427 continue;
3428 }
3429 for (x=0; x < (ssize_t) image->columns; x++)
3430 {
3431 ssize_t
3432 i;
3433
3434 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3435 {
3436 PixelChannel channel = GetPixelChannelChannel(image,i);
3437 PixelTrait traits = GetPixelChannelTraits(image,channel);
3438 if ((traits & UpdatePixelTrait) == 0)
3439 continue;
3440 if (channel_statistics != (const ChannelStatistics *) NULL)
3441 q[i]=(Quantum) (factor*q[i]*QuantumScale*
3442 channel_statistics[channel].standard_deviation);
3443 else
3444 q[i]=(Quantum) (factor*q[i]);
3445 }
3446 q+=(ptrdiff_t) GetPixelChannels(image);
3447 }
3448 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3449 status=MagickFalse;
3450 }
3451 image_view=DestroyCacheView(image_view);
3452 return(status);
3453}
3454
3455static Image *SIMPhaseCorrelationImage(const Image *alpha_image,
3456 const Image *beta_image,const Image *magnitude_image,ExceptionInfo *exception)
3457{
3458 Image
3459 *alpha_fft = (Image *) NULL,
3460 *beta_fft = (Image *) NULL,
3461 *complex_multiplication = (Image *) NULL,
3462 *cross_correlation = (Image *) NULL;
3463
3464 /*
3465 Take the FFT of the beta (reconstruction) image.
3466 */
3467 beta_fft=CloneImage(beta_image,0,0,MagickTrue,exception);
3468 if (beta_fft == NULL)
3469 return((Image *) NULL);
3470 (void) SetImageArtifact(beta_fft,"fourier:normalize","inverse");
3471 beta_fft=ForwardFourierTransformImage(beta_fft,MagickFalse,exception);
3472 if (beta_fft == NULL)
3473 return((Image *) NULL);
3474 /*
3475 Take the FFT of the alpha (test) image.
3476 */
3477 alpha_fft=CloneImage(alpha_image,0,0,MagickTrue,exception);
3478 if (alpha_fft == (Image *) NULL)
3479 {
3480 beta_fft=DestroyImageList(beta_fft);
3481 return((Image *) NULL);
3482 }
3483 (void) SetImageArtifact(alpha_fft,"fourier:normalize","inverse");
3484 alpha_fft=ForwardFourierTransformImage(alpha_fft,MagickFalse,exception);
3485 if (alpha_fft == (Image *) NULL)
3486 {
3487 beta_fft=DestroyImageList(beta_fft);
3488 return((Image *) NULL);
3489 }
3490 /*
3491 Take the complex conjugate of the beta FFT.
3492 */
3493 beta_fft=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
3494 if (beta_fft == (Image *) NULL)
3495 {
3496 alpha_fft=DestroyImageList(alpha_fft);
3497 return((Image *) NULL);
3498 }
3499 /*
3500 Do complex multiplication.
3501 */
3502 AppendImageToList(&beta_fft,alpha_fft);
3503 DisableCompositeClampUnlessSpecified(beta_fft);
3504 DisableCompositeClampUnlessSpecified(beta_fft->next);
3505 complex_multiplication=ComplexImages(beta_fft,MultiplyComplexOperator,
3506 exception);
3507 beta_fft=DestroyImageList(beta_fft);
3508 if (complex_multiplication == (Image *) NULL)
3509 return((Image *) NULL);
3510 /*
3511 Divide the results.
3512 */
3513 CompositeLayers(complex_multiplication,DivideSrcCompositeOp,(Image *)
3514 magnitude_image,0,0,exception);
3515 /*
3516 Do the IFT and return the cross-correlation result.
3517 */
3518 (void) SetImageArtifact(complex_multiplication,"fourier:normalize","inverse");
3519 cross_correlation=InverseFourierTransformImage(complex_multiplication,
3520 complex_multiplication->next,MagickFalse,exception);
3521 complex_multiplication=DestroyImageList(complex_multiplication);
3522 return(cross_correlation);
3523}
3524
3525static MagickBooleanType SIMSetImageMean(Image *image,
3526 const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
3527{
3528 CacheView
3529 *image_view;
3530
3531 MagickBooleanType
3532 status = MagickTrue;
3533
3534 ssize_t
3535 y;
3536
3537 /*
3538 Set image mean.
3539 */
3540 image_view=AcquireAuthenticCacheView(image,exception);
3541#if defined(MAGICKCORE_OPENMP_SUPPORT)
3542 #pragma omp parallel for schedule(static) shared(status) \
3543 magick_number_threads(image,image,image->rows,1)
3544#endif
3545 for (y=0; y < (ssize_t) image->rows; y++)
3546 {
3547 Quantum
3548 *magick_restrict q;
3549
3550 ssize_t
3551 x;
3552
3553 if (status == MagickFalse)
3554 continue;
3555 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3556 if (q == (Quantum *) NULL)
3557 {
3558 status=MagickFalse;
3559 continue;
3560 }
3561 for (x=0; x < (ssize_t) image->columns; x++)
3562 {
3563 ssize_t
3564 i;
3565
3566 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3567 {
3568 PixelChannel channel = GetPixelChannelChannel(image,i);
3569 PixelTrait traits = GetPixelChannelTraits(image,channel);
3570 if ((traits & UpdatePixelTrait) == 0)
3571 continue;
3572 q[i]=(Quantum) channel_statistics[channel].mean;
3573 }
3574 q+=(ptrdiff_t) GetPixelChannels(image);
3575 }
3576 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3577 status=MagickFalse;
3578 }
3579 image_view=DestroyCacheView(image_view);
3580 return(status);
3581}
3582
3583static Image *SIMSubtractImageMean(const Image *alpha_image,
3584 const Image *beta_image,const ChannelStatistics *channel_statistics,
3585 ExceptionInfo *exception)
3586{
3587 CacheView
3588 *beta_view,
3589 *image_view;
3590
3591 Image
3592 *subtract_image;
3593
3594 MagickBooleanType
3595 status = MagickTrue;
3596
3597 ssize_t
3598 y;
3599
3600 /*
3601 Subtract the image mean and pad.
3602 */
3603 subtract_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
3604 MagickTrue,exception);
3605 if (subtract_image == (Image *) NULL)
3606 return(subtract_image);
3607 image_view=AcquireAuthenticCacheView(subtract_image,exception);
3608 beta_view=AcquireVirtualCacheView(beta_image,exception);
3609#if defined(MAGICKCORE_OPENMP_SUPPORT)
3610 #pragma omp parallel for schedule(static) shared(status) \
3611 magick_number_threads(beta_image,subtract_image,subtract_image->rows,1)
3612#endif
3613 for (y=0; y < (ssize_t) subtract_image->rows; y++)
3614 {
3615 const Quantum
3616 *magick_restrict p;
3617
3618 Quantum
3619 *magick_restrict q;
3620
3621 ssize_t
3622 x;
3623
3624 if (status == MagickFalse)
3625 continue;
3626 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,exception);
3627 q=GetCacheViewAuthenticPixels(image_view,0,y,subtract_image->columns,1,
3628 exception);
3629 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3630 {
3631 status=MagickFalse;
3632 continue;
3633 }
3634 for (x=0; x < (ssize_t) subtract_image->columns; x++)
3635 {
3636 ssize_t
3637 i;
3638
3639 for (i=0; i < (ssize_t) GetPixelChannels(subtract_image); i++)
3640 {
3641 PixelChannel channel = GetPixelChannelChannel(subtract_image,i);
3642 PixelTrait traits = GetPixelChannelTraits(subtract_image,channel);
3643 PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
3644 if (((traits & UpdatePixelTrait) == 0) ||
3645 ((beta_traits & UpdatePixelTrait) == 0))
3646 continue;
3647 if ((x >= (ssize_t) beta_image->columns) ||
3648 (y >= (ssize_t) beta_image->rows))
3649 q[i]=(Quantum) 0;
3650 else
3651 q[i]=(Quantum) ((double) GetPixelChannel(beta_image,channel,p)-
3652 channel_statistics[channel].mean);
3653 }
3654 p+=(ptrdiff_t) GetPixelChannels(beta_image);
3655 q+=(ptrdiff_t) GetPixelChannels(subtract_image);
3656 }
3657 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3658 status=MagickFalse;
3659 }
3660 beta_view=DestroyCacheView(beta_view);
3661 image_view=DestroyCacheView(image_view);
3662 if (status == MagickFalse)
3663 subtract_image=DestroyImage(subtract_image);
3664 return(subtract_image);
3665}
3666
3667static Image *SIMUnityImage(const Image *alpha_image,const Image *beta_image,
3668 ExceptionInfo *exception)
3669{
3670 CacheView
3671 *image_view;
3672
3673 Image
3674 *unity_image;
3675
3676 MagickBooleanType
3677 status = MagickTrue;
3678
3679 ssize_t
3680 y;
3681
3682 /*
3683 Create a padded unity image.
3684 */
3685 unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
3686 MagickTrue,exception);
3687 if (unity_image == (Image *) NULL)
3688 return(unity_image);
3689 if (SetImageStorageClass(unity_image,DirectClass,exception) == MagickFalse)
3690 return(DestroyImage(unity_image));
3691 image_view=AcquireAuthenticCacheView(unity_image,exception);
3692#if defined(MAGICKCORE_OPENMP_SUPPORT)
3693 #pragma omp parallel for schedule(static) shared(status) \
3694 magick_number_threads(unity_image,unity_image,unity_image->rows,1)
3695#endif
3696 for (y=0; y < (ssize_t) unity_image->rows; y++)
3697 {
3698 Quantum
3699 *magick_restrict q;
3700
3701 ssize_t
3702 x;
3703
3704 if (status == MagickFalse)
3705 continue;
3706 q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
3707 exception);
3708 if (q == (Quantum *) NULL)
3709 {
3710 status=MagickFalse;
3711 continue;
3712 }
3713 for (x=0; x < (ssize_t) unity_image->columns; x++)
3714 {
3715 ssize_t
3716 i;
3717
3718 for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
3719 {
3720 PixelChannel channel = GetPixelChannelChannel(unity_image,i);
3721 PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
3722 if ((traits & UpdatePixelTrait) == 0)
3723 continue;
3724 if ((x >= (ssize_t) beta_image->columns) ||
3725 (y >= (ssize_t) beta_image->rows))
3726 q[i]=(Quantum) 0;
3727 else
3728 q[i]=QuantumRange;
3729 }
3730 q+=(ptrdiff_t) GetPixelChannels(unity_image);
3731 }
3732 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3733 status=MagickFalse;
3734 }
3735 image_view=DestroyCacheView(image_view);
3736 if (status == MagickFalse)
3737 unity_image=DestroyImage(unity_image);
3738 return(unity_image);
3739}
3740
3741static Image *SIMVarianceImage(Image *alpha_image,const Image *beta_image,
3742 ExceptionInfo *exception)
3743{
3744 CacheView
3745 *beta_view,
3746 *image_view;
3747
3748 Image
3749 *variance_image;
3750
3751 MagickBooleanType
3752 status = MagickTrue;
3753
3754 ssize_t
3755 y;
3756
3757 /*
3758 Compute the variance of the two images.
3759 */
3760 variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
3761 if (variance_image == (Image *) NULL)
3762 return(variance_image);
3763 image_view=AcquireAuthenticCacheView(variance_image,exception);
3764 beta_view=AcquireVirtualCacheView(beta_image,exception);
3765#if defined(MAGICKCORE_OPENMP_SUPPORT)
3766 #pragma omp parallel for schedule(static) shared(status) \
3767 magick_number_threads(beta_image,variance_image,variance_image->rows,1)
3768#endif
3769 for (y=0; y < (ssize_t) variance_image->rows; y++)
3770 {
3771 const Quantum
3772 *magick_restrict p;
3773
3774 Quantum
3775 *magick_restrict q;
3776
3777 ssize_t
3778 x;
3779
3780 if (status == MagickFalse)
3781 continue;
3782 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
3783 exception);
3784 q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
3785 exception);
3786 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3787 {
3788 status=MagickFalse;
3789 continue;
3790 }
3791 for (x=0; x < (ssize_t) variance_image->columns; x++)
3792 {
3793 ssize_t
3794 i;
3795
3796 for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
3797 {
3798 double
3799 error;
3800
3801 PixelChannel channel = GetPixelChannelChannel(variance_image,i);
3802 PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
3803 PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
3804 if (((traits & UpdatePixelTrait) == 0) ||
3805 ((beta_traits & UpdatePixelTrait) == 0))
3806 continue;
3807 error=(double) q[i]-(double) GetPixelChannel(beta_image,channel,p);
3808 q[i]=(Quantum) ((double) ClampToQuantum((double) QuantumRange*
3809 (sqrt(fabs(QuantumScale*error))/sqrt((double) QuantumRange))));
3810 }
3811 p+=(ptrdiff_t) GetPixelChannels(beta_image);
3812 q+=(ptrdiff_t) GetPixelChannels(variance_image);
3813 }
3814 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3815 status=MagickFalse;
3816 }
3817 beta_view=DestroyCacheView(beta_view);
3818 image_view=DestroyCacheView(image_view);
3819 if (status == MagickFalse)
3820 variance_image=DestroyImage(variance_image);
3821 return(variance_image);
3822}
3823
3824static Image *DPCSimilarityImage(const Image *image,const Image *reconstruct,
3825 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3826{
3827#define ThrowDPCSimilarityException() \
3828{ \
3829 if (dot_product_image != (Image *) NULL) \
3830 dot_product_image=DestroyImage(dot_product_image); \
3831 if (magnitude_image != (Image *) NULL) \
3832 magnitude_image=DestroyImage(magnitude_image); \
3833 if (reconstruct_image != (Image *) NULL) \
3834 reconstruct_image=DestroyImage(reconstruct_image); \
3835 if (rx_image != (Image *) NULL) \
3836 rx_image=DestroyImage(rx_image); \
3837 if (ry_image != (Image *) NULL) \
3838 ry_image=DestroyImage(ry_image); \
3839 if (target_image != (Image *) NULL) \
3840 target_image=DestroyImage(target_image); \
3841 if (threshold_image != (Image *) NULL) \
3842 threshold_image=DestroyImage(threshold_image); \
3843 if (trx_image != (Image *) NULL) \
3844 trx_image=DestroyImage(trx_image); \
3845 if (try_image != (Image *) NULL) \
3846 try_image=DestroyImage(try_image); \
3847 if (tx_image != (Image *) NULL) \
3848 tx_image=DestroyImage(tx_image); \
3849 if (ty_image != (Image *) NULL) \
3850 ty_image=DestroyImage(ty_image); \
3851 return((Image *) NULL); \
3852}
3853
3854 double
3855 edge_factor = 0.0,
3856 maxima = 0.0,
3857 mean = 0.0,
3858 standard_deviation = 0.0;
3859
3860 Image
3861 *dot_product_image = (Image *) NULL,
3862 *magnitude_image = (Image *) NULL,
3863 *reconstruct_image = (Image *) NULL,
3864 *rx_image = (Image *) NULL,
3865 *ry_image = (Image *) NULL,
3866 *trx_image = (Image *) NULL,
3867 *target_image = (Image *) NULL,
3868 *threshold_image = (Image *) NULL,
3869 *try_image = (Image *) NULL,
3870 *tx_image = (Image *) NULL,
3871 *ty_image = (Image *) NULL;
3872
3873 MagickBooleanType
3874 status = MagickTrue;
3875
3876 RectangleInfo
3877 geometry;
3878
3879 /*
3880 Dot product correlation-based image similarity using FFT local statistics.
3881 */
3882 target_image=CloneImage(image,0,0,MagickTrue,exception);
3883 if (target_image == (Image *) NULL)
3884 return((Image *) NULL);
3885 /*
3886 Compute the cross correlation of the test and reconstruct magnitudes.
3887 */
3888 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
3889 if (reconstruct_image == (Image *) NULL)
3890 ThrowDPCSimilarityException();
3891 /*
3892 Compute X and Y derivatives of reference image.
3893 */
3894 (void) SetImageVirtualPixelMethod(reconstruct_image,EdgeVirtualPixelMethod,
3895 exception);
3896 rx_image=SIMDerivativeImage(reconstruct_image,"Sobel",exception);
3897 if (rx_image == (Image *) NULL)
3898 ThrowDPCSimilarityException();
3899 ry_image=SIMDerivativeImage(reconstruct_image,"Sobel:90",exception);
3900 reconstruct_image=DestroyImage(reconstruct_image);
3901 if (ry_image == (Image *) NULL)
3902 ThrowDPCSimilarityException();
3903 /*
3904 Compute magnitude of derivatives.
3905 */
3906 magnitude_image=SIMMagnitudeImage(rx_image,ry_image,exception);
3907 if (magnitude_image == (Image *) NULL)
3908 ThrowDPCSimilarityException();
3909 /*
3910 Compute an edge normalization correction.
3911 */
3912 threshold_image=CloneImage(magnitude_image,0,0,MagickTrue,exception);
3913 if (threshold_image == (Image *) NULL)
3914 ThrowDPCSimilarityException();
3915 status=BilevelImage(threshold_image,0.0,exception);
3916 if (status == MagickFalse)
3917 ThrowDPCSimilarityException();
3918 status=GetImageMean(threshold_image,&mean,&standard_deviation,exception);
3919 threshold_image=DestroyImage(threshold_image);
3920 if (status == MagickFalse)
3921 ThrowDPCSimilarityException();
3922 edge_factor=MagickSafeReciprocal(QuantumScale*mean*reconstruct->columns*
3923 reconstruct->rows)+QuantumScale;
3924 /*
3925 Divide X and Y derivitives of reference image by magnitude.
3926 */
3927 trx_image=SIMDivideByMagnitude(rx_image,magnitude_image,image,exception);
3928 rx_image=DestroyImage(rx_image);
3929 if (trx_image == (Image *) NULL)
3930 ThrowDPCSimilarityException();
3931 rx_image=trx_image;
3932 try_image=SIMDivideByMagnitude(ry_image,magnitude_image,image,exception);
3933 magnitude_image=DestroyImage(magnitude_image);
3934 ry_image=DestroyImage(ry_image);
3935 if (try_image == (Image *) NULL)
3936 ThrowDPCSimilarityException();
3937 ry_image=try_image;
3938 /*
3939 Compute X and Y derivatives of image.
3940 */
3941 (void) SetImageVirtualPixelMethod(target_image,EdgeVirtualPixelMethod,
3942 exception);
3943 tx_image=SIMDerivativeImage(target_image,"Sobel",exception);
3944 if (tx_image == (Image *) NULL)
3945 ThrowDPCSimilarityException();
3946 ty_image=SIMDerivativeImage(target_image,"Sobel:90",exception);
3947 target_image=DestroyImage(target_image);
3948 if (ty_image == (Image *) NULL)
3949 ThrowDPCSimilarityException();
3950 /*
3951 Compute magnitude of derivatives.
3952 */
3953 magnitude_image=SIMMagnitudeImage(tx_image,ty_image,exception);
3954 if (magnitude_image == (Image *) NULL)
3955 ThrowDPCSimilarityException();
3956 /*
3957 Divide Lx and Ly by magnitude.
3958 */
3959 trx_image=SIMDivideByMagnitude(tx_image,magnitude_image,image,exception);
3960 tx_image=DestroyImage(tx_image);
3961 if (trx_image == (Image *) NULL)
3962 ThrowDPCSimilarityException();
3963 tx_image=trx_image;
3964 try_image=SIMDivideByMagnitude(ty_image,magnitude_image,image,exception);
3965 ty_image=DestroyImage(ty_image);
3966 magnitude_image=DestroyImage(magnitude_image);
3967 if (try_image == (Image *) NULL)
3968 ThrowDPCSimilarityException();
3969 ty_image=try_image;
3970 /*
3971 Compute the cross correlation of the test and reference images.
3972 */
3973 trx_image=SIMCrossCorrelationImage(tx_image,rx_image,exception);
3974 rx_image=DestroyImage(rx_image);
3975 tx_image=DestroyImage(tx_image);
3976 if (trx_image == (Image *) NULL)
3977 ThrowDPCSimilarityException();
3978 try_image=SIMCrossCorrelationImage(ty_image,ry_image,exception);
3979 ry_image=DestroyImage(ry_image);
3980 ty_image=DestroyImage(ty_image);
3981 if (try_image == (Image *) NULL)
3982 ThrowDPCSimilarityException();
3983 /*
3984 Evaluate dot product correlation image.
3985 */
3986 (void) SetImageArtifact(try_image,"compose:clamp","false");
3987 status=CompositeImage(trx_image,try_image,PlusCompositeOp,MagickTrue,0,0,
3988 exception);
3989 try_image=DestroyImage(try_image);
3990 if (status == MagickFalse)
3991 ThrowDPCSimilarityException();
3992 status=SIMMultiplyImage(trx_image,edge_factor,
3993 (const ChannelStatistics *) NULL,exception);
3994 if (status == MagickFalse)
3995 ThrowDPCSimilarityException();
3996 /*
3997 Crop results.
3998 */
3999 SetGeometry(image,&geometry);
4000 geometry.width=image->columns;
4001 geometry.height=image->rows;
4002 (void) ResetImagePage(trx_image,"0x0+0+0");
4003 dot_product_image=CropImage(trx_image,&geometry,exception);
4004 trx_image=DestroyImage(trx_image);
4005 if (dot_product_image == (Image *) NULL)
4006 ThrowDPCSimilarityException();
4007 (void) ResetImagePage(dot_product_image,"0x0+0+0");
4008 /*
4009 Identify the maxima value in the image and its location.
4010 */
4011 status=GrayscaleImage(dot_product_image,AveragePixelIntensityMethod,
4012 exception);
4013 if (status == MagickFalse)
4014 ThrowDPCSimilarityException();
4015 dot_product_image->depth=32;
4016 dot_product_image->colorspace=GRAYColorspace;
4017 dot_product_image->alpha_trait=UndefinedPixelTrait;
4018 status=SIMFilterImageNaNs(dot_product_image,exception);
4019 if (status == MagickFalse)
4020 ThrowDPCSimilarityException();
4021 status=SIMMaximaImage(dot_product_image,&maxima,offset,exception);
4022 if (status == MagickFalse)
4023 ThrowDPCSimilarityException();
4024 if ((QuantumScale*maxima) > 1.0)
4025 {
4026 status=SIMMultiplyImage(dot_product_image,1.0/(QuantumScale*maxima),
4027 (const ChannelStatistics *) NULL,exception);
4028 maxima=(double) QuantumRange;
4029 }
4030 *similarity_metric=QuantumScale*maxima;
4031 return(dot_product_image);
4032}
4033
4034static Image *MSESimilarityImage(const Image *image,const Image *reconstruct,
4035 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4036{
4037#define ThrowMSESimilarityException() \
4038{ \
4039 if (alpha_image != (Image *) NULL) \
4040 alpha_image=DestroyImage(alpha_image); \
4041 if (beta_image != (Image *) NULL) \
4042 beta_image=DestroyImage(beta_image); \
4043 if (channel_statistics != (ChannelStatistics *) NULL) \
4044 channel_statistics=(ChannelStatistics *) \
4045 RelinquishMagickMemory(channel_statistics); \
4046 if (mean_image != (Image *) NULL) \
4047 mean_image=DestroyImage(mean_image); \
4048 if (mse_image != (Image *) NULL) \
4049 mse_image=DestroyImage(mse_image); \
4050 if (reconstruct_image != (Image *) NULL) \
4051 reconstruct_image=DestroyImage(reconstruct_image); \
4052 if (sum_image != (Image *) NULL) \
4053 sum_image=DestroyImage(sum_image); \
4054 if (alpha_image != (Image *) NULL) \
4055 alpha_image=DestroyImage(alpha_image); \
4056 return((Image *) NULL); \
4057}
4058
4059 ChannelStatistics
4060 *channel_statistics = (ChannelStatistics *) NULL;
4061
4062 double
4063 minima = 0.0;
4064
4065 Image
4066 *alpha_image = (Image *) NULL,
4067 *beta_image = (Image *) NULL,
4068 *mean_image = (Image *) NULL,
4069 *mse_image = (Image *) NULL,
4070 *reconstruct_image = (Image *) NULL,
4071 *sum_image = (Image *) NULL,
4072 *target_image = (Image *) NULL;
4073
4074 MagickBooleanType
4075 status = MagickTrue;
4076
4077 RectangleInfo
4078 geometry;
4079
4080 /*
4081 MSE correlation-based image similarity using FFT local statistics.
4082 */
4083 target_image=SIMSquareImage(image,exception);
4084 if (target_image == (Image *) NULL)
4085 ThrowMSESimilarityException();
4086 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
4087 if (reconstruct_image == (Image *) NULL)
4088 ThrowMSESimilarityException();
4089 /*
4090 Create (U * test)/# pixels.
4091 */
4092 alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
4093 exception);
4094 target_image=DestroyImage(target_image);
4095 if (alpha_image == (Image *) NULL)
4096 ThrowMSESimilarityException();
4097 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(double)
4098 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
4099 if (status == MagickFalse)
4100 ThrowMSESimilarityException();
4101 /*
4102 Create 2*(test * reconstruction)# pixels.
4103 */
4104 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
4105 MagickTrue,0,0,exception);
4106 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
4107 reconstruct_image=DestroyImage(reconstruct_image);
4108 if (beta_image == (Image *) NULL)
4109 ThrowMSESimilarityException();
4110 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(double)
4111 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
4112 if (status == MagickFalse)
4113 ThrowMSESimilarityException();
4114 /*
4115 Mean of reconstruction squared.
4116 */
4117 sum_image=SIMSquareImage(reconstruct,exception);
4118 if (sum_image == (Image *) NULL)
4119 ThrowMSESimilarityException();
4120 channel_statistics=GetImageStatistics(sum_image,exception);
4121 if (channel_statistics == (ChannelStatistics *) NULL)
4122 ThrowMSESimilarityException();
4123 status=SetImageStorageClass(sum_image,DirectClass,exception);
4124 if (status == MagickFalse)
4125 ThrowMSESimilarityException();
4126 status=SIMSetImageMean(sum_image,channel_statistics,exception);
4127 channel_statistics=(ChannelStatistics *)
4128 RelinquishMagickMemory(channel_statistics);
4129 if (status == MagickFalse)
4130 ThrowMSESimilarityException();
4131 /*
4132 Create mean image.
4133 */
4134 AppendImageToList(&sum_image,alpha_image);
4135 AppendImageToList(&sum_image,beta_image);
4136 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
4137 if (mean_image == (Image *) NULL)
4138 ThrowMSESimilarityException();
4139 sum_image=DestroyImage(sum_image);
4140 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
4141 if (status == MagickFalse)
4142 ThrowMSESimilarityException();
4143 /*
4144 Crop to difference of reconstruction and test images.
4145 */
4146 SetGeometry(image,&geometry);
4147 geometry.width=image->columns;
4148 geometry.height=image->rows;
4149 (void) ResetImagePage(mean_image,"0x0+0+0");
4150 mse_image=CropImage(mean_image,&geometry,exception);
4151 mean_image=DestroyImage(mean_image);
4152 if (mse_image == (Image *) NULL)
4153 ThrowMSESimilarityException();
4154 /*
4155 Identify the minima value in the correlation image and its location.
4156 */
4157 (void) ResetImagePage(mse_image,"0x0+0+0");
4158 (void) ClampImage(mse_image,exception);
4159 mse_image->depth=32;
4160 mse_image->colorspace=GRAYColorspace;
4161 mse_image->alpha_trait=UndefinedPixelTrait;
4162 status=SIMMinimaImage(mse_image,&minima,offset,exception);
4163 if (status == MagickFalse)
4164 ThrowMSESimilarityException();
4165 status=NegateImage(mse_image,MagickFalse,exception);
4166 if (status == MagickFalse)
4167 ThrowMSESimilarityException();
4168 alpha_image=DestroyImage(alpha_image);
4169 beta_image=DestroyImage(beta_image);
4170 if ((QuantumScale*minima) < FLT_EPSILON)
4171 minima=0.0;
4172 *similarity_metric=QuantumScale*minima;
4173 return(mse_image);
4174}
4175
4176static Image *NCCSimilarityImage(const Image *image,const Image *reconstruct,
4177 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4178{
4179#define ThrowNCCSimilarityException() \
4180{ \
4181 if (alpha_image != (Image *) NULL) \
4182 alpha_image=DestroyImage(alpha_image); \
4183 if (beta_image != (Image *) NULL) \
4184 beta_image=DestroyImage(beta_image); \
4185 if (channel_statistics != (ChannelStatistics *) NULL) \
4186 channel_statistics=(ChannelStatistics *) \
4187 RelinquishMagickMemory(channel_statistics); \
4188 if (correlation_image != (Image *) NULL) \
4189 correlation_image=DestroyImage(correlation_image); \
4190 if (divide_image != (Image *) NULL) \
4191 divide_image=DestroyImage(divide_image); \
4192 if (ncc_image != (Image *) NULL) \
4193 ncc_image=DestroyImage(ncc_image); \
4194 if (normalize_image != (Image *) NULL) \
4195 normalize_image=DestroyImage(normalize_image); \
4196 if (reconstruct_image != (Image *) NULL) \
4197 reconstruct_image=DestroyImage(reconstruct_image); \
4198 if (target_image != (Image *) NULL) \
4199 target_image=DestroyImage(target_image); \
4200 if (variance_image != (Image *) NULL) \
4201 variance_image=DestroyImage(variance_image); \
4202 return((Image *) NULL); \
4203}
4204
4205 ChannelStatistics
4206 *channel_statistics = (ChannelStatistics *) NULL;
4207
4208 double
4209 maxima = 0.0;
4210
4211 Image
4212 *alpha_image = (Image *) NULL,
4213 *beta_image = (Image *) NULL,
4214 *correlation_image = (Image *) NULL,
4215 *divide_image = (Image *) NULL,
4216 *ncc_image = (Image *) NULL,
4217 *normalize_image = (Image *) NULL,
4218 *reconstruct_image = (Image *) NULL,
4219 *target_image = (Image *) NULL,
4220 *variance_image = (Image *) NULL;
4221
4222 MagickBooleanType
4223 status = MagickTrue;
4224
4225 RectangleInfo
4226 geometry;
4227
4228 /*
4229 NCC correlation-based image similarity with FFT local statistics.
4230 */
4231 target_image=SIMSquareImage(image,exception);
4232 if (target_image == (Image *) NULL)
4233 ThrowNCCSimilarityException();
4234 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
4235 if (reconstruct_image == (Image *) NULL)
4236 ThrowNCCSimilarityException();
4237 /*
4238 Compute the cross correlation of the test and reconstruction images.
4239 */
4240 alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
4241 exception);
4242 target_image=DestroyImage(target_image);
4243 if (alpha_image == (Image *) NULL)
4244 ThrowNCCSimilarityException();
4245 status=SIMMultiplyImage(alpha_image,(double) QuantumRange*
4246 reconstruct->columns*reconstruct->rows,(const ChannelStatistics *) NULL,
4247 exception);
4248 if (status == MagickFalse)
4249 ThrowNCCSimilarityException();
4250 /*
4251 Compute the cross correlation of the source and reconstruction images.
4252 */
4253 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
4254 reconstruct_image=DestroyImage(reconstruct_image);
4255 if (beta_image == (Image *) NULL)
4256 ThrowNCCSimilarityException();
4257 target_image=SIMSquareImage(beta_image,exception);
4258 beta_image=DestroyImage(beta_image);
4259 if (target_image == (Image *) NULL)
4260 ThrowNCCSimilarityException();
4261 status=SIMMultiplyImage(target_image,(double) QuantumRange,
4262 (const ChannelStatistics *) NULL,exception);
4263 if (status == MagickFalse)
4264 ThrowNCCSimilarityException();
4265 /*
4266 Compute the variance of the two images.
4267 */
4268 variance_image=SIMVarianceImage(alpha_image,target_image,exception);
4269 target_image=DestroyImage(target_image);
4270 alpha_image=DestroyImage(alpha_image);
4271 if (variance_image == (Image *) NULL)
4272 ThrowNCCSimilarityException();
4273 /*
4274 Subtract the image mean.
4275 */
4276 channel_statistics=GetImageStatistics(reconstruct,exception);
4277 if (channel_statistics == (ChannelStatistics *) NULL)
4278 ThrowNCCSimilarityException();
4279 status=SIMMultiplyImage(variance_image,1.0,channel_statistics,exception);
4280 if (status == MagickFalse)
4281 ThrowNCCSimilarityException();
4282 normalize_image=SIMSubtractImageMean(image,reconstruct,channel_statistics,
4283 exception);
4284 channel_statistics=(ChannelStatistics *)
4285 RelinquishMagickMemory(channel_statistics);
4286 if (normalize_image == (Image *) NULL)
4287 ThrowNCCSimilarityException();
4288 correlation_image=SIMCrossCorrelationImage(image,normalize_image,exception);
4289 normalize_image=DestroyImage(normalize_image);
4290 if (correlation_image == (Image *) NULL)
4291 ThrowNCCSimilarityException();
4292 /*
4293 Divide the two images.
4294 */
4295 divide_image=SIMDivideImage(correlation_image,variance_image,exception);
4296 correlation_image=DestroyImage(correlation_image);
4297 variance_image=DestroyImage(variance_image);
4298 if (divide_image == (Image *) NULL)
4299 ThrowNCCSimilarityException();
4300 /*
4301 Crop padding.
4302 */
4303 SetGeometry(image,&geometry);
4304 geometry.width=image->columns;
4305 geometry.height=image->rows;
4306 (void) ResetImagePage(divide_image,"0x0+0+0");
4307 ncc_image=CropImage(divide_image,&geometry,exception);
4308 divide_image=DestroyImage(divide_image);
4309 if (ncc_image == (Image *) NULL)
4310 ThrowNCCSimilarityException();
4311 /*
4312 Identify the maxima value in the image and its location.
4313 */
4314 (void) ResetImagePage(ncc_image,"0x0+0+0");
4315 status=GrayscaleImage(ncc_image,AveragePixelIntensityMethod,exception);
4316 if (status == MagickFalse)
4317 ThrowNCCSimilarityException();
4318 ncc_image->depth=32;
4319 ncc_image->colorspace=GRAYColorspace;
4320 ncc_image->alpha_trait=UndefinedPixelTrait;
4321 status=SIMMaximaImage(ncc_image,&maxima,offset,exception);
4322 if (status == MagickFalse)
4323 ThrowNCCSimilarityException();
4324 if ((QuantumScale*maxima) > 1.0)
4325 {
4326 status=SIMMultiplyImage(ncc_image,1.0/(QuantumScale*maxima),
4327 (const ChannelStatistics *) NULL,exception);
4328 maxima=(double) QuantumRange;
4329 }
4330 *similarity_metric=QuantumScale*maxima;
4331 return(ncc_image);
4332}
4333
4334static Image *PhaseSimilarityImage(const Image *image,const Image *reconstruct,
4335 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4336{
4337#define ThrowPhaseSimilarityException() \
4338{ \
4339 if (correlation_image != (Image *) NULL) \
4340 correlation_image=DestroyImage(correlation_image); \
4341 if (fft_images != (Image *) NULL) \
4342 fft_images=DestroyImageList(fft_images); \
4343 if (gamma_image != (Image *) NULL) \
4344 gamma_image=DestroyImage(gamma_image); \
4345 if (magnitude_image != (Image *) NULL) \
4346 magnitude_image=DestroyImage(magnitude_image); \
4347 if (phase_image != (Image *) NULL) \
4348 phase_image=DestroyImage(phase_image); \
4349 if (reconstruct_image != (Image *) NULL) \
4350 reconstruct_image=DestroyImage(reconstruct_image); \
4351 if (reconstruct_magnitude != (Image *) NULL) \
4352 reconstruct_magnitude=DestroyImage(reconstruct_magnitude); \
4353 if (target_image != (Image *) NULL) \
4354 target_image=DestroyImage(target_image); \
4355 if (test_magnitude != (Image *) NULL) \
4356 test_magnitude=DestroyImage(test_magnitude); \
4357 return((Image *) NULL); \
4358}
4359
4360 double
4361 maxima = 0.0;
4362
4363 Image
4364 *correlation_image = (Image *) NULL,
4365 *fft_images = (Image *) NULL,
4366 *gamma_image = (Image *) NULL,
4367 *magnitude_image = (Image *) NULL,
4368 *phase_image = (Image *) NULL,
4369 *reconstruct_image = (Image *) NULL,
4370 *reconstruct_magnitude = (Image *) NULL,
4371 *target_image = (Image *) NULL,
4372 *test_magnitude = (Image *) NULL;
4373
4374 MagickBooleanType
4375 status = MagickTrue;
4376
4377 RectangleInfo
4378 geometry;
4379
4380 /*
4381 Phase correlation-based image similarity using FFT local statistics.
4382 */
4383 target_image=CloneImage(image,0,0,MagickTrue,exception);
4384 if (target_image == (Image *) NULL)
4385 ThrowPhaseSimilarityException();
4386 (void) ResetImagePage(target_image,"0x0+0+0");
4387 GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
4388 &target_image->background_color);
4389 status=SetImageExtent(target_image,2*(size_t) ceil((double) image->columns/
4390 2.0),2*(size_t) ceil((double) image->rows/2.0),exception);
4391 if (status == MagickFalse)
4392 ThrowPhaseSimilarityException();
4393 /*
4394 Compute the cross correlation of the test and reconstruct magnitudes.
4395 */
4396 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
4397 if (reconstruct_image == (Image *) NULL)
4398 ThrowPhaseSimilarityException();
4399 (void) ResetImagePage(reconstruct_image,"0x0+0+0");
4400 GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
4401 &reconstruct_image->background_color);
4402 status=SetImageExtent(reconstruct_image,2*(size_t) ceil((double)
4403 image->columns/2.0),2*(size_t) ceil((double) image->rows/2.0),exception);
4404 if (status == MagickFalse)
4405 ThrowPhaseSimilarityException();
4406 /*
4407 Evaluate phase coorelation image and divide by the product magnitude.
4408 */
4409 (void) SetImageArtifact(target_image,"fourier:normalize","inverse");
4410 fft_images=ForwardFourierTransformImage(target_image,MagickTrue,exception);
4411 if (fft_images == (Image *) NULL)
4412 ThrowPhaseSimilarityException();
4413 test_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
4414 fft_images=DestroyImageList(fft_images);
4415 if (test_magnitude == (Image *) NULL)
4416 ThrowPhaseSimilarityException();
4417 (void) SetImageArtifact(reconstruct_image,"fourier:normalize","inverse");
4418 fft_images=ForwardFourierTransformImage(reconstruct_image,MagickTrue,
4419 exception);
4420 if (fft_images == (Image *) NULL)
4421 ThrowPhaseSimilarityException();
4422 reconstruct_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
4423 fft_images=DestroyImageList(fft_images);
4424 if (reconstruct_magnitude == (Image *) NULL)
4425 ThrowPhaseSimilarityException();
4426 magnitude_image=CloneImage(reconstruct_magnitude,0,0,MagickTrue,exception);
4427 if (magnitude_image == (Image *) NULL)
4428 ThrowPhaseSimilarityException();
4429 DisableCompositeClampUnlessSpecified(magnitude_image);
4430 (void) CompositeImage(magnitude_image,test_magnitude,MultiplyCompositeOp,
4431 MagickTrue,0,0,exception);
4432 /*
4433 Compute the cross correlation of the test and reconstruction images.
4434 */
4435 correlation_image=SIMPhaseCorrelationImage(target_image,reconstruct_image,
4436 magnitude_image,exception);
4437 target_image=DestroyImage(target_image);
4438 reconstruct_image=DestroyImage(reconstruct_image);
4439 test_magnitude=DestroyImage(test_magnitude);
4440 reconstruct_magnitude=DestroyImage(reconstruct_magnitude);
4441 if (correlation_image == (Image *) NULL)
4442 ThrowPhaseSimilarityException();
4443 /*
4444 Identify the maxima value in the image and its location.
4445 */
4446 gamma_image=CloneImage(correlation_image,0,0,MagickTrue,exception);
4447 correlation_image=DestroyImage(correlation_image);
4448 if (gamma_image == (Image *) NULL)
4449 ThrowPhaseSimilarityException();
4450 /*
4451 Crop padding.
4452 */
4453 SetGeometry(image,&geometry);
4454 geometry.width=image->columns;
4455 geometry.height=image->rows;
4456 (void) ResetImagePage(gamma_image,"0x0+0+0");
4457 phase_image=CropImage(gamma_image,&geometry,exception);
4458 gamma_image=DestroyImage(gamma_image);
4459 if (phase_image == (Image *) NULL)
4460 ThrowPhaseSimilarityException();
4461 (void) ResetImagePage(phase_image,"0x0+0+0");
4462 /*
4463 Identify the maxima value in the correlation image and its location.
4464 */
4465 status=GrayscaleImage(phase_image,AveragePixelIntensityMethod,exception);
4466 if (status == MagickFalse)
4467 ThrowPhaseSimilarityException();
4468 phase_image->depth=32;
4469 phase_image->colorspace=GRAYColorspace;
4470 phase_image->alpha_trait=UndefinedPixelTrait;
4471 status=SIMFilterImageNaNs(phase_image,exception);
4472 if (status == MagickFalse)
4473 ThrowPhaseSimilarityException();
4474 status=SIMMaximaImage(phase_image,&maxima,offset,exception);
4475 if (status == MagickFalse)
4476 ThrowPhaseSimilarityException();
4477 magnitude_image=DestroyImage(magnitude_image);
4478 if ((QuantumScale*maxima) > 1.0)
4479 {
4480 status=SIMMultiplyImage(phase_image,1.0/(QuantumScale*maxima),
4481 (const ChannelStatistics *) NULL,exception);
4482 maxima=(double) QuantumRange;
4483 }
4484 *similarity_metric=QuantumScale*maxima;
4485 return(phase_image);
4486}
4487
4488static Image *PSNRSimilarityImage(const Image *image,const Image *reconstruct,
4489 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4490{
4491 Image
4492 *psnr_image = (Image *) NULL;
4493
4494 psnr_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
4495 exception);
4496 if (psnr_image == (Image *) NULL)
4497 return(psnr_image);
4498 *similarity_metric=10.0*MagickSafeLog10(MagickSafeReciprocal(
4499 *similarity_metric))/MagickSafePSNRRecipicol(10.0);
4500 return(psnr_image);
4501}
4502
4503static Image *RMSESimilarityImage(const Image *image,const Image *reconstruct,
4504 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4505{
4506 Image
4507 *rmse_image = (Image *) NULL;
4508
4509 rmse_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
4510 exception);
4511 if (rmse_image == (Image *) NULL)
4512 return(rmse_image);
4513 *similarity_metric=sqrt(*similarity_metric);
4514 return(rmse_image);
4515}
4516#endif
4517
4518static double GetSimilarityMetric(const Image *image,
4519 const Image *reconstruct_image,const MetricType metric,
4520 const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo *exception)
4521{
4522 double
4523 *channel_similarity,
4524 similarity = 0.0;
4525
4526 ExceptionInfo
4527 *sans_exception = AcquireExceptionInfo();
4528
4529 Image
4530 *similarity_image;
4531
4532 MagickBooleanType
4533 status = MagickTrue;
4534
4535 RectangleInfo
4536 geometry;
4537
4538 size_t
4539 length = MaxPixelChannels+1UL;
4540
4541 SetGeometry(reconstruct_image,&geometry);
4542 geometry.x=x_offset;
4543 geometry.y=y_offset;
4544 similarity_image=CropImage(image,&geometry,sans_exception);
4545 sans_exception=DestroyExceptionInfo(sans_exception);
4546 if (similarity_image == (Image *) NULL)
4547 return(NAN);
4548 /*
4549 Get image distortion.
4550 */
4551 channel_similarity=(double *) AcquireQuantumMemory(length,
4552 sizeof(*channel_similarity));
4553 if (channel_similarity == (double *) NULL)
4554 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
4555 (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
4556 switch (metric)
4557 {
4558 case AbsoluteErrorMetric:
4559 {
4560 status=GetAESimilarity(similarity_image,reconstruct_image,
4561 channel_similarity,exception);
4562 break;
4563 }
4564 case DotProductCorrelationErrorMetric:
4565 {
4566 status=GetDPCSimilarity(similarity_image,reconstruct_image,
4567 channel_similarity,exception);
4568 break;
4569 }
4570 case FuzzErrorMetric:
4571 {
4572 status=GetFUZZSimilarity(similarity_image,reconstruct_image,
4573 channel_similarity,exception);
4574 break;
4575 }
4576 case MeanAbsoluteErrorMetric:
4577 {
4578 status=GetMAESimilarity(similarity_image,reconstruct_image,
4579 channel_similarity,exception);
4580 break;
4581 }
4582 case MeanErrorPerPixelErrorMetric:
4583 {
4584 status=GetMEPPSimilarity(similarity_image,reconstruct_image,
4585 channel_similarity,exception);
4586 break;
4587 }
4588 case MeanSquaredErrorMetric:
4589 {
4590 status=GetMSESimilarity(similarity_image,reconstruct_image,
4591 channel_similarity,exception);
4592 break;
4593 }
4594 case NormalizedCrossCorrelationErrorMetric:
4595 {
4596 status=GetNCCSimilarity(similarity_image,reconstruct_image,
4597 channel_similarity,exception);
4598 break;
4599 }
4600 case PeakAbsoluteErrorMetric:
4601 {
4602 status=GetPASimilarity(similarity_image,reconstruct_image,
4603 channel_similarity,exception);
4604 break;
4605 }
4606 case PeakSignalToNoiseRatioErrorMetric:
4607 {
4608 status=GetPSNRSimilarity(similarity_image,reconstruct_image,
4609 channel_similarity,exception);
4610 break;
4611 }
4612 case PerceptualHashErrorMetric:
4613 {
4614 status=GetPHASHSimilarity(similarity_image,reconstruct_image,
4615 channel_similarity,exception);
4616 break;
4617 }
4618 case PhaseCorrelationErrorMetric:
4619 {
4620 status=GetPHASESimilarity(similarity_image,reconstruct_image,
4621 channel_similarity,exception);
4622 break;
4623 }
4624 case PixelDifferenceCountErrorMetric:
4625 {
4626 status=GetPDCSimilarity(similarity_image,reconstruct_image,
4627 channel_similarity,exception);
4628 break;
4629 }
4630 case RootMeanSquaredErrorMetric:
4631 case UndefinedErrorMetric:
4632 default:
4633 {
4634 status=GetRMSESimilarity(similarity_image,reconstruct_image,
4635 channel_similarity,exception);
4636 break;
4637 }
4638 case StructuralDissimilarityErrorMetric:
4639 {
4640 status=GetDSSIMSimilarity(similarity_image,reconstruct_image,
4641 channel_similarity,exception);
4642 break;
4643 }
4644 case StructuralSimilarityErrorMetric:
4645 {
4646 status=GetSSIMSimularity(similarity_image,reconstruct_image,
4647 channel_similarity,exception);
4648 break;
4649 }
4650 }
4651 similarity_image=DestroyImage(similarity_image);
4652 similarity=channel_similarity[CompositePixelChannel];
4653 channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
4654 if (status == MagickFalse)
4655 return(NAN);
4656 return(similarity);
4657}
4658
4659MagickExport Image *SimilarityImage(const Image *image,const Image *reconstruct,
4660 const MetricType metric,const double similarity_threshold,
4661 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4662{
4663#define SimilarityImageTag "Similarity/Image"
4664
4665 typedef struct
4666 {
4667 double
4668 similarity;
4669
4670 ssize_t
4671 x,
4672 y;
4673 } SimilarityInfo;
4674
4675 CacheView
4676 *similarity_view;
4677
4678 Image
4679 *similarity_image = (Image *) NULL;
4680
4681 MagickBooleanType
4682 status = MagickTrue;
4683
4684 MagickOffsetType
4685 progress = 0;
4686
4687 SimilarityInfo
4688 similarity_info = { 0.0, 0, 0 };
4689
4690 ssize_t
4691 y;
4692
4693 assert(image != (const Image *) NULL);
4694 assert(image->signature == MagickCoreSignature);
4695 assert(exception != (ExceptionInfo *) NULL);
4696 assert(exception->signature == MagickCoreSignature);
4697 assert(offset != (RectangleInfo *) NULL);
4698 if (IsEventLogging() != MagickFalse)
4699 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4700 SetGeometry(reconstruct,offset);
4701 *similarity_metric=0.0;
4702 offset->x=0;
4703 offset->y=0;
4704#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
4705{
4706 const char *artifact = GetImageArtifact(image,"compare:frequency-domain");
4707 if (artifact == (const char *) NULL)
4708 artifact=GetImageArtifact(image,"compare:accelerate-ncc");
4709 if (((artifact == (const char *) NULL) ||
4710 (IsStringTrue(artifact) != MagickFalse)) &&
4711 ((image->channels & ReadMaskChannel) == 0))
4712 switch (metric)
4713 {
4714 case DotProductCorrelationErrorMetric:
4715 {
4716 similarity_image=DPCSimilarityImage(image,reconstruct,offset,
4717 similarity_metric,exception);
4718 return(similarity_image);
4719 }
4720 case MeanSquaredErrorMetric:
4721 {
4722 similarity_image=MSESimilarityImage(image,reconstruct,offset,
4723 similarity_metric,exception);
4724 return(similarity_image);
4725 }
4726 case NormalizedCrossCorrelationErrorMetric:
4727 {
4728 similarity_image=NCCSimilarityImage(image,reconstruct,offset,
4729 similarity_metric,exception);
4730 return(similarity_image);
4731 }
4732 case PeakSignalToNoiseRatioErrorMetric:
4733 {
4734 similarity_image=PSNRSimilarityImage(image,reconstruct,offset,
4735 similarity_metric,exception);
4736 return(similarity_image);
4737 }
4738 case PhaseCorrelationErrorMetric:
4739 {
4740 similarity_image=PhaseSimilarityImage(image,reconstruct,offset,
4741 similarity_metric,exception);
4742 return(similarity_image);
4743 }
4744 case RootMeanSquaredErrorMetric:
4745 case UndefinedErrorMetric:
4746 {
4747 similarity_image=RMSESimilarityImage(image,reconstruct,offset,
4748 similarity_metric,exception);
4749 return(similarity_image);
4750 }
4751 default:
4752 break;
4753 }
4754}
4755#endif
4756 if ((image->columns < reconstruct->columns) ||
4757 (image->rows < reconstruct->rows))
4758 {
4759 (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
4760 "GeometryDoesNotContainImage","`%s'",image->filename);
4761 return((Image *) NULL);
4762 }
4763 similarity_image=CloneImage(image,image->columns-reconstruct->columns+1,
4764 image->rows-reconstruct->rows+1,MagickTrue,exception);
4765 if (similarity_image == (Image *) NULL)
4766 return((Image *) NULL);
4767 similarity_image->depth=32;
4768 similarity_image->colorspace=GRAYColorspace;
4769 similarity_image->alpha_trait=UndefinedPixelTrait;
4770 status=SetImageStorageClass(similarity_image,DirectClass,exception);
4771 if (status == MagickFalse)
4772 return(DestroyImage(similarity_image));
4773 /*
4774 Measure similarity of reconstruction image against image.
4775 */
4776 similarity_info.similarity=GetSimilarityMetric(image,reconstruct,metric,
4777 similarity_info.x,similarity_info.y,exception);
4778 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
4779#if defined(MAGICKCORE_OPENMP_SUPPORT)
4780 #pragma omp parallel for schedule(static) shared(similarity_info,status) \
4781 magick_number_threads(image,reconstruct,similarity_image->rows,1)
4782#endif
4783 for (y=0; y < (ssize_t) similarity_image->rows; y++)
4784 {
4785 double
4786 similarity;
4787
4788 MagickBooleanType
4789 threshold_trigger = MagickFalse;
4790
4791 Quantum
4792 *magick_restrict q;
4793
4794 SimilarityInfo
4795 channel_info = similarity_info;
4796
4797 ssize_t
4798 x;
4799
4800 if (status == MagickFalse)
4801 continue;
4802 if (threshold_trigger != MagickFalse)
4803 continue;
4804 q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
4805 similarity_image->columns,1,exception);
4806 if (q == (Quantum *) NULL)
4807 {
4808 status=MagickFalse;
4809 continue;
4810 }
4811 for (x=0; x < (ssize_t) similarity_image->columns; x++)
4812 {
4813 ssize_t
4814 i;
4815
4816 similarity=GetSimilarityMetric((Image *) image,reconstruct,metric,x,y,
4817 exception);
4818 switch (metric)
4819 {
4820 case DotProductCorrelationErrorMetric:
4821 case NormalizedCrossCorrelationErrorMetric:
4822 case PeakSignalToNoiseRatioErrorMetric:
4823 case PhaseCorrelationErrorMetric:
4824 case StructuralSimilarityErrorMetric:
4825 {
4826 if (similarity <= channel_info.similarity)
4827 break;
4828 channel_info.similarity=similarity;
4829 channel_info.x=x;
4830 channel_info.y=y;
4831 break;
4832 }
4833 default:
4834 {
4835 if (similarity >= channel_info.similarity)
4836 break;
4837 channel_info.similarity=similarity;
4838 channel_info.x=x;
4839 channel_info.y=y;
4840 break;
4841 }
4842 }
4843 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4844 {
4845 PixelChannel channel = GetPixelChannelChannel(image,i);
4846 PixelTrait traits = GetPixelChannelTraits(image,channel);
4847 PixelTrait similarity_traits = GetPixelChannelTraits(similarity_image,
4848 channel);
4849 if (((traits & UpdatePixelTrait) == 0) ||
4850 ((similarity_traits & UpdatePixelTrait) == 0))
4851 continue;
4852 switch (metric)
4853 {
4854 case DotProductCorrelationErrorMetric:
4855 case NormalizedCrossCorrelationErrorMetric:
4856 case PeakSignalToNoiseRatioErrorMetric:
4857 case PhaseCorrelationErrorMetric:
4858 case StructuralSimilarityErrorMetric:
4859 {
4860 SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4861 QuantumRange*similarity),q);
4862 break;
4863 }
4864 default:
4865 {
4866 SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4867 QuantumRange*(1.0-similarity)),q);
4868 break;
4869 }
4870 }
4871 }
4872 q+=(ptrdiff_t) GetPixelChannels(similarity_image);
4873 }
4874#if defined(MAGICKCORE_OPENMP_SUPPORT)
4875 #pragma omp critical (MagickCore_GetSimilarityMetric)
4876#endif
4877 switch (metric)
4878 {
4879 case DotProductCorrelationErrorMetric:
4880 case NormalizedCrossCorrelationErrorMetric:
4881 case PeakSignalToNoiseRatioErrorMetric:
4882 case PhaseCorrelationErrorMetric:
4883 case StructuralSimilarityErrorMetric:
4884 {
4885 if (similarity_threshold != DefaultSimilarityThreshold)
4886 if (channel_info.similarity >= similarity_threshold)
4887 threshold_trigger=MagickTrue;
4888 if (channel_info.similarity >= similarity_info.similarity)
4889 similarity_info=channel_info;
4890 break;
4891 }
4892 default:
4893 {
4894 if (similarity_threshold != DefaultSimilarityThreshold)
4895 if (channel_info.similarity < similarity_threshold)
4896 threshold_trigger=MagickTrue;
4897 if (channel_info.similarity < similarity_info.similarity)
4898 similarity_info=channel_info;
4899 break;
4900 }
4901 }
4902 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
4903 status=MagickFalse;
4904 if (image->progress_monitor != (MagickProgressMonitor) NULL)
4905 {
4906 MagickBooleanType
4907 proceed;
4908
4909 progress++;
4910 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
4911 if (proceed == MagickFalse)
4912 status=MagickFalse;
4913 }
4914 }
4915 similarity_view=DestroyCacheView(similarity_view);
4916 if (status == MagickFalse)
4917 similarity_image=DestroyImage(similarity_image);
4918 *similarity_metric=similarity_info.similarity;
4919 if (fabs(*similarity_metric) < MagickEpsilon)
4920 *similarity_metric=0.0;
4921 offset->x=similarity_info.x;
4922 offset->y=similarity_info.y;
4923 (void) FormatImageProperty((Image *) image,"similarity","%.*g",
4924 GetMagickPrecision(),*similarity_metric);
4925 (void) FormatImageProperty((Image *) image,"similarity.offset.x","%.*g",
4926 GetMagickPrecision(),(double) offset->x);
4927 (void) FormatImageProperty((Image *) image,"similarity.offset.y","%.*g",
4928 GetMagickPrecision(),(double) offset->y);
4929 return(similarity_image);
4930}