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