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