MagickCore 7.1.1-46
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{
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
137 highlight,
138 lowlight,
139 masklight;
140
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{
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{
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{
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{
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{
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
941 *image_view,
942 *reconstruct_view;
943
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 rows=MagickMax(image->rows,reconstruct_image->rows);
1115 columns=MagickMax(image->columns,reconstruct_image->columns);
1116 image_view=AcquireVirtualCacheView(image,exception);
1117 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1118#if defined(MAGICKCORE_OPENMP_SUPPORT)
1119 #pragma omp parallel for schedule(static) shared(status) \
1120 magick_number_threads(image,image,rows,1)
1121#endif
1122 for (y=0; y < (ssize_t) rows; y++)
1123 {
1124 double
1125 channel_distortion[MaxPixelChannels+1];
1126
1127 const Quantum
1128 *magick_restrict p,
1129 *magick_restrict q;
1130
1131 ssize_t
1132 j,
1133 x;
1134
1135 if (status == MagickFalse)
1136 continue;
1137 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1138 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1139 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1140 {
1141 status=MagickFalse;
1142 continue;
1143 }
1144 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1145 for (x=0; x < (ssize_t) columns; x++)
1146 {
1147 double
1148 Da,
1149 Sa;
1150
1151 ssize_t
1152 i;
1153
1154 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1155 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1156 {
1157 p+=(ptrdiff_t) GetPixelChannels(image);
1158 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1159 continue;
1160 }
1161 Sa=QuantumScale*(double) GetPixelAlpha(image,p);
1162 Da=QuantumScale*(double) GetPixelAlpha(reconstruct_image,q);
1163 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1164 {
1165 double
1166 distance;
1167
1168 PixelChannel channel = GetPixelChannelChannel(image,i);
1169 PixelTrait traits = GetPixelChannelTraits(image,channel);
1170 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1171 channel);
1172 if ((traits == UndefinedPixelTrait) ||
1173 (reconstruct_traits == UndefinedPixelTrait) ||
1174 ((reconstruct_traits & UpdatePixelTrait) == 0))
1175 continue;
1176 if (channel == AlphaPixelChannel)
1177 distance=QuantumScale*fabs((double) p[i]-(double)
1178 GetPixelChannel(reconstruct_image,channel,q));
1179 else
1180 distance=QuantumScale*fabs(Sa*(double) p[i]-Da*(double)
1181 GetPixelChannel(reconstruct_image,channel,q));
1182 if (distance > channel_distortion[i])
1183 channel_distortion[i]=distance;
1184 if (distance > channel_distortion[CompositePixelChannel])
1185 channel_distortion[CompositePixelChannel]=distance;
1186 }
1187 p+=(ptrdiff_t) GetPixelChannels(image);
1188 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1189 }
1190#if defined(MAGICKCORE_OPENMP_SUPPORT)
1191 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1192#endif
1193 for (j=0; j <= MaxPixelChannels; j++)
1194 if (channel_distortion[j] > distortion[j])
1195 distortion[j]=channel_distortion[j];
1196 }
1197 reconstruct_view=DestroyCacheView(reconstruct_view);
1198 image_view=DestroyCacheView(image_view);
1199 return(status);
1200}
1201
1202static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1203 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1204{
1205 MagickBooleanType
1206 status;
1207
1208 ssize_t
1209 i;
1210
1211 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1212 for (i=0; i <= MaxPixelChannels; i++)
1213 if (fabs(distortion[i]) >= MagickEpsilon)
1214 distortion[i]=(10.0*MagickLog10(distortion[i]))/48.1647;
1215 return(status);
1216}
1217
1218static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1219 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1220{
1222 *channel_phash,
1223 *reconstruct_phash;
1224
1225 const char
1226 *artifact;
1227
1228 ssize_t
1229 channel;
1230
1231 /*
1232 Compute perceptual hash in the sRGB colorspace.
1233 */
1234 channel_phash=GetImagePerceptualHash(image,exception);
1235 if (channel_phash == (ChannelPerceptualHash *) NULL)
1236 return(MagickFalse);
1237 reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1238 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1239 {
1240 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1241 channel_phash);
1242 return(MagickFalse);
1243 }
1244#if defined(MAGICKCORE_OPENMP_SUPPORT)
1245 #pragma omp parallel for schedule(static)
1246#endif
1247 for (channel=0; channel < MaxPixelChannels; channel++)
1248 {
1249 double
1250 difference;
1251
1252 ssize_t
1253 i;
1254
1255 difference=0.0;
1256 for (i=0; i < MaximumNumberOfImageMoments; i++)
1257 {
1258 double
1259 alpha,
1260 beta;
1261
1262 ssize_t
1263 j;
1264
1265 for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1266 {
1267 alpha=channel_phash[channel].phash[j][i];
1268 beta=reconstruct_phash[channel].phash[j][i];
1269 difference+=(beta-alpha)*(beta-alpha);
1270 }
1271 }
1272 distortion[channel]+=difference;
1273#if defined(MAGICKCORE_OPENMP_SUPPORT)
1274 #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1275#endif
1276 distortion[CompositePixelChannel]+=difference;
1277 }
1278 artifact=GetImageArtifact(image,"phash:normalize");
1279 if ((artifact != (const char *) NULL) &&
1280 (IsStringTrue(artifact) != MagickFalse))
1281 {
1282 ssize_t
1283 j;
1284
1285 for (j=0; j <= MaxPixelChannels; j++)
1286 distortion[j]=sqrt(distortion[j]/channel_phash[0].number_channels);
1287 }
1288 /*
1289 Free resources.
1290 */
1291 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1292 reconstruct_phash);
1293 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1294 return(MagickTrue);
1295}
1296
1297static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1298 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1299{
1300 MagickBooleanType
1301 status;
1302
1303 ssize_t
1304 i;
1305
1306 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1307 for (i=0; i <= MaxPixelChannels; i++)
1308 distortion[i]=sqrt(distortion[i]);
1309 return(status);
1310}
1311
1312static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image,
1313 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1314{
1315#define SSIMRadius 5.0
1316#define SSIMSigma 1.5
1317#define SSIMBlocksize 8
1318#define SSIMK1 0.01
1319#define SSIMK2 0.03
1320#define SSIML 1.0
1321
1322 CacheView
1323 *image_view,
1324 *reconstruct_view;
1325
1326 char
1327 geometry[MagickPathExtent];
1328
1329 const char
1330 *artifact;
1331
1332 double
1333 area,
1334 c1,
1335 c2,
1336 radius,
1337 sigma;
1338
1340 *kernel_info;
1341
1342 MagickBooleanType
1343 status;
1344
1345 size_t
1346 columns,
1347 rows;
1348
1349 ssize_t
1350 j,
1351 y;
1352
1353 /*
1354 Compute structural similarity index @
1355 https://en.wikipedia.org/wiki/Structural_similarity.
1356 */
1357 radius=SSIMRadius;
1358 artifact=GetImageArtifact(image,"compare:ssim-radius");
1359 if (artifact != (const char *) NULL)
1360 radius=StringToDouble(artifact,(char **) NULL);
1361 sigma=SSIMSigma;
1362 artifact=GetImageArtifact(image,"compare:ssim-sigma");
1363 if (artifact != (const char *) NULL)
1364 sigma=StringToDouble(artifact,(char **) NULL);
1365 (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1366 radius,sigma);
1367 kernel_info=AcquireKernelInfo(geometry,exception);
1368 if (kernel_info == (KernelInfo *) NULL)
1369 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1370 image->filename);
1371 c1=pow(SSIMK1*SSIML,2.0);
1372 artifact=GetImageArtifact(image,"compare:ssim-k1");
1373 if (artifact != (const char *) NULL)
1374 c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1375 c2=pow(SSIMK2*SSIML,2.0);
1376 artifact=GetImageArtifact(image,"compare:ssim-k2");
1377 if (artifact != (const char *) NULL)
1378 c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1379 status=MagickTrue;
1380 area=0.0;
1381 rows=MagickMax(image->rows,reconstruct_image->rows);
1382 columns=MagickMax(image->columns,reconstruct_image->columns);
1383 image_view=AcquireVirtualCacheView(image,exception);
1384 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1385#if defined(MAGICKCORE_OPENMP_SUPPORT)
1386 #pragma omp parallel for schedule(static,1) shared(area,distortion,status) \
1387 magick_number_threads(image,reconstruct_image,rows,1)
1388#endif
1389 for (y=0; y < (ssize_t) rows; y++)
1390 {
1391 const Quantum
1392 *magick_restrict p,
1393 *magick_restrict q;
1394
1395 double
1396 channel_distortion[MaxPixelChannels+1];
1397
1398 size_t
1399 local_area = 0;
1400
1401 ssize_t
1402 i,
1403 x;
1404
1405 if (status == MagickFalse)
1406 continue;
1407 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1408 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1409 kernel_info->height,exception);
1410 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1411 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1412 kernel_info->height,exception);
1413 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1414 {
1415 status=MagickFalse;
1416 continue;
1417 }
1418 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1419 for (x=0; x < (ssize_t) columns; x++)
1420 {
1421 const Quantum
1422 *magick_restrict reconstruct,
1423 *magick_restrict test;
1424
1425 double
1426 x_pixel_mu[MaxPixelChannels+1],
1427 x_pixel_sigma_squared[MaxPixelChannels+1],
1428 xy_sigma[MaxPixelChannels+1],
1429 y_pixel_mu[MaxPixelChannels+1],
1430 y_pixel_sigma_squared[MaxPixelChannels+1];
1431
1432 MagickRealType
1433 *k;
1434
1435 ssize_t
1436 v;
1437
1438 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1439 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1440 {
1441 p+=(ptrdiff_t) GetPixelChannels(image);
1442 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1443 continue;
1444 }
1445 (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1446 (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1447 (void) memset(xy_sigma,0,sizeof(xy_sigma));
1448 (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1449 (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1450 k=kernel_info->values;
1451 reconstruct=p;
1452 test=q;
1453 for (v=0; v < (ssize_t) kernel_info->height; v++)
1454 {
1455 ssize_t
1456 u;
1457
1458 for (u=0; u < (ssize_t) kernel_info->width; u++)
1459 {
1460 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1461 {
1462 double
1463 x_pixel,
1464 y_pixel;
1465
1466 PixelChannel channel = GetPixelChannelChannel(image,i);
1467 PixelTrait traits = GetPixelChannelTraits(image,channel);
1468 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1469 reconstruct_image,channel);
1470 if ((traits == UndefinedPixelTrait) ||
1471 (reconstruct_traits == UndefinedPixelTrait) ||
1472 ((reconstruct_traits & UpdatePixelTrait) == 0))
1473 continue;
1474 x_pixel=QuantumScale*(double) reconstruct[i];
1475 x_pixel_mu[i]+=(*k)*x_pixel;
1476 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1477 y_pixel=QuantumScale*(double)
1478 GetPixelChannel(reconstruct_image,channel,test);
1479 y_pixel_mu[i]+=(*k)*y_pixel;
1480 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1481 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1482 }
1483 k++;
1484 reconstruct+=(ptrdiff_t) GetPixelChannels(image);
1485 test+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1486 }
1487 reconstruct+=(ptrdiff_t) GetPixelChannels(image)*columns;
1488 test+=(ptrdiff_t) GetPixelChannels(reconstruct_image)*columns;
1489 }
1490 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1491 {
1492 double
1493 ssim,
1494 x_pixel_mu_squared,
1495 x_pixel_sigmas_squared,
1496 xy_mu,
1497 xy_sigmas,
1498 y_pixel_mu_squared,
1499 y_pixel_sigmas_squared;
1500
1501 PixelChannel channel = GetPixelChannelChannel(image,i);
1502 PixelTrait traits = GetPixelChannelTraits(image,channel);
1503 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1504 reconstruct_image,channel);
1505 if ((traits == UndefinedPixelTrait) ||
1506 (reconstruct_traits == UndefinedPixelTrait) ||
1507 ((reconstruct_traits & UpdatePixelTrait) == 0))
1508 continue;
1509 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1510 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1511 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1512 xy_sigmas=xy_sigma[i]-xy_mu;
1513 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1514 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1515 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1516 ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1517 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1518 channel_distortion[i]+=ssim;
1519 channel_distortion[CompositePixelChannel]+=ssim;
1520 }
1521 p+=(ptrdiff_t) GetPixelChannels(image);
1522 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1523 local_area++;
1524 }
1525#if defined(MAGICKCORE_OPENMP_SUPPORT)
1526 #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1527#endif
1528 {
1529 area+=local_area;
1530 for (i=0; i <= MaxPixelChannels; i++)
1531 distortion[i]+=channel_distortion[i];
1532 }
1533 }
1534 image_view=DestroyCacheView(image_view);
1535 reconstruct_view=DestroyCacheView(reconstruct_view);
1536 for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1537 {
1538 PixelChannel channel = GetPixelChannelChannel(image,j);
1539 PixelTrait traits = GetPixelChannelTraits(image,channel);
1540 if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1541 continue;
1542 distortion[j]*=PerceptibleReciprocal(area);
1543 }
1544 distortion[CompositePixelChannel]*=PerceptibleReciprocal(area);
1545 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1546 kernel_info=DestroyKernelInfo(kernel_info);
1547 return(status);
1548}
1549
1550static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image,
1551 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1552{
1553 MagickBooleanType
1554 status;
1555
1556 ssize_t
1557 i;
1558
1559 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1560 distortion,exception);
1561 for (i=0; i <= MaxPixelChannels; i++)
1562 distortion[i]=(1.0-distortion[i])/2.0;
1563 return(status);
1564}
1565
1566MagickExport MagickBooleanType GetImageDistortion(Image *image,
1567 const Image *reconstruct_image,const MetricType metric,double *distortion,
1568 ExceptionInfo *exception)
1569{
1570 double
1571 *channel_distortion;
1572
1573 MagickBooleanType
1574 status;
1575
1576 size_t
1577 length;
1578
1579 assert(image != (Image *) NULL);
1580 assert(image->signature == MagickCoreSignature);
1581 assert(reconstruct_image != (const Image *) NULL);
1582 assert(reconstruct_image->signature == MagickCoreSignature);
1583 assert(distortion != (double *) NULL);
1584 if (IsEventLogging() != MagickFalse)
1585 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1586 *distortion=0.0;
1587 /*
1588 Get image distortion.
1589 */
1590 length=MaxPixelChannels+1UL;
1591 channel_distortion=(double *) AcquireQuantumMemory(length,
1592 sizeof(*channel_distortion));
1593 if (channel_distortion == (double *) NULL)
1594 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1595 (void) memset(channel_distortion,0,length*
1596 sizeof(*channel_distortion));
1597 switch (metric)
1598 {
1599 case AbsoluteErrorMetric:
1600 {
1601 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1602 exception);
1603 break;
1604 }
1605 case FuzzErrorMetric:
1606 {
1607 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1608 exception);
1609 break;
1610 }
1611 case MeanAbsoluteErrorMetric:
1612 {
1613 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1614 channel_distortion,exception);
1615 break;
1616 }
1617 case MeanErrorPerPixelErrorMetric:
1618 {
1619 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1620 exception);
1621 break;
1622 }
1623 case MeanSquaredErrorMetric:
1624 {
1625 status=GetMeanSquaredDistortion(image,reconstruct_image,
1626 channel_distortion,exception);
1627 break;
1628 }
1629 case NormalizedCrossCorrelationErrorMetric:
1630 default:
1631 {
1632 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1633 channel_distortion,exception);
1634 break;
1635 }
1636 case PeakAbsoluteErrorMetric:
1637 {
1638 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1639 channel_distortion,exception);
1640 break;
1641 }
1642 case PeakSignalToNoiseRatioErrorMetric:
1643 {
1644 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1645 channel_distortion,exception);
1646 break;
1647 }
1648 case PerceptualHashErrorMetric:
1649 {
1650 status=GetPerceptualHashDistortion(image,reconstruct_image,
1651 channel_distortion,exception);
1652 break;
1653 }
1654 case RootMeanSquaredErrorMetric:
1655 {
1656 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1657 channel_distortion,exception);
1658 break;
1659 }
1660 case StructuralSimilarityErrorMetric:
1661 {
1662 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1663 channel_distortion,exception);
1664 break;
1665 }
1666 case StructuralDissimilarityErrorMetric:
1667 {
1668 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1669 channel_distortion,exception);
1670 break;
1671 }
1672 }
1673 *distortion=channel_distortion[CompositePixelChannel];
1674 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1675 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1676 *distortion);
1677 return(status);
1678}
1679
1680/*
1681%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1682% %
1683% %
1684% %
1685% G e t I m a g e D i s t o r t i o n s %
1686% %
1687% %
1688% %
1689%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1690%
1691% GetImageDistortions() compares the pixel channels of an image to a
1692% reconstructed image and returns the specified distortion metric for each
1693% channel.
1694%
1695% The format of the GetImageDistortions method is:
1696%
1697% double *GetImageDistortions(const Image *image,
1698% const Image *reconstruct_image,const MetricType metric,
1699% ExceptionInfo *exception)
1700%
1701% A description of each parameter follows:
1702%
1703% o image: the image.
1704%
1705% o reconstruct_image: the reconstruction image.
1706%
1707% o metric: the metric.
1708%
1709% o exception: return any errors or warnings in this structure.
1710%
1711*/
1712MagickExport double *GetImageDistortions(Image *image,
1713 const Image *reconstruct_image,const MetricType metric,
1714 ExceptionInfo *exception)
1715{
1716 double
1717 *channel_distortion;
1718
1719 MagickBooleanType
1720 status;
1721
1722 size_t
1723 length;
1724
1725 assert(image != (Image *) NULL);
1726 assert(image->signature == MagickCoreSignature);
1727 assert(reconstruct_image != (const Image *) NULL);
1728 assert(reconstruct_image->signature == MagickCoreSignature);
1729 if (IsEventLogging() != MagickFalse)
1730 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1731 /*
1732 Get image distortion.
1733 */
1734 length=MaxPixelChannels+1UL;
1735 channel_distortion=(double *) AcquireQuantumMemory(length,
1736 sizeof(*channel_distortion));
1737 if (channel_distortion == (double *) NULL)
1738 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1739 (void) memset(channel_distortion,0,length*
1740 sizeof(*channel_distortion));
1741 status=MagickTrue;
1742 switch (metric)
1743 {
1744 case AbsoluteErrorMetric:
1745 {
1746 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1747 exception);
1748 break;
1749 }
1750 case FuzzErrorMetric:
1751 {
1752 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1753 exception);
1754 break;
1755 }
1756 case MeanAbsoluteErrorMetric:
1757 {
1758 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1759 channel_distortion,exception);
1760 break;
1761 }
1762 case MeanErrorPerPixelErrorMetric:
1763 {
1764 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1765 exception);
1766 break;
1767 }
1768 case MeanSquaredErrorMetric:
1769 {
1770 status=GetMeanSquaredDistortion(image,reconstruct_image,
1771 channel_distortion,exception);
1772 break;
1773 }
1774 case NormalizedCrossCorrelationErrorMetric:
1775 default:
1776 {
1777 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1778 channel_distortion,exception);
1779 break;
1780 }
1781 case PeakAbsoluteErrorMetric:
1782 {
1783 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1784 channel_distortion,exception);
1785 break;
1786 }
1787 case PeakSignalToNoiseRatioErrorMetric:
1788 {
1789 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1790 channel_distortion,exception);
1791 break;
1792 }
1793 case PerceptualHashErrorMetric:
1794 {
1795 status=GetPerceptualHashDistortion(image,reconstruct_image,
1796 channel_distortion,exception);
1797 break;
1798 }
1799 case RootMeanSquaredErrorMetric:
1800 {
1801 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1802 channel_distortion,exception);
1803 break;
1804 }
1805 case StructuralSimilarityErrorMetric:
1806 {
1807 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1808 channel_distortion,exception);
1809 break;
1810 }
1811 case StructuralDissimilarityErrorMetric:
1812 {
1813 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1814 channel_distortion,exception);
1815 break;
1816 }
1817 }
1818 if (status == MagickFalse)
1819 {
1820 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1821 return((double *) NULL);
1822 }
1823 return(channel_distortion);
1824}
1825
1826/*
1827%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1828% %
1829% %
1830% %
1831% I s I m a g e s E q u a l %
1832% %
1833% %
1834% %
1835%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1836%
1837% IsImagesEqual() compare the pixels of two images and returns immediately
1838% if any pixel is not identical.
1839%
1840% The format of the IsImagesEqual method is:
1841%
1842% MagickBooleanType IsImagesEqual(const Image *image,
1843% const Image *reconstruct_image,ExceptionInfo *exception)
1844%
1845% A description of each parameter follows.
1846%
1847% o image: the image.
1848%
1849% o reconstruct_image: the reconstruction image.
1850%
1851% o exception: return any errors or warnings in this structure.
1852%
1853*/
1854MagickExport MagickBooleanType IsImagesEqual(const Image *image,
1855 const Image *reconstruct_image,ExceptionInfo *exception)
1856{
1857 CacheView
1858 *image_view,
1859 *reconstruct_view;
1860
1861 size_t
1862 columns,
1863 rows;
1864
1865 ssize_t
1866 y;
1867
1868 assert(image != (Image *) NULL);
1869 assert(image->signature == MagickCoreSignature);
1870 assert(reconstruct_image != (const Image *) NULL);
1871 assert(reconstruct_image->signature == MagickCoreSignature);
1872 rows=MagickMax(image->rows,reconstruct_image->rows);
1873 columns=MagickMax(image->columns,reconstruct_image->columns);
1874 image_view=AcquireVirtualCacheView(image,exception);
1875 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1876 for (y=0; y < (ssize_t) rows; y++)
1877 {
1878 const Quantum
1879 *magick_restrict p,
1880 *magick_restrict q;
1881
1882 ssize_t
1883 x;
1884
1885 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1886 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1887 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1888 break;
1889 for (x=0; x < (ssize_t) columns; x++)
1890 {
1891 ssize_t
1892 i;
1893
1894 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1895 {
1896 double
1897 distance;
1898
1899 PixelChannel channel = GetPixelChannelChannel(image,i);
1900 PixelTrait traits = GetPixelChannelTraits(image,channel);
1901 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1902 channel);
1903 if ((traits == UndefinedPixelTrait) ||
1904 (reconstruct_traits == UndefinedPixelTrait) ||
1905 ((reconstruct_traits & UpdatePixelTrait) == 0))
1906 continue;
1907 distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
1908 channel,q));
1909 if (distance >= MagickEpsilon)
1910 break;
1911 }
1912 if (i < (ssize_t) GetPixelChannels(image))
1913 break;
1914 p+=(ptrdiff_t) GetPixelChannels(image);
1915 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
1916 }
1917 if (x < (ssize_t) columns)
1918 break;
1919 }
1920 reconstruct_view=DestroyCacheView(reconstruct_view);
1921 image_view=DestroyCacheView(image_view);
1922 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1923}
1924
1925/*
1926%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1927% %
1928% %
1929% %
1930% S e t I m a g e C o l o r M e t r i c %
1931% %
1932% %
1933% %
1934%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1935%
1936% SetImageColorMetric() measures the difference between colors at each pixel
1937% location of two images. A value other than 0 means the colors match
1938% exactly. Otherwise an error measure is computed by summing over all
1939% pixels in an image the distance squared in RGB space between each image
1940% pixel and its corresponding pixel in the reconstruction image. The error
1941% measure is assigned to these image members:
1942%
1943% o mean_error_per_pixel: The mean error for any single pixel in
1944% the image.
1945%
1946% o normalized_mean_error: The normalized mean quantization error for
1947% any single pixel in the image. This distance measure is normalized to
1948% a range between 0 and 1. It is independent of the range of red, green,
1949% and blue values in the image.
1950%
1951% o normalized_maximum_error: The normalized maximum quantization
1952% error for any single pixel in the image. This distance measure is
1953% normalized to a range between 0 and 1. It is independent of the range
1954% of red, green, and blue values in your image.
1955%
1956% A small normalized mean square error, accessed as
1957% image->normalized_mean_error, suggests the images are very similar in
1958% spatial layout and color.
1959%
1960% The format of the SetImageColorMetric method is:
1961%
1962% MagickBooleanType SetImageColorMetric(Image *image,
1963% const Image *reconstruct_image,ExceptionInfo *exception)
1964%
1965% A description of each parameter follows.
1966%
1967% o image: the image.
1968%
1969% o reconstruct_image: the reconstruction image.
1970%
1971% o exception: return any errors or warnings in this structure.
1972%
1973*/
1974MagickExport MagickBooleanType SetImageColorMetric(Image *image,
1975 const Image *reconstruct_image,ExceptionInfo *exception)
1976{
1977 CacheView
1978 *image_view,
1979 *reconstruct_view;
1980
1981 double
1982 area,
1983 maximum_error,
1984 mean_error,
1985 mean_error_per_pixel;
1986
1987 MagickBooleanType
1988 status;
1989
1990 size_t
1991 columns,
1992 rows;
1993
1994 ssize_t
1995 y;
1996
1997 assert(image != (Image *) NULL);
1998 assert(image->signature == MagickCoreSignature);
1999 assert(reconstruct_image != (const Image *) NULL);
2000 assert(reconstruct_image->signature == MagickCoreSignature);
2001 area=0.0;
2002 maximum_error=0.0;
2003 mean_error_per_pixel=0.0;
2004 mean_error=0.0;
2005 rows=MagickMax(image->rows,reconstruct_image->rows);
2006 columns=MagickMax(image->columns,reconstruct_image->columns);
2007 image_view=AcquireVirtualCacheView(image,exception);
2008 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2009 for (y=0; y < (ssize_t) rows; y++)
2010 {
2011 const Quantum
2012 *magick_restrict p,
2013 *magick_restrict q;
2014
2015 ssize_t
2016 x;
2017
2018 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2019 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2020 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2021 break;
2022 for (x=0; x < (ssize_t) columns; x++)
2023 {
2024 ssize_t
2025 i;
2026
2027 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2028 {
2029 double
2030 distance;
2031
2032 PixelChannel channel = GetPixelChannelChannel(image,i);
2033 PixelTrait traits = GetPixelChannelTraits(image,channel);
2034 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2035 channel);
2036 if ((traits == UndefinedPixelTrait) ||
2037 (reconstruct_traits == UndefinedPixelTrait) ||
2038 ((reconstruct_traits & UpdatePixelTrait) == 0))
2039 continue;
2040 distance=fabs((double) p[i]-(double) GetPixelChannel(reconstruct_image,
2041 channel,q));
2042 if (distance >= MagickEpsilon)
2043 {
2044 mean_error_per_pixel+=distance;
2045 mean_error+=distance*distance;
2046 if (distance > maximum_error)
2047 maximum_error=distance;
2048 }
2049 area++;
2050 }
2051 p+=(ptrdiff_t) GetPixelChannels(image);
2052 q+=(ptrdiff_t) GetPixelChannels(reconstruct_image);
2053 }
2054 }
2055 reconstruct_view=DestroyCacheView(reconstruct_view);
2056 image_view=DestroyCacheView(image_view);
2057 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2058 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
2059 mean_error/area);
2060 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2061 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2062 return(status);
2063}
2064
2065/*
2066%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2067% %
2068% %
2069% %
2070% S i m i l a r i t y I m a g e %
2071% %
2072% %
2073% %
2074%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2075%
2076% SimilarityImage() compares the reconstruction of the image and returns the
2077% best match offset. In addition, it returns a similarity image such that an
2078% exact match location is completely white and if none of the pixels match,
2079% black, otherwise some gray level in-between.
2080%
2081% Contributed by Fred Weinhaus.
2082%
2083% The format of the SimilarityImageImage method is:
2084%
2085% Image *SimilarityImage(const Image *image,const Image *reconstruct,
2086% const MetricType metric,const double similarity_threshold,
2087% RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2088%
2089% A description of each parameter follows:
2090%
2091% o image: the image.
2092%
2093% o reconstruct: find an area of the image that closely resembles this image.
2094%
2095% o metric: the metric.
2096%
2097% o similarity_threshold: minimum distortion for (sub)image match.
2098%
2099% o offset: the best match offset of the reconstruction image within the
2100% image.
2101%
2102% o similarity: the computed similarity between the images.
2103%
2104% o exception: return any errors or warnings in this structure.
2105%
2106*/
2107
2108#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
2109static Image *SIMCrossCorrelationImage(const Image *alpha_image,
2110 const Image *beta_image,ExceptionInfo *exception)
2111{
2112 Image
2113 *alpha_fft = (Image *) NULL,
2114 *beta_fft = (Image *) NULL,
2115 *complex_conjugate = (Image *) NULL,
2116 *complex_multiplication = (Image *) NULL,
2117 *cross_correlation = (Image *) NULL,
2118 *temp_image = (Image *) NULL;
2119
2120 /*
2121 Take the FFT of beta (reconstruction) image.
2122 */
2123 temp_image=CloneImage(beta_image,0,0,MagickTrue,exception);
2124 if (temp_image == (Image *) NULL)
2125 return((Image *) NULL);
2126 (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2127 beta_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2128 temp_image=DestroyImageList(temp_image);
2129 if (beta_fft == (Image *) NULL)
2130 return((Image *) NULL);
2131 /*
2132 Take the complex conjugate of beta_fft.
2133 */
2134 complex_conjugate=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
2135 beta_fft=DestroyImageList(beta_fft);
2136 if (complex_conjugate == (Image *) NULL)
2137 return((Image *) NULL);
2138 /*
2139 Take the FFT of the alpha (test) image.
2140 */
2141 temp_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2142 if (temp_image == (Image *) NULL)
2143 {
2144 complex_conjugate=DestroyImageList(complex_conjugate);
2145 return((Image *) NULL);
2146 }
2147 (void) SetImageArtifact(temp_image,"fourier:normalize","inverse");
2148 alpha_fft=ForwardFourierTransformImage(temp_image,MagickFalse,exception);
2149 temp_image=DestroyImageList(temp_image);
2150 if (alpha_fft == (Image *) NULL)
2151 {
2152 complex_conjugate=DestroyImageList(complex_conjugate);
2153 return((Image *) NULL);
2154 }
2155 /*
2156 Do complex multiplication.
2157 */
2158 DisableCompositeClampUnlessSpecified(complex_conjugate);
2159 DisableCompositeClampUnlessSpecified(complex_conjugate->next);
2160 complex_conjugate->next->next=alpha_fft;
2161 complex_multiplication=ComplexImages(complex_conjugate,
2162 MultiplyComplexOperator,exception);
2163 complex_conjugate=DestroyImageList(complex_conjugate);
2164 if (complex_multiplication == (Image *) NULL)
2165 return((Image *) NULL);
2166 /*
2167 Do the IFT and return the cross-correlation result.
2168 */
2169 cross_correlation=InverseFourierTransformImage(complex_multiplication,
2170 complex_multiplication->next,MagickFalse,exception);
2171 complex_multiplication=DestroyImageList(complex_multiplication);
2172 return(cross_correlation);
2173}
2174
2175static Image *SIMDerivativeImage(const Image *image,const char *kernel,
2176 ExceptionInfo *exception)
2177{
2178 Image
2179 *derivative_image;
2180
2182 *kernel_info;
2183
2184 kernel_info=AcquireKernelInfo(kernel,exception);
2185 if (kernel_info == (KernelInfo *) NULL)
2186 return((Image *) NULL);
2187 derivative_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,
2188 exception);
2189 kernel_info=DestroyKernelInfo(kernel_info);
2190 return(derivative_image);
2191}
2192
2193static Image *SIMDivideImage(const Image *numerator_image,
2194 const Image *denominator_image,ExceptionInfo *exception)
2195{
2196 CacheView
2197 *denominator_view,
2198 *numerator_view;
2199
2200 Image
2201 *divide_image;
2202
2203 MagickBooleanType
2204 status;
2205
2206 ssize_t
2207 y;
2208
2209 /*
2210 Divide one image into another.
2211 */
2212 divide_image=CloneImage(numerator_image,0,0,MagickTrue,exception);
2213 if (divide_image == (Image *) NULL)
2214 return(divide_image);
2215 status=MagickTrue;
2216 numerator_view=AcquireAuthenticCacheView(divide_image,exception);
2217 denominator_view=AcquireVirtualCacheView(denominator_image,exception);
2218#if defined(MAGICKCORE_OPENMP_SUPPORT)
2219 #pragma omp parallel for schedule(static) shared(status) \
2220 magick_number_threads(denominator_image,divide_image,divide_image->rows,1)
2221#endif
2222 for (y=0; y < (ssize_t) divide_image->rows; y++)
2223 {
2224 const Quantum
2225 *magick_restrict p;
2226
2227 Quantum
2228 *magick_restrict q;
2229
2230 ssize_t
2231 x;
2232
2233 if (status == MagickFalse)
2234 continue;
2235 p=GetCacheViewVirtualPixels(denominator_view,0,y,
2236 denominator_image->columns,1,exception);
2237 q=GetCacheViewAuthenticPixels(numerator_view,0,y,divide_image->columns,1,
2238 exception);
2239 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2240 {
2241 status=MagickFalse;
2242 continue;
2243 }
2244 for (x=0; x < (ssize_t) divide_image->columns; x++)
2245 {
2246 ssize_t
2247 i;
2248
2249 for (i=0; i < (ssize_t) GetPixelChannels(divide_image); i++)
2250 {
2251 PixelChannel channel = GetPixelChannelChannel(divide_image,i);
2252 PixelTrait traits = GetPixelChannelTraits(divide_image,channel);
2253 if ((traits & UpdatePixelTrait) == 0)
2254 continue;
2255 if (fabs(p[i]) >= MagickEpsilon)
2256 q[i]=(Quantum) ((double) q[i]*PerceptibleReciprocal(QuantumScale*
2257 (double) p[i]));
2258 }
2259 p+=(ptrdiff_t) GetPixelChannels(denominator_image);
2260 q+=(ptrdiff_t) GetPixelChannels(divide_image);
2261 }
2262 if (SyncCacheViewAuthenticPixels(numerator_view,exception) == MagickFalse)
2263 status=MagickFalse;
2264 }
2265 denominator_view=DestroyCacheView(denominator_view);
2266 numerator_view=DestroyCacheView(numerator_view);
2267 if (status == MagickFalse)
2268 divide_image=DestroyImage(divide_image);
2269 return(divide_image);
2270}
2271
2272static Image *SIMDivideByMagnitude(Image *image,Image *magnitude_image,
2273 const Image *source_image,ExceptionInfo *exception)
2274{
2275 Image
2276 *divide_image,
2277 *result_image;
2278
2280 geometry;
2281
2282 divide_image=SIMDivideImage(image,magnitude_image,exception);
2283 if (divide_image == (Image *) NULL)
2284 return((Image *) NULL);
2285 GetPixelInfoRGBA(0,0,0,0,&divide_image->background_color);
2286 SetGeometry(source_image,&geometry);
2287 geometry.width=MagickMax(source_image->columns,divide_image->columns);
2288 geometry.height=MagickMax(source_image->rows,divide_image->rows);
2289 result_image=ExtentImage(divide_image,&geometry,exception);
2290 divide_image=DestroyImage(divide_image);
2291 return(result_image);
2292}
2293
2294static MagickBooleanType SIMLogImage(Image *image,ExceptionInfo *exception)
2295{
2296 CacheView
2297 *image_view;
2298
2299 MagickBooleanType
2300 status;
2301
2302 ssize_t
2303 y;
2304
2305 /*
2306 Take the log of each pixel.
2307 */
2308 status=MagickTrue;
2309 image_view=AcquireAuthenticCacheView(image,exception);
2310#if defined(MAGICKCORE_OPENMP_SUPPORT)
2311 #pragma omp parallel for schedule(static) shared(status) \
2312 magick_number_threads(image,image,image->rows,1)
2313#endif
2314 for (y=0; y < (ssize_t) image->rows; y++)
2315 {
2316 Quantum
2317 *magick_restrict q;
2318
2319 ssize_t
2320 x;
2321
2322 if (status == MagickFalse)
2323 continue;
2324 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2325 if (q == (Quantum *) NULL)
2326 {
2327 status=MagickFalse;
2328 continue;
2329 }
2330 for (x=0; x < (ssize_t) image->columns; x++)
2331 {
2332 ssize_t
2333 i;
2334
2335 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2336 {
2337 PixelChannel channel = GetPixelChannelChannel(image,i);
2338 PixelTrait traits = GetPixelChannelTraits(image,channel);
2339 if ((traits & UpdatePixelTrait) == 0)
2340 continue;
2341 q[i]=(Quantum) (QuantumRange*10.0*MagickLog10((double) q[i]));
2342 }
2343 q+=(ptrdiff_t) GetPixelChannels(image);
2344 }
2345 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2346 status=MagickFalse;
2347 }
2348 image_view=DestroyCacheView(image_view);
2349 return(status);
2350}
2351
2352static Image *SIMSquareImage(const Image *image,ExceptionInfo *exception)
2353{
2354 CacheView
2355 *image_view;
2356
2357 Image
2358 *square_image;
2359
2360 MagickBooleanType
2361 status;
2362
2363 ssize_t
2364 y;
2365
2366 /*
2367 Square each pixel in the image.
2368 */
2369 square_image=CloneImage(image,0,0,MagickTrue,exception);
2370 if (square_image == (Image *) NULL)
2371 return(square_image);
2372 status=MagickTrue;
2373 image_view=AcquireAuthenticCacheView(square_image,exception);
2374#if defined(MAGICKCORE_OPENMP_SUPPORT)
2375 #pragma omp parallel for schedule(static) shared(status) \
2376 magick_number_threads(square_image,square_image,square_image->rows,1)
2377#endif
2378 for (y=0; y < (ssize_t) square_image->rows; y++)
2379 {
2380 Quantum
2381 *magick_restrict q;
2382
2383 ssize_t
2384 x;
2385
2386 if (status == MagickFalse)
2387 continue;
2388 q=GetCacheViewAuthenticPixels(image_view,0,y,square_image->columns,1,
2389 exception);
2390 if (q == (Quantum *) NULL)
2391 {
2392 status=MagickFalse;
2393 continue;
2394 }
2395 for (x=0; x < (ssize_t) square_image->columns; x++)
2396 {
2397 ssize_t
2398 i;
2399
2400 for (i=0; i < (ssize_t) GetPixelChannels(square_image); i++)
2401 {
2402 PixelChannel channel = GetPixelChannelChannel(square_image,i);
2403 PixelTrait traits = GetPixelChannelTraits(square_image,channel);
2404 if ((traits & UpdatePixelTrait) == 0)
2405 continue;
2406 q[i]=(Quantum) ((double) q[i]*QuantumScale*(double) q[i]);
2407 }
2408 q+=(ptrdiff_t) GetPixelChannels(square_image);
2409 }
2410 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2411 status=MagickFalse;
2412 }
2413 image_view=DestroyCacheView(image_view);
2414 if (status == MagickFalse)
2415 square_image=DestroyImage(square_image);
2416 return(square_image);
2417}
2418
2419static Image *SIMMagnitudeImage(Image *alpha_image,Image *beta_image,
2420 ExceptionInfo *exception)
2421{
2422 Image
2423 *magnitude_image,
2424 *xsq_image,
2425 *ysq_image;
2426
2427 MagickBooleanType
2428 status;
2429
2430 (void) SetImageArtifact(alpha_image,"compose:clamp","False");
2431 xsq_image=SIMSquareImage(alpha_image,exception);
2432 if (xsq_image == (Image *) NULL)
2433 return((Image *) NULL);
2434 (void) SetImageArtifact(beta_image,"compose:clamp","False");
2435 ysq_image=SIMSquareImage(beta_image,exception);
2436 if (ysq_image == (Image *) NULL)
2437 {
2438 xsq_image=DestroyImage(xsq_image);
2439 return((Image *) NULL);
2440 }
2441 status=CompositeImage(xsq_image,ysq_image,PlusCompositeOp,MagickTrue,0,0,
2442 exception);
2443 magnitude_image=xsq_image;
2444 ysq_image=DestroyImage(ysq_image);
2445 if (status == MagickFalse)
2446 {
2447 magnitude_image=DestroyImage(magnitude_image);
2448 return((Image *) NULL);
2449 }
2450 status=EvaluateImage(magnitude_image,PowEvaluateOperator,0.5,exception);
2451 if (status == MagickFalse)
2452 {
2453 magnitude_image=DestroyImage(magnitude_image);
2454 return (Image *) NULL;
2455 }
2456 return(magnitude_image);
2457}
2458
2459static MagickBooleanType SIMMaximaImage(const Image *image,double *maxima,
2460 RectangleInfo *offset,ExceptionInfo *exception)
2461{
2462 CacheView
2463 *image_view;
2464
2465 MagickBooleanType
2466 status;
2467
2468 ssize_t
2469 y;
2470
2471 /*
2472 Identify the maxima value in the image and its location.
2473 */
2474 status=MagickTrue;
2475 *maxima=MagickMinimumValue;
2476 offset->x=0;
2477 offset->y=0;
2478 image_view=AcquireVirtualCacheView(image,exception);
2479#if defined(MAGICKCORE_OPENMP_SUPPORT)
2480 #pragma omp parallel for schedule(static) shared(status) \
2481 magick_number_threads(image,image,image->rows,1)
2482#endif
2483 for (y=0; y < (ssize_t) image->rows; y++)
2484 {
2485 const Quantum
2486 *magick_restrict p;
2487
2488 double
2489 row_maxima;
2490
2491 ssize_t
2492 row_x,
2493 x;
2494
2495 if (status == MagickFalse)
2496 continue;
2497 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2498 if (p == (const Quantum *) NULL)
2499 {
2500 status=MagickFalse;
2501 continue;
2502 }
2503 row_maxima=(double) p[0];
2504 row_x=0;
2505 for (x=0; x < (ssize_t) image->columns; x++)
2506 {
2507 ssize_t
2508 i;
2509
2510 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2511 {
2512 PixelChannel channel = GetPixelChannelChannel(image,i);
2513 PixelTrait traits = GetPixelChannelTraits(image,channel);
2514 if ((traits & UpdatePixelTrait) == 0)
2515 continue;
2516 if ((double) p[i] > row_maxima)
2517 {
2518 row_maxima=(double) p[i];
2519 row_x=x;
2520 }
2521 }
2522 p+=(ptrdiff_t) GetPixelChannels(image);
2523 }
2524#if defined(MAGICKCORE_OPENMP_SUPPORT)
2525 #pragma omp critical (MagickCore_SIMMaximaImage)
2526#endif
2527 if (row_maxima > *maxima)
2528 {
2529 *maxima=row_maxima;
2530 offset->x=row_x;
2531 offset->y=y;
2532 }
2533 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2534 status=MagickFalse;
2535 }
2536 image_view=DestroyCacheView(image_view);
2537 return(status);
2538}
2539
2540static MagickBooleanType SIMMinimaImage(const Image *image,double *minima,
2541 RectangleInfo *offset,ExceptionInfo *exception)
2542{
2543 CacheView
2544 *image_view;
2545
2546 MagickBooleanType
2547 status;
2548
2549 ssize_t
2550 y;
2551
2552 /*
2553 Identify the minima value in the image and its location.
2554 */
2555 status=MagickTrue;
2556 *minima=MagickMaximumValue;
2557 offset->x=0;
2558 offset->y=0;
2559 image_view=AcquireVirtualCacheView(image,exception);
2560#if defined(MAGICKCORE_OPENMP_SUPPORT)
2561 #pragma omp parallel for schedule(static) shared(status) \
2562 magick_number_threads(image,image,image->rows,1)
2563#endif
2564 for (y=0; y < (ssize_t) image->rows; y++)
2565 {
2566 const Quantum
2567 *magick_restrict p;
2568
2569 double
2570 row_minima;
2571
2572 ssize_t
2573 row_x,
2574 x;
2575
2576 if (status == MagickFalse)
2577 continue;
2578 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2579 if (p == (const Quantum *) NULL)
2580 {
2581 status=MagickFalse;
2582 continue;
2583 }
2584 row_minima=(double) p[0];
2585 row_x=0;
2586 for (x=0; x < (ssize_t) image->columns; x++)
2587 {
2588 ssize_t
2589 i;
2590
2591 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2592 {
2593 PixelChannel channel = GetPixelChannelChannel(image,i);
2594 PixelTrait traits = GetPixelChannelTraits(image,channel);
2595 if ((traits & UpdatePixelTrait) == 0)
2596 continue;
2597 if ((double) p[i] < row_minima)
2598 {
2599 row_minima=(double) p[i];
2600 row_x=x;
2601 }
2602 }
2603 p+=(ptrdiff_t) GetPixelChannels(image);
2604 }
2605#if defined(MAGICKCORE_OPENMP_SUPPORT)
2606 #pragma omp critical (MagickCore_SIMMinimaImage)
2607#endif
2608 if (row_minima < *minima)
2609 {
2610 *minima=row_minima;
2611 offset->x=row_x;
2612 offset->y=y;
2613 }
2614 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2615 status=MagickFalse;
2616 }
2617 image_view=DestroyCacheView(image_view);
2618 return(status);
2619}
2620
2621static MagickBooleanType SIMMultiplyImage(Image *image,const double factor,
2622 const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
2623{
2624 CacheView
2625 *image_view;
2626
2627 MagickBooleanType
2628 status;
2629
2630 ssize_t
2631 y;
2632
2633 /*
2634 Multiply each pixel by a factor.
2635 */
2636 status=MagickTrue;
2637 image_view=AcquireAuthenticCacheView(image,exception);
2638#if defined(MAGICKCORE_OPENMP_SUPPORT)
2639 #pragma omp parallel for schedule(static) shared(status) \
2640 magick_number_threads(image,image,image->rows,1)
2641#endif
2642 for (y=0; y < (ssize_t) image->rows; y++)
2643 {
2644 Quantum
2645 *magick_restrict q;
2646
2647 ssize_t
2648 x;
2649
2650 if (status == MagickFalse)
2651 continue;
2652 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2653 if (q == (Quantum *) NULL)
2654 {
2655 status=MagickFalse;
2656 continue;
2657 }
2658 for (x=0; x < (ssize_t) image->columns; x++)
2659 {
2660 ssize_t
2661 i;
2662
2663 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2664 {
2665 PixelChannel channel = GetPixelChannelChannel(image,i);
2666 PixelTrait traits = GetPixelChannelTraits(image,channel);
2667 if ((traits & UpdatePixelTrait) == 0)
2668 continue;
2669 if (channel_statistics != (const ChannelStatistics *) NULL)
2670 q[i]=(Quantum) ((double) q[i]*QuantumScale*
2671 channel_statistics[channel].standard_deviation);
2672 q[i]=(Quantum) ((double) q[i]*factor);
2673 }
2674 q+=(ptrdiff_t) GetPixelChannels(image);
2675 }
2676 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2677 status=MagickFalse;
2678 }
2679 image_view=DestroyCacheView(image_view);
2680 return(status);
2681}
2682
2683static Image *SIMPhaseCorrelationImage(const Image *alpha_image,
2684 const Image *beta_image,const Image *magnitude_image,ExceptionInfo *exception)
2685{
2686 Image
2687 *alpha_fft = (Image *) NULL,
2688 *beta_fft = (Image *) NULL,
2689 *complex_multiplication = (Image *) NULL,
2690 *cross_correlation = (Image *) NULL;
2691
2692 /*
2693 Take the FFT of the beta (reconstruction) image.
2694 */
2695 beta_fft=CloneImage(beta_image,0,0,MagickTrue,exception);
2696 if (beta_fft == NULL)
2697 return((Image *) NULL);
2698 (void) SetImageArtifact(beta_fft,"fourier:normalize","inverse");
2699 beta_fft=ForwardFourierTransformImage(beta_fft,MagickFalse,exception);
2700 if (beta_fft == NULL)
2701 return((Image *) NULL);
2702 /*
2703 Take the FFT of the alpha (test) image.
2704 */
2705 alpha_fft=CloneImage(alpha_image,0,0,MagickTrue,exception);
2706 if (alpha_fft == (Image *) NULL)
2707 {
2708 beta_fft=DestroyImageList(beta_fft);
2709 return((Image *) NULL);
2710 }
2711 (void) SetImageArtifact(alpha_fft,"fourier:normalize","inverse");
2712 alpha_fft=ForwardFourierTransformImage(alpha_fft,MagickFalse,exception);
2713 if (alpha_fft == (Image *) NULL)
2714 {
2715 beta_fft=DestroyImageList(beta_fft);
2716 return((Image *) NULL);
2717 }
2718 /*
2719 Take the complex conjugate of the beta FFT.
2720 */
2721 beta_fft=ComplexImages(beta_fft,ConjugateComplexOperator,exception);
2722 if (beta_fft == (Image *) NULL)
2723 {
2724 alpha_fft=DestroyImageList(alpha_fft);
2725 return((Image *) NULL);
2726 }
2727 /*
2728 Do complex multiplication.
2729 */
2730 beta_fft->next->next=alpha_fft;
2731 DisableCompositeClampUnlessSpecified(beta_fft);
2732 DisableCompositeClampUnlessSpecified(beta_fft->next);
2733 complex_multiplication=ComplexImages(beta_fft,MultiplyComplexOperator,
2734 exception);
2735 beta_fft=DestroyImageList(beta_fft);
2736 if (complex_multiplication == (Image *) NULL)
2737 return((Image *) NULL);
2738 /*
2739 Divide the results.
2740 */
2741 CompositeLayers(complex_multiplication,DivideSrcCompositeOp,(Image *)
2742 magnitude_image,0,0,exception);
2743 /*
2744 Do the IFT and return the cross-correlation result.
2745 */
2746 (void) SetImageArtifact(complex_multiplication,"fourier:normalize","inverse");
2747 cross_correlation=InverseFourierTransformImage(complex_multiplication,
2748 complex_multiplication->next,MagickFalse,exception);
2749 complex_multiplication=DestroyImageList(complex_multiplication);
2750 return(cross_correlation);
2751}
2752
2753static MagickBooleanType SIMSetImageMean(const Image *image,
2754 const ChannelStatistics *channel_statistics,ExceptionInfo *exception)
2755{
2756 CacheView
2757 *image_view;
2758
2759 MagickBooleanType
2760 status;
2761
2762 ssize_t
2763 y;
2764
2765 /*
2766 Set image mean.
2767 */
2768 status=MagickTrue;
2769 image_view=AcquireAuthenticCacheView(image,exception);
2770#if defined(MAGICKCORE_OPENMP_SUPPORT)
2771 #pragma omp parallel for schedule(static) shared(status) \
2772 magick_number_threads(image,image,image->rows,1)
2773#endif
2774 for (y=0; y < (ssize_t) image->rows; y++)
2775 {
2776 Quantum
2777 *magick_restrict q;
2778
2779 ssize_t
2780 x;
2781
2782 if (status == MagickFalse)
2783 continue;
2784 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2785 if (q == (Quantum *) NULL)
2786 {
2787 status=MagickFalse;
2788 continue;
2789 }
2790 for (x=0; x < (ssize_t) image->columns; x++)
2791 {
2792 ssize_t
2793 i;
2794
2795 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2796 {
2797 PixelChannel channel = GetPixelChannelChannel(image,i);
2798 PixelTrait traits = GetPixelChannelTraits(image,channel);
2799 if ((traits & UpdatePixelTrait) == 0)
2800 continue;
2801 q[i]=(Quantum) channel_statistics[channel].mean;
2802 }
2803 q+=(ptrdiff_t) GetPixelChannels(image);
2804 }
2805 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2806 status=MagickFalse;
2807 }
2808 image_view=DestroyCacheView(image_view);
2809 return(status);
2810}
2811
2812static Image *SIMSubtractImageMean(const Image *alpha_image,
2813 const Image *beta_image,const ChannelStatistics *channel_statistics,
2814 ExceptionInfo *exception)
2815{
2816 CacheView
2817 *beta_view,
2818 *image_view;
2819
2820 Image
2821 *subtract_image;
2822
2823 MagickBooleanType
2824 status;
2825
2826 ssize_t
2827 y;
2828
2829 /*
2830 Subtract the image mean and pad.
2831 */
2832 subtract_image=CloneImage(beta_image,alpha_image->columns,alpha_image->rows,
2833 MagickTrue,exception);
2834 if (subtract_image == (Image *) NULL)
2835 return(subtract_image);
2836 status=MagickTrue;
2837 image_view=AcquireAuthenticCacheView(subtract_image,exception);
2838 beta_view=AcquireVirtualCacheView(beta_image,exception);
2839#if defined(MAGICKCORE_OPENMP_SUPPORT)
2840 #pragma omp parallel for schedule(static) shared(status) \
2841 magick_number_threads(beta_image,subtract_image,subtract_image->rows,1)
2842#endif
2843 for (y=0; y < (ssize_t) subtract_image->rows; y++)
2844 {
2845 const Quantum
2846 *magick_restrict p;
2847
2848 Quantum
2849 *magick_restrict q;
2850
2851 ssize_t
2852 x;
2853
2854 if (status == MagickFalse)
2855 continue;
2856 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,exception);
2857 q=GetCacheViewAuthenticPixels(image_view,0,y,subtract_image->columns,1,
2858 exception);
2859 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2860 {
2861 status=MagickFalse;
2862 continue;
2863 }
2864 for (x=0; x < (ssize_t) subtract_image->columns; x++)
2865 {
2866 ssize_t
2867 i;
2868
2869 for (i=0; i < (ssize_t) GetPixelChannels(subtract_image); i++)
2870 {
2871 PixelChannel channel = GetPixelChannelChannel(subtract_image,i);
2872 PixelTrait traits = GetPixelChannelTraits(subtract_image,channel);
2873 if ((traits & UpdatePixelTrait) == 0)
2874 continue;
2875 if ((x >= (ssize_t) beta_image->columns) ||
2876 (y >= (ssize_t) beta_image->rows))
2877 q[i]=(Quantum) 0;
2878 else
2879 q[i]=(Quantum) ((double) p[i]-channel_statistics[channel].mean);
2880 }
2881 p+=(ptrdiff_t) GetPixelChannels(beta_image);
2882 q+=(ptrdiff_t) GetPixelChannels(subtract_image);
2883 }
2884 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2885 status=MagickFalse;
2886 }
2887 beta_view=DestroyCacheView(beta_view);
2888 image_view=DestroyCacheView(image_view);
2889 if (status == MagickFalse)
2890 subtract_image=DestroyImage(subtract_image);
2891 return(subtract_image);
2892}
2893
2894static Image *SIMUnityImage(const Image *alpha_image,const Image *beta_image,
2895 ExceptionInfo *exception)
2896{
2897 CacheView
2898 *image_view;
2899
2900 Image
2901 *unity_image;
2902
2903 MagickBooleanType
2904 status;
2905
2906 ssize_t
2907 y;
2908
2909 /*
2910 Create a padded unity image.
2911 */
2912 unity_image=CloneImage(alpha_image,alpha_image->columns,alpha_image->rows,
2913 MagickTrue,exception);
2914 if (unity_image == (Image *) NULL)
2915 return(unity_image);
2916 if (SetImageStorageClass(unity_image,DirectClass,exception) == MagickFalse)
2917 return(DestroyImage(unity_image));
2918 status=MagickTrue;
2919 image_view=AcquireAuthenticCacheView(unity_image,exception);
2920#if defined(MAGICKCORE_OPENMP_SUPPORT)
2921 #pragma omp parallel for schedule(static) shared(status) \
2922 magick_number_threads(unity_image,unity_image,unity_image->rows,1)
2923#endif
2924 for (y=0; y < (ssize_t) unity_image->rows; y++)
2925 {
2926 Quantum
2927 *magick_restrict q;
2928
2929 ssize_t
2930 x;
2931
2932 if (status == MagickFalse)
2933 continue;
2934 q=GetCacheViewAuthenticPixels(image_view,0,y,unity_image->columns,1,
2935 exception);
2936 if (q == (Quantum *) NULL)
2937 {
2938 status=MagickFalse;
2939 continue;
2940 }
2941 for (x=0; x < (ssize_t) unity_image->columns; x++)
2942 {
2943 ssize_t
2944 i;
2945
2946 for (i=0; i < (ssize_t) GetPixelChannels(unity_image); i++)
2947 {
2948 PixelChannel channel = GetPixelChannelChannel(unity_image,i);
2949 PixelTrait traits = GetPixelChannelTraits(unity_image,channel);
2950 if ((traits & UpdatePixelTrait) == 0)
2951 continue;
2952 if ((x >= (ssize_t) beta_image->columns) ||
2953 (y >= (ssize_t) beta_image->rows))
2954 q[i]=(Quantum) 0;
2955 else
2956 q[i]=QuantumRange;
2957 }
2958 q+=(ptrdiff_t) GetPixelChannels(unity_image);
2959 }
2960 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2961 status=MagickFalse;
2962 }
2963 image_view=DestroyCacheView(image_view);
2964 if (status == MagickFalse)
2965 unity_image=DestroyImage(unity_image);
2966 return(unity_image);
2967}
2968
2969static Image *SIMVarianceImage(Image *alpha_image,const Image *beta_image,
2970 ExceptionInfo *exception)
2971{
2972 CacheView
2973 *beta_view,
2974 *image_view;
2975
2976 Image
2977 *variance_image;
2978
2979 MagickBooleanType
2980 status;
2981
2982 ssize_t
2983 y;
2984
2985 /*
2986 Compute the variance of the two images.
2987 */
2988 variance_image=CloneImage(alpha_image,0,0,MagickTrue,exception);
2989 if (variance_image == (Image *) NULL)
2990 return(variance_image);
2991 status=MagickTrue;
2992 image_view=AcquireAuthenticCacheView(variance_image,exception);
2993 beta_view=AcquireVirtualCacheView(beta_image,exception);
2994#if defined(MAGICKCORE_OPENMP_SUPPORT)
2995 #pragma omp parallel for schedule(static) shared(status) \
2996 magick_number_threads(beta_image,variance_image,variance_image->rows,1)
2997#endif
2998 for (y=0; y < (ssize_t) variance_image->rows; y++)
2999 {
3000 const Quantum
3001 *magick_restrict p;
3002
3003 Quantum
3004 *magick_restrict q;
3005
3006 ssize_t
3007 x;
3008
3009 if (status == MagickFalse)
3010 continue;
3011 p=GetCacheViewVirtualPixels(beta_view,0,y,beta_image->columns,1,
3012 exception);
3013 q=GetCacheViewAuthenticPixels(image_view,0,y,variance_image->columns,1,
3014 exception);
3015 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3016 {
3017 status=MagickFalse;
3018 continue;
3019 }
3020 for (x=0; x < (ssize_t) variance_image->columns; x++)
3021 {
3022 ssize_t
3023 i;
3024
3025 for (i=0; i < (ssize_t) GetPixelChannels(variance_image); i++)
3026 {
3027 PixelChannel channel = GetPixelChannelChannel(variance_image,i);
3028 PixelTrait traits = GetPixelChannelTraits(variance_image,channel);
3029 if ((traits & UpdatePixelTrait) == 0)
3030 continue;
3031 q[i]=(Quantum) ((double) ClampToQuantum((double) QuantumRange*
3032 sqrt(fabs(QuantumScale*((double) q[i]-(double) p[i]))))/
3033 sqrt((double) QuantumRange));
3034 }
3035 p+=(ptrdiff_t) GetPixelChannels(beta_image);
3036 q+=(ptrdiff_t) GetPixelChannels(variance_image);
3037 }
3038 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3039 status=MagickFalse;
3040 }
3041 beta_view=DestroyCacheView(beta_view);
3042 image_view=DestroyCacheView(image_view);
3043 if (status == MagickFalse)
3044 variance_image=DestroyImage(variance_image);
3045 return(variance_image);
3046}
3047
3048static Image *DPCSimilarityImage(const Image *image,const Image *reconstruct,
3049 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3050{
3051#define ThrowDPCSimilarityException() \
3052{ \
3053 if (dot_product_image != (Image *) NULL) \
3054 dot_product_image=DestroyImage(dot_product_image); \
3055 if (magnitude_image != (Image *) NULL) \
3056 magnitude_image=DestroyImage(magnitude_image); \
3057 if (reconstruct_image != (Image *) NULL) \
3058 reconstruct_image=DestroyImage(reconstruct_image); \
3059 if (rx_image != (Image *) NULL) \
3060 rx_image=DestroyImage(rx_image); \
3061 if (ry_image != (Image *) NULL) \
3062 ry_image=DestroyImage(ry_image); \
3063 if (test_image != (Image *) NULL) \
3064 test_image=DestroyImage(test_image); \
3065 if (threshold_image != (Image *) NULL) \
3066 threshold_image=DestroyImage(threshold_image); \
3067 if (trx_image != (Image *) NULL) \
3068 trx_image=DestroyImage(trx_image); \
3069 if (try_image != (Image *) NULL) \
3070 try_image=DestroyImage(try_image); \
3071 if (tx_image != (Image *) NULL) \
3072 tx_image=DestroyImage(tx_image); \
3073 if (ty_image != (Image *) NULL) \
3074 ty_image=DestroyImage(ty_image); \
3075 return((Image *) NULL); \
3076}
3077
3078 double
3079 edge_factor = 0.0,
3080 maxima = 0.0,
3081 mean = 0.0,
3082 standard_deviation = 0.0;
3083
3084 Image
3085 *dot_product_image = (Image *) NULL,
3086 *magnitude_image = (Image *) NULL,
3087 *reconstruct_image = (Image *) NULL,
3088 *rx_image = (Image *) NULL,
3089 *ry_image = (Image *) NULL,
3090 *trx_image = (Image *) NULL,
3091 *temp_image = (Image *) NULL,
3092 *test_image = (Image *) NULL,
3093 *threshold_image = (Image *) NULL,
3094 *try_image = (Image *) NULL,
3095 *tx_image = (Image *) NULL,
3096 *ty_image = (Image *) NULL;
3097
3098 MagickBooleanType
3099 status;
3100
3102 geometry;
3103
3104 /*
3105 Dot product correlation-based image similarity using FFT local statistics.
3106 */
3107 test_image=CloneImage(image,0,0,MagickTrue,exception);
3108 if (test_image == (Image *) NULL)
3109 return((Image *) NULL);
3110 GetPixelInfoRGBA(0,0,0,0,&test_image->background_color);
3111 (void) ResetImagePage(test_image,"0x0+0+0");
3112 status=SetImageExtent(test_image,2*(size_t) ceil((double) image->columns/2.0),
3113 2*(size_t) ceil((double) image->rows/2.0),exception);
3114 if (status == MagickFalse)
3115 ThrowDPCSimilarityException();
3116 (void) SetImageAlphaChannel(test_image,OffAlphaChannel,exception);
3117 /*
3118 Compute the cross correlation of the test and reconstruct magnitudes.
3119 */
3120 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
3121 if (reconstruct_image == (Image *) NULL)
3122 ThrowDPCSimilarityException();
3123 GetPixelInfoRGBA(0,0,0,0,&reconstruct_image->background_color);
3124 (void) ResetImagePage(reconstruct_image,"0x0+0+0");
3125 (void) SetImageAlphaChannel(reconstruct_image,OffAlphaChannel,exception);
3126 /*
3127 Compute X and Y derivatives of reference image.
3128 */
3129 (void) SetImageVirtualPixelMethod(reconstruct_image,EdgeVirtualPixelMethod,
3130 exception);
3131 rx_image=SIMDerivativeImage(reconstruct_image,"Sobel",exception);
3132 if (rx_image == (Image *) NULL)
3133 ThrowDPCSimilarityException();
3134 ry_image=SIMDerivativeImage(reconstruct_image,"Sobel:90",exception);
3135 reconstruct_image=DestroyImage(reconstruct_image);
3136 if (ry_image == (Image *) NULL)
3137 ThrowDPCSimilarityException();
3138 /*
3139 Compute magnitude of derivatives.
3140 */
3141 magnitude_image=SIMMagnitudeImage(rx_image,ry_image,exception);
3142 if (magnitude_image == (Image *) NULL)
3143 ThrowDPCSimilarityException();
3144 /*
3145 Compute an edge normalization correction.
3146 */
3147 threshold_image=CloneImage(magnitude_image,0,0,MagickTrue,exception);
3148 if (threshold_image == (Image *) NULL)
3149 ThrowDPCSimilarityException();
3150 status=BilevelImage(threshold_image,0.0,exception);
3151 if (status == MagickFalse)
3152 ThrowDPCSimilarityException();
3153 status=GetImageMean(threshold_image,&mean,&standard_deviation,exception);
3154 threshold_image=DestroyImage(threshold_image);
3155 if (status == MagickFalse)
3156 ThrowDPCSimilarityException();
3157 edge_factor=1.0/(QuantumScale*mean)/reconstruct->columns/reconstruct->rows;
3158 /*
3159 Divide X and Y derivitives of reference image by magnitude.
3160 */
3161 temp_image=SIMDivideByMagnitude(rx_image,magnitude_image,image,exception);
3162 rx_image=DestroyImage(rx_image);
3163 if (temp_image == (Image *) NULL)
3164 ThrowDPCSimilarityException();
3165 rx_image=temp_image;
3166 try_image=SIMDivideByMagnitude(ry_image,magnitude_image,image,exception);
3167 magnitude_image=DestroyImage(magnitude_image);
3168 ry_image=DestroyImage(ry_image);
3169 if (try_image == (Image *) NULL)
3170 ThrowDPCSimilarityException();
3171 ry_image=try_image;
3172 /*
3173 Compute X and Y derivatives of image.
3174 */
3175 (void) SetImageVirtualPixelMethod(test_image,EdgeVirtualPixelMethod,
3176 exception);
3177 tx_image=SIMDerivativeImage(test_image,"Sobel",exception);
3178 if (tx_image == (Image *) NULL)
3179 ThrowDPCSimilarityException();
3180 ty_image=SIMDerivativeImage(test_image,"Sobel:90",exception);
3181 test_image=DestroyImage(test_image);
3182 if (ty_image == (Image *) NULL)
3183 ThrowDPCSimilarityException();
3184 /*
3185 Compute magnitude of derivatives.
3186 */
3187 magnitude_image=SIMMagnitudeImage(tx_image,ty_image,exception);
3188 if (magnitude_image == (Image *) NULL)
3189 ThrowDPCSimilarityException();
3190 /*
3191 Divide Lx and Ly by magnitude.
3192 */
3193 temp_image=SIMDivideByMagnitude(tx_image,magnitude_image,image,exception);
3194 tx_image=DestroyImage(tx_image);
3195 if (temp_image == (Image *) NULL)
3196 ThrowDPCSimilarityException();
3197 tx_image=temp_image;
3198 try_image=SIMDivideByMagnitude(ty_image,magnitude_image,image,exception);
3199 ty_image=DestroyImage(ty_image);
3200 magnitude_image=DestroyImage(magnitude_image);
3201 if (try_image == (Image *) NULL)
3202 ThrowDPCSimilarityException();
3203 ty_image=try_image;
3204 /*
3205 Compute the cross correlation of the test and reference images.
3206 */
3207 trx_image=SIMCrossCorrelationImage(tx_image,rx_image,exception);
3208 rx_image=DestroyImage(rx_image);
3209 tx_image=DestroyImage(tx_image);
3210 if (trx_image == (Image *) NULL)
3211 ThrowDPCSimilarityException();
3212 try_image=SIMCrossCorrelationImage(ty_image,ry_image,exception);
3213 ry_image=DestroyImage(ry_image);
3214 ty_image=DestroyImage(ty_image);
3215 if (try_image == (Image *) NULL)
3216 ThrowDPCSimilarityException();
3217 /*
3218 Evaluate dot product correlation image.
3219 */
3220 (void) SetImageArtifact(try_image,"compose:clamp","False");
3221 status=CompositeImage(trx_image,try_image,PlusCompositeOp,MagickTrue,0,0,
3222 exception);
3223 try_image=DestroyImage(try_image);
3224 if (status == MagickFalse)
3225 ThrowDPCSimilarityException();
3226 status=SIMMultiplyImage(trx_image,edge_factor,
3227 (const ChannelStatistics *) NULL,exception);
3228 if (status == MagickFalse)
3229 ThrowDPCSimilarityException();
3230 /*
3231 Crop results.
3232 */
3233 SetGeometry(image,&geometry);
3234 geometry.width=image->columns;
3235 geometry.height=image->rows;
3236 (void) ResetImagePage(trx_image,"0x0+0+0");
3237puts("a");
3238 dot_product_image=CropImage(trx_image,&geometry,exception);
3239 trx_image=DestroyImage(trx_image);
3240 if (dot_product_image == (Image *) NULL)
3241 ThrowDPCSimilarityException();
3242 (void) ResetImagePage(dot_product_image,"0x0+0+0");
3243 /*
3244 Identify the maxima value in the image and its location.
3245 */
3246 status=GrayscaleImage(dot_product_image,AveragePixelIntensityMethod,
3247 exception);
3248 if (status == MagickFalse)
3249 ThrowDPCSimilarityException();
3250 dot_product_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3251 status=SIMMaximaImage(dot_product_image,&maxima,offset,exception);
3252 if (status == MagickFalse)
3253 ThrowDPCSimilarityException();
3254 *similarity_metric=1.0-QuantumScale*maxima;
3255 return(dot_product_image);
3256}
3257
3258static Image *MSESimilarityImage(const Image *image,const Image *reconstruct,
3259 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3260{
3261#define ThrowMSESimilarityException() \
3262{ \
3263 if (alpha_image != (Image *) NULL) \
3264 alpha_image=DestroyImage(alpha_image); \
3265 if (beta_image != (Image *) NULL) \
3266 beta_image=DestroyImage(beta_image); \
3267 if (channel_statistics != (ChannelStatistics *) NULL) \
3268 channel_statistics=(ChannelStatistics *) \
3269 RelinquishMagickMemory(channel_statistics); \
3270 if (mean_image != (Image *) NULL) \
3271 mean_image=DestroyImage(mean_image); \
3272 if (mse_image != (Image *) NULL) \
3273 mse_image=DestroyImage(mse_image); \
3274 if (reconstruct_image != (Image *) NULL) \
3275 reconstruct_image=DestroyImage(reconstruct_image); \
3276 if (sum_image != (Image *) NULL) \
3277 sum_image=DestroyImage(sum_image); \
3278 if (alpha_image != (Image *) NULL) \
3279 alpha_image=DestroyImage(alpha_image); \
3280 return((Image *) NULL); \
3281}
3282
3284 *channel_statistics = (ChannelStatistics *) NULL;
3285
3286 double
3287 minima = 0.0;
3288
3289 Image
3290 *alpha_image = (Image *) NULL,
3291 *beta_image = (Image *) NULL,
3292 *mean_image = (Image *) NULL,
3293 *mse_image = (Image *) NULL,
3294 *reconstruct_image = (Image *) NULL,
3295 *sum_image = (Image *) NULL,
3296 *test_image = (Image *) NULL;
3297
3298 MagickBooleanType
3299 status;
3300
3302 geometry;
3303
3304 /*
3305 MSE correlation-based image similarity using FFT local statistics.
3306 */
3307 test_image=SIMSquareImage(image,exception);
3308 if (test_image == (Image *) NULL)
3309 ThrowMSESimilarityException();
3310 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3311 if (reconstruct_image == (Image *) NULL)
3312 ThrowMSESimilarityException();
3313 /*
3314 Create (U * test)/# pixels.
3315 */
3316 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3317 test_image=DestroyImage(test_image);
3318 if (alpha_image == (Image *) NULL)
3319 ThrowMSESimilarityException();
3320 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(double)
3321 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
3322 if (status == MagickFalse)
3323 ThrowMSESimilarityException();
3324 /*
3325 Create 2*(text * reconstruction)# pixels.
3326 */
3327 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
3328 MagickTrue,0,0,exception);
3329 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3330 reconstruct_image=DestroyImage(reconstruct_image);
3331 if (beta_image == (Image *) NULL)
3332 ThrowMSESimilarityException();
3333 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(double)
3334 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
3335 if (status == MagickFalse)
3336 ThrowMSESimilarityException();
3337 /*
3338 Mean of reconstruction squared.
3339 */
3340 sum_image=SIMSquareImage(reconstruct,exception);
3341 if (sum_image == (Image *) NULL)
3342 ThrowMSESimilarityException();
3343 channel_statistics=GetImageStatistics(sum_image,exception);
3344 if (channel_statistics == (ChannelStatistics *) NULL)
3345 ThrowMSESimilarityException();
3346 status=SetImageExtent(sum_image,image->columns,image->rows,exception);
3347 if (status == MagickFalse)
3348 ThrowMSESimilarityException();
3349 status=SetImageStorageClass(sum_image,DirectClass,exception);
3350 if (status == MagickFalse)
3351 ThrowMSESimilarityException();
3352 status=SIMSetImageMean(sum_image,channel_statistics,exception);
3353 channel_statistics=(ChannelStatistics *)
3354 RelinquishMagickMemory(channel_statistics);
3355 if (status == MagickFalse)
3356 ThrowMSESimilarityException();
3357 /*
3358 Create mean image.
3359 */
3360 AppendImageToList(&sum_image,alpha_image);
3361 AppendImageToList(&sum_image,beta_image);
3362 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
3363 if (mean_image == (Image *) NULL)
3364 ThrowMSESimilarityException();
3365 sum_image=DestroyImage(sum_image);
3366 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
3367 if (status == MagickFalse)
3368 ThrowMSESimilarityException();
3369 mean_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3370 /*
3371 Crop to difference of reconstruction and test images.
3372 */
3373 SetGeometry(image,&geometry);
3374 geometry.width=image->columns;
3375 geometry.height=image->rows;
3376 (void) ResetImagePage(mean_image,"0x0+0+0");
3377 mse_image=CropImage(mean_image,&geometry,exception);
3378 mean_image=DestroyImage(mean_image);
3379 if (mse_image == (Image *) NULL)
3380 ThrowMSESimilarityException();
3381 /*
3382 Locate minimum.
3383 */
3384 (void) ResetImagePage(mse_image,"0x0+0+0");
3385 (void) ClampImage(mse_image,exception);
3386 status=SIMMinimaImage(mse_image,&minima,offset,exception);
3387 if (status == MagickFalse)
3388 ThrowMSESimilarityException();
3389 status=NegateImage(mse_image,MagickFalse,exception);
3390 if (status == MagickFalse)
3391 ThrowMSESimilarityException();
3392 *similarity_metric=QuantumScale*minima;
3393 alpha_image=DestroyImage(alpha_image);
3394 beta_image=DestroyImage(beta_image);
3395 return(mse_image);
3396}
3397
3398static Image *NCCSimilarityImage(const Image *image,const Image *reconstruct,
3399 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3400{
3401#define ThrowNCCSimilarityException() \
3402{ \
3403 if (alpha_image != (Image *) NULL) \
3404 alpha_image=DestroyImage(alpha_image); \
3405 if (beta_image != (Image *) NULL) \
3406 beta_image=DestroyImage(beta_image); \
3407 if (channel_statistics != (ChannelStatistics *) NULL) \
3408 channel_statistics=(ChannelStatistics *) \
3409 RelinquishMagickMemory(channel_statistics); \
3410 if (correlation_image != (Image *) NULL) \
3411 correlation_image=DestroyImage(correlation_image); \
3412 if (divide_image != (Image *) NULL) \
3413 divide_image=DestroyImage(divide_image); \
3414 if (ncc_image != (Image *) NULL) \
3415 ncc_image=DestroyImage(ncc_image); \
3416 if (normalize_image != (Image *) NULL) \
3417 normalize_image=DestroyImage(normalize_image); \
3418 if (reconstruct_image != (Image *) NULL) \
3419 reconstruct_image=DestroyImage(reconstruct_image); \
3420 if (test_image != (Image *) NULL) \
3421 test_image=DestroyImage(test_image); \
3422 if (variance_image != (Image *) NULL) \
3423 variance_image=DestroyImage(variance_image); \
3424 return((Image *) NULL); \
3425}
3426
3428 *channel_statistics = (ChannelStatistics *) NULL;
3429
3430 double
3431 maxima = 0.0;
3432
3433 Image
3434 *alpha_image = (Image *) NULL,
3435 *beta_image = (Image *) NULL,
3436 *correlation_image = (Image *) NULL,
3437 *divide_image = (Image *) NULL,
3438 *ncc_image = (Image *) NULL,
3439 *normalize_image = (Image *) NULL,
3440 *reconstruct_image = (Image *) NULL,
3441 *test_image = (Image *) NULL,
3442 *variance_image = (Image *) NULL;
3443
3444 MagickBooleanType
3445 status;
3446
3448 geometry;
3449
3450 /*
3451 NCC correlation-based image similarity with FFT local statistics.
3452 */
3453 test_image=SIMSquareImage(image,exception);
3454 if (test_image == (Image *) NULL)
3455 ThrowNCCSimilarityException();
3456 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3457 if (reconstruct_image == (Image *) NULL)
3458 ThrowNCCSimilarityException();
3459 /*
3460 Compute the cross correlation of the test and reconstruction images.
3461 */
3462 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3463 test_image=DestroyImage(test_image);
3464 if (alpha_image == (Image *) NULL)
3465 ThrowNCCSimilarityException();
3466 status=SIMMultiplyImage(alpha_image,(double) QuantumRange*
3467 reconstruct->columns*reconstruct->rows,(const ChannelStatistics *) NULL,
3468 exception);
3469 if (status == MagickFalse)
3470 ThrowNCCSimilarityException();
3471 /*
3472 Compute the cross correlation of the source and reconstruction images.
3473 */
3474 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3475 reconstruct_image=DestroyImage(reconstruct_image);
3476 if (beta_image == (Image *) NULL)
3477 ThrowNCCSimilarityException();
3478 test_image=SIMSquareImage(beta_image,exception);
3479 beta_image=DestroyImage(beta_image);
3480 if (test_image == (Image *) NULL)
3481 ThrowNCCSimilarityException();
3482 status=SIMMultiplyImage(test_image,(double) QuantumRange,
3483 (const ChannelStatistics *) NULL,exception);
3484 if (status == MagickFalse)
3485 ThrowNCCSimilarityException();
3486 /*
3487 Compute the variance of the two images.
3488 */
3489 variance_image=SIMVarianceImage(alpha_image,test_image,exception);
3490 test_image=DestroyImage(test_image);
3491 alpha_image=DestroyImage(alpha_image);
3492 if (variance_image == (Image *) NULL)
3493 ThrowNCCSimilarityException();
3494 /*
3495 Subtract the image mean.
3496 */
3497 channel_statistics=GetImageStatistics(reconstruct,exception);
3498 if (channel_statistics == (ChannelStatistics *) NULL)
3499 ThrowNCCSimilarityException();
3500 status=SIMMultiplyImage(variance_image,1.0,channel_statistics,exception);
3501 if (status == MagickFalse)
3502 ThrowNCCSimilarityException();
3503 normalize_image=SIMSubtractImageMean(image,reconstruct,channel_statistics,
3504 exception);
3505 channel_statistics=(ChannelStatistics *)
3506 RelinquishMagickMemory(channel_statistics);
3507 if (normalize_image == (Image *) NULL)
3508 ThrowNCCSimilarityException();
3509 correlation_image=SIMCrossCorrelationImage(image,normalize_image,exception);
3510 normalize_image=DestroyImage(normalize_image);
3511 if (correlation_image == (Image *) NULL)
3512 ThrowNCCSimilarityException();
3513 /*
3514 Divide the two images.
3515 */
3516 divide_image=SIMDivideImage(correlation_image,variance_image,exception);
3517 correlation_image=DestroyImage(correlation_image);
3518 variance_image=DestroyImage(variance_image);
3519 if (divide_image == (Image *) NULL)
3520 ThrowNCCSimilarityException();
3521 /*
3522 Crop padding.
3523 */
3524 SetGeometry(image,&geometry);
3525 geometry.width=image->columns;
3526 geometry.height=image->rows;
3527 (void) ResetImagePage(divide_image,"0x0+0+0");
3528 ncc_image=CropImage(divide_image,&geometry,exception);
3529 divide_image=DestroyImage(divide_image);
3530 if (ncc_image == (Image *) NULL)
3531 ThrowNCCSimilarityException();
3532 /*
3533 Identify the maxima value in the image and its location.
3534 */
3535 (void) ResetImagePage(ncc_image,"0x0+0+0");
3536 status=GrayscaleImage(ncc_image,AveragePixelIntensityMethod,exception);
3537 if (status == MagickFalse)
3538 ThrowNCCSimilarityException();
3539 ncc_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3540 status=SIMMaximaImage(ncc_image,&maxima,offset,exception);
3541 if (status == MagickFalse)
3542 ThrowNCCSimilarityException();
3543 *similarity_metric=1.0-QuantumScale*maxima;
3544 return(ncc_image);
3545}
3546
3547static Image *PhaseSimilarityImage(const Image *image,const Image *reconstruct,
3548 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3549{
3550#define ThrowPhaseSimilarityException() \
3551{ \
3552 if (correlation_image != (Image *) NULL) \
3553 correlation_image=DestroyImage(correlation_image); \
3554 if (fft_images != (Image *) NULL) \
3555 fft_images=DestroyImageList(fft_images); \
3556 if (gamma_image != (Image *) NULL) \
3557 gamma_image=DestroyImage(gamma_image); \
3558 if (magnitude_image != (Image *) NULL) \
3559 magnitude_image=DestroyImage(magnitude_image); \
3560 if (phase_image != (Image *) NULL) \
3561 phase_image=DestroyImage(phase_image); \
3562 if (reconstruct_image != (Image *) NULL) \
3563 reconstruct_image=DestroyImage(reconstruct_image); \
3564 if (reconstruct_magnitude != (Image *) NULL) \
3565 reconstruct_magnitude=DestroyImage(reconstruct_magnitude); \
3566 if (test_image != (Image *) NULL) \
3567 test_image=DestroyImage(test_image); \
3568 if (test_magnitude != (Image *) NULL) \
3569 test_magnitude=DestroyImage(test_magnitude); \
3570 return((Image *) NULL); \
3571}
3572
3573 double
3574 maxima = 0.0;
3575
3576 Image
3577 *correlation_image = (Image *) NULL,
3578 *fft_images = (Image *) NULL,
3579 *gamma_image = (Image *) NULL,
3580 *magnitude_image = (Image *) NULL,
3581 *phase_image = (Image *) NULL,
3582 *reconstruct_image = (Image *) NULL,
3583 *reconstruct_magnitude = (Image *) NULL,
3584 *test_image = (Image *) NULL,
3585 *test_magnitude = (Image *) NULL;
3586
3587 MagickBooleanType
3588 status;
3589
3591 geometry;
3592
3593 /*
3594 Phase correlation-based image similarity using FFT local statistics.
3595 */
3596 test_image=CloneImage(image,0,0,MagickTrue,exception);
3597 if (test_image == (Image *) NULL)
3598 ThrowPhaseSimilarityException();
3599 GetPixelInfoRGBA(0,0,0,0,&test_image->background_color);
3600 (void) ResetImagePage(test_image,"0x0+0+0");
3601 status=SetImageExtent(test_image,2*(size_t) ceil((double) image->columns/2.0),
3602 2*(size_t) ceil((double) image->rows/2.0),exception);
3603 if (status == MagickFalse)
3604 ThrowPhaseSimilarityException();
3605 (void) SetImageAlphaChannel(test_image,OffAlphaChannel,exception);
3606 /*
3607 Compute the cross correlation of the test and reconstruct magnitudes.
3608 */
3609 reconstruct_image=CloneImage(reconstruct,0,0,MagickTrue,exception);
3610 if (reconstruct_image == (Image *) NULL)
3611 ThrowPhaseSimilarityException();
3612 GetPixelInfoRGBA(0,0,0,0,&reconstruct_image->background_color);
3613 (void) ResetImagePage(reconstruct_image,"0x0+0+0");
3614 status=SetImageExtent(reconstruct_image,2*(size_t) ceil((double)
3615 image->columns/2.0),2*(size_t) ceil((double) image->rows/2.0),exception);
3616 if (status == MagickFalse)
3617 ThrowPhaseSimilarityException();
3618 (void) SetImageAlphaChannel(reconstruct_image,OffAlphaChannel,exception);
3619 (void) SetImageArtifact(test_image,"fourier:normalize","inverse");
3620 fft_images=ForwardFourierTransformImage(test_image,MagickTrue,exception);
3621 if (fft_images == (Image *) NULL)
3622 ThrowPhaseSimilarityException();
3623 test_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
3624 fft_images=DestroyImageList(fft_images);
3625 if (test_magnitude == (Image *) NULL)
3626 ThrowPhaseSimilarityException();
3627 (void) SetImageArtifact(reconstruct_image,"fourier:normalize","inverse");
3628 fft_images=ForwardFourierTransformImage(reconstruct_image,MagickTrue,
3629 exception);
3630 if (fft_images == (Image *) NULL)
3631 ThrowPhaseSimilarityException();
3632 reconstruct_magnitude=CloneImage(fft_images,0,0,MagickTrue,exception);
3633 fft_images=DestroyImageList(fft_images);
3634 if (reconstruct_magnitude == (Image *) NULL)
3635 ThrowPhaseSimilarityException();
3636 magnitude_image=CloneImage(reconstruct_magnitude,0,0,MagickTrue,exception);
3637 if (magnitude_image == (Image *) NULL)
3638 ThrowPhaseSimilarityException();
3639 DisableCompositeClampUnlessSpecified(magnitude_image);
3640 (void) CompositeImage(magnitude_image,test_magnitude,MultiplyCompositeOp,
3641 MagickTrue,0,0,exception);
3642 /*
3643 Compute the cross correlation of the test and reconstruction images.
3644 */
3645 correlation_image=SIMPhaseCorrelationImage(test_image,reconstruct_image,
3646 magnitude_image,exception);
3647 test_image=DestroyImage(test_image);
3648 reconstruct_image=DestroyImage(reconstruct_image);
3649 test_magnitude=DestroyImage(test_magnitude);
3650 reconstruct_magnitude=DestroyImage(reconstruct_magnitude);
3651 if (correlation_image == (Image *) NULL)
3652 ThrowPhaseSimilarityException();
3653 /*
3654 Identify the maxima value in the image and its location.
3655 */
3656 gamma_image=CloneImage(correlation_image,0,0,MagickTrue,exception);
3657 correlation_image=DestroyImage(correlation_image);
3658 if (gamma_image == (Image *) NULL)
3659 ThrowPhaseSimilarityException();
3660 /*
3661 Crop padding.
3662 */
3663 SetGeometry(image,&geometry);
3664 geometry.width=image->columns;
3665 geometry.height=image->rows;
3666 (void) ResetImagePage(gamma_image,"0x0+0+0");
3667 phase_image=CropImage(gamma_image,&geometry,exception);
3668 gamma_image=DestroyImage(gamma_image);
3669 if (phase_image == (Image *) NULL)
3670 ThrowPhaseSimilarityException();
3671 (void) ResetImagePage(phase_image,"0x0+0+0");
3672 /*
3673 Identify the maxima value in the image and its location.
3674 */
3675 status=GrayscaleImage(phase_image,AveragePixelIntensityMethod,exception);
3676 if (status == MagickFalse)
3677 ThrowPhaseSimilarityException();
3678 phase_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3679 status=SIMMaximaImage(phase_image,&maxima,offset,exception);
3680 if (status == MagickFalse)
3681 ThrowPhaseSimilarityException();
3682 *similarity_metric=QuantumScale*maxima;
3683 magnitude_image=DestroyImage(magnitude_image);
3684 return(phase_image);
3685}
3686
3687static Image *PSNRSimilarityImage(const Image *image,const Image *reconstruct,
3688 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3689{
3690#define ThrowPSNRSimilarityException() \
3691{ \
3692 if (alpha_image != (Image *) NULL) \
3693 alpha_image=DestroyImage(alpha_image); \
3694 if (beta_image != (Image *) NULL) \
3695 beta_image=DestroyImage(beta_image); \
3696 if (channel_statistics != (ChannelStatistics *) NULL) \
3697 channel_statistics=(ChannelStatistics *) \
3698 RelinquishMagickMemory(channel_statistics); \
3699 if (mean_image != (Image *) NULL) \
3700 mean_image=DestroyImage(mean_image); \
3701 if (psnr_image != (Image *) NULL) \
3702 psnr_image=DestroyImage(psnr_image); \
3703 if (reconstruct_image != (Image *) NULL) \
3704 reconstruct_image=DestroyImage(reconstruct_image); \
3705 if (sum_image != (Image *) NULL) \
3706 sum_image=DestroyImage(sum_image); \
3707 if (test_image != (Image *) NULL) \
3708 test_image=DestroyImage(test_image); \
3709 return((Image *) NULL); \
3710}
3711
3713 *channel_statistics = (ChannelStatistics *) NULL;
3714
3715 double
3716 minima = 0.0;
3717
3718 Image
3719 *alpha_image = (Image *) NULL,
3720 *beta_image = (Image *) NULL,
3721 *mean_image = (Image *) NULL,
3722 *psnr_image = (Image *) NULL,
3723 *reconstruct_image = (Image *) NULL,
3724 *sum_image = (Image *) NULL,
3725 *test_image = (Image *) NULL;
3726
3727 MagickBooleanType
3728 status;
3729
3731 geometry;
3732
3733 /*
3734 PSNR correlation-based image similarity using FFT local statistics.
3735 */
3736 test_image=SIMSquareImage(image,exception);
3737 if (test_image == (Image *) NULL)
3738 ThrowPSNRSimilarityException();
3739 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3740 if (reconstruct_image == (Image *) NULL)
3741 ThrowPSNRSimilarityException();
3742 /*
3743 Create (U * test)/# pixels.
3744 */
3745 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3746 test_image=DestroyImage(test_image);
3747 if (alpha_image == (Image *) NULL)
3748 ThrowPSNRSimilarityException();
3749 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(double)
3750 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
3751 if (status == MagickFalse)
3752 ThrowPSNRSimilarityException();
3753 /*
3754 Create 2*(text * reconstruction)# pixels.
3755 */
3756 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
3757 MagickTrue,0,0,exception);
3758 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3759 reconstruct_image=DestroyImage(reconstruct_image);
3760 if (beta_image == (Image *) NULL)
3761 ThrowPSNRSimilarityException();
3762 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(double)
3763 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
3764 if (status == MagickFalse)
3765 ThrowPSNRSimilarityException();
3766 /*
3767 Mean of reconstruction squared.
3768 */
3769 sum_image=SIMSquareImage(reconstruct,exception);
3770 if (sum_image == (Image *) NULL)
3771 ThrowPSNRSimilarityException();
3772 channel_statistics=GetImageStatistics(sum_image,exception);
3773 if (channel_statistics == (ChannelStatistics *) NULL)
3774 ThrowPSNRSimilarityException();
3775 status=SetImageExtent(sum_image,image->columns,image->rows,exception);
3776 if (status == MagickFalse)
3777 ThrowPSNRSimilarityException();
3778 status=SetImageStorageClass(sum_image,DirectClass,exception);
3779 if (status == MagickFalse)
3780 ThrowPSNRSimilarityException();
3781 status=SIMSetImageMean(sum_image,channel_statistics,exception);
3782 channel_statistics=(ChannelStatistics *)
3783 RelinquishMagickMemory(channel_statistics);
3784 if (status == MagickFalse)
3785 ThrowPSNRSimilarityException();
3786 /*
3787 Create mean image.
3788 */
3789 AppendImageToList(&sum_image,alpha_image);
3790 AppendImageToList(&sum_image,beta_image);
3791 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
3792 if (mean_image == (Image *) NULL)
3793 ThrowPSNRSimilarityException();
3794 sum_image=DestroyImage(sum_image);
3795 status=SIMLogImage(mean_image,exception);
3796 if (status == MagickFalse)
3797 ThrowPSNRSimilarityException();
3798 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
3799 if (status == MagickFalse)
3800 ThrowPSNRSimilarityException();
3801 mean_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3802 status=SIMMultiplyImage(mean_image,1.0/48.1647,
3803 (const ChannelStatistics *) NULL,exception);
3804 if (status == MagickFalse)
3805 ThrowPSNRSimilarityException();
3806 /*
3807 Crop to difference of reconstruction and test images.
3808 */
3809 SetGeometry(image,&geometry);
3810 geometry.width=image->columns;
3811 geometry.height=image->rows;
3812 (void) ResetImagePage(mean_image,"0x0+0+0");
3813 psnr_image=CropImage(mean_image,&geometry,exception);
3814 mean_image=DestroyImage(mean_image);
3815 if (psnr_image == (Image *) NULL)
3816 ThrowPSNRSimilarityException();
3817 /*
3818 Locate minimum.
3819 */
3820 (void) ResetImagePage(psnr_image,"0x0+0+0");
3821 (void) EvaluateImage(psnr_image,MaxEvaluateOperator,0.0,exception);
3822 status=SIMMinimaImage(psnr_image,&minima,offset,exception);
3823 if (status == MagickFalse)
3824 ThrowPSNRSimilarityException();
3825 status=NegateImage(psnr_image,MagickFalse,exception);
3826 if (status == MagickFalse)
3827 ThrowPSNRSimilarityException();
3828 *similarity_metric=(-1.0+QuantumScale*minima);
3829 alpha_image=DestroyImage(alpha_image);
3830 beta_image=DestroyImage(beta_image);
3831 return(psnr_image);
3832}
3833
3834static Image *RMSESimilarityImage(const Image *image,const Image *reconstruct,
3835 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
3836{
3837#define ThrowRMSESimilarityException() \
3838{ \
3839 if (alpha_image != (Image *) NULL) \
3840 alpha_image=DestroyImage(alpha_image); \
3841 if (beta_image != (Image *) NULL) \
3842 beta_image=DestroyImage(beta_image); \
3843 if (channel_statistics != (ChannelStatistics *) NULL) \
3844 channel_statistics=(ChannelStatistics *) \
3845 RelinquishMagickMemory(channel_statistics); \
3846 if (mean_image != (Image *) NULL) \
3847 mean_image=DestroyImage(mean_image); \
3848 if (rmse_image != (Image *) NULL) \
3849 rmse_image=DestroyImage(rmse_image); \
3850 if (reconstruct_image != (Image *) NULL) \
3851 reconstruct_image=DestroyImage(reconstruct_image); \
3852 if (sum_image != (Image *) NULL) \
3853 sum_image=DestroyImage(sum_image); \
3854 if (alpha_image != (Image *) NULL) \
3855 alpha_image=DestroyImage(alpha_image); \
3856 return((Image *) NULL); \
3857}
3858
3860 *channel_statistics = (ChannelStatistics *) NULL;
3861
3862 double
3863 minima = 0.0;
3864
3865 Image
3866 *alpha_image = (Image *) NULL,
3867 *beta_image = (Image *) NULL,
3868 *mean_image = (Image *) NULL,
3869 *reconstruct_image = (Image *) NULL,
3870 *rmse_image = (Image *) NULL,
3871 *sum_image = (Image *) NULL,
3872 *test_image = (Image *) NULL;
3873
3874 MagickBooleanType
3875 status;
3876
3878 geometry;
3879
3880 /*
3881 RMSE correlation-based image similarity using FFT local statistics.
3882 */
3883 test_image=SIMSquareImage(image,exception);
3884 if (test_image == (Image *) NULL)
3885 ThrowRMSESimilarityException();
3886 reconstruct_image=SIMUnityImage(image,reconstruct,exception);
3887 if (reconstruct_image == (Image *) NULL)
3888 ThrowRMSESimilarityException();
3889 /*
3890 Create (U * test)/# pixels.
3891 */
3892 alpha_image=SIMCrossCorrelationImage(test_image,reconstruct_image,exception);
3893 test_image=DestroyImage(test_image);
3894 if (alpha_image == (Image *) NULL)
3895 ThrowRMSESimilarityException();
3896 status=SIMMultiplyImage(alpha_image,1.0/reconstruct->columns/(double)
3897 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
3898 if (status == MagickFalse)
3899 ThrowRMSESimilarityException();
3900 /*
3901 Create 2*(text * reconstruction)# pixels.
3902 */
3903 (void) CompositeImage(reconstruct_image,reconstruct,CopyCompositeOp,
3904 MagickTrue,0,0,exception);
3905 beta_image=SIMCrossCorrelationImage(image,reconstruct_image,exception);
3906 reconstruct_image=DestroyImage(reconstruct_image);
3907 if (beta_image == (Image *) NULL)
3908 ThrowRMSESimilarityException();
3909 status=SIMMultiplyImage(beta_image,-2.0/reconstruct->columns/(double)
3910 reconstruct->rows,(const ChannelStatistics *) NULL,exception);
3911 if (status == MagickFalse)
3912 ThrowRMSESimilarityException();
3913 /*
3914 Mean of reconstruction squared.
3915 */
3916 sum_image=SIMSquareImage(reconstruct,exception);
3917 if (sum_image == (Image *) NULL)
3918 ThrowRMSESimilarityException();
3919 channel_statistics=GetImageStatistics(sum_image,exception);
3920 if (channel_statistics == (ChannelStatistics *) NULL)
3921 ThrowRMSESimilarityException();
3922 status=SetImageExtent(sum_image,image->columns,image->rows,exception);
3923 if (status == MagickFalse)
3924 ThrowRMSESimilarityException();
3925 status=SetImageStorageClass(sum_image,DirectClass,exception);
3926 if (status == MagickFalse)
3927 ThrowRMSESimilarityException();
3928 status=SIMSetImageMean(sum_image,channel_statistics,exception);
3929 channel_statistics=(ChannelStatistics *)
3930 RelinquishMagickMemory(channel_statistics);
3931 if (status == MagickFalse)
3932 ThrowRMSESimilarityException();
3933 /*
3934 Create mean image.
3935 */
3936 AppendImageToList(&sum_image,alpha_image);
3937 AppendImageToList(&sum_image,beta_image);
3938 mean_image=EvaluateImages(sum_image,SumEvaluateOperator,exception);
3939 if (mean_image == (Image *) NULL)
3940 ThrowRMSESimilarityException();
3941 status=EvaluateImage(mean_image,PowEvaluateOperator,0.5,exception);
3942 if (mean_image == (Image *) NULL)
3943 ThrowRMSESimilarityException();
3944 sum_image=DestroyImage(sum_image);
3945 status=GrayscaleImage(mean_image,AveragePixelIntensityMethod,exception);
3946 if (status == MagickFalse)
3947 ThrowRMSESimilarityException();
3948 mean_image->depth=MAGICKCORE_QUANTUM_DEPTH;
3949 /*
3950 Crop to difference of reconstruction and test images.
3951 */
3952 SetGeometry(image,&geometry);
3953 geometry.width=image->columns;
3954 geometry.height=image->rows;
3955 (void) ResetImagePage(mean_image,"0x0+0+0");
3956 rmse_image=CropImage(mean_image,&geometry,exception);
3957 mean_image=DestroyImage(mean_image);
3958 if (rmse_image == (Image *) NULL)
3959 ThrowRMSESimilarityException();
3960 /*
3961 Locate minimum.
3962 */
3963 (void) ResetImagePage(rmse_image,"0x0+0+0");
3964 status=SIMMinimaImage(rmse_image,&minima,offset,exception);
3965 if (status == MagickFalse)
3966 ThrowRMSESimilarityException();
3967 status=NegateImage(rmse_image,MagickFalse,exception);
3968 if (status == MagickFalse)
3969 ThrowRMSESimilarityException();
3970 *similarity_metric=QuantumScale*minima;
3971 alpha_image=DestroyImage(alpha_image);
3972 beta_image=DestroyImage(beta_image);
3973 return(rmse_image);
3974}
3975
3976#endif
3977
3978static double GetSimilarityMetric(const Image *image,const Image *reconstruct,
3979 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
3980 ExceptionInfo *exception)
3981{
3982 double
3983 distortion;
3984
3986 *sans_exception = AcquireExceptionInfo();
3987
3988 Image
3989 *similarity_image;
3990
3991 MagickBooleanType
3992 status;
3993
3995 geometry;
3996
3997 SetGeometry(reconstruct,&geometry);
3998 geometry.x=x_offset;
3999 geometry.y=y_offset;
4000 similarity_image=CropImage(image,&geometry,sans_exception);
4001 sans_exception=DestroyExceptionInfo(sans_exception);
4002 if (similarity_image == (Image *) NULL)
4003 return(0.0);
4004 distortion=0.0;
4005 status=GetImageDistortion(similarity_image,reconstruct,metric,&distortion,
4006 exception);
4007 similarity_image=DestroyImage(similarity_image);
4008 if (status == MagickFalse)
4009 return(0.0);
4010 return(distortion);
4011}
4012
4013MagickExport Image *SimilarityImage(const Image *image,const Image *reconstruct,
4014 const MetricType metric,const double similarity_threshold,
4015 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
4016{
4017#define SimilarityImageTag "Similarity/Image"
4018
4019 CacheView
4020 *similarity_view;
4021
4022 Image
4023 *similarity_image = (Image *) NULL;
4024
4025 MagickBooleanType
4026 status;
4027
4028 MagickOffsetType
4029 progress;
4030
4031 size_t
4032 rows;
4033
4034 ssize_t
4035 y;
4036
4037 assert(image != (const Image *) NULL);
4038 assert(image->signature == MagickCoreSignature);
4039 assert(exception != (ExceptionInfo *) NULL);
4040 assert(exception->signature == MagickCoreSignature);
4041 assert(offset != (RectangleInfo *) NULL);
4042 if (IsEventLogging() != MagickFalse)
4043 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4044 SetGeometry(reconstruct,offset);
4045 *similarity_metric=MagickMaximumValue;
4046#if defined(MAGICKCORE_HDRI_SUPPORT) && defined(MAGICKCORE_FFTW_DELEGATE)
4047{
4048 const char *artifact = GetImageArtifact(image,"compare:frequency-domain");
4049 if (artifact == (const char *) NULL)
4050 artifact=GetImageArtifact(image,"compare:accelerate-ncc");
4051 if ((image->channels & ReadMaskChannel) == 0)
4052 switch (metric)
4053 {
4054 case DotProductCorrelationErrorMetric:
4055 {
4056 similarity_image=DPCSimilarityImage(image,reconstruct,offset,
4057 similarity_metric,exception);
4058 return(similarity_image);
4059 }
4060 case MeanSquaredErrorMetric:
4061 {
4062 if ((artifact != (const char *) NULL) &&
4063 (IsStringTrue(artifact) == MagickFalse))
4064 break;
4065 similarity_image=MSESimilarityImage(image,reconstruct,offset,
4066 similarity_metric,exception);
4067 return(similarity_image);
4068 }
4069 case NormalizedCrossCorrelationErrorMetric:
4070 {
4071 if ((artifact != (const char *) NULL) &&
4072 (IsStringTrue(artifact) == MagickFalse))
4073 break;
4074 similarity_image=NCCSimilarityImage(image,reconstruct,offset,
4075 similarity_metric,exception);
4076 return(similarity_image);
4077 }
4078 case PeakSignalToNoiseRatioErrorMetric:
4079 {
4080 if ((artifact != (const char *) NULL) &&
4081 (IsStringTrue(artifact) == MagickFalse))
4082 break;
4083 similarity_image=PSNRSimilarityImage(image,reconstruct,offset,
4084 similarity_metric,exception);
4085 return(similarity_image);
4086 }
4087 case PhaseCorrelationErrorMetric:
4088 {
4089 similarity_image=PhaseSimilarityImage(image,reconstruct,offset,
4090 similarity_metric,exception);
4091 return(similarity_image);
4092 }
4093 case RootMeanSquaredErrorMetric:
4094 {
4095 if ((artifact != (const char *) NULL) &&
4096 (IsStringTrue(artifact) == MagickFalse))
4097 break;
4098 similarity_image=RMSESimilarityImage(image,reconstruct,offset,
4099 similarity_metric,exception);
4100 return(similarity_image);
4101 }
4102 default: break;
4103 }
4104}
4105#else
4106 if ((metric == DotProductCorrelationErrorMetric) ||
4107 (metric == PhaseCorrelationErrorMetric))
4108 {
4109 (void) ThrowMagickException(exception,GetMagickModule(),
4110 MissingDelegateError,"DelegateLibrarySupportNotBuiltIn",
4111 "'%s' (HDRI, FFT)",image->filename);
4112 return((Image *) NULL);
4113 }
4114#endif
4115 if ((image->columns < reconstruct->columns) ||
4116 (image->rows < reconstruct->rows))
4117 {
4118 (void) ThrowMagickException(exception,GetMagickModule(),OptionWarning,
4119 "GeometryDoesNotContainImage","`%s'",image->filename);
4120 return((Image *) NULL);
4121 }
4122 similarity_image=CloneImage(image,image->columns,image->rows,MagickTrue,
4123 exception);
4124 if (similarity_image == (Image *) NULL)
4125 return((Image *) NULL);
4126 similarity_image->depth=MAGICKCORE_QUANTUM_DEPTH;
4127 similarity_image->alpha_trait=UndefinedPixelTrait;
4128 similarity_image->type=GrayscaleType;
4129 status=SetImageStorageClass(similarity_image,DirectClass,exception);
4130 if (status == MagickFalse)
4131 return(DestroyImage(similarity_image));
4132 /*
4133 Measure similarity of reconstruction image against image.
4134 */
4135 status=MagickTrue;
4136 progress=0;
4137 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
4138 rows=similarity_image->rows;
4139#if defined(MAGICKCORE_OPENMP_SUPPORT)
4140 #pragma omp parallel for schedule(static,1) \
4141 shared(progress,similarity_metric,status) \
4142 magick_number_threads(similarity_image,similarity_image,rows << 3,1)
4143#endif
4144 for (y=0; y < (ssize_t) rows; y++)
4145 {
4146 double
4147 similarity;
4148
4149 Quantum
4150 *magick_restrict q;
4151
4152 ssize_t
4153 x;
4154
4155 if (status == MagickFalse)
4156 continue;
4157#if defined(MAGICKCORE_OPENMP_SUPPORT)
4158 #pragma omp flush(similarity_metric)
4159#endif
4160 if (*similarity_metric <= similarity_threshold)
4161 continue;
4162 q=QueueCacheViewAuthenticPixels(similarity_view,0,y,
4163 similarity_image->columns,1,exception);
4164 if (q == (Quantum *) NULL)
4165 {
4166 status=MagickFalse;
4167 continue;
4168 }
4169 for (x=0; x < (ssize_t) similarity_image->columns; x++)
4170 {
4171 ssize_t
4172 i;
4173
4174#if defined(MAGICKCORE_OPENMP_SUPPORT)
4175 #pragma omp flush(similarity_metric)
4176#endif
4177 if (*similarity_metric <= similarity_threshold)
4178 break;
4179 similarity=GetSimilarityMetric(image,reconstruct,metric,x,y,exception);
4180 if ((metric == DotProductCorrelationErrorMetric) ||
4181 (metric == PhaseCorrelationErrorMetric) ||
4182 (metric == NormalizedCrossCorrelationErrorMetric) ||
4183 (metric == StructuralSimilarityErrorMetric) ||
4184 (metric == UndefinedErrorMetric))
4185 similarity=1.0-similarity;
4186 if (metric == PerceptualHashErrorMetric)
4187 similarity=MagickMin(0.01*similarity,1.0);
4188#if defined(MAGICKCORE_OPENMP_SUPPORT)
4189 #pragma omp critical (MagickCore_SimilarityImage)
4190#endif
4191 if (similarity < *similarity_metric)
4192 {
4193 offset->x=x;
4194 offset->y=y;
4195 *similarity_metric=similarity;
4196 }
4197 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4198 {
4199 PixelChannel channel = GetPixelChannelChannel(image,i);
4200 PixelTrait traits = GetPixelChannelTraits(image,channel);
4201 PixelTrait similarity_traits = GetPixelChannelTraits(similarity_image,
4202 channel);
4203 if ((traits == UndefinedPixelTrait) ||
4204 (similarity_traits == UndefinedPixelTrait) ||
4205 ((similarity_traits & UpdatePixelTrait) == 0))
4206 continue;
4207 switch (metric)
4208 {
4209 case FuzzErrorMetric:
4210 case MeanAbsoluteErrorMetric:
4211 case MeanSquaredErrorMetric:
4212 case NormalizedCrossCorrelationErrorMetric:
4213 case PerceptualHashErrorMetric:
4214 case RootMeanSquaredErrorMetric:
4215 case StructuralSimilarityErrorMetric:
4216 case StructuralDissimilarityErrorMetric:
4217 {
4218 SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4219 QuantumRange-QuantumRange*similarity),q);
4220 break;
4221 }
4222 default:
4223 {
4224 SetPixelChannel(similarity_image,channel,ClampToQuantum((double)
4225 QuantumRange*similarity),q);
4226 break;
4227 }
4228 }
4229 }
4230 q+=(ptrdiff_t) GetPixelChannels(similarity_image);
4231 }
4232 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
4233 status=MagickFalse;
4234 if (image->progress_monitor != (MagickProgressMonitor) NULL)
4235 {
4236 MagickBooleanType
4237 proceed;
4238
4239#if defined(MAGICKCORE_OPENMP_SUPPORT)
4240 #pragma omp atomic
4241#endif
4242 progress++;
4243 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
4244 if (proceed == MagickFalse)
4245 status=MagickFalse;
4246 }
4247 }
4248 similarity_view=DestroyCacheView(similarity_view);
4249 if (status == MagickFalse)
4250 similarity_image=DestroyImage(similarity_image);
4251 return(similarity_image);
4252}