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-2021 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/exception-private.h"
57 #include "MagickCore/geometry.h"
58 #include "MagickCore/image-private.h"
59 #include "MagickCore/list.h"
60 #include "MagickCore/log.h"
61 #include "MagickCore/memory_.h"
62 #include "MagickCore/monitor.h"
63 #include "MagickCore/monitor-private.h"
64 #include "MagickCore/option.h"
65 #include "MagickCore/pixel-accessor.h"
66 #include "MagickCore/property.h"
67 #include "MagickCore/resource_.h"
68 #include "MagickCore/string_.h"
69 #include "MagickCore/statistic.h"
70 #include "MagickCore/string-private.h"
71 #include "MagickCore/thread-private.h"
72 #include "MagickCore/transform.h"
73 #include "MagickCore/utility.h"
74 #include "MagickCore/version.h"
75
76 /*
77 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78 % %
79 % %
80 % %
81 % C o m p a r e I m a g e s %
82 % %
83 % %
84 % %
85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86 %
87 % CompareImages() compares one or more pixel channels of an image to a
88 % reconstructed image and returns the difference image.
89 %
90 % The format of the CompareImages method is:
91 %
92 % Image *CompareImages(const Image *image,const Image *reconstruct_image,
93 % const MetricType metric,double *distortion,ExceptionInfo *exception)
94 %
95 % A description of each parameter follows:
96 %
97 % o image: the image.
98 %
99 % o reconstruct_image: the reconstruct image.
100 %
101 % o metric: the metric.
102 %
103 % o distortion: the computed distortion between the images.
104 %
105 % o exception: return any errors or warnings in this structure.
106 %
107 */
108
GetImageChannels(const Image * image)109 static size_t GetImageChannels(const Image *image)
110 {
111 ssize_t
112 i;
113
114 size_t
115 channels;
116
117 channels=0;
118 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
119 {
120 PixelChannel channel = GetPixelChannelChannel(image,i);
121 PixelTrait traits = GetPixelChannelTraits(image,channel);
122 if ((traits & UpdatePixelTrait) != 0)
123 channels++;
124 }
125 return(channels == 0 ? (size_t) 1 : channels);
126 }
127
CompareImages(Image * image,const Image * reconstruct_image,const MetricType metric,double * distortion,ExceptionInfo * exception)128 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
129 const MetricType metric,double *distortion,ExceptionInfo *exception)
130 {
131 CacheView
132 *highlight_view,
133 *image_view,
134 *reconstruct_view;
135
136 const char
137 *artifact;
138
139 double
140 fuzz;
141
142 Image
143 *clone_image,
144 *difference_image,
145 *highlight_image;
146
147 MagickBooleanType
148 status;
149
150 PixelInfo
151 highlight,
152 lowlight,
153 masklight;
154
155 RectangleInfo
156 geometry;
157
158 size_t
159 columns,
160 rows;
161
162 ssize_t
163 y;
164
165 assert(image != (Image *) NULL);
166 assert(image->signature == MagickCoreSignature);
167 if (image->debug != MagickFalse)
168 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
169 assert(reconstruct_image != (const Image *) NULL);
170 assert(reconstruct_image->signature == MagickCoreSignature);
171 assert(distortion != (double *) NULL);
172 *distortion=0.0;
173 if (image->debug != MagickFalse)
174 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
175 status=GetImageDistortion(image,reconstruct_image,metric,distortion,
176 exception);
177 if (status == MagickFalse)
178 return((Image *) NULL);
179 columns=MagickMax(image->columns,reconstruct_image->columns);
180 rows=MagickMax(image->rows,reconstruct_image->rows);
181 SetGeometry(image,&geometry);
182 geometry.width=columns;
183 geometry.height=rows;
184 clone_image=CloneImage(image,0,0,MagickTrue,exception);
185 if (clone_image == (Image *) NULL)
186 return((Image *) NULL);
187 (void) SetImageMask(clone_image,ReadPixelMask,(Image *) NULL,exception);
188 difference_image=ExtentImage(clone_image,&geometry,exception);
189 clone_image=DestroyImage(clone_image);
190 if (difference_image == (Image *) NULL)
191 return((Image *) NULL);
192 (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel,exception);
193 highlight_image=CloneImage(image,columns,rows,MagickTrue,exception);
194 if (highlight_image == (Image *) NULL)
195 {
196 difference_image=DestroyImage(difference_image);
197 return((Image *) NULL);
198 }
199 status=SetImageStorageClass(highlight_image,DirectClass,exception);
200 if (status == MagickFalse)
201 {
202 difference_image=DestroyImage(difference_image);
203 highlight_image=DestroyImage(highlight_image);
204 return((Image *) NULL);
205 }
206 (void) SetImageMask(highlight_image,ReadPixelMask,(Image *) NULL,exception);
207 (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel,exception);
208 (void) QueryColorCompliance("#f1001ecc",AllCompliance,&highlight,exception);
209 artifact=GetImageArtifact(image,"compare:highlight-color");
210 if (artifact != (const char *) NULL)
211 (void) QueryColorCompliance(artifact,AllCompliance,&highlight,exception);
212 (void) QueryColorCompliance("#ffffffcc",AllCompliance,&lowlight,exception);
213 artifact=GetImageArtifact(image,"compare:lowlight-color");
214 if (artifact != (const char *) NULL)
215 (void) QueryColorCompliance(artifact,AllCompliance,&lowlight,exception);
216 (void) QueryColorCompliance("#888888cc",AllCompliance,&masklight,exception);
217 artifact=GetImageArtifact(image,"compare:masklight-color");
218 if (artifact != (const char *) NULL)
219 (void) QueryColorCompliance(artifact,AllCompliance,&masklight,exception);
220 /*
221 Generate difference image.
222 */
223 status=MagickTrue;
224 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
225 image_view=AcquireVirtualCacheView(image,exception);
226 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
227 highlight_view=AcquireAuthenticCacheView(highlight_image,exception);
228 #if defined(MAGICKCORE_OPENMP_SUPPORT)
229 #pragma omp parallel for schedule(static) shared(status) \
230 magick_number_threads(image,highlight_image,rows,1)
231 #endif
232 for (y=0; y < (ssize_t) rows; y++)
233 {
234 MagickBooleanType
235 sync;
236
237 const Quantum
238 *magick_restrict p,
239 *magick_restrict q;
240
241 Quantum
242 *magick_restrict r;
243
244 ssize_t
245 x;
246
247 if (status == MagickFalse)
248 continue;
249 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
250 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
251 r=QueueCacheViewAuthenticPixels(highlight_view,0,y,columns,1,exception);
252 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL) ||
253 (r == (Quantum *) NULL))
254 {
255 status=MagickFalse;
256 continue;
257 }
258 for (x=0; x < (ssize_t) columns; x++)
259 {
260 double
261 Da,
262 Sa;
263
264 MagickStatusType
265 difference;
266
267 ssize_t
268 i;
269
270 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
271 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
272 {
273 SetPixelViaPixelInfo(highlight_image,&masklight,r);
274 p+=GetPixelChannels(image);
275 q+=GetPixelChannels(reconstruct_image);
276 r+=GetPixelChannels(highlight_image);
277 continue;
278 }
279 difference=MagickFalse;
280 Sa=QuantumScale*GetPixelAlpha(image,p);
281 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
282 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
283 {
284 double
285 distance,
286 pixel;
287
288 PixelChannel channel = GetPixelChannelChannel(image,i);
289 PixelTrait traits = GetPixelChannelTraits(image,channel);
290 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
291 channel);
292 if ((traits == UndefinedPixelTrait) ||
293 (reconstruct_traits == UndefinedPixelTrait) ||
294 ((reconstruct_traits & UpdatePixelTrait) == 0))
295 continue;
296 if (channel == AlphaPixelChannel)
297 pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
298 else
299 pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
300 distance=pixel*pixel;
301 if (distance >= fuzz)
302 {
303 difference=MagickTrue;
304 break;
305 }
306 }
307 if (difference == MagickFalse)
308 SetPixelViaPixelInfo(highlight_image,&lowlight,r);
309 else
310 SetPixelViaPixelInfo(highlight_image,&highlight,r);
311 p+=GetPixelChannels(image);
312 q+=GetPixelChannels(reconstruct_image);
313 r+=GetPixelChannels(highlight_image);
314 }
315 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
316 if (sync == MagickFalse)
317 status=MagickFalse;
318 }
319 highlight_view=DestroyCacheView(highlight_view);
320 reconstruct_view=DestroyCacheView(reconstruct_view);
321 image_view=DestroyCacheView(image_view);
322 (void) CompositeImage(difference_image,highlight_image,image->compose,
323 MagickTrue,0,0,exception);
324 highlight_image=DestroyImage(highlight_image);
325 if (status == MagickFalse)
326 difference_image=DestroyImage(difference_image);
327 return(difference_image);
328 }
329
330 /*
331 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
332 % %
333 % %
334 % %
335 % G e t I m a g e D i s t o r t i o n %
336 % %
337 % %
338 % %
339 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
340 %
341 % GetImageDistortion() compares one or more pixel channels of an image to a
342 % reconstructed image and returns the specified distortion metric.
343 %
344 % The format of the GetImageDistortion method is:
345 %
346 % MagickBooleanType GetImageDistortion(const Image *image,
347 % const Image *reconstruct_image,const MetricType metric,
348 % double *distortion,ExceptionInfo *exception)
349 %
350 % A description of each parameter follows:
351 %
352 % o image: the image.
353 %
354 % o reconstruct_image: the reconstruct image.
355 %
356 % o metric: the metric.
357 %
358 % o distortion: the computed distortion between the images.
359 %
360 % o exception: return any errors or warnings in this structure.
361 %
362 */
363
GetAbsoluteDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)364 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
365 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
366 {
367 CacheView
368 *image_view,
369 *reconstruct_view;
370
371 double
372 fuzz;
373
374 MagickBooleanType
375 status;
376
377 size_t
378 columns,
379 rows;
380
381 ssize_t
382 y;
383
384 /*
385 Compute the absolute difference in pixels between two images.
386 */
387 status=MagickTrue;
388 fuzz=GetFuzzyColorDistance(image,reconstruct_image);
389 rows=MagickMax(image->rows,reconstruct_image->rows);
390 columns=MagickMax(image->columns,reconstruct_image->columns);
391 image_view=AcquireVirtualCacheView(image,exception);
392 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
393 #if defined(MAGICKCORE_OPENMP_SUPPORT)
394 #pragma omp parallel for schedule(static) shared(status) \
395 magick_number_threads(image,image,rows,1)
396 #endif
397 for (y=0; y < (ssize_t) rows; y++)
398 {
399 double
400 channel_distortion[MaxPixelChannels+1];
401
402 const Quantum
403 *magick_restrict p,
404 *magick_restrict q;
405
406 ssize_t
407 j,
408 x;
409
410 if (status == MagickFalse)
411 continue;
412 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
413 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
414 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
415 {
416 status=MagickFalse;
417 continue;
418 }
419 (void) memset(channel_distortion,0,sizeof(channel_distortion));
420 for (x=0; x < (ssize_t) columns; x++)
421 {
422 double
423 Da,
424 Sa;
425
426 MagickBooleanType
427 difference;
428
429 ssize_t
430 i;
431
432 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
433 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
434 {
435 p+=GetPixelChannels(image);
436 q+=GetPixelChannels(reconstruct_image);
437 continue;
438 }
439 difference=MagickFalse;
440 Sa=QuantumScale*GetPixelAlpha(image,p);
441 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
442 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
443 {
444 double
445 distance,
446 pixel;
447
448 PixelChannel channel = GetPixelChannelChannel(image,i);
449 PixelTrait traits = GetPixelChannelTraits(image,channel);
450 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
451 channel);
452 if ((traits == UndefinedPixelTrait) ||
453 (reconstruct_traits == UndefinedPixelTrait) ||
454 ((reconstruct_traits & UpdatePixelTrait) == 0))
455 continue;
456 if (channel == AlphaPixelChannel)
457 pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
458 else
459 pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
460 distance=pixel*pixel;
461 if (distance >= fuzz)
462 {
463 channel_distortion[i]++;
464 difference=MagickTrue;
465 }
466 }
467 if (difference != MagickFalse)
468 channel_distortion[CompositePixelChannel]++;
469 p+=GetPixelChannels(image);
470 q+=GetPixelChannels(reconstruct_image);
471 }
472 #if defined(MAGICKCORE_OPENMP_SUPPORT)
473 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
474 #endif
475 for (j=0; j <= MaxPixelChannels; j++)
476 distortion[j]+=channel_distortion[j];
477 }
478 reconstruct_view=DestroyCacheView(reconstruct_view);
479 image_view=DestroyCacheView(image_view);
480 return(status);
481 }
482
GetFuzzDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)483 static MagickBooleanType GetFuzzDistortion(const Image *image,
484 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
485 {
486 CacheView
487 *image_view,
488 *reconstruct_view;
489
490 double
491 area;
492
493 MagickBooleanType
494 status;
495
496 ssize_t
497 j;
498
499 size_t
500 columns,
501 rows;
502
503 ssize_t
504 y;
505
506 status=MagickTrue;
507 rows=MagickMax(image->rows,reconstruct_image->rows);
508 columns=MagickMax(image->columns,reconstruct_image->columns);
509 area=0.0;
510 image_view=AcquireVirtualCacheView(image,exception);
511 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
512 #if defined(MAGICKCORE_OPENMP_SUPPORT)
513 #pragma omp parallel for schedule(static) shared(status) \
514 magick_number_threads(image,image,rows,1) reduction(+:area)
515 #endif
516 for (y=0; y < (ssize_t) rows; y++)
517 {
518 double
519 channel_distortion[MaxPixelChannels+1];
520
521 const Quantum
522 *magick_restrict p,
523 *magick_restrict q;
524
525 ssize_t
526 x;
527
528 if (status == MagickFalse)
529 continue;
530 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
531 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
532 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
533 {
534 status=MagickFalse;
535 continue;
536 }
537 (void) memset(channel_distortion,0,sizeof(channel_distortion));
538 for (x=0; x < (ssize_t) columns; x++)
539 {
540 double
541 Da,
542 Sa;
543
544 ssize_t
545 i;
546
547 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
548 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
549 {
550 p+=GetPixelChannels(image);
551 q+=GetPixelChannels(reconstruct_image);
552 continue;
553 }
554 Sa=QuantumScale*GetPixelAlpha(image,p);
555 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
556 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
557 {
558 double
559 distance;
560
561 PixelChannel channel = GetPixelChannelChannel(image,i);
562 PixelTrait traits = GetPixelChannelTraits(image,channel);
563 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
564 channel);
565 if ((traits == UndefinedPixelTrait) ||
566 (reconstruct_traits == UndefinedPixelTrait) ||
567 ((reconstruct_traits & UpdatePixelTrait) == 0))
568 continue;
569 if (channel == AlphaPixelChannel)
570 distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
571 channel,q));
572 else
573 distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
574 channel,q));
575 channel_distortion[i]+=distance*distance;
576 channel_distortion[CompositePixelChannel]+=distance*distance;
577 }
578 area++;
579 p+=GetPixelChannels(image);
580 q+=GetPixelChannels(reconstruct_image);
581 }
582 #if defined(MAGICKCORE_OPENMP_SUPPORT)
583 #pragma omp critical (MagickCore_GetFuzzDistortion)
584 #endif
585 for (j=0; j <= MaxPixelChannels; j++)
586 distortion[j]+=channel_distortion[j];
587 }
588 reconstruct_view=DestroyCacheView(reconstruct_view);
589 image_view=DestroyCacheView(image_view);
590 area=PerceptibleReciprocal(area);
591 for (j=0; j <= MaxPixelChannels; j++)
592 distortion[j]*=area;
593 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
594 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
595 return(status);
596 }
597
GetMeanAbsoluteDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)598 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
599 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
600 {
601 CacheView
602 *image_view,
603 *reconstruct_view;
604
605 double
606 area;
607
608 MagickBooleanType
609 status;
610
611 ssize_t
612 j;
613
614 size_t
615 columns,
616 rows;
617
618 ssize_t
619 y;
620
621 status=MagickTrue;
622 rows=MagickMax(image->rows,reconstruct_image->rows);
623 columns=MagickMax(image->columns,reconstruct_image->columns);
624 area=0.0;
625 image_view=AcquireVirtualCacheView(image,exception);
626 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
627 #if defined(MAGICKCORE_OPENMP_SUPPORT)
628 #pragma omp parallel for schedule(static) shared(status) \
629 magick_number_threads(image,image,rows,1) reduction(+:area)
630 #endif
631 for (y=0; y < (ssize_t) rows; y++)
632 {
633 double
634 channel_distortion[MaxPixelChannels+1];
635
636 const Quantum
637 *magick_restrict p,
638 *magick_restrict q;
639
640 ssize_t
641 x;
642
643 if (status == MagickFalse)
644 continue;
645 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
646 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
647 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
648 {
649 status=MagickFalse;
650 continue;
651 }
652 (void) memset(channel_distortion,0,sizeof(channel_distortion));
653 for (x=0; x < (ssize_t) columns; x++)
654 {
655 double
656 Da,
657 Sa;
658
659 ssize_t
660 i;
661
662 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
663 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
664 {
665 p+=GetPixelChannels(image);
666 q+=GetPixelChannels(reconstruct_image);
667 continue;
668 }
669 Sa=QuantumScale*GetPixelAlpha(image,p);
670 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
671 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
672 {
673 double
674 distance;
675
676 PixelChannel channel = GetPixelChannelChannel(image,i);
677 PixelTrait traits = GetPixelChannelTraits(image,channel);
678 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
679 channel);
680 if ((traits == UndefinedPixelTrait) ||
681 (reconstruct_traits == UndefinedPixelTrait) ||
682 ((reconstruct_traits & UpdatePixelTrait) == 0))
683 continue;
684 if (channel == AlphaPixelChannel)
685 distance=QuantumScale*fabs((double) (p[i]-(double)
686 GetPixelChannel(reconstruct_image,channel,q)));
687 else
688 distance=QuantumScale*fabs((double) (Sa*p[i]-Da*
689 GetPixelChannel(reconstruct_image,channel,q)));
690 channel_distortion[i]+=distance;
691 channel_distortion[CompositePixelChannel]+=distance;
692 }
693 area++;
694 p+=GetPixelChannels(image);
695 q+=GetPixelChannels(reconstruct_image);
696 }
697 #if defined(MAGICKCORE_OPENMP_SUPPORT)
698 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
699 #endif
700 for (j=0; j <= MaxPixelChannels; j++)
701 distortion[j]+=channel_distortion[j];
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
GetMeanErrorPerPixel(Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)712 static 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+=GetPixelChannels(image);
771 q+=GetPixelChannels(reconstruct_image);
772 continue;
773 }
774 Sa=QuantumScale*GetPixelAlpha(image,p);
775 Da=QuantumScale*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((double) (Sa*p[i]-Da*
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+=GetPixelChannels(image);
803 q+=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
GetMeanSquaredDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)815 static 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 ssize_t
829 j;
830
831 size_t
832 columns,
833 rows;
834
835 ssize_t
836 y;
837
838 status=MagickTrue;
839 rows=MagickMax(image->rows,reconstruct_image->rows);
840 columns=MagickMax(image->columns,reconstruct_image->columns);
841 area=0.0;
842 image_view=AcquireVirtualCacheView(image,exception);
843 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
844 #if defined(MAGICKCORE_OPENMP_SUPPORT)
845 #pragma omp parallel for schedule(static) shared(status) \
846 magick_number_threads(image,image,rows,1) reduction(+:area)
847 #endif
848 for (y=0; y < (ssize_t) rows; y++)
849 {
850 double
851 channel_distortion[MaxPixelChannels+1];
852
853 const Quantum
854 *magick_restrict p,
855 *magick_restrict q;
856
857 ssize_t
858 x;
859
860 if (status == MagickFalse)
861 continue;
862 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
863 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
864 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
865 {
866 status=MagickFalse;
867 continue;
868 }
869 (void) memset(channel_distortion,0,sizeof(channel_distortion));
870 for (x=0; x < (ssize_t) columns; x++)
871 {
872 double
873 Da,
874 Sa;
875
876 ssize_t
877 i;
878
879 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
880 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
881 {
882 p+=GetPixelChannels(image);
883 q+=GetPixelChannels(reconstruct_image);
884 continue;
885 }
886 Sa=QuantumScale*GetPixelAlpha(image,p);
887 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
888 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
889 {
890 double
891 distance;
892
893 PixelChannel channel = GetPixelChannelChannel(image,i);
894 PixelTrait traits = GetPixelChannelTraits(image,channel);
895 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
896 channel);
897 if ((traits == UndefinedPixelTrait) ||
898 (reconstruct_traits == UndefinedPixelTrait) ||
899 ((reconstruct_traits & UpdatePixelTrait) == 0))
900 continue;
901 if (channel == AlphaPixelChannel)
902 distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
903 channel,q));
904 else
905 distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
906 channel,q));
907 channel_distortion[i]+=distance*distance;
908 channel_distortion[CompositePixelChannel]+=distance*distance;
909 }
910 area++;
911 p+=GetPixelChannels(image);
912 q+=GetPixelChannels(reconstruct_image);
913 }
914 #if defined(MAGICKCORE_OPENMP_SUPPORT)
915 #pragma omp critical (MagickCore_GetMeanSquaredError)
916 #endif
917 for (j=0; j <= MaxPixelChannels; j++)
918 distortion[j]+=channel_distortion[j];
919 }
920 reconstruct_view=DestroyCacheView(reconstruct_view);
921 image_view=DestroyCacheView(image_view);
922 area=PerceptibleReciprocal(area);
923 for (j=0; j <= MaxPixelChannels; j++)
924 distortion[j]*=area;
925 distortion[CompositePixelChannel]/=GetImageChannels(image);
926 return(status);
927 }
928
GetNormalizedCrossCorrelationDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)929 static MagickBooleanType GetNormalizedCrossCorrelationDistortion(
930 const Image *image,const Image *reconstruct_image,double *distortion,
931 ExceptionInfo *exception)
932 {
933 #define SimilarityImageTag "Similarity/Image"
934
935 CacheView
936 *image_view,
937 *reconstruct_view;
938
939 ChannelStatistics
940 *image_statistics,
941 *reconstruct_statistics;
942
943 double
944 area;
945
946 MagickBooleanType
947 status;
948
949 MagickOffsetType
950 progress;
951
952 ssize_t
953 i;
954
955 size_t
956 columns,
957 rows;
958
959 ssize_t
960 y;
961
962 /*
963 Normalize to account for variation due to lighting and exposure condition.
964 */
965 image_statistics=GetImageStatistics(image,exception);
966 reconstruct_statistics=GetImageStatistics(reconstruct_image,exception);
967 if ((image_statistics == (ChannelStatistics *) NULL) ||
968 (reconstruct_statistics == (ChannelStatistics *) NULL))
969 {
970 if (image_statistics != (ChannelStatistics *) NULL)
971 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
972 image_statistics);
973 if (reconstruct_statistics != (ChannelStatistics *) NULL)
974 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
975 reconstruct_statistics);
976 return(MagickFalse);
977 }
978 status=MagickTrue;
979 progress=0;
980 for (i=0; i <= MaxPixelChannels; i++)
981 distortion[i]=0.0;
982 rows=MagickMax(image->rows,reconstruct_image->rows);
983 columns=MagickMax(image->columns,reconstruct_image->columns);
984 area=0.0;
985 image_view=AcquireVirtualCacheView(image,exception);
986 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
987 for (y=0; y < (ssize_t) rows; y++)
988 {
989 const Quantum
990 *magick_restrict p,
991 *magick_restrict q;
992
993 ssize_t
994 x;
995
996 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
997 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
998 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
999 {
1000 status=MagickFalse;
1001 break;
1002 }
1003 for (x=0; x < (ssize_t) columns; x++)
1004 {
1005 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1006 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1007 {
1008 p+=GetPixelChannels(image);
1009 q+=GetPixelChannels(reconstruct_image);
1010 continue;
1011 }
1012 area++;
1013 p+=GetPixelChannels(image);
1014 q+=GetPixelChannels(reconstruct_image);
1015 }
1016 }
1017 area=PerceptibleReciprocal(area);
1018 for (y=0; y < (ssize_t) rows; y++)
1019 {
1020 const Quantum
1021 *magick_restrict p,
1022 *magick_restrict q;
1023
1024 ssize_t
1025 x;
1026
1027 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1028 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1029 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1030 {
1031 status=MagickFalse;
1032 break;
1033 }
1034 for (x=0; x < (ssize_t) columns; x++)
1035 {
1036 double
1037 Da,
1038 Sa;
1039
1040 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1041 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1042 {
1043 p+=GetPixelChannels(image);
1044 q+=GetPixelChannels(reconstruct_image);
1045 continue;
1046 }
1047 Sa=QuantumScale*GetPixelAlpha(image,p);
1048 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1049 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1050 {
1051 PixelChannel channel = GetPixelChannelChannel(image,i);
1052 PixelTrait traits = GetPixelChannelTraits(image,channel);
1053 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1054 channel);
1055 if ((traits == UndefinedPixelTrait) ||
1056 (reconstruct_traits == UndefinedPixelTrait) ||
1057 ((reconstruct_traits & UpdatePixelTrait) == 0))
1058 continue;
1059 if (channel == AlphaPixelChannel)
1060 {
1061 distortion[i]+=area*QuantumScale*(p[i]-
1062 image_statistics[channel].mean)*(GetPixelChannel(
1063 reconstruct_image,channel,q)-
1064 reconstruct_statistics[channel].mean);
1065 }
1066 else
1067 {
1068 distortion[i]+=area*QuantumScale*(Sa*p[i]-
1069 image_statistics[channel].mean)*(Da*GetPixelChannel(
1070 reconstruct_image,channel,q)-
1071 reconstruct_statistics[channel].mean);
1072 }
1073 }
1074 p+=GetPixelChannels(image);
1075 q+=GetPixelChannels(reconstruct_image);
1076 }
1077 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1078 {
1079 MagickBooleanType
1080 proceed;
1081
1082 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1083 #pragma omp atomic
1084 #endif
1085 progress++;
1086 proceed=SetImageProgress(image,SimilarityImageTag,progress,rows);
1087 if (proceed == MagickFalse)
1088 {
1089 status=MagickFalse;
1090 break;
1091 }
1092 }
1093 }
1094 reconstruct_view=DestroyCacheView(reconstruct_view);
1095 image_view=DestroyCacheView(image_view);
1096 /*
1097 Divide by the standard deviation.
1098 */
1099 distortion[CompositePixelChannel]=0.0;
1100 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1101 {
1102 double
1103 gamma;
1104
1105 PixelChannel channel = GetPixelChannelChannel(image,i);
1106 gamma=image_statistics[channel].standard_deviation*
1107 reconstruct_statistics[channel].standard_deviation;
1108 gamma=PerceptibleReciprocal(gamma);
1109 distortion[i]=QuantumRange*gamma*distortion[i];
1110 distortion[CompositePixelChannel]+=distortion[i]*distortion[i];
1111 }
1112 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]/
1113 GetImageChannels(image));
1114 /*
1115 Free resources.
1116 */
1117 reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1118 reconstruct_statistics);
1119 image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
1120 image_statistics);
1121 return(status);
1122 }
1123
GetPeakAbsoluteDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)1124 static MagickBooleanType GetPeakAbsoluteDistortion(const Image *image,
1125 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1126 {
1127 CacheView
1128 *image_view,
1129 *reconstruct_view;
1130
1131 MagickBooleanType
1132 status;
1133
1134 size_t
1135 columns,
1136 rows;
1137
1138 ssize_t
1139 y;
1140
1141 status=MagickTrue;
1142 rows=MagickMax(image->rows,reconstruct_image->rows);
1143 columns=MagickMax(image->columns,reconstruct_image->columns);
1144 image_view=AcquireVirtualCacheView(image,exception);
1145 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1146 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1147 #pragma omp parallel for schedule(static) shared(status) \
1148 magick_number_threads(image,image,rows,1)
1149 #endif
1150 for (y=0; y < (ssize_t) rows; y++)
1151 {
1152 double
1153 channel_distortion[MaxPixelChannels+1];
1154
1155 const Quantum
1156 *magick_restrict p,
1157 *magick_restrict q;
1158
1159 ssize_t
1160 j,
1161 x;
1162
1163 if (status == MagickFalse)
1164 continue;
1165 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1166 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1167 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1168 {
1169 status=MagickFalse;
1170 continue;
1171 }
1172 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1173 for (x=0; x < (ssize_t) columns; x++)
1174 {
1175 double
1176 Da,
1177 Sa;
1178
1179 ssize_t
1180 i;
1181
1182 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1183 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1184 {
1185 p+=GetPixelChannels(image);
1186 q+=GetPixelChannels(reconstruct_image);
1187 continue;
1188 }
1189 Sa=QuantumScale*GetPixelAlpha(image,p);
1190 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
1191 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1192 {
1193 double
1194 distance;
1195
1196 PixelChannel channel = GetPixelChannelChannel(image,i);
1197 PixelTrait traits = GetPixelChannelTraits(image,channel);
1198 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1199 channel);
1200 if ((traits == UndefinedPixelTrait) ||
1201 (reconstruct_traits == UndefinedPixelTrait) ||
1202 ((reconstruct_traits & UpdatePixelTrait) == 0))
1203 continue;
1204 if (channel == AlphaPixelChannel)
1205 distance=QuantumScale*fabs((double) (p[i]-(double)
1206 GetPixelChannel(reconstruct_image,channel,q)));
1207 else
1208 distance=QuantumScale*fabs((double) (Sa*p[i]-Da*
1209 GetPixelChannel(reconstruct_image,channel,q)));
1210 if (distance > channel_distortion[i])
1211 channel_distortion[i]=distance;
1212 if (distance > channel_distortion[CompositePixelChannel])
1213 channel_distortion[CompositePixelChannel]=distance;
1214 }
1215 p+=GetPixelChannels(image);
1216 q+=GetPixelChannels(reconstruct_image);
1217 }
1218 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1219 #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1220 #endif
1221 for (j=0; j <= MaxPixelChannels; j++)
1222 if (channel_distortion[j] > distortion[j])
1223 distortion[j]=channel_distortion[j];
1224 }
1225 reconstruct_view=DestroyCacheView(reconstruct_view);
1226 image_view=DestroyCacheView(image_view);
1227 return(status);
1228 }
1229
MagickLog10(const double x)1230 static inline double MagickLog10(const double x)
1231 {
1232 #define Log10Epsilon (1.0e-11)
1233
1234 if (fabs(x) < Log10Epsilon)
1235 return(log10(Log10Epsilon));
1236 return(log10(fabs(x)));
1237 }
1238
GetPeakSignalToNoiseRatio(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)1239 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1240 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1241 {
1242 MagickBooleanType
1243 status;
1244
1245 ssize_t
1246 i;
1247
1248 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1249 for (i=0; i <= MaxPixelChannels; i++)
1250 if (fabs(distortion[i]) < MagickEpsilon)
1251 distortion[i]=INFINITY;
1252 else
1253 distortion[i]=10.0*MagickLog10(1.0)-10.0*MagickLog10(distortion[i]);
1254 return(status);
1255 }
1256
GetPerceptualHashDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)1257 static MagickBooleanType GetPerceptualHashDistortion(const Image *image,
1258 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1259 {
1260 ChannelPerceptualHash
1261 *channel_phash,
1262 *reconstruct_phash;
1263
1264 const char
1265 *artifact;
1266
1267 MagickBooleanType
1268 normalize;
1269
1270 ssize_t
1271 channel;
1272
1273 /*
1274 Compute perceptual hash in the sRGB colorspace.
1275 */
1276 channel_phash=GetImagePerceptualHash(image,exception);
1277 if (channel_phash == (ChannelPerceptualHash *) NULL)
1278 return(MagickFalse);
1279 reconstruct_phash=GetImagePerceptualHash(reconstruct_image,exception);
1280 if (reconstruct_phash == (ChannelPerceptualHash *) NULL)
1281 {
1282 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1283 channel_phash);
1284 return(MagickFalse);
1285 }
1286 artifact=GetImageArtifact(image,"phash:normalize");
1287 normalize=(artifact == (const char *) NULL) ||
1288 (IsStringTrue(artifact) == MagickFalse) ? MagickFalse : MagickTrue;
1289 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1290 #pragma omp parallel for schedule(static)
1291 #endif
1292 for (channel=0; channel < MaxPixelChannels; channel++)
1293 {
1294 double
1295 difference;
1296
1297 ssize_t
1298 i;
1299
1300 difference=0.0;
1301 for (i=0; i < MaximumNumberOfImageMoments; i++)
1302 {
1303 double
1304 alpha,
1305 beta;
1306
1307 ssize_t
1308 j;
1309
1310 for (j=0; j < (ssize_t) channel_phash[0].number_colorspaces; j++)
1311 {
1312 alpha=channel_phash[channel].phash[j][i];
1313 beta=reconstruct_phash[channel].phash[j][i];
1314 if (normalize == MagickFalse)
1315 difference+=(beta-alpha)*(beta-alpha);
1316 else
1317 difference=sqrt((beta-alpha)*(beta-alpha)/
1318 channel_phash[0].number_channels);
1319 }
1320 }
1321 distortion[channel]+=difference;
1322 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1323 #pragma omp critical (MagickCore_GetPerceptualHashDistortion)
1324 #endif
1325 distortion[CompositePixelChannel]+=difference;
1326 }
1327 /*
1328 Free resources.
1329 */
1330 reconstruct_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(
1331 reconstruct_phash);
1332 channel_phash=(ChannelPerceptualHash *) RelinquishMagickMemory(channel_phash);
1333 return(MagickTrue);
1334 }
1335
GetRootMeanSquaredDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)1336 static MagickBooleanType GetRootMeanSquaredDistortion(const Image *image,
1337 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1338 {
1339 MagickBooleanType
1340 status;
1341
1342 ssize_t
1343 i;
1344
1345 status=GetMeanSquaredDistortion(image,reconstruct_image,distortion,exception);
1346 for (i=0; i <= MaxPixelChannels; i++)
1347 distortion[i]=sqrt(distortion[i]);
1348 return(status);
1349 }
1350
GetStructuralSimilarityDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)1351 static MagickBooleanType GetStructuralSimilarityDistortion(const Image *image,
1352 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1353 {
1354 #define SSIMRadius 5.0
1355 #define SSIMSigma 1.5
1356 #define SSIMBlocksize 8
1357 #define SSIMK1 0.01
1358 #define SSIMK2 0.03
1359 #define SSIML 1.0
1360
1361 CacheView
1362 *image_view,
1363 *reconstruct_view;
1364
1365 char
1366 geometry[MagickPathExtent];
1367
1368 const char
1369 *artifact;
1370
1371 double
1372 area,
1373 c1,
1374 c2,
1375 radius,
1376 sigma;
1377
1378 KernelInfo
1379 *kernel_info;
1380
1381 MagickBooleanType
1382 status;
1383
1384 ssize_t
1385 i;
1386
1387 size_t
1388 columns,
1389 rows;
1390
1391 ssize_t
1392 y;
1393
1394 /*
1395 Compute structural similarity index @
1396 https://en.wikipedia.org/wiki/Structural_similarity.
1397 */
1398 radius=SSIMRadius;
1399 artifact=GetImageArtifact(image,"compare:ssim-radius");
1400 if (artifact != (const char *) NULL)
1401 radius=StringToDouble(artifact,(char **) NULL);
1402 sigma=SSIMSigma;
1403 artifact=GetImageArtifact(image,"compare:ssim-sigma");
1404 if (artifact != (const char *) NULL)
1405 sigma=StringToDouble(artifact,(char **) NULL);
1406 (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1407 radius,sigma);
1408 kernel_info=AcquireKernelInfo(geometry,exception);
1409 if (kernel_info == (KernelInfo *) NULL)
1410 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1411 image->filename);
1412 c1=pow(SSIMK1*SSIML,2.0);
1413 artifact=GetImageArtifact(image,"compare:ssim-k1");
1414 if (artifact != (const char *) NULL)
1415 c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1416 c2=pow(SSIMK2*SSIML,2.0);
1417 artifact=GetImageArtifact(image,"compare:ssim-k2");
1418 if (artifact != (const char *) NULL)
1419 c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1420 status=MagickTrue;
1421 area=0.0;
1422 rows=MagickMax(image->rows,reconstruct_image->rows);
1423 columns=MagickMax(image->columns,reconstruct_image->columns);
1424 image_view=AcquireVirtualCacheView(image,exception);
1425 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1426 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1427 #pragma omp parallel for schedule(static) shared(status) \
1428 magick_number_threads(image,reconstruct_image,rows,1)
1429 #endif
1430 for (y=0; y < (ssize_t) rows; y++)
1431 {
1432 double
1433 channel_distortion[MaxPixelChannels+1];
1434
1435 const Quantum
1436 *magick_restrict p,
1437 *magick_restrict q;
1438
1439 ssize_t
1440 i,
1441 x;
1442
1443 if (status == MagickFalse)
1444 continue;
1445 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1446 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1447 kernel_info->height,exception);
1448 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1449 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1450 kernel_info->height,exception);
1451 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1452 {
1453 status=MagickFalse;
1454 continue;
1455 }
1456 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1457 for (x=0; x < (ssize_t) columns; x++)
1458 {
1459 double
1460 x_pixel_mu[MaxPixelChannels+1],
1461 x_pixel_sigma_squared[MaxPixelChannels+1],
1462 xy_sigma[MaxPixelChannels+1],
1463 y_pixel_mu[MaxPixelChannels+1],
1464 y_pixel_sigma_squared[MaxPixelChannels+1];
1465
1466 const Quantum
1467 *magick_restrict reference,
1468 *magick_restrict target;
1469
1470 MagickRealType
1471 *k;
1472
1473 ssize_t
1474 v;
1475
1476 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
1477 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
1478 {
1479 p+=GetPixelChannels(image);
1480 q+=GetPixelChannels(reconstruct_image);
1481 continue;
1482 }
1483 (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1484 (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1485 (void) memset(xy_sigma,0,sizeof(xy_sigma));
1486 (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1487 (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1488 (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1489 k=kernel_info->values;
1490 reference=p;
1491 target=q;
1492 for (v=0; v < (ssize_t) kernel_info->height; v++)
1493 {
1494 ssize_t
1495 u;
1496
1497 for (u=0; u < (ssize_t) kernel_info->width; u++)
1498 {
1499 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1500 {
1501 double
1502 x_pixel,
1503 y_pixel;
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=QuantumScale*reference[i];
1514 x_pixel_mu[i]+=(*k)*x_pixel;
1515 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1516 y_pixel=QuantumScale*
1517 GetPixelChannel(reconstruct_image,channel,target);
1518 y_pixel_mu[i]+=(*k)*y_pixel;
1519 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1520 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1521 }
1522 k++;
1523 reference+=GetPixelChannels(image);
1524 target+=GetPixelChannels(reconstruct_image);
1525 }
1526 reference+=GetPixelChannels(image)*columns;
1527 target+=GetPixelChannels(reconstruct_image)*columns;
1528 }
1529 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1530 {
1531 double
1532 ssim,
1533 x_pixel_mu_squared,
1534 x_pixel_sigmas_squared,
1535 xy_mu,
1536 xy_sigmas,
1537 y_pixel_mu_squared,
1538 y_pixel_sigmas_squared;
1539
1540 PixelChannel channel = GetPixelChannelChannel(image,i);
1541 PixelTrait traits = GetPixelChannelTraits(image,channel);
1542 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1543 reconstruct_image,channel);
1544 if ((traits == UndefinedPixelTrait) ||
1545 (reconstruct_traits == UndefinedPixelTrait) ||
1546 ((reconstruct_traits & UpdatePixelTrait) == 0))
1547 continue;
1548 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1549 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1550 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1551 xy_sigmas=xy_sigma[i]-xy_mu;
1552 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1553 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1554 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1555 ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1556 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1557 channel_distortion[i]+=ssim;
1558 channel_distortion[CompositePixelChannel]+=ssim;
1559 }
1560 area++;
1561 p+=GetPixelChannels(image);
1562 q+=GetPixelChannels(reconstruct_image);
1563 }
1564 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1565 #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1566 #endif
1567 for (i=0; i <= MaxPixelChannels; i++)
1568 distortion[i]+=channel_distortion[i];
1569 }
1570 image_view=DestroyCacheView(image_view);
1571 reconstruct_view=DestroyCacheView(reconstruct_view);
1572 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1573 {
1574 PixelChannel channel = GetPixelChannelChannel(image,i);
1575 PixelTrait traits = GetPixelChannelTraits(image,channel);
1576 if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1577 continue;
1578 distortion[i]/=area;
1579 }
1580 distortion[CompositePixelChannel]/=area;
1581 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1582 kernel_info=DestroyKernelInfo(kernel_info);
1583 return(status);
1584 }
1585
GetStructuralDisimilarityDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)1586 static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image,
1587 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1588 {
1589 MagickBooleanType
1590 status;
1591
1592 ssize_t
1593 i;
1594
1595 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1596 distortion,exception);
1597 for (i=0; i <= MaxPixelChannels; i++)
1598 distortion[i]=(1.0-(distortion[i]))/2.0;
1599 return(status);
1600 }
1601
GetImageDistortion(Image * image,const Image * reconstruct_image,const MetricType metric,double * distortion,ExceptionInfo * exception)1602 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1603 const Image *reconstruct_image,const MetricType metric,double *distortion,
1604 ExceptionInfo *exception)
1605 {
1606 double
1607 *channel_distortion;
1608
1609 MagickBooleanType
1610 status;
1611
1612 size_t
1613 length;
1614
1615 assert(image != (Image *) NULL);
1616 assert(image->signature == MagickCoreSignature);
1617 if (image->debug != MagickFalse)
1618 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1619 assert(reconstruct_image != (const Image *) NULL);
1620 assert(reconstruct_image->signature == MagickCoreSignature);
1621 assert(distortion != (double *) NULL);
1622 *distortion=0.0;
1623 if (image->debug != MagickFalse)
1624 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1625 /*
1626 Get image distortion.
1627 */
1628 length=MaxPixelChannels+1UL;
1629 channel_distortion=(double *) AcquireQuantumMemory(length,
1630 sizeof(*channel_distortion));
1631 if (channel_distortion == (double *) NULL)
1632 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1633 (void) memset(channel_distortion,0,length*
1634 sizeof(*channel_distortion));
1635 switch (metric)
1636 {
1637 case AbsoluteErrorMetric:
1638 {
1639 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1640 exception);
1641 break;
1642 }
1643 case FuzzErrorMetric:
1644 {
1645 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1646 exception);
1647 break;
1648 }
1649 case MeanAbsoluteErrorMetric:
1650 {
1651 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1652 channel_distortion,exception);
1653 break;
1654 }
1655 case MeanErrorPerPixelErrorMetric:
1656 {
1657 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1658 exception);
1659 break;
1660 }
1661 case MeanSquaredErrorMetric:
1662 {
1663 status=GetMeanSquaredDistortion(image,reconstruct_image,
1664 channel_distortion,exception);
1665 break;
1666 }
1667 case NormalizedCrossCorrelationErrorMetric:
1668 default:
1669 {
1670 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1671 channel_distortion,exception);
1672 break;
1673 }
1674 case PeakAbsoluteErrorMetric:
1675 {
1676 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1677 channel_distortion,exception);
1678 break;
1679 }
1680 case PeakSignalToNoiseRatioErrorMetric:
1681 {
1682 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1683 channel_distortion,exception);
1684 break;
1685 }
1686 case PerceptualHashErrorMetric:
1687 {
1688 status=GetPerceptualHashDistortion(image,reconstruct_image,
1689 channel_distortion,exception);
1690 break;
1691 }
1692 case RootMeanSquaredErrorMetric:
1693 {
1694 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1695 channel_distortion,exception);
1696 break;
1697 }
1698 case StructuralSimilarityErrorMetric:
1699 {
1700 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1701 channel_distortion,exception);
1702 break;
1703 }
1704 case StructuralDissimilarityErrorMetric:
1705 {
1706 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1707 channel_distortion,exception);
1708 break;
1709 }
1710 }
1711 *distortion=channel_distortion[CompositePixelChannel];
1712 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1713 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1714 *distortion);
1715 return(status);
1716 }
1717
1718 /*
1719 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1720 % %
1721 % %
1722 % %
1723 % G e t I m a g e D i s t o r t i o n s %
1724 % %
1725 % %
1726 % %
1727 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1728 %
1729 % GetImageDistortions() compares the pixel channels of an image to a
1730 % reconstructed image and returns the specified distortion metric for each
1731 % channel.
1732 %
1733 % The format of the GetImageDistortions method is:
1734 %
1735 % double *GetImageDistortions(const Image *image,
1736 % const Image *reconstruct_image,const MetricType metric,
1737 % ExceptionInfo *exception)
1738 %
1739 % A description of each parameter follows:
1740 %
1741 % o image: the image.
1742 %
1743 % o reconstruct_image: the reconstruct image.
1744 %
1745 % o metric: the metric.
1746 %
1747 % o exception: return any errors or warnings in this structure.
1748 %
1749 */
GetImageDistortions(Image * image,const Image * reconstruct_image,const MetricType metric,ExceptionInfo * exception)1750 MagickExport double *GetImageDistortions(Image *image,
1751 const Image *reconstruct_image,const MetricType metric,
1752 ExceptionInfo *exception)
1753 {
1754 double
1755 *channel_distortion;
1756
1757 MagickBooleanType
1758 status;
1759
1760 size_t
1761 length;
1762
1763 assert(image != (Image *) NULL);
1764 assert(image->signature == MagickCoreSignature);
1765 if (image->debug != MagickFalse)
1766 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1767 assert(reconstruct_image != (const Image *) NULL);
1768 assert(reconstruct_image->signature == MagickCoreSignature);
1769 if (image->debug != MagickFalse)
1770 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1771 /*
1772 Get image distortion.
1773 */
1774 length=MaxPixelChannels+1UL;
1775 channel_distortion=(double *) AcquireQuantumMemory(length,
1776 sizeof(*channel_distortion));
1777 if (channel_distortion == (double *) NULL)
1778 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1779 (void) memset(channel_distortion,0,length*
1780 sizeof(*channel_distortion));
1781 status=MagickTrue;
1782 switch (metric)
1783 {
1784 case AbsoluteErrorMetric:
1785 {
1786 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1787 exception);
1788 break;
1789 }
1790 case FuzzErrorMetric:
1791 {
1792 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1793 exception);
1794 break;
1795 }
1796 case MeanAbsoluteErrorMetric:
1797 {
1798 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1799 channel_distortion,exception);
1800 break;
1801 }
1802 case MeanErrorPerPixelErrorMetric:
1803 {
1804 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1805 exception);
1806 break;
1807 }
1808 case MeanSquaredErrorMetric:
1809 {
1810 status=GetMeanSquaredDistortion(image,reconstruct_image,
1811 channel_distortion,exception);
1812 break;
1813 }
1814 case NormalizedCrossCorrelationErrorMetric:
1815 default:
1816 {
1817 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1818 channel_distortion,exception);
1819 break;
1820 }
1821 case PeakAbsoluteErrorMetric:
1822 {
1823 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1824 channel_distortion,exception);
1825 break;
1826 }
1827 case PeakSignalToNoiseRatioErrorMetric:
1828 {
1829 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1830 channel_distortion,exception);
1831 break;
1832 }
1833 case PerceptualHashErrorMetric:
1834 {
1835 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1836 channel_distortion,exception);
1837 break;
1838 }
1839 case RootMeanSquaredErrorMetric:
1840 {
1841 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1842 channel_distortion,exception);
1843 break;
1844 }
1845 case StructuralSimilarityErrorMetric:
1846 {
1847 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1848 channel_distortion,exception);
1849 break;
1850 }
1851 case StructuralDissimilarityErrorMetric:
1852 {
1853 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1854 channel_distortion,exception);
1855 break;
1856 }
1857 }
1858 if (status == MagickFalse)
1859 {
1860 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1861 return((double *) NULL);
1862 }
1863 return(channel_distortion);
1864 }
1865
1866 /*
1867 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1868 % %
1869 % %
1870 % %
1871 % I s I m a g e s E q u a l %
1872 % %
1873 % %
1874 % %
1875 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1876 %
1877 % IsImagesEqual() compare the pixels of two images and returns immediately
1878 % if any pixel is not identical.
1879 %
1880 % The format of the IsImagesEqual method is:
1881 %
1882 % MagickBooleanType IsImagesEqual(const Image *image,
1883 % const Image *reconstruct_image,ExceptionInfo *exception)
1884 %
1885 % A description of each parameter follows.
1886 %
1887 % o image: the image.
1888 %
1889 % o reconstruct_image: the reconstruct image.
1890 %
1891 % o exception: return any errors or warnings in this structure.
1892 %
1893 */
IsImagesEqual(const Image * image,const Image * reconstruct_image,ExceptionInfo * exception)1894 MagickExport MagickBooleanType IsImagesEqual(const Image *image,
1895 const Image *reconstruct_image,ExceptionInfo *exception)
1896 {
1897 CacheView
1898 *image_view,
1899 *reconstruct_view;
1900
1901 size_t
1902 columns,
1903 rows;
1904
1905 ssize_t
1906 y;
1907
1908 assert(image != (Image *) NULL);
1909 assert(image->signature == MagickCoreSignature);
1910 assert(reconstruct_image != (const Image *) NULL);
1911 assert(reconstruct_image->signature == MagickCoreSignature);
1912 rows=MagickMax(image->rows,reconstruct_image->rows);
1913 columns=MagickMax(image->columns,reconstruct_image->columns);
1914 image_view=AcquireVirtualCacheView(image,exception);
1915 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1916 for (y=0; y < (ssize_t) rows; y++)
1917 {
1918 const Quantum
1919 *magick_restrict p,
1920 *magick_restrict q;
1921
1922 ssize_t
1923 x;
1924
1925 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1926 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1927 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1928 break;
1929 for (x=0; x < (ssize_t) columns; x++)
1930 {
1931 ssize_t
1932 i;
1933
1934 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1935 {
1936 double
1937 distance;
1938
1939 PixelChannel channel = GetPixelChannelChannel(image,i);
1940 PixelTrait traits = GetPixelChannelTraits(image,channel);
1941 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1942 channel);
1943 if ((traits == UndefinedPixelTrait) ||
1944 (reconstruct_traits == UndefinedPixelTrait) ||
1945 ((reconstruct_traits & UpdatePixelTrait) == 0))
1946 continue;
1947 distance=fabs((double) (p[i]-(double) GetPixelChannel(reconstruct_image,
1948 channel,q)));
1949 if (distance >= MagickEpsilon)
1950 break;
1951 }
1952 if (i < (ssize_t) GetPixelChannels(image))
1953 break;
1954 p+=GetPixelChannels(image);
1955 q+=GetPixelChannels(reconstruct_image);
1956 }
1957 if (x < (ssize_t) columns)
1958 break;
1959 }
1960 reconstruct_view=DestroyCacheView(reconstruct_view);
1961 image_view=DestroyCacheView(image_view);
1962 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1963 }
1964
1965 /*
1966 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1967 % %
1968 % %
1969 % %
1970 % S e t I m a g e C o l o r M e t r i c %
1971 % %
1972 % %
1973 % %
1974 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1975 %
1976 % SetImageColorMetric() measures the difference between colors at each pixel
1977 % location of two images. A value other than 0 means the colors match
1978 % exactly. Otherwise an error measure is computed by summing over all
1979 % pixels in an image the distance squared in RGB space between each image
1980 % pixel and its corresponding pixel in the reconstruct image. The error
1981 % measure is assigned to these image members:
1982 %
1983 % o mean_error_per_pixel: The mean error for any single pixel in
1984 % the image.
1985 %
1986 % o normalized_mean_error: The normalized mean quantization error for
1987 % any single pixel in the image. This distance measure is normalized to
1988 % a range between 0 and 1. It is independent of the range of red, green,
1989 % and blue values in the image.
1990 %
1991 % o normalized_maximum_error: The normalized maximum quantization
1992 % error for any single pixel in the image. This distance measure is
1993 % normalized to a range between 0 and 1. It is independent of the range
1994 % of red, green, and blue values in your image.
1995 %
1996 % A small normalized mean square error, accessed as
1997 % image->normalized_mean_error, suggests the images are very similar in
1998 % spatial layout and color.
1999 %
2000 % The format of the SetImageColorMetric method is:
2001 %
2002 % MagickBooleanType SetImageColorMetric(Image *image,
2003 % const Image *reconstruct_image,ExceptionInfo *exception)
2004 %
2005 % A description of each parameter follows.
2006 %
2007 % o image: the image.
2008 %
2009 % o reconstruct_image: the reconstruct image.
2010 %
2011 % o exception: return any errors or warnings in this structure.
2012 %
2013 */
SetImageColorMetric(Image * image,const Image * reconstruct_image,ExceptionInfo * exception)2014 MagickExport MagickBooleanType SetImageColorMetric(Image *image,
2015 const Image *reconstruct_image,ExceptionInfo *exception)
2016 {
2017 CacheView
2018 *image_view,
2019 *reconstruct_view;
2020
2021 double
2022 area,
2023 maximum_error,
2024 mean_error,
2025 mean_error_per_pixel;
2026
2027 MagickBooleanType
2028 status;
2029
2030 size_t
2031 columns,
2032 rows;
2033
2034 ssize_t
2035 y;
2036
2037 assert(image != (Image *) NULL);
2038 assert(image->signature == MagickCoreSignature);
2039 assert(reconstruct_image != (const Image *) NULL);
2040 assert(reconstruct_image->signature == MagickCoreSignature);
2041 area=0.0;
2042 maximum_error=0.0;
2043 mean_error_per_pixel=0.0;
2044 mean_error=0.0;
2045 rows=MagickMax(image->rows,reconstruct_image->rows);
2046 columns=MagickMax(image->columns,reconstruct_image->columns);
2047 image_view=AcquireVirtualCacheView(image,exception);
2048 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2049 for (y=0; y < (ssize_t) rows; y++)
2050 {
2051 const Quantum
2052 *magick_restrict p,
2053 *magick_restrict q;
2054
2055 ssize_t
2056 x;
2057
2058 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2059 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2060 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2061 break;
2062 for (x=0; x < (ssize_t) columns; x++)
2063 {
2064 ssize_t
2065 i;
2066
2067 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2068 {
2069 double
2070 distance;
2071
2072 PixelChannel channel = GetPixelChannelChannel(image,i);
2073 PixelTrait traits = GetPixelChannelTraits(image,channel);
2074 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2075 channel);
2076 if ((traits == UndefinedPixelTrait) ||
2077 (reconstruct_traits == UndefinedPixelTrait) ||
2078 ((reconstruct_traits & UpdatePixelTrait) == 0))
2079 continue;
2080 distance=fabs((double) (p[i]-(double) GetPixelChannel(reconstruct_image,
2081 channel,q)));
2082 if (distance >= MagickEpsilon)
2083 {
2084 mean_error_per_pixel+=distance;
2085 mean_error+=distance*distance;
2086 if (distance > maximum_error)
2087 maximum_error=distance;
2088 }
2089 area++;
2090 }
2091 p+=GetPixelChannels(image);
2092 q+=GetPixelChannels(reconstruct_image);
2093 }
2094 }
2095 reconstruct_view=DestroyCacheView(reconstruct_view);
2096 image_view=DestroyCacheView(image_view);
2097 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2098 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
2099 mean_error/area);
2100 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2101 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2102 return(status);
2103 }
2104
2105 /*
2106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2107 % %
2108 % %
2109 % %
2110 % S i m i l a r i t y I m a g e %
2111 % %
2112 % %
2113 % %
2114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2115 %
2116 % SimilarityImage() compares the reference image of the image and returns the
2117 % best match offset. In addition, it returns a similarity image such that an
2118 % exact match location is completely white and if none of the pixels match,
2119 % black, otherwise some gray level in-between.
2120 %
2121 % The format of the SimilarityImageImage method is:
2122 %
2123 % Image *SimilarityImage(const Image *image,const Image *reference,
2124 % const MetricType metric,const double similarity_threshold,
2125 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2126 %
2127 % A description of each parameter follows:
2128 %
2129 % o image: the image.
2130 %
2131 % o reference: find an area of the image that closely resembles this image.
2132 %
2133 % o metric: the metric.
2134 %
2135 % o similarity_threshold: minimum distortion for (sub)image match.
2136 %
2137 % o offset: the best match offset of the reference image within the image.
2138 %
2139 % o similarity: the computed similarity between the images.
2140 %
2141 % o exception: return any errors or warnings in this structure.
2142 %
2143 */
2144
GetSimilarityMetric(const Image * image,const Image * reference,const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo * exception)2145 static double GetSimilarityMetric(const Image *image,const Image *reference,
2146 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2147 ExceptionInfo *exception)
2148 {
2149 double
2150 distortion;
2151
2152 Image
2153 *similarity_image;
2154
2155 MagickBooleanType
2156 status;
2157
2158 RectangleInfo
2159 geometry;
2160
2161 SetGeometry(reference,&geometry);
2162 geometry.x=x_offset;
2163 geometry.y=y_offset;
2164 similarity_image=CropImage(image,&geometry,exception);
2165 if (similarity_image == (Image *) NULL)
2166 return(0.0);
2167 distortion=0.0;
2168 status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2169 exception);
2170 similarity_image=DestroyImage(similarity_image);
2171 if (status == MagickFalse)
2172 return(0.0);
2173 return(distortion);
2174 }
2175
SimilarityImage(const Image * image,const Image * reference,const MetricType metric,const double similarity_threshold,RectangleInfo * offset,double * similarity_metric,ExceptionInfo * exception)2176 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2177 const MetricType metric,const double similarity_threshold,
2178 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2179 {
2180 #define SimilarityImageTag "Similarity/Image"
2181
2182 CacheView
2183 *similarity_view;
2184
2185 Image
2186 *similarity_image;
2187
2188 MagickBooleanType
2189 status;
2190
2191 MagickOffsetType
2192 progress;
2193
2194 ssize_t
2195 y;
2196
2197 assert(image != (const Image *) NULL);
2198 assert(image->signature == MagickCoreSignature);
2199 if (image->debug != MagickFalse)
2200 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2201 assert(exception != (ExceptionInfo *) NULL);
2202 assert(exception->signature == MagickCoreSignature);
2203 assert(offset != (RectangleInfo *) NULL);
2204 SetGeometry(reference,offset);
2205 *similarity_metric=MagickMaximumValue;
2206 similarity_image=CloneImage(image,image->columns-reference->columns+1,
2207 image->rows-reference->rows+1,MagickTrue,exception);
2208 if (similarity_image == (Image *) NULL)
2209 return((Image *) NULL);
2210 status=SetImageStorageClass(similarity_image,DirectClass,exception);
2211 if (status == MagickFalse)
2212 {
2213 similarity_image=DestroyImage(similarity_image);
2214 return((Image *) NULL);
2215 }
2216 (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2217 exception);
2218 /*
2219 Measure similarity of reference image against image.
2220 */
2221 status=MagickTrue;
2222 progress=0;
2223 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2224 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2225 #pragma omp parallel for schedule(static) \
2226 shared(progress,status,similarity_metric) \
2227 magick_number_threads(image,image,image->rows-reference->rows+1,1)
2228 #endif
2229 for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2230 {
2231 double
2232 similarity;
2233
2234 Quantum
2235 *magick_restrict q;
2236
2237 ssize_t
2238 x;
2239
2240 if (status == MagickFalse)
2241 continue;
2242 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2243 #pragma omp flush(similarity_metric)
2244 #endif
2245 if (*similarity_metric <= similarity_threshold)
2246 continue;
2247 q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2248 1,exception);
2249 if (q == (Quantum *) NULL)
2250 {
2251 status=MagickFalse;
2252 continue;
2253 }
2254 for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2255 {
2256 ssize_t
2257 i;
2258
2259 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2260 #pragma omp flush(similarity_metric)
2261 #endif
2262 if (*similarity_metric <= similarity_threshold)
2263 break;
2264 similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2265 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2266 #pragma omp critical (MagickCore_SimilarityImage)
2267 #endif
2268 if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2269 (metric == UndefinedErrorMetric))
2270 similarity=1.0-similarity;
2271 if (similarity < *similarity_metric)
2272 {
2273 offset->x=x;
2274 offset->y=y;
2275 *similarity_metric=similarity;
2276 }
2277 if (metric == PerceptualHashErrorMetric)
2278 similarity=MagickMin(0.01*similarity,1.0);
2279 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2280 {
2281 PixelChannel channel = GetPixelChannelChannel(image,i);
2282 PixelTrait traits = GetPixelChannelTraits(image,channel);
2283 PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2284 channel);
2285 if ((traits == UndefinedPixelTrait) ||
2286 (similarity_traits == UndefinedPixelTrait) ||
2287 ((similarity_traits & UpdatePixelTrait) == 0))
2288 continue;
2289 SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2290 QuantumRange*similarity),q);
2291 }
2292 q+=GetPixelChannels(similarity_image);
2293 }
2294 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2295 status=MagickFalse;
2296 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2297 {
2298 MagickBooleanType
2299 proceed;
2300
2301 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2302 #pragma omp atomic
2303 #endif
2304 progress++;
2305 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2306 if (proceed == MagickFalse)
2307 status=MagickFalse;
2308 }
2309 }
2310 similarity_view=DestroyCacheView(similarity_view);
2311 if (status == MagickFalse)
2312 similarity_image=DestroyImage(similarity_image);
2313 return(similarity_image);
2314 }
2315