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