granicus.if.org Git - imagemagick/blob - MagickCore/accelerate-private.h

   1 /*
   2   Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization
   3   dedicated to making software imaging solutions freely available.
   4
   5   You may not use this file except in compliance with the License.
   6   obtain a copy of the License at
   7
   8     http://www.imagemagick.org/script/license.php
   9
  10   Unless required by applicable law or agreed to in writing, software
  11   distributed under the License is distributed on an "AS IS" BASIS,
  12   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13   See the License for the specific language governing permissions and
  14   limitations under the License.
  15
  16   MagickCore private methods for accelerated functions.
  17 */
  18
  19 #ifndef _MAGICKCORE_ACCELERATE_PRIVATE_H
  20 #define _MAGICKCORE_ACCELERATE_PRIVATE_H
  21
  22 #if defined(__cplusplus) || defined(c_plusplus)
  23 extern "C" {
  24 #endif
  25
  26
  27 #if defined(MAGICKCORE_OPENCL_SUPPORT)
  28
  29 #define OPENCL_DEFINE(VAR,...)  "\n #""define " #VAR " " #__VA_ARGS__ " \n"
  30 #define OPENCL_ELIF(...)        "\n #""elif " #__VA_ARGS__ " \n"
  31 #define OPENCL_ELSE()           "\n #""else " " \n"
  32 #define OPENCL_ENDIF()          "\n #""endif " " \n"
  33 #define OPENCL_IF(...)          "\n #""if " #__VA_ARGS__ " \n"
  34 #define STRINGIFY(...) #__VA_ARGS__ "\n"
  35
  36 typedef struct _FloatPixelPacket
  37 {
  38 #ifdef MAGICK_PIXEL_RGBA
  39   MagickRealType
  40     red,
  41     green,
  42     blue,
  43     opacity;
  44 #endif
  45 #ifdef MAGICK_PIXEL_BGRA
  46   MagickRealType
  47     blue,
  48     green,
  49     red,
  50     opacity;
  51 #endif
  52 } FloatPixelPacket;
  53
  54 const char* accelerateKernels =
  55   STRINGIFY(
  56      typedef enum
  57      {
  58        UndefinedChannel,
  59        RedChannel = 0x0001,
  60        GrayChannel = 0x0001,
  61        CyanChannel = 0x0001,
  62        GreenChannel = 0x0002,
  63        MagentaChannel = 0x0002,
  64        BlueChannel = 0x0004,
  65        YellowChannel = 0x0004,
  66        AlphaChannel = 0x0008,
  67        OpacityChannel = 0x0008,
  68        MatteChannel = 0x0008,     /* deprecated */
  69        BlackChannel = 0x0020,
  70        IndexChannel = 0x0020,
  71        CompositeChannels = 0x002F,
  72        AllChannels = 0x7ffffff,
  73        /*
  74        Special purpose channel types.
  75        */
  76        TrueAlphaChannel = 0x0040, /* extract actual alpha channel from opacity */
  77        RGBChannels = 0x0080,      /* set alpha from  grayscale mask in RGB */
  78        GrayChannels = 0x0080,
  79        SyncChannels = 0x0100,     /* channels should be modified equally */
  80        DefaultChannels = ((AllChannels | SyncChannels) &~ OpacityChannel)
  81      } ChannelType;
  82   )
  83
  84   OPENCL_IF((MAGICKCORE_QUANTUM_DEPTH == 8))
  85
  86   STRINGIFY(
  87     inline CLQuantum ScaleCharToQuantum(const unsigned char value)
  88     {
  89       return((CLQuantum) value);
  90     }
  91   )
  92
  93   OPENCL_ELIF((MAGICKCORE_QUANTUM_DEPTH == 16))
  94
  95   STRINGIFY(
  96     inline CLQuantum ScaleCharToQuantum(const unsigned char value)
  97     {
  98       return((CLQuantum) (257.0f*value));
  99     }
 100   )
 101
 102   OPENCL_ELIF((MAGICKCORE_QUANTUM_DEPTH == 32))
 103
 104   STRINGIFY(
 105     inline CLQuantum ScaleCharToQuantum(const unsigned char value)
 106     {
 107       return((Quantum) (16843009.0*value));
 108     }
 109   )
 110
 111   OPENCL_ENDIF()
 112
 113
 114   STRINGIFY(
 115     inline int ClampToCanvas(const int offset,const int range)
 116       {
 117         return clamp(offset, (int)0, range-1);
 118       }
 119   )
 120
 121   STRINGIFY(
 122     inline int ClampToCanvasWithHalo(const int offset,const int range, const int edge, const int section)
 123       {
 124         return clamp(offset, section?(int)(0-edge):(int)0, section?(range-1):(range-1+edge));
 125       }
 126   )
 127
 128   STRINGIFY(
 129     inline CLQuantum ClampToQuantum(const float value)
 130       {
 131         return (CLQuantum) (clamp(value, 0.0f, (float) QuantumRange) + 0.5f);
 132       }
 133   )
 134
 135   STRINGIFY(
 136     inline uint ScaleQuantumToMap(CLQuantum value)
 137       {
 138         if (value >= (CLQuantum) MaxMap)
 139           return ((uint)MaxMap);
 140         else
 141           return ((uint)value);
 142       }
 143   )
 144
 145   STRINGIFY(
 146     inline float PerceptibleReciprocal(const float x)
 147     {
 148       float sign = x < (float) 0.0 ? (float) -1.0 : (float) 1.0;
 149       return((sign*x) >= MagickEpsilon ? (float) 1.0/x : sign*((float) 1.0/MagickEpsilon));
 150     }
 151   )
 152
 153   OPENCL_DEFINE(GetPixelAlpha(pixel),(QuantumRange-(pixel).w))
 154
 155   STRINGIFY(
 156   typedef enum
 157   {
 158     UndefinedPixelIntensityMethod = 0,
 159     AveragePixelIntensityMethod,
 160     BrightnessPixelIntensityMethod,
 161     LightnessPixelIntensityMethod,
 162     Rec601LumaPixelIntensityMethod,
 163     Rec601LuminancePixelIntensityMethod,
 164     Rec709LumaPixelIntensityMethod,
 165     Rec709LuminancePixelIntensityMethod,
 166     RMSPixelIntensityMethod,
 167     MSPixelIntensityMethod
 168   } PixelIntensityMethod;
 169   )
 170
 171   STRINGIFY(
 172   typedef enum
 173   {
 174     UndefinedColorspace,
 175     RGBColorspace,            /* Linear RGB colorspace */
 176     GRAYColorspace,           /* greyscale (linear) image (faked 1 channel) */
 177     TransparentColorspace,
 178     OHTAColorspace,
 179     LabColorspace,
 180     XYZColorspace,
 181     YCbCrColorspace,
 182     YCCColorspace,
 183     YIQColorspace,
 184     YPbPrColorspace,
 185     YUVColorspace,
 186     CMYKColorspace,           /* negared linear RGB with black separated */
 187     sRGBColorspace,           /* Default: non-lienar sRGB colorspace */
 188     HSBColorspace,
 189     HSLColorspace,
 190     HWBColorspace,
 191     Rec601LumaColorspace,
 192     Rec601YCbCrColorspace,
 193     Rec709LumaColorspace,
 194     Rec709YCbCrColorspace,
 195     LogColorspace,
 196     CMYColorspace,            /* negated linear RGB colorspace */
 197     LuvColorspace,
 198     HCLColorspace,
 199     LCHColorspace,            /* alias for LCHuv */
 200     LMSColorspace,
 201     LCHabColorspace,          /* Cylindrical (Polar) Lab */
 202     LCHuvColorspace,          /* Cylindrical (Polar) Luv */
 203     scRGBColorspace,
 204     HSIColorspace,
 205     HSVColorspace,            /* alias for HSB */
 206     HCLpColorspace,
 207     YDbDrColorspace
 208   } ColorspaceType;
 209   )
 210
 211   STRINGIFY(
 212   inline float RoundToUnity(const float value)
 213    {
 214      return clamp(value,0.0f,1.0f);
 215    }
 216   )
 217
 218   STRINGIFY(
 219
 220   inline CLQuantum getBlue(CLPixelType p)                   { return p.x; }
 221   inline void setBlue(CLPixelType* p, CLQuantum value)      { (*p).x = value; }
 222   inline float getBlueF4(float4 p)                          { return p.x; }
 223   inline void setBlueF4(float4* p, float value)             { (*p).x = value; }
 224
 225   inline CLQuantum getGreen(CLPixelType p)                  { return p.y; }
 226   inline void setGreen(CLPixelType* p, CLQuantum value)     { (*p).y = value; }
 227   inline float getGreenF4(float4 p)                         { return p.y; }
 228   inline void setGreenF4(float4* p, float value)            { (*p).y = value; }
 229
 230   inline CLQuantum getRed(CLPixelType p)                    { return p.z; }
 231   inline void setRed(CLPixelType* p, CLQuantum value)       { (*p).z = value; }
 232   inline float getRedF4(float4 p)                           { return p.z; }
 233   inline void setRedF4(float4* p, float value)              { (*p).z = value; }
 234
 235   inline CLQuantum getOpacity(CLPixelType p)                { return p.w; }
 236   inline void setOpacity(CLPixelType* p, CLQuantum value)   { (*p).w = value; }
 237   inline float getOpacityF4(float4 p)                       { return p.w; }
 238   inline void setOpacityF4(float4* p, float value)          { (*p).w = value; }
 239
 240   inline void setGray(CLPixelType* p, CLQuantum value)      { (*p).z = value; (*p).y = value; (*p).x = value; }
 241
 242   inline float GetPixelIntensity(const int method, const int colorspace, CLPixelType p)
 243   {
 244     float red = getRed(p);
 245     float green = getGreen(p);
 246     float blue = getBlue(p);
 247
 248     float intensity;
 249
 250     if (colorspace == GRAYColorspace)
 251       return red;
 252
 253     switch (method)
 254     {
 255       case AveragePixelIntensityMethod:
 256         {
 257           intensity=(red+green+blue)/3.0;
 258           break;
 259         }
 260       case BrightnessPixelIntensityMethod:
 261         {
 262           intensity=max(max(red,green),blue);
 263           break;
 264         }
 265       case LightnessPixelIntensityMethod:
 266         {
 267           intensity=(min(min(red,green),blue)+
 268               max(max(red,green),blue))/2.0;
 269           break;
 270         }
 271       case MSPixelIntensityMethod:
 272         {
 273           intensity=(float) (((float) red*red+green*green+blue*blue)/
 274               (3.0*QuantumRange));
 275           break;
 276         }
 277       case Rec601LumaPixelIntensityMethod:
 278         {
 279           /*
 280           if (image->colorspace == RGBColorspace)
 281           {
 282             red=EncodePixelGamma(red);
 283             green=EncodePixelGamma(green);
 284             blue=EncodePixelGamma(blue);
 285           }
 286           */
 287           intensity=0.298839*red+0.586811*green+0.114350*blue;
 288           break;
 289         }
 290       case Rec601LuminancePixelIntensityMethod:
 291         {
 292           /*
 293           if (image->colorspace == sRGBColorspace)
 294           {
 295             red=DecodePixelGamma(red);
 296             green=DecodePixelGamma(green);
 297             blue=DecodePixelGamma(blue);
 298           }
 299           */
 300           intensity=0.298839*red+0.586811*green+0.114350*blue;
 301           break;
 302         }
 303       case Rec709LumaPixelIntensityMethod:
 304       default:
 305         {
 306           /*
 307           if (image->colorspace == RGBColorspace)
 308           {
 309             red=EncodePixelGamma(red);
 310             green=EncodePixelGamma(green);
 311             blue=EncodePixelGamma(blue);
 312           }
 313           */
 314           intensity=0.212656*red+0.715158*green+0.072186*blue;
 315           break;
 316         }
 317       case Rec709LuminancePixelIntensityMethod:
 318         {
 319           /*
 320           if (image->colorspace == sRGBColorspace)
 321           {
 322             red=DecodePixelGamma(red);
 323             green=DecodePixelGamma(green);
 324             blue=DecodePixelGamma(blue);
 325           }
 326           */
 327           intensity=0.212656*red+0.715158*green+0.072186*blue;
 328           break;
 329         }
 330       case RMSPixelIntensityMethod:
 331         {
 332           intensity=(float) (sqrt((float) red*red+green*green+blue*blue)/
 333               sqrt(3.0));
 334           break;
 335         }
 336     }
 337
 338     return intensity;
 339
 340   }
 341   )
 342
 343   STRINGIFY(
 344     __kernel
 345     void ConvolveOptimized(const __global CLPixelType *input, __global CLPixelType *output,
 346     const unsigned int imageWidth, const unsigned int imageHeight,
 347     __constant float *filter, const unsigned int filterWidth, const unsigned int filterHeight,
 348     const uint matte, const ChannelType channel, __local CLPixelType *pixelLocalCache, __local float* filterCache) {
 349
 350       int2 blockID;
 351       blockID.x = get_group_id(0);
 352       blockID.y = get_group_id(1);
 353
 354       // image area processed by this workgroup
 355       int2 imageAreaOrg;
 356       imageAreaOrg.x = blockID.x * get_local_size(0);
 357       imageAreaOrg.y = blockID.y * get_local_size(1);
 358
 359       int2 midFilterDimen;
 360       midFilterDimen.x = (filterWidth-1)/2;
 361       midFilterDimen.y = (filterHeight-1)/2;
 362
 363       int2 cachedAreaOrg = imageAreaOrg - midFilterDimen;
 364
 365       // dimension of the local cache
 366       int2 cachedAreaDimen;
 367       cachedAreaDimen.x = get_local_size(0) + filterWidth - 1;
 368       cachedAreaDimen.y = get_local_size(1) + filterHeight - 1;
 369
 370       // cache the pixels accessed by this workgroup in local memory
 371       int localID = get_local_id(1)*get_local_size(0)+get_local_id(0);
 372       int cachedAreaNumPixels = cachedAreaDimen.x * cachedAreaDimen.y;
 373       int groupSize = get_local_size(0) * get_local_size(1);
 374       for (int i = localID; i < cachedAreaNumPixels; i+=groupSize) {
 375
 376         int2 cachedAreaIndex;
 377         cachedAreaIndex.x = i % cachedAreaDimen.x;
 378         cachedAreaIndex.y = i / cachedAreaDimen.x;
 379
 380         int2 imagePixelIndex;
 381         imagePixelIndex = cachedAreaOrg + cachedAreaIndex;
 382
 383         // only support EdgeVirtualPixelMethod through ClampToCanvas
 384         // TODO: implement other virtual pixel method
 385         imagePixelIndex.x = ClampToCanvas(imagePixelIndex.x, imageWidth);
 386         imagePixelIndex.y = ClampToCanvas(imagePixelIndex.y, imageHeight);
 387
 388         pixelLocalCache[i] = input[imagePixelIndex.y * imageWidth + imagePixelIndex.x];
 389       }
 390
 391       // cache the filter
 392       for (int i = localID; i < filterHeight*filterWidth; i+=groupSize) {
 393         filterCache[i] = filter[i];
 394       }
 395       barrier(CLK_LOCAL_MEM_FENCE);
 396
 397
 398       int2 imageIndex;
 399       imageIndex.x = imageAreaOrg.x + get_local_id(0);
 400       imageIndex.y = imageAreaOrg.y + get_local_id(1);
 401
 402       // if out-of-range, stops here and quit
 403       if (imageIndex.x >= imageWidth
 404         || imageIndex.y >= imageHeight) {
 405           return;
 406       }
 407
 408       int filterIndex = 0;
 409       float4 sum = (float4)0.0f;
 410       float gamma = 0.0f;
 411       if (((channel & OpacityChannel) == 0) || (matte == 0)) {
 412         int cacheIndexY = get_local_id(1);
 413         for (int j = 0; j < filterHeight; j++) {
 414           int cacheIndexX = get_local_id(0);
 415           for (int i = 0; i < filterWidth; i++) {
 416             CLPixelType p = pixelLocalCache[cacheIndexY*cachedAreaDimen.x + cacheIndexX];
 417             float f = filterCache[filterIndex];
 418
 419             sum.x += f * p.x;
 420             sum.y += f * p.y;
 421             sum.z += f * p.z;
 422             sum.w += f * p.w;
 423
 424             gamma += f;
 425             filterIndex++;
 426             cacheIndexX++;
 427           }
 428           cacheIndexY++;
 429         }
 430       }
 431       else {
 432         int cacheIndexY = get_local_id(1);
 433         for (int j = 0; j < filterHeight; j++) {
 434           int cacheIndexX = get_local_id(0);
 435           for (int i = 0; i < filterWidth; i++) {
 436
 437             CLPixelType p = pixelLocalCache[cacheIndexY*cachedAreaDimen.x + cacheIndexX];
 438             float alpha = QuantumScale*(QuantumRange-p.w);
 439             float f = filterCache[filterIndex];
 440             float g = alpha * f;
 441
 442             sum.x += g*p.x;
 443             sum.y += g*p.y;
 444             sum.z += g*p.z;
 445             sum.w += f*p.w;
 446
 447             gamma += g;
 448             filterIndex++;
 449             cacheIndexX++;
 450           }
 451           cacheIndexY++;
 452         }
 453         gamma = PerceptibleReciprocal(gamma);
 454         sum.xyz = gamma*sum.xyz;
 455       }
 456       CLPixelType outputPixel;
 457       outputPixel.x = ClampToQuantum(sum.x);
 458       outputPixel.y = ClampToQuantum(sum.y);
 459       outputPixel.z = ClampToQuantum(sum.z);
 460       outputPixel.w = ((channel & OpacityChannel)!=0)?ClampToQuantum(sum.w):input[imageIndex.y * imageWidth + imageIndex.x].w;
 461
 462       output[imageIndex.y * imageWidth + imageIndex.x] = outputPixel;
 463     }
 464   )
 465
 466   STRINGIFY(
 467     __kernel
 468     void Convolve(const __global CLPixelType *input, __global CLPixelType *output,
 469                   const uint imageWidth, const uint imageHeight,
 470                   __constant float *filter, const unsigned int filterWidth, const unsigned int filterHeight,
 471                   const uint matte, const ChannelType channel) {
 472
 473       int2 imageIndex;
 474       imageIndex.x = get_global_id(0);
 475       imageIndex.y = get_global_id(1);
 476
 477       /*
 478       unsigned int imageWidth = get_global_size(0);
 479       unsigned int imageHeight = get_global_size(1);
 480       */
 481       if (imageIndex.x >= imageWidth
 482           || imageIndex.y >= imageHeight)
 483           return;
 484
 485       int2 midFilterDimen;
 486       midFilterDimen.x = (filterWidth-1)/2;
 487       midFilterDimen.y = (filterHeight-1)/2;
 488
 489       int filterIndex = 0;
 490       float4 sum = (float4)0.0f;
 491       float gamma = 0.0f;
 492       if (((channel & OpacityChannel) == 0) || (matte == 0)) {
 493         for (int j = 0; j < filterHeight; j++) {
 494           int2 inputPixelIndex;
 495           inputPixelIndex.y = imageIndex.y - midFilterDimen.y + j;
 496           inputPixelIndex.y = ClampToCanvas(inputPixelIndex.y, imageHeight);
 497           for (int i = 0; i < filterWidth; i++) {
 498             inputPixelIndex.x = imageIndex.x - midFilterDimen.x + i;
 499             inputPixelIndex.x = ClampToCanvas(inputPixelIndex.x, imageWidth);
 500
 501             CLPixelType p = input[inputPixelIndex.y * imageWidth + inputPixelIndex.x];
 502             float f = filter[filterIndex];
 503
 504             sum.x += f * p.x;
 505             sum.y += f * p.y;
 506             sum.z += f * p.z;
 507             sum.w += f * p.w;
 508
 509             gamma += f;
 510
 511             filterIndex++;
 512           }
 513         }
 514       }
 515       else {
 516
 517         for (int j = 0; j < filterHeight; j++) {
 518           int2 inputPixelIndex;
 519           inputPixelIndex.y = imageIndex.y - midFilterDimen.y + j;
 520           inputPixelIndex.y = ClampToCanvas(inputPixelIndex.y, imageHeight);
 521           for (int i = 0; i < filterWidth; i++) {
 522             inputPixelIndex.x = imageIndex.x - midFilterDimen.x + i;
 523             inputPixelIndex.x = ClampToCanvas(inputPixelIndex.x, imageWidth);
 524
 525             CLPixelType p = input[inputPixelIndex.y * imageWidth + inputPixelIndex.x];
 526             float alpha = QuantumScale*(QuantumRange-p.w);
 527             float f = filter[filterIndex];
 528             float g = alpha * f;
 529
 530             sum.x += g*p.x;
 531             sum.y += g*p.y;
 532             sum.z += g*p.z;
 533             sum.w += f*p.w;
 534
 535             gamma += g;
 536
 537
 538             filterIndex++;
 539           }
 540         }
 541         gamma = PerceptibleReciprocal(gamma);
 542         sum.xyz = gamma*sum.xyz;
 543       }
 544
 545       CLPixelType outputPixel;
 546       outputPixel.x = ClampToQuantum(sum.x);
 547       outputPixel.y = ClampToQuantum(sum.y);
 548       outputPixel.z = ClampToQuantum(sum.z);
 549       outputPixel.w = ((channel & OpacityChannel)!=0)?ClampToQuantum(sum.w):input[imageIndex.y * imageWidth + imageIndex.x].w;
 550
 551       output[imageIndex.y * imageWidth + imageIndex.x] = outputPixel;
 552     }
 553   )
 554
 555   STRINGIFY(
 556      typedef enum
 557      {
 558        UndefinedFunction,
 559        PolynomialFunction,
 560        SinusoidFunction,
 561        ArcsinFunction,
 562        ArctanFunction
 563      } MagickFunction;
 564   )
 565
 566   STRINGIFY(
 567
 568     /*
 569     apply FunctionImageChannel(braightness-contrast)
 570     */
 571     CLPixelType ApplyFunction(CLPixelType pixel,const MagickFunction function,
 572         const unsigned int number_parameters,
 573         __constant float *parameters)
 574       {
 575         float4 result = (float4) 0.0f;
 576         switch (function)
 577         {
 578         case PolynomialFunction:
 579           {
 580             for (unsigned int i=0; i < number_parameters; i++)
 581               result = result*(float4)QuantumScale*convert_float4(pixel) + parameters[i];
 582             result *= (float4)QuantumRange;
 583             break;
 584           }
 585         case SinusoidFunction:
 586           {
 587             float  freq,phase,ampl,bias;
 588             freq  = ( number_parameters >= 1 ) ? parameters[0] : 1.0f;
 589             phase = ( number_parameters >= 2 ) ? parameters[1] : 0.0f;
 590             ampl  = ( number_parameters >= 3 ) ? parameters[2] : 0.5f;
 591             bias  = ( number_parameters >= 4 ) ? parameters[3] : 0.5f;
 592             result.x = QuantumRange*(ampl*sin(2.0f*MagickPI*
 593               (freq*QuantumScale*(float)pixel.x + phase/360.0f)) + bias);
 594             result.y = QuantumRange*(ampl*sin(2.0f*MagickPI*
 595               (freq*QuantumScale*(float)pixel.y + phase/360.0f)) + bias);
 596             result.z = QuantumRange*(ampl*sin(2.0f*MagickPI*
 597               (freq*QuantumScale*(float)pixel.z + phase/360.0f)) + bias);
 598             result.w = QuantumRange*(ampl*sin(2.0f*MagickPI*
 599               (freq*QuantumScale*(float)pixel.w + phase/360.0f)) + bias);
 600             break;
 601           }
 602         case ArcsinFunction:
 603           {
 604             float  width,range,center,bias;
 605             width  = ( number_parameters >= 1 ) ? parameters[0] : 1.0f;
 606             center = ( number_parameters >= 2 ) ? parameters[1] : 0.5f;
 607             range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0f;
 608             bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5f;
 609
 610             result.x = 2.0f/width*(QuantumScale*(float)pixel.x - center);
 611             result.x = range/MagickPI*asin(result.x)+bias;
 612             result.x = ( result.x <= -1.0f ) ? bias - range/2.0f : result.x;
 613             result.x = ( result.x >= 1.0f ) ? bias + range/2.0f : result.x;
 614
 615             result.y = 2.0f/width*(QuantumScale*(float)pixel.y - center);
 616             result.y = range/MagickPI*asin(result.y)+bias;
 617             result.y = ( result.y <= -1.0f ) ? bias - range/2.0f : result.y;
 618             result.y = ( result.y >= 1.0f ) ? bias + range/2.0f : result.y;
 619
 620             result.z = 2.0f/width*(QuantumScale*(float)pixel.z - center);
 621             result.z = range/MagickPI*asin(result.z)+bias;
 622             result.z = ( result.z <= -1.0f ) ? bias - range/2.0f : result.x;
 623             result.z = ( result.z >= 1.0f ) ? bias + range/2.0f : result.x;
 624
 625
 626             result.w = 2.0f/width*(QuantumScale*(float)pixel.w - center);
 627             result.w = range/MagickPI*asin(result.w)+bias;
 628             result.w = ( result.w <= -1.0f ) ? bias - range/2.0f : result.w;
 629             result.w = ( result.w >= 1.0f ) ? bias + range/2.0f : result.w;
 630
 631             result *= (float4)QuantumRange;
 632             break;
 633           }
 634         case ArctanFunction:
 635           {
 636             float slope,range,center,bias;
 637             slope  = ( number_parameters >= 1 ) ? parameters[0] : 1.0f;
 638             center = ( number_parameters >= 2 ) ? parameters[1] : 0.5f;
 639             range  = ( number_parameters >= 3 ) ? parameters[2] : 1.0f;
 640             bias   = ( number_parameters >= 4 ) ? parameters[3] : 0.5f;
 641             result = (float4)MagickPI*(float4)slope*((float4)QuantumScale*convert_float4(pixel)-(float4)center);
 642             result = (float4)QuantumRange*((float4)range/(float4)MagickPI*atan(result) + (float4)bias);
 643             break;
 644           }
 645         case UndefinedFunction:
 646           break;
 647         }
 648         return (CLPixelType) (ClampToQuantum(result.x), ClampToQuantum(result.y),
 649           ClampToQuantum(result.z), ClampToQuantum(result.w));
 650       }
 651     )
 652
 653     STRINGIFY(
 654     /*
 655     Improve brightness / contrast of the image
 656     channel : define which channel is improved
 657     function : the function called to enchance the brightness contrast
 658     number_parameters : numbers of parameters
 659     parameters : the parameter
 660     */
 661     __kernel void FunctionImage(__global CLPixelType *im,
 662                                         const ChannelType channel, const MagickFunction function,
 663                                         const unsigned int number_parameters, __constant float *parameters)
 664       {
 665         const int x = get_global_id(0);
 666         const int y = get_global_id(1);
 667         const int columns = get_global_size(0);
 668         const int c = x + y * columns;
 669         im[c] = ApplyFunction(im[c], function, number_parameters, parameters);
 670       }
 671     )
 672
 673     STRINGIFY(
 674     /*
 675     */
 676     __kernel void Stretch(__global CLPixelType * restrict im,
 677       const ChannelType channel,
 678       __global CLPixelType * restrict stretch_map,
 679       const float4 white, const float4 black)
 680       {
 681         const int x = get_global_id(0);
 682         const int y = get_global_id(1);
 683         const int columns = get_global_size(0);
 684         const int c = x + y * columns;
 685
 686         uint ePos;
 687         CLPixelType oValue, eValue;
 688         CLQuantum red, green, blue, opacity;
 689
 690         //read from global
 691         oValue=im[c];
 692
 693         if ((channel & RedChannel) != 0)
 694         {
 695           if (getRedF4(white) != getRedF4(black))
 696           {
 697             ePos = ScaleQuantumToMap(getRed(oValue));
 698             eValue = stretch_map[ePos];
 699             red = getRed(eValue);
 700           }
 701         }
 702
 703         if ((channel & GreenChannel) != 0)
 704         {
 705           if (getGreenF4(white) != getGreenF4(black))
 706           {
 707             ePos = ScaleQuantumToMap(getGreen(oValue));
 708             eValue = stretch_map[ePos];
 709             green = getGreen(eValue);
 710           }
 711         }
 712
 713         if ((channel & BlueChannel) != 0)
 714         {
 715           if (getBlueF4(white) != getBlueF4(black))
 716           {
 717             ePos = ScaleQuantumToMap(getBlue(oValue));
 718             eValue = stretch_map[ePos];
 719             blue = getBlue(eValue);
 720           }
 721         }
 722
 723         if ((channel & OpacityChannel) != 0)
 724         {
 725           if (getOpacityF4(white) != getOpacityF4(black))
 726           {
 727             ePos = ScaleQuantumToMap(getOpacity(oValue));
 728             eValue = stretch_map[ePos];
 729             opacity = getOpacity(eValue);
 730           }
 731         }
 732
 733         //write back
 734         im[c]=(CLPixelType)(blue, green, red, opacity);
 735
 736       }
 737     )
 738
 739
 740     STRINGIFY(
 741     /*
 742     */
 743     __kernel void Equalize(__global CLPixelType * restrict im,
 744       const ChannelType channel,
 745       __global CLPixelType * restrict equalize_map,
 746       const float4 white, const float4 black)
 747       {
 748         const int x = get_global_id(0);
 749         const int y = get_global_id(1);
 750         const int columns = get_global_size(0);
 751         const int c = x + y * columns;
 752
 753         uint ePos;
 754         CLPixelType oValue, eValue;
 755         CLQuantum red, green, blue, opacity;
 756
 757         //read from global
 758         oValue=im[c];
 759
 760         if ((channel & SyncChannels) != 0)
 761         {
 762           if (getRedF4(white) != getRedF4(black))
 763           {
 764             ePos = ScaleQuantumToMap(getRed(oValue));
 765             eValue = equalize_map[ePos];
 766             red = getRed(eValue);
 767             ePos = ScaleQuantumToMap(getGreen(oValue));
 768             eValue = equalize_map[ePos];
 769             green = getRed(eValue);
 770             ePos = ScaleQuantumToMap(getBlue(oValue));
 771             eValue = equalize_map[ePos];
 772             blue = getRed(eValue);
 773             ePos = ScaleQuantumToMap(getOpacity(oValue));
 774             eValue = equalize_map[ePos];
 775             opacity = getRed(eValue);
 776
 777             //write back
 778             im[c]=(CLPixelType)(blue, green, red, opacity);
 779           }
 780
 781         }
 782
 783         // for equalizing, we always need all channels?
 784         // otherwise something more
 785
 786      }
 787     )
 788
 789     STRINGIFY(
 790     /*
 791     */
 792     __kernel void Histogram(__global CLPixelType * restrict im,
 793       const ChannelType channel,
 794       const int method,
 795       const int colorspace,
 796       __global uint4 * restrict histogram)
 797       {
 798         const int x = get_global_id(0);
 799         const int y = get_global_id(1);
 800         const int columns = get_global_size(0);
 801         const int c = x + y * columns;
 802         if ((channel & SyncChannels) != 0)
 803         {
 804           float intensity = GetPixelIntensity(method, colorspace,im[c]);
 805           uint pos = ScaleQuantumToMap(ClampToQuantum(intensity));
 806           atomic_inc((__global uint *)(&(histogram[pos]))+2); //red position
 807         }
 808         else
 809         {
 810           // for equalizing, we always need all channels?
 811           // otherwise something more
 812         }
 813       }
 814     )
 815
 816     STRINGIFY(
 817       /*
 818       Reduce image noise and reduce detail levels by row
 819       im: input pixels filtered_in  filtered_im: output pixels
 820       filter : convolve kernel  width: convolve kernel size
 821       channel : define which channel is blured
 822       is_RGBA_BGRA : define the input is RGBA or BGRA
 823       */
 824       __kernel void BlurRow(__global CLPixelType *im, __global float4 *filtered_im,
 825                          const ChannelType channel, __constant float *filter,
 826                          const unsigned int width,
 827                          const unsigned int imageColumns, const unsigned int imageRows,
 828                          __local CLPixelType *temp)
 829       {
 830         const int x = get_global_id(0);
 831         const int y = get_global_id(1);
 832
 833         const int columns = imageColumns;
 834
 835         const unsigned int radius = (width-1)/2;
 836         const int wsize = get_local_size(0);
 837         const unsigned int loadSize = wsize+width;
 838
 839         //load chunk only for now
 840         //event_t e = async_work_group_copy(temp+radius, im+x+y*columns, wsize, 0);
 841         //wait_group_events(1,&e);
 842
 843         //parallel load and clamp
 844         /*
 845         int count = 0;
 846         for (int i=0; i < loadSize; i=i+wsize)
 847         {
 848           int currentX = x + wsize*(count++);
 849
 850           int localId = get_local_id(0);
 851
 852           if ((localId+i) > loadSize)
 853             break;
 854
 855           temp[localId+i] = im[y*columns+ClampToCanvas(currentX-radius, columns)];
 856
 857           if (y==0 && get_group_id(0) == 0)
 858           {
 859             printf("(%d %d) temp %d load %d currentX %d\n", x, y, localId+i, ClampToCanvas(currentX-radius, columns), currentX);
 860           }
 861         }
 862         */
 863
 864         //group coordinate
 865         const int groupX=get_local_size(0)*get_group_id(0);
 866         const int groupY=get_local_size(1)*get_group_id(1);
 867
 868         //parallel load and clamp
 869         for (int i=get_local_id(0); i < loadSize; i=i+get_local_size(0))
 870         {
 871           //int cx = ClampToCanvas(groupX+i, columns);
 872           temp[i] = im[y * columns + ClampToCanvas(i+groupX-radius, columns)];
 873
 874           if (0 && y==0 && get_group_id(1) == 0)
 875           {
 876             printf("(%d %d) temp %d load %d groupX %d\n", x, y, i, ClampToCanvas(groupX+i, columns), groupX);
 877           }
 878         }
 879
 880         // barrier
 881         barrier(CLK_LOCAL_MEM_FENCE);
 882
 883         // only do the work if this is not a patched item
 884         if (get_global_id(0) < columns)
 885         {
 886           // compute
 887           float4 result = (float4) 0;
 888
 889           int i = 0;
 890
 891           \n #ifndef UFACTOR   \n
 892           \n #define UFACTOR 8 \n
 893           \n #endif                  \n
 894
 895           for ( ; i+UFACTOR < width; )
 896           {
 897             \n #pragma unroll UFACTOR\n
 898             for (int j=0; j < UFACTOR; j++, i++)
 899             {
 900               result+=filter[i]*convert_float4(temp[i+get_local_id(0)]);
 901             }
 902           }
 903
 904           for ( ; i < width; i++)
 905           {
 906             result+=filter[i]*convert_float4(temp[i+get_local_id(0)]);
 907           }
 908
 909           result.x = ClampToQuantum(result.x);
 910           result.y = ClampToQuantum(result.y);
 911           result.z = ClampToQuantum(result.z);
 912           result.w = ClampToQuantum(result.w);
 913
 914           // write back to global
 915           filtered_im[y*columns+x] = result;
 916         }
 917       }
 918     )
 919
 920     STRINGIFY(
 921       /*
 922       Reduce image noise and reduce detail levels by row
 923       im: input pixels filtered_in  filtered_im: output pixels
 924       filter : convolve kernel  width: convolve kernel size
 925       channel : define which channel is blured
 926       is_RGBA_BGRA : define the input is RGBA or BGRA
 927       */
 928       __kernel void BlurRowSection(__global CLPixelType *im, __global float4 *filtered_im,
 929                          const ChannelType channel, __constant float *filter,
 930                          const unsigned int width,
 931                          const unsigned int imageColumns, const unsigned int imageRows,
 932                          __local CLPixelType *temp,
 933                          const unsigned int offsetRows, const unsigned int section)
 934       {
 935         const int x = get_global_id(0);
 936         const int y = get_global_id(1);
 937
 938         const int columns = imageColumns;
 939
 940         const unsigned int radius = (width-1)/2;
 941         const int wsize = get_local_size(0);
 942         const unsigned int loadSize = wsize+width;
 943
 944         //group coordinate
 945         const int groupX=get_local_size(0)*get_group_id(0);
 946         const int groupY=get_local_size(1)*get_group_id(1);
 947
 948         //offset the input data, assuming section is 0, 1
 949         im += imageColumns * (offsetRows - radius * section);
 950
 951         //parallel load and clamp
 952         for (int i=get_local_id(0); i < loadSize; i=i+get_local_size(0))
 953         {
 954           //int cx = ClampToCanvas(groupX+i, columns);
 955           temp[i] = im[y * columns + ClampToCanvas(i+groupX-radius, columns)];
 956
 957           if (0 && y==0 && get_group_id(1) == 0)
 958           {
 959             printf("(%d %d) temp %d load %d groupX %d\n", x, y, i, ClampToCanvas(groupX+i, columns), groupX);
 960           }
 961         }
 962
 963         // barrier
 964         barrier(CLK_LOCAL_MEM_FENCE);
 965
 966         // only do the work if this is not a patched item
 967         if (get_global_id(0) < columns)
 968         {
 969           // compute
 970           float4 result = (float4) 0;
 971
 972           int i = 0;
 973
 974           \n #ifndef UFACTOR   \n
 975           \n #define UFACTOR 8 \n
 976           \n #endif                  \n
 977
 978           for ( ; i+UFACTOR < width; )
 979           {
 980             \n #pragma unroll UFACTOR\n
 981             for (int j=0; j < UFACTOR; j++, i++)
 982             {
 983               result+=filter[i]*convert_float4(temp[i+get_local_id(0)]);
 984             }
 985           }
 986
 987           for ( ; i < width; i++)
 988           {
 989             result+=filter[i]*convert_float4(temp[i+get_local_id(0)]);
 990           }
 991
 992           result.x = ClampToQuantum(result.x);
 993           result.y = ClampToQuantum(result.y);
 994           result.z = ClampToQuantum(result.z);
 995           result.w = ClampToQuantum(result.w);
 996
 997           // write back to global
 998           filtered_im[y*columns+x] = result;
 999         }
1000
1001       }
1002     )
1003
1004     STRINGIFY(
1005       /*
1006       Reduce image noise and reduce detail levels by line
1007       im: input pixels filtered_in  filtered_im: output pixels
1008       filter : convolve kernel  width: convolve kernel size
1009       channel : define which channel is blured\
1010       is_RGBA_BGRA : define the input is RGBA or BGRA
1011       */
1012       __kernel void BlurColumn(const __global float4 *blurRowData, __global CLPixelType *filtered_im,
1013                                 const ChannelType channel, __constant float *filter,
1014                                 const unsigned int width,
1015                                 const unsigned int imageColumns, const unsigned int imageRows,
1016                                 __local float4 *temp)
1017       {
1018         const int x = get_global_id(0);
1019         const int y = get_global_id(1);
1020
1021         //const int columns = get_global_size(0);
1022         //const int rows = get_global_size(1);
1023         const int columns = imageColumns;
1024         const int rows = imageRows;
1025
1026         unsigned int radius = (width-1)/2;
1027         const int wsize = get_local_size(1);
1028         const unsigned int loadSize = wsize+width;
1029
1030         //group coordinate
1031         const int groupX=get_local_size(0)*get_group_id(0);
1032         const int groupY=get_local_size(1)*get_group_id(1);
1033         //notice that get_local_size(0) is 1, so
1034         //groupX=get_group_id(0);
1035
1036         //parallel load and clamp
1037         for (int i = get_local_id(1); i < loadSize; i=i+get_local_size(1))
1038         {
1039           temp[i] = blurRowData[ClampToCanvas(i+groupY-radius, rows) * columns + groupX];
1040         }
1041
1042         // barrier
1043         barrier(CLK_LOCAL_MEM_FENCE);
1044
1045         // only do the work if this is not a patched item
1046         if (get_global_id(1) < rows)
1047         {
1048           // compute
1049           float4 result = (float4) 0;
1050
1051           int i = 0;
1052
1053           \n #ifndef UFACTOR   \n
1054           \n #define UFACTOR 8 \n
1055           \n #endif                  \n
1056
1057           for ( ; i+UFACTOR < width; )
1058           {
1059             \n #pragma unroll UFACTOR \n
1060             for (int j=0; j < UFACTOR; j++, i++)
1061             {
1062               result+=filter[i]*temp[i+get_local_id(1)];
1063             }
1064           }
1065
1066           for ( ; i < width; i++)
1067           {
1068             result+=filter[i]*temp[i+get_local_id(1)];
1069           }
1070
1071           result.x = ClampToQuantum(result.x);
1072           result.y = ClampToQuantum(result.y);
1073           result.z = ClampToQuantum(result.z);
1074           result.w = ClampToQuantum(result.w);
1075
1076           // write back to global
1077           filtered_im[y*columns+x] = (CLPixelType) (result.x,result.y,result.z,result.w);
1078         }
1079
1080       }
1081     )
1082
1083
1084     STRINGIFY(
1085       /*
1086       Reduce image noise and reduce detail levels by line
1087       im: input pixels filtered_in  filtered_im: output pixels
1088       filter : convolve kernel  width: convolve kernel size
1089       channel : define which channel is blured\
1090       is_RGBA_BGRA : define the input is RGBA or BGRA
1091       */
1092       __kernel void BlurColumnSection(const __global float4 *blurRowData, __global CLPixelType *filtered_im,
1093                                 const ChannelType channel, __constant float *filter,
1094                                 const unsigned int width,
1095                                 const unsigned int imageColumns, const unsigned int imageRows,
1096                                 __local float4 *temp,
1097                                 const unsigned int offsetRows, const unsigned int section)
1098       {
1099         const int x = get_global_id(0);
1100         const int y = get_global_id(1);
1101
1102         //const int columns = get_global_size(0);
1103         //const int rows = get_global_size(1);
1104         const int columns = imageColumns;
1105         const int rows = imageRows;
1106
1107         unsigned int radius = (width-1)/2;
1108         const int wsize = get_local_size(1);
1109         const unsigned int loadSize = wsize+width;
1110
1111         //group coordinate
1112         const int groupX=get_local_size(0)*get_group_id(0);
1113         const int groupY=get_local_size(1)*get_group_id(1);
1114         //notice that get_local_size(0) is 1, so
1115         //groupX=get_group_id(0);
1116
1117         // offset the input data
1118         blurRowData += imageColumns * radius * section;
1119
1120         //parallel load and clamp
1121         for (int i = get_local_id(1); i < loadSize; i=i+get_local_size(1))
1122         {
1123           int pos = ClampToCanvasWithHalo(i+groupY-radius, rows, radius, section) * columns + groupX;
1124           temp[i] = *(blurRowData+pos);
1125         }
1126
1127         // barrier
1128         barrier(CLK_LOCAL_MEM_FENCE);
1129
1130         // only do the work if this is not a patched item
1131         if (get_global_id(1) < rows)
1132         {
1133           // compute
1134           float4 result = (float4) 0;
1135
1136           int i = 0;
1137
1138           \n #ifndef UFACTOR   \n
1139           \n #define UFACTOR 8 \n
1140           \n #endif                  \n
1141
1142           for ( ; i+UFACTOR < width; )
1143           {
1144             \n #pragma unroll UFACTOR \n
1145             for (int j=0; j < UFACTOR; j++, i++)
1146             {
1147               result+=filter[i]*temp[i+get_local_id(1)];
1148             }
1149           }
1150           for ( ; i < width; i++)
1151           {
1152             result+=filter[i]*temp[i+get_local_id(1)];
1153           }
1154
1155           result.x = ClampToQuantum(result.x);
1156           result.y = ClampToQuantum(result.y);
1157           result.z = ClampToQuantum(result.z);
1158           result.w = ClampToQuantum(result.w);
1159
1160           // offset the output data
1161           filtered_im += imageColumns * offsetRows;
1162
1163           // write back to global
1164           filtered_im[y*columns+x] = (CLPixelType) (result.x,result.y,result.z,result.w);
1165         }
1166
1167       }
1168     )
1169
1170
1171     STRINGIFY(
1172     __kernel void UnsharpMaskBlurColumn(const __global CLPixelType* inputImage,
1173           const __global float4 *blurRowData, __global CLPixelType *filtered_im,
1174           const unsigned int imageColumns, const unsigned int imageRows,
1175           __local float4* cachedData, __local float* cachedFilter,
1176           const ChannelType channel, const __global float *filter, const unsigned int width,
1177           const float gain, const float threshold)
1178     {
1179       const unsigned int radius = (width-1)/2;
1180
1181       // cache the pixel shared by the workgroup
1182       const int groupX = get_group_id(0);
1183       const int groupStartY = get_group_id(1)*get_local_size(1) - radius;
1184       const int groupStopY = (get_group_id(1)+1)*get_local_size(1) + radius;
1185
1186       if (groupStartY >= 0
1187           && groupStopY < imageRows) {
1188         event_t e = async_work_group_strided_copy(cachedData
1189                                                 ,blurRowData+groupStartY*imageColumns+groupX
1190                                                 ,groupStopY-groupStartY,imageColumns,0);
1191         wait_group_events(1,&e);
1192       }
1193       else {
1194         for (int i = get_local_id(1); i < (groupStopY - groupStartY); i+=get_local_size(1)) {
1195           cachedData[i] = blurRowData[ClampToCanvas(groupStartY+i,imageRows)*imageColumns+ groupX];
1196         }
1197         barrier(CLK_LOCAL_MEM_FENCE);
1198       }
1199       // cache the filter as well
1200       event_t e = async_work_group_copy(cachedFilter,filter,width,0);
1201       wait_group_events(1,&e);
1202
1203       // only do the work if this is not a patched item
1204       //const int cy = get_group_id(1)*get_local_size(1)+get_local_id(1);
1205       const int cy = get_global_id(1);
1206
1207       if (cy < imageRows) {
1208         float4 blurredPixel = (float4) 0.0f;
1209
1210         int i = 0;
1211
1212         \n #ifndef UFACTOR   \n
1213           \n #define UFACTOR 8 \n
1214           \n #endif                  \n
1215
1216           for ( ; i+UFACTOR < width; )
1217           {
1218             \n #pragma unroll UFACTOR \n
1219               for (int j=0; j < UFACTOR; j++, i++)
1220               {
1221                 blurredPixel+=cachedFilter[i]*cachedData[i+get_local_id(1)];
1222               }
1223           }
1224
1225         for ( ; i < width; i++)
1226         {
1227           blurredPixel+=cachedFilter[i]*cachedData[i+get_local_id(1)];
1228         }
1229
1230         blurredPixel = floor((float4)(ClampToQuantum(blurredPixel.x), ClampToQuantum(blurredPixel.y)
1231                                       ,ClampToQuantum(blurredPixel.z), ClampToQuantum(blurredPixel.w)));
1232
1233         float4 inputImagePixel = convert_float4(inputImage[cy*imageColumns+groupX]);
1234         float4 outputPixel = inputImagePixel - blurredPixel;
1235
1236         float quantumThreshold = QuantumRange*threshold;
1237
1238         int4 mask = isless(fabs(2.0f*outputPixel), (float4)quantumThreshold);
1239         outputPixel = select(inputImagePixel + outputPixel * gain, inputImagePixel, mask);
1240
1241         //write back
1242         filtered_im[cy*imageColumns+groupX] = (CLPixelType) (ClampToQuantum(outputPixel.x), ClampToQuantum(outputPixel.y)
1243                                                             ,ClampToQuantum(outputPixel.z), ClampToQuantum(outputPixel.w));
1244
1245       }
1246     }
1247
1248     __kernel void UnsharpMaskBlurColumnSection(const __global CLPixelType* inputImage,
1249           const __global float4 *blurRowData, __global CLPixelType *filtered_im,
1250           const unsigned int imageColumns, const unsigned int imageRows,
1251           __local float4* cachedData, __local float* cachedFilter,
1252           const ChannelType channel, const __global float *filter, const unsigned int width,
1253           const float gain, const float threshold,
1254           const unsigned int offsetRows, const unsigned int section)
1255     {
1256       const unsigned int radius = (width-1)/2;
1257
1258       // cache the pixel shared by the workgroup
1259       const int groupX = get_group_id(0);
1260       const int groupStartY = get_group_id(1)*get_local_size(1) - radius;
1261       const int groupStopY = (get_group_id(1)+1)*get_local_size(1) + radius;
1262
1263       // offset the input data
1264       blurRowData += imageColumns * radius * section;
1265
1266       if (groupStartY >= 0
1267           && groupStopY < imageRows) {
1268         event_t e = async_work_group_strided_copy(cachedData
1269                                                 ,blurRowData+groupStartY*imageColumns+groupX
1270                                                 ,groupStopY-groupStartY,imageColumns,0);
1271         wait_group_events(1,&e);
1272       }
1273       else {
1274         for (int i = get_local_id(1); i < (groupStopY - groupStartY); i+=get_local_size(1)) {
1275           int pos = ClampToCanvasWithHalo(groupStartY+i,imageRows, radius, section)*imageColumns+ groupX;
1276           cachedData[i] = *(blurRowData + pos);
1277         }
1278         barrier(CLK_LOCAL_MEM_FENCE);
1279       }
1280       // cache the filter as well
1281       event_t e = async_work_group_copy(cachedFilter,filter,width,0);
1282       wait_group_events(1,&e);
1283
1284       // only do the work if this is not a patched item
1285       //const int cy = get_group_id(1)*get_local_size(1)+get_local_id(1);
1286       const int cy = get_global_id(1);
1287
1288       if (cy < imageRows) {
1289         float4 blurredPixel = (float4) 0.0f;
1290
1291         int i = 0;
1292
1293         \n #ifndef UFACTOR   \n
1294           \n #define UFACTOR 8 \n
1295           \n #endif                  \n
1296
1297           for ( ; i+UFACTOR < width; )
1298           {
1299             \n #pragma unroll UFACTOR \n
1300               for (int j=0; j < UFACTOR; j++, i++)
1301               {
1302                 blurredPixel+=cachedFilter[i]*cachedData[i+get_local_id(1)];
1303               }
1304           }
1305
1306         for ( ; i < width; i++)
1307         {
1308           blurredPixel+=cachedFilter[i]*cachedData[i+get_local_id(1)];
1309         }
1310
1311         blurredPixel = floor((float4)(ClampToQuantum(blurredPixel.x), ClampToQuantum(blurredPixel.y)
1312                                       ,ClampToQuantum(blurredPixel.z), ClampToQuantum(blurredPixel.w)));
1313
1314         // offset the output data
1315         inputImage += imageColumns * offsetRows;
1316         filtered_im += imageColumns * offsetRows;
1317
1318         float4 inputImagePixel = convert_float4(inputImage[cy*imageColumns+groupX]);
1319         float4 outputPixel = inputImagePixel - blurredPixel;
1320
1321         float quantumThreshold = QuantumRange*threshold;
1322
1323         int4 mask = isless(fabs(2.0f*outputPixel), (float4)quantumThreshold);
1324         outputPixel = select(inputImagePixel + outputPixel * gain, inputImagePixel, mask);
1325
1326         //write back
1327         filtered_im[cy*imageColumns+groupX] = (CLPixelType) (ClampToQuantum(outputPixel.x), ClampToQuantum(outputPixel.y)
1328                                                             ,ClampToQuantum(outputPixel.z), ClampToQuantum(outputPixel.w));
1329
1330       }
1331
1332     }
1333     )
1334
1335
1336
1337   STRINGIFY(
1338
1339   __kernel void HullPass1(const __global CLPixelType *inputImage, __global CLPixelType *outputImage
1340   , const unsigned int imageWidth, const unsigned int imageHeight
1341   , const int2 offset, const int polarity, const int matte) {
1342
1343     int x = get_global_id(0);
1344     int y = get_global_id(1);
1345
1346     CLPixelType v = inputImage[y*imageWidth+x];
1347
1348     int2 neighbor;
1349     neighbor.y = y + offset.y;
1350     neighbor.x = x + offset.x;
1351
1352     int2 clampedNeighbor;
1353     clampedNeighbor.x = ClampToCanvas(neighbor.x, imageWidth);
1354     clampedNeighbor.y = ClampToCanvas(neighbor.y, imageHeight);
1355
1356     CLPixelType r = (clampedNeighbor.x == neighbor.x
1357                      && clampedNeighbor.y == neighbor.y)?inputImage[clampedNeighbor.y*imageWidth+clampedNeighbor.x]
1358     :(CLPixelType)0;
1359
1360     int sv[4];
1361     sv[0] = (int)v.x;
1362     sv[1] = (int)v.y;
1363     sv[2] = (int)v.z;
1364     sv[3] = (int)v.w;
1365
1366     int sr[4];
1367     sr[0] = (int)r.x;
1368     sr[1] = (int)r.y;
1369     sr[2] = (int)r.z;
1370     sr[3] = (int)r.w;
1371
1372     if (polarity > 0) {
1373       \n #pragma unroll 4\n
1374       for (unsigned int i = 0; i < 4; i++) {
1375         sv[i] = (sr[i] >= (sv[i]+ScaleCharToQuantum(2)))?(sv[i]+ScaleCharToQuantum(1)):sv[i];
1376       }
1377     }
1378     else {
1379       \n #pragma unroll 4\n
1380       for (unsigned int i = 0; i < 4; i++) {
1381         sv[i] = (sr[i] <= (sv[i]-ScaleCharToQuantum(2)))?(sv[i]-ScaleCharToQuantum(1)):sv[i];
1382       }
1383
1384     }
1385
1386     v.x = (CLQuantum)sv[0];
1387     v.y = (CLQuantum)sv[1];
1388     v.z = (CLQuantum)sv[2];
1389
1390     if (matte!=0)
1391       v.w = (CLQuantum)sv[3];
1392
1393     outputImage[y*imageWidth+x] = v;
1394
1395     }
1396
1397
1398   )
1399
1400
1401
1402   STRINGIFY(
1403
1404   __kernel void HullPass2(const __global CLPixelType *inputImage, __global CLPixelType *outputImage
1405   , const unsigned int imageWidth, const unsigned int imageHeight
1406   , const int2 offset, const int polarity, const int matte) {
1407
1408     int x = get_global_id(0);
1409     int y = get_global_id(1);
1410
1411     CLPixelType v = inputImage[y*imageWidth+x];
1412
1413     int2 neighbor, clampedNeighbor;
1414
1415     neighbor.y = y + offset.y;
1416     neighbor.x = x + offset.x;
1417     clampedNeighbor.x = ClampToCanvas(neighbor.x, imageWidth);
1418     clampedNeighbor.y = ClampToCanvas(neighbor.y, imageHeight);
1419
1420     CLPixelType r = (clampedNeighbor.x == neighbor.x
1421       && clampedNeighbor.y == neighbor.y)?inputImage[clampedNeighbor.y*imageWidth+clampedNeighbor.x]
1422     :(CLPixelType)0;
1423
1424
1425     neighbor.y = y - offset.y;
1426     neighbor.x = x - offset.x;
1427     clampedNeighbor.x = ClampToCanvas(neighbor.x, imageWidth);
1428     clampedNeighbor.y = ClampToCanvas(neighbor.y, imageHeight);
1429
1430     CLPixelType s = (clampedNeighbor.x == neighbor.x
1431       && clampedNeighbor.y == neighbor.y)?inputImage[clampedNeighbor.y*imageWidth+clampedNeighbor.x]
1432     :(CLPixelType)0;
1433
1434
1435     int sv[4];
1436     sv[0] = (int)v.x;
1437     sv[1] = (int)v.y;
1438     sv[2] = (int)v.z;
1439     sv[3] = (int)v.w;
1440
1441     int sr[4];
1442     sr[0] = (int)r.x;
1443     sr[1] = (int)r.y;
1444     sr[2] = (int)r.z;
1445     sr[3] = (int)r.w;
1446
1447     int ss[4];
1448     ss[0] = (int)s.x;
1449     ss[1] = (int)s.y;
1450     ss[2] = (int)s.z;
1451     ss[3] = (int)s.w;
1452
1453     if (polarity > 0) {
1454       \n #pragma unroll 4\n
1455       for (unsigned int i = 0; i < 4; i++) {
1456         //sv[i] = (ss[i] >= (sv[i]+ScaleCharToQuantum(2)) && sr[i] > sv[i] )   ? (sv[i]+ScaleCharToQuantum(1)):sv[i];
1457         //
1458         //sv[i] =(!( (int)(ss[i] >= (sv[i]+ScaleCharToQuantum(2))) && (int) (sr[i] > sv[i] ) ))  ? sv[i]:(sv[i]+ScaleCharToQuantum(1));
1459         //sv[i] =(( (int)( ss[i] < (sv[i]+ScaleCharToQuantum(2))) || (int) ( sr[i] <= sv[i] ) ))  ? sv[i]:(sv[i]+ScaleCharToQuantum(1));
1460         sv[i] =(( (int)( ss[i] < (sv[i]+ScaleCharToQuantum(2))) + (int) ( sr[i] <= sv[i] ) ) !=0)  ? sv[i]:(sv[i]+ScaleCharToQuantum(1));
1461       }
1462     }
1463     else {
1464       \n #pragma unroll 4\n
1465       for (unsigned int i = 0; i < 4; i++) {
1466         //sv[i] = (ss[i] <= (sv[i]-ScaleCharToQuantum(2)) && sr[i] < sv[i] )   ? (sv[i]-ScaleCharToQuantum(1)):sv[i];
1467         //
1468         //sv[i] = ( (int)(ss[i] <= (sv[i]-ScaleCharToQuantum(2)) ) + (int)( sr[i] < sv[i] ) ==0)   ? sv[i]:(sv[i]-ScaleCharToQuantum(1));
1469         sv[i] = (( (int)(ss[i] > (sv[i]-ScaleCharToQuantum(2))) + (int)( sr[i] >= sv[i] )) !=0)   ? sv[i]:(sv[i]-ScaleCharToQuantum(1));
1470       }
1471     }
1472
1473     v.x = (CLQuantum)sv[0];
1474     v.y = (CLQuantum)sv[1];
1475     v.z = (CLQuantum)sv[2];
1476
1477     if (matte!=0)
1478       v.w = (CLQuantum)sv[3];
1479
1480     outputImage[y*imageWidth+x] = v;
1481
1482     }
1483
1484
1485   )
1486
1487
1488   STRINGIFY(
1489     __kernel void RotationalBlur(const __global CLPixelType *im, __global CLPixelType *filtered_im,
1490                                  const float4 bias,
1491                                  const unsigned int channel, const unsigned int matte,
1492                                  const float2 blurCenter,
1493                                  __constant float *cos_theta, __constant float *sin_theta,
1494                                  const unsigned int cossin_theta_size)
1495       {
1496         const int x = get_global_id(0);
1497         const int y = get_global_id(1);
1498         const int columns = get_global_size(0);
1499         const int rows = get_global_size(1);
1500         unsigned int step = 1;
1501         float center_x = (float) x - blurCenter.x;
1502         float center_y = (float) y - blurCenter.y;
1503         float radius = hypot(center_x, center_y);
1504
1505         //float blur_radius = hypot((float) columns/2.0f, (float) rows/2.0f);
1506         float blur_radius = hypot(blurCenter.x, blurCenter.y);
1507
1508         if (radius > MagickEpsilon)
1509         {
1510           step = (unsigned int) (blur_radius / radius);
1511           if (step == 0)
1512             step = 1;
1513           if (step >= cossin_theta_size)
1514             step = cossin_theta_size-1;
1515         }
1516
1517         float4 result;
1518         result.x = (float)bias.x;
1519         result.y = (float)bias.y;
1520         result.z = (float)bias.z;
1521         result.w = (float)bias.w;
1522         float normalize = 0.0f;
1523
1524         if (((channel & OpacityChannel) == 0) || (matte == 0)) {
1525           for (unsigned int i=0; i<cossin_theta_size; i+=step)
1526           {
1527             result += convert_float4(im[
1528               ClampToCanvas(blurCenter.x+center_x*cos_theta[i]-center_y*sin_theta[i]+0.5f,columns)+
1529                 ClampToCanvas(blurCenter.y+center_x*sin_theta[i]+center_y*cos_theta[i]+0.5f, rows)*columns]);
1530               normalize += 1.0f;
1531           }
1532           normalize = PerceptibleReciprocal(normalize);
1533           result = result * normalize;
1534         }
1535         else {
1536           float gamma = 0.0f;
1537           for (unsigned int i=0; i<cossin_theta_size; i+=step)
1538           {
1539             float4 p = convert_float4(im[
1540               ClampToCanvas(blurCenter.x+center_x*cos_theta[i]-center_y*sin_theta[i]+0.5f,columns)+
1541                 ClampToCanvas(blurCenter.y+center_x*sin_theta[i]+center_y*cos_theta[i]+0.5f, rows)*columns]);
1542
1543             float alpha = (float)(QuantumScale*(QuantumRange-p.w));
1544             result.x += alpha * p.x;
1545             result.y += alpha * p.y;
1546             result.z += alpha * p.z;
1547             result.w += p.w;
1548             gamma+=alpha;
1549             normalize += 1.0f;
1550           }
1551           gamma = PerceptibleReciprocal(gamma);
1552           normalize = PerceptibleReciprocal(normalize);
1553           result.x = gamma*result.x;
1554           result.y = gamma*result.y;
1555           result.z = gamma*result.z;
1556           result.w = normalize*result.w;
1557         }
1558         filtered_im[y * columns + x] = (CLPixelType) (ClampToQuantum(result.x), ClampToQuantum(result.y),
1559           ClampToQuantum(result.z), ClampToQuantum(result.w));
1560       }
1561   )
1562
1563   STRINGIFY(
1564
1565   inline float3 ConvertRGBToHSB(CLPixelType pixel) {
1566     float3 HueSaturationBrightness;
1567     HueSaturationBrightness.x = 0.0f; // Hue
1568     HueSaturationBrightness.y = 0.0f; // Saturation
1569     HueSaturationBrightness.z = 0.0f; // Brightness
1570
1571     float r=(float) getRed(pixel);
1572     float g=(float) getGreen(pixel);
1573     float b=(float) getBlue(pixel);
1574
1575     float tmin=min(min(r,g),b);
1576     float tmax=max(max(r,g),b);
1577
1578     if (tmax!=0.0f) {
1579       float delta=tmax-tmin;
1580       HueSaturationBrightness.y=delta/tmax;
1581       HueSaturationBrightness.z=QuantumScale*tmax;
1582
1583       if (delta != 0.0f) {
1584   HueSaturationBrightness.x = ((r == tmax)?0.0f:((g == tmax)?2.0f:4.0f));
1585   HueSaturationBrightness.x += ((r == tmax)?(g-b):((g == tmax)?(b-r):(r-g)))/delta;
1586         HueSaturationBrightness.x/=6.0f;
1587         HueSaturationBrightness.x += (HueSaturationBrightness.x < 0.0f)?0.0f:1.0f;
1588       }
1589     }
1590     return HueSaturationBrightness;
1591   }
1592
1593   inline CLPixelType ConvertHSBToRGB(float3 HueSaturationBrightness) {
1594
1595     float hue = HueSaturationBrightness.x;
1596     float brightness = HueSaturationBrightness.z;
1597     float saturation = HueSaturationBrightness.y;
1598
1599     CLPixelType rgb;
1600
1601     if (saturation == 0.0f) {
1602       setRed(&rgb,ClampToQuantum(QuantumRange*brightness));
1603       setGreen(&rgb,getRed(rgb));
1604       setBlue(&rgb,getRed(rgb));
1605     }
1606     else {
1607
1608       float h=6.0f*(hue-floor(hue));
1609       float f=h-floor(h);
1610       float p=brightness*(1.0f-saturation);
1611       float q=brightness*(1.0f-saturation*f);
1612       float t=brightness*(1.0f-(saturation*(1.0f-f)));
1613
1614       float clampedBrightness = ClampToQuantum(QuantumRange*brightness);
1615       float clamped_t = ClampToQuantum(QuantumRange*t);
1616       float clamped_p = ClampToQuantum(QuantumRange*p);
1617       float clamped_q = ClampToQuantum(QuantumRange*q);
1618       int ih = (int)h;
1619       setRed(&rgb, (ih == 1)?clamped_q:
1620         (ih == 2 || ih == 3)?clamped_p:
1621         (ih == 4)?clamped_t:
1622                  clampedBrightness);
1623
1624       setGreen(&rgb, (ih == 1 || ih == 2)?clampedBrightness:
1625         (ih == 3)?clamped_q:
1626         (ih == 4 || ih == 5)?clamped_p:
1627                  clamped_t);
1628
1629       setBlue(&rgb, (ih == 2)?clamped_t:
1630         (ih == 3 || ih == 4)?clampedBrightness:
1631         (ih == 5)?clamped_q:
1632                  clamped_p);
1633     }
1634     return rgb;
1635   }
1636
1637   __kernel void Contrast(__global CLPixelType *im, const unsigned int sharpen)
1638   {
1639
1640     const int sign = sharpen!=0?1:-1;
1641     const int x = get_global_id(0);
1642     const int y = get_global_id(1);
1643     const int columns = get_global_size(0);
1644     const int c = x + y * columns;
1645
1646     CLPixelType pixel = im[c];
1647     float3 HueSaturationBrightness = ConvertRGBToHSB(pixel);
1648     float brightness = HueSaturationBrightness.z;
1649     brightness+=0.5f*sign*(0.5f*(sinpi(brightness-0.5f)+1.0f)-brightness);
1650     brightness = clamp(brightness,0.0f,1.0f);
1651     HueSaturationBrightness.z = brightness;
1652
1653     CLPixelType filteredPixel = ConvertHSBToRGB(HueSaturationBrightness);
1654     filteredPixel.w = pixel.w;
1655     im[c] = filteredPixel;
1656   }
1657
1658
1659   )
1660
1661   STRINGIFY(
1662
1663   inline void ConvertRGBToHSL(const CLQuantum red,const CLQuantum green, const CLQuantum blue,
1664     float *hue, float *saturation, float *lightness)
1665   {
1666   float
1667     c,
1668     tmax,
1669     tmin;
1670
1671   /*
1672      Convert RGB to HSL colorspace.
1673      */
1674   tmax=max(QuantumScale*red,max(QuantumScale*green, QuantumScale*blue));
1675   tmin=min(QuantumScale*red,min(QuantumScale*green, QuantumScale*blue));
1676
1677   c=tmax-tmin;
1678
1679   *lightness=(tmax+tmin)/2.0;
1680   if (c <= 0.0)
1681   {
1682     *hue=0.0;
1683     *saturation=0.0;
1684     return;
1685   }
1686
1687   if (tmax == (QuantumScale*red))
1688   {
1689     *hue=(QuantumScale*green-QuantumScale*blue)/c;
1690     if ((QuantumScale*green) < (QuantumScale*blue))
1691       *hue+=6.0;
1692   }
1693   else
1694     if (tmax == (QuantumScale*green))
1695       *hue=2.0+(QuantumScale*blue-QuantumScale*red)/c;
1696     else
1697       *hue=4.0+(QuantumScale*red-QuantumScale*green)/c;
1698
1699   *hue*=60.0/360.0;
1700   if (*lightness <= 0.5)
1701     *saturation=c/(2.0*(*lightness));
1702   else
1703     *saturation=c/(2.0-2.0*(*lightness));
1704   }
1705
1706   inline void ConvertHSLToRGB(const float hue,const float saturation, const float lightness,
1707       CLQuantum *red,CLQuantum *green,CLQuantum *blue)
1708   {
1709     float
1710       b,
1711       c,
1712       g,
1713       h,
1714       tmin,
1715       r,
1716       x;
1717
1718     /*
1719        Convert HSL to RGB colorspace.
1720        */
1721     h=hue*360.0;
1722     if (lightness <= 0.5)
1723       c=2.0*lightness*saturation;
1724     else
1725       c=(2.0-2.0*lightness)*saturation;
1726     tmin=lightness-0.5*c;
1727     h-=360.0*floor(h/360.0);
1728     h/=60.0;
1729     x=c*(1.0-fabs(h-2.0*floor(h/2.0)-1.0));
1730     switch ((int) floor(h))
1731     {
1732       case 0:
1733         {
1734           r=tmin+c;
1735           g=tmin+x;
1736           b=tmin;
1737           break;
1738         }
1739       case 1:
1740         {
1741           r=tmin+x;
1742           g=tmin+c;
1743           b=tmin;
1744           break;
1745         }
1746       case 2:
1747         {
1748           r=tmin;
1749           g=tmin+c;
1750           b=tmin+x;
1751           break;
1752         }
1753       case 3:
1754         {
1755           r=tmin;
1756           g=tmin+x;
1757           b=tmin+c;
1758           break;
1759         }
1760       case 4:
1761         {
1762           r=tmin+x;
1763           g=tmin;
1764           b=tmin+c;
1765           break;
1766         }
1767       case 5:
1768         {
1769           r=tmin+c;
1770           g=tmin;
1771           b=tmin+x;
1772           break;
1773         }
1774       default:
1775         {
1776           r=0.0;
1777           g=0.0;
1778           b=0.0;
1779         }
1780     }
1781     *red=ClampToQuantum(QuantumRange*r);
1782     *green=ClampToQuantum(QuantumRange*g);
1783     *blue=ClampToQuantum(QuantumRange*b);
1784   }
1785
1786   inline void ModulateHSL(const float percent_hue, const float percent_saturation,const float percent_lightness,
1787     CLQuantum *red,CLQuantum *green,CLQuantum *blue)
1788   {
1789     float
1790       hue,
1791       lightness,
1792       saturation;
1793
1794     /*
1795     Increase or decrease color lightness, saturation, or hue.
1796     */
1797     ConvertRGBToHSL(*red,*green,*blue,&hue,&saturation,&lightness);
1798     hue+=0.5*(0.01*percent_hue-1.0);
1799     while (hue < 0.0)
1800       hue+=1.0;
1801     while (hue >= 1.0)
1802       hue-=1.0;
1803     saturation*=0.01*percent_saturation;
1804     lightness*=0.01*percent_lightness;
1805     ConvertHSLToRGB(hue,saturation,lightness,red,green,blue);
1806   }
1807
1808   __kernel void Modulate(__global CLPixelType *im,
1809     const float percent_brightness,
1810     const float percent_hue,
1811     const float percent_saturation,
1812     const int colorspace)
1813   {
1814
1815     const int x = get_global_id(0);
1816     const int y = get_global_id(1);
1817     const int columns = get_global_size(0);
1818     const int c = x + y * columns;
1819
1820     CLPixelType pixel = im[c];
1821
1822     CLQuantum
1823         blue,
1824         green,
1825         red;
1826
1827     red=getRed(pixel);
1828     green=getGreen(pixel);
1829     blue=getBlue(pixel);
1830
1831     switch (colorspace)
1832     {
1833       case HSLColorspace:
1834       default:
1835         {
1836           ModulateHSL(percent_hue, percent_saturation, percent_brightness,
1837               &red, &green, &blue);
1838         }
1839
1840     }
1841
1842     CLPixelType filteredPixel;
1843
1844     setRed(&filteredPixel, red);
1845     setGreen(&filteredPixel, green);
1846     setBlue(&filteredPixel, blue);
1847     filteredPixel.w = pixel.w;
1848
1849     im[c] = filteredPixel;
1850   }
1851   )
1852
1853   STRINGIFY(
1854   __kernel void Negate(__global CLPixelType *im,
1855     const ChannelType channel)
1856   {
1857
1858     const int x = get_global_id(0);
1859     const int y = get_global_id(1);
1860     const int columns = get_global_size(0);
1861     const int c = x + y * columns;
1862
1863     CLPixelType pixel = im[c];
1864
1865     CLQuantum
1866         blue,
1867         green,
1868         red;
1869
1870     red=getRed(pixel);
1871     green=getGreen(pixel);
1872     blue=getBlue(pixel);
1873
1874     CLPixelType filteredPixel;
1875
1876     if ((channel & RedChannel) !=0)
1877       setRed(&filteredPixel, QuantumRange-red);
1878     if ((channel & GreenChannel) !=0)
1879       setGreen(&filteredPixel, QuantumRange-green);
1880     if ((channel & BlueChannel) !=0)
1881       setBlue(&filteredPixel, QuantumRange-blue);
1882
1883     filteredPixel.w = pixel.w;
1884
1885     im[c] = filteredPixel;
1886   }
1887   )
1888
1889   STRINGIFY(
1890   __kernel void Grayscale(__global CLPixelType *im,
1891     const int method, const int colorspace)
1892   {
1893
1894     const int x = get_global_id(0);
1895     const int y = get_global_id(1);
1896     const int columns = get_global_size(0);
1897     const int c = x + y * columns;
1898
1899     CLPixelType pixel = im[c];
1900
1901     float
1902         blue,
1903         green,
1904         intensity,
1905         red;
1906
1907     red=(float)getRed(pixel);
1908     green=(float)getGreen(pixel);
1909     blue=(float)getBlue(pixel);
1910
1911     intensity=0.0;
1912
1913     CLPixelType filteredPixel;
1914
1915     switch (method)
1916     {
1917       case AveragePixelIntensityMethod:
1918         {
1919           intensity=(red+green+blue)/3.0;
1920           break;
1921         }
1922       case BrightnessPixelIntensityMethod:
1923         {
1924           intensity=max(max(red,green),blue);
1925           break;
1926         }
1927       case LightnessPixelIntensityMethod:
1928         {
1929           intensity=(min(min(red,green),blue)+
1930               max(max(red,green),blue))/2.0;
1931           break;
1932         }
1933       case MSPixelIntensityMethod:
1934         {
1935           intensity=(float) (((float) red*red+green*green+
1936                 blue*blue)/(3.0*QuantumRange));
1937           break;
1938         }
1939       case Rec601LumaPixelIntensityMethod:
1940         {
1941           /*
1942           if (colorspace == RGBColorspace)
1943           {
1944             red=EncodePixelGamma(red);
1945             green=EncodePixelGamma(green);
1946             blue=EncodePixelGamma(blue);
1947           }
1948           */
1949           intensity=0.298839*red+0.586811*green+0.114350*blue;
1950           break;
1951         }
1952       case Rec601LuminancePixelIntensityMethod:
1953         {
1954           /*
1955           if (image->colorspace == sRGBColorspace)
1956           {
1957             red=DecodePixelGamma(red);
1958             green=DecodePixelGamma(green);
1959             blue=DecodePixelGamma(blue);
1960           }
1961           */
1962           intensity=0.298839*red+0.586811*green+0.114350*blue;
1963           break;
1964         }
1965       case Rec709LumaPixelIntensityMethod:
1966       default:
1967         {
1968           /*
1969           if (image->colorspace == RGBColorspace)
1970           {
1971             red=EncodePixelGamma(red);
1972             green=EncodePixelGamma(green);
1973             blue=EncodePixelGamma(blue);
1974           }
1975           */
1976           intensity=0.212656*red+0.715158*green+0.072186*blue;
1977           break;
1978         }
1979       case Rec709LuminancePixelIntensityMethod:
1980         {
1981           /*
1982           if (image->colorspace == sRGBColorspace)
1983           {
1984             red=DecodePixelGamma(red);
1985             green=DecodePixelGamma(green);
1986             blue=DecodePixelGamma(blue);
1987           }
1988           */
1989           intensity=0.212656*red+0.715158*green+0.072186*blue;
1990           break;
1991         }
1992       case RMSPixelIntensityMethod:
1993         {
1994           intensity=(float) (sqrt((float) red*red+green*green+
1995                 blue*blue)/sqrt(3.0));
1996           break;
1997         }
1998
1999     }
2000
2001     setGray(&filteredPixel, ClampToQuantum(intensity));
2002
2003     filteredPixel.w = pixel.w;
2004
2005     im[c] = filteredPixel;
2006   }
2007   )
2008
2009   STRINGIFY(
2010   // Based on Box from resize.c
2011   float BoxResizeFilter(const float x)
2012   {
2013     return 1.0f;
2014   }
2015   )
2016
2017   STRINGIFY(
2018   // Based on CubicBC from resize.c
2019   float CubicBC(const float x,const __global float* resizeFilterCoefficients)
2020   {
2021     /*
2022     Cubic Filters using B,C determined values:
2023     Mitchell-Netravali  B = 1/3 C = 1/3  "Balanced" cubic spline filter
2024     Catmull-Rom         B = 0   C = 1/2  Interpolatory and exact on linears
2025     Spline              B = 1   C = 0    B-Spline Gaussian approximation
2026     Hermite             B = 0   C = 0    B-Spline interpolator
2027
2028     See paper by Mitchell and Netravali, Reconstruction Filters in Computer
2029     Graphics Computer Graphics, Volume 22, Number 4, August 1988
2030     http://www.cs.utexas.edu/users/fussell/courses/cs384g/lectures/mitchell/
2031     Mitchell.pdf.
2032
2033     Coefficents are determined from B,C values:
2034     P0 = (  6 - 2*B       )/6 = coeff[0]
2035     P1 =         0
2036     P2 = (-18 +12*B + 6*C )/6 = coeff[1]
2037     P3 = ( 12 - 9*B - 6*C )/6 = coeff[2]
2038     Q0 = (      8*B +24*C )/6 = coeff[3]
2039     Q1 = (    -12*B -48*C )/6 = coeff[4]
2040     Q2 = (      6*B +30*C )/6 = coeff[5]
2041     Q3 = (    - 1*B - 6*C )/6 = coeff[6]
2042
2043     which are used to define the filter:
2044
2045     P0 + P1*x + P2*x^2 + P3*x^3      0 <= x < 1
2046     Q0 + Q1*x + Q2*x^2 + Q3*x^3      1 <= x < 2
2047
2048     which ensures function is continuous in value and derivative (slope).
2049     */
2050     if (x < 1.0)
2051       return(resizeFilterCoefficients[0]+x*(x*
2052       (resizeFilterCoefficients[1]+x*resizeFilterCoefficients[2])));
2053     if (x < 2.0)
2054       return(resizeFilterCoefficients[3]+x*(resizeFilterCoefficients[4]+x*
2055       (resizeFilterCoefficients[5]+x*resizeFilterCoefficients[6])));
2056     return(0.0);
2057   }
2058   )
2059
2060   STRINGIFY(
2061   float Sinc(const float x)
2062   {
2063     if (x != 0.0f)
2064     {
2065       const float alpha=(float) (MagickPI*x);
2066       return sinpi(x)/alpha;
2067     }
2068     return(1.0f);
2069   }
2070   )
2071
2072   STRINGIFY(
2073   float Triangle(const float x)
2074   {
2075     /*
2076     1st order (linear) B-Spline, bilinear interpolation, Tent 1D filter, or
2077     a Bartlett 2D Cone filter.  Also used as a Bartlett Windowing function
2078     for Sinc().
2079     */
2080     return ((x<1.0f)?(1.0f-x):0.0f);
2081   }
2082   )
2083
2084
2085   STRINGIFY(
2086   float Hanning(const float x)
2087   {
2088     /*
2089     Cosine window function:
2090       0.5+0.5*cos(pi*x).
2091     */
2092     const float cosine=cos((MagickPI*x));
2093     return(0.5f+0.5f*cosine);
2094   }
2095   )
2096
2097   STRINGIFY(
2098   float Hamming(const float x)
2099   {
2100     /*
2101       Offset cosine window function:
2102        .54 + .46 cos(pi x).
2103     */
2104     const float cosine=cos((MagickPI*x));
2105     return(0.54f+0.46f*cosine);
2106   }
2107   )
2108
2109   STRINGIFY(
2110   float Blackman(const float x)
2111   {
2112     /*
2113       Blackman: 2nd order cosine windowing function:
2114         0.42 + 0.5 cos(pi x) + 0.08 cos(2pi x)
2115
2116       Refactored by Chantal Racette and Nicolas Robidoux to one trig call and
2117       five flops.
2118     */
2119     const float cosine=cos((MagickPI*x));
2120     return(0.34f+cosine*(0.5f+cosine*0.16f));
2121   }
2122   )
2123
2124
2125   STRINGIFY(
2126   typedef enum {
2127     BoxWeightingFunction = 0,
2128     TriangleWeightingFunction,
2129     CubicBCWeightingFunction,
2130     HanningWeightingFunction,
2131     HammingWeightingFunction,
2132     BlackmanWeightingFunction,
2133     GaussianWeightingFunction,
2134     QuadraticWeightingFunction,
2135     JincWeightingFunction,
2136     SincWeightingFunction,
2137     SincFastWeightingFunction,
2138     KaiserWeightingFunction,
2139     WelshWeightingFunction,
2140     BohmanWeightingFunction,
2141     LagrangeWeightingFunction,
2142     CosineWeightingFunction,
2143   } ResizeWeightingFunctionType;
2144   )
2145
2146   STRINGIFY(
2147   inline float applyResizeFilter(const float x, const ResizeWeightingFunctionType filterType, const __global float* filterCoefficients)
2148   {
2149     switch (filterType)
2150     {
2151     /* Call Sinc even for SincFast to get better precision on GPU
2152        and to avoid thread divergence.  Sinc is pretty fast on GPU anyway...*/
2153     case SincWeightingFunction:
2154     case SincFastWeightingFunction:
2155       return Sinc(x);
2156     case CubicBCWeightingFunction:
2157       return CubicBC(x,filterCoefficients);
2158     case BoxWeightingFunction:
2159       return BoxResizeFilter(x);
2160     case TriangleWeightingFunction:
2161       return Triangle(x);
2162     case HanningWeightingFunction:
2163       return Hanning(x);
2164     case HammingWeightingFunction:
2165       return Hamming(x);
2166     case BlackmanWeightingFunction:
2167       return Blackman(x);
2168
2169     default:
2170       return 0.0f;
2171     }
2172   }
2173   )
2174
2175
2176   STRINGIFY(
2177   inline float getResizeFilterWeight(const __global float* resizeFilterCubicCoefficients, const ResizeWeightingFunctionType resizeFilterType
2178            , const ResizeWeightingFunctionType resizeWindowType
2179            , const float resizeFilterScale, const float resizeWindowSupport, const float resizeFilterBlur, const float x)
2180   {
2181     float scale;
2182     float xBlur = fabs(x/resizeFilterBlur);
2183     if (resizeWindowSupport < MagickEpsilon
2184         || resizeWindowType == BoxWeightingFunction)
2185     {
2186       scale = 1.0f;
2187     }
2188     else
2189     {
2190       scale = resizeFilterScale;
2191       scale = applyResizeFilter(xBlur*scale, resizeWindowType, resizeFilterCubicCoefficients);
2192     }
2193     float weight = scale * applyResizeFilter(xBlur, resizeFilterType, resizeFilterCubicCoefficients);
2194     return weight;
2195   }
2196
2197   )
2198
2199   ;
2200   const char* accelerateKernels2 =
2201
2202   STRINGIFY(
2203
2204   inline unsigned int getNumWorkItemsPerPixel(const unsigned int pixelPerWorkgroup, const unsigned int numWorkItems) {
2205     return (numWorkItems/pixelPerWorkgroup);
2206   }
2207
2208   // returns the index of the pixel for the current workitem to compute.
2209   // returns -1 if this workitem doesn't need to participate in any computation
2210   inline int pixelToCompute(const unsigned itemID, const unsigned int pixelPerWorkgroup, const unsigned int numWorkItems) {
2211     const unsigned int numWorkItemsPerPixel = getNumWorkItemsPerPixel(pixelPerWorkgroup, numWorkItems);
2212     int pixelIndex = itemID/numWorkItemsPerPixel;
2213     pixelIndex = (pixelIndex<pixelPerWorkgroup)?pixelIndex:-1;
2214     return pixelIndex;
2215   }
2216
2217   )
2218
2219   STRINGIFY(
2220  __kernel __attribute__((reqd_work_group_size(256, 1, 1)))
2221  void ResizeHorizontalFilter(const __global CLPixelType* inputImage, const unsigned int inputColumns, const unsigned int inputRows, const unsigned int matte
2222   , const float xFactor, __global CLPixelType* filteredImage, const unsigned int filteredColumns, const unsigned int filteredRows
2223   , const int resizeFilterType, const int resizeWindowType
2224   , const __global float* resizeFilterCubicCoefficients
2225   , const float resizeFilterScale, const float resizeFilterSupport, const float resizeFilterWindowSupport, const float resizeFilterBlur
2226   , __local CLPixelType* inputImageCache, const int numCachedPixels, const unsigned int pixelPerWorkgroup, const unsigned int pixelChunkSize
2227   , __local float4* outputPixelCache, __local float* densityCache, __local float* gammaCache) {
2228
2229
2230     // calculate the range of resized image pixels computed by this workgroup
2231     const unsigned int startX = get_group_id(0)*pixelPerWorkgroup;
2232     const unsigned int stopX = min(startX + pixelPerWorkgroup,filteredColumns);
2233     const unsigned int actualNumPixelToCompute = stopX - startX;
2234
2235     // calculate the range of input image pixels to cache
2236     float scale = max(1.0f/xFactor+MagickEpsilon ,1.0f);
2237     const float support = max(scale*resizeFilterSupport,0.5f);
2238     scale = PerceptibleReciprocal(scale);
2239
2240     const int cacheRangeStartX = max((int)((startX+0.5f)/xFactor+MagickEpsilon-support+0.5f),(int)(0));
2241     const int cacheRangeEndX = min((int)(cacheRangeStartX + numCachedPixels), (int)inputColumns);
2242
2243     // cache the input pixels into local memory
2244     const unsigned int y = get_global_id(1);
2245     event_t e = async_work_group_copy(inputImageCache,inputImage+y*inputColumns+cacheRangeStartX,cacheRangeEndX-cacheRangeStartX,0);
2246     wait_group_events(1,&e);
2247
2248     unsigned int totalNumChunks = (actualNumPixelToCompute+pixelChunkSize-1)/pixelChunkSize;
2249     for (unsigned int chunk = 0; chunk < totalNumChunks; chunk++)
2250     {
2251
2252       const unsigned int chunkStartX = startX + chunk*pixelChunkSize;
2253       const unsigned int chunkStopX = min(chunkStartX + pixelChunkSize, stopX);
2254       const unsigned int actualNumPixelInThisChunk = chunkStopX - chunkStartX;
2255
2256       // determine which resized pixel computed by this workitem
2257       const unsigned int itemID = get_local_id(0);
2258       const unsigned int numItems = getNumWorkItemsPerPixel(actualNumPixelInThisChunk, get_local_size(0));
2259
2260       const int pixelIndex = pixelToCompute(itemID, actualNumPixelInThisChunk, get_local_size(0));
2261
2262       float4 filteredPixel = (float4)0.0f;
2263       float density = 0.0f;
2264       float gamma = 0.0f;
2265       // -1 means this workitem doesn't participate in the computation
2266       if (pixelIndex != -1) {
2267
2268         // x coordinated of the resized pixel computed by this workitem
2269         const int x = chunkStartX + pixelIndex;
2270
2271         // calculate how many steps required for this pixel
2272         const float bisect = (x+0.5)/xFactor+MagickEpsilon;
2273         const unsigned int start = (unsigned int)max(bisect-support+0.5f,0.0f);
2274         const unsigned int stop  = (unsigned int)min(bisect+support+0.5f,(float)inputColumns);
2275         const unsigned int n = stop - start;
2276
2277         // calculate how many steps this workitem will contribute
2278         unsigned int numStepsPerWorkItem = n / numItems;
2279         numStepsPerWorkItem += ((numItems*numStepsPerWorkItem)==n?0:1);
2280
2281         const unsigned int startStep = (itemID%numItems)*numStepsPerWorkItem;
2282         if (startStep < n) {
2283           const unsigned int stopStep = min(startStep+numStepsPerWorkItem, n);
2284
2285           unsigned int cacheIndex = start+startStep-cacheRangeStartX;
2286           if (matte == 0) {
2287
2288             for (unsigned int i = startStep; i < stopStep; i++,cacheIndex++) {
2289               float4 cp = convert_float4(inputImageCache[cacheIndex]);
2290
2291               float weight = getResizeFilterWeight(resizeFilterCubicCoefficients,(ResizeWeightingFunctionType)resizeFilterType
2292                 , (ResizeWeightingFunctionType)resizeWindowType
2293                 , resizeFilterScale, resizeFilterWindowSupport, resizeFilterBlur,scale*(start+i-bisect+0.5));
2294
2295               filteredPixel += ((float4)weight)*cp;
2296               density+=weight;
2297             }
2298
2299
2300           }
2301           else {
2302             for (unsigned int i = startStep; i < stopStep; i++,cacheIndex++) {
2303               CLPixelType p = inputImageCache[cacheIndex];
2304
2305               float weight = getResizeFilterWeight(resizeFilterCubicCoefficients,(ResizeWeightingFunctionType)resizeFilterType
2306                 , (ResizeWeightingFunctionType)resizeWindowType
2307                 , resizeFilterScale, resizeFilterWindowSupport, resizeFilterBlur,scale*(start+i-bisect+0.5));
2308
2309               float alpha = weight * QuantumScale * GetPixelAlpha(p);
2310               float4 cp = convert_float4(p);
2311
2312               filteredPixel.x += alpha * cp.x;
2313               filteredPixel.y += alpha * cp.y;
2314               filteredPixel.z += alpha * cp.z;
2315               filteredPixel.w += weight * cp.w;
2316
2317               density+=weight;
2318               gamma+=alpha;
2319             }
2320          }
2321       }
2322     }
2323
2324     // initialize the accumulators to zero
2325     if (itemID < actualNumPixelInThisChunk) {
2326       outputPixelCache[itemID] = (float4)0.0f;
2327       densityCache[itemID] = 0.0f;
2328       if (matte != 0)
2329         gammaCache[itemID] = 0.0f;
2330     }
2331     barrier(CLK_LOCAL_MEM_FENCE);
2332
2333     // accumulatte the filtered pixel value and the density
2334     for (unsigned int i = 0; i < numItems; i++) {
2335       if (pixelIndex != -1) {
2336         if (itemID%numItems == i) {
2337           outputPixelCache[pixelIndex]+=filteredPixel;
2338           densityCache[pixelIndex]+=density;
2339           if (matte!=0) {
2340             gammaCache[pixelIndex]+=gamma;
2341           }
2342         }
2343       }
2344       barrier(CLK_LOCAL_MEM_FENCE);
2345     }
2346
2347     if (itemID < actualNumPixelInThisChunk) {
2348       if (matte==0) {
2349         float density = densityCache[itemID];
2350         float4 filteredPixel = outputPixelCache[itemID];
2351         if (density!= 0.0f && density != 1.0)
2352         {
2353           density = PerceptibleReciprocal(density);
2354           filteredPixel *= (float4)density;
2355         }
2356         filteredImage[y*filteredColumns+chunkStartX+itemID] = (CLPixelType) (ClampToQuantum(filteredPixel.x)
2357                                                                        , ClampToQuantum(filteredPixel.y)
2358                                                                        , ClampToQuantum(filteredPixel.z)
2359                                                                        , ClampToQuantum(filteredPixel.w));
2360       }
2361       else {
2362         float density = densityCache[itemID];
2363         float gamma = gammaCache[itemID];
2364         float4 filteredPixel = outputPixelCache[itemID];
2365
2366         if (density!= 0.0f && density != 1.0) {
2367           density = PerceptibleReciprocal(density);
2368           filteredPixel *= (float4)density;
2369           gamma *= density;
2370         }
2371         gamma = PerceptibleReciprocal(gamma);
2372
2373         CLPixelType fp;
2374         fp = (CLPixelType) ( ClampToQuantum(gamma*filteredPixel.x)
2375           , ClampToQuantum(gamma*filteredPixel.y)
2376           , ClampToQuantum(gamma*filteredPixel.z)
2377           , ClampToQuantum(filteredPixel.w));
2378
2379         filteredImage[y*filteredColumns+chunkStartX+itemID] = fp;
2380
2381       }
2382     }
2383
2384     } // end of chunking loop
2385   }
2386   )
2387
2388
2389
2390   STRINGIFY(
2391  __kernel __attribute__((reqd_work_group_size(256, 1, 1)))
2392  void ResizeHorizontalFilterSinc(const __global CLPixelType* inputImage, const unsigned int inputColumns, const unsigned int inputRows, const unsigned int matte
2393   , const float xFactor, __global CLPixelType* filteredImage, const unsigned int filteredColumns, const unsigned int filteredRows
2394   , const int resizeFilterType, const int resizeWindowType
2395   , const __global float* resizeFilterCubicCoefficients
2396   , const float resizeFilterScale, const float resizeFilterSupport, const float resizeFilterWindowSupport, const float resizeFilterBlur
2397   , __local CLPixelType* inputImageCache, const int numCachedPixels, const unsigned int pixelPerWorkgroup, const unsigned int pixelChunkSize
2398   , __local float4* outputPixelCache, __local float* densityCache, __local float* gammaCache) {
2399
2400     ResizeHorizontalFilter(inputImage,inputColumns,inputRows,matte
2401     ,xFactor, filteredImage, filteredColumns, filteredRows
2402     ,SincWeightingFunction, SincWeightingFunction
2403     ,resizeFilterCubicCoefficients
2404     ,resizeFilterScale, resizeFilterSupport, resizeFilterWindowSupport, resizeFilterBlur
2405     ,inputImageCache, numCachedPixels, pixelPerWorkgroup, pixelChunkSize
2406     ,outputPixelCache, densityCache, gammaCache);
2407
2408   }
2409   )
2410
2411
2412   STRINGIFY(
2413  __kernel __attribute__((reqd_work_group_size(1, 256, 1)))
2414  void ResizeVerticalFilter(const __global CLPixelType* inputImage, const unsigned int inputColumns, const unsigned int inputRows, const unsigned int matte
2415   , const float yFactor, __global CLPixelType* filteredImage, const unsigned int filteredColumns, const unsigned int filteredRows
2416   , const int resizeFilterType, const int resizeWindowType
2417   , const __global float* resizeFilterCubicCoefficients
2418   , const float resizeFilterScale, const float resizeFilterSupport, const float resizeFilterWindowSupport, const float resizeFilterBlur
2419   , __local CLPixelType* inputImageCache, const int numCachedPixels, const unsigned int pixelPerWorkgroup, const unsigned int pixelChunkSize
2420   , __local float4* outputPixelCache, __local float* densityCache, __local float* gammaCache) {
2421
2422
2423     // calculate the range of resized image pixels computed by this workgroup
2424     const unsigned int startY = get_group_id(1)*pixelPerWorkgroup;
2425     const unsigned int stopY = min(startY + pixelPerWorkgroup,filteredRows);
2426     const unsigned int actualNumPixelToCompute = stopY - startY;
2427
2428     // calculate the range of input image pixels to cache
2429     float scale = max(1.0f/yFactor+MagickEpsilon ,1.0f);
2430     const float support = max(scale*resizeFilterSupport,0.5f);
2431     scale = PerceptibleReciprocal(scale);
2432
2433     const int cacheRangeStartY = max((int)((startY+0.5f)/yFactor+MagickEpsilon-support+0.5f),(int)(0));
2434     const int cacheRangeEndY = min((int)(cacheRangeStartY + numCachedPixels), (int)inputRows);
2435
2436     // cache the input pixels into local memory
2437     const unsigned int x = get_global_id(0);
2438     event_t e = async_work_group_strided_copy(inputImageCache, inputImage+cacheRangeStartY*inputColumns+x, cacheRangeEndY-cacheRangeStartY, inputColumns, 0);
2439     wait_group_events(1,&e);
2440
2441     unsigned int totalNumChunks = (actualNumPixelToCompute+pixelChunkSize-1)/pixelChunkSize;
2442     for (unsigned int chunk = 0; chunk < totalNumChunks; chunk++)
2443     {
2444
2445       const unsigned int chunkStartY = startY + chunk*pixelChunkSize;
2446       const unsigned int chunkStopY = min(chunkStartY + pixelChunkSize, stopY);
2447       const unsigned int actualNumPixelInThisChunk = chunkStopY - chunkStartY;
2448
2449       // determine which resized pixel computed by this workitem
2450       const unsigned int itemID = get_local_id(1);
2451       const unsigned int numItems = getNumWorkItemsPerPixel(actualNumPixelInThisChunk, get_local_size(1));
2452
2453       const int pixelIndex = pixelToCompute(itemID, actualNumPixelInThisChunk, get_local_size(1));
2454
2455       float4 filteredPixel = (float4)0.0f;
2456       float density = 0.0f;
2457       float gamma = 0.0f;
2458       // -1 means this workitem doesn't participate in the computation
2459       if (pixelIndex != -1) {
2460
2461         // x coordinated of the resized pixel computed by this workitem
2462         const int y = chunkStartY + pixelIndex;
2463
2464         // calculate how many steps required for this pixel
2465         const float bisect = (y+0.5)/yFactor+MagickEpsilon;
2466         const unsigned int start = (unsigned int)max(bisect-support+0.5f,0.0f);
2467         const unsigned int stop  = (unsigned int)min(bisect+support+0.5f,(float)inputRows);
2468         const unsigned int n = stop - start;
2469
2470         // calculate how many steps this workitem will contribute
2471         unsigned int numStepsPerWorkItem = n / numItems;
2472         numStepsPerWorkItem += ((numItems*numStepsPerWorkItem)==n?0:1);
2473
2474         const unsigned int startStep = (itemID%numItems)*numStepsPerWorkItem;
2475         if (startStep < n) {
2476           const unsigned int stopStep = min(startStep+numStepsPerWorkItem, n);
2477
2478           unsigned int cacheIndex = start+startStep-cacheRangeStartY;
2479           if (matte == 0) {
2480
2481             for (unsigned int i = startStep; i < stopStep; i++,cacheIndex++) {
2482               float4 cp = convert_float4(inputImageCache[cacheIndex]);
2483
2484               float weight = getResizeFilterWeight(resizeFilterCubicCoefficients,(ResizeWeightingFunctionType)resizeFilterType
2485                 , (ResizeWeightingFunctionType)resizeWindowType
2486                 , resizeFilterScale, resizeFilterWindowSupport, resizeFilterBlur,scale*(start+i-bisect+0.5));
2487
2488               filteredPixel += ((float4)weight)*cp;
2489               density+=weight;
2490             }
2491
2492
2493           }
2494           else {
2495             for (unsigned int i = startStep; i < stopStep; i++,cacheIndex++) {
2496               CLPixelType p = inputImageCache[cacheIndex];
2497
2498               float weight = getResizeFilterWeight(resizeFilterCubicCoefficients,(ResizeWeightingFunctionType)resizeFilterType
2499                 , (ResizeWeightingFunctionType)resizeWindowType
2500                 , resizeFilterScale, resizeFilterWindowSupport, resizeFilterBlur,scale*(start+i-bisect+0.5));
2501
2502               float alpha = weight * QuantumScale * GetPixelAlpha(p);
2503               float4 cp = convert_float4(p);
2504
2505               filteredPixel.x += alpha * cp.x;
2506               filteredPixel.y += alpha * cp.y;
2507               filteredPixel.z += alpha * cp.z;
2508               filteredPixel.w += weight * cp.w;
2509
2510               density+=weight;
2511               gamma+=alpha;
2512             }
2513          }
2514       }
2515     }
2516
2517     // initialize the accumulators to zero
2518     if (itemID < actualNumPixelInThisChunk) {
2519       outputPixelCache[itemID] = (float4)0.0f;
2520       densityCache[itemID] = 0.0f;
2521       if (matte != 0)
2522         gammaCache[itemID] = 0.0f;
2523     }
2524     barrier(CLK_LOCAL_MEM_FENCE);
2525
2526     // accumulatte the filtered pixel value and the density
2527     for (unsigned int i = 0; i < numItems; i++) {
2528       if (pixelIndex != -1) {
2529         if (itemID%numItems == i) {
2530           outputPixelCache[pixelIndex]+=filteredPixel;
2531           densityCache[pixelIndex]+=density;
2532           if (matte!=0) {
2533             gammaCache[pixelIndex]+=gamma;
2534           }
2535         }
2536       }
2537       barrier(CLK_LOCAL_MEM_FENCE);
2538     }
2539
2540     if (itemID < actualNumPixelInThisChunk) {
2541       if (matte==0) {
2542         float density = densityCache[itemID];
2543         float4 filteredPixel = outputPixelCache[itemID];
2544         if (density!= 0.0f && density != 1.0)
2545         {
2546           density = PerceptibleReciprocal(density);
2547           filteredPixel *= (float4)density;
2548         }
2549         filteredImage[(chunkStartY+itemID)*filteredColumns+x] = (CLPixelType) (ClampToQuantum(filteredPixel.x)
2550                                                                        , ClampToQuantum(filteredPixel.y)
2551                                                                        , ClampToQuantum(filteredPixel.z)
2552                                                                        , ClampToQuantum(filteredPixel.w));
2553       }
2554       else {
2555         float density = densityCache[itemID];
2556         float gamma = gammaCache[itemID];
2557         float4 filteredPixel = outputPixelCache[itemID];
2558
2559         if (density!= 0.0f && density != 1.0) {
2560           density = PerceptibleReciprocal(density);
2561           filteredPixel *= (float4)density;
2562           gamma *= density;
2563         }
2564         gamma = PerceptibleReciprocal(gamma);
2565
2566         CLPixelType fp;
2567         fp = (CLPixelType) ( ClampToQuantum(gamma*filteredPixel.x)
2568           , ClampToQuantum(gamma*filteredPixel.y)
2569           , ClampToQuantum(gamma*filteredPixel.z)
2570           , ClampToQuantum(filteredPixel.w));
2571
2572         filteredImage[(chunkStartY+itemID)*filteredColumns+x] = fp;
2573
2574       }
2575     }
2576
2577     } // end of chunking loop
2578   }
2579   )
2580
2581
2582
2583   STRINGIFY(
2584  __kernel __attribute__((reqd_work_group_size(1, 256, 1)))
2585  void ResizeVerticalFilterSinc(const __global CLPixelType* inputImage, const unsigned int inputColumns, const unsigned int inputRows, const unsigned int matte
2586   , const float yFactor, __global CLPixelType* filteredImage, const unsigned int filteredColumns, const unsigned int filteredRows
2587   , const int resizeFilterType, const int resizeWindowType
2588   , const __global float* resizeFilterCubicCoefficients
2589   , const float resizeFilterScale, const float resizeFilterSupport, const float resizeFilterWindowSupport, const float resizeFilterBlur
2590   , __local CLPixelType* inputImageCache, const int numCachedPixels, const unsigned int pixelPerWorkgroup, const unsigned int pixelChunkSize
2591   , __local float4* outputPixelCache, __local float* densityCache, __local float* gammaCache) {
2592     ResizeVerticalFilter(inputImage,inputColumns,inputRows,matte
2593       ,yFactor,filteredImage,filteredColumns,filteredRows
2594       ,SincWeightingFunction, SincWeightingFunction
2595       ,resizeFilterCubicCoefficients
2596       ,resizeFilterScale,resizeFilterSupport,resizeFilterWindowSupport,resizeFilterBlur
2597       ,inputImageCache,numCachedPixels,pixelPerWorkgroup,pixelChunkSize
2598       ,outputPixelCache,densityCache,gammaCache);
2599   }
2600   )
2601
2602   STRINGIFY(
2603
2604   inline float GetPseudoRandomValue(uint4* seed, const float normalizeRand) {
2605     uint4 s = *seed;
2606     do {
2607       unsigned int alpha = (unsigned int) (s.y ^ (s.y << 11));
2608       s.y=s.z;
2609       s.z=s.w;
2610       s.w=s.x;
2611       s.x = (s.x ^ (s.x >> 19)) ^ (alpha ^ (alpha >> 8));
2612     } while (s.x == ~0UL);
2613     *seed = s;
2614     return (normalizeRand*s.x);
2615   }
2616
2617   __kernel void randomNumberGeneratorKernel(__global uint* seeds, const float normalizeRand
2618                                            , __global float* randomNumbers, const uint init
2619                                            ,const uint numRandomNumbers) {
2620
2621     unsigned int id = get_global_id(0);
2622     unsigned int seed[4];
2623
2624     if (init!=0) {
2625       seed[0] = seeds[id*4];
2626       seed[1] = 0x50a7f451;
2627       seed[2] = 0x5365417e;
2628       seed[3] = 0xc3a4171a;
2629     }
2630     else {
2631       seed[0] = seeds[id*4];
2632       seed[1] = seeds[id*4+1];
2633       seed[2] = seeds[id*4+2];
2634       seed[3] = seeds[id*4+3];
2635     }
2636
2637     unsigned int numRandomNumbersPerItem = (numRandomNumbers+get_global_size(0)-1)/get_global_size(0);
2638     for (unsigned int i = 0; i < numRandomNumbersPerItem; i++) {
2639       do
2640       {
2641         unsigned int alpha=(unsigned int) (seed[1] ^ (seed[1] << 11));
2642         seed[1]=seed[2];
2643         seed[2]=seed[3];
2644         seed[3]=seed[0];
2645         seed[0]=(seed[0] ^ (seed[0] >> 19)) ^ (alpha ^ (alpha >> 8));
2646       } while (seed[0] == ~0UL);
2647       unsigned int pos = (get_group_id(0)*get_local_size(0)*numRandomNumbersPerItem)
2648                           + get_local_size(0) * i + get_local_id(0);
2649
2650       if (pos >= numRandomNumbers)
2651         break;
2652       randomNumbers[pos] = normalizeRand*seed[0];
2653     }
2654
2655     /* save the seeds for the time*/
2656     seeds[id*4]   = seed[0];
2657     seeds[id*4+1] = seed[1];
2658     seeds[id*4+2] = seed[2];
2659     seeds[id*4+3] = seed[3];
2660   }
2661
2662   )
2663
2664
2665   STRINGIFY(
2666
2667   typedef enum
2668   {
2669     UndefinedNoise,
2670     UniformNoise,
2671     GaussianNoise,
2672     MultiplicativeGaussianNoise,
2673     ImpulseNoise,
2674     LaplacianNoise,
2675     PoissonNoise,
2676     RandomNoise
2677   } NoiseType;
2678
2679   typedef struct {
2680     const global float* rns;
2681   } RandomNumbers;
2682
2683
2684   float ReadPseudoRandomValue(RandomNumbers* r) {
2685     float v = *r->rns;
2686     r->rns++;
2687     return v;
2688   }
2689   )
2690
2691   OPENCL_DEFINE(SigmaUniform, (attenuate*0.015625f))
2692   OPENCL_DEFINE(SigmaGaussian,(attenuate*0.015625f))
2693   OPENCL_DEFINE(SigmaImpulse,  (attenuate*0.1f))
2694   OPENCL_DEFINE(SigmaLaplacian, (attenuate*0.0390625f))
2695   OPENCL_DEFINE(SigmaMultiplicativeGaussian,  (attenuate*0.5f))
2696   OPENCL_DEFINE(SigmaPoisson,  (attenuate*12.5f))
2697   OPENCL_DEFINE(SigmaRandom,  (attenuate))
2698   OPENCL_DEFINE(TauGaussian,  (attenuate*0.078125f))
2699
2700   STRINGIFY(
2701   float GenerateDifferentialNoise(RandomNumbers* r, CLQuantum pixel, NoiseType noise_type, float attenuate) {
2702
2703     float
2704       alpha,
2705       beta,
2706       noise,
2707       sigma;
2708
2709     noise = 0.0f;
2710     alpha=ReadPseudoRandomValue(r);
2711     switch(noise_type) {
2712     case UniformNoise:
2713     default:
2714       {
2715         noise=(pixel+QuantumRange*SigmaUniform*(alpha-0.5f));
2716         break;
2717       }
2718     case GaussianNoise:
2719       {
2720         float
2721           gamma,
2722           tau;
2723
2724         if (alpha == 0.0f)
2725           alpha=1.0f;
2726         beta=ReadPseudoRandomValue(r);
2727         gamma=sqrt(-2.0f*log(alpha));
2728         sigma=gamma*cospi((2.0f*beta));
2729         tau=gamma*sinpi((2.0f*beta));
2730         noise=(float)(pixel+sqrt((float) pixel)*SigmaGaussian*sigma+
2731                       QuantumRange*TauGaussian*tau);
2732         break;
2733       }
2734
2735
2736     case ImpulseNoise:
2737     {
2738       if (alpha < (SigmaImpulse/2.0f))
2739         noise=0.0f;
2740       else
2741         if (alpha >= (1.0f-(SigmaImpulse/2.0f)))
2742           noise=(float)QuantumRange;
2743         else
2744           noise=(float)pixel;
2745       break;
2746     }
2747     case LaplacianNoise:
2748     {
2749       if (alpha <= 0.5f)
2750         {
2751           if (alpha <= MagickEpsilon)
2752             noise=(float) (pixel-QuantumRange);
2753           else
2754             noise=(float) (pixel+QuantumRange*SigmaLaplacian*log(2.0f*alpha)+
2755               0.5f);
2756           break;
2757         }
2758       beta=1.0f-alpha;
2759       if (beta <= (0.5f*MagickEpsilon))
2760         noise=(float) (pixel+QuantumRange);
2761       else
2762         noise=(float) (pixel-QuantumRange*SigmaLaplacian*log(2.0f*beta)+0.5f);
2763       break;
2764     }
2765     case MultiplicativeGaussianNoise:
2766     {
2767       sigma=1.0f;
2768       if (alpha > MagickEpsilon)
2769         sigma=sqrt(-2.0f*log(alpha));
2770       beta=ReadPseudoRandomValue(r);
2771       noise=(float) (pixel+pixel*SigmaMultiplicativeGaussian*sigma*
2772         cospi((float) (2.0f*beta))/2.0f);
2773       break;
2774     }
2775     case PoissonNoise:
2776     {
2777       float
2778         poisson;
2779       unsigned int i;
2780       poisson=exp(-SigmaPoisson*QuantumScale*pixel);
2781       for (i=0; alpha > poisson; i++)
2782       {
2783         beta=ReadPseudoRandomValue(r);
2784         alpha*=beta;
2785       }
2786       noise=(float) (QuantumRange*i/SigmaPoisson);
2787       break;
2788     }
2789     case RandomNoise:
2790     {
2791       noise=(float) (QuantumRange*SigmaRandom*alpha);
2792       break;
2793     }
2794
2795     };
2796     return noise;
2797   }
2798
2799   __kernel
2800   void AddNoiseImage(const __global CLPixelType* inputImage, __global CLPixelType* filteredImage
2801                     ,const unsigned int inputColumns, const unsigned int inputRows
2802                     ,const ChannelType channel
2803                     ,const NoiseType noise_type, const float attenuate
2804                     ,const __global float* randomNumbers, const unsigned int numRandomNumbersPerPixel
2805                     ,const unsigned int rowOffset) {
2806
2807     unsigned int x = get_global_id(0);
2808     unsigned int y = get_global_id(1) + rowOffset;
2809     RandomNumbers r;
2810     r.rns = randomNumbers + (get_global_id(1) * inputColumns + get_global_id(0))*numRandomNumbersPerPixel;
2811
2812     CLPixelType p = inputImage[y*inputColumns+x];
2813     CLPixelType q = filteredImage[y*inputColumns+x];
2814
2815     if ((channel&RedChannel)!=0) {
2816       setRed(&q,ClampToQuantum(GenerateDifferentialNoise(&r,getRed(p),noise_type,attenuate)));
2817     }
2818
2819     if ((channel&GreenChannel)!=0) {
2820       setGreen(&q,ClampToQuantum(GenerateDifferentialNoise(&r,getGreen(p),noise_type,attenuate)));
2821     }
2822
2823     if ((channel&BlueChannel)!=0) {
2824       setBlue(&q,ClampToQuantum(GenerateDifferentialNoise(&r,getBlue(p),noise_type,attenuate)));
2825     }
2826
2827     if ((channel & OpacityChannel) != 0) {
2828       setOpacity(&q,ClampToQuantum(GenerateDifferentialNoise(&r,getOpacity(p),noise_type,attenuate)));
2829     }
2830
2831     filteredImage[y*inputColumns+x] = q;
2832   }
2833
2834   )
2835
2836   STRINGIFY(
2837   __kernel
2838   void RandomImage(__global CLPixelType* inputImage,
2839                    const uint imageColumns, const uint imageRows,
2840                    __global uint* seeds,
2841                    const float randNormNumerator,
2842                    const uint randNormDenominator) {
2843
2844     unsigned int numGenerators = get_global_size(0);
2845     unsigned numRandPixelsPerWorkItem = ((imageColumns*imageRows) + (numGenerators-1))
2846                                         / numGenerators;
2847
2848     uint4 s;
2849     s.x = seeds[get_global_id(0)*4];
2850     s.y = seeds[get_global_id(0)*4+1];
2851     s.z = seeds[get_global_id(0)*4+2];
2852     s.w = seeds[get_global_id(0)*4+3];
2853
2854     unsigned int offset = get_group_id(0) * get_local_size(0) * numRandPixelsPerWorkItem;
2855     for (unsigned int n = 0; n < numRandPixelsPerWorkItem; n++)
2856     {
2857       int i = offset + n*get_local_size(0) + get_local_id(0);
2858       if (i >= imageColumns*imageRows)
2859         break;
2860
2861       float rand = GetPseudoRandomValue(&s,randNormNumerator/randNormDenominator);
2862       CLQuantum v = ClampToQuantum(QuantumRange*rand);
2863
2864       CLPixelType p;
2865       setRed(&p,v);
2866       setGreen(&p,v);
2867       setBlue(&p,v);
2868       setOpacity(&p,0);
2869
2870       inputImage[i] = p;
2871     }
2872
2873     seeds[get_global_id(0)*4]   = s.x;
2874     seeds[get_global_id(0)*4+1] = s.y;
2875     seeds[get_global_id(0)*4+2] = s.z;
2876     seeds[get_global_id(0)*4+3] = s.w;
2877   }
2878   )
2879
2880   STRINGIFY(
2881     __kernel
2882     void MotionBlur(const __global CLPixelType *input, __global CLPixelType *output,
2883                     const unsigned int imageWidth, const unsigned int imageHeight,
2884                     const __global float *filter, const unsigned int width, const __global int2* offset,
2885                     const float4 bias,
2886                     const ChannelType channel, const unsigned int matte) {
2887
2888       int2 currentPixel;
2889       currentPixel.x = get_global_id(0);
2890       currentPixel.y = get_global_id(1);
2891
2892       if (currentPixel.x >= imageWidth
2893           || currentPixel.y >= imageHeight)
2894           return;
2895
2896       float4 pixel;
2897       pixel.x = (float)bias.x;
2898       pixel.y = (float)bias.y;
2899       pixel.z = (float)bias.z;
2900       pixel.w = (float)bias.w;
2901
2902       if (((channel & OpacityChannel) == 0) || (matte == 0)) {
2903
2904         for (int i = 0; i < width; i++) {
2905           // only support EdgeVirtualPixelMethod through ClampToCanvas
2906           // TODO: implement other virtual pixel method
2907           int2 samplePixel = currentPixel + offset[i];
2908           samplePixel.x = ClampToCanvas(samplePixel.x, imageWidth);
2909           samplePixel.y = ClampToCanvas(samplePixel.y, imageHeight);
2910           CLPixelType samplePixelValue = input[ samplePixel.y * imageWidth + samplePixel.x];
2911
2912           pixel.x += (filter[i] * (float)samplePixelValue.x);
2913           pixel.y += (filter[i] * (float)samplePixelValue.y);
2914           pixel.z += (filter[i] * (float)samplePixelValue.z);
2915           pixel.w += (filter[i] * (float)samplePixelValue.w);
2916         }
2917
2918         CLPixelType outputPixel;
2919         outputPixel.x = ClampToQuantum(pixel.x);
2920         outputPixel.y = ClampToQuantum(pixel.y);
2921         outputPixel.z = ClampToQuantum(pixel.z);
2922         outputPixel.w = ClampToQuantum(pixel.w);
2923         output[currentPixel.y * imageWidth + currentPixel.x] = outputPixel;
2924       }
2925       else {
2926
2927         float gamma = 0.0f;
2928         for (int i = 0; i < width; i++) {
2929           // only support EdgeVirtualPixelMethod through ClampToCanvas
2930           // TODO: implement other virtual pixel method
2931           int2 samplePixel = currentPixel + offset[i];
2932           samplePixel.x = ClampToCanvas(samplePixel.x, imageWidth);
2933           samplePixel.y = ClampToCanvas(samplePixel.y, imageHeight);
2934
2935           CLPixelType samplePixelValue = input[ samplePixel.y * imageWidth + samplePixel.x];
2936
2937           float alpha = QuantumScale*(QuantumRange-samplePixelValue.w);
2938           float k = filter[i];
2939           pixel.x = pixel.x + k * alpha * samplePixelValue.x;
2940           pixel.y = pixel.y + k * alpha * samplePixelValue.y;
2941           pixel.z = pixel.z + k * alpha * samplePixelValue.z;
2942
2943           pixel.w += k * alpha * samplePixelValue.w;
2944
2945           gamma+=k*alpha;
2946         }
2947         gamma = PerceptibleReciprocal(gamma);
2948         pixel.xyz = gamma*pixel.xyz;
2949
2950         CLPixelType outputPixel;
2951         outputPixel.x = ClampToQuantum(pixel.x);
2952         outputPixel.y = ClampToQuantum(pixel.y);
2953         outputPixel.z = ClampToQuantum(pixel.z);
2954         outputPixel.w = ClampToQuantum(pixel.w);
2955         output[currentPixel.y * imageWidth + currentPixel.x] = outputPixel;
2956       }
2957     }
2958   )
2959
2960   STRINGIFY(
2961     typedef enum
2962     {
2963       UndefinedCompositeOp,
2964       NoCompositeOp,
2965       ModulusAddCompositeOp,
2966       AtopCompositeOp,
2967       BlendCompositeOp,
2968       BumpmapCompositeOp,
2969       ChangeMaskCompositeOp,
2970       ClearCompositeOp,
2971       ColorBurnCompositeOp,
2972       ColorDodgeCompositeOp,
2973       ColorizeCompositeOp,
2974       CopyBlackCompositeOp,
2975       CopyBlueCompositeOp,
2976       CopyCompositeOp,
2977       CopyCyanCompositeOp,
2978       CopyGreenCompositeOp,
2979       CopyMagentaCompositeOp,
2980       CopyOpacityCompositeOp,
2981       CopyRedCompositeOp,
2982       CopyYellowCompositeOp,
2983       DarkenCompositeOp,
2984       DstAtopCompositeOp,
2985       DstCompositeOp,
2986       DstInCompositeOp,
2987       DstOutCompositeOp,
2988       DstOverCompositeOp,
2989       DifferenceCompositeOp,
2990       DisplaceCompositeOp,
2991       DissolveCompositeOp,
2992       ExclusionCompositeOp,
2993       HardLightCompositeOp,
2994       HueCompositeOp,
2995       InCompositeOp,
2996       LightenCompositeOp,
2997       LinearLightCompositeOp,
2998       LuminizeCompositeOp,
2999       MinusDstCompositeOp,
3000       ModulateCompositeOp,
3001       MultiplyCompositeOp,
3002       OutCompositeOp,
3003       OverCompositeOp,
3004       OverlayCompositeOp,
3005       PlusCompositeOp,
3006       ReplaceCompositeOp,
3007       SaturateCompositeOp,
3008       ScreenCompositeOp,
3009       SoftLightCompositeOp,
3010       SrcAtopCompositeOp,
3011       SrcCompositeOp,
3012       SrcInCompositeOp,
3013       SrcOutCompositeOp,
3014       SrcOverCompositeOp,
3015       ModulusSubtractCompositeOp,
3016       ThresholdCompositeOp,
3017       XorCompositeOp,
3018       /* These are new operators, added after the above was last sorted.
3019        * The list should be re-sorted only when a new library version is
3020        * created.
3021        */
3022       DivideDstCompositeOp,
3023       DistortCompositeOp,
3024       BlurCompositeOp,
3025       PegtopLightCompositeOp,
3026       VividLightCompositeOp,
3027       PinLightCompositeOp,
3028       LinearDodgeCompositeOp,
3029       LinearBurnCompositeOp,
3030       MathematicsCompositeOp,
3031       DivideSrcCompositeOp,
3032       MinusSrcCompositeOp,
3033       DarkenIntensityCompositeOp,
3034       LightenIntensityCompositeOp
3035     } CompositeOperator;
3036   )
3037
3038   STRINGIFY(
3039     inline float ColorDodge(const float Sca,
3040       const float Sa,const float Dca,const float Da)
3041     {
3042       /*
3043         Oct 2004 SVG specification.
3044       */
3045       if ((Sca*Da+Dca*Sa) >= Sa*Da)
3046         return(Sa*Da+Sca*(1.0-Da)+Dca*(1.0-Sa));
3047       return(Dca*Sa*Sa/(Sa-Sca)+Sca*(1.0-Da)+Dca*(1.0-Sa));
3048
3049
3050       /*
3051         New specification, March 2009 SVG specification.  This specification was
3052         also wrong of non-overlap cases.
3053       */
3054       /*
3055       if ((fabs(Sca-Sa) < MagickEpsilon) && (fabs(Dca) < MagickEpsilon))
3056         return(Sca*(1.0-Da));
3057       if (fabs(Sca-Sa) < MagickEpsilon)
3058         return(Sa*Da+Sca*(1.0-Da)+Dca*(1.0-Sa));
3059       return(Sa*MagickMin(Da,Dca*Sa/(Sa-Sca)));
3060       */
3061
3062       /*
3063         Working from first principles using the original formula:
3064
3065            f(Sc,Dc) = Dc/(1-Sc)
3066
3067         This works correctly! Looks like the 2004 model was right but just
3068         required a extra condition for correct handling.
3069       */
3070
3071       /*
3072       if ((fabs(Sca-Sa) < MagickEpsilon) && (fabs(Dca) < MagickEpsilon))
3073         return(Sca*(1.0-Da)+Dca*(1.0-Sa));
3074       if (fabs(Sca-Sa) < MagickEpsilon)
3075         return(Sa*Da+Sca*(1.0-Da)+Dca*(1.0-Sa));
3076       return(Dca*Sa*Sa/(Sa-Sca)+Sca*(1.0-Da)+Dca*(1.0-Sa));
3077       */
3078     }
3079
3080     inline void CompositeColorDodge(const float4 *p,
3081       const float4 *q,float4 *composite) {
3082
3083       float
3084       Da,
3085       gamma,
3086       Sa;
3087
3088       Sa=1.0f-QuantumScale*getOpacityF4(*p);  /* simplify and speed up equations */
3089       Da=1.0f-QuantumScale*getOpacityF4(*q);
3090       gamma=RoundToUnity(Sa+Da-Sa*Da); /* over blend, as per SVG doc */
3091       setOpacityF4(composite, QuantumRange*(1.0-gamma));
3092       gamma=QuantumRange/(fabs(gamma) < MagickEpsilon ? MagickEpsilon : gamma);
3093       setRedF4(composite,gamma*ColorDodge(QuantumScale*getRedF4(*p)*Sa,Sa,QuantumScale*
3094         getRedF4(*q)*Da,Da));
3095       setGreenF4(composite,gamma*ColorDodge(QuantumScale*getGreenF4(*p)*Sa,Sa,QuantumScale*
3096         getGreenF4(*q)*Da,Da));
3097       setBlueF4(composite,gamma*ColorDodge(QuantumScale*getBlueF4(*p)*Sa,Sa,QuantumScale*
3098         getBlueF4(*q)*Da,Da));
3099     }
3100   )
3101
3102   STRINGIFY(
3103     inline void MagickPixelCompositePlus(const float4 *p,
3104       const float alpha,const float4 *q,
3105       const float beta,float4 *composite)
3106     {
3107       float
3108         gamma;
3109
3110       float
3111         Da,
3112         Sa;
3113       /*
3114         Add two pixels with the given opacities.
3115       */
3116       Sa=1.0-QuantumScale*alpha;
3117       Da=1.0-QuantumScale*beta;
3118       gamma=RoundToUnity(Sa+Da);  /* 'Plus' blending -- not 'Over' blending */
3119       setOpacityF4(composite,(float) QuantumRange*(1.0-gamma));
3120       gamma=PerceptibleReciprocal(gamma);
3121       setRedF4(composite,gamma*(Sa*getRedF4(*p)+Da*getRedF4(*q)));
3122       setGreenF4(composite,gamma*(Sa*getGreenF4(*p)+Da*getGreenF4(*q)));
3123       setBlueF4(composite,gamma*(Sa*getBlueF4(*p)+Da*getBlueF4(*q)));
3124     }
3125   )
3126
3127   STRINGIFY(
3128     inline void MagickPixelCompositeBlend(const float4 *p,
3129       const float alpha,const float4 *q,
3130       const float beta,float4 *composite)
3131     {
3132       MagickPixelCompositePlus(p,(float) (QuantumRange-alpha*
3133       (QuantumRange-getOpacityF4(*p))),q,(float) (QuantumRange-beta*
3134       (QuantumRange-getOpacityF4(*q))),composite);
3135     }
3136   )
3137
3138   STRINGIFY(
3139     __kernel
3140     void Composite(__global CLPixelType *image,
3141                    const unsigned int imageWidth,
3142                    const unsigned int imageHeight,
3143                    const __global CLPixelType *compositeImage,
3144                    const unsigned int compositeWidth,
3145                    const unsigned int compositeHeight,
3146                    const unsigned int compose,
3147                    const ChannelType channel,
3148                    const unsigned int matte,
3149                    const float destination_dissolve,
3150                    const float source_dissolve) {
3151
3152       uint2 index;
3153       index.x = get_global_id(0);
3154       index.y = get_global_id(1);
3155
3156
3157       if (index.x >= imageWidth
3158         || index.y >= imageHeight) {
3159           return;
3160       }
3161       const CLPixelType inputPixel = image[index.y*imageWidth+index.x];
3162       float4 destination;
3163       setRedF4(&destination,getRed(inputPixel));
3164       setGreenF4(&destination,getGreen(inputPixel));
3165       setBlueF4(&destination,getBlue(inputPixel));
3166
3167
3168       const CLPixelType compositePixel
3169         = compositeImage[index.y*imageWidth+index.x];
3170       float4 source;
3171       setRedF4(&source,getRed(compositePixel));
3172       setGreenF4(&source,getGreen(compositePixel));
3173       setBlueF4(&source,getBlue(compositePixel));
3174
3175       if (matte != 0) {
3176         setOpacityF4(&destination,getOpacity(inputPixel));
3177         setOpacityF4(&source,getOpacity(compositePixel));
3178       }
3179       else {
3180         setOpacityF4(&destination,0.0f);
3181         setOpacityF4(&source,0.0f);
3182       }
3183
3184       float4 composite=destination;
3185
3186       CompositeOperator op = (CompositeOperator)compose;
3187       switch (op) {
3188       case ColorDodgeCompositeOp:
3189         CompositeColorDodge(&source,&destination,&composite);
3190         break;
3191       case BlendCompositeOp:
3192         MagickPixelCompositeBlend(&source,source_dissolve,&destination,
3193             destination_dissolve,&composite);
3194         break;
3195       default:
3196         // unsupported operators
3197         break;
3198       };
3199
3200       CLPixelType outputPixel;
3201       setRed(&outputPixel, ClampToQuantum(getRedF4(composite)));
3202       setGreen(&outputPixel, ClampToQuantum(getGreenF4(composite)));
3203       setBlue(&outputPixel, ClampToQuantum(getBlueF4(composite)));
3204       setOpacity(&outputPixel, ClampToQuantum(getOpacityF4(composite)));
3205       image[index.y*imageWidth+index.x] = outputPixel;
3206     }
3207   )
3208
3209   ;
3210
3211 #endif // MAGICKCORE_OPENCL_SUPPORT
3212
3213 #if defined(__cplusplus) || defined(c_plusplus)
3214 }
3215 #endif
3216
3217 #endif // _MAGICKCORE_ACCELERATE_PRIVATE_H