)
- STRINGIFY(\r
- __kernel void UnsharpMask(__global CLPixelType *im, __global CLPixelType *filtered_im,\r
- __constant float *filter,\r
- const unsigned int width, \r
- const unsigned int imageColumns, const unsigned int imageRows,\r
- __local float4 *pixels, \r
- const float gain, const float threshold, const unsigned int justBlur)\r
- {\r
- const int x = get_global_id(0);\r
- const int y = get_global_id(1);\r
-\r
- const unsigned int radius = (width - 1) / 2;\r
- \r
- int row = y - radius;\r
- int baseRow = get_group_id(1) * get_local_size(1) - radius;\r
- int endRow = (get_group_id(1) + 1) * get_local_size(1) + radius;\r
- \r
- while (row < endRow) {\r
- int srcy = (row < 0) ? -row : row; // mirror pad\r
- srcy = (srcy >= imageRows) ? (2 * imageRows - srcy - 1) : srcy;\r
- \r
- float4 value = 0.0f;\r
- \r
- int ix = x - radius;\r
- int i = 0;\r
-\r
- while (i + 7 < width) {\r
- for (int j = 0; j < 8; ++j) { // unrolled\r
- int srcx = ix + j;\r
- srcx = (srcx < 0) ? -srcx : srcx;\r
- srcx = (srcx >= imageColumns) ? (2 * imageColumns - srcx - 1) : srcx;\r
- value += filter[i + j] * convert_float4(im[srcx + srcy * imageColumns]);\r
- }\r
- ix += 8;\r
- i += 8;\r
- }\r
-\r
- while (i < width) {\r
- int srcx = (ix < 0) ? -ix : ix; // mirror pad\r
- srcx = (srcx >= imageColumns) ? (2 * imageColumns - srcx - 1) : srcx;\r
- value += filter[i] * convert_float4(im[srcx + srcy * imageColumns]);\r
- ++i;\r
- ++ix;\r
- } \r
- pixels[(row - baseRow) * get_local_size(0) + get_local_id(0)] = value;\r
- row += get_local_size(1);\r
- }\r
- \r
- \r
- barrier(CLK_LOCAL_MEM_FENCE);\r
-\r
- \r
- const int px = get_local_id(0);\r
- const int py = get_local_id(1);\r
- const int prp = get_local_size(0);\r
- float4 value = (float4)(0.0f);\r
- \r
- int i = 0;\r
- while (i + 7 < width) { // unrolled\r
- value += (float4)(filter[i]) * pixels[px + (py + i) * prp];\r
- value += (float4)(filter[i]) * pixels[px + (py + i + 1) * prp];\r
- value += (float4)(filter[i]) * pixels[px + (py + i + 2) * prp];\r
- value += (float4)(filter[i]) * pixels[px + (py + i + 3) * prp];\r
- value += (float4)(filter[i]) * pixels[px + (py + i + 4) * prp];\r
- value += (float4)(filter[i]) * pixels[px + (py + i + 5) * prp];\r
- value += (float4)(filter[i]) * pixels[px + (py + i + 6) * prp];\r
- value += (float4)(filter[i]) * pixels[px + (py + i + 7) * prp];\r
- i += 8;\r
- }\r
- while (i < width) {\r
- value += (float4)(filter[i]) * pixels[px + (py + i) * prp];\r
- ++i;\r
- }\r
-\r
- if (justBlur == 0) { // apply sharpening\r
- float4 srcPixel = convert_float4(im[x + y * imageColumns]);\r
- float4 diff = srcPixel - value;\r
-\r
- float quantumThreshold = QuantumRange*threshold;\r
-\r
- int4 mask = isless(fabs(2.0f * diff), (float4)quantumThreshold);\r
- value = select(srcPixel + diff * gain, srcPixel, mask);\r
- }\r
- \r
- if ((x < imageColumns) && (y < imageRows))\r
- filtered_im[x + y * imageColumns] = (CLPixelType)(ClampToQuantum(value.s0), ClampToQuantum(value.s1), ClampToQuantum(value.s2), ClampToQuantum(value.s3));\r
- } \r
+ STRINGIFY(
+ __kernel void UnsharpMask(__global CLPixelType *im, __global CLPixelType *filtered_im,
+ __constant float *filter,
+ const unsigned int width,
+ const unsigned int imageColumns, const unsigned int imageRows,
+ __local float4 *pixels,
+ const float gain, const float threshold, const unsigned int justBlur)
+ {
+ const int x = get_global_id(0);
+ const int y = get_global_id(1);
+
+ const unsigned int radius = (width - 1) / 2;
+
+ int row = y - radius;
+ int baseRow = get_group_id(1) * get_local_size(1) - radius;
+ int endRow = (get_group_id(1) + 1) * get_local_size(1) + radius;
+
+ while (row < endRow) {
+ int srcy = (row < 0) ? -row : row; // mirror pad
+ srcy = (srcy >= imageRows) ? (2 * imageRows - srcy - 1) : srcy;
+
+ float4 value = 0.0f;
+
+ int ix = x - radius;
+ int i = 0;
+
+ while (i + 7 < width) {
+ for (int j = 0; j < 8; ++j) { // unrolled
+ int srcx = ix + j;
+ srcx = (srcx < 0) ? -srcx : srcx;
+ srcx = (srcx >= imageColumns) ? (2 * imageColumns - srcx - 1) : srcx;
+ value += filter[i + j] * convert_float4(im[srcx + srcy * imageColumns]);
+ }
+ ix += 8;
+ i += 8;
+ }
+
+ while (i < width) {
+ int srcx = (ix < 0) ? -ix : ix; // mirror pad
+ srcx = (srcx >= imageColumns) ? (2 * imageColumns - srcx - 1) : srcx;
+ value += filter[i] * convert_float4(im[srcx + srcy * imageColumns]);
+ ++i;
+ ++ix;
+ }
+ pixels[(row - baseRow) * get_local_size(0) + get_local_id(0)] = value;
+ row += get_local_size(1);
+ }
+
+
+ barrier(CLK_LOCAL_MEM_FENCE);
+
+
+ const int px = get_local_id(0);
+ const int py = get_local_id(1);
+ const int prp = get_local_size(0);
+ float4 value = (float4)(0.0f);
+
+ int i = 0;
+ while (i + 7 < width) { // unrolled
+ value += (float4)(filter[i]) * pixels[px + (py + i) * prp];
+ value += (float4)(filter[i]) * pixels[px + (py + i + 1) * prp];
+ value += (float4)(filter[i]) * pixels[px + (py + i + 2) * prp];
+ value += (float4)(filter[i]) * pixels[px + (py + i + 3) * prp];
+ value += (float4)(filter[i]) * pixels[px + (py + i + 4) * prp];
+ value += (float4)(filter[i]) * pixels[px + (py + i + 5) * prp];
+ value += (float4)(filter[i]) * pixels[px + (py + i + 6) * prp];
+ value += (float4)(filter[i]) * pixels[px + (py + i + 7) * prp];
+ i += 8;
+ }
+ while (i < width) {
+ value += (float4)(filter[i]) * pixels[px + (py + i) * prp];
+ ++i;
+ }
+
+ if (justBlur == 0) { // apply sharpening
+ float4 srcPixel = convert_float4(im[x + y * imageColumns]);
+ float4 diff = srcPixel - value;
+
+ float quantumThreshold = QuantumRange*threshold;
+
+ int4 mask = isless(fabs(2.0f * diff), (float4)quantumThreshold);
+ value = select(srcPixel + diff * gain, srcPixel, mask);
+ }
+
+ if ((x < imageColumns) && (y < imageRows))
+ filtered_im[x + y * imageColumns] = (CLPixelType)(ClampToQuantum(value.s0), ClampToQuantum(value.s1), ClampToQuantum(value.s2), ClampToQuantum(value.s3));
+ }
)
Part of MWC64X by David Thomas, dt10@imperial.ac.uk
This is provided under BSD, full license is with the main package.
See http://www.doc.ic.ac.uk/~dt10/research
-*/\r
-\r
+*/
+
// Pre: a<M, b<M
// Post: r=(a+b) mod M
ulong MWC_AddMod64(ulong a, ulong b, ulong M)
ulong x=MWC_MulMod64(MWC_BASEID, m, M);
return (uint2)((uint)(x/A), (uint)(x%A));
-}\r
-\r
-//! Represents the state of a particular generator\r
+}
+
+//! Represents the state of a particular generator
typedef struct{ uint x; uint c; } mwc64x_state_t;
enum{ MWC64X_A = 4294883355U };
uint res=s->x ^ s->c;
MWC64X_Step(s);
return res;
-}\r
-\r
-//\r
-// End of MWC64X excerpt\r
-//\r
-\r
-\r
- typedef enum\r
- {\r
- UndefinedNoise,\r
- UniformNoise,\r
- GaussianNoise,\r
- MultiplicativeGaussianNoise,\r
- ImpulseNoise,\r
- LaplacianNoise,\r
- PoissonNoise,\r
- RandomNoise\r
- } NoiseType;\r
-\r
-\r
- float mwcReadPseudoRandomValue(mwc64x_state_t* rng) {\r
- return (1.0f * MWC64X_NextUint(rng)) / (float)(0xffffffff); // normalized to 1.0\r
- }\r
-\r
- \r
- float mwcGenerateDifferentialNoise(mwc64x_state_t* r, CLQuantum pixel, NoiseType noise_type, float attenuate) {\r
- \r
- float \r
- alpha,\r
- beta,\r
- noise,\r
- sigma;\r
-\r
- noise = 0.0f;\r
- alpha=mwcReadPseudoRandomValue(r);\r
- switch(noise_type) {\r
- case UniformNoise:\r
- default:\r
- {\r
- noise=(pixel+QuantumRange*SigmaUniform*(alpha-0.5f));\r
- break;\r
- }\r
- case GaussianNoise:\r
- {\r
- float\r
- gamma,\r
- tau;\r
-\r
- if (alpha == 0.0f)\r
- alpha=1.0f;\r
- beta=mwcReadPseudoRandomValue(r);\r
- gamma=sqrt(-2.0f*log(alpha));\r
- sigma=gamma*cospi((2.0f*beta));\r
- tau=gamma*sinpi((2.0f*beta));\r
- noise=(float)(pixel+sqrt((float) pixel)*SigmaGaussian*sigma+\r
- QuantumRange*TauGaussian*tau); \r
- break;\r
- }\r
-\r
-\r
- case ImpulseNoise:\r
- {\r
- if (alpha < (SigmaImpulse/2.0f))\r
- noise=0.0f;\r
- else\r
- if (alpha >= (1.0f-(SigmaImpulse/2.0f)))\r
- noise=(float)QuantumRange;\r
- else\r
- noise=(float)pixel;\r
- break;\r
- }\r
- case LaplacianNoise:\r
- {\r
- if (alpha <= 0.5f)\r
- {\r
- if (alpha <= MagickEpsilon)\r
- noise=(float) (pixel-QuantumRange);\r
- else\r
- noise=(float) (pixel+QuantumRange*SigmaLaplacian*log(2.0f*alpha)+\r
- 0.5f);\r
- break;\r
- }\r
- beta=1.0f-alpha;\r
- if (beta <= (0.5f*MagickEpsilon))\r
- noise=(float) (pixel+QuantumRange);\r
- else\r
- noise=(float) (pixel-QuantumRange*SigmaLaplacian*log(2.0f*beta)+0.5f);\r
- break;\r
- }\r
- case MultiplicativeGaussianNoise:\r
- {\r
- sigma=1.0f;\r
- if (alpha > MagickEpsilon)\r
- sigma=sqrt(-2.0f*log(alpha));\r
- beta=mwcReadPseudoRandomValue(r);\r
- noise=(float) (pixel+pixel*SigmaMultiplicativeGaussian*sigma*\r
- cospi((float) (2.0f*beta))/2.0f);\r
- break;\r
- }\r
- case PoissonNoise:\r
- {\r
- float \r
- poisson;\r
- unsigned int i;\r
- poisson=exp(-SigmaPoisson*QuantumScale*pixel);\r
- for (i=0; alpha > poisson; i++)\r
- {\r
- beta=mwcReadPseudoRandomValue(r);\r
- alpha*=beta;\r
- }\r
- noise=(float) (QuantumRange*i/SigmaPoisson);\r
- break;\r
- }\r
- case RandomNoise:\r
- {\r
- noise=(float) (QuantumRange*SigmaRandom*alpha);\r
- break;\r
- }\r
-\r
- };\r
- return noise;\r
- }\r
-\r
-\r
-\r
-\r
-\r
- __kernel\r
- void GenerateNoiseImage(const __global CLPixelType* inputImage, __global CLPixelType* filteredImage\r
- ,const unsigned int inputPixelCount, const unsigned int pixelsPerWorkItem\r
- ,const ChannelType channel \r
- ,const NoiseType noise_type, const float attenuate\r
- ,const unsigned int seed0, const unsigned int seed1\r
- ,const unsigned int numRandomNumbersPerPixel) {\r
-\r
+}
+
+//
+// End of MWC64X excerpt
+//
+
+
+ typedef enum
+ {
+ UndefinedNoise,
+ UniformNoise,
+ GaussianNoise,
+ MultiplicativeGaussianNoise,
+ ImpulseNoise,
+ LaplacianNoise,
+ PoissonNoise,
+ RandomNoise
+ } NoiseType;
+
+
+ float mwcReadPseudoRandomValue(mwc64x_state_t* rng) {
+ return (1.0f * MWC64X_NextUint(rng)) / (float)(0xffffffff); // normalized to 1.0
+ }
+
+
+ float mwcGenerateDifferentialNoise(mwc64x_state_t* r, CLQuantum pixel, NoiseType noise_type, float attenuate) {
+
+ float
+ alpha,
+ beta,
+ noise,
+ sigma;
+
+ noise = 0.0f;
+ alpha=mwcReadPseudoRandomValue(r);
+ switch(noise_type) {
+ case UniformNoise:
+ default:
+ {
+ noise=(pixel+QuantumRange*SigmaUniform*(alpha-0.5f));
+ break;
+ }
+ case GaussianNoise:
+ {
+ float
+ gamma,
+ tau;
+
+ if (alpha == 0.0f)
+ alpha=1.0f;
+ beta=mwcReadPseudoRandomValue(r);
+ gamma=sqrt(-2.0f*log(alpha));
+ sigma=gamma*cospi((2.0f*beta));
+ tau=gamma*sinpi((2.0f*beta));
+ noise=(float)(pixel+sqrt((float) pixel)*SigmaGaussian*sigma+
+ QuantumRange*TauGaussian*tau);
+ break;
+ }
+
+
+ case ImpulseNoise:
+ {
+ if (alpha < (SigmaImpulse/2.0f))
+ noise=0.0f;
+ else
+ if (alpha >= (1.0f-(SigmaImpulse/2.0f)))
+ noise=(float)QuantumRange;
+ else
+ noise=(float)pixel;
+ break;
+ }
+ case LaplacianNoise:
+ {
+ if (alpha <= 0.5f)
+ {
+ if (alpha <= MagickEpsilon)
+ noise=(float) (pixel-QuantumRange);
+ else
+ noise=(float) (pixel+QuantumRange*SigmaLaplacian*log(2.0f*alpha)+
+ 0.5f);
+ break;
+ }
+ beta=1.0f-alpha;
+ if (beta <= (0.5f*MagickEpsilon))
+ noise=(float) (pixel+QuantumRange);
+ else
+ noise=(float) (pixel-QuantumRange*SigmaLaplacian*log(2.0f*beta)+0.5f);
+ break;
+ }
+ case MultiplicativeGaussianNoise:
+ {
+ sigma=1.0f;
+ if (alpha > MagickEpsilon)
+ sigma=sqrt(-2.0f*log(alpha));
+ beta=mwcReadPseudoRandomValue(r);
+ noise=(float) (pixel+pixel*SigmaMultiplicativeGaussian*sigma*
+ cospi((float) (2.0f*beta))/2.0f);
+ break;
+ }
+ case PoissonNoise:
+ {
+ float
+ poisson;
+ unsigned int i;
+ poisson=exp(-SigmaPoisson*QuantumScale*pixel);
+ for (i=0; alpha > poisson; i++)
+ {
+ beta=mwcReadPseudoRandomValue(r);
+ alpha*=beta;
+ }
+ noise=(float) (QuantumRange*i/SigmaPoisson);
+ break;
+ }
+ case RandomNoise:
+ {
+ noise=(float) (QuantumRange*SigmaRandom*alpha);
+ break;
+ }
+
+ };
+ return noise;
+ }
+
+
+
+
+
+ __kernel
+ void GenerateNoiseImage(const __global CLPixelType* inputImage, __global CLPixelType* filteredImage
+ ,const unsigned int inputPixelCount, const unsigned int pixelsPerWorkItem
+ ,const ChannelType channel
+ ,const NoiseType noise_type, const float attenuate
+ ,const unsigned int seed0, const unsigned int seed1
+ ,const unsigned int numRandomNumbersPerPixel) {
+
mwc64x_state_t rng;
rng.x = seed0;
- rng.c = seed1;\r
-\r
- uint span = pixelsPerWorkItem * numRandomNumbersPerPixel; // length of RNG substream each workitem will use\r
- uint offset = span * get_local_size(0) * get_group_id(0); // offset of this workgroup's RNG substream (in master stream);\r
-\r
- MWC64X_SeedStreams(&rng, offset, span); // Seed the RNG streams\r
-\r
- uint pos = get_local_size(0) * get_group_id(0) * pixelsPerWorkItem + get_local_id(0); // pixel to process\r
-\r
- uint count = pixelsPerWorkItem;\r
-\r
- while (count > 0) {\r
- if (pos < inputPixelCount) {\r
- CLPixelType p = inputImage[pos];\r
-\r
+ rng.c = seed1;
+
+ uint span = pixelsPerWorkItem * numRandomNumbersPerPixel; // length of RNG substream each workitem will use
+ uint offset = span * get_local_size(0) * get_group_id(0); // offset of this workgroup's RNG substream (in master stream);
+
+ MWC64X_SeedStreams(&rng, offset, span); // Seed the RNG streams
+
+ uint pos = get_local_size(0) * get_group_id(0) * pixelsPerWorkItem + get_local_id(0); // pixel to process
+
+ uint count = pixelsPerWorkItem;
+
+ while (count > 0) {
+ if (pos < inputPixelCount) {
+ CLPixelType p = inputImage[pos];
+
if ((channel&RedChannel)!=0) {
setRed(&p,ClampToQuantum(mwcGenerateDifferentialNoise(&rng,getRed(p),noise_type,attenuate)));
}
if ((channel & OpacityChannel) != 0) {
setOpacity(&p,ClampToQuantum(mwcGenerateDifferentialNoise(&rng,getOpacity(p),noise_type,attenuate)));
- }\r
-\r
- filteredImage[pos] = p;\r
- //filteredImage[pos] = (CLPixelType)(MWC64X_NextUint(&rng) % 256, MWC64X_NextUint(&rng) % 256, MWC64X_NextUint(&rng) % 256, 255);\r
- }\r
- pos += get_local_size(0);\r
- --count;\r
- }\r
- }\r
+ }
+
+ filteredImage[pos] = p;
+ //filteredImage[pos] = (CLPixelType)(MWC64X_NextUint(&rng) % 256, MWC64X_NextUint(&rng) % 256, MWC64X_NextUint(&rng) % 256, 255);
+ }
+ pos += get_local_size(0);
+ --count;
+ }
+ }
)
projects / imagemagick / commitdiff