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