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