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