compare.c

Go to the documentation of this file.
00001 /*
00002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00003 %                                                                             %
00004 %                                                                             %
00005 %                                                                             %
00006 %               CCCC   OOO   M   M  PPPP    AAA   RRRR    EEEEE               %
00007 %              C      O   O  MM MM  P   P  A   A  R   R   E                   %
00008 %              C      O   O  M M M  PPPP   AAAAA  RRRR    EEE                 %
00009 %              C      O   O  M   M  P      A   A  R R     E                   %
00010 %               CCCC   OOO   M   M  P      A   A  R  R    EEEEE               %
00011 %                                                                             %
00012 %                                                                             %
00013 %                      MagickCore Image Comparison Methods                    %
00014 %                                                                             %
00015 %                              Software Design                                %
00016 %                                John Cristy                                  %
00017 %                               December 2003                                 %
00018 %                                                                             %
00019 %                                                                             %
00020 %  Copyright 1999-2009 ImageMagick Studio LLC, a non-profit organization      %
00021 %  dedicated to making software imaging solutions freely available.           %
00022 %                                                                             %
00023 %  You may not use this file except in compliance with the License.  You may  %
00024 %  obtain a copy of the License at                                            %
00025 %                                                                             %
00026 %    http://www.imagemagick.org/script/license.php                            %
00027 %                                                                             %
00028 %  Unless required by applicable law or agreed to in writing, software        %
00029 %  distributed under the License is distributed on an "AS IS" BASIS,          %
00030 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
00031 %  See the License for the specific language governing permissions and        %
00032 %  limitations under the License.                                             %
00033 %                                                                             %
00034 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00035 %
00036 %
00037 %
00038 */
00039 
00040 /*
00041   Include declarations.
00042 */
00043 #include "magick/studio.h"
00044 #include "magick/artifact.h"
00045 #include "magick/cache-view.h"
00046 #include "magick/client.h"
00047 #include "magick/color.h"
00048 #include "magick/color-private.h"
00049 #include "magick/colorspace.h"
00050 #include "magick/colorspace-private.h"
00051 #include "magick/compare.h"
00052 #include "magick/composite-private.h"
00053 #include "magick/constitute.h"
00054 #include "magick/exception-private.h"
00055 #include "magick/geometry.h"
00056 #include "magick/image-private.h"
00057 #include "magick/list.h"
00058 #include "magick/log.h"
00059 #include "magick/memory_.h"
00060 #include "magick/monitor.h"
00061 #include "magick/monitor-private.h"
00062 #include "magick/option.h"
00063 #include "magick/pixel-private.h"
00064 #include "magick/resource_.h"
00065 #include "magick/string_.h"
00066 #include "magick/utility.h"
00067 #include "magick/version.h"
00068 
00069 /*
00070 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00071 %                                                                             %
00072 %                                                                             %
00073 %                                                                             %
00074 %   C o m p a r e I m a g e C h a n n e l s                                   %
00075 %                                                                             %
00076 %                                                                             %
00077 %                                                                             %
00078 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00079 %
00080 %  CompareImageChannels() compares one or more image channels of an image
00081 %  to a reconstructed image and returns the difference image.
00082 %
00083 %  The format of the CompareImageChannels method is:
00084 %
00085 %      Image *CompareImageChannels(const Image *image,
00086 %        const Image *reconstruct_image,const ChannelType channel,
00087 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
00088 %
00089 %  A description of each parameter follows:
00090 %
00091 %    o image: the image.
00092 %
00093 %    o reconstruct_image: the reconstruct image.
00094 %
00095 %    o channel: the channel.
00096 %
00097 %    o metric: the metric.
00098 %
00099 %    o distortion: the computed distortion between the images.
00100 %
00101 %    o exception: return any errors or warnings in this structure.
00102 %
00103 */
00104 
00105 MagickExport Image *CompareImages(Image *image,const Image *reconstruct_image,
00106   const MetricType metric,double *distortion,ExceptionInfo *exception)
00107 {
00108   Image
00109     *highlight_image;
00110 
00111   highlight_image=CompareImageChannels(image,reconstruct_image,AllChannels,
00112     metric,distortion,exception);
00113   return(highlight_image);
00114 }
00115 
00116 MagickExport Image *CompareImageChannels(Image *image,
00117   const Image *reconstruct_image,const ChannelType channel,
00118   const MetricType metric,double *distortion,ExceptionInfo *exception)
00119 {
00120   const char
00121     *artifact;
00122 
00123   Image
00124     *difference_image,
00125     *highlight_image;
00126 
00127   long
00128     y;
00129 
00130   MagickBooleanType
00131     status;
00132 
00133   MagickPixelPacket
00134     highlight,
00135     lowlight,
00136     zero;
00137 
00138   CacheView
00139     *highlight_view,
00140     *image_view,
00141     *reconstruct_view;
00142 
00143   assert(image != (Image *) NULL);
00144   assert(image->signature == MagickSignature);
00145   if (image->debug != MagickFalse)
00146     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00147   assert(reconstruct_image != (const Image *) NULL);
00148   assert(reconstruct_image->signature == MagickSignature);
00149   assert(distortion != (double *) NULL);
00150   *distortion=0.0;
00151   if (image->debug != MagickFalse)
00152     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00153   if ((reconstruct_image->columns != image->columns) ||
00154       (reconstruct_image->rows != image->rows))
00155     ThrowImageException(ImageError,"ImageSizeDiffers");
00156   status=GetImageChannelDistortion(image,reconstruct_image,channel,metric,
00157     distortion,exception);
00158   if (status == MagickFalse)
00159     return((Image *) NULL);
00160   difference_image=CloneImage(image,0,0,MagickTrue,exception);
00161   if (difference_image == (Image *) NULL)
00162     return((Image *) NULL);
00163   (void) SetImageAlphaChannel(difference_image,OpaqueAlphaChannel);
00164   highlight_image=CloneImage(image,image->columns,image->rows,MagickTrue,
00165     exception);
00166   if (highlight_image == (Image *) NULL)
00167     {
00168       difference_image=DestroyImage(difference_image);
00169       return((Image *) NULL);
00170     }
00171   if (SetImageStorageClass(highlight_image,DirectClass) == MagickFalse)
00172     {
00173       InheritException(exception,&highlight_image->exception);
00174       difference_image=DestroyImage(difference_image);
00175       highlight_image=DestroyImage(highlight_image);
00176       return((Image *) NULL);
00177     }
00178   (void) SetImageAlphaChannel(highlight_image,OpaqueAlphaChannel);
00179   (void) QueryMagickColor("#f1001ecc",&highlight,exception);
00180   artifact=GetImageArtifact(image,"highlight-color");
00181   if (artifact != (const char *) NULL)
00182     (void) QueryMagickColor(artifact,&highlight,exception);
00183   (void) QueryMagickColor("#ffffffcc",&lowlight,exception);
00184   artifact=GetImageArtifact(image,"lowlight-color");
00185   if (artifact != (const char *) NULL)
00186     (void) QueryMagickColor(artifact,&lowlight,exception);
00187   if (highlight_image->colorspace == CMYKColorspace)
00188     {
00189       ConvertRGBToCMYK(&highlight);
00190       ConvertRGBToCMYK(&lowlight);
00191     }
00192   /*
00193     Generate difference image.
00194   */
00195   status=MagickTrue;
00196   GetMagickPixelPacket(image,&zero);
00197   image_view=AcquireCacheView(image);
00198   reconstruct_view=AcquireCacheView(reconstruct_image);
00199   highlight_view=AcquireCacheView(highlight_image);
00200 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00201   #pragma omp parallel for schedule(dynamic,4) shared(status)
00202 #endif
00203   for (y=0; y < (long) image->rows; y++)
00204   {
00205     MagickBooleanType
00206       sync;
00207 
00208     MagickPixelPacket
00209       pixel,
00210       reconstruct_pixel;
00211 
00212     register const IndexPacket
00213       *__restrict indexes,
00214       *__restrict reconstruct_indexes;
00215 
00216     register const PixelPacket
00217       *__restrict p,
00218       *__restrict q;
00219 
00220     register IndexPacket
00221       *__restrict highlight_indexes;
00222 
00223     register long
00224       x;
00225 
00226     register PixelPacket
00227       *__restrict r;
00228 
00229     if (status == MagickFalse)
00230       continue;
00231     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
00232     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
00233       1,exception);
00234     r=QueueCacheViewAuthenticPixels(highlight_view,0,y,highlight_image->columns,
00235       1,exception);
00236     if ((p == (const PixelPacket *) NULL) ||
00237         (q == (const PixelPacket *) NULL) || (r == (PixelPacket *) NULL))
00238       {
00239         status=MagickFalse;
00240         continue;
00241       }
00242     indexes=GetCacheViewVirtualIndexQueue(image_view);
00243     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
00244     highlight_indexes=GetCacheViewAuthenticIndexQueue(highlight_view);
00245     pixel=zero;
00246     reconstruct_pixel=zero;
00247     for (x=0; x < (long) image->columns; x++)
00248     {
00249       MagickStatusType
00250         difference;
00251 
00252       SetMagickPixelPacket(image,p,indexes+x,&pixel);
00253       SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
00254         &reconstruct_pixel);
00255       difference=MagickFalse;
00256       if (channel == AllChannels)
00257         {
00258           if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
00259             difference=MagickTrue;
00260         }
00261       else
00262         {
00263           if (((channel & RedChannel) != 0) && (p->red != q->red))
00264             difference=MagickTrue;
00265           if (((channel & GreenChannel) != 0) && (p->green != q->green))
00266             difference=MagickTrue;
00267           if (((channel & BlueChannel) != 0) && (p->blue != q->blue))
00268             difference=MagickTrue;
00269           if (((channel & OpacityChannel) != 0) &&
00270               (image->matte != MagickFalse) && (p->opacity != q->opacity))
00271             difference=MagickTrue;
00272           if ((((channel & IndexChannel) != 0) &&
00273                (image->colorspace == CMYKColorspace) &&
00274                (reconstruct_image->colorspace == CMYKColorspace)) &&
00275               (indexes[x] != reconstruct_indexes[x]))
00276             difference=MagickTrue;
00277         }
00278       if (difference != MagickFalse)
00279         SetPixelPacket(highlight_image,&highlight,r,highlight_indexes+x);
00280       else
00281         SetPixelPacket(highlight_image,&lowlight,r,highlight_indexes+x);
00282       p++;
00283       q++;
00284       r++;
00285     }
00286     sync=SyncCacheViewAuthenticPixels(highlight_view,exception);
00287     if (sync == MagickFalse)
00288       status=MagickFalse;
00289   }
00290   highlight_view=DestroyCacheView(highlight_view);
00291   reconstruct_view=DestroyCacheView(reconstruct_view);
00292   image_view=DestroyCacheView(image_view);
00293   (void) CompositeImage(difference_image,image->compose,highlight_image,0,0);
00294   highlight_image=DestroyImage(highlight_image);
00295   if (status == MagickFalse)
00296     difference_image=DestroyImage(difference_image);
00297   return(difference_image);
00298 }
00299 
00300 /*
00301 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00302 %                                                                             %
00303 %                                                                             %
00304 %                                                                             %
00305 %   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                         %
00306 %                                                                             %
00307 %                                                                             %
00308 %                                                                             %
00309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00310 %
00311 %  GetImageChannelDistortion() compares one or more image channels of an image
00312 %  to a reconstructed image and returns the specified distortion metric.
00313 %
00314 %  The format of the CompareImageChannels method is:
00315 %
00316 %      MagickBooleanType GetImageChannelDistortion(const Image *image,
00317 %        const Image *reconstruct_image,const ChannelType channel,
00318 %        const MetricType metric,double *distortion,ExceptionInfo *exception)
00319 %
00320 %  A description of each parameter follows:
00321 %
00322 %    o image: the image.
00323 %
00324 %    o reconstruct_image: the reconstruct image.
00325 %
00326 %    o channel: the channel.
00327 %
00328 %    o metric: the metric.
00329 %
00330 %    o distortion: the computed distortion between the images.
00331 %
00332 %    o exception: return any errors or warnings in this structure.
00333 %
00334 */
00335 
00336 MagickExport MagickBooleanType GetImageDistortion(Image *image,
00337   const Image *reconstruct_image,const MetricType metric,double *distortion,
00338   ExceptionInfo *exception)
00339 {
00340   MagickBooleanType
00341     status;
00342 
00343   status=GetImageChannelDistortion(image,reconstruct_image,AllChannels,
00344     metric,distortion,exception);
00345   return(status);
00346 }
00347 
00348 static MagickBooleanType GetAbsoluteError(const Image *image,
00349   const Image *reconstruct_image,const ChannelType channel,double *distortion,
00350   ExceptionInfo *exception)
00351 {
00352   long
00353     y;
00354 
00355   MagickBooleanType
00356     status;
00357 
00358   MagickPixelPacket
00359     zero;
00360 
00361   CacheView
00362     *image_view,
00363     *reconstruct_view;
00364 
00365   /*
00366     Compute the absolute difference in pixels between two images.
00367   */
00368   status=MagickTrue;
00369   GetMagickPixelPacket(image,&zero);
00370   image_view=AcquireCacheView(image);
00371   reconstruct_view=AcquireCacheView(reconstruct_image);
00372 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00373   #pragma omp parallel for schedule(dynamic,4) shared(status)
00374 #endif
00375   for (y=0; y < (long) image->rows; y++)
00376   {
00377     double
00378       channel_distortion[AllChannels+1];
00379 
00380     MagickPixelPacket
00381       pixel,
00382       reconstruct_pixel;
00383 
00384     register const IndexPacket
00385       *__restrict indexes,
00386       *__restrict reconstruct_indexes;
00387 
00388     register const PixelPacket
00389       *__restrict p,
00390       *__restrict q;
00391 
00392     register long
00393       i,
00394       x;
00395 
00396     if (status == MagickFalse)
00397       continue;
00398     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
00399     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
00400       1,exception);
00401     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
00402       {
00403         status=MagickFalse;
00404         continue;
00405       }
00406     indexes=GetCacheViewVirtualIndexQueue(image_view);
00407     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
00408     pixel=zero;
00409     reconstruct_pixel=pixel;
00410     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
00411     for (x=0; x < (long) image->columns; x++)
00412     {
00413       SetMagickPixelPacket(image,p,indexes+x,&pixel);
00414       SetMagickPixelPacket(reconstruct_image,q,reconstruct_indexes+x,
00415         &reconstruct_pixel);
00416       if (IsMagickColorSimilar(&pixel,&reconstruct_pixel) == MagickFalse)
00417         {
00418           if ((channel & RedChannel) != 0)
00419             channel_distortion[RedChannel]++;
00420           if ((channel & GreenChannel) != 0)
00421             channel_distortion[GreenChannel]++;
00422           if ((channel & BlueChannel) != 0)
00423             channel_distortion[BlueChannel]++;
00424           if (((channel & OpacityChannel) != 0) &&
00425               (image->matte != MagickFalse))
00426             channel_distortion[OpacityChannel]++;
00427           if (((channel & IndexChannel) != 0) &&
00428               (image->colorspace == CMYKColorspace))
00429             channel_distortion[BlackChannel]++;
00430           channel_distortion[AllChannels]++;
00431         }
00432       p++;
00433       q++;
00434     }
00435 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00436   #pragma omp critical (MagickCore_GetAbsoluteError)
00437 #endif
00438     for (i=0; i <= (long) AllChannels; i++)
00439       distortion[i]+=channel_distortion[i];
00440   }
00441   reconstruct_view=DestroyCacheView(reconstruct_view);
00442   image_view=DestroyCacheView(image_view);
00443   return(status);
00444 }
00445 
00446 static unsigned long GetNumberChannels(const Image *image,
00447   const ChannelType channel)
00448 {
00449   unsigned long
00450     channels;
00451 
00452   channels=0;
00453   if ((channel & RedChannel) != 0)
00454     channels++;
00455   if ((channel & GreenChannel) != 0)
00456     channels++;
00457   if ((channel & BlueChannel) != 0)
00458     channels++;
00459   if (((channel & OpacityChannel) != 0) &&
00460        (image->matte != MagickFalse))
00461     channels++;
00462   if (((channel & IndexChannel) != 0) &&
00463       (image->colorspace == CMYKColorspace))
00464     channels++;
00465   return(channels);
00466 }
00467 
00468 static MagickBooleanType GetMeanAbsoluteError(const Image *image,
00469   const Image *reconstruct_image,const ChannelType channel,
00470   double *distortion,ExceptionInfo *exception)
00471 {
00472   long
00473     y;
00474 
00475   MagickBooleanType
00476     status;
00477 
00478   register long
00479     i;
00480 
00481   CacheView
00482     *image_view,
00483     *reconstruct_view;
00484 
00485   status=MagickTrue;
00486   image_view=AcquireCacheView(image);
00487   reconstruct_view=AcquireCacheView(reconstruct_image);
00488 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00489   #pragma omp parallel for schedule(dynamic,4) shared(status)
00490 #endif
00491   for (y=0; y < (long) image->rows; y++)
00492   {
00493     double
00494       channel_distortion[AllChannels+1];
00495 
00496     register const IndexPacket
00497       *__restrict indexes,
00498       *__restrict reconstruct_indexes;
00499 
00500     register const PixelPacket
00501       *__restrict p,
00502       *__restrict q;
00503 
00504     register long
00505       i,
00506       x;
00507 
00508     if (status == MagickFalse)
00509       continue;
00510     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
00511     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
00512       reconstruct_image->columns,1,exception);
00513     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
00514       {
00515         status=MagickFalse;
00516         continue;
00517       }
00518     indexes=GetCacheViewVirtualIndexQueue(image_view);
00519     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
00520     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
00521     for (x=0; x < (long) image->columns; x++)
00522     {
00523       MagickRealType
00524         distance;
00525 
00526       if ((channel & RedChannel) != 0)
00527         {
00528           distance=QuantumScale*fabs(p->red-(double) q->red);
00529           channel_distortion[RedChannel]+=distance;
00530           channel_distortion[AllChannels]+=distance;
00531         }
00532       if ((channel & GreenChannel) != 0)
00533         {
00534           distance=QuantumScale*fabs(p->green-(double) q->green);
00535           channel_distortion[GreenChannel]+=distance;
00536           channel_distortion[AllChannels]+=distance;
00537         }
00538       if ((channel & BlueChannel) != 0)
00539         {
00540           distance=QuantumScale*fabs(p->blue-(double) q->blue);
00541           channel_distortion[BlueChannel]+=distance;
00542           channel_distortion[AllChannels]+=distance;
00543         }
00544       if (((channel & OpacityChannel) != 0) &&
00545           (image->matte != MagickFalse))
00546         {
00547           distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
00548           channel_distortion[OpacityChannel]+=distance;
00549           channel_distortion[AllChannels]+=distance;
00550         }
00551       if (((channel & IndexChannel) != 0) &&
00552           (image->colorspace == CMYKColorspace))
00553         {
00554           distance=QuantumScale*fabs(indexes[x]-(double)
00555             reconstruct_indexes[x]);
00556           channel_distortion[BlackChannel]+=distance;
00557           channel_distortion[AllChannels]+=distance;
00558         }
00559       p++;
00560       q++;
00561     }
00562 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00563   #pragma omp critical (MagickCore_GetMeanAbsoluteError)
00564 #endif
00565     for (i=0; i <= (long) AllChannels; i++)
00566       distortion[i]+=channel_distortion[i];
00567   }
00568   reconstruct_view=DestroyCacheView(reconstruct_view);
00569   image_view=DestroyCacheView(image_view);
00570   for (i=0; i <= (long) AllChannels; i++)
00571     distortion[i]/=((double) image->columns*image->rows);
00572   distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
00573   return(status);
00574 }
00575 
00576 static MagickBooleanType GetMeanErrorPerPixel(Image *image,
00577   const Image *reconstruct_image,const ChannelType channel,double *distortion,
00578   ExceptionInfo *exception)
00579 {
00580   long
00581     y;
00582 
00583   MagickBooleanType
00584     status;
00585 
00586   MagickRealType
00587     alpha,
00588     area,
00589     beta,
00590     maximum_error,
00591     mean_error;
00592 
00593   CacheView
00594     *image_view,
00595     *reconstruct_view;
00596 
00597   status=MagickTrue;
00598   alpha=1.0;
00599   beta=1.0;
00600   area=0.0;
00601   maximum_error=0.0;
00602   mean_error=0.0;
00603   image_view=AcquireCacheView(image);
00604   reconstruct_view=AcquireCacheView(reconstruct_image);
00605   for (y=0; y < (long) image->rows; y++)
00606   {
00607     register const IndexPacket
00608       *__restrict indexes,
00609       *__restrict reconstruct_indexes;
00610 
00611     register const PixelPacket
00612       *__restrict p,
00613       *__restrict q;
00614 
00615     register long
00616       x;
00617 
00618     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
00619     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
00620       1,exception);
00621     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
00622       {
00623         status=MagickFalse;
00624         break;
00625       }
00626     indexes=GetCacheViewVirtualIndexQueue(image_view);
00627     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
00628     for (x=0; x < (long) image->columns; x++)
00629     {
00630       MagickRealType
00631         distance;
00632 
00633       if ((channel & OpacityChannel) != 0)
00634         {
00635           if (image->matte != MagickFalse)
00636             alpha=(MagickRealType) (QuantumScale*(QuantumRange-p->opacity));
00637           if (reconstruct_image->matte != MagickFalse)
00638             beta=(MagickRealType) (QuantumScale*(QuantumRange-q->opacity));
00639         }
00640       if ((channel & RedChannel) != 0)
00641         {
00642           distance=fabs(alpha*p->red-beta*q->red);
00643           distortion[RedChannel]+=distance;
00644           distortion[AllChannels]+=distance;
00645           mean_error+=distance*distance;
00646           if (distance > maximum_error)
00647             maximum_error=distance;
00648           area++;
00649         }
00650       if ((channel & GreenChannel) != 0)
00651         {
00652           distance=fabs(alpha*p->green-beta*q->green);
00653           distortion[GreenChannel]+=distance;
00654           distortion[AllChannels]+=distance;
00655           mean_error+=distance*distance;
00656           if (distance > maximum_error)
00657             maximum_error=distance;
00658           area++;
00659         }
00660       if ((channel & BlueChannel) != 0)
00661         {
00662           distance=fabs(alpha*p->blue-beta*q->blue);
00663           distortion[BlueChannel]+=distance;
00664           distortion[AllChannels]+=distance;
00665           mean_error+=distance*distance;
00666           if (distance > maximum_error)
00667             maximum_error=distance;
00668           area++;
00669         }
00670       if (((channel & OpacityChannel) != 0) &&
00671           (image->matte != MagickFalse))
00672         {
00673           distance=fabs((double) p->opacity-q->opacity);
00674           distortion[OpacityChannel]+=distance;
00675           distortion[AllChannels]+=distance;
00676           mean_error+=distance*distance;
00677           if (distance > maximum_error)
00678             maximum_error=distance;
00679           area++;
00680         }
00681       if (((channel & IndexChannel) != 0) &&
00682           (image->colorspace == CMYKColorspace) &&
00683           (reconstruct_image->colorspace == CMYKColorspace))
00684         {
00685           distance=fabs(alpha*indexes[x]-beta*reconstruct_indexes[x]);
00686           distortion[BlackChannel]+=distance;
00687           distortion[AllChannels]+=distance;
00688           mean_error+=distance*distance;
00689           if (distance > maximum_error)
00690             maximum_error=distance;
00691           area++;
00692         }
00693       p++;
00694       q++;
00695     }
00696   }
00697   reconstruct_view=DestroyCacheView(reconstruct_view);
00698   image_view=DestroyCacheView(image_view);
00699   image->error.mean_error_per_pixel=distortion[AllChannels]/area;
00700   image->error.normalized_mean_error=QuantumScale*QuantumScale*mean_error/area;
00701   image->error.normalized_maximum_error=QuantumScale*maximum_error;
00702   return(status);
00703 }
00704 
00705 static MagickBooleanType GetMeanSquaredError(const Image *image,
00706   const Image *reconstruct_image,const ChannelType channel,
00707   double *distortion,ExceptionInfo *exception)
00708 {
00709   long
00710     y;
00711 
00712   MagickBooleanType
00713     status;
00714 
00715   register long
00716     i;
00717 
00718   CacheView
00719     *image_view,
00720     *reconstruct_view;
00721 
00722   status=MagickTrue;
00723   image_view=AcquireCacheView(image);
00724   reconstruct_view=AcquireCacheView(reconstruct_image);
00725 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00726   #pragma omp parallel for schedule(dynamic,4) shared(status)
00727 #endif
00728   for (y=0; y < (long) image->rows; y++)
00729   {
00730     double
00731       channel_distortion[AllChannels+1];
00732 
00733     register const IndexPacket
00734       *__restrict indexes,
00735       *__restrict reconstruct_indexes;
00736 
00737     register const PixelPacket
00738       *__restrict p,
00739       *__restrict q;
00740 
00741     register long
00742       i,
00743       x;
00744 
00745     if (status == MagickFalse)
00746       continue;
00747     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
00748     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
00749       reconstruct_image->columns,1,exception);
00750     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
00751       {
00752         status=MagickFalse;
00753         continue;
00754       }
00755     indexes=GetCacheViewVirtualIndexQueue(image_view);
00756     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
00757     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
00758     for (x=0; x < (long) image->columns; x++)
00759     {
00760       MagickRealType
00761         distance;
00762 
00763       if ((channel & RedChannel) != 0)
00764         {
00765           distance=QuantumScale*(p->red-(MagickRealType) q->red);
00766           channel_distortion[RedChannel]+=distance*distance;
00767           channel_distortion[AllChannels]+=distance*distance;
00768         }
00769       if ((channel & GreenChannel) != 0)
00770         {
00771           distance=QuantumScale*(p->green-(MagickRealType) q->green);
00772           channel_distortion[GreenChannel]+=distance*distance;
00773           channel_distortion[AllChannels]+=distance*distance;
00774         }
00775       if ((channel & BlueChannel) != 0)
00776         {
00777           distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
00778           channel_distortion[BlueChannel]+=distance*distance;
00779           channel_distortion[AllChannels]+=distance*distance;
00780         }
00781       if (((channel & OpacityChannel) != 0) &&
00782           (image->matte != MagickFalse))
00783         {
00784           distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
00785           channel_distortion[OpacityChannel]+=distance*distance;
00786           channel_distortion[AllChannels]+=distance*distance;
00787         }
00788       if (((channel & IndexChannel) != 0) &&
00789           (image->colorspace == CMYKColorspace) &&
00790           (reconstruct_image->colorspace == CMYKColorspace))
00791         {
00792           distance=QuantumScale*(indexes[x]-(MagickRealType)
00793             reconstruct_indexes[x]);
00794           channel_distortion[BlackChannel]+=distance*distance;
00795           channel_distortion[AllChannels]+=distance*distance;
00796         }
00797       p++;
00798       q++;
00799     }
00800 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00801   #pragma omp critical (MagickCore_GetMeanSquaredError)
00802 #endif
00803     for (i=0; i <= (long) AllChannels; i++)
00804       distortion[i]+=channel_distortion[i];
00805   }
00806   reconstruct_view=DestroyCacheView(reconstruct_view);
00807   image_view=DestroyCacheView(image_view);
00808   for (i=0; i <= (long) AllChannels; i++)
00809     distortion[i]/=((double) image->columns*image->rows);
00810   distortion[AllChannels]/=(double) GetNumberChannels(image,channel);
00811   return(status);
00812 }
00813 
00814 static MagickBooleanType GetPeakAbsoluteError(const Image *image,
00815   const Image *reconstruct_image,const ChannelType channel,
00816   double *distortion,ExceptionInfo *exception)
00817 {
00818   long
00819     y;
00820 
00821   MagickBooleanType
00822     status;
00823 
00824   CacheView
00825     *image_view,
00826     *reconstruct_view;
00827 
00828   status=MagickTrue;
00829   image_view=AcquireCacheView(image);
00830   reconstruct_view=AcquireCacheView(reconstruct_image);
00831 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00832   #pragma omp parallel for schedule(dynamic,4) shared(status)
00833 #endif
00834   for (y=0; y < (long) image->rows; y++)
00835   {
00836     double
00837       channel_distortion[AllChannels+1];
00838 
00839     register const IndexPacket
00840       *__restrict indexes,
00841       *__restrict reconstruct_indexes;
00842 
00843     register const PixelPacket
00844       *__restrict p,
00845       *__restrict q;
00846 
00847     register long
00848       i,
00849       x;
00850 
00851     if (status == MagickFalse)
00852       continue;
00853     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
00854     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,
00855       reconstruct_image->columns,1,exception);
00856     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
00857       {
00858         status=MagickFalse;
00859         continue;
00860       }
00861     indexes=GetCacheViewVirtualIndexQueue(image_view);
00862     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
00863     (void) ResetMagickMemory(channel_distortion,0,sizeof(channel_distortion));
00864     for (x=0; x < (long) image->columns; x++)
00865     {
00866       MagickRealType
00867         distance;
00868 
00869       if ((channel & RedChannel) != 0)
00870         {
00871           distance=QuantumScale*fabs(p->red-(double) q->red);
00872           if (distance > channel_distortion[RedChannel])
00873             channel_distortion[RedChannel]=distance;
00874           if (distance > channel_distortion[AllChannels])
00875             channel_distortion[AllChannels]=distance;
00876         }
00877       if ((channel & GreenChannel) != 0)
00878         {
00879           distance=QuantumScale*fabs(p->green-(double) q->green);
00880           if (distance > channel_distortion[GreenChannel])
00881             channel_distortion[GreenChannel]=distance;
00882           if (distance > channel_distortion[AllChannels])
00883             channel_distortion[AllChannels]=distance;
00884         }
00885       if ((channel & BlueChannel) != 0)
00886         {
00887           distance=QuantumScale*fabs(p->blue-(double) q->blue);
00888           if (distance > channel_distortion[BlueChannel])
00889             channel_distortion[BlueChannel]=distance;
00890           if (distance > channel_distortion[AllChannels])
00891             channel_distortion[AllChannels]=distance;
00892         }
00893       if (((channel & OpacityChannel) != 0) &&
00894           (image->matte != MagickFalse))
00895         {
00896           distance=QuantumScale*fabs(p->opacity-(double) q->opacity);
00897           if (distance > channel_distortion[OpacityChannel])
00898             channel_distortion[OpacityChannel]=distance;
00899           if (distance > channel_distortion[AllChannels])
00900             channel_distortion[AllChannels]=distance;
00901         }
00902       if (((channel & IndexChannel) != 0) &&
00903           (image->colorspace == CMYKColorspace) &&
00904           (reconstruct_image->colorspace == CMYKColorspace))
00905         {
00906           distance=QuantumScale*fabs(indexes[x]-(double)
00907             reconstruct_indexes[x]);
00908           if (distance > channel_distortion[BlackChannel])
00909             channel_distortion[BlackChannel]=distance;
00910           if (distance > channel_distortion[AllChannels])
00911             channel_distortion[AllChannels]=distance;
00912         }
00913       p++;
00914       q++;
00915     }
00916 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00917   #pragma omp critical (MagickCore_GetPeakAbsoluteError)
00918 #endif
00919     for (i=0; i <= (long) AllChannels; i++)
00920       if (channel_distortion[i] > distortion[i])
00921         distortion[i]=channel_distortion[i];
00922   }
00923   reconstruct_view=DestroyCacheView(reconstruct_view);
00924   image_view=DestroyCacheView(image_view);
00925   return(status);
00926 }
00927 
00928 static MagickBooleanType GetPeakSignalToNoiseRatio(const Image *image,
00929   const Image *reconstruct_image,const ChannelType channel,
00930   double *distortion,ExceptionInfo *exception)
00931 {
00932   MagickBooleanType
00933     status;
00934 
00935   status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
00936     exception);
00937   if ((channel & RedChannel) != 0)
00938     distortion[RedChannel]=20.0*log10((double) 1.0/sqrt(
00939       distortion[RedChannel]));
00940   if ((channel & GreenChannel) != 0)
00941     distortion[GreenChannel]=20.0*log10((double) 1.0/sqrt(
00942       distortion[GreenChannel]));
00943   if ((channel & BlueChannel) != 0)
00944     distortion[BlueChannel]=20.0*log10((double) 1.0/sqrt(
00945       distortion[BlueChannel]));
00946   if (((channel & OpacityChannel) != 0) &&
00947       (image->matte != MagickFalse))
00948     distortion[OpacityChannel]=20.0*log10((double) 1.0/sqrt(
00949       distortion[OpacityChannel]));
00950   if (((channel & IndexChannel) != 0) &&
00951       (image->colorspace == CMYKColorspace))
00952     distortion[BlackChannel]=20.0*log10((double) 1.0/sqrt(
00953       distortion[BlackChannel]));
00954   distortion[AllChannels]=20.0*log10((double) 1.0/sqrt(
00955     distortion[AllChannels]));
00956   return(status);
00957 }
00958 
00959 static MagickBooleanType GetRootMeanSquaredError(const Image *image,
00960   const Image *reconstruct_image,const ChannelType channel,
00961   double *distortion,ExceptionInfo *exception)
00962 {
00963   MagickBooleanType
00964     status;
00965 
00966   status=GetMeanSquaredError(image,reconstruct_image,channel,distortion,
00967     exception);
00968   if ((channel & RedChannel) != 0)
00969     distortion[RedChannel]=sqrt(distortion[RedChannel]);
00970   if ((channel & GreenChannel) != 0)
00971     distortion[GreenChannel]=sqrt(distortion[GreenChannel]);
00972   if ((channel & BlueChannel) != 0)
00973     distortion[BlueChannel]=sqrt(distortion[BlueChannel]);
00974   if (((channel & OpacityChannel) != 0) &&
00975       (image->matte != MagickFalse))
00976     distortion[OpacityChannel]=sqrt(distortion[OpacityChannel]);
00977   if (((channel & IndexChannel) != 0) &&
00978       (image->colorspace == CMYKColorspace))
00979     distortion[BlackChannel]=sqrt(distortion[BlackChannel]);
00980   distortion[AllChannels]=sqrt(distortion[AllChannels]);
00981   return(status);
00982 }
00983 
00984 MagickExport MagickBooleanType GetImageChannelDistortion(Image *image,
00985   const Image *reconstruct_image,const ChannelType channel,
00986   const MetricType metric,double *distortion,ExceptionInfo *exception)
00987 {
00988   double
00989     *channel_distortion;
00990 
00991   MagickBooleanType
00992     status;
00993 
00994   size_t
00995     length;
00996 
00997   assert(image != (Image *) NULL);
00998   assert(image->signature == MagickSignature);
00999   if (image->debug != MagickFalse)
01000     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01001   assert(reconstruct_image != (const Image *) NULL);
01002   assert(reconstruct_image->signature == MagickSignature);
01003   assert(distortion != (double *) NULL);
01004   *distortion=0.0;
01005   if (image->debug != MagickFalse)
01006     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01007   if ((reconstruct_image->columns != image->columns) ||
01008       (reconstruct_image->rows != image->rows))
01009     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
01010   /*
01011     Get image distortion.
01012   */
01013   length=AllChannels+1UL;
01014   channel_distortion=(double *) AcquireQuantumMemory(length,
01015     sizeof(*channel_distortion));
01016   if (channel_distortion == (double *) NULL)
01017     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
01018   (void) ResetMagickMemory(channel_distortion,0,length*
01019     sizeof(*channel_distortion));
01020   switch (metric)
01021   {
01022     case AbsoluteErrorMetric:
01023     {
01024       status=GetAbsoluteError(image,reconstruct_image,channel,
01025         channel_distortion,exception);
01026       break;
01027     }
01028     case MeanAbsoluteErrorMetric:
01029     {
01030       status=GetMeanAbsoluteError(image,reconstruct_image,channel,
01031         channel_distortion,exception);
01032       break;
01033     }
01034     case MeanErrorPerPixelMetric:
01035     {
01036       status=GetMeanErrorPerPixel(image,reconstruct_image,channel,
01037         channel_distortion,exception);
01038       break;
01039     }
01040     case MeanSquaredErrorMetric:
01041     {
01042       status=GetMeanSquaredError(image,reconstruct_image,channel,
01043         channel_distortion,exception);
01044       break;
01045     }
01046     case PeakAbsoluteErrorMetric:
01047     default:
01048     {
01049       status=GetPeakAbsoluteError(image,reconstruct_image,channel,
01050         channel_distortion,exception);
01051       break;
01052     }
01053     case PeakSignalToNoiseRatioMetric:
01054     {
01055       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,channel,
01056         channel_distortion,exception);
01057       break;
01058     }
01059     case RootMeanSquaredErrorMetric:
01060     {
01061       status=GetRootMeanSquaredError(image,reconstruct_image,channel,
01062         channel_distortion,exception);
01063       break;
01064     }
01065   }
01066   *distortion=channel_distortion[AllChannels];
01067   channel_distortion=(double *) RelinquishMagickMemory(channel_distortion);
01068   return(status);
01069 }
01070 
01071 /*
01072 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01073 %                                                                             %
01074 %                                                                             %
01075 %                                                                             %
01076 %   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                       %
01077 %                                                                             %
01078 %                                                                             %
01079 %                                                                             %
01080 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01081 %
01082 %  GetImageChannelDistrortion() compares the image channels of an image to a
01083 %  reconstructed image and returns the specified distortion metric for each
01084 %  channel.
01085 %
01086 %  The format of the CompareImageChannels method is:
01087 %
01088 %      double *GetImageChannelDistortions(const Image *image,
01089 %        const Image *reconstruct_image,const MetricType metric,
01090 %        ExceptionInfo *exception)
01091 %
01092 %  A description of each parameter follows:
01093 %
01094 %    o image: the image.
01095 %
01096 %    o reconstruct_image: the reconstruct image.
01097 %
01098 %    o metric: the metric.
01099 %
01100 %    o exception: return any errors or warnings in this structure.
01101 %
01102 */
01103 MagickExport double *GetImageChannelDistortions(Image *image,
01104   const Image *reconstruct_image,const MetricType metric,
01105   ExceptionInfo *exception)
01106 {
01107   double
01108     *channel_distortion;
01109 
01110   MagickBooleanType
01111     status;
01112 
01113   size_t
01114     length;
01115 
01116   assert(image != (Image *) NULL);
01117   assert(image->signature == MagickSignature);
01118   if (image->debug != MagickFalse)
01119     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01120   assert(reconstruct_image != (const Image *) NULL);
01121   assert(reconstruct_image->signature == MagickSignature);
01122   if (image->debug != MagickFalse)
01123     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01124   if ((reconstruct_image->columns != image->columns) ||
01125       (reconstruct_image->rows != image->rows))
01126     {
01127       (void) ThrowMagickException(&image->exception,GetMagickModule(),
01128         ImageError,"ImageSizeDiffers","`%s'",image->filename);
01129       return((double *) NULL);
01130     }
01131   /*
01132     Get image distortion.
01133   */
01134   length=AllChannels+1UL;
01135   channel_distortion=(double *) AcquireQuantumMemory(length,
01136     sizeof(*channel_distortion));
01137   if (channel_distortion == (double *) NULL)
01138     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
01139   (void) ResetMagickMemory(channel_distortion,0,length*
01140     sizeof(*channel_distortion));
01141   switch (metric)
01142   {
01143     case AbsoluteErrorMetric:
01144     {
01145       status=GetAbsoluteError(image,reconstruct_image,AllChannels,
01146         channel_distortion,exception);
01147       break;
01148     }
01149     case MeanAbsoluteErrorMetric:
01150     {
01151       status=GetMeanAbsoluteError(image,reconstruct_image,AllChannels,
01152         channel_distortion,exception);
01153       break;
01154     }
01155     case MeanErrorPerPixelMetric:
01156     {
01157       status=GetMeanErrorPerPixel(image,reconstruct_image,AllChannels,
01158         channel_distortion,exception);
01159       break;
01160     }
01161     case MeanSquaredErrorMetric:
01162     {
01163       status=GetMeanSquaredError(image,reconstruct_image,AllChannels,
01164         channel_distortion,exception);
01165       break;
01166     }
01167     case PeakAbsoluteErrorMetric:
01168     default:
01169     {
01170       status=GetPeakAbsoluteError(image,reconstruct_image,AllChannels,
01171         channel_distortion,exception);
01172       break;
01173     }
01174     case PeakSignalToNoiseRatioMetric:
01175     {
01176       status=GetPeakSignalToNoiseRatio(image,reconstruct_image,AllChannels,
01177         channel_distortion,exception);
01178       break;
01179     }
01180     case RootMeanSquaredErrorMetric:
01181     {
01182       status=GetRootMeanSquaredError(image,reconstruct_image,AllChannels,
01183         channel_distortion,exception);
01184       break;
01185     }
01186   }
01187   return(channel_distortion);
01188 }
01189 
01190 /*
01191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01192 %                                                                             %
01193 %                                                                             %
01194 %                                                                             %
01195 %  I s I m a g e s E q u a l                                                  %
01196 %                                                                             %
01197 %                                                                             %
01198 %                                                                             %
01199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01200 %
01201 %  IsImagesEqual() measures the difference between colors at each pixel
01202 %  location of two images.  A value other than 0 means the colors match
01203 %  exactly.  Otherwise an error measure is computed by summing over all
01204 %  pixels in an image the distance squared in RGB space between each image
01205 %  pixel and its corresponding pixel in the reconstruct image.  The error
01206 %  measure is assigned to these image members:
01207 %
01208 %    o mean_error_per_pixel:  The mean error for any single pixel in
01209 %      the image.
01210 %
01211 %    o normalized_mean_error:  The normalized mean quantization error for
01212 %      any single pixel in the image.  This distance measure is normalized to
01213 %      a range between 0 and 1.  It is independent of the range of red, green,
01214 %      and blue values in the image.
01215 %
01216 %    o normalized_maximum_error:  The normalized maximum quantization
01217 %      error for any single pixel in the image.  This distance measure is
01218 %      normalized to a range between 0 and 1.  It is independent of the range
01219 %      of red, green, and blue values in your image.
01220 %
01221 %  A small normalized mean square error, accessed as
01222 %  image->normalized_mean_error, suggests the images are very similar in
01223 %  spatial layout and color.
01224 %
01225 %  The format of the IsImagesEqual method is:
01226 %
01227 %      MagickBooleanType IsImagesEqual(Image *image,
01228 %        const Image *reconstruct_image)
01229 %
01230 %  A description of each parameter follows.
01231 %
01232 %    o image: the image.
01233 %
01234 %    o reconstruct_image: the reconstruct image.
01235 %
01236 */
01237 MagickExport MagickBooleanType IsImagesEqual(Image *image,
01238   const Image *reconstruct_image)
01239 {
01240   ExceptionInfo
01241     *exception;
01242 
01243   long
01244     y;
01245 
01246   MagickBooleanType
01247     status;
01248 
01249   MagickRealType
01250     area,
01251     maximum_error,
01252     mean_error,
01253     mean_error_per_pixel;
01254 
01255   CacheView
01256     *image_view,
01257     *reconstruct_view;
01258 
01259   assert(image != (Image *) NULL);
01260   assert(image->signature == MagickSignature);
01261   assert(reconstruct_image != (const Image *) NULL);
01262   assert(reconstruct_image->signature == MagickSignature);
01263   if ((reconstruct_image->columns != image->columns) ||
01264       (reconstruct_image->rows != image->rows))
01265     ThrowBinaryException(ImageError,"ImageSizeDiffers",image->filename);
01266   area=0.0;
01267   maximum_error=0.0;
01268   mean_error_per_pixel=0.0;
01269   mean_error=0.0;
01270   exception=(&image->exception);
01271   image_view=AcquireCacheView(image);
01272   reconstruct_view=AcquireCacheView(reconstruct_image);
01273   for (y=0; y < (long) image->rows; y++)
01274   {
01275     register const IndexPacket
01276       *__restrict indexes,
01277       *__restrict reconstruct_indexes;
01278 
01279     register const PixelPacket
01280       *__restrict p,
01281       *__restrict q;
01282 
01283     register long
01284       x;
01285 
01286     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
01287     q=GetCacheViewVirtualPixels(reconstruct_view,0,y,reconstruct_image->columns,
01288       1,exception);
01289     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
01290       break;
01291     indexes=GetCacheViewVirtualIndexQueue(image_view);
01292     reconstruct_indexes=GetCacheViewVirtualIndexQueue(reconstruct_view);
01293     for (x=0; x < (long) image->columns; x++)
01294     {
01295       MagickRealType
01296         distance;
01297 
01298       distance=fabs(p->red-(double) q->red);
01299       mean_error_per_pixel+=distance;
01300       mean_error+=distance*distance;
01301       if (distance > maximum_error)
01302         maximum_error=distance;
01303       area++;
01304       distance=fabs(p->green-(double) q->green);
01305       mean_error_per_pixel+=distance;
01306       mean_error+=distance*distance;
01307       if (distance > maximum_error)
01308         maximum_error=distance;
01309       area++;
01310       distance=fabs(p->blue-(double) q->blue);
01311       mean_error_per_pixel+=distance;
01312       mean_error+=distance*distance;
01313       if (distance > maximum_error)
01314         maximum_error=distance;
01315       area++;
01316       if (image->matte != MagickFalse)
01317         {
01318           distance=fabs(p->opacity-(double) q->opacity);
01319           mean_error_per_pixel+=distance;
01320           mean_error+=distance*distance;
01321           if (distance > maximum_error)
01322             maximum_error=distance;
01323           area++;
01324         }
01325       if ((image->colorspace == CMYKColorspace) &&
01326           (reconstruct_image->colorspace == CMYKColorspace))
01327         {
01328           distance=fabs(indexes[x]-(double) reconstruct_indexes[x]);
01329           mean_error_per_pixel+=distance;
01330           mean_error+=distance*distance;
01331           if (distance > maximum_error)
01332             maximum_error=distance;
01333           area++;
01334         }
01335       p++;
01336       q++;
01337     }
01338   }
01339   reconstruct_view=DestroyCacheView(reconstruct_view);
01340   image_view=DestroyCacheView(image_view);
01341   image->error.mean_error_per_pixel=(double) (mean_error_per_pixel/area);
01342   image->error.normalized_mean_error=(double) (QuantumScale*QuantumScale*
01343     mean_error/area);
01344   image->error.normalized_maximum_error=(double) (QuantumScale*maximum_error);
01345   status=image->error.mean_error_per_pixel == 0.0 ? MagickTrue : MagickFalse;
01346   return(status);
01347 }
01348 
01349 /*
01350 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01351 %                                                                             %
01352 %                                                                             %
01353 %                                                                             %
01354 %   S i m i l a r i t y I m a g e                                             %
01355 %                                                                             %
01356 %                                                                             %
01357 %                                                                             %
01358 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01359 %
01360 %  SimilarityImage() compares the reference image of the image and returns the
01361 %  best match offset.  In addition, it returns a similarity image such that an
01362 %  exact match location is completely white and if none of the pixels match,
01363 %  black, otherwise some gray level in-between.
01364 %
01365 %  The format of the SimilarityImageImage method is:
01366 %
01367 %      Image *SimilarityImage(const Image *image,const Image *reference,
01368 %        RectangleInfo *offset,double *similarity,ExceptionInfo *exception)
01369 %
01370 %  A description of each parameter follows:
01371 %
01372 %    o image: the image.
01373 %
01374 %    o reference: find an area of the image that closely resembles this image.
01375 %
01376 %    o the best match offset of the reference image within the image.
01377 %
01378 %    o similarity: the computed similarity between the images.
01379 %
01380 %    o exception: return any errors or warnings in this structure.
01381 %
01382 */
01383 
01384 static double GetSimilarityMetric(const Image *image,const Image *reference,
01385   const long x_offset,const long y_offset,ExceptionInfo *exception)
01386 {
01387   double
01388     similarity;
01389 
01390   long
01391     y;
01392 
01393   CacheView
01394     *image_view,
01395     *reference_view;
01396 
01397   MagickBooleanType
01398     status;
01399 
01400   /*
01401     Compute the similarity in pixels between two images.
01402   */
01403   status=MagickTrue;
01404   similarity=0.0;
01405   image_view=AcquireCacheView(image);
01406   reference_view=AcquireCacheView(reference);
01407 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01408   #pragma omp parallel for schedule(dynamic,4) shared(status)
01409 #endif
01410   for (y=0; y < (long) reference->rows; y++)
01411   {
01412     register const IndexPacket
01413       *__restrict indexes,
01414       *__restrict reference_indexes;
01415 
01416     register const PixelPacket
01417       *__restrict p,
01418       *__restrict q;
01419 
01420     register long
01421       x;
01422 
01423     if (status == MagickFalse)
01424       continue;
01425     p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset+y,
01426       reference->columns,1,exception);
01427     q=GetCacheViewVirtualPixels(reference_view,0,y,reference->columns,1,
01428       exception);
01429     if ((p == (const PixelPacket *) NULL) || (q == (const PixelPacket *) NULL))
01430       {
01431         status=MagickFalse;
01432         continue;
01433       }
01434     indexes=GetCacheViewVirtualIndexQueue(image_view);
01435     reference_indexes=GetCacheViewVirtualIndexQueue(reference_view);
01436     for (x=0; x < (long) reference->columns; x++)
01437     {
01438       double
01439         thread_similarity;
01440 
01441       MagickRealType
01442         distance;
01443 
01444       thread_similarity=0.0;
01445       distance=QuantumScale*(p->red-(MagickRealType) q->red);
01446       thread_similarity+=distance*distance;
01447       distance=QuantumScale*(p->green-(MagickRealType) q->green);
01448       thread_similarity+=distance*distance;
01449       distance=QuantumScale*(p->blue-(MagickRealType) q->blue);
01450       thread_similarity+=distance*distance;
01451       if ((image->matte != MagickFalse) && (reference->matte != MagickFalse))
01452         {
01453           distance=QuantumScale*(p->opacity-(MagickRealType) q->opacity);
01454           thread_similarity+=distance*distance;
01455         }
01456       if ((image->colorspace == CMYKColorspace) &&
01457           (reference->colorspace == CMYKColorspace))
01458         {
01459           distance=QuantumScale*(indexes[x]-(MagickRealType)
01460             reference_indexes[x]);
01461           thread_similarity+=distance*distance;
01462         }
01463 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01464   #pragma omp critical (MagickCore_GetSimilarityMetric)
01465 #endif
01466       similarity+=thread_similarity;
01467       p++;
01468       q++;
01469     }
01470   }
01471   reference_view=DestroyCacheView(reference_view);
01472   image_view=DestroyCacheView(image_view);
01473   if (status == MagickFalse)
01474     return(0.0);
01475   similarity/=((double) reference->columns*reference->rows);
01476   similarity/=(double) GetNumberChannels(reference,AllChannels);
01477   return(sqrt(similarity));
01478 }
01479 
01480 MagickExport Image *SimilarityImage(Image *image,const Image *reference,
01481   RectangleInfo *offset,double *similarity_metric,ExceptionInfo *exception)
01482 {
01483 #define SimilarityImageTag  "Similarity/Image"
01484 
01485   long
01486     progress,
01487     y;
01488 
01489   Image
01490     *similarity_image;
01491 
01492   MagickBooleanType
01493     status;
01494 
01495   CacheView
01496     *similarity_view;
01497 
01498   assert(image != (const Image *) NULL);
01499   assert(image->signature == MagickSignature);
01500   if (image->debug != MagickFalse)
01501     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01502   assert(exception != (ExceptionInfo *) NULL);
01503   assert(exception->signature == MagickSignature);
01504   assert(offset != (RectangleInfo *) NULL);
01505   SetGeometry(reference,offset);
01506   *similarity_metric=1.0;
01507   if ((reference->columns > image->columns) ||
01508       (reference->rows > image->rows))
01509     ThrowImageException(ImageError,"ImageSizeDiffers");
01510   similarity_image=CloneImage(image,image->columns-reference->columns+1,
01511     image->rows-reference->rows+1,MagickTrue,exception);
01512   if (similarity_image == (Image *) NULL)
01513     return((Image *) NULL);
01514   if (SetImageStorageClass(similarity_image,DirectClass) == MagickFalse)
01515     {
01516       InheritException(exception,&similarity_image->exception);
01517       similarity_image=DestroyImage(similarity_image);
01518       return((Image *) NULL);
01519     }
01520   /*
01521     Measure similarity of reference image against image.
01522   */
01523   status=MagickTrue;
01524   progress=0;
01525   similarity_view=AcquireCacheView(similarity_image);
01526 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01527   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
01528 #endif
01529   for (y=0; y < (long) (image->rows-reference->rows+1); y++)
01530   {
01531     double
01532       similarity;
01533 
01534     register long
01535       x;
01536 
01537     register PixelPacket
01538       *__restrict q;
01539 
01540     if (status == MagickFalse)
01541       continue;
01542     q=GetCacheViewAuthenticPixels(similarity_view,0,y,
01543       similarity_image->columns,1,exception);
01544     if (q == (const PixelPacket *) NULL)
01545       {
01546         status=MagickFalse;
01547         continue;
01548       }
01549     for (x=0; x < (long) (image->columns-reference->columns+1); x++)
01550     {
01551       similarity=GetSimilarityMetric(image,reference,x,y,exception);
01552 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01553   #pragma omp critical (MagickCore_SimilarityImage)
01554 #endif
01555       if (similarity < *similarity_metric)
01556         {
01557           *similarity_metric=similarity;
01558           offset->x=x;
01559           offset->y=y;
01560         }
01561       q->red=RoundToQuantum(QuantumRange-QuantumRange*similarity);
01562       q->green=q->red;
01563       q->blue=q->red;
01564       q++;
01565     }
01566     if (SyncCacheViewAuthenticPixels(similarity_view,exception) == MagickFalse)
01567       status=MagickFalse;
01568     if (image->progress_monitor != (MagickProgressMonitor) NULL)
01569       {
01570         MagickBooleanType
01571           proceed;
01572 
01573 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01574   #pragma omp critical (MagickCore_SimilarityImage)
01575 #endif
01576         proceed=SetImageProgress(image,SimilarityImageTag,progress++,
01577           image->rows);
01578         if (proceed == MagickFalse)
01579           status=MagickFalse;
01580       }
01581   }
01582   similarity_view=DestroyCacheView(similarity_view);
01583   return(similarity_image);
01584 }

Generated on 19 Nov 2009 for MagickCore by  doxygen 1.6.1