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