MagickCore 7.1.2-21
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#define PHASHNormalizationFactor 389.373723242
1913
1914 ChannelPerceptualHash
1915 *channel_phash,
1916 *reconstruct_phash;
1917
1918 const char
1919 *artifact;
1920
1921 ssize_t
1922 k;
1923
1924 /*
1925 Compute the perceptual hash similarity.
1926 */
1927 channel_phash=GetImagePerceptualHash(image,exception);
1928 if (channel_phash == (ChannelPerceptualHash *) NULL)
1929 return(MagickFalse);
1930 reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1931 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1932 {
1933 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1934 channel_phash);
1935 return(MagickFalse);
1936 }
1937 for (k=0; k < MaxPixelChannels; k++)
1938 {
1939 double
1940 difference;
1941
1942 ssize_t
1943 i;
1944
1945 PixelChannel channel = GetPixelChannelChannel(image,k);
1946 PixelTrait traits = GetPixelChannelTraits(image,channel);
1947 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1948 channel);
1949 if (((traits & UpdatePixelTrait) == 0) ||
1950 ((reconstruct_traits & UpdatePixelTrait) == 0))
1951 continue;
1952 difference=0.0;
1953 for (i=0; i < MaximumNumberOfImageMoments; i++)
1954 {
1955 double
1956 alpha,
1957 beta;
1958
1959 ssize_t
1960 j;
1961
1962 for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1963 {
1964 double
1965 error;
1966
1967 alpha=channel_phash[k].phash[j][i];
1968 beta=reconstruct_phash[k].phash[j][i];
1969 error=beta-alpha;
1970 if (IsNaN(error) != 0)
1971 error=0.0;
1972 difference+=error*error/PHASHNormalizationFactor;
1973 }
1974 }
1975 similarity[k]+=difference;
1976 similarity[CompositePixelChannel]+=difference;
1977 }
1978 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
1979 artifact=GetImageArtifact(image,"phash:normalize");
1980 if (IsStringTrue(artifact) != MagickFalse)
1981 {
1982 ssize_t
1983 j;
1984
1985 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1986 {
1987 PixelChannel channel = GetPixelChannelChannel(image,j);
1988 PixelTrait traits = GetPixelChannelTraits(image,channel);
1989 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1990 channel);
1991 if (((traits & UpdatePixelTrait) == 0) ||
1992 ((reconstruct_traits & UpdatePixelTrait) == 0))
1993 continue;
1994 similarity[j]=sqrt(similarity[j]/channel_phash[0].number_colorspaces);
1995 }
1996 similarity[CompositePixelChannel]=sqrt(similarity[CompositePixelChannel]/
1997 channel_phash[0].number_colorspaces);
1998 }
1999 /*
2000 Free resources.
2001 */
2002 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
2003 reconstruct_phash);
2004 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
2005 return(MagickTrue);
2006}
2007
2008static MagickBooleanType GetRMSESimilarity(const Image *image,
2009 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2010{
2011#define RMSESquareRoot(x) sqrt((x) < 0.0 ? 0.0 : (x))
2012
2013 MagickBooleanType
2014 status = MagickTrue;
2015
2016 ssize_t
2017 i;
2018
2019 /*
2020 Compute the root mean-squared error similarity.
2021 */
2022 status=GetMSESimilarity(image,reconstruct_image,similarity,exception);
2023 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2024 {
2025 PixelChannel channel = GetPixelChannelChannel(image,i);
2026 PixelTrait traits = GetPixelChannelTraits(image,channel);
2027 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2028 channel);
2029 if (((traits & UpdatePixelTrait) == 0) ||
2030 ((reconstruct_traits & UpdatePixelTrait) == 0))
2031 continue;
2032 similarity[i]=RMSESquareRoot(similarity[i]);
2033 }
2034 similarity[CompositePixelChannel]=RMSESquareRoot(
2035 similarity[CompositePixelChannel]);
2036 return(status);
2037}
2038
2039static MagickBooleanType GetSSIMSimularity(const Image *image,
2040 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2041{
2042#define SSIMRadius 5.0
2043#define SSIMSigma 1.5
2044#define SSIMK1 0.01
2045#define SSIMK2 0.03
2046#define SSIML 1.0
2047
2048 CacheView
2049 *image_view,
2050 *reconstruct_view;
2051
2052 char
2053 geometry[MagickPathExtent];
2054
2055 const char
2056 *artifact;
2057
2058 double
2059 area = 0.0,
2060 c1,
2061 c2,
2062 radius,
2063 sigma;
2064
2065 KernelInfo
2066 *kernel_info;
2067
2068 MagickBooleanType
2069 status = MagickTrue;
2070
2071 size_t
2072 columns,
2073 rows;
2074
2075 ssize_t
2076 l,
2077 y;
2078
2079 /*
2080 Compute the structual similarity index similarity.
2081 */
2082 radius=SSIMRadius;
2083 artifact=GetImageArtifact(image,"compare:ssim-radius");
2084 if (artifact != (const char *) NULL)
2085 radius=StringToDouble(artifact,(char **) NULL);
2086 sigma=SSIMSigma;
2087 artifact=GetImageArtifact(image,"compare:ssim-sigma");
2088 if (artifact != (const char *) NULL)
2089 sigma=StringToDouble(artifact,(char **) NULL);
2090 (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
2091 radius,sigma);
2092 kernel_info=AcquireKernelInfo(geometry,exception);
2093 if (kernel_info == (KernelInfo *) NULL)
2094 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
2095 image->filename);
2096 c1=pow(SSIMK1*SSIML,2.0);
2097 artifact=GetImageArtifact(image,"compare:ssim-k1");
2098 if (artifact != (const char *) NULL)
2099 c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
2100 c2=pow(SSIMK2*SSIML,2.0);
2101 artifact=GetImageArtifact(image,"compare:ssim-k2");
2102 if (artifact != (const char *) NULL)
2103 c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
2104 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
2105 image_view=AcquireVirtualCacheView(image,exception);
2106 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2107#if defined(MAGICKCORE_OPENMP_SUPPORT)
2108 #pragma omp parallel for schedule(static) shared(area,similarity,status) \
2109 magick_number_threads(image,reconstruct_image,rows,1)
2110#endif
2111 for (y=0; y < (ssize_t) rows; y++)
2112 {
2113 const Quantum
2114 *magick_restrict p,
2115 *magick_restrict q;
2116
2117 double
2118 channel_area = 0.0,
2119 channel_similarity[MaxPixelChannels+1] = { 0.0 };
2120
2121 ssize_t
2122 i,
2123 x;
2124
2125 if (status == MagickFalse)
2126 continue;
2127 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
2128 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
2129 kernel_info->height,exception);
2130 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
2131 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
2132 kernel_info->height,exception);
2133 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
2134 {
2135 status=MagickFalse;
2136 continue;
2137 }
2138 for (x=0; x < (ssize_t) columns; x++)
2139 {
2140 const Quantum
2141 *magick_restrict reconstruct,
2142 *magick_restrict test;
2143
2144 double
2145 x_pixel_mu[MaxPixelChannels+1] = { 0.0 },
2146 x_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 },
2147 xy_sigma[MaxPixelChannels+1] = { 0.0 },
2148 y_pixel_mu[MaxPixelChannels+1] = { 0.0 },
2149 y_pixel_sigma_squared[MaxPixelChannels+1] = { 0.0 };
2150
2151 MagickRealType
2152 *k;
2153
2154 ssize_t
2155 v;
2156
2157 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
2158 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
2159 {
2160 p+=(ptrdiff_t) GetPixelChannels(image);
2161 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2162 continue;
2163 }
2164 k=kernel_info->values;
2165 test=p;
2166 reconstruct=q;
2167 for (v=0; v < (ssize_t) kernel_info->height; v++)
2168 {
2169 ssize_t
2170 u;
2171
2172 for (u=0; u < (ssize_t) kernel_info->width; u++)
2173 {
2174 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2175 {
2176 double
2177 x_pixel,
2178 y_pixel;
2179
2180 PixelChannel channel = GetPixelChannelChannel(image,i);
2181 PixelTrait traits = GetPixelChannelTraits(image,channel);
2182 PixelTrait reconstruct_traits = GetPixelChannelTraits(
2183 reconstruct_image,channel);
2184 if (((traits & UpdatePixelTrait) == 0) ||
2185 ((reconstruct_traits & UpdatePixelTrait) == 0))
2186 continue;
2187 x_pixel=QuantumScale*(double) test[i];
2188 x_pixel_mu[i]+=(*k)*x_pixel;
2189 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
2190 y_pixel=QuantumScale*(double)
2191 GetPixelChannel(reconstruct_image,channel,reconstruct);
2192 y_pixel_mu[i]+=(*k)*y_pixel;
2193 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
2194 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
2195 }
2196 k++;
2197 test+=(ptrdiff_t) GetPixelChannels(image);
2198 reconstruct+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2199 }
2200 test+=(ptrdiff_t) GetPixelChannels(image)*columns;
2201 reconstruct+=(ptrdiff_t) GetPixelChannels(reconstruct_image)*columns;
2202 }
2203 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2204 {
2205 double
2206 ssim,
2207 x_pixel_mu_squared,
2208 x_pixel_sigmas_squared,
2209 xy_mu,
2210 xy_sigmas,
2211 y_pixel_mu_squared,
2212 y_pixel_sigmas_squared;
2213
2214 PixelChannel channel = GetPixelChannelChannel(image,i);
2215 PixelTrait traits = GetPixelChannelTraits(image,channel);
2216 PixelTrait reconstruct_traits = GetPixelChannelTraits(
2217 reconstruct_image,channel);
2218 if (((traits & UpdatePixelTrait) == 0) ||
2219 ((reconstruct_traits & UpdatePixelTrait) == 0))
2220 continue;
2221 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
2222 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
2223 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
2224 xy_sigmas=xy_sigma[i]-xy_mu;
2225 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
2226 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
2227 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))*
2228 MagickSafeReciprocal((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
2229 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
2230 channel_similarity[i]+=ssim;
2231 channel_similarity[CompositePixelChannel]+=ssim;
2232 }
2233 p+=(ptrdiff_t) GetPixelChannels(image);
2234 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2235 channel_area++;
2236 }
2237#if defined(MAGICKCORE_OPENMP_SUPPORT)
2238 #pragma omp critical (MagickCore_GetSSIMSimularity)
2239#endif
2240 {
2241 ssize_t
2242 j;
2243
2244 area+=channel_area;
2245 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
2246 {
2247 PixelChannel channel = GetPixelChannelChannel(image,j);
2248 PixelTrait traits = GetPixelChannelTraits(image,channel);
2249 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2250 channel);
2251 if (((traits & UpdatePixelTrait) == 0) ||
2252 ((reconstruct_traits & UpdatePixelTrait) == 0))
2253 continue;
2254 similarity[j]+=channel_similarity[j];
2255 }
2256 similarity[CompositePixelChannel]+=
2257 channel_similarity[CompositePixelChannel];
2258 }
2259 }
2260 image_view=DestroyCacheView(image_view);
2261 reconstruct_view=DestroyCacheView(reconstruct_view);
2262 area=MagickSafeReciprocal(area);
2263 for (l=0; l < (ssize_t) GetPixelChannels(image); l++)
2264 {
2265 PixelChannel channel = GetPixelChannelChannel(image,l);
2266 PixelTrait traits = GetPixelChannelTraits(image,channel);
2267 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2268 channel);
2269 if (((traits & UpdatePixelTrait) == 0) ||
2270 ((reconstruct_traits & UpdatePixelTrait) == 0))
2271 continue;
2272 similarity[l]*=area;
2273 }
2274 similarity[CompositePixelChannel]*=area;
2275 similarity[CompositePixelChannel]/=(double) GetImageChannels(image);
2276 kernel_info=DestroyKernelInfo(kernel_info);
2277 return(status);
2278}
2279
2280static MagickBooleanType GetDSSIMSimilarity(const Image *image,
2281 const Image *reconstruct_image,double *similarity,ExceptionInfo *exception)
2282{
2283 MagickBooleanType
2284 status = MagickTrue;
2285
2286 ssize_t
2287 i;
2288
2289 /*
2290 Compute the structual dissimilarity index similarity.
2291 */
2292 status=GetSSIMSimularity(image,reconstruct_image,similarity,exception);
2293 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2294 {
2295 PixelChannel channel = GetPixelChannelChannel(image,i);
2296 PixelTrait traits = GetPixelChannelTraits(image,channel);
2297 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2298 channel);
2299 if (((traits & UpdatePixelTrait) == 0) ||
2300 ((reconstruct_traits & UpdatePixelTrait) == 0))
2301 continue;
2302 similarity[i]=(1.0-similarity[i])/2.0;
2303 }
2304 similarity[CompositePixelChannel]=(1.0-similarity[CompositePixelChannel])/2.0;
2305 return(status);
2306}
2307
2308MagickExport MagickBooleanType GetImageDistortion(Image *image,
2309 const Image *reconstruct_image,const MetricType metric,double *distortion,
2310 ExceptionInfo *exception)
2311{
2312#define CompareMetricNotSupportedException "metric not supported"
2313
2314 double
2315 *channel_similarity;
2316
2317 MagickBooleanType
2318 status = MagickTrue;
2319
2320 size_t
2321 length;
2322
2323 assert(image != (Image *) NULL);
2324 assert(image->signature == MagickCoreSignature);
2325 assert(reconstruct_image != (const Image *) NULL);
2326 assert(reconstruct_image->signature == MagickCoreSignature);
2327 assert(distortion != (double *) NULL);
2328 if (IsEventLogging() != MagickFalse)
2329 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2330 /*
2331 Get image distortion.
2332 */
2333 *distortion=0.0;
2334 length=MaxPixelChannels+1UL;
2335 channel_similarity=(double *) AcquireQuantumMemory(length,
2336 sizeof(*channel_similarity));
2337 if (channel_similarity == (double *) NULL)
2338 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2339 (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
2340 switch (metric)
2341 {
2342 case AbsoluteErrorMetric:
2343 {
2344 status=GetAESimilarity(image,reconstruct_image,channel_similarity,
2345 exception);
2346 break;
2347 }
2348 case DotProductCorrelationErrorMetric:
2349 {
2350 status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
2351 exception);
2352 break;
2353 }
2354 case FuzzErrorMetric:
2355 {
2356 status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
2357 exception);
2358 break;
2359 }
2360 case MeanAbsoluteErrorMetric:
2361 {
2362 status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
2363 exception);
2364 break;
2365 }
2366 case MeanErrorPerPixelErrorMetric:
2367 {
2368 status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2369 exception);
2370 break;
2371 }
2372 case MeanSquaredErrorMetric:
2373 {
2374 status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
2375 exception);
2376 break;
2377 }
2378 case NormalizedCrossCorrelationErrorMetric:
2379 {
2380 status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
2381 exception);
2382 break;
2383 }
2384 case PeakAbsoluteErrorMetric:
2385 {
2386 status=GetPASimilarity(image,reconstruct_image,channel_similarity,
2387 exception);
2388 break;
2389 }
2390 case PeakSignalToNoiseRatioErrorMetric:
2391 {
2392 status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
2393 exception);
2394 break;
2395 }
2396 case PerceptualHashErrorMetric:
2397 {
2398 status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
2399 exception);
2400 break;
2401 }
2402 case PhaseCorrelationErrorMetric:
2403 {
2404 status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
2405 exception);
2406 break;
2407 }
2408 case PixelDifferenceCountErrorMetric:
2409 {
2410 status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
2411 exception);
2412 break;
2413 }
2414 case RootMeanSquaredErrorMetric:
2415 case UndefinedErrorMetric:
2416 default:
2417 {
2418 status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
2419 exception);
2420 break;
2421 }
2422 case StructuralDissimilarityErrorMetric:
2423 {
2424 status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
2425 exception);
2426 break;
2427 }
2428 case StructuralSimilarityErrorMetric:
2429 {
2430 status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
2431 exception);
2432 break;
2433 }
2434 }
2435 *distortion=channel_similarity[CompositePixelChannel];
2436 switch (metric)
2437 {
2438 case DotProductCorrelationErrorMetric:
2439 case NormalizedCrossCorrelationErrorMetric:
2440 case PhaseCorrelationErrorMetric:
2441 case StructuralSimilarityErrorMetric:
2442 {
2443 *distortion=(1.0-(*distortion))/2.0;
2444 break;
2445 }
2446 default: break;
2447 }
2448 channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
2449 if (fabs(*distortion) < MagickEpsilon)
2450 *distortion=0.0;
2451 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
2452 *distortion);
2453 return(status);
2454}
2455␌
2456/*
2457%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2458% %
2459% %
2460% %
2461% G e t I m a g e D i s t o r t i o n s %
2462% %
2463% %
2464% %
2465%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2466%
2467% GetImageDistortions() compares the pixel channels of an image to a
2468% reconstructed image and returns the specified metric for each channel.
2469%
2470% The format of the GetImageDistortions method is:
2471%
2472% double *GetImageDistortions(const Image *image,
2473% const Image *reconstruct_image,const MetricType metric,
2474% ExceptionInfo *exception)
2475%
2476% A description of each parameter follows:
2477%
2478% o image: the image.
2479%
2480% o reconstruct_image: the reconstruction image.
2481%
2482% o metric: the metric.
2483%
2484% o exception: return any errors or warnings in this structure.
2485%
2486*/
2487MagickExport double *GetImageDistortions(Image *image,
2488 const Image *reconstruct_image,const MetricType metric,
2489 ExceptionInfo *exception)
2490{
2491 double
2492 *distortion,
2493 *channel_similarity;
2494
2495 MagickBooleanType
2496 status = MagickTrue;
2497
2498 size_t
2499 length;
2500
2501 ssize_t
2502 i;
2503
2504 assert(image != (Image *) NULL);
2505 assert(image->signature == MagickCoreSignature);
2506 assert(reconstruct_image != (const Image *) NULL);
2507 assert(reconstruct_image->signature == MagickCoreSignature);
2508 if (IsEventLogging() != MagickFalse)
2509 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2510 /*
2511 Get image distortion.
2512 */
2513 length=MaxPixelChannels+1UL;
2514 channel_similarity=(double *) AcquireQuantumMemory(length,
2515 sizeof(*channel_similarity));
2516 if (channel_similarity == (double *) NULL)
2517 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2518 (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
2519 switch (metric)
2520 {
2521 case AbsoluteErrorMetric:
2522 {
2523 status=GetAESimilarity(image,reconstruct_image,channel_similarity,
2524 exception);
2525 break;
2526 }
2527 case DotProductCorrelationErrorMetric:
2528 {
2529 status=GetDPCSimilarity(image,reconstruct_image,channel_similarity,
2530 exception);
2531 break;
2532 }
2533 case FuzzErrorMetric:
2534 {
2535 status=GetFUZZSimilarity(image,reconstruct_image,channel_similarity,
2536 exception);
2537 break;
2538 }
2539 case MeanAbsoluteErrorMetric:
2540 {
2541 status=GetMAESimilarity(image,reconstruct_image,channel_similarity,
2542 exception);
2543 break;
2544 }
2545 case MeanErrorPerPixelErrorMetric:
2546 {
2547 status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2548 exception);
2549 break;
2550 }
2551 case MeanSquaredErrorMetric:
2552 {
2553 status=GetMSESimilarity(image,reconstruct_image,channel_similarity,
2554 exception);
2555 break;
2556 }
2557 case NormalizedCrossCorrelationErrorMetric:
2558 {
2559 status=GetNCCSimilarity(image,reconstruct_image,channel_similarity,
2560 exception);
2561 break;
2562 }
2563 case PeakAbsoluteErrorMetric:
2564 {
2565 status=GetPASimilarity(image,reconstruct_image,channel_similarity,
2566 exception);
2567 break;
2568 }
2569 case PeakSignalToNoiseRatioErrorMetric:
2570 {
2571 status=GetPSNRSimilarity(image,reconstruct_image,channel_similarity,
2572 exception);
2573 break;
2574 }
2575 case PerceptualHashErrorMetric:
2576 {
2577 status=GetPHASHSimilarity(image,reconstruct_image,channel_similarity,
2578 exception);
2579 break;
2580 }
2581 case PhaseCorrelationErrorMetric:
2582 {
2583 status=GetPHASESimilarity(image,reconstruct_image,channel_similarity,
2584 exception);
2585 break;
2586 }
2587 case PixelDifferenceCountErrorMetric:
2588 {
2589 status=GetPDCSimilarity(image,reconstruct_image,channel_similarity,
2590 exception);
2591 break;
2592 }
2593 case RootMeanSquaredErrorMetric:
2594 case UndefinedErrorMetric:
2595 default:
2596 {
2597 status=GetRMSESimilarity(image,reconstruct_image,channel_similarity,
2598 exception);
2599 break;
2600 }
2601 case StructuralDissimilarityErrorMetric:
2602 {
2603 status=GetDSSIMSimilarity(image,reconstruct_image,channel_similarity,
2604 exception);
2605 break;
2606 }
2607 case StructuralSimilarityErrorMetric:
2608 {
2609 status=GetSSIMSimularity(image,reconstruct_image,channel_similarity,
2610 exception);
2611 break;
2612 }
2613 }
2614 if (status == MagickFalse)
2615 {
2616 channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
2617 return((double *) NULL);
2618 }
2619 distortion=channel_similarity;
2620 switch (metric)
2621 {
2622 case DotProductCorrelationErrorMetric:
2623 case NormalizedCrossCorrelationErrorMetric:
2624 case PhaseCorrelationErrorMetric:
2625 case StructuralSimilarityErrorMetric:
2626 {
2627 for (i=0; i <= MaxPixelChannels; i++)
2628 distortion[i]=(1.0-distortion[i])/2.0;
2629 break;
2630 }
2631 default: break;
2632 }
2633 for (i=0; i <= MaxPixelChannels; i++)
2634 if (fabs(distortion[i]) < MagickEpsilon)
2635 distortion[i]=0.0;
2636 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
2637 distortion[CompositePixelChannel]);
2638 return(distortion);
2639}
2640␌
2641/*
2642%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2643% %
2644% %
2645% %
2646% I s I m a g e s E q u a l %
2647% %
2648% %
2649% %
2650%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2651%
2652% IsImagesEqual() compare the pixels of two images and returns immediately
2653% if any pixel is not identical.
2654%
2655% The format of the IsImagesEqual method is:
2656%
2657% MagickBooleanType IsImagesEqual(const Image *image,
2658% const Image *reconstruct_image,ExceptionInfo *exception)
2659%
2660% A description of each parameter follows.
2661%
2662% o image: the image.
2663%
2664% o reconstruct_image: the reconstruction image.
2665%
2666% o exception: return any errors or warnings in this structure.
2667%
2668*/
2669MagickExport MagickBooleanType IsImagesEqual(const Image *image,
2670 const Image *reconstruct_image,ExceptionInfo *exception)
2671{
2672 CacheView
2673 *image_view,
2674 *reconstruct_view;
2675
2676 size_t
2677 columns,
2678 rows;
2679
2680 ssize_t
2681 y;
2682
2683 assert(image != (Image *) NULL);
2684 assert(image->signature == MagickCoreSignature);
2685 assert(reconstruct_image != (const Image *) NULL);
2686 assert(reconstruct_image->signature == MagickCoreSignature);
2687 SetImageCompareBounds(image,reconstruct_image,&columns,&rows);
2688 image_view=AcquireVirtualCacheView(image,exception);
2689 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2690 for (y=0; y < (ssize_t) rows; y++)
2691 {
2692 const Quantum
2693 *magick_restrict p,
2694 *magick_restrict q;
2695
2696 ssize_t
2697 x;
2698
2699 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2700 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2701 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
2702 break;
2703 for (x=0; x < (ssize_t) columns; x++)
2704 {
2705 ssize_t
2706 i;
2707
2708 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2709 {
2710 double
2711 distance;
2712
2713 PixelChannel channel = GetPixelChannelChannel(image,i);
2714 PixelTrait traits = GetPixelChannelTraits(image,channel);
2715 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2716 channel);
2717 if (((traits & UpdatePixelTrait) == 0) ||
2718 ((reconstruct_traits & UpdatePixelTrait) == 0))
2719 continue;
2720 distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
2721 channel,q));
2722 if (distance >= MagickEpsilon)
2723 break;
2724 }
2725 if (i < (ssize_t) GetPixelChannels(image))
2726 break;
2727 p+=(ptrdiff_t) GetPixelChannels(image);
2728 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2729 }
2730 if (x < (ssize_t) columns)
2731 break;
2732 }
2733 reconstruct_view=DestroyCacheView(reconstruct_view);
2734 image_view=DestroyCacheView(image_view);
2735 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
2736}
2737␌
2738/*
2739%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2740% %
2741% %
2742% %
2743% S e t I m a g e C o l o r M e t r i c %
2744% %
2745% %
2746% %
2747%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2748%
2749% SetImageColorMetric() measures the difference between colors at each pixel
2750% location of two images. A value other than 0 means the colors match
2751% exactly. Otherwise an error measure is computed by summing over all
2752% pixels in an image the distance squared in RGB space between each image
2753% pixel and its corresponding pixel in the reconstruction image. The error
2754% measure is assigned to these image members:
2755%
2756% o mean_error_per_pixel: The mean error for any single pixel in
2757% the image.
2758%
2759% o normalized_mean_error: The normalized mean quantization error for
2760% any single pixel in the image. This distance measure is normalized to
2761% a range between 0 and 1. It is independent of the range of red, green,
2762% and blue values in the image.
2763%
2764% o normalized_maximum_error: The normalized maximum quantization
2765% error for any single pixel in the image. This distance measure is
2766% normalized to a range between 0 and 1. It is independent of the range
2767% of red, green, and blue values in your image.
2768%
2769% A small normalized mean square error, accessed as
2770% image->normalized_mean_error, suggests the images are very similar in
2771% spatial layout and color.
2772%
2773% The format of the SetImageColorMetric method is:
2774%
2775% MagickBooleanType SetImageColorMetric(Image *image,
2776% const Image *reconstruct_image,ExceptionInfo *exception)
2777%
2778% A description of each parameter follows.
2779%
2780% o image: the image.
2781%
2782% o reconstruct_image: the reconstruction image.
2783%
2784% o exception: return any errors or warnings in this structure.
2785%
2786*/
2787MagickExport MagickBooleanType SetImageColorMetric(Image *image,
2788 const Image *reconstruct_image,ExceptionInfo *exception)
2789{
2790 double
2791 channel_similarity[MaxPixelChannels+1] = { 0.0 };
2792
2793 MagickBooleanType
2794 status;
2795
2796 status=GetMEPPSimilarity(image,reconstruct_image,channel_similarity,
2797 exception);
2798 if (status == MagickFalse)
2799 return(MagickFalse);
2800 status=fabs(image->error.mean_error_per_pixel) < MagickEpsilon ?
2801 MagickTrue : MagickFalse;
2802 return(status);
2803}
2804␌
2805/*
2806%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2807% %
2808% %
2809% %
2810% S i m i l a r i t y I m a g e %
2811% %
2812% %
2813% %
2814%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2815%
2816% SimilarityImage() compares the reconstruction of the image and returns the
2817% best match offset. In addition, it returns a similarity image such that an
2818% exact match location is completely white and if none of the pixels match,
2819% black, otherwise some gray level in-between.
2820%
2821% Contributed by Fred Weinhaus.
2822%
2823% The format of the SimilarityImageImage method is:
2824%
2825% Image *SimilarityImage(const Image *image,const Image *reconstruct,
2826% const MetricType metric,const double similarity_threshold,
2827% RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2828%
2829% A description of each parameter follows:
2830%
2831% o image: the image.
2832%
2833% o reconstruct: find an area of the image that closely resembles this image.
2834%
2835% o metric: the metric.
2836%
2837% o similarity_threshold: minimum similarity for (sub)image match.
2838%
2839% o offset: the best match offset of the reconstruction image within the
2840% image.
2841%
2842% o similarity: the computed similarity between the images.
2843%
2844% o exception: return any errors or warnings in this structure.
2845%
2846*/
2847
2848#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2849static Image *SIMCrossCorrelationImage(const Image *alpha_image,
2850 const Image *beta_image,ExceptionInfo *exception)
2851{
2852 Image
2853 *alpha_fft = (Image *) NULL,
2854 *beta_fft = (Image *) NULL,
2855 *complex_conjugate = (Image *) NULL,
2856 *complex_multiplication = (Image *) NULL,
2857 *cross_correlation = (Image *) NULL,
2858 *temp_image = (Image *) NULL;
2859
2860 /*
2861 Take the FFT of beta (reconstruction) image.
2862 */
2863 temp_image=CloneImage(beta_image,0,0,MagickTrue,exception);
2864 if (temp_image == (Image *) NULL)
2865 return((Image *) NULL);
2866 (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2867 beta_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2868 temp_image=DestroyImageList(temp_image);
2869 if (beta_fft == (Image *) NULL)
2870 return((Image *) NULL);
2871 /*
2872 Take the complex conjugate of beta_fft.
2873 */
2874 complex_conjugate=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
2875 beta_fft=DestroyImageList(beta_fft);
2876 if (complex_conjugate == (Image *) NULL)
2877 return((Image *) NULL);
2878 /*
2879 Take the FFT of the alpha (test) image.
2880 */
2881 temp_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2882 if (temp_image == (Image *) NULL)
2883 {
2884 complex_conjugate=DestroyImageList(complex_conjugate);
2885 return((Image *) NULL);
2886 }
2887 (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2888 alpha_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2889 temp_image=DestroyImageList(temp_image);
2890 if (alpha_fft == (Image *) NULL)
2891 {
2892 complex_conjugate=DestroyImageList(complex_conjugate);
2893 return((Image *) NULL);
2894 }
2895 /*
2896 Do complex multiplication.
2897 */
2898 DisableCompositeClampUnlessSpecified(complex_conjugate);
2899 DisableCompositeClampUnlessSpecified(complex_conjugate->next);
2900 AppendImageToList(&complex_conjugate,alpha_fft);
2901 complex_multiplication=ComplexImages(complex_conjugate,
2902 MultiplyComplexOperator,exception);
2903 complex_conjugate=DestroyImageList(complex_conjugate);
2904 if (complex_multiplication == (Image *) NULL)
2905 return((Image *) NULL);
2906 /*
2907 Do the IFT and return the cross-correlation result.
2908 */
2909 cross_correlation=InverseFourierTransformImage(complex_multiplication,
2910 complex_multiplication->next,MagickFalse,exception);
2911 complex_multiplication=DestroyImageList(complex_multiplication);
2912 return(cross_correlation);
2913}
2914
2915static Image *SIMDerivativeImage(const Image *image,const char *kernel,
2916 ExceptionInfo *exception)
2917{
2918 Image
2919 *derivative_image;
2920
2921 KernelInfo
2922 *kernel_info;
2923
2924 kernel_info=AcquireKernelInfo(kernel,exception);
2925 if (kernel_info == (KernelInfo *) NULL)
2926 return((Image *) NULL);
2927 derivative_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,
2928 exception);
2929 kernel_info=DestroyKernelInfo(kernel_info);
2930 return(derivative_image);
2931}
2932
2933static Image *SIMDivideImage(const Image *numerator_image,
2934 const Image *denominator_image,ExceptionInfo *exception)
2935{
2936 CacheView
2937 *denominator_view,
2938 *numerator_view;
2939
2940 Image
2941 *divide_image;
2942
2943 MagickBooleanType
2944 status = MagickTrue;
2945
2946 ssize_t
2947 y;
2948
2949 /*
2950 Divide one image into another.
2951 */
2952 divide_image=CloneImage(numerator_image,0,0,MagickTrue,exception);
2953 if (divide_image == (Image *) NULL)
2954 return(divide_image);
2955 numerator_view=AcquireAuthenticCacheView(divide_image,exception);
2956 denominator_view=AcquireVirtualCacheView(denominator_image,exception);
2957#if defined(MAGICKCORE_OPENMP_SUPPORT)
2958 #pragma omp parallel for schedule(static) shared(status) \
2959 magick_number_threads(denominator_image,divide_image,divide_image->rows,1)
2960#endif
2961 for (y=0; y < (ssize_t) divide_image->rows; y++)
2962 {
2963 const Quantum
2964 *magick_restrict p;
2965
2966 Quantum
2967 *magick_restrict q;
2968
2969 ssize_t
2970 x;
2971
2972 if (status == MagickFalse)
2973 continue;
2974 p=GetCacheViewVirtualPixels(denominator_view,0,y,
2975 denominator_image->columns,1,exception);
2976 q=GetCacheViewAuthenticPixels(numerator_view,0,y,divide_image->columns,1,
2977 exception);
2978 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2979 {
2980 status=MagickFalse;
2981 continue;
2982 }
2983 for (x=0; x < (ssize_t) divide_image->columns; x++)
2984 {
2985 ssize_t
2986 i;
2987
2988 for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
2989 {
2990 PixelChannel channel = GetPixelChannelChannel(divide_image,i);
2991 PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
2992 PixelTrait denominator_traits = GetPixelChannelTraits(denominator_image,
2993 channel);
2994 if (((traits & UpdatePixelTrait) == 0) ||
2995 ((denominator_traits & UpdatePixelTrait) == 0))
2996 continue;
2997 q[i]=(Quantum) ((double) q[i]*MagickSafeReciprocal(QuantumScale*
2998 (double) GetPixelChannel(denominator_image,channel,p)));
2999 }
3000 p+=(ptrdiff_t) GetPixelChannels(denominator_image);
3001 q+=(ptrdiff_t) GetPixelChannels(divide_image);
3002 }
3003 if (SyncCacheViewAuthenticPixels(numerator_view,exception) == MagickFalse)
3004 status=MagickFalse;
3005 }
3006 denominator_view=DestroyCacheView(denominator_view);
3007 numerator_view=DestroyCacheView(numerator_view);
3008 if (status == MagickFalse)
3009 divide_image=DestroyImage(divide_image);
3010 return(divide_image);
3011}
3012
3013static Image *SIMDivideByMagnitude(Image *image,Image *magnitude_image,
3014 const Image *source_image,ExceptionInfo *exception)
3015{
3016 Image
3017 *divide_image,
3018 *result_image;
3019
3020 RectangleInfo
3021 geometry;
3022
3023 divide_image=SIMDivideImage(image,magnitude_image,exception);
3024 if (divide_image == (Image *) NULL)
3025 return((Image *) NULL);
3026 GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
3027 &divide_image->background_color);
3028 SetGeometry(source_image,&geometry);
3029 geometry.width=MagickMax(source_image->columns,divide_image->columns);
3030 geometry.height=MagickMax(source_image->rows,divide_image->rows);
3031 result_image=ExtentImage(divide_image,&geometry,exception);
3032 divide_image=DestroyImage(divide_image);
3033 return(result_image);
3034}
3035
3036static MagickBooleanType SIMFilterImageNaNs(Image *image,
3037 ExceptionInfo *exception)
3038{
3039 CacheView
3040 *image_view;
3041
3042 MagickBooleanType
3043 status = MagickTrue;
3044
3045 ssize_t
3046 y;
3047
3048 /*
3049 Square each pixel in the image.
3050 */
3051 image_view=AcquireAuthenticCacheView(image,exception);
3052#if defined(MAGICKCORE_OPENMP_SUPPORT)
3053 #pragma omp parallel for schedule(static) shared(status) \
3054 magick_number_threads(image,image,image->rows,1)
3055#endif
3056 for (y=0; y < (ssize_t) image->rows; y++)
3057 {
3058 Quantum
3059 *magick_restrict q;
3060
3061 ssize_t
3062 x;
3063
3064 if (status == MagickFalse)
3065 continue;
3066 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3067 if (q == (Quantum *) NULL)
3068 {
3069 status=MagickFalse;
3070 continue;
3071 }
3072 for (x=0; x < (ssize_t) image->columns; x++)
3073 {
3074 ssize_t
3075 i;
3076
3077 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3078 {
3079 PixelChannel channel = GetPixelChannelChannel(image,i);
3080 PixelTrait traits = GetPixelChannelTraits(image,channel);
3081 if ((traits & UpdatePixelTrait) == 0)
3082 continue;
3083 if (IsNaN((double) q[i]) != 0)
3084 q[i]=(Quantum) 0;
3085 }
3086 q+=(ptrdiff_t) GetPixelChannels(image);
3087 }
3088 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3089 status=MagickFalse;
3090 }
3091 image_view=DestroyCacheView(image_view);
3092 return(status);
3093}
3094
3095static Image *SIMSquareImage(const Image *image,ExceptionInfo *exception)
3096{
3097 CacheView
3098 *image_view;
3099
3100 Image
3101 *square_image;
3102
3103 MagickBooleanType
3104 status = MagickTrue;
3105
3106 ssize_t
3107 y;
3108
3109 /*
3110 Square each pixel in the image.
3111 */
3112 square_image=CloneImage(image,0,0,MagickTrue,exception);
3113 if (square_image == (Image *) NULL)
3114 return(square_image);
3115 image_view=AcquireAuthenticCacheView(square_image,exception);
3116#if defined(MAGICKCORE_OPENMP_SUPPORT)
3117 #pragma omp parallel for schedule(static) shared(status) \
3118 magick_number_threads(square_image,square_image,square_image->rows,1)
3119#endif
3120 for (y=0; y < (ssize_t) square_image->rows; y++)
3121 {
3122 Quantum
3123 *magick_restrict q;
3124
3125 ssize_t
3126 x;
3127
3128 if (status == MagickFalse)
3129 continue;
3130 q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
3131 exception);
3132 if (q == (Quantum *) NULL)
3133 {
3134 status=MagickFalse;
3135 continue;
3136 }
3137 for (x=0; x < (ssize_t) square_image->columns; x++)
3138 {
3139 ssize_t
3140 i;
3141
3142 for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
3143 {
3144 PixelChannel channel = GetPixelChannelChannel(square_image,i);
3145 PixelTrait traits = GetPixelChannelTraits(square_image,channel);
3146 if ((traits & UpdatePixelTrait) == 0)
3147 continue;
3148 q[i]=(Quantum) (QuantumScale*q[i]*q[i]);
3149 }
3150 q+=(ptrdiff_t) GetPixelChannels(square_image);
3151 }
3152 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3153 status=MagickFalse;
3154 }
3155 image_view=DestroyCacheView(image_view);
3156 if (status == MagickFalse)
3157 square_image=DestroyImage(square_image);
3158 return(square_image);
3159}
3160
3161static Image *SIMMagnitudeImage(Image *alpha_image,Image *beta_image,
3162 ExceptionInfo *exception)
3163{
3164 Image
3165 *magnitude_image,
3166 *xsq_image,
3167 *ysq_image;
3168
3169 MagickBooleanType
3170 status = MagickTrue;
3171
3172 (void) SetImageArtifact(alpha_image,"compose:clamp","False");
3173 xsq_image=SIMSquareImage(alpha_image,exception);
3174 if (xsq_image == (Image *) NULL)
3175 return((Image *) NULL);
3176 (void) SetImageArtifact(beta_image,"compose:clamp","False");
3177 ysq_image=SIMSquareImage(beta_image,exception);
3178 if (ysq_image == (Image *) NULL)
3179 {
3180 xsq_image=DestroyImage(xsq_image);
3181 return((Image *) NULL);
3182 }
3183 status=CompositeImage(xsq_image,ysq_image,PlusCompositeOp,MagickTrue,0,0,
3184 exception);
3185 magnitude_image=xsq_image;
3186 ysq_image=DestroyImage(ysq_image);
3187 if (status == MagickFalse)
3188 {
3189 magnitude_image=DestroyImage(magnitude_image);
3190 return((Image *) NULL);
3191 }
3192 status=EvaluateImage(magnitude_image,PowEvaluateOperator,0.5,exception);
3193 if (status == MagickFalse)
3194 {
3195 magnitude_image=DestroyImage(magnitude_image);
3196 return (Image *) NULL;
3197 }
3198 return(magnitude_image);
3199}
3200
3201static MagickBooleanType SIMMaximaImage(const Image *image,double *maxima,
3202 RectangleInfo *offset,ExceptionInfo *exception)
3203{
3204 typedef struct
3205 {
3206 double
3207 maxima;
3208
3209 ssize_t
3210 x,
3211 y;
3212 } MaximaInfo;
3213
3214 CacheView
3215 *image_view;
3216
3217 const Quantum
3218 *magick_restrict q;
3219
3220 MagickBooleanType
3221 status = MagickTrue;
3222
3223 MaximaInfo
3224 maxima_info = { -MagickMaximumValue, 0, 0 };
3225
3226 ssize_t
3227 y;
3228
3229 /*
3230 Identify the maxima value in the image and its location.
3231 */
3232 image_view=AcquireVirtualCacheView(image,exception);
3233 q=GetCacheViewVirtualPixels(image_view,maxima_info.x,maxima_info.y,1,1,
3234 exception);
3235 if (q != (const Quantum *) NULL)
3236 maxima_info.maxima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
3237#if defined(MAGICKCORE_OPENMP_SUPPORT)
3238 #pragma omp parallel for schedule(static) shared(maxima_info,status) \
3239 magick_number_threads(image,image,image->rows,1)
3240#endif
3241 for (y=0; y < (ssize_t) image->rows; y++)
3242 {
3243 const Quantum
3244 *magick_restrict p;
3245
3246 MaximaInfo
3247 channel_maxima = { -MagickMaximumValue, 0, 0 };
3248
3249 ssize_t
3250 x;
3251
3252 if (status == MagickFalse)
3253 continue;
3254 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3255 if (p == (const Quantum *) NULL)
3256 {
3257 status=MagickFalse;
3258 continue;
3259 }
3260 channel_maxima=maxima_info;
3261 for (x=0; x < (ssize_t) image->columns; x++)
3262 {
3263 ssize_t
3264 i;
3265
3266 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3267 {
3268 double
3269 pixel;
3270
3271 PixelChannel channel = GetPixelChannelChannel(image,i);
3272 PixelTrait traits = GetPixelChannelTraits(image,channel);
3273 if ((traits & UpdatePixelTrait) == 0)
3274 continue;
3275 pixel=(double) p[i];
3276 if (IsNaN(pixel) != 0)
3277 pixel=0.0;
3278 if (pixel > channel_maxima.maxima)
3279 {
3280 channel_maxima.maxima=(double) p[i];
3281 channel_maxima.x=x;
3282 channel_maxima.y=y;
3283 }
3284 }
3285 p+=(ptrdiff_t) GetPixelChannels(image);
3286 }
3287#if defined(MAGICKCORE_OPENMP_SUPPORT)
3288 #pragma omp critical (MagickCore_SIMMaximaImage)
3289#endif
3290 if (channel_maxima.maxima > maxima_info.maxima)
3291 maxima_info=channel_maxima;
3292 }
3293 image_view=DestroyCacheView(image_view);
3294 *maxima=maxima_info.maxima;
3295 offset->x=maxima_info.x;
3296 offset->y=maxima_info.y;
3297 return(status);
3298}
3299
3300static MagickBooleanType SIMMinimaImage(const Image *image,double *minima,
3301 RectangleInfo *offset,ExceptionInfo *exception)
3302{
3303 typedef struct
3304 {
3305 double
3306 minima;
3307
3308 ssize_t
3309 x,
3310 y;
3311 } MinimaInfo;
3312
3313 CacheView
3314 *image_view;
3315
3316 const Quantum
3317 *magick_restrict q;
3318
3319 MagickBooleanType
3320 status = MagickTrue;
3321
3322 MinimaInfo
3323 minima_info = { MagickMaximumValue, 0, 0 };
3324
3325 ssize_t
3326 y;
3327
3328 /*
3329 Identify the minima value in the image and its location.
3330 */
3331 image_view=AcquireVirtualCacheView(image,exception);
3332 q=GetCacheViewVirtualPixels(image_view,minima_info.x,minima_info.y,1,1,
3333 exception);
3334 if (q != (const Quantum *) NULL)
3335 minima_info.minima=IsNaN((double) q[0]) != 0 ? 0.0 : (double) q[0];
3336#if defined(MAGICKCORE_OPENMP_SUPPORT)
3337 #pragma omp parallel for schedule(static) shared(minima_info,status) \
3338 magick_number_threads(image,image,image->rows,1)
3339#endif
3340 for (y=0; y < (ssize_t) image->rows; y++)
3341 {
3342 const Quantum
3343 *magick_restrict p;
3344
3345 MinimaInfo
3346 channel_minima = { MagickMaximumValue, 0, 0 };
3347
3348 ssize_t
3349 x;
3350
3351 if (status == MagickFalse)
3352 continue;
3353 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3354 if (p == (const Quantum *) NULL)
3355 {
3356 status=MagickFalse;
3357 continue;
3358 }
3359 channel_minima=minima_info;
3360 for (x=0; x < (ssize_t) image->columns; x++)
3361 {
3362 ssize_t
3363 i;
3364
3365 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3366 {
3367 double
3368 pixel;
3369
3370 PixelChannel channel = GetPixelChannelChannel(image,i);
3371 PixelTrait traits = GetPixelChannelTraits(image,channel);
3372 if ((traits & UpdatePixelTrait) == 0)
3373 continue;
3374 pixel=(double) p[i];
3375 if (IsNaN(pixel) != 0)
3376 pixel=0.0;
3377 if (pixel < channel_minima.minima)
3378 {
3379 channel_minima.minima=pixel;
3380 channel_minima.x=x;
3381 channel_minima.y=y;
3382 }
3383 }
3384 p+=(ptrdiff_t) GetPixelChannels(image);
3385 }
3386#if defined(MAGICKCORE_OPENMP_SUPPORT)
3387 #pragma omp critical (MagickCore_SIMMinimaImage)
3388#endif
3389 if (channel_minima.minima < minima_info.minima)
3390 minima_info=channel_minima;
3391 }
3392 image_view=DestroyCacheView(image_view);
3393 *minima=minima_info.minima;
3394 offset->x=minima_info.x;
3395 offset->y=minima_info.y;
3396 return(status);
3397}
3398
3399static MagickBooleanType SIMMultiplyImage(Image *image,const double factor,
3400 const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
3401{
3402 CacheView
3403 *image_view;
3404
3405 MagickBooleanType
3406 status = MagickTrue;
3407
3408 ssize_t
3409 y;
3410
3411 /*
3412 Multiply each pixel by a factor.
3413 */
3414 image_view=AcquireAuthenticCacheView(image,exception);
3415#if defined(MAGICKCORE_OPENMP_SUPPORT)
3416 #pragma omp parallel for schedule(static) shared(status) \
3417 magick_number_threads(image,image,image->rows,1)
3418#endif
3419 for (y=0; y < (ssize_t) image->rows; y++)
3420 {
3421 Quantum
3422 *magick_restrict q;
3423
3424 ssize_t
3425 x;
3426
3427 if (status == MagickFalse)
3428 continue;
3429 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3430 if (q == (Quantum *) NULL)
3431 {
3432 status=MagickFalse;
3433 continue;
3434 }
3435 for (x=0; x < (ssize_t) image->columns; x++)
3436 {
3437 ssize_t
3438 i;
3439
3440 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3441 {
3442 PixelChannel channel = GetPixelChannelChannel(image,i);
3443 PixelTrait traits = GetPixelChannelTraits(image,channel);
3444 if ((traits & UpdatePixelTrait) == 0)
3445 continue;
3446 if (channel_statistics != (const ChannelStatistics *) NULL)
3447 q[i]=(Quantum) (factor*q[i]*QuantumScale*
3448 channel_statistics[channel].standard_deviation);
3449 else
3450 q[i]=(Quantum) (factor*q[i]);
3451 }
3452 q+=(ptrdiff_t) GetPixelChannels(image);
3453 }
3454 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3455 status=MagickFalse;
3456 }
3457 image_view=DestroyCacheView(image_view);
3458 return(status);
3459}
3460
3461static Image *SIMPhaseCorrelationImage(const Image *alpha_image,
3462 const Image *beta_image,const Image *magnitude_image,ExceptionInfo *exception)
3463{
3464 Image
3465 *alpha_fft = (Image *) NULL,
3466 *beta_fft = (Image *) NULL,
3467 *complex_multiplication = (Image *) NULL,
3468 *cross_correlation = (Image *) NULL;
3469
3470 /*
3471 Take the FFT of the beta (reconstruction) image.
3472 */
3473 beta_fft=CloneImage(beta_image,0,0,MagickTrue,exception);
3474 if (beta_fft == NULL)
3475 return((Image *) NULL);
3476 (void) SetImageArtifact(beta_fft,"fourier:normalize","inverse");
3477 beta_fft=ForwardFourierTransformImage(beta_fft,MagickFalse,exception);
3478 if (beta_fft == NULL)
3479 return((Image *) NULL);
3480 /*
3481 Take the FFT of the alpha (test) image.
3482 */
3483 alpha_fft=CloneImage(alpha_image,0,0,MagickTrue,exception);
3484 if (alpha_fft == (Image *) NULL)
3485 {
3486 beta_fft=DestroyImageList(beta_fft);
3487 return((Image *) NULL);
3488 }
3489 (void) SetImageArtifact(alpha_fft,"fourier:normalize","inverse");
3490 alpha_fft=ForwardFourierTransformImage(alpha_fft,MagickFalse,exception);
3491 if (alpha_fft == (Image *) NULL)
3492 {
3493 beta_fft=DestroyImageList(beta_fft);
3494 return((Image *) NULL);
3495 }
3496 /*
3497 Take the complex conjugate of the beta FFT.
3498 */
3499 beta_fft=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
3500 if (beta_fft == (Image *) NULL)
3501 {
3502 alpha_fft=DestroyImageList(alpha_fft);
3503 return((Image *) NULL);
3504 }
3505 /*
3506 Do complex multiplication.
3507 */
3508 AppendImageToList(&beta_fft,alpha_fft);
3509 DisableCompositeClampUnlessSpecified(beta_fft);
3510 DisableCompositeClampUnlessSpecified(beta_fft->next);
3511 complex_multiplication=ComplexImages(beta_fft,MultiplyComplexOperator,
3512 exception);
3513 beta_fft=DestroyImageList(beta_fft);
3514 if (complex_multiplication == (Image *) NULL)
3515 return((Image *) NULL);
3516 /*
3517 Divide the results.
3518 */
3519 CompositeLayers(complex_multiplication,DivideSrcCompositeOp,(Image *)
3520 magnitude_image,0,0,exception);
3521 /*
3522 Do the IFT and return the cross-correlation result.
3523 */
3524 (void) SetImageArtifact(complex_multiplication,"fourier:normalize","inverse");
3525 cross_correlation=InverseFourierTransformImage(complex_multiplication,
3526 complex_multiplication->next,MagickFalse,exception);
3527 complex_multiplication=DestroyImageList(complex_multiplication);
3528 return(cross_correlation);
3529}
3530
3531static MagickBooleanType SIMSetImageMean(Image *image,
3532 const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
3533{
3534 CacheView
3535 *image_view;
3536
3537 MagickBooleanType
3538 status = MagickTrue;
3539
3540 ssize_t
3541 y;
3542
3543 /*
3544 Set image mean.
3545 */
3546 image_view=AcquireAuthenticCacheView(image,exception);
3547#if defined(MAGICKCORE_OPENMP_SUPPORT)
3548 #pragma omp parallel for schedule(static) shared(status) \
3549 magick_number_threads(image,image,image->rows,1)
3550#endif
3551 for (y=0; y < (ssize_t) image->rows; y++)
3552 {
3553 Quantum
3554 *magick_restrict q;
3555
3556 ssize_t
3557 x;
3558
3559 if (status == MagickFalse)
3560 continue;
3561 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3562 if (q == (Quantum *) NULL)
3563 {
3564 status=MagickFalse;
3565 continue;
3566 }
3567 for (x=0; x < (ssize_t) image->columns; x++)
3568 {
3569 ssize_t
3570 i;
3571
3572 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3573 {
3574 PixelChannel channel = GetPixelChannelChannel(image,i);
3575 PixelTrait traits = GetPixelChannelTraits(image,channel);
3576 if ((traits & UpdatePixelTrait) == 0)
3577 continue;
3578 q[i]=(Quantum) channel_statistics[channel].mean;
3579 }
3580 q+=(ptrdiff_t) GetPixelChannels(image);
3581 }
3582 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3583 status=MagickFalse;
3584 }
3585 image_view=DestroyCacheView(image_view);
3586 return(status);
3587}
3588
3589static Image *SIMSubtractImageMean(const Image *alpha_image,
3590 const Image *beta_image,const ChannelStatistics *channel_statistics,
3591 ExceptionInfo *exception)
3592{
3593 CacheView
3594 *beta_view,
3595 *image_view;
3596
3597 Image
3598 *subtract_image;
3599
3600 MagickBooleanType
3601 status = MagickTrue;
3602
3603 ssize_t
3604 y;
3605
3606 /*
3607 Subtract the image mean and pad.
3608 */
3609 subtract_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
3610 MagickTrue,exception);
3611 if (subtract_image == (Image *) NULL)
3612 return(subtract_image);
3613 image_view=AcquireAuthenticCacheView(subtract_image,exception);
3614 beta_view=AcquireVirtualCacheView(beta_image,exception);
3615#if defined(MAGICKCORE_OPENMP_SUPPORT)
3616 #pragma omp parallel for schedule(static) shared(status) \
3617 magick_number_threads(beta_image,subtract_image,subtract_image->rows,1)
3618#endif
3619 for (y=0; y < (ssize_t) subtract_image->rows; y++)
3620 {
3621 const Quantum
3622 *magick_restrict p;
3623
3624 Quantum
3625 *magick_restrict q;
3626
3627 ssize_t
3628 x;
3629
3630 if (status == MagickFalse)
3631 continue;
3632 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,exception);
3633 q=GetCacheViewAuthenticPixels(image_view,0,y,subtract_image->columns,1,
3634 exception);
3635 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3636 {
3637 status=MagickFalse;
3638 continue;
3639 }
3640 for (x=0; x < (ssize_t) subtract_image->columns; x++)
3641 {
3642 ssize_t
3643 i;
3644
3645 for (i=0; i < (ssize_t) GetPixelChannels(subtract_image); i++)
3646 {
3647 PixelChannel channel = GetPixelChannelChannel(subtract_image,i);
3648 PixelTrait traits = GetPixelChannelTraits(subtract_image,channel);
3649 PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
3650 if (((traits & UpdatePixelTrait) == 0) ||
3651 ((beta_traits & UpdatePixelTrait) == 0))
3652 continue;
3653 if ((x >= (ssize_t) beta_image->columns) ||
3654 (y >= (ssize_t) beta_image->rows))
3655 q[i]=(Quantum) 0;
3656 else
3657 q[i]=(Quantum) ((double) GetPixelChannel(beta_image,channel,p)-
3658 channel_statistics[channel].mean);
3659 }
3660 p+=(ptrdiff_t) GetPixelChannels(beta_image);
3661 q+=(ptrdiff_t) GetPixelChannels(subtract_image);
3662 }
3663 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3664 status=MagickFalse;
3665 }
3666 beta_view=DestroyCacheView(beta_view);
3667 image_view=DestroyCacheView(image_view);
3668 if (status == MagickFalse)
3669 subtract_image=DestroyImage(subtract_image);
3670 return(subtract_image);
3671}
3672
3673static Image *SIMUnityImage(const Image *alpha_image,const Image *beta_image,
3674 ExceptionInfo *exception)
3675{
3676 CacheView
3677 *image_view;
3678
3679 Image
3680 *unity_image;
3681
3682 MagickBooleanType
3683 status = MagickTrue;
3684
3685 ssize_t
3686 y;
3687
3688 /*
3689 Create a padded unity image.
3690 */
3691 unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
3692 MagickTrue,exception);
3693 if (unity_image == (Image *) NULL)
3694 return(unity_image);
3695 if (SetImageStorageClass(unity_image,DirectClass,exception) == MagickFalse)
3696 return(DestroyImage(unity_image));
3697 image_view=AcquireAuthenticCacheView(unity_image,exception);
3698#if defined(MAGICKCORE_OPENMP_SUPPORT)
3699 #pragma omp parallel for schedule(static) shared(status) \
3700 magick_number_threads(unity_image,unity_image,unity_image->rows,1)
3701#endif
3702 for (y=0; y < (ssize_t) unity_image->rows; y++)
3703 {
3704 Quantum
3705 *magick_restrict q;
3706
3707 ssize_t
3708 x;
3709
3710 if (status == MagickFalse)
3711 continue;
3712 q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
3713 exception);
3714 if (q == (Quantum *) NULL)
3715 {
3716 status=MagickFalse;
3717 continue;
3718 }
3719 for (x=0; x < (ssize_t) unity_image->columns; x++)
3720 {
3721 ssize_t
3722 i;
3723
3724 for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
3725 {
3726 PixelChannel channel = GetPixelChannelChannel(unity_image,i);
3727 PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
3728 if ((traits & UpdatePixelTrait) == 0)
3729 continue;
3730 if ((x >= (ssize_t) beta_image->columns) ||
3731 (y >= (ssize_t) beta_image->rows))
3732 q[i]=(Quantum) 0;
3733 else
3734 q[i]=QuantumRange;
3735 }
3736 q+=(ptrdiff_t) GetPixelChannels(unity_image);
3737 }
3738 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3739 status=MagickFalse;
3740 }
3741 image_view=DestroyCacheView(image_view);
3742 if (status == MagickFalse)
3743 unity_image=DestroyImage(unity_image);
3744 return(unity_image);
3745}
3746
3747static Image *SIMVarianceImage(Image *alpha_image,const Image *beta_image,
3748 ExceptionInfo *exception)
3749{
3750 CacheView
3751 *beta_view,
3752 *image_view;
3753
3754 Image
3755 *variance_image;
3756
3757 MagickBooleanType
3758 status = MagickTrue;
3759
3760 ssize_t
3761 y;
3762
3763 /*
3764 Compute the variance of the two images.
3765 */
3766 variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
3767 if (variance_image == (Image *) NULL)
3768 return(variance_image);
3769 image_view=AcquireAuthenticCacheView(variance_image,exception);
3770 beta_view=AcquireVirtualCacheView(beta_image,exception);
3771#if defined(MAGICKCORE_OPENMP_SUPPORT)
3772 #pragma omp parallel for schedule(static) shared(status) \
3773 magick_number_threads(beta_image,variance_image,variance_image->rows,1)
3774#endif
3775 for (y=0; y < (ssize_t) variance_image->rows; y++)
3776 {
3777 const Quantum
3778 *magick_restrict p;
3779
3780 Quantum
3781 *magick_restrict q;
3782
3783 ssize_t
3784 x;
3785
3786 if (status == MagickFalse)
3787 continue;
3788 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
3789 exception);
3790 q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
3791 exception);
3792 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3793 {
3794 status=MagickFalse;
3795 continue;
3796 }
3797 for (x=0; x < (ssize_t) variance_image->columns; x++)
3798 {
3799 ssize_t
3800 i;
3801
3802 for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
3803 {
3804 double
3805 error;
3806
3807 PixelChannel channel = GetPixelChannelChannel(variance_image,i);
3808 PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
3809 PixelTrait beta_traits = GetPixelChannelTraits(beta_image,channel);
3810 if (((traits & UpdatePixelTrait) == 0) ||
3811 ((beta_traits & UpdatePixelTrait) == 0))
3812 continue;
3813 error=(double) q[i]-(double) GetPixelChannel(beta_image,channel,p);
3814 q[i]=(Quantum) ((double) ClampToQuantum((double) QuantumRange*
3815 (sqrt(fabs(QuantumScale*error))/sqrt((double) QuantumRange))));
3816 }
3817 p+=(ptrdiff_t) GetPixelChannels(beta_image);
3818 q+=(ptrdiff_t) GetPixelChannels(variance_image);
3819 }
3820 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3821 status=MagickFalse;
3822 }
3823 beta_view=DestroyCacheView(beta_view);
3824 image_view=DestroyCacheView(image_view);
3825 if (status == MagickFalse)
3826 variance_image=DestroyImage(variance_image);
3827 return(variance_image);
3828}
3829
3830static Image *DPCSimilarityImage(const Image *image,const Image *reconstruct,
3831 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3832{
3833#define ThrowDPCSimilarityException() \
3834{ \
3835 if (dot_product_image != (Image *) NULL) \
3836 dot_product_image=DestroyImage(dot_product_image); \
3837 if (magnitude_image != (Image *) NULL) \
3838 magnitude_image=DestroyImage(magnitude_image); \
3839 if (reconstruct_image != (Image *) NULL) \
3840 reconstruct_image=DestroyImage(reconstruct_image); \
3841 if (rx_image != (Image *) NULL) \
3842 rx_image=DestroyImage(rx_image); \
3843 if (ry_image != (Image *) NULL) \
3844 ry_image=DestroyImage(ry_image); \
3845 if (target_image != (Image *) NULL) \
3846 target_image=DestroyImage(target_image); \
3847 if (threshold_image != (Image *) NULL) \
3848 threshold_image=DestroyImage(threshold_image); \
3849 if (trx_image != (Image *) NULL) \
3850 trx_image=DestroyImage(trx_image); \
3851 if (try_image != (Image *) NULL) \
3852 try_image=DestroyImage(try_image); \
3853 if (tx_image != (Image *) NULL) \
3854 tx_image=DestroyImage(tx_image); \
3855 if (ty_image != (Image *) NULL) \
3856 ty_image=DestroyImage(ty_image); \
3857 return((Image *) NULL); \
3858}
3859
3860 double
3861 edge_factor = 0.0,
3862 maxima = 0.0,
3863 mean = 0.0,
3864 standard_deviation = 0.0;
3865
3866 Image
3867 *dot_product_image = (Image *) NULL,
3868 *magnitude_image = (Image *) NULL,
3869 *reconstruct_image = (Image *) NULL,
3870 *rx_image = (Image *) NULL,
3871 *ry_image = (Image *) NULL,
3872 *trx_image = (Image *) NULL,
3873 *target_image = (Image *) NULL,
3874 *threshold_image = (Image *) NULL,
3875 *try_image = (Image *) NULL,
3876 *tx_image = (Image *) NULL,
3877 *ty_image = (Image *) NULL;
3878
3879 MagickBooleanType
3880 status = MagickTrue;
3881
3882 RectangleInfo
3883 geometry;
3884
3885 /*
3886 Dot product correlation-based image similarity using FFT local statistics.
3887 */
3888 target_image=CloneImage(image,0,0,MagickTrue,exception);
3889 if (target_image == (Image *) NULL)
3890 return((Image *) NULL);
3891 /*
3892 Compute the cross correlation of the test and reconstruct magnitudes.
3893 */
3894 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
3895 if (reconstruct_image == (Image *) NULL)
3896 ThrowDPCSimilarityException();
3897 /*
3898 Compute X and Y derivatives of reference image.
3899 */
3900 (void) SetImageVirtualPixelMethod(reconstruct_image,EdgeVirtualPixelMethod,
3901 exception);
3902 rx_image=SIMDerivativeImage(reconstruct_image,"Sobel",exception);
3903 if (rx_image == (Image *) NULL)
3904 ThrowDPCSimilarityException();
3905 ry_image=SIMDerivativeImage(reconstruct_image,"Sobel:90",exception);
3906 reconstruct_image=DestroyImage(reconstruct_image);
3907 if (ry_image == (Image *) NULL)
3908 ThrowDPCSimilarityException();
3909 /*
3910 Compute magnitude of derivatives.
3911 */
3912 magnitude_image=SIMMagnitudeImage(rx_image,ry_image,exception);
3913 if (magnitude_image == (Image *) NULL)
3914 ThrowDPCSimilarityException();
3915 /*
3916 Compute an edge normalization correction.
3917 */
3918 threshold_image=CloneImage(magnitude_image,0,0,MagickTrue,exception);
3919 if (threshold_image == (Image *) NULL)
3920 ThrowDPCSimilarityException();
3921 status=BilevelImage(threshold_image,0.0,exception);
3922 if (status == MagickFalse)
3923 ThrowDPCSimilarityException();
3924 status=GetImageMean(threshold_image,&mean,&standard_deviation,exception);
3925 threshold_image=DestroyImage(threshold_image);
3926 if (status == MagickFalse)
3927 ThrowDPCSimilarityException();
3928 edge_factor=MagickSafeReciprocal(QuantumScale*mean*reconstruct->columns*
3929 reconstruct->rows)+QuantumScale;
3930 /*
3931 Divide X and Y derivitives of reference image by magnitude.
3932 */
3933 trx_image=SIMDivideByMagnitude(rx_image,magnitude_image,image,exception);
3934 rx_image=DestroyImage(rx_image);
3935 if (trx_image == (Image *) NULL)
3936 ThrowDPCSimilarityException();
3937 rx_image=trx_image;
3938 try_image=SIMDivideByMagnitude(ry_image,magnitude_image,image,exception);
3939 magnitude_image=DestroyImage(magnitude_image);
3940 ry_image=DestroyImage(ry_image);
3941 if (try_image == (Image *) NULL)
3942 ThrowDPCSimilarityException();
3943 ry_image=try_image;
3944 /*
3945 Compute X and Y derivatives of image.
3946 */
3947 (void) SetImageVirtualPixelMethod(target_image,EdgeVirtualPixelMethod,
3948 exception);
3949 tx_image=SIMDerivativeImage(target_image,"Sobel",exception);
3950 if (tx_image == (Image *) NULL)
3951 ThrowDPCSimilarityException();
3952 ty_image=SIMDerivativeImage(target_image,"Sobel:90",exception);
3953 target_image=DestroyImage(target_image);
3954 if (ty_image == (Image *) NULL)
3955 ThrowDPCSimilarityException();
3956 /*
3957 Compute magnitude of derivatives.
3958 */
3959 magnitude_image=SIMMagnitudeImage(tx_image,ty_image,exception);
3960 if (magnitude_image == (Image *) NULL)
3961 ThrowDPCSimilarityException();
3962 /*
3963 Divide Lx and Ly by magnitude.
3964 */
3965 trx_image=SIMDivideByMagnitude(tx_image,magnitude_image,image,exception);
3966 tx_image=DestroyImage(tx_image);
3967 if (trx_image == (Image *) NULL)
3968 ThrowDPCSimilarityException();
3969 tx_image=trx_image;
3970 try_image=SIMDivideByMagnitude(ty_image,magnitude_image,image,exception);
3971 ty_image=DestroyImage(ty_image);
3972 magnitude_image=DestroyImage(magnitude_image);
3973 if (try_image == (Image *) NULL)
3974 ThrowDPCSimilarityException();
3975 ty_image=try_image;
3976 /*
3977 Compute the cross correlation of the test and reference images.
3978 */
3979 trx_image=SIMCrossCorrelationImage(tx_image,rx_image,exception);
3980 rx_image=DestroyImage(rx_image);
3981 tx_image=DestroyImage(tx_image);
3982 if (trx_image == (Image *) NULL)
3983 ThrowDPCSimilarityException();
3984 try_image=SIMCrossCorrelationImage(ty_image,ry_image,exception);
3985 ry_image=DestroyImage(ry_image);
3986 ty_image=DestroyImage(ty_image);
3987 if (try_image == (Image *) NULL)
3988 ThrowDPCSimilarityException();
3989 /*
3990 Evaluate dot product correlation image.
3991 */
3992 (void) SetImageArtifact(try_image,"compose:clamp","false");
3993 status=CompositeImage(trx_image,try_image,PlusCompositeOp,MagickTrue,0,0,
3994 exception);
3995 try_image=DestroyImage(try_image);
3996 if (status == MagickFalse)
3997 ThrowDPCSimilarityException();
3998 status=SIMMultiplyImage(trx_image,edge_factor,
3999 (const ChannelStatistics *) NULL,exception);
4000 if (status == MagickFalse)
4001 ThrowDPCSimilarityException();
4002 /*
4003 Crop results.
4004 */
4005 SetGeometry(image,&geometry);
4006 geometry.width=image->columns;
4007 geometry.height=image->rows;
4008 (void) ResetImagePage(trx_image,"0x0+0+0");
4009 dot_product_image=CropImage(trx_image,&geometry,exception);
4010 trx_image=DestroyImage(trx_image);
4011 if (dot_product_image == (Image *) NULL)
4012 ThrowDPCSimilarityException();
4013 (void) ResetImagePage(dot_product_image,"0x0+0+0");
4014 /*
4015 Identify the maxima value in the image and its location.
4016 */
4017 status=GrayscaleImage(dot_product_image,AveragePixelIntensityMethod,
4018 exception);
4019 if (status == MagickFalse)
4020 ThrowDPCSimilarityException();
4021 dot_product_image->depth=32;
4022 dot_product_image->colorspace=GRAYColorspace;
4023 dot_product_image->alpha_trait=UndefinedPixelTrait;
4024 status=SIMFilterImageNaNs(dot_product_image,exception);
4025 if (status == MagickFalse)
4026 ThrowDPCSimilarityException();
4027 status=SIMMaximaImage(dot_product_image,&maxima,offset,exception);
4028 if (status == MagickFalse)
4029 ThrowDPCSimilarityException();
4030 if ((QuantumScale*maxima) > 1.0)
4031 {
4032 status=SIMMultiplyImage(dot_product_image,1.0/(QuantumScale*maxima),
4033 (const ChannelStatistics *) NULL,exception);
4034 maxima=(double) QuantumRange;
4035 }
4036 *similarity_metric=QuantumScale*maxima;
4037 return(dot_product_image);
4038}
4039
4040static Image *MSESimilarityImage(const Image *image,const Image *reconstruct,
4041 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4042{
4043#define ThrowMSESimilarityException() \
4044{ \
4045 if (alpha_image != (Image *) NULL) \
4046 alpha_image=DestroyImage(alpha_image); \
4047 if (beta_image != (Image *) NULL) \
4048 beta_image=DestroyImage(beta_image); \
4049 if (channel_statistics != (ChannelStatistics *) NULL) \
4050 channel_statistics=(ChannelStatistics *) \
4051 RelinquishMagickMemory(channel_statistics); \
4052 if (mean_image != (Image *) NULL) \
4053 mean_image=DestroyImage(mean_image); \
4054 if (mse_image != (Image *) NULL) \
4055 mse_image=DestroyImage(mse_image); \
4056 if (reconstruct_image != (Image *) NULL) \
4057 reconstruct_image=DestroyImage(reconstruct_image); \
4058 if (sum_image != (Image *) NULL) \
4059 sum_image=DestroyImage(sum_image); \
4060 if (alpha_image != (Image *) NULL) \
4061 alpha_image=DestroyImage(alpha_image); \
4062 return((Image *) NULL); \
4063}
4064
4065 ChannelStatistics
4066 *channel_statistics = (ChannelStatistics *) NULL;
4067
4068 double
4069 minima = 0.0;
4070
4071 Image
4072 *alpha_image = (Image *) NULL,
4073 *beta_image = (Image *) NULL,
4074 *mean_image = (Image *) NULL,
4075 *mse_image = (Image *) NULL,
4076 *reconstruct_image = (Image *) NULL,
4077 *sum_image = (Image *) NULL,
4078 *target_image = (Image *) NULL;
4079
4080 MagickBooleanType
4081 status = MagickTrue;
4082
4083 RectangleInfo
4084 geometry;
4085
4086 /*
4087 MSE correlation-based image similarity using FFT local statistics.
4088 */
4089 target_image=SIMSquareImage(image,exception);
4090 if (target_image == (Image *) NULL)
4091 ThrowMSESimilarityException();
4092 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
4093 if (reconstruct_image == (Image *) NULL)
4094 ThrowMSESimilarityException();
4095 /*
4096 Create (U * test)/# pixels.
4097 */
4098 alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
4099 exception);
4100 target_image=DestroyImage(target_image);
4101 if (alpha_image == (Image *) NULL)
4102 ThrowMSESimilarityException();
4103 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(double)
4104 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
4105 if (status == MagickFalse)
4106 ThrowMSESimilarityException();
4107 /*
4108 Create 2*(test * reconstruction)# pixels.
4109 */
4110 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
4111 MagickTrue,0,0,exception);
4112 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
4113 reconstruct_image=DestroyImage(reconstruct_image);
4114 if (beta_image == (Image *) NULL)
4115 ThrowMSESimilarityException();
4116 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(double)
4117 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
4118 if (status == MagickFalse)
4119 ThrowMSESimilarityException();
4120 /*
4121 Mean of reconstruction squared.
4122 */
4123 sum_image=SIMSquareImage(reconstruct,exception);
4124 if (sum_image == (Image *) NULL)
4125 ThrowMSESimilarityException();
4126 channel_statistics=GetImageStatistics(sum_image,exception);
4127 if (channel_statistics == (ChannelStatistics *) NULL)
4128 ThrowMSESimilarityException();
4129 status=SetImageStorageClass(sum_image,DirectClass,exception);
4130 if (status == MagickFalse)
4131 ThrowMSESimilarityException();
4132 status=SIMSetImageMean(sum_image,channel_statistics,exception);
4133 channel_statistics=(ChannelStatistics *)
4134 RelinquishMagickMemory(channel_statistics);
4135 if (status == MagickFalse)
4136 ThrowMSESimilarityException();
4137 /*
4138 Create mean image.
4139 */
4140 AppendImageToList(&sum_image,alpha_image);
4141 AppendImageToList(&sum_image,beta_image);
4142 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
4143 if (mean_image == (Image *) NULL)
4144 ThrowMSESimilarityException();
4145 sum_image=DestroyImage(sum_image);
4146 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
4147 if (status == MagickFalse)
4148 ThrowMSESimilarityException();
4149 /*
4150 Crop to difference of reconstruction and test images.
4151 */
4152 SetGeometry(image,&geometry);
4153 geometry.width=image->columns;
4154 geometry.height=image->rows;
4155 (void) ResetImagePage(mean_image,"0x0+0+0");
4156 mse_image=CropImage(mean_image,&geometry,exception);
4157 mean_image=DestroyImage(mean_image);
4158 if (mse_image == (Image *) NULL)
4159 ThrowMSESimilarityException();
4160 /*
4161 Identify the minima value in the correlation image and its location.
4162 */
4163 (void) ResetImagePage(mse_image,"0x0+0+0");
4164 (void) ClampImage(mse_image,exception);
4165 mse_image->depth=32;
4166 mse_image->colorspace=GRAYColorspace;
4167 mse_image->alpha_trait=UndefinedPixelTrait;
4168 status=SIMMinimaImage(mse_image,&minima,offset,exception);
4169 if (status == MagickFalse)
4170 ThrowMSESimilarityException();
4171 status=NegateImage(mse_image,MagickFalse,exception);
4172 if (status == MagickFalse)
4173 ThrowMSESimilarityException();
4174 alpha_image=DestroyImage(alpha_image);
4175 beta_image=DestroyImage(beta_image);
4176 if ((QuantumScale*minima) < FLT_EPSILON)
4177 minima=0.0;
4178 *similarity_metric=QuantumScale*minima;
4179 return(mse_image);
4180}
4181
4182static Image *NCCSimilarityImage(const Image *image,const Image *reconstruct,
4183 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4184{
4185#define ThrowNCCSimilarityException() \
4186{ \
4187 if (alpha_image != (Image *) NULL) \
4188 alpha_image=DestroyImage(alpha_image); \
4189 if (beta_image != (Image *) NULL) \
4190 beta_image=DestroyImage(beta_image); \
4191 if (channel_statistics != (ChannelStatistics *) NULL) \
4192 channel_statistics=(ChannelStatistics *) \
4193 RelinquishMagickMemory(channel_statistics); \
4194 if (correlation_image != (Image *) NULL) \
4195 correlation_image=DestroyImage(correlation_image); \
4196 if (divide_image != (Image *) NULL) \
4197 divide_image=DestroyImage(divide_image); \
4198 if (ncc_image != (Image *) NULL) \
4199 ncc_image=DestroyImage(ncc_image); \
4200 if (normalize_image != (Image *) NULL) \
4201 normalize_image=DestroyImage(normalize_image); \
4202 if (reconstruct_image != (Image *) NULL) \
4203 reconstruct_image=DestroyImage(reconstruct_image); \
4204 if (target_image != (Image *) NULL) \
4205 target_image=DestroyImage(target_image); \
4206 if (variance_image != (Image *) NULL) \
4207 variance_image=DestroyImage(variance_image); \
4208 return((Image *) NULL); \
4209}
4210
4211 ChannelStatistics
4212 *channel_statistics = (ChannelStatistics *) NULL;
4213
4214 double
4215 maxima = 0.0;
4216
4217 Image
4218 *alpha_image = (Image *) NULL,
4219 *beta_image = (Image *) NULL,
4220 *correlation_image = (Image *) NULL,
4221 *divide_image = (Image *) NULL,
4222 *ncc_image = (Image *) NULL,
4223 *normalize_image = (Image *) NULL,
4224 *reconstruct_image = (Image *) NULL,
4225 *target_image = (Image *) NULL,
4226 *variance_image = (Image *) NULL;
4227
4228 MagickBooleanType
4229 status = MagickTrue;
4230
4231 RectangleInfo
4232 geometry;
4233
4234 /*
4235 NCC correlation-based image similarity with FFT local statistics.
4236 */
4237 target_image=SIMSquareImage(image,exception);
4238 if (target_image == (Image *) NULL)
4239 ThrowNCCSimilarityException();
4240 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
4241 if (reconstruct_image == (Image *) NULL)
4242 ThrowNCCSimilarityException();
4243 /*
4244 Compute the cross correlation of the test and reconstruction images.
4245 */
4246 alpha_image=SIMCrossCorrelationImage(target_image,reconstruct_image,
4247 exception);
4248 target_image=DestroyImage(target_image);
4249 if (alpha_image == (Image *) NULL)
4250 ThrowNCCSimilarityException();
4251 status=SIMMultiplyImage(alpha_image,(double) QuantumRange*
4252 reconstruct->columns*reconstruct->rows,(const ChannelStatistics *) NULL,
4253 exception);
4254 if (status == MagickFalse)
4255 ThrowNCCSimilarityException();
4256 /*
4257 Compute the cross correlation of the source and reconstruction images.
4258 */
4259 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
4260 reconstruct_image=DestroyImage(reconstruct_image);
4261 if (beta_image == (Image *) NULL)
4262 ThrowNCCSimilarityException();
4263 target_image=SIMSquareImage(beta_image,exception);
4264 beta_image=DestroyImage(beta_image);
4265 if (target_image == (Image *) NULL)
4266 ThrowNCCSimilarityException();
4267 status=SIMMultiplyImage(target_image,(double) QuantumRange,
4268 (const ChannelStatistics *) NULL,exception);
4269 if (status == MagickFalse)
4270 ThrowNCCSimilarityException();
4271 /*
4272 Compute the variance of the two images.
4273 */
4274 variance_image=SIMVarianceImage(alpha_image,target_image,exception);
4275 target_image=DestroyImage(target_image);
4276 alpha_image=DestroyImage(alpha_image);
4277 if (variance_image == (Image *) NULL)
4278 ThrowNCCSimilarityException();
4279 /*
4280 Subtract the image mean.
4281 */
4282 channel_statistics=GetImageStatistics(reconstruct,exception);
4283 if (channel_statistics == (ChannelStatistics *) NULL)
4284 ThrowNCCSimilarityException();
4285 status=SIMMultiplyImage(variance_image,1.0,channel_statistics,exception);
4286 if (status == MagickFalse)
4287 ThrowNCCSimilarityException();
4288 normalize_image=SIMSubtractImageMean(image,reconstruct,channel_statistics,
4289 exception);
4290 channel_statistics=(ChannelStatistics *)
4291 RelinquishMagickMemory(channel_statistics);
4292 if (normalize_image == (Image *) NULL)
4293 ThrowNCCSimilarityException();
4294 correlation_image=SIMCrossCorrelationImage(image,normalize_image,exception);
4295 normalize_image=DestroyImage(normalize_image);
4296 if (correlation_image == (Image *) NULL)
4297 ThrowNCCSimilarityException();
4298 /*
4299 Divide the two images.
4300 */
4301 divide_image=SIMDivideImage(correlation_image,variance_image,exception);
4302 correlation_image=DestroyImage(correlation_image);
4303 variance_image=DestroyImage(variance_image);
4304 if (divide_image == (Image *) NULL)
4305 ThrowNCCSimilarityException();
4306 /*
4307 Crop padding.
4308 */
4309 SetGeometry(image,&geometry);
4310 geometry.width=image->columns;
4311 geometry.height=image->rows;
4312 (void) ResetImagePage(divide_image,"0x0+0+0");
4313 ncc_image=CropImage(divide_image,&geometry,exception);
4314 divide_image=DestroyImage(divide_image);
4315 if (ncc_image == (Image *) NULL)
4316 ThrowNCCSimilarityException();
4317 /*
4318 Identify the maxima value in the image and its location.
4319 */
4320 (void) ResetImagePage(ncc_image,"0x0+0+0");
4321 status=GrayscaleImage(ncc_image,AveragePixelIntensityMethod,exception);
4322 if (status == MagickFalse)
4323 ThrowNCCSimilarityException();
4324 ncc_image->depth=32;
4325 ncc_image->colorspace=GRAYColorspace;
4326 ncc_image->alpha_trait=UndefinedPixelTrait;
4327 status=SIMMaximaImage(ncc_image,&maxima,offset,exception);
4328 if (status == MagickFalse)
4329 ThrowNCCSimilarityException();
4330 if ((QuantumScale*maxima) > 1.0)
4331 {
4332 status=SIMMultiplyImage(ncc_image,1.0/(QuantumScale*maxima),
4333 (const ChannelStatistics *) NULL,exception);
4334 maxima=(double) QuantumRange;
4335 }
4336 *similarity_metric=QuantumScale*maxima;
4337 return(ncc_image);
4338}
4339
4340static Image *PhaseSimilarityImage(const Image *image,const Image *reconstruct,
4341 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4342{
4343#define ThrowPhaseSimilarityException() \
4344{ \
4345 if (correlation_image != (Image *) NULL) \
4346 correlation_image=DestroyImage(correlation_image); \
4347 if (fft_images != (Image *) NULL) \
4348 fft_images=DestroyImageList(fft_images); \
4349 if (gamma_image != (Image *) NULL) \
4350 gamma_image=DestroyImage(gamma_image); \
4351 if (magnitude_image != (Image *) NULL) \
4352 magnitude_image=DestroyImage(magnitude_image); \
4353 if (phase_image != (Image *) NULL) \
4354 phase_image=DestroyImage(phase_image); \
4355 if (reconstruct_image != (Image *) NULL) \
4356 reconstruct_image=DestroyImage(reconstruct_image); \
4357 if (reconstruct_magnitude != (Image *) NULL) \
4358 reconstruct_magnitude=DestroyImage(reconstruct_magnitude); \
4359 if (target_image != (Image *) NULL) \
4360 target_image=DestroyImage(target_image); \
4361 if (test_magnitude != (Image *) NULL) \
4362 test_magnitude=DestroyImage(test_magnitude); \
4363 return((Image *) NULL); \
4364}
4365
4366 double
4367 maxima = 0.0;
4368
4369 Image
4370 *correlation_image = (Image *) NULL,
4371 *fft_images = (Image *) NULL,
4372 *gamma_image = (Image *) NULL,
4373 *magnitude_image = (Image *) NULL,
4374 *phase_image = (Image *) NULL,
4375 *reconstruct_image = (Image *) NULL,
4376 *reconstruct_magnitude = (Image *) NULL,
4377 *target_image = (Image *) NULL,
4378 *test_magnitude = (Image *) NULL;
4379
4380 MagickBooleanType
4381 status = MagickTrue;
4382
4383 RectangleInfo
4384 geometry;
4385
4386 /*
4387 Phase correlation-based image similarity using FFT local statistics.
4388 */
4389 target_image=CloneImage(image,0,0,MagickTrue,exception);
4390 if (target_image == (Image *) NULL)
4391 ThrowPhaseSimilarityException();
4392 (void) ResetImagePage(target_image,"0x0+0+0");
4393 GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
4394 &target_image->background_color);
4395 status=SetImageExtent(target_image,2*(size_t) ceil((double) image->columns/
4396 2.0),2*(size_t) ceil((double) image->rows/2.0),exception);
4397 if (status == MagickFalse)
4398 ThrowPhaseSimilarityException();
4399 /*
4400 Compute the cross correlation of the test and reconstruct magnitudes.
4401 */
4402 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
4403 if (reconstruct_image == (Image *) NULL)
4404 ThrowPhaseSimilarityException();
4405 (void) ResetImagePage(reconstruct_image,"0x0+0+0");
4406 GetPixelInfoRGBA((Quantum) 0,(Quantum) 0,(Quantum) 0,(Quantum) 0,
4407 &reconstruct_image->background_color);
4408 status=SetImageExtent(reconstruct_image,2*(size_t) ceil((double)
4409 image->columns/2.0),2*(size_t) ceil((double) image->rows/2.0),exception);
4410 if (status == MagickFalse)
4411 ThrowPhaseSimilarityException();
4412 /*
4413 Evaluate phase coorelation image and divide by the product magnitude.
4414 */
4415 (void) SetImageArtifact(target_image,"fourier:normalize","inverse");
4416 fft_images=ForwardFourierTransformImage(target_image,MagickTrue,exception);
4417 if (fft_images == (Image *) NULL)
4418 ThrowPhaseSimilarityException();
4419 test_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
4420 fft_images=DestroyImageList(fft_images);
4421 if (test_magnitude == (Image *) NULL)
4422 ThrowPhaseSimilarityException();
4423 (void) SetImageArtifact(reconstruct_image,"fourier:normalize","inverse");
4424 fft_images=ForwardFourierTransformImage(reconstruct_image,MagickTrue,
4425 exception);
4426 if (fft_images == (Image *) NULL)
4427 ThrowPhaseSimilarityException();
4428 reconstruct_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
4429 fft_images=DestroyImageList(fft_images);
4430 if (reconstruct_magnitude == (Image *) NULL)
4431 ThrowPhaseSimilarityException();
4432 magnitude_image=CloneImage(reconstruct_magnitude,0,0,MagickTrue,exception);
4433 if (magnitude_image == (Image *) NULL)
4434 ThrowPhaseSimilarityException();
4435 DisableCompositeClampUnlessSpecified(magnitude_image);
4436 (void) CompositeImage(magnitude_image,test_magnitude,MultiplyCompositeOp,
4437 MagickTrue,0,0,exception);
4438 /*
4439 Compute the cross correlation of the test and reconstruction images.
4440 */
4441 correlation_image=SIMPhaseCorrelationImage(target_image,reconstruct_image,
4442 magnitude_image,exception);
4443 target_image=DestroyImage(target_image);
4444 reconstruct_image=DestroyImage(reconstruct_image);
4445 test_magnitude=DestroyImage(test_magnitude);
4446 reconstruct_magnitude=DestroyImage(reconstruct_magnitude);
4447 if (correlation_image == (Image *) NULL)
4448 ThrowPhaseSimilarityException();
4449 /*
4450 Identify the maxima value in the image and its location.
4451 */
4452 gamma_image=CloneImage(correlation_image,0,0,MagickTrue,exception);
4453 correlation_image=DestroyImage(correlation_image);
4454 if (gamma_image == (Image *) NULL)
4455 ThrowPhaseSimilarityException();
4456 /*
4457 Crop padding.
4458 */
4459 SetGeometry(image,&geometry);
4460 geometry.width=image->columns;
4461 geometry.height=image->rows;
4462 (void) ResetImagePage(gamma_image,"0x0+0+0");
4463 phase_image=CropImage(gamma_image,&geometry,exception);
4464 gamma_image=DestroyImage(gamma_image);
4465 if (phase_image == (Image *) NULL)
4466 ThrowPhaseSimilarityException();
4467 (void) ResetImagePage(phase_image,"0x0+0+0");
4468 /*
4469 Identify the maxima value in the correlation image and its location.
4470 */
4471 status=GrayscaleImage(phase_image,AveragePixelIntensityMethod,exception);
4472 if (status == MagickFalse)
4473 ThrowPhaseSimilarityException();
4474 phase_image->depth=32;
4475 phase_image->colorspace=GRAYColorspace;
4476 phase_image->alpha_trait=UndefinedPixelTrait;
4477 status=SIMFilterImageNaNs(phase_image,exception);
4478 if (status == MagickFalse)
4479 ThrowPhaseSimilarityException();
4480 status=SIMMaximaImage(phase_image,&maxima,offset,exception);
4481 if (status == MagickFalse)
4482 ThrowPhaseSimilarityException();
4483 magnitude_image=DestroyImage(magnitude_image);
4484 if ((QuantumScale*maxima) > 1.0)
4485 {
4486 status=SIMMultiplyImage(phase_image,1.0/(QuantumScale*maxima),
4487 (const ChannelStatistics *) NULL,exception);
4488 maxima=(double) QuantumRange;
4489 }
4490 *similarity_metric=QuantumScale*maxima;
4491 return(phase_image);
4492}
4493
4494static Image *PSNRSimilarityImage(const Image *image,const Image *reconstruct,
4495 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4496{
4497 Image
4498 *psnr_image = (Image *) NULL;
4499
4500 psnr_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
4501 exception);
4502 if (psnr_image == (Image *) NULL)
4503 return(psnr_image);
4504 *similarity_metric=10.0*MagickSafeLog10(MagickSafeReciprocal(
4505 *similarity_metric))/MagickSafePSNRRecipicol(10.0);
4506 return(psnr_image);
4507}
4508
4509static Image *RMSESimilarityImage(const Image *image,const Image *reconstruct,
4510 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4511{
4512 Image
4513 *rmse_image = (Image *) NULL;
4514
4515 rmse_image=MSESimilarityImage(image,reconstruct,offset,similarity_metric,
4516 exception);
4517 if (rmse_image == (Image *) NULL)
4518 return(rmse_image);
4519 *similarity_metric=sqrt(*similarity_metric);
4520 return(rmse_image);
4521}
4522#endif
4523
4524static double GetSimilarityMetric(const Image *image,
4525 const Image *reconstruct_image,const MetricType metric,
4526 const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo *exception)
4527{
4528 double
4529 *channel_similarity,
4530 similarity = 0.0;
4531
4532 ExceptionInfo
4533 *sans_exception = AcquireExceptionInfo();
4534
4535 Image
4536 *similarity_image;
4537
4538 MagickBooleanType
4539 status = MagickTrue;
4540
4541 RectangleInfo
4542 geometry;
4543
4544 size_t
4545 length = MaxPixelChannels+1UL;
4546
4547 SetGeometry(reconstruct_image,&geometry);
4548 geometry.x=x_offset;
4549 geometry.y=y_offset;
4550 similarity_image=CropImage(image,&geometry,sans_exception);
4551 sans_exception=DestroyExceptionInfo(sans_exception);
4552 if (similarity_image == (Image *) NULL)
4553 return(NAN);
4554 /*
4555 Get image distortion.
4556 */
4557 channel_similarity=(double *) AcquireQuantumMemory(length,
4558 sizeof(*channel_similarity));
4559 if (channel_similarity == (double *) NULL)
4560 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
4561 (void) memset(channel_similarity,0,length*sizeof(*channel_similarity));
4562 switch (metric)
4563 {
4564 case AbsoluteErrorMetric:
4565 {
4566 status=GetAESimilarity(similarity_image,reconstruct_image,
4567 channel_similarity,exception);
4568 break;
4569 }
4570 case DotProductCorrelationErrorMetric:
4571 {
4572 status=GetDPCSimilarity(similarity_image,reconstruct_image,
4573 channel_similarity,exception);
4574 break;
4575 }
4576 case FuzzErrorMetric:
4577 {
4578 status=GetFUZZSimilarity(similarity_image,reconstruct_image,
4579 channel_similarity,exception);
4580 break;
4581 }
4582 case MeanAbsoluteErrorMetric:
4583 {
4584 status=GetMAESimilarity(similarity_image,reconstruct_image,
4585 channel_similarity,exception);
4586 break;
4587 }
4588 case MeanErrorPerPixelErrorMetric:
4589 {
4590 status=GetMEPPSimilarity(similarity_image,reconstruct_image,
4591 channel_similarity,exception);
4592 break;
4593 }
4594 case MeanSquaredErrorMetric:
4595 {
4596 status=GetMSESimilarity(similarity_image,reconstruct_image,
4597 channel_similarity,exception);
4598 break;
4599 }
4600 case NormalizedCrossCorrelationErrorMetric:
4601 {
4602 status=GetNCCSimilarity(similarity_image,reconstruct_image,
4603 channel_similarity,exception);
4604 break;
4605 }
4606 case PeakAbsoluteErrorMetric:
4607 {
4608 status=GetPASimilarity(similarity_image,reconstruct_image,
4609 channel_similarity,exception);
4610 break;
4611 }
4612 case PeakSignalToNoiseRatioErrorMetric:
4613 {
4614 status=GetPSNRSimilarity(similarity_image,reconstruct_image,
4615 channel_similarity,exception);
4616 break;
4617 }
4618 case PerceptualHashErrorMetric:
4619 {
4620 status=GetPHASHSimilarity(similarity_image,reconstruct_image,
4621 channel_similarity,exception);
4622 break;
4623 }
4624 case PhaseCorrelationErrorMetric:
4625 {
4626 status=GetPHASESimilarity(similarity_image,reconstruct_image,
4627 channel_similarity,exception);
4628 break;
4629 }
4630 case PixelDifferenceCountErrorMetric:
4631 {
4632 status=GetPDCSimilarity(similarity_image,reconstruct_image,
4633 channel_similarity,exception);
4634 break;
4635 }
4636 case RootMeanSquaredErrorMetric:
4637 case UndefinedErrorMetric:
4638 default:
4639 {
4640 status=GetRMSESimilarity(similarity_image,reconstruct_image,
4641 channel_similarity,exception);
4642 break;
4643 }
4644 case StructuralDissimilarityErrorMetric:
4645 {
4646 status=GetDSSIMSimilarity(similarity_image,reconstruct_image,
4647 channel_similarity,exception);
4648 break;
4649 }
4650 case StructuralSimilarityErrorMetric:
4651 {
4652 status=GetSSIMSimularity(similarity_image,reconstruct_image,
4653 channel_similarity,exception);
4654 break;
4655 }
4656 }
4657 similarity_image=DestroyImage(similarity_image);
4658 similarity=channel_similarity[CompositePixelChannel];
4659 channel_similarity=(double *) RelinquishMagickMemory(channel_similarity);
4660 if (status == MagickFalse)
4661 return(NAN);
4662 return(similarity);
4663}
4664
4665MagickExport Image *SimilarityImage(const Image *image,const Image *reconstruct,
4666 const MetricType metric,const double similarity_threshold,
4667 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4668{
4669#define SimilarityImageTag "Similarity/Image"
4670
4671 typedef struct
4672 {
4673 double
4674 similarity;
4675
4676 ssize_t
4677 x,
4678 y;
4679 } SimilarityInfo;
4680
4681 CacheView
4682 *similarity_view;
4683
4684 Image
4685 *similarity_image = (Image *) NULL;
4686
4687 MagickBooleanType
4688 status = MagickTrue;
4689
4690 MagickOffsetType
4691 progress = 0;
4692
4693 SimilarityInfo
4694 similarity_info = { 0.0, 0, 0 };
4695
4696 ssize_t
4697 y;
4698
4699 assert(image != (const Image *) NULL);
4700 assert(image->signature == MagickCoreSignature);
4701 assert(exception != (ExceptionInfo *) NULL);
4702 assert(exception->signature == MagickCoreSignature);
4703 assert(offset != (RectangleInfo *) NULL);
4704 if (IsEventLogging() != MagickFalse)
4705 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4706 SetGeometry(reconstruct,offset);
4707 *similarity_metric=0.0;
4708 offset->x=0;
4709 offset->y=0;
4710#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
4711{
4712 const char *artifact = GetImageArtifact(image,"compare:frequency-domain");
4713 if (artifact == (const char *) NULL)
4714 artifact=GetImageArtifact(image,"compare:accelerate-ncc");
4715 if (((artifact == (const char *) NULL) ||
4716 (IsStringTrue(artifact) != MagickFalse)) &&
4717 ((image->channels & ReadMaskChannel) == 0))
4718 switch (metric)
4719 {
4720 case DotProductCorrelationErrorMetric:
4721 {
4722 similarity_image=DPCSimilarityImage(image,reconstruct,offset,
4723 similarity_metric,exception);
4724 return(similarity_image);
4725 }
4726 case MeanSquaredErrorMetric:
4727 {
4728 similarity_image=MSESimilarityImage(image,reconstruct,offset,
4729 similarity_metric,exception);
4730 return(similarity_image);
4731 }
4732 case NormalizedCrossCorrelationErrorMetric:
4733 {
4734 similarity_image=NCCSimilarityImage(image,reconstruct,offset,
4735 similarity_metric,exception);
4736 return(similarity_image);
4737 }
4738 case PeakSignalToNoiseRatioErrorMetric:
4739 {
4740 similarity_image=PSNRSimilarityImage(image,reconstruct,offset,
4741 similarity_metric,exception);
4742 return(similarity_image);
4743 }
4744 case PhaseCorrelationErrorMetric:
4745 {
4746 similarity_image=PhaseSimilarityImage(image,reconstruct,offset,
4747 similarity_metric,exception);
4748 return(similarity_image);
4749 }
4750 case RootMeanSquaredErrorMetric:
4751 case UndefinedErrorMetric:
4752 {
4753 similarity_image=RMSESimilarityImage(image,reconstruct,offset,
4754 similarity_metric,exception);
4755 return(similarity_image);
4756 }
4757 default:
4758 break;
4759 }
4760}
4761#endif
4762 if ((image->columns < reconstruct->columns) ||
4763 (image->rows < reconstruct->rows))
4764 {
4765 (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
4766 "GeometryDoesNotContainImage","`%s'",image->filename);
4767 return((Image *) NULL);
4768 }
4769 similarity_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4770 exception);
4771 if (similarity_image == (Image *) NULL)
4772 return((Image *) NULL);
4773 similarity_image->depth=32;
4774 similarity_image->colorspace=GRAYColorspace;
4775 similarity_image->alpha_trait=UndefinedPixelTrait;
4776 status=SetImageStorageClass(similarity_image,DirectClass,exception);
4777 if (status == MagickFalse)
4778 return(DestroyImage(similarity_image));
4779 /*
4780 Measure similarity of reconstruction image against image.
4781 */
4782 similarity_info.similarity=GetSimilarityMetric(image,reconstruct,metric,
4783 similarity_info.x,similarity_info.y,exception);
4784 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
4785#if defined(MAGICKCORE_OPENMP_SUPPORT)
4786 #pragma omp parallel for schedule(static) shared(similarity_info,status) \
4787 magick_number_threads(image,reconstruct,similarity_image->rows,1)
4788#endif
4789 for (y=0; y < (ssize_t) similarity_image->rows; y++)
4790 {
4791 double
4792 similarity;
4793
4794 MagickBooleanType
4795 threshold_trigger = MagickFalse;
4796
4797 Quantum
4798 *magick_restrict q;
4799
4800 SimilarityInfo
4801 channel_info = similarity_info;
4802
4803 ssize_t
4804 x;
4805
4806 if (status == MagickFalse)
4807 continue;
4808 if (threshold_trigger != MagickFalse)
4809 continue;
4810 q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
4811 similarity_image->columns,1,exception);
4812 if (q == (Quantum *) NULL)
4813 {
4814 status=MagickFalse;
4815 continue;
4816 }
4817 for (x=0; x < (ssize_t) similarity_image->columns; x++)
4818 {
4819 ssize_t
4820 i;
4821
4822 similarity=GetSimilarityMetric((Image *) image,reconstruct,metric,x,y,
4823 exception);
4824 switch (metric)
4825 {
4826 case DotProductCorrelationErrorMetric:
4827 case NormalizedCrossCorrelationErrorMetric:
4828 case PeakSignalToNoiseRatioErrorMetric:
4829 case PhaseCorrelationErrorMetric:
4830 case StructuralSimilarityErrorMetric:
4831 {
4832 if (similarity <= channel_info.similarity)
4833 break;
4834 channel_info.similarity=similarity;
4835 channel_info.x=x;
4836 channel_info.y=y;
4837 break;
4838 }
4839 default:
4840 {
4841 if (similarity >= channel_info.similarity)
4842 break;
4843 channel_info.similarity=similarity;
4844 channel_info.x=x;
4845 channel_info.y=y;
4846 break;
4847 }
4848 }
4849 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4850 {
4851 PixelChannel channel = GetPixelChannelChannel(image,i);
4852 PixelTrait traits = GetPixelChannelTraits(image,channel);
4853 PixelTrait similarity_traits = GetPixelChannelTraits(similarity_image,
4854 channel);
4855 if (((traits & UpdatePixelTrait) == 0) ||
4856 ((similarity_traits & UpdatePixelTrait) == 0))
4857 continue;
4858 switch (metric)
4859 {
4860 case DotProductCorrelationErrorMetric:
4861 case NormalizedCrossCorrelationErrorMetric:
4862 case PeakSignalToNoiseRatioErrorMetric:
4863 case PhaseCorrelationErrorMetric:
4864 case StructuralSimilarityErrorMetric:
4865 {
4866 SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4867 QuantumRange*similarity),q);
4868 break;
4869 }
4870 default:
4871 {
4872 SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4873 QuantumRange*(1.0-similarity)),q);
4874 break;
4875 }
4876 }
4877 }
4878 q+=(ptrdiff_t) GetPixelChannels(similarity_image);
4879 }
4880#if defined(MAGICKCORE_OPENMP_SUPPORT)
4881 #pragma omp critical (MagickCore_GetSimilarityMetric)
4882#endif
4883 switch (metric)
4884 {
4885 case DotProductCorrelationErrorMetric:
4886 case NormalizedCrossCorrelationErrorMetric:
4887 case PeakSignalToNoiseRatioErrorMetric:
4888 case PhaseCorrelationErrorMetric:
4889 case StructuralSimilarityErrorMetric:
4890 {
4891 if (similarity_threshold != DefaultSimilarityThreshold)
4892 if (channel_info.similarity >= similarity_threshold)
4893 threshold_trigger=MagickTrue;
4894 if (channel_info.similarity >= similarity_info.similarity)
4895 similarity_info=channel_info;
4896 break;
4897 }
4898 default:
4899 {
4900 if (similarity_threshold != DefaultSimilarityThreshold)
4901 if (channel_info.similarity < similarity_threshold)
4902 threshold_trigger=MagickTrue;
4903 if (channel_info.similarity < similarity_info.similarity)
4904 similarity_info=channel_info;
4905 break;
4906 }
4907 }
4908 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
4909 status=MagickFalse;
4910 if (image->progress_monitor != (MagickProgressMonitor) NULL)
4911 {
4912 MagickBooleanType
4913 proceed;
4914
4915 progress++;
4916 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
4917 if (proceed == MagickFalse)
4918 status=MagickFalse;
4919 }
4920 }
4921 similarity_view=DestroyCacheView(similarity_view);
4922 if (status == MagickFalse)
4923 similarity_image=DestroyImage(similarity_image);
4924 *similarity_metric=similarity_info.similarity;
4925 if (fabs(*similarity_metric) < MagickEpsilon)
4926 *similarity_metric=0.0;
4927 offset->x=similarity_info.x;
4928 offset->y=similarity_info.y;
4929 (void) FormatImageProperty((Image *) image,"similarity","%.*g",
4930 GetMagickPrecision(),*similarity_metric);
4931 (void) FormatImageProperty((Image *) image,"similarity.offset.x","%.*g",
4932 GetMagickPrecision(),(double) offset->x);
4933 (void) FormatImageProperty((Image *) image,"similarity.offset.y","%.*g",
4934 GetMagickPrecision(),(double) offset->y);
4935 return(similarity_image);
4936}