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