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