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