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