data/kernels/demosaic_amaze.cl
Rama master · réplica de lectura del repositorio original.
/*
This file is part of darktable,
Copyright (C) 2026 darktable developers.
darktable is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
darktable is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with darktable. If not, see <http://www.gnu.org/licenses/>.
*/
// AMaZE demosaic (src/iop/demosaicing/amaze.cc, RawTherapee code by Emil
// Martinec and Ingo Weyrich).
//
// The CPU code works on 160x160 tiles, 128 apart, with a 16 pixel border
// mirrored at the image edges. A tile computes its border pixels itself
// and its later passes read those, so neighboring tiles see different
// values where they overlap. To match that, every buffer here stores each
// tile on its own (TS x TS), and every pass runs per tile with the CPU
// margin to the tile edge (2, 4, 6, ...). Sharing the overlap between
// tiles saves about a third of the memory but differs from the CPU by up to
// 2e-3 near the tile seams.
//
// Three CPU loops update a buffer in place while reading neighbors that the
// same loop already updated (the color difference bounding, the hvwt and pmwt
// neighbor votes), so their result depends on the scan order inside a tile.
// Those run here one CPU tile per work group, in the CPU order (amaze_cd_bound,
// amaze_vote). The Nyquist test has no per-tile bounding box, which changes
// nothing visible.
//
// The CPU stores the R/B-only arrays at half width, indexed by indx >> 1.
// Here they are full size, at the R/B site of each pair; the one place where
// the half indexing reaches further than the neighboring site is rbint in
// the G correction (indx1 -/+ v1 is two rows away), kept as is.
//
// Many decisions compare two nearly equal floats, so contracting to fma
// flips some of them against the CPU, which is built without it for
// amaze.cc (src/iop/CMakeLists.txt).
#include "common.h"
#pragma OPENCL FP_CONTRACT OFF
#define eps 1e-5f
#define epssq 1e-10f
#define arthresh 0.75f
#define nyqthresh 0.5f
// gaussian on 5x5 quincunx, sigma=1.2
constant float gaussodd[4] = { 0.14659727707323927f, 0.103592713382435f, 0.0732036125103057f,
0.0365543548389495f };
// gaussian on 5x5, multiplied with nyqthresh
constant float gaussgrad[6] = { 0.5f * 0.07384411893421103f, 0.5f * 0.06207511968171489f,
0.5f * 0.0521818194747806f, 0.5f * 0.03687419286733595f,
0.5f * 0.03099732204057846f, 0.5f * 0.018413194161458882f };
// gaussian on 5x5 alt quincunx, sigma=1.5
constant float gausseven[2] = { 0.13719494435797422f, 0.05640252782101291f };
// gaussian on quincunx grid
constant float gquinc[4] = { 0.169917f, 0.108947f, 0.069855f, 0.0287182f };
#define TS 160
// a pixel of one CPU tile: the tile at (tx, ty) covers the padded image from
// (128 tx, 128 ty), cc1 x rr1 pixels, and is stored on its own, TS x TS
typedef struct
{
int i; // index into the per-tile buffers
int c, r; // position in the tile
int cc1, rr1; // tile width and height
int col, row; // position in the padded image
bool slot; // false beyond the last tile, where the work size is rounded up
} AmazePix;
AmazePix _pix(const int pw, const int ph)
{
AmazePix p;
const int tx = get_global_id(0) / TS;
const int ty = get_global_id(1) / TS;
p.c = get_global_id(0) - tx * TS;
p.r = get_global_id(1) - ty * TS;
const int ntx = (pw - 16 + 127) / 128;
const int nty = (ph - 16 + 127) / 128;
// the work size may be rounded up beyond the last tile
p.cc1 = tx < ntx && ty < nty ? min(TS, pw - tx * 128) : 0;
p.rr1 = tx < ntx && ty < nty ? min(TS, ph - ty * 128) : 0;
p.col = tx * 128 + p.c;
p.row = ty * 128 + p.r;
p.slot = tx < ntx && ty < nty;
p.i = (ty * ntx + tx) * TS * TS + p.r * TS + p.c;
return p;
}
// at least m pixels inside the tile edge
bool _in(const AmazePix p, const int m)
{
return p.c >= m && p.c < p.cc1 - m && p.r >= m && p.r < p.rr1 - m;
}
// the 16 pixel border mirrors the image without repeating its edge pixel,
// which keeps the CFA phase
int _mirror(int i, int n)
{
i = i < 0 ? -i : i;
return i >= n ? 2 * n - 2 - i : i;
}
float _ulim(float x, float y, float z)
{
return y < z ? clamp(x, y, z) : clamp(x, z, y);
}
float _clampnan(float x)
{
if(isnan(x)) return 0.5f;
return isinf(x) ? clamp(x, 0.0f, 1.0f) : x;
}
// the gradient based directional weights, computed where they are read
// instead of stored
float _dirwts0(global const float *cfa, int i, int v1)
{
return eps + fabs(cfa[i + 2 * v1] - cfa[i]) + fabs(cfa[i] - cfa[i - 2 * v1]) + fabs(cfa[i + v1] - cfa[i - v1]);
}
float _dirwts1(global const float *cfa, int i)
{
return eps + fabs(cfa[i + 2] - cfa[i]) + fabs(cfa[i] - cfa[i - 2]) + fabs(cfa[i + 1] - cfa[i - 1]);
}
// the CPU's delp/delm at an R/B site and Dgrbsq1p/m at a G site, also
// computed where they are read
float _delp(global const float *cfa, int i, int v1) { return fabs(cfa[i - v1 + 1] - cfa[i + v1 - 1]); }
float _delm(global const float *cfa, int i, int v1) { return fabs(cfa[i + v1 + 1] - cfa[i - v1 - 1]); }
float _sq1p(global const float *cfa, int i, int v1)
{
return fsquare(cfa[i] - cfa[i + v1 - 1]) + fsquare(cfa[i] - cfa[i - v1 + 1]);
}
float _sq1m(global const float *cfa, int i, int v1)
{
return fsquare(cfa[i] - cfa[i - v1 - 1]) + fsquare(cfa[i] - cfa[i + v1 + 1]);
}
// the raw data with its mirrored border, and the gradient term of the
// Nyquist test
kernel void amaze_init(
read_only image2d_t in, global float *cfa,
const int width, const int height, const int pw, const int ph)
{
const AmazePix px = _pix(pw, ph);
if(px.c >= px.cc1 || px.r >= px.rr1)
{
if(px.slot) cfa[px.i] = 0.0f;
return;
}
const int x = px.col - 16, y = px.row - 16;
const bool xout = x < 0 || x >= width, yout = y < 0 || y >= height;
// the CPU fills the top left, top right and bottom left corners of the
// border from 16 pixels further in than the edges (amaze.cc:349)
const int sx = x < 0 && yout ? 16 - x : _mirror(x, width);
const int sy = y < 0 && xout ? 16 - y : _mirror(y, height);
cfa[px.i] = Areadsingle(in, sx, sy);
}
kernel void amaze_grad(
global float *cfa, global float *delhvsqsum, const int pw, const int ph)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 2))
{
if(px.slot) { delhvsqsum[px.i] = 0.0f; }
return;
}
const int i = px.i;
delhvsqsum[i] = fsquare(fabs(cfa[i + 1] - cfa[i - 1])) + fsquare(fabs(cfa[i + TS] - cfa[i - TS]));
}
// interpolate vertical and horizontal color differences
kernel void amaze_cd(
global float *cfa, global float *vcd, global float *hcd, global float *vcdalt, global float *hcdalt,
global float *dgintv, global float *dginth,
const int pw, const int ph, const unsigned int filters, const float clip_pt8)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 4))
{
if(px.slot) { vcd[px.i] = 0.0f; hcd[px.i] = 0.0f; vcdalt[px.i] = 0.0f; hcdalt[px.i] = 0.0f; dgintv[px.i] = 0.0f; dginth[px.i] = 0.0f; }
return;
}
const int v1 = TS, v2 = 2 * TS;
const int i = px.i;
const float c0 = cfa[i];
const float d0m2 = _dirwts0(cfa, i - v2, v1), d0 = _dirwts0(cfa, i, v1), d0p2 = _dirwts0(cfa, i + v2, v1);
const float d1m2 = _dirwts1(cfa, i - 2), d1 = _dirwts1(cfa, i), d1p2 = _dirwts1(cfa, i + 2);
// color ratios in each cardinal direction
const float cru = fp32_div(cfa[i - v1] * (d0m2 + d0), d0m2 * (eps + c0) + d0 * (eps + cfa[i - v2]));
const float crd = fp32_div(cfa[i + v1] * (d0p2 + d0), d0p2 * (eps + c0) + d0 * (eps + cfa[i + v2]));
const float crl = fp32_div(cfa[i - 1] * (d1m2 + d1), d1m2 * (eps + c0) + d1 * (eps + cfa[i - 2]));
const float crr = fp32_div(cfa[i + 1] * (d1p2 + d1), d1p2 * (eps + c0) + d1 * (eps + cfa[i + 2]));
// G interpolated in vert/hor directions using Hamilton-Adams method
const float guha = cfa[i - v1] + 0.5f * (c0 - cfa[i - v2]);
const float gdha = cfa[i + v1] + 0.5f * (c0 - cfa[i + v2]);
const float glha = cfa[i - 1] + 0.5f * (c0 - cfa[i - 2]);
const float grha = cfa[i + 1] + 0.5f * (c0 - cfa[i + 2]);
// G interpolated in vert/hor directions using adaptive ratios
float guar = fabs(1.0f - cru) < arthresh ? c0 * cru : guha;
float gdar = fabs(1.0f - crd) < arthresh ? c0 * crd : gdha;
float glar = fabs(1.0f - crl) < arthresh ? c0 * crl : glha;
float grar = fabs(1.0f - crr) < arthresh ? c0 * crr : grha;
// adaptive weights for vertical/horizontal directions
const float hwt = fp32_div(_dirwts1(cfa, i - 1), _dirwts1(cfa, i - 1) + _dirwts1(cfa, i + 1));
const float vwt = fp32_div(_dirwts0(cfa, i - v1, v1), _dirwts0(cfa, i + v1, v1) + _dirwts0(cfa, i - v1, v1));
// interpolated G via adaptive weights of cardinal evaluations
const float Gintvha = vwt * gdha + (1.0f - vwt) * guha;
const float Ginthha = hwt * grha + (1.0f - hwt) * glha;
float v, h, va, ha;
if((FC(row, col, filters) & 1) != 0)
{
v = c0 - (vwt * gdar + (1.0f - vwt) * guar);
h = c0 - (hwt * grar + (1.0f - hwt) * glar);
va = c0 - Gintvha;
ha = c0 - Ginthha;
}
else
{
v = (vwt * gdar + (1.0f - vwt) * guar) - c0;
h = (hwt * grar + (1.0f - hwt) * glar) - c0;
va = Gintvha - c0;
ha = Ginthha - c0;
}
if(c0 > clip_pt8 || Gintvha > clip_pt8 || Ginthha > clip_pt8)
{
// use HA if highlights are (nearly) clipped
guar = guha;
gdar = gdha;
glar = glha;
grar = grha;
v = va;
h = ha;
}
vcd[i] = v;
hcd[i] = h;
vcdalt[i] = va;
hcdalt[i] = ha;
// differences of interpolations in opposite directions
dgintv[i] = fmin(fsquare(guha - gdha), fsquare(guar - gdar));
dginth[i] = fmin(fsquare(glha - grha), fsquare(glar - grar));
}
// the CPU tile of the work group
typedef struct
{
int base; // its first index in the per-tile buffers
int ox, oy; // origin in the padded image
int rr1, cc1; // tile height and width
} AmazeTile;
AmazeTile _tile(int pw, int ph)
{
AmazeTile t;
const int gx = get_group_id(0);
const int gy = get_group_id(1);
t.base = (gy * ((pw - 16 + 127) / 128) + gx) * TS * TS;
t.ox = gx * 128;
t.oy = gy * 128;
t.rr1 = min(TS, ph - t.oy);
t.cc1 = min(TS, pw - t.ox);
return t;
}
// the bounding of one color difference in regions of high saturation, with
// n1 and n2 the raw neighbors on its axis
float _cdbound(float d, float c0, float n1, float n2, bool gsite, float clip_pt)
{
if(gsite)
{
const float Gint = -d + c0; // R or B
if(d > 0.0f)
{
if(3.0f * d > (Gint + c0))
d = -_ulim(Gint, n1, n2) + c0;
else
{
const float wt = 1.0f - fp32_div(3.0f * d, eps + Gint + c0);
d = wt * d + (1.0f - wt) * (-_ulim(Gint, n1, n2) + c0);
}
}
if(Gint > clip_pt) d = -_ulim(Gint, n1, n2) + c0;
}
else
{
const float Gint = d + c0; // interpolated G
if(d < 0.0f)
{
if(3.0f * d < -(Gint + c0))
d = _ulim(Gint, n1, n2) - c0;
else
{
const float wt = 1.0f + fp32_div(3.0f * d, eps + Gint + c0);
d = wt * d + (1.0f - wt) * (_ulim(Gint, n1, n2) - c0);
}
}
if(Gint > clip_pt) d = _ulim(Gint, n1, n2) - c0;
}
return d;
}
float _var3v(float a, float b, float c)
{
return 3.0f * (fsquare(a) + fsquare(b) + fsquare(c)) - fsquare(a + b + c);
}
// choose the smoother color difference and bound it in regions of high
// saturation, one CPU tile per work group. The CPU updates hcd and vcd in
// place, so each value sees the one two pixels before it on its axis
// already done. That makes every row of hcd and every column of vcd two
// chains (one per CFA phase), run here one per thread from the tile edge as
// on the CPU. A pass that reads only the old values differs by
// up to 36/255 on fine detail.
kernel void amaze_cd_bound(
global float *cfa, global float *vcd, global float *hcd, global float *vcdalt, global float *hcdalt,
global float *vcd2, global float *hcd2,
const int pw, const int ph, const unsigned int filters, const float clip_pt)
{
const AmazeTile t = _tile(pw, ph);
const int l = get_local_id(1) * 16 + get_local_id(0);
const int v1 = TS, v2 = 2 * TS;
const int nh = 2 * (t.rr1 - 8); // hcd chains: rows 4 .. rr1-4, two phases each
const int nv = 2 * (t.cc1 - 8); // vcd chains
// the chains cover rows and columns 4 to rr1/cc1 - 5 of the tile
for(int k = l; k < TS * TS; k += 256)
{
const int r = k / TS, c = k - r * TS;
if(r < 4 || r >= t.rr1 - 4 || c < 4 || c >= t.cc1 - 4)
{
hcd2[t.base + k] = 0.0f;
vcd2[t.base + k] = 0.0f;
}
}
for(int job = l; job < nh + nv; job += 256)
{
const bool horiz = job < nh;
const int k = horiz ? job : job - nh;
const int line = 4 + (k >> 1); // tile row (hcd) or column (vcd)
const int len = horiz ? t.cc1 : t.rr1;
const int step = horiz ? 2 : v2;
const int nstep = horiz ? 1 : v1;
// the CPU computes the color differences from 4 pixels inside the tile
// edge and reads zeros beyond that
float prev = 0.0f, aprev = 0.0f;
int p = 4 + (k & 1);
int i = horiz ? t.base + line * TS + p : t.base + p * TS + line;
float cur = horiz ? hcd[i] : vcd[i];
float acur = horiz ? hcdalt[i] : vcdalt[i];
for(; p < len - 4; p += 2, i += step)
{
const bool in = p + 2 < len - 4;
const float next = in ? (horiz ? hcd[i + step] : vcd[i + step]) : 0.0f;
const float anext = in ? (horiz ? hcdalt[i + step] : vcdalt[i + step]) : 0.0f;
float x = _var3v(aprev, acur, anext) < _var3v(prev, cur, next) ? acur : cur;
const int row = horiz ? t.oy + line : t.oy + p;
const int col = horiz ? t.ox + p : t.ox + line;
x = _cdbound(x, cfa[i], cfa[i - nstep], cfa[i + nstep], (FC(row, col, filters) & 1) != 0, clip_pt);
if(horiz) hcd2[i] = x;
else vcd2[i] = x;
prev = x;
aprev = acur;
cur = next;
acur = anext;
}
}
}
// the CPU's hvwt and pmwt votes: an R/B site takes the average of its four
// diagonal neighbors when that discriminates the direction better. They
// update in place, so the two neighbors above are already done; here one
// work group walks the rows of one CPU tile, from row m to rr1 - m. A
// pass that reads only the old values differs by up to 42/255.
kernel void amaze_vote(
global float *src, global float *dst,
const int pw, const int ph, const unsigned int filters, const int m)
{
const AmazeTile t = _tile(pw, ph);
const int c = get_local_id(1) * 16 + get_local_id(0); // tile column
local float vote_row[160];
const bool active = c < t.cc1;
if(active) vote_row[c] = src[t.base + (m - 1) * TS + c];
barrier(CLK_LOCAL_MEM_FENCE);
for(int r = 0; r < t.rr1; r++)
{
const int i = t.base + r * TS + c;
const bool rb = active && (FC(t.oy + r, t.ox + c, filters) & 1) == 0;
float w = rb ? src[i] : 0.0f;
if(rb && r >= m && r < t.rr1 - m && c >= m && c < t.cc1 - m)
{
const float alt = 0.25f * (vote_row[c - 1] + vote_row[c + 1] + src[i + TS - 1] + src[i + TS + 1]);
if(fabs(0.5f - w) < fabs(0.5f - alt)) w = alt;
}
if(c < TS) dst[i] = w;
barrier(CLK_LOCAL_MEM_FENCE);
// the R/B sites of this row are the other phase of the previous one
if(rb && r >= m - 1) vote_row[c] = w;
barrier(CLK_LOCAL_MEM_FENCE);
}
}
// the square of the color difference, computed where it is read
float _cdd(global const float *vcd, global const float *hcd, int i) { return fsquare(vcd[i] - hcd[i]); }
// R/B sites: the adaptive weight of the vertical against the horizontal
// interpolation, and the Nyquist texture test
kernel void amaze_hvwt(
global float *cfa, global float *vcd, global float *hcd, global float *dgintv, global float *dginth,
global float *delhvsqsum, global float *hvwt, global float *nyquist,
const int pw, const int ph, const unsigned int filters)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 6) || (FC(row, col, filters) & 1) != 0)
{
if(px.slot) { hvwt[px.i] = 0.0f; nyquist[px.i] = 0.0f; }
return;
}
const int ts = TS;
const int v1 = TS, v2 = 2 * TS, v3 = 3 * TS;
const int p1 = -ts + 1, p2 = -2 * ts + 2, m1 = ts + 1, m2 = 2 * ts + 2;
const int i = px.i;
// color difference variances in cardinal directions
const float uave = vcd[i] + vcd[i - v1] + vcd[i - v2] + vcd[i - v3];
const float dave = vcd[i] + vcd[i + v1] + vcd[i + v2] + vcd[i + v3];
const float lave = hcd[i] + hcd[i - 1] + hcd[i - 2] + hcd[i - 3];
const float rave = hcd[i] + hcd[i + 1] + hcd[i + 2] + hcd[i + 3];
float Dgrbvvaru = fsquare(vcd[i] - uave) + fsquare(vcd[i - v1] - uave) + fsquare(vcd[i - v2] - uave)
+ fsquare(vcd[i - v3] - uave);
float Dgrbvvard = fsquare(vcd[i] - dave) + fsquare(vcd[i + v1] - dave) + fsquare(vcd[i + v2] - dave)
+ fsquare(vcd[i + v3] - dave);
float Dgrbhvarl = fsquare(hcd[i] - lave) + fsquare(hcd[i - 1] - lave) + fsquare(hcd[i - 2] - lave)
+ fsquare(hcd[i - 3] - lave);
float Dgrbhvarr = fsquare(hcd[i] - rave) + fsquare(hcd[i + 1] - rave) + fsquare(hcd[i + 2] - rave)
+ fsquare(hcd[i + 3] - rave);
const float hwt = fp32_div(_dirwts1(cfa, i - 1), _dirwts1(cfa, i - 1) + _dirwts1(cfa, i + 1));
const float vwt = fp32_div(_dirwts0(cfa, i - v1, v1), _dirwts0(cfa, i + v1, v1) + _dirwts0(cfa, i - v1, v1));
const float vcdvar = epssq + vwt * Dgrbvvard + (1.0f - vwt) * Dgrbvvaru;
const float hcdvar = epssq + hwt * Dgrbhvarr + (1.0f - hwt) * Dgrbhvarl;
// fluctuations in up/down and left/right interpolations of colors
Dgrbvvaru = dgintv[i] + dgintv[i - v1] + dgintv[i - v2];
Dgrbvvard = dgintv[i] + dgintv[i + v1] + dgintv[i + v2];
Dgrbhvarl = dginth[i] + dginth[i - 1] + dginth[i - 2];
Dgrbhvarr = dginth[i] + dginth[i + 1] + dginth[i + 2];
const float vcdvar1 = epssq + vwt * Dgrbvvard + (1.0f - vwt) * Dgrbvvaru;
const float hcdvar1 = epssq + hwt * Dgrbhvarr + (1.0f - hwt) * Dgrbhvarl;
// adaptive weights for G interpolation
const float varwt = fp32_div(hcdvar, vcdvar + hcdvar);
const float diffwt = fp32_div(hcdvar1, vcdvar1 + hcdvar1);
// if both agree on the direction, choose the one with the strongest
// discrimination, otherwise the u/d and l/r difference fluctuation weights
hvwt[i] = ((0.5f - varwt) * (0.5f - diffwt) > 0.0f && fabs(0.5f - diffwt) < fabs(0.5f - varwt)) ? varwt : diffwt;
const float nyqutest
= (gaussodd[0] * _cdd(vcd, hcd, i)
+ gaussodd[1] * (_cdd(vcd, hcd, i - m1) + _cdd(vcd, hcd, i + p1) + _cdd(vcd, hcd, i - p1) + _cdd(vcd, hcd, i + m1))
+ gaussodd[2] * (_cdd(vcd, hcd, i - v2) + _cdd(vcd, hcd, i - 2) + _cdd(vcd, hcd, i + 2) + _cdd(vcd, hcd, i + v2))
+ gaussodd[3] * (_cdd(vcd, hcd, i - m2) + _cdd(vcd, hcd, i + p2) + _cdd(vcd, hcd, i - p2) + _cdd(vcd, hcd, i + m2)))
- (gaussgrad[0] * delhvsqsum[i]
+ gaussgrad[1] * (delhvsqsum[i - v1] + delhvsqsum[i + 1] + delhvsqsum[i - 1] + delhvsqsum[i + v1])
+ gaussgrad[2] * (delhvsqsum[i - m1] + delhvsqsum[i + p1] + delhvsqsum[i - p1] + delhvsqsum[i + m1])
+ gaussgrad[3] * (delhvsqsum[i - v2] + delhvsqsum[i - 2] + delhvsqsum[i + 2] + delhvsqsum[i + v2])
+ gaussgrad[4] * (delhvsqsum[i - v2 - 1] + delhvsqsum[i - v2 + 1] + delhvsqsum[i - ts - 2]
+ delhvsqsum[i - ts + 2] + delhvsqsum[i + ts - 2] + delhvsqsum[i + ts + 2]
+ delhvsqsum[i + v2 - 1] + delhvsqsum[i + v2 + 1])
+ gaussgrad[5] * (delhvsqsum[i - m2] + delhvsqsum[i + p2] + delhvsqsum[i - p2] + delhvsqsum[i + m2]));
nyquist[i] = nyqutest > 0.0f ? 1.0f : 0.0f;
}
// if most of your neighbors are named Nyquist, it's likely that you're one
// too, or not
kernel void amaze_nyquist2(
global float *nyquist, global float *nyquist2, const int pw, const int ph, const unsigned int filters)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 8) || (FC(row, col, filters) & 1) != 0)
{
if(px.slot) nyquist2[px.i] = 0.0f;
return;
}
const int ts = TS;
const int v2 = 2 * TS, p1 = -ts + 1, m1 = ts + 1;
const int i = px.i;
const float n = nyquist[i - v2] + nyquist[i - m1] + nyquist[i + p1] + nyquist[i - 2] + nyquist[i + 2]
+ nyquist[i - p1] + nyquist[i + m1] + nyquist[i + v2];
nyquist2[i] = n > 4.0f ? 1.0f : (n < 4.0f ? 0.0f : nyquist[i]);
}
// in areas of Nyquist texture, do area interpolation
kernel void amaze_area(
global float *cfa, global float *nyquist2, global float *hvwt,
const int pw, const int ph, const unsigned int filters)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 8) || (FC(row, col, filters) & 1) != 0) return;
const int v1 = TS;
const int i = px.i;
if(nyquist2[i] == 0.0f) return;
float sumcfa = 0.0f, sumh = 0.0f, sumv = 0.0f, sumsqh = 0.0f, sumsqv = 0.0f, areawt = 0.0f;
for(int y = -6; y < 7; y += 2)
for(int x = -6; x < 7; x += 2)
{
const int k = i + y * TS + x;
if(nyquist2[k] != 0.0f)
{
const float c = cfa[k];
sumcfa += c;
sumh += cfa[k - 1] + cfa[k + 1];
sumv += cfa[k - v1] + cfa[k + v1];
sumsqh += fsquare(c - cfa[k - 1]) + fsquare(c - cfa[k + 1]);
sumsqv += fsquare(c - cfa[k - v1]) + fsquare(c - cfa[k + v1]);
areawt += 1.0f;
}
}
// horizontal and vertical color differences, and adaptive weight
sumh = sumcfa - 0.5f * sumh;
sumv = sumcfa - 0.5f * sumv;
areawt = 0.5f * areawt;
const float hcdvar = epssq + fabs(areawt * sumsqh - sumh * sumh);
const float vcdvar = epssq + fabs(areawt * sumsqv - sumv * sumv);
hvwt[i] = fp32_div(hcdvar, vcdvar + hcdvar);
}
// populate G at R/B sites
kernel void amaze_green(
global float *cfa, global float *vcd, global float *hcd, global float *hvwt, global float *nyquist2,
global float *dgrb, global float *green, global float *dgrb2h, global float *dgrb2v,
const int pw, const int ph, const unsigned int filters)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 8) || (FC(row, col, filters) & 1) != 0)
{
if(px.slot) { dgrb[px.i] = 0.0f; green[px.i] = 0.0f; dgrb2h[px.i] = 0.0f; dgrb2v[px.i] = 0.0f; }
return;
}
const int v1 = TS;
const int i = px.i;
const float w = hvwt[i]; // after amaze_vote
const float d = hcd[i] + w * (vcd[i] - hcd[i]); // evaluate color differences
const float g = cfa[i] + d; // evaluate G (finally!)
dgrb[i] = d;
green[i] = g;
// local curvature in G (preparation for nyquist refinement step); the
// neighbors are G sites, where G is the raw value
const bool nyq = nyquist2[i] != 0.0f;
dgrb2h[i] = nyq ? fsquare(g - 0.5f * (cfa[i - 1] + cfa[i + 1])) : 0.0f;
dgrb2v[i] = nyq ? fsquare(g - 0.5f * (cfa[i - v1] + cfa[i + v1])) : 0.0f;
}
float _gvar(global const float *b, int i, int ts)
{
const int v2 = 2 * ts, p1 = -ts + 1, p2 = -2 * ts + 2, m1 = ts + 1, m2 = 2 * ts + 2;
return epssq + (gquinc[0] * b[i]
+ gquinc[1] * (b[i - m1] + b[i + p1] + b[i - p1] + b[i + m1])
+ gquinc[2] * (b[i - v2] + b[i - 2] + b[i + 2] + b[i + v2])
+ gquinc[3] * (b[i - m2] + b[i + p2] + b[i - p2] + b[i + m2]));
}
// refine Nyquist areas using G curvatures
kernel void amaze_nyquist_refine(
global float *cfa, global float *vcd, global float *hcd, global float *nyquist2,
global float *dgrb2h, global float *dgrb2v, global float *dgrb, global float *green,
const int pw, const int ph, const unsigned int filters)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 8) || (FC(row, col, filters) & 1) != 0) return;
const int i = px.i;
if(nyquist2[i] == 0.0f) return;
// local averages (over Nyquist pixels only) of G curvature squared
const float gvarh = _gvar(dgrb2h, i, TS);
const float gvarv = _gvar(dgrb2v, i, TS);
// use the results as weights for refined G interpolation
const float d = fp32_div((hcd[i] * gvarv + vcd[i] * gvarh), gvarv + gvarh);
dgrb[i] = d;
green[i] = cfa[i] + d;
}
float _rbdiag(float c, float cr, float n1, float n2)
{
return fabs(1.0f - cr) < arthresh ? c * cr : n1 + 0.5f * (c - n2);
}
float _rbvar(global const float *cfa, int i, int ts, bool plus)
{
const int v1 = ts, v2 = 2 * ts;
const int k[12] = { i - v1, i - 1, i + 1, i + v1,
i - v2 - 1, i - v2 + 1, i - 2 - v1, i + 2 - v1, i - 2 + v1, i + 2 + v1, i + v2 - 1, i + v2 + 1 };
float s0 = 0.0f, s1 = 0.0f;
for(int n = 0; n < 4; n++) s0 += plus ? _sq1p(cfa, k[n], ts) : _sq1m(cfa, k[n], ts);
for(int n = 4; n < 12; n++) s1 += plus ? _sq1p(cfa, k[n], ts) : _sq1m(cfa, k[n], ts);
return epssq + (gausseven[0] * s0 + gausseven[1] * s1);
}
// diagonal interpolation correction
kernel void amaze_diag(
global float *cfa, global float *rbm, global float *rbp, global float *pmwt,
const int pw, const int ph, const unsigned int filters, const float clip_pt)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 8) || (FC(row, col, filters) & 1) != 0)
{
if(px.slot) { rbm[px.i] = 0.0f; rbp[px.i] = 0.0f; pmwt[px.i] = 0.0f; }
return;
}
const int ts = TS;
const int p1 = -ts + 1, p2 = -2 * ts + 2, m1 = ts + 1, m2 = 2 * ts + 2;
const int i = px.i;
const float c0 = cfa[i];
// diagonal color ratios
const float crse = fp32_div(2.0f * cfa[i + m1], eps + c0 + cfa[i + m2]);
const float crnw = fp32_div(2.0f * cfa[i - m1], eps + c0 + cfa[i - m2]);
const float crne = fp32_div(2.0f * cfa[i + p1], eps + c0 + cfa[i + p2]);
const float crsw = fp32_div(2.0f * cfa[i - p1], eps + c0 + cfa[i - p2]);
// assign B/R at R/B sites
const float rbse = _rbdiag(c0, crse, cfa[i + m1], cfa[i + m2]);
const float rbnw = _rbdiag(c0, crnw, cfa[i - m1], cfa[i - m2]);
const float rbne = _rbdiag(c0, crne, cfa[i + p1], cfa[i + p2]);
const float rbsw = _rbdiag(c0, crsw, cfa[i - p1], cfa[i - p2]);
const float dm = _delm(cfa, i, ts), dp = _delp(cfa, i, ts);
const float wtse = eps + dm + _delm(cfa, i + m1, ts) + _delm(cfa, i + m2, ts);
const float wtnw = eps + dm + _delm(cfa, i - m1, ts) + _delm(cfa, i - m2, ts);
const float wtne = eps + dp + _delp(cfa, i + p1, ts) + _delp(cfa, i + p2, ts);
const float wtsw = eps + dp + _delp(cfa, i - p1, ts) + _delp(cfa, i - p2, ts);
float m = fp32_div((wtse * rbnw + wtnw * rbse), wtse + wtnw);
float p = fp32_div((wtne * rbsw + wtsw * rbne), wtne + wtsw);
// variance of R-B in plus/minus directions
const float rbvarm = _rbvar(cfa, i, ts, false);
pmwt[i] = fp32_div(rbvarm, _rbvar(cfa, i, ts, true) + rbvarm);
// bound the interpolation in regions of high saturation
if(p < c0)
{
if(2.0f * p < c0)
p = _ulim(p, cfa[i - p1], cfa[i + p1]);
else
{
const float pwt = fp32_div(2.0f * (c0 - p), eps + p + c0);
p = pwt * p + (1.0f - pwt) * _ulim(p, cfa[i - p1], cfa[i + p1]);
}
}
if(m < c0)
{
if(2.0f * m < c0)
m = _ulim(m, cfa[i - m1], cfa[i + m1]);
else
{
const float mwt = fp32_div(2.0f * (c0 - m), eps + m + c0);
m = mwt * m + (1.0f - mwt) * _ulim(m, cfa[i - m1], cfa[i + m1]);
}
}
if(p > clip_pt) p = _ulim(p, cfa[i - p1], cfa[i + p1]);
if(m > clip_pt) m = _ulim(m, cfa[i - m1], cfa[i + m1]);
rbm[i] = m;
rbp[i] = p;
}
// R+B interpolated
kernel void amaze_rbint(
global float *cfa, global float *rbm, global float *rbp, global float *pmwt,
global float *rbint, const int pw, const int ph, const unsigned int filters)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 10) || (FC(row, col, filters) & 1) != 0)
{
if(px.slot) { rbint[px.i] = 0.0f; }
return;
}
const int i = px.i;
const float w = pmwt[i]; // after amaze_vote
rbint[i] = 0.5f * (cfa[i] + rbm[i] * (1.0f - w) + rbp[i] * w);
}
float _gdir(float rb0, float cr, float c, float rbn)
{
return fabs(1.0f - cr) < arthresh ? rb0 * cr : c + 0.5f * (rb0 - rbn);
}
// interpolate G vertically/horizontally using R+B values where the
// diagonal discrimination is stronger
kernel void amaze_green_diag(
global float *cfa, global float *hvwt, global float *pmwt, global float *rbint,
global float *dgrb, global float *green,
const int pw, const int ph, const unsigned int filters, const float clip_pt)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
if(!_in(px, 12) || (FC(row, col, filters) & 1) != 0) return;
const int v1 = TS;
const int i = px.i;
if(fabs(0.5f - pmwt[i]) < fabs(0.5f - hvwt[i])) return;
// the CPU's rbint[indx1 -/+ v1] is two rows away in its half width
// array, and rbint[indx1 -/+ 1] two columns
const float rb0 = rbint[i];
const float rbu = rbint[i - 2 * v1], rbd = rbint[i + 2 * v1], rbl = rbint[i - 2], rbr = rbint[i + 2];
// color ratios for G interpolation
const float cru = fp32_div(cfa[i - v1] * 2.0f, eps + rb0 + rbu);
const float crd = fp32_div(cfa[i + v1] * 2.0f, eps + rb0 + rbd);
const float crl = fp32_div(cfa[i - 1] * 2.0f, eps + rb0 + rbl);
const float crr = fp32_div(cfa[i + 1] * 2.0f, eps + rb0 + rbr);
// interpolated G via adaptive ratios or Hamilton-Adams in each cardinal direction
const float gu = _gdir(rb0, cru, cfa[i - v1], rbu);
const float gd = _gdir(rb0, crd, cfa[i + v1], rbd);
const float gl = _gdir(rb0, crl, cfa[i - 1], rbl);
const float gr = _gdir(rb0, crr, cfa[i + 1], rbr);
// interpolated G via adaptive weights of cardinal evaluations
const float d0u = _dirwts0(cfa, i - v1, v1), d0d = _dirwts0(cfa, i + v1, v1);
const float d1l = _dirwts1(cfa, i - 1), d1r = _dirwts1(cfa, i + 1);
float Gintv = fp32_div((d0u * gd + d0d * gu), d0d + d0u);
float Ginth = fp32_div((d1l * gr + d1r * gl), d1l + d1r);
// bound the interpolation in regions of high saturation
if(Gintv < rb0)
{
if(2.0f * Gintv < rb0)
Gintv = _ulim(Gintv, cfa[i - v1], cfa[i + v1]);
else
{
const float vwt = fp32_div(2.0f * (rb0 - Gintv), eps + Gintv + rb0);
Gintv = vwt * Gintv + (1.0f - vwt) * _ulim(Gintv, cfa[i - v1], cfa[i + v1]);
}
}
if(Ginth < rb0)
{
if(2.0f * Ginth < rb0)
Ginth = _ulim(Ginth, cfa[i - 1], cfa[i + 1]);
else
{
const float hwt = fp32_div(2.0f * (rb0 - Ginth), eps + Ginth + rb0);
Ginth = hwt * Ginth + (1.0f - hwt) * _ulim(Ginth, cfa[i - 1], cfa[i + 1]);
}
}
if(Ginth > clip_pt) Ginth = _ulim(Ginth, cfa[i - 1], cfa[i + 1]);
if(Gintv > clip_pt) Gintv = _ulim(Gintv, cfa[i - v1], cfa[i + v1]);
const float g = Ginth * (1.0f - hvwt[i]) + Gintv * hvwt[i];
green[i] = g;
dgrb[i] = g - cfa[i];
}
float _chroma(global const float *d, int i, int ts)
{
const int v2 = 2 * ts, p1 = -ts + 1, p3 = -3 * ts + 3, m1 = ts + 1, m3 = 3 * ts + 3;
const float wtnw = fp32_div(1.0f, eps + fabs(d[i - m1] - d[i + m1]) + fabs(d[i - m1] - d[i - m3]) + fabs(d[i + m1] - d[i - m3]));
const float wtne = fp32_div(1.0f, eps + fabs(d[i + p1] - d[i - p1]) + fabs(d[i + p1] - d[i + p3]) + fabs(d[i - p1] - d[i + p3]));
const float wtsw = fp32_div(1.0f, eps + fabs(d[i - p1] - d[i + p1]) + fabs(d[i - p1] - d[i + m3]) + fabs(d[i + p1] - d[i - p3]));
const float wtse = fp32_div(1.0f, eps + fabs(d[i + m1] - d[i - m1]) + fabs(d[i + m1] - d[i - p3]) + fabs(d[i - m1] - d[i + m3]));
return fp32_div(wtnw * (1.325f * d[i - m1] - 0.175f * d[i - m3] - 0.075f * d[i - m1 - 2] - 0.075f * d[i - m1 - v2])
+ wtne * (1.325f * d[i + p1] - 0.175f * d[i + p3] - 0.075f * d[i + p1 + 2] - 0.075f * d[i + p1 + v2])
+ wtsw * (1.325f * d[i - p1] - 0.175f * d[i - p3] - 0.075f * d[i - p1 - 2] - 0.075f * d[i - p1 - v2])
+ wtse * (1.325f * d[i + m1] - 0.175f * d[i + m3] - 0.075f * d[i + m1 + 2] - 0.075f * d[i + m1 + v2]),
wtnw + wtne + wtsw + wtse);
}
// fancy chrominance interpolation: G-R and G-B at every R/B site, the own
// color difference from dgrb and the other from the diagonal neighbors,
// which are of the other color
kernel void amaze_chroma(
global float *dgrb, global float *dgr, global float *dgb,
const int pw, const int ph, const unsigned int filters)
{
const AmazePix px = _pix(pw, ph);
const int col = px.col, row = px.row;
const int c = FC(row, col, filters);
if(!_in(px, 14) || (c & 1) != 0)
{
if(px.slot) { dgr[px.i] = 0.0f; dgb[px.i] = 0.0f; }
return;
}
const int i = px.i;
const float own = dgrb[i];
const float other = _chroma(dgrb, i, TS);
dgr[i] = c == 0 ? own : other;
dgb[i] = c == 0 ? other : own;
}
kernel void amaze_output(
global float *cfa, global float *green, global float *hvwt, global float *dgr, global float *dgb,
write_only image2d_t out, const int width, const int height, const int pw, const unsigned int filters)
{
const int x = get_global_id(0);
const int y = get_global_id(1);
if(x >= width || y >= height) return;
const int row = y + 16, col = x + 16;
// the tile whose output this pixel is: its rows and columns 16 to 143
const int ntx = (pw - 16 + 127) / 128;
const int tx = x / 128, ty = y / 128;
const int v1 = TS;
const int i = (ty * ntx + tx) * TS * TS + (row - ty * 128) * TS + col - tx * 128;
float r, g, b;
if((FC(row, col, filters) & 1) == 0)
{ // R/B site
g = green[i];
r = g - dgr[i];
b = g - dgb[i];
}
else
{ // G site: the color differences of the four R/B neighbors
g = cfa[i];
const float wu = hvwt[i - v1], wd = hvwt[i + v1], wl = 1.0f - hvwt[i - 1], wr = 1.0f - hvwt[i + 1];
const float temp = fp32_div(1.0f, wu + 2.0f - hvwt[i + 1] - hvwt[i - 1] + wd);
r = g - (wu * dgr[i - v1] + wr * dgr[i + 1] + wl * dgr[i - 1] + wd * dgr[i + v1]) * temp;
b = g - (wu * dgb[i - v1] + wr * dgb[i + 1] + wl * dgb[i - 1] + wd * dgb[i + v1]) * temp;
}
write_imagef(out, (int2)(x, y), (float4)(_clampnan(r), _clampnan(g), _clampnan(b), 0.0f));
}