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