]> granicus.if.org Git - imagemagick/blob - magick/compare.c
(no commit message)
[imagemagick] / magick / compare.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %                                                                             %
4 %                                                                             %
5 %                                                                             %
6 %               CCCC   OOO   M   M  PPPP    AAA   RRRR    EEEEE               %
7 %              C      O   O  MM MM  P   P  A   A  R   R   E                   %
8 %              C      O   O  M M M  PPPP   AAAAA  RRRR    EEE                 %
9 %              C      O   O  M   M  P      A   A  R R     E                   %
10 %               CCCC   OOO   M   M  P      A   A  R  R    EEEEE               %
11 %                                                                             %
12 %                                                                             %
13 %                      MagickCore Image Comparison Methods                    %
14 %                                                                             %
15 %                              Software Design                                %
16 %                                John Cristy                                  %
17 %                               December 2003                                 %
18 %                                                                             %
19 %                                                                             %
20 %  Copyright 1999-2010 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 %    http://www.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 \f
40 /*
41   Include declarations.
42 */
43 #include "magick/studio.h"
44 #include "magick/artifact.h"
45 #include "magick/cache-view.h"
46 #include "magick/client.h"
47 #include "magick/color.h"
48 #include "magick/color-private.h"
49 #include "magick/colorspace.h"
50 #include "magick/colorspace-private.h"
51 #include "magick/compare.h"
52 #include "magick/composite-private.h"
53 #include "magick/constitute.h"
54 #include "magick/exception-private.h"
55 #include "magick/geometry.h"
56 #include "magick/image-private.h"
57 #include "magick/list.h"
58 #include "magick/log.h"
59 #include "magick/memory_.h"
60 #include "magick/monitor.h"
61 #include "magick/monitor-private.h"
62 #include "magick/option.h"
63 #include "magick/pixel-private.h"
64 #include "magick/resource_.h"
65 #include "magick/string_.h"
66 #include "magick/statistic.h"
67 #include "magick/utility.h"
68 #include "magick/version.h"
69 \f
70 /*
71 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72 %                                                                             %
73 %                                                                             %
74 %                                                                             %
75 %   C o m p a r e I m a g e C h a n n e l s                                   %
76 %                                                                             %
77 %                                                                             %
78 %                                                                             %
79 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 %
81 %  CompareImageChannels() compares one or more image channels of an image
82 %  to a reconstructed image and returns the difference image.
83 %
84 %  The format of the CompareImageChannels method is:
85 %
86 %      Image *CompareImageChannels(const Image *image,
87 %        const Image *reconstruct_image,const ChannelType channel,
88 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
89 %
90 %  A description of each parameter follows:
91 %
92 %    o image: the image.
93 %
94 %    o reconstruct_image: the reconstruct image.
95 %
96 %    o channel: the channel.
97 %
98 %    o metric: the metric.
99 %
100 %    o distortion: the computed distortion between the images.
101 %
102 %    o exception: return any errors or warnings in this structure.
103 %
104 */
105
106 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
107   const MetricType metric,double *distortion,ExceptionInfo *exception)
108 {
109   Image
110     *highlight_image;
111
112   highlight_image=CompareImageChannels(image,reconstruct_image,AllChannels,
113     metric,distortion,exception);
114   return(highlight_image);
115 }
116
117 MagickExport Image *CompareImageChannels(Image *image,
118   const Image *reconstruct_image,const ChannelType channel,
119   const MetricType metric,double *distortion,ExceptionInfo *exception)
120 {
121   CacheView
122     *highlight_view,
123     *image_view,
124     *reconstruct_view;
125
126   const char
127     *artifact;
128
129   Image
130     *difference_image,
131     *highlight_image;
132
133   ssize_t
134     y;
135
136   MagickBooleanType
137     status;
138
139   MagickPixelPacket
140     highlight,
141     lowlight,
142     zero;
143
144   assert(image != (Image *) NULL);
145   assert(image->signature == MagickSignature);
146   if (image->debug != MagickFalse)
147     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
148   assert(reconstruct_image != (const Image *) NULL);
149   assert(reconstruct_image->signature == MagickSignature);
150   assert(distortion != (double *) NULL);
151   *distortion=0.0;
152   if (image->debug != MagickFalse)
153     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
154   if ((reconstruct_image->columns != image->columns) ||
155       (reconstruct_image->rows != image->rows))
156     ThrowImageException(ImageError,"ImageSizeDiffers");
157   status=GetImageChannelDistortion(image,reconstruct_image,channel,metric,
158     distortion,exception);
159   if (status == MagickFalse)
160     return((Image *) NULL);
161   difference_image=CloneImage(image,0,0,MagickTrue,exception);
162   if (difference_image == (Image *) NULL)
163     return((Image *) NULL);
164   (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel);
165   highlight_image=CloneImage(image,image->columns,image->rows,MagickTrue,
166     exception);
167   if (highlight_image == (Image *) NULL)
168     {
169       difference_image=DestroyImage(difference_image);
170       return((Image *) NULL);
171     }
172   if (SetImageStorageClass(highlight_image,DirectClass) == MagickFalse)
173     {
174       InheritException(exception,&highlight_image->exception);
175       difference_image=DestroyImage(difference_image);
176       highlight_image=DestroyImage(highlight_image);
177       return((Image *) NULL);
178     }
179   (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel);
180   (void) QueryMagickColor("#f1001ecc",&highlight,exception);
181   artifact=GetImageArtifact(image,"highlight-color");
182   if (artifact != (const char *) NULL)
183     (void) QueryMagickColor(artifact,&highlight,exception);
184   (void) QueryMagickColor("#ffffffcc",&lowlight,exception);
185   artifact=GetImageArtifact(image,"lowlight-color");
186   if (artifact != (const char *) NULL)
187     (void) QueryMagickColor(artifact,&lowlight,exception);
188   if (highlight_image->colorspace == CMYKColorspace)
189     {
190       ConvertRGBToCMYK(&highlight);
191       ConvertRGBToCMYK(&lowlight);
192     }
193   /*
194     Generate difference image.
195   */
196   status=MagickTrue;
197   GetMagickPixelPacket(image,&zero);
198   image_view=AcquireCacheView(image);
199   reconstruct_view=AcquireCacheView(reconstruct_image);
200   highlight_view=AcquireCacheView(highlight_image);
201 #if defined(MAGICKCORE_OPENMP_SUPPORT)
202   #pragma omp parallel for schedule(dynamic,4) shared(status)
203 #endif
204   for (y=0; y < (ssize_t) image->rows; y++)
205   {
206     MagickBooleanType
207       sync;
208
209     MagickPixelPacket
210       pixel,
211       reconstruct_pixel;
212
213     register const IndexPacket
214       *restrict indexes,
215       *restrict reconstruct_indexes;
216
217     register const PixelPacket
218       *restrict p,
219       *restrict q;
220
221     register IndexPacket
222       *restrict highlight_indexes;
223
224     register ssize_t
225       x;
226
227     register PixelPacket
228       *restrict r;
229
230     if (status == MagickFalse)
231       continue;
232     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
233     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
234       1,exception);
235     r=QueueCacheViewAuthenticPixels(highlight_view,0,y,highlight_image->columns,
236       1,exception);
237     if ((p == (const PixelPacket *) NULL) ||
238         (q == (const PixelPacket *) NULL) || (r == (PixelPacket *) NULL))
239       {
240         status=MagickFalse;
241         continue;
242       }
243     indexes=GetCacheViewVirtualIndexQueue(image_view);
244     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
245     highlight_indexes=GetCacheViewAuthenticIndexQueue(highlight_view);
246     pixel=zero;
247     reconstruct_pixel=zero;
248     for (x=0; x < (ssize_t) image->columns; x++)
249     {
250       MagickStatusType
251         difference;
252
253       SetMagickPixelPacket(image,p,indexes+x,&pixel);
254       SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
255         &reconstruct_pixel);
256       difference=MagickFalse;
257       if (channel == AllChannels)
258         {
259           if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
260             difference=MagickTrue;
261         }
262       else
263         {
264           if (((channel & RedChannel) != 0) && (p->red != q->red))
265             difference=MagickTrue;
266           if (((channel & GreenChannel) != 0) && (p->green != q->green))
267             difference=MagickTrue;
268           if (((channel & BlueChannel) != 0) && (p->blue != q->blue))
269             difference=MagickTrue;
270           if (((channel & OpacityChannel) != 0) &&
271               (image->matte != MagickFalse) && (p->opacity != q->opacity))
272             difference=MagickTrue;
273           if ((((channel & IndexChannel) != 0) &&
274                (image->colorspace == CMYKColorspace) &&
275                (reconstruct_image->colorspace == CMYKColorspace)) &&
276               (indexes[x] != reconstruct_indexes[x]))
277             difference=MagickTrue;
278         }
279       if (difference != MagickFalse)
280         SetPixelPacket(highlight_image,&highlight,r,highlight_indexes+x);
281       else
282         SetPixelPacket(highlight_image,&lowlight,r,highlight_indexes+x);
283       p++;
284       q++;
285       r++;
286     }
287     sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
288     if (sync == MagickFalse)
289       status=MagickFalse;
290   }
291   highlight_view=DestroyCacheView(highlight_view);
292   reconstruct_view=DestroyCacheView(reconstruct_view);
293   image_view=DestroyCacheView(image_view);
294   (void) CompositeImage(difference_image,image->compose,highlight_image,0,0);
295   highlight_image=DestroyImage(highlight_image);
296   if (status == MagickFalse)
297     difference_image=DestroyImage(difference_image);
298   return(difference_image);
299 }
300 \f
301 /*
302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
303 %                                                                             %
304 %                                                                             %
305 %                                                                             %
306 %   G e t I m a g e C h a n n e l D i s t o r t i o n                         %
307 %                                                                             %
308 %                                                                             %
309 %                                                                             %
310 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
311 %
312 %  GetImageChannelDistortion() compares one or more image channels of an image
313 %  to a reconstructed image and returns the specified distortion metric.
314 %
315 %  The format of the CompareImageChannels method is:
316 %
317 %      MagickBooleanType GetImageChannelDistortion(const Image *image,
318 %        const Image *reconstruct_image,const ChannelType channel,
319 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
320 %
321 %  A description of each parameter follows:
322 %
323 %    o image: the image.
324 %
325 %    o reconstruct_image: the reconstruct image.
326 %
327 %    o channel: the channel.
328 %
329 %    o metric: the metric.
330 %
331 %    o distortion: the computed distortion between the images.
332 %
333 %    o exception: return any errors or warnings in this structure.
334 %
335 */
336
337 MagickExport MagickBooleanType GetImageDistortion(Image *image,
338   const Image *reconstruct_image,const MetricType metric,double *distortion,
339   ExceptionInfo *exception)
340 {
341   MagickBooleanType
342     status;
343
344   status=GetImageChannelDistortion(image,reconstruct_image,AllChannels,
345     metric,distortion,exception);
346   return(status);
347 }
348
349 static MagickBooleanType GetAbsoluteError(const Image *image,
350   const Image *reconstruct_image,const ChannelType channel,double *distortion,
351   ExceptionInfo *exception)
352 {
353   CacheView
354     *image_view,
355     *reconstruct_view;
356
357   ssize_t
358     y;
359
360   MagickBooleanType
361     status;
362
363   MagickPixelPacket
364     zero;
365
366   /*
367     Compute the absolute difference in pixels between two images.
368   */
369   status=MagickTrue;
370   GetMagickPixelPacket(image,&zero);
371   image_view=AcquireCacheView(image);
372   reconstruct_view=AcquireCacheView(reconstruct_image);
373 #if defined(MAGICKCORE_OPENMP_SUPPORT)
374   #pragma omp parallel for schedule(dynamic,4) shared(status)
375 #endif
376   for (y=0; y < (ssize_t) image->rows; y++)
377   {
378     double
379       channel_distortion[AllChannels+1];
380
381     MagickPixelPacket
382       pixel,
383       reconstruct_pixel;
384
385     register const IndexPacket
386       *restrict indexes,
387       *restrict reconstruct_indexes;
388
389     register const PixelPacket
390       *restrict p,
391       *restrict q;
392
393     register ssize_t
394       i,
395       x;
396
397     if (status == MagickFalse)
398       continue;
399     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
400     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
401       1,exception);
402     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
403       {
404         status=MagickFalse;
405         continue;
406       }
407     indexes=GetCacheViewVirtualIndexQueue(image_view);
408     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
409     pixel=zero;
410     reconstruct_pixel=pixel;
411     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
412     for (x=0; x < (ssize_t) image->columns; x++)
413     {
414       SetMagickPixelPacket(image,p,indexes+x,&pixel);
415       SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
416         &reconstruct_pixel);
417       if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
418         {
419           if ((channel & RedChannel) != 0)
420             channel_distortion[RedChannel]++;
421           if ((channel & GreenChannel) != 0)
422             channel_distortion[GreenChannel]++;
423           if ((channel & BlueChannel) != 0)
424             channel_distortion[BlueChannel]++;
425           if (((channel & OpacityChannel) != 0) &&
426               (image->matte != MagickFalse))
427             channel_distortion[OpacityChannel]++;
428           if (((channel & IndexChannel) != 0) &&
429               (image->colorspace == CMYKColorspace))
430             channel_distortion[BlackChannel]++;
431           channel_distortion[AllChannels]++;
432         }
433       p++;
434       q++;
435     }
436 #if defined(MAGICKCORE_OPENMP_SUPPORT)
437   #pragma omp critical (MagickCore_GetAbsoluteError)
438 #endif
439     for (i=0; i <= (ssize_t) AllChannels; i++)
440       distortion[i]+=channel_distortion[i];
441   }
442   reconstruct_view=DestroyCacheView(reconstruct_view);
443   image_view=DestroyCacheView(image_view);
444   return(status);
445 }
446
447 static size_t GetNumberChannels(const Image *image,
448   const ChannelType channel)
449 {
450   size_t
451     channels;
452
453   channels=0;
454   if ((channel & RedChannel) != 0)
455     channels++;
456   if ((channel & GreenChannel) != 0)
457     channels++;
458   if ((channel & BlueChannel) != 0)
459     channels++;
460   if (((channel & OpacityChannel) != 0) &&
461        (image->matte != MagickFalse))
462     channels++;
463   if (((channel & IndexChannel) != 0) &&
464       (image->colorspace == CMYKColorspace))
465     channels++;
466   return(channels);
467 }
468
469 static MagickBooleanType GetMeanAbsoluteError(const Image *image,
470   const Image *reconstruct_image,const ChannelType channel,
471   double *distortion,ExceptionInfo *exception)
472 {
473   CacheView
474     *image_view,
475     *reconstruct_view;
476
477   ssize_t
478     y;
479
480   MagickBooleanType
481     status;
482
483   register ssize_t
484     i;
485
486   status=MagickTrue;
487   image_view=AcquireCacheView(image);
488   reconstruct_view=AcquireCacheView(reconstruct_image);
489 #if defined(MAGICKCORE_OPENMP_SUPPORT)
490   #pragma omp parallel for schedule(dynamic,4) shared(status)
491 #endif
492   for (y=0; y < (ssize_t) image->rows; y++)
493   {
494     double
495       channel_distortion[AllChannels+1];
496
497     register const IndexPacket
498       *restrict indexes,
499       *restrict reconstruct_indexes;
500
501     register const PixelPacket
502       *restrict p,
503       *restrict q;
504
505     register ssize_t
506       i,
507       x;
508
509     if (status == MagickFalse)
510       continue;
511     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
512     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
513       reconstruct_image->columns,1,exception);
514     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
515       {
516         status=MagickFalse;
517         continue;
518       }
519     indexes=GetCacheViewVirtualIndexQueue(image_view);
520     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
521     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
522     for (x=0; x < (ssize_t) image->columns; x++)
523     {
524       MagickRealType
525         distance;
526
527       if ((channel & RedChannel) != 0)
528         {
529           distance=QuantumScale*fabs(p->red-(double) q->red);
530           channel_distortion[RedChannel]+=distance;
531           channel_distortion[AllChannels]+=distance;
532         }
533       if ((channel & GreenChannel) != 0)
534         {
535           distance=QuantumScale*fabs(p->green-(double) q->green);
536           channel_distortion[GreenChannel]+=distance;
537           channel_distortion[AllChannels]+=distance;
538         }
539       if ((channel & BlueChannel) != 0)
540         {
541           distance=QuantumScale*fabs(p->blue-(double) q->blue);
542           channel_distortion[BlueChannel]+=distance;
543           channel_distortion[AllChannels]+=distance;
544         }
545       if (((channel & OpacityChannel) != 0) &&
546           (image->matte != MagickFalse))
547         {
548           distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
549           channel_distortion[OpacityChannel]+=distance;
550           channel_distortion[AllChannels]+=distance;
551         }
552       if (((channel & IndexChannel) != 0) &&
553           (image->colorspace == CMYKColorspace))
554         {
555           distance=QuantumScale*fabs(indexes[x]-(double)
556             reconstruct_indexes[x]);
557           channel_distortion[BlackChannel]+=distance;
558           channel_distortion[AllChannels]+=distance;
559         }
560       p++;
561       q++;
562     }
563 #if defined(MAGICKCORE_OPENMP_SUPPORT)
564   #pragma omp critical (MagickCore_GetMeanAbsoluteError)
565 #endif
566     for (i=0; i <= (ssize_t) AllChannels; i++)
567       distortion[i]+=channel_distortion[i];
568   }
569   reconstruct_view=DestroyCacheView(reconstruct_view);
570   image_view=DestroyCacheView(image_view);
571   for (i=0; i <= (ssize_t) AllChannels; i++)
572     distortion[i]/=((double) image->columns*image->rows);
573   distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
574   return(status);
575 }
576
577 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
578   const Image *reconstruct_image,const ChannelType channel,double *distortion,
579   ExceptionInfo *exception)
580 {
581   CacheView
582     *image_view,
583     *reconstruct_view;
584
585   ssize_t
586     y;
587
588   MagickBooleanType
589     status;
590
591   MagickRealType
592     alpha,
593     area,
594     beta,
595     maximum_error,
596     mean_error;
597
598   status=MagickTrue;
599   alpha=1.0;
600   beta=1.0;
601   area=0.0;
602   maximum_error=0.0;
603   mean_error=0.0;
604   image_view=AcquireCacheView(image);
605   reconstruct_view=AcquireCacheView(reconstruct_image);
606   for (y=0; y < (ssize_t) image->rows; y++)
607   {
608     register const IndexPacket
609       *restrict indexes,
610       *restrict reconstruct_indexes;
611
612     register const PixelPacket
613       *restrict p,
614       *restrict q;
615
616     register ssize_t
617       x;
618
619     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
620     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
621       1,exception);
622     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
623       {
624         status=MagickFalse;
625         break;
626       }
627     indexes=GetCacheViewVirtualIndexQueue(image_view);
628     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
629     for (x=0; x < (ssize_t) image->columns; x++)
630     {
631       MagickRealType
632         distance;
633
634       if ((channel & OpacityChannel) != 0)
635         {
636           if (image->matte != MagickFalse)
637             alpha=(MagickRealType) (QuantumScale*(GetAlphaPixelComponent(p)));
638           if (reconstruct_image->matte != MagickFalse)
639             beta=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(q));
640         }
641       if ((channel & RedChannel) != 0)
642         {
643           distance=fabs(alpha*p->red-beta*q->red);
644           distortion[RedChannel]+=distance;
645           distortion[AllChannels]+=distance;
646           mean_error+=distance*distance;
647           if (distance > maximum_error)
648             maximum_error=distance;
649           area++;
650         }
651       if ((channel & GreenChannel) != 0)
652         {
653           distance=fabs(alpha*p->green-beta*q->green);
654           distortion[GreenChannel]+=distance;
655           distortion[AllChannels]+=distance;
656           mean_error+=distance*distance;
657           if (distance > maximum_error)
658             maximum_error=distance;
659           area++;
660         }
661       if ((channel & BlueChannel) != 0)
662         {
663           distance=fabs(alpha*p->blue-beta*q->blue);
664           distortion[BlueChannel]+=distance;
665           distortion[AllChannels]+=distance;
666           mean_error+=distance*distance;
667           if (distance > maximum_error)
668             maximum_error=distance;
669           area++;
670         }
671       if (((channel & OpacityChannel) != 0) &&
672           (image->matte != MagickFalse))
673         {
674           distance=fabs((double) p->opacity-q->opacity);
675           distortion[OpacityChannel]+=distance;
676           distortion[AllChannels]+=distance;
677           mean_error+=distance*distance;
678           if (distance > maximum_error)
679             maximum_error=distance;
680           area++;
681         }
682       if (((channel & IndexChannel) != 0) &&
683           (image->colorspace == CMYKColorspace) &&
684           (reconstruct_image->colorspace == CMYKColorspace))
685         {
686           distance=fabs(alpha*indexes[x]-beta*reconstruct_indexes[x]);
687           distortion[BlackChannel]+=distance;
688           distortion[AllChannels]+=distance;
689           mean_error+=distance*distance;
690           if (distance > maximum_error)
691             maximum_error=distance;
692           area++;
693         }
694       p++;
695       q++;
696     }
697   }
698   reconstruct_view=DestroyCacheView(reconstruct_view);
699   image_view=DestroyCacheView(image_view);
700   image->error.mean_error_per_pixel=distortion[AllChannels]/area;
701   image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
702   image->error.normalized_maximum_error=QuantumScale*maximum_error;
703   return(status);
704 }
705
706 static MagickBooleanType GetMeanSquaredError(const Image *image,
707   const Image *reconstruct_image,const ChannelType channel,
708   double *distortion,ExceptionInfo *exception)
709 {
710   CacheView
711     *image_view,
712     *reconstruct_view;
713
714   ssize_t
715     y;
716
717   MagickBooleanType
718     status;
719
720   register ssize_t
721     i;
722
723   status=MagickTrue;
724   image_view=AcquireCacheView(image);
725   reconstruct_view=AcquireCacheView(reconstruct_image);
726 #if defined(MAGICKCORE_OPENMP_SUPPORT)
727   #pragma omp parallel for schedule(dynamic,4) shared(status)
728 #endif
729   for (y=0; y < (ssize_t) image->rows; y++)
730   {
731     double
732       channel_distortion[AllChannels+1];
733
734     register const IndexPacket
735       *restrict indexes,
736       *restrict reconstruct_indexes;
737
738     register const PixelPacket
739       *restrict p,
740       *restrict q;
741
742     register ssize_t
743       i,
744       x;
745
746     if (status == MagickFalse)
747       continue;
748     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
749     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
750       reconstruct_image->columns,1,exception);
751     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
752       {
753         status=MagickFalse;
754         continue;
755       }
756     indexes=GetCacheViewVirtualIndexQueue(image_view);
757     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
758     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
759     for (x=0; x < (ssize_t) image->columns; x++)
760     {
761       MagickRealType
762         distance;
763
764       if ((channel & RedChannel) != 0)
765         {
766           distance=QuantumScale*(p->red-(MagickRealType) q->red);
767           channel_distortion[RedChannel]+=distance*distance;
768           channel_distortion[AllChannels]+=distance*distance;
769         }
770       if ((channel & GreenChannel) != 0)
771         {
772           distance=QuantumScale*(p->green-(MagickRealType) q->green);
773           channel_distortion[GreenChannel]+=distance*distance;
774           channel_distortion[AllChannels]+=distance*distance;
775         }
776       if ((channel & BlueChannel) != 0)
777         {
778           distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
779           channel_distortion[BlueChannel]+=distance*distance;
780           channel_distortion[AllChannels]+=distance*distance;
781         }
782       if (((channel & OpacityChannel) != 0) &&
783           (image->matte != MagickFalse))
784         {
785           distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
786           channel_distortion[OpacityChannel]+=distance*distance;
787           channel_distortion[AllChannels]+=distance*distance;
788         }
789       if (((channel & IndexChannel) != 0) &&
790           (image->colorspace == CMYKColorspace) &&
791           (reconstruct_image->colorspace == CMYKColorspace))
792         {
793           distance=QuantumScale*(indexes[x]-(MagickRealType)
794             reconstruct_indexes[x]);
795           channel_distortion[BlackChannel]+=distance*distance;
796           channel_distortion[AllChannels]+=distance*distance;
797         }
798       p++;
799       q++;
800     }
801 #if defined(MAGICKCORE_OPENMP_SUPPORT)
802   #pragma omp critical (MagickCore_GetMeanSquaredError)
803 #endif
804     for (i=0; i <= (ssize_t) AllChannels; i++)
805       distortion[i]+=channel_distortion[i];
806   }
807   reconstruct_view=DestroyCacheView(reconstruct_view);
808   image_view=DestroyCacheView(image_view);
809   for (i=0; i <= (ssize_t) AllChannels; i++)
810     distortion[i]/=((double) image->columns*image->rows);
811   distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
812   return(status);
813 }
814
815 static MagickBooleanType GetNormalizedCrossCorrelationError(const Image *image,
816   const Image *reconstruct_image,const ChannelType channel,
817   double *distortion,ExceptionInfo *exception)
818 {
819 #define SimilarityImageTag  "Similarity/Image"
820
821   CacheView
822     *correlation_view,
823     *image_view,
824     *reconstruct_view;
825
826   ChannelStatistics
827     *correlation_statistics,
828     *image_statistics,
829     *reconstruct_statistics;
830
831   Image
832     *correlation_image;
833
834   MagickOffsetType
835     progress;
836
837   ssize_t
838     y;
839
840   MagickBooleanType
841     status;
842
843   register ssize_t
844     i;
845
846   correlation_image=CloneImage(image,image->columns,image->rows,MagickTrue,
847     exception);
848   if (correlation_image == (Image *) NULL)
849     return(MagickFalse);
850   if (SetImageStorageClass(correlation_image,DirectClass) == MagickFalse)
851     {
852       InheritException(exception,&correlation_image->exception);
853       correlation_image=DestroyImage(correlation_image);
854       return(MagickFalse);
855     }
856   image_statistics=GetImageChannelStatistics(image,exception);
857   reconstruct_statistics=GetImageChannelStatistics(reconstruct_image,exception);
858   status=MagickTrue;
859   progress=0;
860   image_view=AcquireCacheView(image);
861   reconstruct_view=AcquireCacheView(reconstruct_image);
862   correlation_view=AcquireCacheView(correlation_image);
863 #if defined(MAGICKCORE_OPENMP_SUPPORT)
864   #pragma omp parallel for schedule(dynamic,4) shared(status,progress)
865 #endif
866   for (y=0; y < (ssize_t) image->rows; y++)
867   {
868     register const IndexPacket
869       *restrict indexes,
870       *restrict reconstruct_indexes;
871
872     register const PixelPacket
873       *restrict p,
874       *restrict q;
875
876     register IndexPacket
877       *restrict correlation_indexes;
878
879     register PixelPacket
880       *restrict c;
881
882     register ssize_t
883       x;
884
885     if (status == MagickFalse)
886       continue;
887     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
888     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
889       reconstruct_image->columns,1,exception);
890     c=GetCacheViewAuthenticPixels(correlation_view,0,y,
891       correlation_image->columns,1,exception);
892     if ((p == (const PixelPacket *) NULL) ||
893         (q == (const PixelPacket *) NULL) || (c == (PixelPacket *) NULL))
894       {
895         status=MagickFalse;
896         continue;
897       }
898     indexes=GetCacheViewVirtualIndexQueue(image_view);
899     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
900     correlation_indexes=GetCacheViewAuthenticIndexQueue(correlation_view);
901     for (x=0; x < (ssize_t) image->columns; x++)
902     {
903       if ((channel & RedChannel) != 0)
904         c->red=(p->red-image_statistics[RedChannel].mean)*
905           (q->red-reconstruct_statistics[RedChannel].mean);
906       if ((channel & GreenChannel) != 0)
907         c->green=(p->green-image_statistics[GreenChannel].mean)*
908           (q->green-reconstruct_statistics[GreenChannel].mean);
909       if ((channel & BlueChannel) != 0)
910         c->blue=(p->blue-image_statistics[BlueChannel].mean)*
911           (q->blue-reconstruct_statistics[BlueChannel].mean);
912       if (((channel & OpacityChannel) != 0) &&
913           (image->matte != MagickFalse))
914         c->opacity=(p->opacity-image_statistics[OpacityChannel].mean)*
915           (q->opacity-reconstruct_statistics[OpacityChannel].mean);
916       if (((channel & IndexChannel) != 0) &&
917           (image->colorspace == CMYKColorspace) &&
918           (reconstruct_image->colorspace == CMYKColorspace))
919         correlation_indexes[x]=(indexes[x]-
920           image_statistics[OpacityChannel].mean)*(reconstruct_indexes[x]-
921           reconstruct_statistics[OpacityChannel].mean);
922       p++;
923       q++;
924       c++;
925     }
926     if (SyncCacheViewAuthenticPixels(correlation_view,exception) == MagickFalse)
927       status=MagickFalse;
928     if (image->progress_monitor != (MagickProgressMonitor) NULL)
929       {
930         MagickBooleanType
931           proceed;
932
933 #if defined(MAGICKCORE_OPENMP_SUPPORT)
934   #pragma omp critical (MagickCore_GetNormalizedCrossCorrelationError)
935 #endif
936         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
937           image->rows);
938         if (proceed == MagickFalse)
939           status=MagickFalse;
940       }
941   }
942   correlation_view=DestroyCacheView(correlation_view);
943   reconstruct_view=DestroyCacheView(reconstruct_view);
944   image_view=DestroyCacheView(image_view);
945   correlation_statistics=GetImageChannelStatistics(correlation_image,exception);
946   correlation_image=DestroyImage(correlation_image);
947   for (i=0; i < (ssize_t) AllChannels; i++)
948     distortion[i]=correlation_statistics[i].mean/(
949       image_statistics[i].standard_deviation*
950       reconstruct_statistics[i].standard_deviation);
951   distortion[AllChannels]=0.0;
952   if ((channel & RedChannel) != 0)
953     distortion[AllChannels]+=distortion[RedChannel]*distortion[RedChannel];
954   if ((channel & GreenChannel) != 0)
955     distortion[AllChannels]+=distortion[GreenChannel]*distortion[GreenChannel];
956   if ((channel & BlueChannel) != 0)
957     distortion[AllChannels]+=distortion[BlueChannel]*distortion[BlueChannel];
958   if (((channel & OpacityChannel) != 0) && (image->matte != MagickFalse))
959     distortion[AllChannels]+=distortion[OpacityChannel]*
960       distortion[OpacityChannel];
961   if (((channel & IndexChannel) != 0) &&
962       (image->colorspace == CMYKColorspace))
963     distortion[AllChannels]+=distortion[BlackChannel]*distortion[BlackChannel];
964   distortion[AllChannels]=sqrt(distortion[AllChannels])/
965     GetNumberChannels(image,AllChannels);
966   reconstruct_statistics=(ChannelStatistics *) RelinquishMagickMemory(
967     reconstruct_statistics);
968   image_statistics=(ChannelStatistics *) RelinquishMagickMemory(
969     image_statistics);
970   correlation_statistics=(ChannelStatistics *) RelinquishMagickMemory(
971     correlation_statistics);
972   return(status);
973 }
974
975 static MagickBooleanType GetPeakAbsoluteError(const Image *image,
976   const Image *reconstruct_image,const ChannelType channel,
977   double *distortion,ExceptionInfo *exception)
978 {
979   CacheView
980     *image_view,
981     *reconstruct_view;
982
983   ssize_t
984     y;
985
986   MagickBooleanType
987     status;
988
989   status=MagickTrue;
990   image_view=AcquireCacheView(image);
991   reconstruct_view=AcquireCacheView(reconstruct_image);
992 #if defined(MAGICKCORE_OPENMP_SUPPORT)
993   #pragma omp parallel for schedule(dynamic,4) shared(status)
994 #endif
995   for (y=0; y < (ssize_t) image->rows; y++)
996   {
997     double
998       channel_distortion[AllChannels+1];
999
1000     register const IndexPacket
1001       *restrict indexes,
1002       *restrict reconstruct_indexes;
1003
1004     register const PixelPacket
1005       *restrict p,
1006       *restrict q;
1007
1008     register ssize_t
1009       i,
1010       x;
1011
1012     if (status == MagickFalse)
1013       continue;
1014     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1015     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
1016       reconstruct_image->columns,1,exception);
1017     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1018       {
1019         status=MagickFalse;
1020         continue;
1021       }
1022     indexes=GetCacheViewVirtualIndexQueue(image_view);
1023     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1024     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
1025     for (x=0; x < (ssize_t) image->columns; x++)
1026     {
1027       MagickRealType
1028         distance;
1029
1030       if ((channel & RedChannel) != 0)
1031         {
1032           distance=QuantumScale*fabs(p->red-(double) q->red);
1033           if (distance > channel_distortion[RedChannel])
1034             channel_distortion[RedChannel]=distance;
1035           if (distance > channel_distortion[AllChannels])
1036             channel_distortion[AllChannels]=distance;
1037         }
1038       if ((channel & GreenChannel) != 0)
1039         {
1040           distance=QuantumScale*fabs(p->green-(double) q->green);
1041           if (distance > channel_distortion[GreenChannel])
1042             channel_distortion[GreenChannel]=distance;
1043           if (distance > channel_distortion[AllChannels])
1044             channel_distortion[AllChannels]=distance;
1045         }
1046       if ((channel & BlueChannel) != 0)
1047         {
1048           distance=QuantumScale*fabs(p->blue-(double) q->blue);
1049           if (distance > channel_distortion[BlueChannel])
1050             channel_distortion[BlueChannel]=distance;
1051           if (distance > channel_distortion[AllChannels])
1052             channel_distortion[AllChannels]=distance;
1053         }
1054       if (((channel & OpacityChannel) != 0) &&
1055           (image->matte != MagickFalse))
1056         {
1057           distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
1058           if (distance > channel_distortion[OpacityChannel])
1059             channel_distortion[OpacityChannel]=distance;
1060           if (distance > channel_distortion[AllChannels])
1061             channel_distortion[AllChannels]=distance;
1062         }
1063       if (((channel & IndexChannel) != 0) &&
1064           (image->colorspace == CMYKColorspace) &&
1065           (reconstruct_image->colorspace == CMYKColorspace))
1066         {
1067           distance=QuantumScale*fabs(indexes[x]-(double)
1068             reconstruct_indexes[x]);
1069           if (distance > channel_distortion[BlackChannel])
1070             channel_distortion[BlackChannel]=distance;
1071           if (distance > channel_distortion[AllChannels])
1072             channel_distortion[AllChannels]=distance;
1073         }
1074       p++;
1075       q++;
1076     }
1077 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1078   #pragma omp critical (MagickCore_GetPeakAbsoluteError)
1079 #endif
1080     for (i=0; i <= (ssize_t) AllChannels; i++)
1081       if (channel_distortion[i] > distortion[i])
1082         distortion[i]=channel_distortion[i];
1083   }
1084   reconstruct_view=DestroyCacheView(reconstruct_view);
1085   image_view=DestroyCacheView(image_view);
1086   return(status);
1087 }
1088
1089 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
1090   const Image *reconstruct_image,const ChannelType channel,
1091   double *distortion,ExceptionInfo *exception)
1092 {
1093   MagickBooleanType
1094     status;
1095
1096   status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
1097     exception);
1098   if ((channel & RedChannel) != 0)
1099     distortion[RedChannel]=20.0*log10((double) 1.0/sqrt(
1100       distortion[RedChannel]));
1101   if ((channel & GreenChannel) != 0)
1102     distortion[GreenChannel]=20.0*log10((double) 1.0/sqrt(
1103       distortion[GreenChannel]));
1104   if ((channel & BlueChannel) != 0)
1105     distortion[BlueChannel]=20.0*log10((double) 1.0/sqrt(
1106       distortion[BlueChannel]));
1107   if (((channel & OpacityChannel) != 0) &&
1108       (image->matte != MagickFalse))
1109     distortion[OpacityChannel]=20.0*log10((double) 1.0/sqrt(
1110       distortion[OpacityChannel]));
1111   if (((channel & IndexChannel) != 0) &&
1112       (image->colorspace == CMYKColorspace))
1113     distortion[BlackChannel]=20.0*log10((double) 1.0/sqrt(
1114       distortion[BlackChannel]));
1115   distortion[AllChannels]=20.0*log10((double) 1.0/sqrt(
1116     distortion[AllChannels]));
1117   return(status);
1118 }
1119
1120 static MagickBooleanType GetRootMeanSquaredError(const Image *image,
1121   const Image *reconstruct_image,const ChannelType channel,
1122   double *distortion,ExceptionInfo *exception)
1123 {
1124   MagickBooleanType
1125     status;
1126
1127   status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
1128     exception);
1129   if ((channel & RedChannel) != 0)
1130     distortion[RedChannel]=sqrt(distortion[RedChannel]);
1131   if ((channel & GreenChannel) != 0)
1132     distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
1133   if ((channel & BlueChannel) != 0)
1134     distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
1135   if (((channel & OpacityChannel) != 0) &&
1136       (image->matte != MagickFalse))
1137     distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
1138   if (((channel & IndexChannel) != 0) &&
1139       (image->colorspace == CMYKColorspace))
1140     distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
1141   distortion[AllChannels]=sqrt(distortion[AllChannels]);
1142   return(status);
1143 }
1144
1145 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
1146   const Image *reconstruct_image,const ChannelType channel,
1147   const MetricType metric,double *distortion,ExceptionInfo *exception)
1148 {
1149   double
1150     *channel_distortion;
1151
1152   MagickBooleanType
1153     status;
1154
1155   size_t
1156     length;
1157
1158   assert(image != (Image *) NULL);
1159   assert(image->signature == MagickSignature);
1160   if (image->debug != MagickFalse)
1161     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1162   assert(reconstruct_image != (const Image *) NULL);
1163   assert(reconstruct_image->signature == MagickSignature);
1164   assert(distortion != (double *) NULL);
1165   *distortion=0.0;
1166   if (image->debug != MagickFalse)
1167     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1168   if ((reconstruct_image->columns != image->columns) ||
1169       (reconstruct_image->rows != image->rows))
1170     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1171   /*
1172     Get image distortion.
1173   */
1174   length=AllChannels+1UL;
1175   channel_distortion=(double *) AcquireQuantumMemory(length,
1176     sizeof(*channel_distortion));
1177   if (channel_distortion == (double *) NULL)
1178     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1179   (void) ResetMagickMemory(channel_distortion,0,length*
1180     sizeof(*channel_distortion));
1181   switch (metric)
1182   {
1183     case AbsoluteErrorMetric:
1184     {
1185       status=GetAbsoluteError(image,reconstruct_image,channel,
1186         channel_distortion,exception);
1187       break;
1188     }
1189     case MeanAbsoluteErrorMetric:
1190     {
1191       status=GetMeanAbsoluteError(image,reconstruct_image,channel,
1192         channel_distortion,exception);
1193       break;
1194     }
1195     case MeanErrorPerPixelMetric:
1196     {
1197       status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
1198         channel_distortion,exception);
1199       break;
1200     }
1201     case MeanSquaredErrorMetric:
1202     {
1203       status=GetMeanSquaredError(image,reconstruct_image,channel,
1204         channel_distortion,exception);
1205       break;
1206     }
1207     case NormalizedCrossCorrelationErrorMetric:
1208     {
1209       status=GetNormalizedCrossCorrelationError(image,reconstruct_image,channel,
1210         channel_distortion,exception);
1211       break;
1212     }
1213     case PeakAbsoluteErrorMetric:
1214     default:
1215     {
1216       status=GetPeakAbsoluteError(image,reconstruct_image,channel,
1217         channel_distortion,exception);
1218       break;
1219     }
1220     case PeakSignalToNoiseRatioMetric:
1221     {
1222       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
1223         channel_distortion,exception);
1224       break;
1225     }
1226     case RootMeanSquaredErrorMetric:
1227     {
1228       status=GetRootMeanSquaredError(image,reconstruct_image,channel,
1229         channel_distortion,exception);
1230       break;
1231     }
1232   }
1233   *distortion=channel_distortion[AllChannels];
1234   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
1235   return(status);
1236 }
1237 \f
1238 /*
1239 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1240 %                                                                             %
1241 %                                                                             %
1242 %                                                                             %
1243 %   G e t I m a g e C h a n n e l D i s t o r t i o n s                       %
1244 %                                                                             %
1245 %                                                                             %
1246 %                                                                             %
1247 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1248 %
1249 %  GetImageChannelDistrortion() compares the image channels of an image to a
1250 %  reconstructed image and returns the specified distortion metric for each
1251 %  channel.
1252 %
1253 %  The format of the CompareImageChannels method is:
1254 %
1255 %      double *GetImageChannelDistortions(const Image *image,
1256 %        const Image *reconstruct_image,const MetricType metric,
1257 %        ExceptionInfo *exception)
1258 %
1259 %  A description of each parameter follows:
1260 %
1261 %    o image: the image.
1262 %
1263 %    o reconstruct_image: the reconstruct image.
1264 %
1265 %    o metric: the metric.
1266 %
1267 %    o exception: return any errors or warnings in this structure.
1268 %
1269 */
1270 MagickExport double *GetImageChannelDistortions(Image *image,
1271   const Image *reconstruct_image,const MetricType metric,
1272   ExceptionInfo *exception)
1273 {
1274   double
1275     *channel_distortion;
1276
1277   MagickBooleanType
1278     status;
1279
1280   size_t
1281     length;
1282
1283   assert(image != (Image *) NULL);
1284   assert(image->signature == MagickSignature);
1285   if (image->debug != MagickFalse)
1286     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1287   assert(reconstruct_image != (const Image *) NULL);
1288   assert(reconstruct_image->signature == MagickSignature);
1289   if (image->debug != MagickFalse)
1290     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1291   if ((reconstruct_image->columns != image->columns) ||
1292       (reconstruct_image->rows != image->rows))
1293     {
1294       (void) ThrowMagickException(&image->exception,GetMagickModule(),
1295         ImageError,"ImageSizeDiffers","`%s'",image->filename);
1296       return((double *) NULL);
1297     }
1298   /*
1299     Get image distortion.
1300   */
1301   length=AllChannels+1UL;
1302   channel_distortion=(double *) AcquireQuantumMemory(length,
1303     sizeof(*channel_distortion));
1304   if (channel_distortion == (double *) NULL)
1305     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
1306   (void) ResetMagickMemory(channel_distortion,0,length*
1307     sizeof(*channel_distortion));
1308   switch (metric)
1309   {
1310     case AbsoluteErrorMetric:
1311     {
1312       status=GetAbsoluteError(image,reconstruct_image,AllChannels,
1313         channel_distortion,exception);
1314       break;
1315     }
1316     case MeanAbsoluteErrorMetric:
1317     {
1318       status=GetMeanAbsoluteError(image,reconstruct_image,AllChannels,
1319         channel_distortion,exception);
1320       break;
1321     }
1322     case MeanErrorPerPixelMetric:
1323     {
1324       status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
1325         channel_distortion,exception);
1326       break;
1327     }
1328     case MeanSquaredErrorMetric:
1329     {
1330       status=GetMeanSquaredError(image,reconstruct_image,AllChannels,
1331         channel_distortion,exception);
1332       break;
1333     }
1334     case NormalizedCrossCorrelationErrorMetric:
1335     {
1336       status=GetNormalizedCrossCorrelationError(image,reconstruct_image,
1337         AllChannels,channel_distortion,exception);
1338       break;
1339     }
1340     case PeakAbsoluteErrorMetric:
1341     default:
1342     {
1343       status=GetPeakAbsoluteError(image,reconstruct_image,AllChannels,
1344         channel_distortion,exception);
1345       break;
1346     }
1347     case PeakSignalToNoiseRatioMetric:
1348     {
1349       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
1350         channel_distortion,exception);
1351       break;
1352     }
1353     case RootMeanSquaredErrorMetric:
1354     {
1355       status=GetRootMeanSquaredError(image,reconstruct_image,AllChannels,
1356         channel_distortion,exception);
1357       break;
1358     }
1359   }
1360   return(channel_distortion);
1361 }
1362 \f
1363 /*
1364 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1365 %                                                                             %
1366 %                                                                             %
1367 %                                                                             %
1368 %  I s I m a g e s E q u a l                                                  %
1369 %                                                                             %
1370 %                                                                             %
1371 %                                                                             %
1372 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1373 %
1374 %  IsImagesEqual() measures the difference between colors at each pixel
1375 %  location of two images.  A value other than 0 means the colors match
1376 %  exactly.  Otherwise an error measure is computed by summing over all
1377 %  pixels in an image the distance squared in RGB space between each image
1378 %  pixel and its corresponding pixel in the reconstruct image.  The error
1379 %  measure is assigned to these image members:
1380 %
1381 %    o mean_error_per_pixel:  The mean error for any single pixel in
1382 %      the image.
1383 %
1384 %    o normalized_mean_error:  The normalized mean quantization error for
1385 %      any single pixel in the image.  This distance measure is normalized to
1386 %      a range between 0 and 1.  It is independent of the range of red, green,
1387 %      and blue values in the image.
1388 %
1389 %    o normalized_maximum_error:  The normalized maximum quantization
1390 %      error for any single pixel in the image.  This distance measure is
1391 %      normalized to a range between 0 and 1.  It is independent of the range
1392 %      of red, green, and blue values in your image.
1393 %
1394 %  A small normalized mean square error, accessed as
1395 %  image->normalized_mean_error, suggests the images are very similar in
1396 %  spatial layout and color.
1397 %
1398 %  The format of the IsImagesEqual method is:
1399 %
1400 %      MagickBooleanType IsImagesEqual(Image *image,
1401 %        const Image *reconstruct_image)
1402 %
1403 %  A description of each parameter follows.
1404 %
1405 %    o image: the image.
1406 %
1407 %    o reconstruct_image: the reconstruct image.
1408 %
1409 */
1410 MagickExport MagickBooleanType IsImagesEqual(Image *image,
1411   const Image *reconstruct_image)
1412 {
1413   CacheView
1414     *image_view,
1415     *reconstruct_view;
1416
1417   ExceptionInfo
1418     *exception;
1419
1420   ssize_t
1421     y;
1422
1423   MagickBooleanType
1424     status;
1425
1426   MagickRealType
1427     area,
1428     maximum_error,
1429     mean_error,
1430     mean_error_per_pixel;
1431
1432   assert(image != (Image *) NULL);
1433   assert(image->signature == MagickSignature);
1434   assert(reconstruct_image != (const Image *) NULL);
1435   assert(reconstruct_image->signature == MagickSignature);
1436   if ((reconstruct_image->columns != image->columns) ||
1437       (reconstruct_image->rows != image->rows))
1438     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
1439   area=0.0;
1440   maximum_error=0.0;
1441   mean_error_per_pixel=0.0;
1442   mean_error=0.0;
1443   exception=(&image->exception);
1444   image_view=AcquireCacheView(image);
1445   reconstruct_view=AcquireCacheView(reconstruct_image);
1446   for (y=0; y < (ssize_t) image->rows; y++)
1447   {
1448     register const IndexPacket
1449       *restrict indexes,
1450       *restrict reconstruct_indexes;
1451
1452     register const PixelPacket
1453       *restrict p,
1454       *restrict q;
1455
1456     register ssize_t
1457       x;
1458
1459     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1460     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
1461       1,exception);
1462     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1463       break;
1464     indexes=GetCacheViewVirtualIndexQueue(image_view);
1465     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
1466     for (x=0; x < (ssize_t) image->columns; x++)
1467     {
1468       MagickRealType
1469         distance;
1470
1471       distance=fabs(p->red-(double) q->red);
1472       mean_error_per_pixel+=distance;
1473       mean_error+=distance*distance;
1474       if (distance > maximum_error)
1475         maximum_error=distance;
1476       area++;
1477       distance=fabs(p->green-(double) q->green);
1478       mean_error_per_pixel+=distance;
1479       mean_error+=distance*distance;
1480       if (distance > maximum_error)
1481         maximum_error=distance;
1482       area++;
1483       distance=fabs(p->blue-(double) q->blue);
1484       mean_error_per_pixel+=distance;
1485       mean_error+=distance*distance;
1486       if (distance > maximum_error)
1487         maximum_error=distance;
1488       area++;
1489       if (image->matte != MagickFalse)
1490         {
1491           distance=fabs(p->opacity-(double) q->opacity);
1492           mean_error_per_pixel+=distance;
1493           mean_error+=distance*distance;
1494           if (distance > maximum_error)
1495             maximum_error=distance;
1496           area++;
1497         }
1498       if ((image->colorspace == CMYKColorspace) &&
1499           (reconstruct_image->colorspace == CMYKColorspace))
1500         {
1501           distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
1502           mean_error_per_pixel+=distance;
1503           mean_error+=distance*distance;
1504           if (distance > maximum_error)
1505             maximum_error=distance;
1506           area++;
1507         }
1508       p++;
1509       q++;
1510     }
1511   }
1512   reconstruct_view=DestroyCacheView(reconstruct_view);
1513   image_view=DestroyCacheView(image_view);
1514   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
1515   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
1516     mean_error/area);
1517   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
1518   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
1519   return(status);
1520 }
1521 \f
1522 /*
1523 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1524 %                                                                             %
1525 %                                                                             %
1526 %                                                                             %
1527 %   S i m i l a r i t y I m a g e                                             %
1528 %                                                                             %
1529 %                                                                             %
1530 %                                                                             %
1531 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1532 %
1533 %  SimilarityImage() compares the reference image of the image and returns the
1534 %  best match offset.  In addition, it returns a similarity image such that an
1535 %  exact match location is completely white and if none of the pixels match,
1536 %  black, otherwise some gray level in-between.
1537 %
1538 %  The format of the SimilarityImageImage method is:
1539 %
1540 %      Image *SimilarityImage(const Image *image,const Image *reference,
1541 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
1542 %
1543 %  A description of each parameter follows:
1544 %
1545 %    o image: the image.
1546 %
1547 %    o reference: find an area of the image that closely resembles this image.
1548 %
1549 %    o the best match offset of the reference image within the image.
1550 %
1551 %    o similarity: the computed similarity between the images.
1552 %
1553 %    o exception: return any errors or warnings in this structure.
1554 %
1555 */
1556
1557 static double GetSimilarityMetric(const Image *image,const Image *reference,
1558   const ssize_t x_offset,const ssize_t y_offset,ExceptionInfo *exception)
1559 {
1560   CacheView
1561     *image_view,
1562     *reference_view;
1563
1564   double
1565     similarity;
1566
1567   ssize_t
1568     y;
1569
1570   MagickBooleanType
1571     status;
1572
1573   /*
1574     Compute the similarity in pixels between two images.
1575   */
1576   status=MagickTrue;
1577   similarity=0.0;
1578   image_view=AcquireCacheView(image);
1579   reference_view=AcquireCacheView(reference);
1580   for (y=0; y < (ssize_t) reference->rows; y++)
1581   {
1582     register const IndexPacket
1583       *restrict indexes,
1584       *restrict reference_indexes;
1585
1586     register const PixelPacket
1587       *restrict p,
1588       *restrict q;
1589
1590     register ssize_t
1591       x;
1592
1593     if (status == MagickFalse)
1594       continue;
1595     p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset+y,
1596       reference->columns,1,exception);
1597     q=GetCacheViewVirtualPixels(reference_view,0,y,reference->columns,1,
1598       exception);
1599     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
1600       {
1601         status=MagickFalse;
1602         continue;
1603       }
1604     indexes=GetCacheViewVirtualIndexQueue(image_view);
1605     reference_indexes=GetCacheViewVirtualIndexQueue(reference_view);
1606     for (x=0; x < (ssize_t) reference->columns; x++)
1607     {
1608       double
1609         thread_similarity;
1610
1611       MagickRealType
1612         distance;
1613
1614       thread_similarity=0.0;
1615       distance=QuantumScale*(p->red-(MagickRealType) q->red);
1616       thread_similarity+=distance*distance;
1617       distance=QuantumScale*(p->green-(MagickRealType) q->green);
1618       thread_similarity+=distance*distance;
1619       distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
1620       thread_similarity+=distance*distance;
1621       if ((image->matte != MagickFalse) && (reference->matte != MagickFalse))
1622         {
1623           distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
1624           thread_similarity+=distance*distance;
1625         }
1626       if ((image->colorspace == CMYKColorspace) &&
1627           (reference->colorspace == CMYKColorspace))
1628         {
1629           distance=QuantumScale*(indexes[x]-(MagickRealType)
1630             reference_indexes[x]);
1631           thread_similarity+=distance*distance;
1632         }
1633       similarity+=thread_similarity;
1634       p++;
1635       q++;
1636     }
1637   }
1638   reference_view=DestroyCacheView(reference_view);
1639   image_view=DestroyCacheView(image_view);
1640   if (status == MagickFalse)
1641     return(0.0);
1642   similarity/=((double) reference->columns*reference->rows);
1643   similarity/=(double) GetNumberChannels(reference,AllChannels);
1644   return(sqrt(similarity));
1645 }
1646
1647 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
1648   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
1649 {
1650 #define SimilarityImageTag  "Similarity/Image"
1651
1652   CacheView
1653     *similarity_view;
1654
1655   Image
1656     *similarity_image;
1657
1658   MagickBooleanType
1659     status;
1660
1661   MagickOffsetType
1662     progress;
1663
1664   ssize_t
1665     y;
1666
1667   assert(image != (const Image *) NULL);
1668   assert(image->signature == MagickSignature);
1669   if (image->debug != MagickFalse)
1670     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1671   assert(exception != (ExceptionInfo *) NULL);
1672   assert(exception->signature == MagickSignature);
1673   assert(offset != (RectangleInfo *) NULL);
1674   SetGeometry(reference,offset);
1675   *similarity_metric=1.0;
1676   if ((reference->columns > image->columns) || (reference->rows > image->rows))
1677     ThrowImageException(ImageError,"ImageSizeDiffers");
1678   similarity_image=CloneImage(image,image->columns-reference->columns+1,
1679     image->rows-reference->rows+1,MagickTrue,exception);
1680   if (similarity_image == (Image *) NULL)
1681     return((Image *) NULL);
1682   if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
1683     {
1684       InheritException(exception,&similarity_image->exception);
1685       similarity_image=DestroyImage(similarity_image);
1686       return((Image *) NULL);
1687     }
1688   /*
1689     Measure similarity of reference image against image.
1690   */
1691   status=MagickTrue;
1692   progress=0;
1693   similarity_view=AcquireCacheView(similarity_image);
1694 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1695   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1696 #endif
1697   for (y=0; y < (ssize_t) (image->rows-reference->rows+1); y++)
1698   {
1699     double
1700       similarity;
1701
1702     register ssize_t
1703       x;
1704
1705     register PixelPacket
1706       *restrict q;
1707
1708     if (status == MagickFalse)
1709       continue;
1710     q=GetCacheViewAuthenticPixels(similarity_view,0,y,
1711       similarity_image->columns,1,exception);
1712     if (q == (const PixelPacket *) NULL)
1713       {
1714         status=MagickFalse;
1715         continue;
1716       }
1717     for (x=0; x < (ssize_t) (image->columns-reference->columns+1); x++)
1718     {
1719       similarity=GetSimilarityMetric(image,reference,x,y,exception);
1720 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1721   #pragma omp critical (MagickCore_SimilarityImage)
1722 #endif
1723       if (similarity < *similarity_metric)
1724         {
1725           *similarity_metric=similarity;
1726           offset->x=x;
1727           offset->y=y;
1728         }
1729       q->red=ClampToQuantum(QuantumRange-QuantumRange*similarity);
1730       q->green=q->red;
1731       q->blue=q->red;
1732       q++;
1733     }
1734     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
1735       status=MagickFalse;
1736     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1737       {
1738         MagickBooleanType
1739           proceed;
1740
1741 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1742   #pragma omp critical (MagickCore_SimilarityImage)
1743 #endif
1744         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
1745           image->rows);
1746         if (proceed == MagickFalse)
1747           status=MagickFalse;
1748       }
1749   }
1750   similarity_view=DestroyCacheView(similarity_view);
1751   return(similarity_image);
1752 }