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-2019 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 %
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 register 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 register const Quantum
238 *magick_restrict p,
239 *magick_restrict q;
240
241 register Quantum
242 *magick_restrict r;
243
244 register 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 register 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
287 PixelChannel channel = GetPixelChannelChannel(image,i);
288 PixelTrait traits = GetPixelChannelTraits(image,channel);
289 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
290 channel);
291 if ((traits == UndefinedPixelTrait) ||
292 (reconstruct_traits == UndefinedPixelTrait) ||
293 ((reconstruct_traits & UpdatePixelTrait) == 0))
294 continue;
295 if (channel == AlphaPixelChannel)
296 distance=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
297 else
298 distance=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
299 if ((distance*distance) > fuzz)
300 {
301 difference=MagickTrue;
302 break;
303 }
304 }
305 if (difference == MagickFalse)
306 SetPixelViaPixelInfo(highlight_image,&lowlight,r);
307 else
308 SetPixelViaPixelInfo(highlight_image,&highlight,r);
309 p+=GetPixelChannels(image);
310 q+=GetPixelChannels(reconstruct_image);
311 r+=GetPixelChannels(highlight_image);
312 }
313 sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
314 if (sync == MagickFalse)
315 status=MagickFalse;
316 }
317 highlight_view=DestroyCacheView(highlight_view);
318 reconstruct_view=DestroyCacheView(reconstruct_view);
319 image_view=DestroyCacheView(image_view);
320 (void) CompositeImage(difference_image,highlight_image,image->compose,
321 MagickTrue,0,0,exception);
322 highlight_image=DestroyImage(highlight_image);
323 if (status == MagickFalse)
324 difference_image=DestroyImage(difference_image);
325 return(difference_image);
326 }
327
328 /*
329 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
330 % %
331 % %
332 % %
333 % G e t I m a g e D i s t o r t i o n %
334 % %
335 % %
336 % %
337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
338 %
339 % GetImageDistortion() compares one or more pixel channels of an image to a
340 % reconstructed image and returns the specified distortion metric.
341 %
342 % The format of the GetImageDistortion method is:
343 %
344 % MagickBooleanType GetImageDistortion(const Image *image,
345 % const Image *reconstruct_image,const MetricType metric,
346 % double *distortion,ExceptionInfo *exception)
347 %
348 % A description of each parameter follows:
349 %
350 % o image: the image.
351 %
352 % o reconstruct_image: the reconstruct image.
353 %
354 % o metric: the metric.
355 %
356 % o distortion: the computed distortion between the images.
357 %
358 % o exception: return any errors or warnings in this structure.
359 %
360 */
361
GetAbsoluteDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)362 static MagickBooleanType GetAbsoluteDistortion(const Image *image,
363 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
364 {
365 CacheView
366 *image_view,
367 *reconstruct_view;
368
369 double
370 fuzz;
371
372 MagickBooleanType
373 status;
374
375 size_t
376 columns,
377 rows;
378
379 ssize_t
380 y;
381
382 /*
383 Compute the absolute difference in pixels between two images.
384 */
385 status=MagickTrue;
386 fuzz=(double) MagickMin(GetPixelChannels(image),
387 GetPixelChannels(reconstruct_image))*
388 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 register const Quantum
403 *magick_restrict p,
404 *magick_restrict q;
405
406 register 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 distance,
425 Sa;
426
427 MagickBooleanType
428 difference;
429
430 register ssize_t
431 i;
432
433 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
434 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
435 {
436 p+=GetPixelChannels(image);
437 q+=GetPixelChannels(reconstruct_image);
438 continue;
439 }
440 difference=MagickFalse;
441 distance=0.0;
442 Sa=QuantumScale*GetPixelAlpha(image,p);
443 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
444 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
445 {
446 double
447 pixel;
448
449 PixelChannel channel = GetPixelChannelChannel(image,i);
450 PixelTrait traits = GetPixelChannelTraits(image,channel);
451 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
452 channel);
453 if ((traits == UndefinedPixelTrait) ||
454 (reconstruct_traits == UndefinedPixelTrait) ||
455 ((reconstruct_traits & UpdatePixelTrait) == 0))
456 continue;
457 if (channel == AlphaPixelChannel)
458 pixel=(double) p[i]-GetPixelChannel(reconstruct_image,channel,q);
459 else
460 pixel=Sa*p[i]-Da*GetPixelChannel(reconstruct_image,channel,q);
461 distance+=pixel*pixel;
462 if (distance > fuzz)
463 {
464 channel_distortion[i]++;
465 difference=MagickTrue;
466 }
467 }
468 if (difference != MagickFalse)
469 channel_distortion[CompositePixelChannel]++;
470 p+=GetPixelChannels(image);
471 q+=GetPixelChannels(reconstruct_image);
472 }
473 #if defined(MAGICKCORE_OPENMP_SUPPORT)
474 #pragma omp critical (MagickCore_GetAbsoluteDistortion)
475 #endif
476 for (j=0; j <= MaxPixelChannels; j++)
477 distortion[j]+=channel_distortion[j];
478 }
479 reconstruct_view=DestroyCacheView(reconstruct_view);
480 image_view=DestroyCacheView(image_view);
481 return(status);
482 }
483
GetFuzzDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)484 static MagickBooleanType GetFuzzDistortion(const Image *image,
485 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
486 {
487 CacheView
488 *image_view,
489 *reconstruct_view;
490
491 double
492 area;
493
494 MagickBooleanType
495 status;
496
497 register ssize_t
498 j;
499
500 size_t
501 columns,
502 rows;
503
504 ssize_t
505 y;
506
507 status=MagickTrue;
508 rows=MagickMax(image->rows,reconstruct_image->rows);
509 columns=MagickMax(image->columns,reconstruct_image->columns);
510 area=0.0;
511 image_view=AcquireVirtualCacheView(image,exception);
512 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
513 #if defined(MAGICKCORE_OPENMP_SUPPORT)
514 #pragma omp parallel for schedule(static) shared(status) \
515 magick_number_threads(image,image,rows,1) reduction(+:area)
516 #endif
517 for (y=0; y < (ssize_t) rows; y++)
518 {
519 double
520 channel_distortion[MaxPixelChannels+1];
521
522 register const Quantum
523 *magick_restrict p,
524 *magick_restrict q;
525
526 register ssize_t
527 x;
528
529 if (status == MagickFalse)
530 continue;
531 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
532 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
533 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
534 {
535 status=MagickFalse;
536 continue;
537 }
538 (void) memset(channel_distortion,0,sizeof(channel_distortion));
539 for (x=0; x < (ssize_t) columns; x++)
540 {
541 double
542 Da,
543 Sa;
544
545 register ssize_t
546 i;
547
548 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
549 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
550 {
551 p+=GetPixelChannels(image);
552 q+=GetPixelChannels(reconstruct_image);
553 continue;
554 }
555 Sa=QuantumScale*GetPixelAlpha(image,p);
556 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
557 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
558 {
559 double
560 distance;
561
562 PixelChannel channel = GetPixelChannelChannel(image,i);
563 PixelTrait traits = GetPixelChannelTraits(image,channel);
564 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
565 channel);
566 if ((traits == UndefinedPixelTrait) ||
567 (reconstruct_traits == UndefinedPixelTrait) ||
568 ((reconstruct_traits & UpdatePixelTrait) == 0))
569 continue;
570 if (channel == AlphaPixelChannel)
571 distance=QuantumScale*(p[i]-GetPixelChannel(reconstruct_image,
572 channel,q));
573 else
574 distance=QuantumScale*(Sa*p[i]-Da*GetPixelChannel(reconstruct_image,
575 channel,q));
576 channel_distortion[i]+=distance*distance;
577 channel_distortion[CompositePixelChannel]+=distance*distance;
578 }
579 area++;
580 p+=GetPixelChannels(image);
581 q+=GetPixelChannels(reconstruct_image);
582 }
583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
584 #pragma omp critical (MagickCore_GetFuzzDistortion)
585 #endif
586 for (j=0; j <= MaxPixelChannels; j++)
587 distortion[j]+=channel_distortion[j];
588 }
589 reconstruct_view=DestroyCacheView(reconstruct_view);
590 image_view=DestroyCacheView(image_view);
591 area=PerceptibleReciprocal(area);
592 for (j=0; j <= MaxPixelChannels; j++)
593 distortion[j]*=area;
594 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
595 distortion[CompositePixelChannel]=sqrt(distortion[CompositePixelChannel]);
596 return(status);
597 }
598
GetMeanAbsoluteDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)599 static MagickBooleanType GetMeanAbsoluteDistortion(const Image *image,
600 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
601 {
602 CacheView
603 *image_view,
604 *reconstruct_view;
605
606 double
607 area;
608
609 MagickBooleanType
610 status;
611
612 register ssize_t
613 j;
614
615 size_t
616 columns,
617 rows;
618
619 ssize_t
620 y;
621
622 status=MagickTrue;
623 rows=MagickMax(image->rows,reconstruct_image->rows);
624 columns=MagickMax(image->columns,reconstruct_image->columns);
625 area=0.0;
626 image_view=AcquireVirtualCacheView(image,exception);
627 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
628 #if defined(MAGICKCORE_OPENMP_SUPPORT)
629 #pragma omp parallel for schedule(static) shared(status) \
630 magick_number_threads(image,image,rows,1) reduction(+:area)
631 #endif
632 for (y=0; y < (ssize_t) rows; y++)
633 {
634 double
635 channel_distortion[MaxPixelChannels+1];
636
637 register const Quantum
638 *magick_restrict p,
639 *magick_restrict q;
640
641 register ssize_t
642 x;
643
644 if (status == MagickFalse)
645 continue;
646 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
647 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
648 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
649 {
650 status=MagickFalse;
651 continue;
652 }
653 (void) memset(channel_distortion,0,sizeof(channel_distortion));
654 for (x=0; x < (ssize_t) columns; x++)
655 {
656 double
657 Da,
658 Sa;
659
660 register ssize_t
661 i;
662
663 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
664 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
665 {
666 p+=GetPixelChannels(image);
667 q+=GetPixelChannels(reconstruct_image);
668 continue;
669 }
670 Sa=QuantumScale*GetPixelAlpha(image,p);
671 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
672 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
673 {
674 double
675 distance;
676
677 PixelChannel channel = GetPixelChannelChannel(image,i);
678 PixelTrait traits = GetPixelChannelTraits(image,channel);
679 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
680 channel);
681 if ((traits == UndefinedPixelTrait) ||
682 (reconstruct_traits == UndefinedPixelTrait) ||
683 ((reconstruct_traits & UpdatePixelTrait) == 0))
684 continue;
685 if (channel == AlphaPixelChannel)
686 distance=QuantumScale*fabs((double) p[i]-
687 GetPixelChannel(reconstruct_image,channel,q));
688 else
689 distance=QuantumScale*fabs(Sa*p[i]-Da*
690 GetPixelChannel(reconstruct_image,channel,q));
691 channel_distortion[i]+=distance;
692 channel_distortion[CompositePixelChannel]+=distance;
693 }
694 area++;
695 p+=GetPixelChannels(image);
696 q+=GetPixelChannels(reconstruct_image);
697 }
698 #if defined(MAGICKCORE_OPENMP_SUPPORT)
699 #pragma omp critical (MagickCore_GetMeanAbsoluteError)
700 #endif
701 for (j=0; j <= MaxPixelChannels; j++)
702 distortion[j]+=channel_distortion[j];
703 }
704 reconstruct_view=DestroyCacheView(reconstruct_view);
705 image_view=DestroyCacheView(image_view);
706 area=PerceptibleReciprocal(area);
707 for (j=0; j <= MaxPixelChannels; j++)
708 distortion[j]*=area;
709 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
710 return(status);
711 }
712
GetMeanErrorPerPixel(Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)713 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
714 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
715 {
716 CacheView
717 *image_view,
718 *reconstruct_view;
719
720 MagickBooleanType
721 status;
722
723 double
724 area,
725 maximum_error,
726 mean_error;
727
728 size_t
729 columns,
730 rows;
731
732 ssize_t
733 y;
734
735 status=MagickTrue;
736 area=0.0;
737 maximum_error=0.0;
738 mean_error=0.0;
739 rows=MagickMax(image->rows,reconstruct_image->rows);
740 columns=MagickMax(image->columns,reconstruct_image->columns);
741 image_view=AcquireVirtualCacheView(image,exception);
742 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
743 for (y=0; y < (ssize_t) rows; y++)
744 {
745 register const Quantum
746 *magick_restrict p,
747 *magick_restrict q;
748
749 register ssize_t
750 x;
751
752 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
753 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
754 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
755 {
756 status=MagickFalse;
757 break;
758 }
759 for (x=0; x < (ssize_t) columns; x++)
760 {
761 double
762 Da,
763 Sa;
764
765 register ssize_t
766 i;
767
768 if ((GetPixelReadMask(image,p) <= (QuantumRange/2)) ||
769 (GetPixelReadMask(reconstruct_image,q) <= (QuantumRange/2)))
770 {
771 p+=GetPixelChannels(image);
772 q+=GetPixelChannels(reconstruct_image);
773 continue;
774 }
775 Sa=QuantumScale*GetPixelAlpha(image,p);
776 Da=QuantumScale*GetPixelAlpha(reconstruct_image,q);
777 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
778 {
779 double
780 distance;
781
782 PixelChannel channel = GetPixelChannelChannel(image,i);
783 PixelTrait traits = GetPixelChannelTraits(image,channel);
784 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
785 channel);
786 if ((traits == UndefinedPixelTrait) ||
787 (reconstruct_traits == UndefinedPixelTrait) ||
788 ((reconstruct_traits & UpdatePixelTrait) == 0))
789 continue;
790 if (channel == AlphaPixelChannel)
791 distance=fabs((double) p[i]-
792 GetPixelChannel(reconstruct_image,channel,q));
793 else
794 distance=fabs(Sa*p[i]-Da*
795 GetPixelChannel(reconstruct_image,channel,q));
796 distortion[i]+=distance;
797 distortion[CompositePixelChannel]+=distance;
798 mean_error+=distance*distance;
799 if (distance > maximum_error)
800 maximum_error=distance;
801 area++;
802 }
803 p+=GetPixelChannels(image);
804 q+=GetPixelChannels(reconstruct_image);
805 }
806 }
807 reconstruct_view=DestroyCacheView(reconstruct_view);
808 image_view=DestroyCacheView(image_view);
809 image->error.mean_error_per_pixel=distortion[CompositePixelChannel]/area;
810 image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
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 register 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 register const Quantum
854 *magick_restrict p,
855 *magick_restrict q;
856
857 register 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 register 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 register 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 register const Quantum
990 *magick_restrict p,
991 *magick_restrict q;
992
993 register 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 register const Quantum
1021 *magick_restrict p,
1022 *magick_restrict q;
1023
1024 register 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 register const Quantum
1156 *magick_restrict p,
1157 *magick_restrict q;
1158
1159 register 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 register 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]-
1206 GetPixelChannel(reconstruct_image,channel,q));
1207 else
1208 distance=QuantumScale*fabs(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 register 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 register 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 register 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 register 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 c1,
1373 c2,
1374 radius,
1375 sigma;
1376
1377 KernelInfo
1378 *kernel_info;
1379
1380 MagickBooleanType
1381 status;
1382
1383 register ssize_t
1384 i;
1385
1386 size_t
1387 columns,
1388 rows;
1389
1390 ssize_t
1391 y;
1392
1393 /*
1394 Compute structural similarity index @
1395 https://en.wikipedia.org/wiki/Structural_similarity.
1396 */
1397 radius=SSIMRadius;
1398 artifact=GetImageArtifact(image,"compare:ssim-radius");
1399 if (artifact != (const char *) NULL)
1400 radius=StringToDouble(artifact,(char **) NULL);
1401 sigma=SSIMSigma;
1402 artifact=GetImageArtifact(image,"compare:ssim-sigma");
1403 if (artifact != (const char *) NULL)
1404 sigma=StringToDouble(artifact,(char **) NULL);
1405 (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1406 radius,sigma);
1407 kernel_info=AcquireKernelInfo(geometry,exception);
1408 if (kernel_info == (KernelInfo *) NULL)
1409 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1410 image->filename);
1411 c1=pow(SSIMK1*SSIML,2.0);
1412 artifact=GetImageArtifact(image,"compare:ssim-k1");
1413 if (artifact != (const char *) NULL)
1414 c1=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1415 c2=pow(SSIMK2*SSIML,2.0);
1416 artifact=GetImageArtifact(image,"compare:ssim-k2");
1417 if (artifact != (const char *) NULL)
1418 c2=pow(StringToDouble(artifact,(char **) NULL)*SSIML,2.0);
1419 status=MagickTrue;
1420 rows=MagickMax(image->rows,reconstruct_image->rows);
1421 columns=MagickMax(image->columns,reconstruct_image->columns);
1422 image_view=AcquireVirtualCacheView(image,exception);
1423 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1424 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1425 #pragma omp parallel for schedule(static) shared(status) \
1426 magick_number_threads(image,reconstruct_image,rows,1)
1427 #endif
1428 for (y=0; y < (ssize_t) rows; y++)
1429 {
1430 double
1431 channel_distortion[MaxPixelChannels+1];
1432
1433 register const Quantum
1434 *magick_restrict p,
1435 *magick_restrict q;
1436
1437 register ssize_t
1438 i,
1439 x;
1440
1441 if (status == MagickFalse)
1442 continue;
1443 p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
1444 ((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1445 kernel_info->height,exception);
1446 q=GetCacheViewVirtualPixels(reconstruct_view,-((ssize_t) kernel_info->width/
1447 2L),y-((ssize_t) kernel_info->height/2L),columns+kernel_info->width,
1448 kernel_info->height,exception);
1449 if ((p == (const Quantum *) NULL) || (q == (const Quantum *) NULL))
1450 {
1451 status=MagickFalse;
1452 continue;
1453 }
1454 (void) memset(channel_distortion,0,sizeof(channel_distortion));
1455 for (x=0; x < (ssize_t) columns; x++)
1456 {
1457 double
1458 x_pixel_mu[MaxPixelChannels+1],
1459 x_pixel_sigma_squared[MaxPixelChannels+1],
1460 xy_sigma[MaxPixelChannels+1],
1461 y_pixel_mu[MaxPixelChannels+1],
1462 y_pixel_sigma_squared[MaxPixelChannels+1];
1463
1464 register const Quantum
1465 *magick_restrict reference,
1466 *magick_restrict target;
1467
1468 register MagickRealType
1469 *k;
1470
1471 ssize_t
1472 v;
1473
1474 (void) memset(x_pixel_mu,0,sizeof(x_pixel_mu));
1475 (void) memset(x_pixel_sigma_squared,0,sizeof(x_pixel_sigma_squared));
1476 (void) memset(xy_sigma,0,sizeof(xy_sigma));
1477 (void) memset(x_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1478 (void) memset(y_pixel_mu,0,sizeof(y_pixel_mu));
1479 (void) memset(y_pixel_sigma_squared,0,sizeof(y_pixel_sigma_squared));
1480 k=kernel_info->values;
1481 reference=p;
1482 target=q;
1483 for (v=0; v < (ssize_t) kernel_info->height; v++)
1484 {
1485 register ssize_t
1486 u;
1487
1488 for (u=0; u < (ssize_t) kernel_info->width; u++)
1489 {
1490 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1491 {
1492 double
1493 x_pixel,
1494 y_pixel;
1495
1496 PixelChannel channel = GetPixelChannelChannel(image,i);
1497 PixelTrait traits = GetPixelChannelTraits(image,channel);
1498 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1499 reconstruct_image,channel);
1500 if ((traits == UndefinedPixelTrait) ||
1501 (reconstruct_traits == UndefinedPixelTrait) ||
1502 ((reconstruct_traits & UpdatePixelTrait) == 0))
1503 continue;
1504 x_pixel=QuantumScale*reference[i];
1505 x_pixel_mu[i]+=(*k)*x_pixel;
1506 x_pixel_sigma_squared[i]+=(*k)*x_pixel*x_pixel;
1507 y_pixel=QuantumScale*
1508 GetPixelChannel(reconstruct_image,channel,target);
1509 y_pixel_mu[i]+=(*k)*y_pixel;
1510 y_pixel_sigma_squared[i]+=(*k)*y_pixel*y_pixel;
1511 xy_sigma[i]+=(*k)*x_pixel*y_pixel;
1512 }
1513 k++;
1514 reference+=GetPixelChannels(image);
1515 target+=GetPixelChannels(reconstruct_image);
1516 }
1517 reference+=GetPixelChannels(image)*columns;
1518 target+=GetPixelChannels(reconstruct_image)*columns;
1519 }
1520 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1521 {
1522 double
1523 ssim,
1524 x_pixel_mu_squared,
1525 x_pixel_sigmas_squared,
1526 xy_mu,
1527 xy_sigmas,
1528 y_pixel_mu_squared,
1529 y_pixel_sigmas_squared;
1530
1531 PixelChannel channel = GetPixelChannelChannel(image,i);
1532 PixelTrait traits = GetPixelChannelTraits(image,channel);
1533 PixelTrait reconstruct_traits = GetPixelChannelTraits(
1534 reconstruct_image,channel);
1535 if ((traits == UndefinedPixelTrait) ||
1536 (reconstruct_traits == UndefinedPixelTrait) ||
1537 ((reconstruct_traits & UpdatePixelTrait) == 0))
1538 continue;
1539 x_pixel_mu_squared=x_pixel_mu[i]*x_pixel_mu[i];
1540 y_pixel_mu_squared=y_pixel_mu[i]*y_pixel_mu[i];
1541 xy_mu=x_pixel_mu[i]*y_pixel_mu[i];
1542 xy_sigmas=xy_sigma[i]-xy_mu;
1543 x_pixel_sigmas_squared=x_pixel_sigma_squared[i]-x_pixel_mu_squared;
1544 y_pixel_sigmas_squared=y_pixel_sigma_squared[i]-y_pixel_mu_squared;
1545 ssim=((2.0*xy_mu+c1)*(2.0*xy_sigmas+c2))/
1546 ((x_pixel_mu_squared+y_pixel_mu_squared+c1)*
1547 (x_pixel_sigmas_squared+y_pixel_sigmas_squared+c2));
1548 channel_distortion[i]+=ssim;
1549 channel_distortion[CompositePixelChannel]+=ssim;
1550 }
1551 p+=GetPixelChannels(image);
1552 q+=GetPixelChannels(reconstruct_image);
1553 }
1554 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1555 #pragma omp critical (MagickCore_GetStructuralSimilarityDistortion)
1556 #endif
1557 for (i=0; i <= MaxPixelChannels; i++)
1558 distortion[i]+=channel_distortion[i];
1559 }
1560 image_view=DestroyCacheView(image_view);
1561 reconstruct_view=DestroyCacheView(reconstruct_view);
1562 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1563 {
1564 PixelChannel channel = GetPixelChannelChannel(image,i);
1565 PixelTrait traits = GetPixelChannelTraits(image,channel);
1566 if ((traits == UndefinedPixelTrait) || ((traits & UpdatePixelTrait) == 0))
1567 continue;
1568 distortion[i]/=((double) columns*rows);
1569 }
1570 distortion[CompositePixelChannel]/=((double) columns*rows);
1571 distortion[CompositePixelChannel]/=(double) GetImageChannels(image);
1572 kernel_info=DestroyKernelInfo(kernel_info);
1573 return(status);
1574 }
1575
GetStructuralDisimilarityDistortion(const Image * image,const Image * reconstruct_image,double * distortion,ExceptionInfo * exception)1576 static MagickBooleanType GetStructuralDisimilarityDistortion(const Image *image,
1577 const Image *reconstruct_image,double *distortion,ExceptionInfo *exception)
1578 {
1579 MagickBooleanType
1580 status;
1581
1582 register ssize_t
1583 i;
1584
1585 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1586 distortion,exception);
1587 for (i=0; i <= MaxPixelChannels; i++)
1588 distortion[i]=(1.0-(distortion[i]))/2.0;
1589 return(status);
1590 }
1591
GetImageDistortion(Image * image,const Image * reconstruct_image,const MetricType metric,double * distortion,ExceptionInfo * exception)1592 MagickExport MagickBooleanType GetImageDistortion(Image *image,
1593 const Image *reconstruct_image,const MetricType metric,double *distortion,
1594 ExceptionInfo *exception)
1595 {
1596 double
1597 *channel_distortion;
1598
1599 MagickBooleanType
1600 status;
1601
1602 size_t
1603 length;
1604
1605 assert(image != (Image *) NULL);
1606 assert(image->signature == MagickCoreSignature);
1607 if (image->debug != MagickFalse)
1608 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1609 assert(reconstruct_image != (const Image *) NULL);
1610 assert(reconstruct_image->signature == MagickCoreSignature);
1611 assert(distortion != (double *) NULL);
1612 *distortion=0.0;
1613 if (image->debug != MagickFalse)
1614 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1615 /*
1616 Get image distortion.
1617 */
1618 length=MaxPixelChannels+1;
1619 channel_distortion=(double *) AcquireQuantumMemory(length,
1620 sizeof(*channel_distortion));
1621 if (channel_distortion == (double *) NULL)
1622 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1623 (void) memset(channel_distortion,0,length*
1624 sizeof(*channel_distortion));
1625 switch (metric)
1626 {
1627 case AbsoluteErrorMetric:
1628 {
1629 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1630 exception);
1631 break;
1632 }
1633 case FuzzErrorMetric:
1634 {
1635 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1636 exception);
1637 break;
1638 }
1639 case MeanAbsoluteErrorMetric:
1640 {
1641 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1642 channel_distortion,exception);
1643 break;
1644 }
1645 case MeanErrorPerPixelErrorMetric:
1646 {
1647 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1648 exception);
1649 break;
1650 }
1651 case MeanSquaredErrorMetric:
1652 {
1653 status=GetMeanSquaredDistortion(image,reconstruct_image,
1654 channel_distortion,exception);
1655 break;
1656 }
1657 case NormalizedCrossCorrelationErrorMetric:
1658 default:
1659 {
1660 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1661 channel_distortion,exception);
1662 break;
1663 }
1664 case PeakAbsoluteErrorMetric:
1665 {
1666 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1667 channel_distortion,exception);
1668 break;
1669 }
1670 case PeakSignalToNoiseRatioErrorMetric:
1671 {
1672 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1673 channel_distortion,exception);
1674 break;
1675 }
1676 case PerceptualHashErrorMetric:
1677 {
1678 status=GetPerceptualHashDistortion(image,reconstruct_image,
1679 channel_distortion,exception);
1680 break;
1681 }
1682 case RootMeanSquaredErrorMetric:
1683 {
1684 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1685 channel_distortion,exception);
1686 break;
1687 }
1688 case StructuralSimilarityErrorMetric:
1689 {
1690 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1691 channel_distortion,exception);
1692 break;
1693 }
1694 case StructuralDissimilarityErrorMetric:
1695 {
1696 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1697 channel_distortion,exception);
1698 break;
1699 }
1700 }
1701 *distortion=channel_distortion[CompositePixelChannel];
1702 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1703 (void) FormatImageProperty(image,"distortion","%.*g",GetMagickPrecision(),
1704 *distortion);
1705 return(status);
1706 }
1707
1708 /*
1709 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1710 % %
1711 % %
1712 % %
1713 % G e t I m a g e D i s t o r t i o n s %
1714 % %
1715 % %
1716 % %
1717 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1718 %
1719 % GetImageDistortions() compares the pixel channels of an image to a
1720 % reconstructed image and returns the specified distortion metric for each
1721 % channel.
1722 %
1723 % The format of the GetImageDistortions method is:
1724 %
1725 % double *GetImageDistortions(const Image *image,
1726 % const Image *reconstruct_image,const MetricType metric,
1727 % ExceptionInfo *exception)
1728 %
1729 % A description of each parameter follows:
1730 %
1731 % o image: the image.
1732 %
1733 % o reconstruct_image: the reconstruct image.
1734 %
1735 % o metric: the metric.
1736 %
1737 % o exception: return any errors or warnings in this structure.
1738 %
1739 */
GetImageDistortions(Image * image,const Image * reconstruct_image,const MetricType metric,ExceptionInfo * exception)1740 MagickExport double *GetImageDistortions(Image *image,
1741 const Image *reconstruct_image,const MetricType metric,
1742 ExceptionInfo *exception)
1743 {
1744 double
1745 *channel_distortion;
1746
1747 MagickBooleanType
1748 status;
1749
1750 size_t
1751 length;
1752
1753 assert(image != (Image *) NULL);
1754 assert(image->signature == MagickCoreSignature);
1755 if (image->debug != MagickFalse)
1756 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1757 assert(reconstruct_image != (const Image *) NULL);
1758 assert(reconstruct_image->signature == MagickCoreSignature);
1759 if (image->debug != MagickFalse)
1760 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1761 /*
1762 Get image distortion.
1763 */
1764 length=MaxPixelChannels+1UL;
1765 channel_distortion=(double *) AcquireQuantumMemory(length,
1766 sizeof(*channel_distortion));
1767 if (channel_distortion == (double *) NULL)
1768 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1769 (void) memset(channel_distortion,0,length*
1770 sizeof(*channel_distortion));
1771 status=MagickTrue;
1772 switch (metric)
1773 {
1774 case AbsoluteErrorMetric:
1775 {
1776 status=GetAbsoluteDistortion(image,reconstruct_image,channel_distortion,
1777 exception);
1778 break;
1779 }
1780 case FuzzErrorMetric:
1781 {
1782 status=GetFuzzDistortion(image,reconstruct_image,channel_distortion,
1783 exception);
1784 break;
1785 }
1786 case MeanAbsoluteErrorMetric:
1787 {
1788 status=GetMeanAbsoluteDistortion(image,reconstruct_image,
1789 channel_distortion,exception);
1790 break;
1791 }
1792 case MeanErrorPerPixelErrorMetric:
1793 {
1794 status=GetMeanErrorPerPixel(image,reconstruct_image,channel_distortion,
1795 exception);
1796 break;
1797 }
1798 case MeanSquaredErrorMetric:
1799 {
1800 status=GetMeanSquaredDistortion(image,reconstruct_image,
1801 channel_distortion,exception);
1802 break;
1803 }
1804 case NormalizedCrossCorrelationErrorMetric:
1805 default:
1806 {
1807 status=GetNormalizedCrossCorrelationDistortion(image,reconstruct_image,
1808 channel_distortion,exception);
1809 break;
1810 }
1811 case PeakAbsoluteErrorMetric:
1812 {
1813 status=GetPeakAbsoluteDistortion(image,reconstruct_image,
1814 channel_distortion,exception);
1815 break;
1816 }
1817 case PeakSignalToNoiseRatioErrorMetric:
1818 {
1819 status=GetPeakSignalToNoiseRatio(image,reconstruct_image,
1820 channel_distortion,exception);
1821 break;
1822 }
1823 case PerceptualHashErrorMetric:
1824 {
1825 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1826 channel_distortion,exception);
1827 break;
1828 }
1829 case RootMeanSquaredErrorMetric:
1830 {
1831 status=GetRootMeanSquaredDistortion(image,reconstruct_image,
1832 channel_distortion,exception);
1833 break;
1834 }
1835 case StructuralSimilarityErrorMetric:
1836 {
1837 status=GetStructuralSimilarityDistortion(image,reconstruct_image,
1838 channel_distortion,exception);
1839 break;
1840 }
1841 case StructuralDissimilarityErrorMetric:
1842 {
1843 status=GetStructuralDisimilarityDistortion(image,reconstruct_image,
1844 channel_distortion,exception);
1845 break;
1846 }
1847 }
1848 if (status == MagickFalse)
1849 {
1850 channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1851 return((double *) NULL);
1852 }
1853 return(channel_distortion);
1854 }
1855
1856 /*
1857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1858 % %
1859 % %
1860 % %
1861 % I s I m a g e s E q u a l %
1862 % %
1863 % %
1864 % %
1865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1866 %
1867 % IsImagesEqual() compare the pixels of two images and returns immediately
1868 % if any pixel is not identical.
1869 %
1870 % The format of the IsImagesEqual method is:
1871 %
1872 % MagickBooleanType IsImagesEqual(const Image *image,
1873 % const Image *reconstruct_image,ExceptionInfo *exception)
1874 %
1875 % A description of each parameter follows.
1876 %
1877 % o image: the image.
1878 %
1879 % o reconstruct_image: the reconstruct image.
1880 %
1881 % o exception: return any errors or warnings in this structure.
1882 %
1883 */
IsImagesEqual(const Image * image,const Image * reconstruct_image,ExceptionInfo * exception)1884 MagickExport MagickBooleanType IsImagesEqual(const Image *image,
1885 const Image *reconstruct_image,ExceptionInfo *exception)
1886 {
1887 CacheView
1888 *image_view,
1889 *reconstruct_view;
1890
1891 size_t
1892 columns,
1893 rows;
1894
1895 ssize_t
1896 y;
1897
1898 assert(image != (Image *) NULL);
1899 assert(image->signature == MagickCoreSignature);
1900 assert(reconstruct_image != (const Image *) NULL);
1901 assert(reconstruct_image->signature == MagickCoreSignature);
1902 rows=MagickMax(image->rows,reconstruct_image->rows);
1903 columns=MagickMax(image->columns,reconstruct_image->columns);
1904 image_view=AcquireVirtualCacheView(image,exception);
1905 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
1906 for (y=0; y < (ssize_t) rows; y++)
1907 {
1908 register const Quantum
1909 *magick_restrict p,
1910 *magick_restrict q;
1911
1912 register ssize_t
1913 x;
1914
1915 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
1916 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
1917 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1918 break;
1919 for (x=0; x < (ssize_t) columns; x++)
1920 {
1921 register ssize_t
1922 i;
1923
1924 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1925 {
1926 double
1927 distance;
1928
1929 PixelChannel channel = GetPixelChannelChannel(image,i);
1930 PixelTrait traits = GetPixelChannelTraits(image,channel);
1931 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
1932 channel);
1933 if ((traits == UndefinedPixelTrait) ||
1934 (reconstruct_traits == UndefinedPixelTrait) ||
1935 ((reconstruct_traits & UpdatePixelTrait) == 0))
1936 continue;
1937 distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
1938 channel,q));
1939 if (distance >= MagickEpsilon)
1940 break;
1941 }
1942 if (i < (ssize_t) GetPixelChannels(image))
1943 break;
1944 p+=GetPixelChannels(image);
1945 q+=GetPixelChannels(reconstruct_image);
1946 }
1947 if (x < (ssize_t) columns)
1948 break;
1949 }
1950 reconstruct_view=DestroyCacheView(reconstruct_view);
1951 image_view=DestroyCacheView(image_view);
1952 return(y < (ssize_t) rows ? MagickFalse : MagickTrue);
1953 }
1954
1955 /*
1956 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1957 % %
1958 % %
1959 % %
1960 % S e t I m a g e C o l o r M e t r i c %
1961 % %
1962 % %
1963 % %
1964 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1965 %
1966 % SetImageColorMetric() measures the difference between colors at each pixel
1967 % location of two images. A value other than 0 means the colors match
1968 % exactly. Otherwise an error measure is computed by summing over all
1969 % pixels in an image the distance squared in RGB space between each image
1970 % pixel and its corresponding pixel in the reconstruct image. The error
1971 % measure is assigned to these image members:
1972 %
1973 % o mean_error_per_pixel: The mean error for any single pixel in
1974 % the image.
1975 %
1976 % o normalized_mean_error: The normalized mean quantization error for
1977 % any single pixel in the image. This distance measure is normalized to
1978 % a range between 0 and 1. It is independent of the range of red, green,
1979 % and blue values in the image.
1980 %
1981 % o normalized_maximum_error: The normalized maximum quantization
1982 % error for any single pixel in the image. This distance measure is
1983 % normalized to a range between 0 and 1. It is independent of the range
1984 % of red, green, and blue values in your image.
1985 %
1986 % A small normalized mean square error, accessed as
1987 % image->normalized_mean_error, suggests the images are very similar in
1988 % spatial layout and color.
1989 %
1990 % The format of the SetImageColorMetric method is:
1991 %
1992 % MagickBooleanType SetImageColorMetric(Image *image,
1993 % const Image *reconstruct_image,ExceptionInfo *exception)
1994 %
1995 % A description of each parameter follows.
1996 %
1997 % o image: the image.
1998 %
1999 % o reconstruct_image: the reconstruct image.
2000 %
2001 % o exception: return any errors or warnings in this structure.
2002 %
2003 */
SetImageColorMetric(Image * image,const Image * reconstruct_image,ExceptionInfo * exception)2004 MagickExport MagickBooleanType SetImageColorMetric(Image *image,
2005 const Image *reconstruct_image,ExceptionInfo *exception)
2006 {
2007 CacheView
2008 *image_view,
2009 *reconstruct_view;
2010
2011 double
2012 area,
2013 maximum_error,
2014 mean_error,
2015 mean_error_per_pixel;
2016
2017 MagickBooleanType
2018 status;
2019
2020 size_t
2021 columns,
2022 rows;
2023
2024 ssize_t
2025 y;
2026
2027 assert(image != (Image *) NULL);
2028 assert(image->signature == MagickCoreSignature);
2029 assert(reconstruct_image != (const Image *) NULL);
2030 assert(reconstruct_image->signature == MagickCoreSignature);
2031 area=0.0;
2032 maximum_error=0.0;
2033 mean_error_per_pixel=0.0;
2034 mean_error=0.0;
2035 rows=MagickMax(image->rows,reconstruct_image->rows);
2036 columns=MagickMax(image->columns,reconstruct_image->columns);
2037 image_view=AcquireVirtualCacheView(image,exception);
2038 reconstruct_view=AcquireVirtualCacheView(reconstruct_image,exception);
2039 for (y=0; y < (ssize_t) rows; y++)
2040 {
2041 register const Quantum
2042 *magick_restrict p,
2043 *magick_restrict q;
2044
2045 register ssize_t
2046 x;
2047
2048 p=GetCacheViewVirtualPixels(image_view,0,y,columns,1,exception);
2049 q=GetCacheViewVirtualPixels(reconstruct_view,0,y,columns,1,exception);
2050 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2051 break;
2052 for (x=0; x < (ssize_t) columns; x++)
2053 {
2054 register ssize_t
2055 i;
2056
2057 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2058 {
2059 double
2060 distance;
2061
2062 PixelChannel channel = GetPixelChannelChannel(image,i);
2063 PixelTrait traits = GetPixelChannelTraits(image,channel);
2064 PixelTrait reconstruct_traits = GetPixelChannelTraits(reconstruct_image,
2065 channel);
2066 if ((traits == UndefinedPixelTrait) ||
2067 (reconstruct_traits == UndefinedPixelTrait) ||
2068 ((reconstruct_traits & UpdatePixelTrait) == 0))
2069 continue;
2070 distance=fabs(p[i]-(double) GetPixelChannel(reconstruct_image,
2071 channel,q));
2072 if (distance >= MagickEpsilon)
2073 {
2074 mean_error_per_pixel+=distance;
2075 mean_error+=distance*distance;
2076 if (distance > maximum_error)
2077 maximum_error=distance;
2078 }
2079 area++;
2080 }
2081 p+=GetPixelChannels(image);
2082 q+=GetPixelChannels(reconstruct_image);
2083 }
2084 }
2085 reconstruct_view=DestroyCacheView(reconstruct_view);
2086 image_view=DestroyCacheView(image_view);
2087 image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
2088 image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
2089 mean_error/area);
2090 image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
2091 status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
2092 return(status);
2093 }
2094
2095 /*
2096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2097 % %
2098 % %
2099 % %
2100 % S i m i l a r i t y I m a g e %
2101 % %
2102 % %
2103 % %
2104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2105 %
2106 % SimilarityImage() compares the reference image of the image and returns the
2107 % best match offset. In addition, it returns a similarity image such that an
2108 % exact match location is completely white and if none of the pixels match,
2109 % black, otherwise some gray level in-between.
2110 %
2111 % The format of the SimilarityImageImage method is:
2112 %
2113 % Image *SimilarityImage(const Image *image,const Image *reference,
2114 % const MetricType metric,const double similarity_threshold,
2115 % RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
2116 %
2117 % A description of each parameter follows:
2118 %
2119 % o image: the image.
2120 %
2121 % o reference: find an area of the image that closely resembles this image.
2122 %
2123 % o metric: the metric.
2124 %
2125 % o similarity_threshold: minimum distortion for (sub)image match.
2126 %
2127 % o offset: the best match offset of the reference image within the image.
2128 %
2129 % o similarity: the computed similarity between the images.
2130 %
2131 % o exception: return any errors or warnings in this structure.
2132 %
2133 */
2134
GetSimilarityMetric(const Image * image,const Image * reference,const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo * exception)2135 static double GetSimilarityMetric(const Image *image,const Image *reference,
2136 const MetricType metric,const ssize_t x_offset,const ssize_t y_offset,
2137 ExceptionInfo *exception)
2138 {
2139 double
2140 distortion;
2141
2142 Image
2143 *similarity_image;
2144
2145 MagickBooleanType
2146 status;
2147
2148 RectangleInfo
2149 geometry;
2150
2151 SetGeometry(reference,&geometry);
2152 geometry.x=x_offset;
2153 geometry.y=y_offset;
2154 similarity_image=CropImage(image,&geometry,exception);
2155 if (similarity_image == (Image *) NULL)
2156 return(0.0);
2157 distortion=0.0;
2158 status=GetImageDistortion(similarity_image,reference,metric,&distortion,
2159 exception);
2160 similarity_image=DestroyImage(similarity_image);
2161 if (status == MagickFalse)
2162 return(0.0);
2163 return(distortion);
2164 }
2165
SimilarityImage(const Image * image,const Image * reference,const MetricType metric,const double similarity_threshold,RectangleInfo * offset,double * similarity_metric,ExceptionInfo * exception)2166 MagickExport Image *SimilarityImage(const Image *image,const Image *reference,
2167 const MetricType metric,const double similarity_threshold,
2168 RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
2169 {
2170 #define SimilarityImageTag "Similarity/Image"
2171
2172 CacheView
2173 *similarity_view;
2174
2175 Image
2176 *similarity_image;
2177
2178 MagickBooleanType
2179 status;
2180
2181 MagickOffsetType
2182 progress;
2183
2184 ssize_t
2185 y;
2186
2187 assert(image != (const Image *) NULL);
2188 assert(image->signature == MagickCoreSignature);
2189 if (image->debug != MagickFalse)
2190 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2191 assert(exception != (ExceptionInfo *) NULL);
2192 assert(exception->signature == MagickCoreSignature);
2193 assert(offset != (RectangleInfo *) NULL);
2194 SetGeometry(reference,offset);
2195 *similarity_metric=MagickMaximumValue;
2196 similarity_image=CloneImage(image,image->columns-reference->columns+1,
2197 image->rows-reference->rows+1,MagickTrue,exception);
2198 if (similarity_image == (Image *) NULL)
2199 return((Image *) NULL);
2200 status=SetImageStorageClass(similarity_image,DirectClass,exception);
2201 if (status == MagickFalse)
2202 {
2203 similarity_image=DestroyImage(similarity_image);
2204 return((Image *) NULL);
2205 }
2206 (void) SetImageAlphaChannel(similarity_image,DeactivateAlphaChannel,
2207 exception);
2208 /*
2209 Measure similarity of reference image against image.
2210 */
2211 status=MagickTrue;
2212 progress=0;
2213 similarity_view=AcquireAuthenticCacheView(similarity_image,exception);
2214 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2215 #pragma omp parallel for schedule(static) \
2216 shared(progress,status,similarity_metric) \
2217 magick_number_threads(image,image,image->rows-reference->rows+1,1)
2218 #endif
2219 for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
2220 {
2221 double
2222 similarity;
2223
2224 register Quantum
2225 *magick_restrict q;
2226
2227 register ssize_t
2228 x;
2229
2230 if (status == MagickFalse)
2231 continue;
2232 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2233 #pragma omp flush(similarity_metric)
2234 #endif
2235 if (*similarity_metric <= similarity_threshold)
2236 continue;
2237 q=GetCacheViewAuthenticPixels(similarity_view,0,y,similarity_image->columns,
2238 1,exception);
2239 if (q == (Quantum *) NULL)
2240 {
2241 status=MagickFalse;
2242 continue;
2243 }
2244 for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
2245 {
2246 register ssize_t
2247 i;
2248
2249 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2250 #pragma omp flush(similarity_metric)
2251 #endif
2252 if (*similarity_metric <= similarity_threshold)
2253 break;
2254 similarity=GetSimilarityMetric(image,reference,metric,x,y,exception);
2255 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2256 #pragma omp critical (MagickCore_SimilarityImage)
2257 #endif
2258 if ((metric == NormalizedCrossCorrelationErrorMetric) ||
2259 (metric == UndefinedErrorMetric))
2260 similarity=1.0-similarity;
2261 if (similarity < *similarity_metric)
2262 {
2263 offset->x=x;
2264 offset->y=y;
2265 *similarity_metric=similarity;
2266 }
2267 if (metric == PerceptualHashErrorMetric)
2268 similarity=MagickMin(0.01*similarity,1.0);
2269 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2270 {
2271 PixelChannel channel = GetPixelChannelChannel(image,i);
2272 PixelTrait traits = GetPixelChannelTraits(image,channel);
2273 PixelTrait similarity_traits=GetPixelChannelTraits(similarity_image,
2274 channel);
2275 if ((traits == UndefinedPixelTrait) ||
2276 (similarity_traits == UndefinedPixelTrait) ||
2277 ((similarity_traits & UpdatePixelTrait) == 0))
2278 continue;
2279 SetPixelChannel(similarity_image,channel,ClampToQuantum(QuantumRange-
2280 QuantumRange*similarity),q);
2281 }
2282 q+=GetPixelChannels(similarity_image);
2283 }
2284 if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
2285 status=MagickFalse;
2286 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2287 {
2288 MagickBooleanType
2289 proceed;
2290
2291 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2292 #pragma omp atomic
2293 #endif
2294 progress++;
2295 proceed=SetImageProgress(image,SimilarityImageTag,progress,image->rows);
2296 if (proceed == MagickFalse)
2297 status=MagickFalse;
2298 }
2299 }
2300 similarity_view=DestroyCacheView(similarity_view);
2301 if (status == MagickFalse)
2302 similarity_image=DestroyImage(similarity_image);
2303 return(similarity_image);
2304 }
2305