diff options
13 files changed, 326 insertions, 276 deletions
diff --git a/source/blender/blenlib/BLI_math_geom.h b/source/blender/blenlib/BLI_math_geom.h index 9073fbb8b80..ad7dab0d7a6 100644 --- a/source/blender/blenlib/BLI_math_geom.h +++ b/source/blender/blenlib/BLI_math_geom.h @@ -222,6 +222,8 @@ int barycentric_inside_triangle_v2(const float w[3]); void resolve_tri_uv(float r_uv[2], const float st[2], const float st0[2], const float st1[2], const float st2[2]); void resolve_quad_uv(float uv[2], const float st[2], const float st0[2], const float st1[2], const float st2[2], const float st3[2]); +void resolve_quad_uv_deriv(float r_uv[2], float r_deriv[2][2], + const float st[2], const float st0[2], const float st1[2], const float st2[2], const float st3[2]); /* use to find the point of a UV on a face */ void interp_bilinear_quad_v3(float data[4][3], float u, float v, float res[3]); diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c index 45612335323..c4c0adf2666 100644 --- a/source/blender/blenlib/intern/math_geom.c +++ b/source/blender/blenlib/intern/math_geom.c @@ -2685,6 +2685,13 @@ void resolve_tri_uv(float r_uv[2], const float st[2], const float st0[2], const /* bilinear reverse */ void resolve_quad_uv(float r_uv[2], const float st[2], const float st0[2], const float st1[2], const float st2[2], const float st3[2]) { + resolve_quad_uv_deriv(r_uv, NULL, st, st0, st1, st2, st3); +} + +/* bilinear reverse with derivatives */ +void resolve_quad_uv_deriv(float r_uv[2], float r_deriv[2][2], + const float st[2], const float st0[2], const float st1[2], const float st2[2], const float st3[2]) +{ const double signed_area = (st0[0] * st1[1] - st0[1] * st1[0]) + (st1[0] * st2[1] - st1[1] * st2[0]) + (st2[0] * st3[1] - st2[1] * st3[0]) + (st3[0] * st0[1] - st3[1] * st0[0]); @@ -2732,6 +2739,32 @@ void resolve_quad_uv(float r_uv[2], const float st[2], const float st0[2], const if (IS_ZERO(denom) == 0) r_uv[1] = (float)((double)((1.0f - r_uv[0]) * (st0[i] - st[i]) + r_uv[0] * (st1[i] - st[i])) / denom); } + + if (r_deriv) { + float tmp1[2], tmp2[2], s[2], t[2]; + double denom; + + /* clear outputs */ + zero_v2(r_deriv[0]); + zero_v2(r_deriv[1]); + + sub_v2_v2v2(tmp1, st1, st0); + sub_v2_v2v2(tmp2, st2, st3); + interp_v2_v2v2(s, tmp1, tmp2, r_uv[1]); + sub_v2_v2v2(tmp1, st3, st0); + sub_v2_v2v2(tmp2, st2, st1); + interp_v2_v2v2(t, tmp1, tmp2, r_uv[0]); + + denom = t[0]*s[1] - t[1]*s[0]; + + if (!IS_ZERO(denom)) { + double inv_denom = 1.0 / denom; + r_deriv[0][0] = -t[1] * inv_denom; + r_deriv[0][1] = t[0] * inv_denom; + r_deriv[1][0] = s[1] * inv_denom; + r_deriv[1][1] = -s[0] * inv_denom; + } + } } #undef IS_ZERO diff --git a/source/blender/compositor/intern/COM_MemoryBuffer.cpp b/source/blender/compositor/intern/COM_MemoryBuffer.cpp index f10e6696c6a..8a010561adf 100644 --- a/source/blender/compositor/intern/COM_MemoryBuffer.cpp +++ b/source/blender/compositor/intern/COM_MemoryBuffer.cpp @@ -214,136 +214,108 @@ static const float EWA_WTS[EWA_MAXIDX + 1] = { 0.00754159f, 0.00625989f, 0.00498819f, 0.00372644f, 0.00247454f, 0.00123242f, 0.f }; -static void radangle2imp(float a2, float b2, float th, float *A, float *B, float *C, float *F) +static void ellipse_bounds(float A, float B, float C, float F, float &xmax, float &ymax) { - float ct2 = cosf(th); - const float st2 = 1.f - ct2 * ct2; // <- sin(th)^2 - ct2 *= ct2; - *A = a2 * st2 + b2 * ct2; - *B = (b2 - a2) * sinf(2.f * th); - *C = a2 * ct2 + b2 * st2; - *F = a2 * b2; -} - -// all tests here are done to make sure possible overflows are hopefully minimized -static void imp2radangle(float A, float B, float C, float F, float *a, float *b, float *th, float *ecc) -{ - if (F <= 1e-5f) { // use arbitrary major radius, zero minor, infinite eccentricity - *a = sqrtf(A > C ? A : C); - *b = 0.f; - *ecc = 1e10f; - *th = 0.5f * (atan2f(B, A - C) + (float)M_PI); + float denom = 4.0f*A*C - B*B; + if (denom > 0.0f && A != 0.0f && C != 0.0f) { + xmax = sqrt(F)/(2.0f*A) * (sqrt(F*(4.0f*A - B*B/C)) + B*B*sqrt(F/(C*denom))); + ymax = sqrt(F)/(2.0f*C) * (sqrt(F*(4.0f*C - B*B/A)) + B*B*sqrt(F/(A*denom))); } else { - const float AmC = A - C, ApC = A + C, F2 = F * 2.f; - const float r = sqrtf(AmC * AmC + B * B); - float d = ApC - r; - *a = (d <= 0.f) ? sqrtf(A > C ? A : C) : sqrtf(F2 / d); - d = ApC + r; - if (d <= 0.f) { - *b = 0.f; - *ecc = 1e10f; - } - else { - *b = sqrtf(F2 / d); - *ecc = *a / *b; - } - /* incr theta by 0.5 * pi (angle of major axis) */ - *th = 0.5f * (atan2f(B, AmC) + (float)M_PI); + xmax = 0.0f; + ymax = 0.0f; } } -static float clipuv(float x, float limit) +static void ellipse_params(float Ux, float Uy, float Vx, float Vy, float &A, float &B, float &C, float &F, float &umax, float &vmax) { - x = (x < 0) ? 0 : ((x >= limit) ? (limit - 1) : x); - return x; + A = Vx*Vx + Vy*Vy; + B = -2.0f * (Ux*Vx + Uy*Vy); + C = Ux*Ux + Uy*Uy; + F = A*C - B*B * 0.25f; + + float factor = (F != 0.0f ? (float)(EWA_MAXIDX+1) / F : 0.0f); + A *= factor; + B *= factor; + C *= factor; + F = (float)(EWA_MAXIDX+1); + + ellipse_bounds(A, B, C, sqrtf(F), umax, vmax); } /** - * \note \a sampler at the moment is either 'COM_PS_NEAREST' or not, other values won't matter. + * Filtering method based on + * "Creating raster omnimax images from multiple perspective views using the elliptical weighted average filter" + * by Ned Greene and Paul S. Heckbert (1986) */ -void MemoryBuffer::readEWA(float result[4], float fx, float fy, float dx, float dy, PixelSampler sampler) +void MemoryBuffer::readEWA(float result[4], const float uv[2], const float derivatives[2][2], PixelSampler sampler) { - const int width = this->getWidth(), height = this->getHeight(); - - // scaling dxt/dyt by full resolution can cause overflow because of huge A/B/C and esp. F values, - // scaling by aspect ratio alone does the opposite, so try something in between instead... - const float ff2 = width, ff = sqrtf(ff2), q = height / ff; - const float Ux = dx * ff, Vx = dx * q, Uy = dy * ff, Vy = dy * q; - float A = Vx * Vx + Vy * Vy; - float B = -2.f * (Ux * Vx + Uy * Vy); - float C = Ux * Ux + Uy * Uy; - float F = A * C - B * B * 0.25f; - float a, b, th, ecc, a2, b2, ue, ve, U0, V0, DDQ, U, ac1, ac2, BU, d; - int u, v, u1, u2, v1, v2; - // The so-called 'high' quality ewa method simply adds a constant of 1 to both A & C, - // so the ellipse always covers at least some texels. But since the filter is now always larger, - // it also means that everywhere else it's also more blurry then ideally should be the case. - // So instead here the ellipse radii are modified instead whenever either is too low. - // Use a different radius based on interpolation switch, just enough to anti-alias when interpolation is off, - // and slightly larger to make result a bit smoother than bilinear interpolation when interpolation is on - // (minimum values: const float rmin = intpol ? 1.f : 0.5f;) - - /* note: 0.765625f is too sharp, 1.0 will not blur with an exact pixel sample - * useful to avoid blurring when there is no distortion */ -#if 0 - const float rmin = ((sampler != COM_PS_NEAREST) ? 1.5625f : 0.765625f) / ff2; -#else - const float rmin = ((sampler != COM_PS_NEAREST) ? 1.5625f : 1.0f ) / ff2; -#endif - imp2radangle(A, B, C, F, &a, &b, &th, &ecc); - if ((b2 = b * b) < rmin) { - if ((a2 = a * a) < rmin) { - B = 0.f; - A = C = rmin; - F = A * C; - } - else { - b2 = rmin; - radangle2imp(a2, b2, th, &A, &B, &C, &F); - } - } + zero_v4(result); + int width = this->getWidth(), height = this->getHeight(); + if (width == 0 || height == 0) + return; - ue = ff * sqrtf(C); - ve = ff * sqrtf(A); - d = (float)(EWA_MAXIDX + 1) / (F * ff2); - A *= d; - B *= d; - C *= d; - - U0 = fx; - V0 = fy; - u1 = (int)(floorf(U0 - ue)); - u2 = (int)(ceilf(U0 + ue)); - v1 = (int)(floorf(V0 - ve)); - v2 = (int)(ceilf(V0 + ve)); - U0 -= 0.5f; - V0 -= 0.5f; - DDQ = 2.f * A; - U = u1 - U0; - ac1 = A * (2.f * U + 1.f); - ac2 = A * U * U; - BU = B * U; - - d = result[0] = result[1] = result[2] = result[3] = 0.f; - for (v = v1; v <= v2; ++v) { - const float V = v - V0; - float DQ = ac1 + B * V; - float Q = (C * V + BU) * V + ac2; - for (u = u1; u <= u2; ++u) { - if (Q < (float)(EWA_MAXIDX + 1)) { + float u = uv[0], v = uv[1]; + float Ux = derivatives[0][0], Vx = derivatives[1][0], Uy = derivatives[0][1], Vy = derivatives[1][1]; + float A, B, C, F, ue, ve; + ellipse_params(Ux, Uy, Vx, Vy, A, B, C, F, ue, ve); + + /* Note: highly eccentric ellipses can lead to large texture space areas to filter! + * This is limited somewhat by the EWA_WTS size in the loop, but a nicer approach + * could be the one found in + * "High Quality Elliptical Texture Filtering on GPU" + * by Pavlos Mavridis and Georgios Papaioannou + * in which the eccentricity of the ellipse is clamped. + */ + + int U0 = (int)u; + int V0 = (int)v; + /* pixel offset for interpolation */ + float ufac = u - floor(u), vfac = v - floor(v); + /* filter size */ + int u1 = (int)(u - ue); + int u2 = (int)(u + ue); + int v1 = (int)(v - ve); + int v2 = (int)(v + ve); + + /* sane clamping to avoid unnecessarily huge loops */ + /* note: if eccentricity gets clamped (see above), + * the ue/ve limits can also be lowered accordingly + */ + if (U0-u1 > EWA_MAXIDX) u1 = U0 - EWA_MAXIDX; + if (u2-U0 > EWA_MAXIDX) u2 = U0 + EWA_MAXIDX; + if (V0-v1 > EWA_MAXIDX) v1 = V0 - EWA_MAXIDX; + if (v2-V0 > EWA_MAXIDX) v2 = V0 + EWA_MAXIDX; + + float DDQ = 2.f * A; + float U = u1 - U0; + float ac1 = A * (2.f*U + 1.f); + float ac2 = A * U*U; + float BU = B * U; + + float sum = 0.0f; + for (int v = v1; v <= v2; ++v) { + float V = v - V0; + + float DQ = ac1 + B*V; + float Q = (C*V + BU)*V + ac2; + for (int u = u1; u <= u2; ++u) { + if (Q < F) { float tc[4]; - const float wt = EWA_WTS[(Q < 0.f) ? 0 : (unsigned int)Q]; - read(tc, clipuv(u, width), clipuv(v, height)); + const float wt = EWA_WTS[CLAMPIS((int)Q, 0, EWA_MAXIDX)]; + switch (sampler) { + case COM_PS_NEAREST: read(tc, u, v); break; + case COM_PS_BILINEAR: readBilinear(tc, (float)u+ufac, (float)v+vfac); break; + case COM_PS_BICUBIC: readBilinear(tc, (float)u+ufac, (float)v+vfac); break; /* XXX no readBicubic method yet */ + default: zero_v4(tc); break; + } madd_v4_v4fl(result, tc, wt); - d += wt; + sum += wt; } Q += DQ; DQ += DDQ; } } - // d should hopefully never be zero anymore - d = 1.f / d; - mul_v4_fl(result, d); + mul_v4_fl(result, (sum != 0.0f ? 1.0f / sum : 0.0f)); } diff --git a/source/blender/compositor/intern/COM_MemoryBuffer.h b/source/blender/compositor/intern/COM_MemoryBuffer.h index 548744f6660..f9cde79c07f 100644 --- a/source/blender/compositor/intern/COM_MemoryBuffer.h +++ b/source/blender/compositor/intern/COM_MemoryBuffer.h @@ -245,7 +245,7 @@ public: result[3] = color1[3] * mvaluex + color3[3] * valuex; } - void readEWA(float result[4], float fx, float fy, float dx, float dy, PixelSampler sampler); + void readEWA(float result[4], const float uv[2], const float derivatives[2][2], PixelSampler sampler); /** * @brief is this MemoryBuffer a temporarily buffer (based on an area, not on a chunk) diff --git a/source/blender/compositor/intern/COM_SocketReader.h b/source/blender/compositor/intern/COM_SocketReader.h index 2168611d0c0..c996ef5bbeb 100644 --- a/source/blender/compositor/intern/COM_SocketReader.h +++ b/source/blender/compositor/intern/COM_SocketReader.h @@ -88,7 +88,7 @@ protected: * @param dy * @param inputBuffers chunks that can be read by their ReadBufferOperation. */ - virtual void executePixelFiltered(float output[4], float x, float y, float dx, float dy, PixelSampler sampler) {} + virtual void executePixelFiltered(float output[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler) {} public: inline void readSampled(float result[4], float x, float y, PixelSampler sampler) { @@ -97,7 +97,7 @@ public: inline void read(float result[4], int x, int y, void *chunkData) { executePixel(result, x, y, chunkData); } - inline void readFiltered(float result[4], float x, float y, float dx, float dy, PixelSampler sampler) { + inline void readFiltered(float result[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler) { executePixelFiltered(result, x, y, dx, dy, sampler); } diff --git a/source/blender/compositor/operations/COM_DisplaceOperation.cpp b/source/blender/compositor/operations/COM_DisplaceOperation.cpp index 17692bd18ef..96031e7acc4 100644 --- a/source/blender/compositor/operations/COM_DisplaceOperation.cpp +++ b/source/blender/compositor/operations/COM_DisplaceOperation.cpp @@ -49,54 +49,92 @@ void DisplaceOperation::initExecution() this->m_height_x4 = this->getHeight() * 4; } - -/* minimum distance (in pixels) a pixel has to be displaced - * in order to take effect */ -#define DISPLACE_EPSILON 0.01f - -void DisplaceOperation::executePixel(float output[4], int x, int y, void *data) +void DisplaceOperation::executePixelSampled(float output[4], float x, float y, PixelSampler sampler) { - float inVector[4]; - float inScale[4]; + float xy[2] = { x, y }; + float uv[2], deriv[2][2]; - float p_dx, p_dy; /* main displacement in pixel space */ - float d_dx, d_dy; - float dxt, dyt; - float u, v; + pixelTransform(xy, uv, deriv); - this->m_inputScaleXProgram->readSampled(inScale, x, y, COM_PS_NEAREST); - float xs = inScale[0]; - this->m_inputScaleYProgram->readSampled(inScale, x, y, COM_PS_NEAREST); - float ys = inScale[0]; + /* EWA filtering (without nearest it gets blurry with NO distortion) */ + this->m_inputColorProgram->readFiltered(output, uv[0], uv[1], deriv[0], deriv[1], COM_PS_BILINEAR); +} +bool DisplaceOperation::read_displacement(float x, float y, float xscale, float yscale, const float origin[2], float &r_u, float &r_v) +{ + float width = m_inputVectorProgram->getWidth(); + float height = m_inputVectorProgram->getHeight(); + if (x < 0.0f || x >= width || y < 0.0f || y >= height) { + r_u = 0.0f; + r_v = 0.0f; + return false; + } + else { + float col[4]; + m_inputVectorProgram->readSampled(col, x, y, COM_PS_BILINEAR); + r_u = origin[0] - col[0]*xscale + 0.5f; + r_v = origin[1] - col[1]*yscale + 0.5f; + return true; + } +} + +void DisplaceOperation::pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2]) +{ + float col[4]; + float uv[2]; /* temporary variables for derivative estimation */ + int num; + + m_inputScaleXProgram->readSampled(col, xy[0], xy[1], COM_PS_NEAREST); + float xs = col[0]; + m_inputScaleYProgram->readSampled(col, xy[0], xy[1], COM_PS_NEAREST); + float ys = col[0]; /* clamp x and y displacement to triple image resolution - * to prevent hangs from huge values mistakenly plugged in eg. z buffers */ - CLAMP(xs, -this->m_width_x4, this->m_width_x4); - CLAMP(ys, -this->m_height_x4, this->m_height_x4); - - this->m_inputVectorProgram->readSampled(inVector, x, y, COM_PS_NEAREST); - p_dx = inVector[0] * xs; - p_dy = inVector[1] * ys; + CLAMP(xs, -m_width_x4, m_width_x4); + CLAMP(ys, -m_height_x4, m_height_x4); /* displaced pixel in uv coords, for image sampling */ - u = x - p_dx + 0.5f; - v = y - p_dy + 0.5f; - - /* calc derivatives */ - this->m_inputVectorProgram->readSampled(inVector, x + 1, y, COM_PS_NEAREST); - d_dx = inVector[0] * xs; - this->m_inputVectorProgram->readSampled(inVector, x, y + 1, COM_PS_NEAREST); - d_dy = inVector[1] * ys; - - /* clamp derivatives to minimum displacement distance in UV space */ - dxt = p_dx - d_dx; - dyt = p_dy - d_dy; - - dxt = signf(dxt) * max_ff(fabsf(dxt), DISPLACE_EPSILON) / this->getWidth(); - dyt = signf(dyt) * max_ff(fabsf(dyt), DISPLACE_EPSILON) / this->getHeight(); + read_displacement(xy[0], xy[1], xs, ys, xy, r_uv[0], r_uv[1]); + + /* Estimate partial derivatives using 1-pixel offsets */ + const float epsilon[2] = { 1.0f, 1.0f }; + + zero_v2(r_deriv[0]); + zero_v2(r_deriv[1]); + + num = 0; + if (read_displacement(xy[0] + epsilon[0], xy[1], xs, ys, xy, uv[0], uv[1])) { + r_deriv[0][0] += uv[0] - r_uv[0]; + r_deriv[1][0] += uv[1] - r_uv[1]; + ++num; + } + if (read_displacement(xy[0] - epsilon[0], xy[1], xs, ys, xy, uv[0], uv[1])) { + r_deriv[0][0] += r_uv[0] - uv[0]; + r_deriv[1][0] += r_uv[1] - uv[1]; + ++num; + } + if (num > 0) { + float numinv = 1.0f / (float)num; + r_deriv[0][0] *= numinv; + r_deriv[1][0] *= numinv; + } - /* EWA filtering (without nearest it gets blurry with NO distortion) */ - this->m_inputColorProgram->readFiltered(output, u, v, dxt, dyt, COM_PS_NEAREST); + num = 0; + if (read_displacement(xy[0], xy[1] + epsilon[1], xs, ys, xy, uv[0], uv[1])) { + r_deriv[0][1] += uv[0] - r_uv[0]; + r_deriv[1][1] += uv[1] - r_uv[1]; + ++num; + } + if (read_displacement(xy[0], xy[1] - epsilon[1], xs, ys, xy, uv[0], uv[1])) { + r_deriv[0][1] += r_uv[0] - uv[0]; + r_deriv[1][1] += r_uv[1] - uv[1]; + ++num; + } + if (num > 0) { + float numinv = 1.0f / (float)num; + r_deriv[0][1] *= numinv; + r_deriv[1][1] *= numinv; + } } void DisplaceOperation::deinitExecution() @@ -126,10 +164,10 @@ bool DisplaceOperation::determineDependingAreaOfInterest(rcti *input, ReadBuffer /* vector */ operation = getInputOperation(1); - vectorInput.xmax = input->xmax + 2; - vectorInput.xmin = input->xmin; - vectorInput.ymax = input->ymax + 2; - vectorInput.ymin = input->ymin; + vectorInput.xmax = input->xmax + 1; + vectorInput.xmin = input->xmin - 1; + vectorInput.ymax = input->ymax + 1; + vectorInput.ymin = input->ymin - 1; if (operation->determineDependingAreaOfInterest(&vectorInput, readOperation, output)) { return true; } diff --git a/source/blender/compositor/operations/COM_DisplaceOperation.h b/source/blender/compositor/operations/COM_DisplaceOperation.h index daf54294aa1..cec7937d9d6 100644 --- a/source/blender/compositor/operations/COM_DisplaceOperation.h +++ b/source/blender/compositor/operations/COM_DisplaceOperation.h @@ -48,7 +48,9 @@ public: /** * the inner loop of this program */ - void executePixel(float output[4], int x, int y, void *data); + void executePixelSampled(float output[4], float x, float y, PixelSampler sampler); + + void pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2]); /** * Initialize the execution @@ -59,5 +61,8 @@ public: * Deinitialize the execution */ void deinitExecution(); + +private: + bool read_displacement(float x, float y, float xscale, float yscale, const float origin[2], float &r_u, float &r_v); }; #endif diff --git a/source/blender/compositor/operations/COM_MapUVOperation.cpp b/source/blender/compositor/operations/COM_MapUVOperation.cpp index 289b447dec7..292f073548a 100644 --- a/source/blender/compositor/operations/COM_MapUVOperation.cpp +++ b/source/blender/compositor/operations/COM_MapUVOperation.cpp @@ -42,46 +42,28 @@ void MapUVOperation::initExecution() void MapUVOperation::executePixelSampled(float output[4], float x, float y, PixelSampler sampler) { - float inputUV[4]; - float uv_a[4], uv_b[4]; - float u, v; + float xy[2] = { x, y }; + float uv[2], deriv[2][2], alpha; - float dx, dy; - float uv_l, uv_r; - float uv_u, uv_d; - - this->m_inputUVProgram->readSampled(inputUV, x, y, sampler); - if (inputUV[2] == 0.f) { + pixelTransform(xy, uv, deriv, alpha); + if (alpha == 0.0f) { zero_v4(output); return; } - /* adaptive sampling, red (U) channel */ - this->m_inputUVProgram->readSampled(uv_a, x - 1, y, COM_PS_NEAREST); - this->m_inputUVProgram->readSampled(uv_b, x + 1, y, COM_PS_NEAREST); - uv_l = uv_a[2] != 0.f ? fabsf(inputUV[0] - uv_a[0]) : 0.f; - uv_r = uv_b[2] != 0.f ? fabsf(inputUV[0] - uv_b[0]) : 0.f; - - dx = 0.5f * (uv_l + uv_r); - - /* adaptive sampling, green (V) channel */ - this->m_inputUVProgram->readSampled(uv_a, x, y - 1, COM_PS_NEAREST); - this->m_inputUVProgram->readSampled(uv_b, x, y + 1, COM_PS_NEAREST); - uv_u = uv_a[2] != 0.f ? fabsf(inputUV[1] - uv_a[1]) : 0.f; - uv_d = uv_b[2] != 0.f ? fabsf(inputUV[1] - uv_b[1]) : 0.f; - - dy = 0.5f * (uv_u + uv_d); + /* EWA filtering */ + this->m_inputColorProgram->readFiltered(output, uv[0], uv[1], deriv[0], deriv[1], COM_PS_BILINEAR); + /* UV to alpha threshold */ const float threshold = this->m_alpha * 0.05f; - float alpha = 1.0f - threshold * (dx + dy); - if (alpha < 0.f) alpha = 0.f; - else alpha *= inputUV[2]; - - /* EWA filtering */ - u = inputUV[0] * this->m_inputColorProgram->getWidth(); - v = inputUV[1] * this->m_inputColorProgram->getHeight(); - - this->m_inputColorProgram->readFiltered(output, u, v, dx, dy, COM_PS_NEAREST); + /* XXX alpha threshold is used to fade out pixels on boundaries with invalid derivatives. + * this calculation is not very well defined, should be looked into if it becomes a problem ... + */ + float du = len_v2(deriv[0]); + float dv = len_v2(deriv[1]); + float factor = 1.0f - threshold * (du + dv); + if (factor < 0.f) alpha = 0.f; + else alpha *= factor; /* "premul" */ if (alpha < 1.0f) { @@ -89,6 +71,74 @@ void MapUVOperation::executePixelSampled(float output[4], float x, float y, Pixe } } +bool MapUVOperation::read_uv(float x, float y, float &r_u, float &r_v, float &r_alpha) +{ + float width = m_inputUVProgram->getWidth(); + float height = m_inputUVProgram->getHeight(); + if (x < 0.0f || x >= width || y < 0.0f || y >= height) { + r_u = 0.0f; + r_v = 0.0f; + r_alpha = 0.0f; + return false; + } + else { + float col[4]; + m_inputUVProgram->readSampled(col, x, y, COM_PS_BILINEAR); + r_u = col[0] * width; + r_v = col[1] * height; + r_alpha = col[2]; + return true; + } +} + +void MapUVOperation::pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2], float &r_alpha) +{ + float uv[2], alpha; /* temporary variables for derivative estimation */ + int num; + + read_uv(xy[0], xy[1], r_uv[0], r_uv[1], r_alpha); + + /* Estimate partial derivatives using 1-pixel offsets */ + const float epsilon[2] = { 1.0f, 1.0f }; + + zero_v2(r_deriv[0]); + zero_v2(r_deriv[1]); + + num = 0; + if (read_uv(xy[0] + epsilon[0], xy[1], uv[0], uv[1], alpha)) { + r_deriv[0][0] += uv[0] - r_uv[0]; + r_deriv[1][0] += uv[1] - r_uv[1]; + ++num; + } + if (read_uv(xy[0] - epsilon[0], xy[1], uv[0], uv[1], alpha)) { + r_deriv[0][0] += r_uv[0] - uv[0]; + r_deriv[1][0] += r_uv[1] - uv[1]; + ++num; + } + if (num > 0) { + float numinv = 1.0f / (float)num; + r_deriv[0][0] *= numinv; + r_deriv[1][0] *= numinv; + } + + num = 0; + if (read_uv(xy[0], xy[1] + epsilon[1], uv[0], uv[1], alpha)) { + r_deriv[0][1] += uv[0] - r_uv[0]; + r_deriv[1][1] += uv[1] - r_uv[1]; + ++num; + } + if (read_uv(xy[0], xy[1] - epsilon[1], uv[0], uv[1], alpha)) { + r_deriv[0][1] += r_uv[0] - uv[0]; + r_deriv[1][1] += r_uv[1] - uv[1]; + ++num; + } + if (num > 0) { + float numinv = 1.0f / (float)num; + r_deriv[0][1] *= numinv; + r_deriv[1][1] *= numinv; + } +} + void MapUVOperation::deinitExecution() { this->m_inputUVProgram = NULL; diff --git a/source/blender/compositor/operations/COM_MapUVOperation.h b/source/blender/compositor/operations/COM_MapUVOperation.h index fe8bfd2a9ac..796ee952607 100644 --- a/source/blender/compositor/operations/COM_MapUVOperation.h +++ b/source/blender/compositor/operations/COM_MapUVOperation.h @@ -46,7 +46,9 @@ public: * the inner loop of this program */ void executePixelSampled(float output[4], float x, float y, PixelSampler sampler); - + + void pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2], float &r_alpha); + /** * Initialize the execution */ @@ -58,5 +60,9 @@ public: void deinitExecution(); void setAlpha(float alpha) { this->m_alpha = alpha; } + +private: + bool read_uv(float x, float y, float &r_u, float &r_v, float &r_alpha); }; + #endif diff --git a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp index dba3a6f3505..7780023c465 100644 --- a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp +++ b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.cpp @@ -36,60 +36,17 @@ extern "C" { #include "BKE_tracking.h" } -BLI_INLINE bool isPointInsideQuad(const float x, const float y, const float corners[4][2]) -{ - float point[2]; - - point[0] = x; - point[1] = y; - - return isect_point_tri_v2(point, corners[0], corners[1], corners[2]) || - isect_point_tri_v2(point, corners[0], corners[2], corners[3]); -} - -BLI_INLINE void warpCoord(float x, float y, float matrix[3][3], float uv[2]) +BLI_INLINE void warpCoord(float x, float y, float matrix[3][3], float uv[2], float deriv[2][2]) { float vec[3] = {x, y, 1.0f}; mul_m3_v3(matrix, vec); - vec[0] /= vec[2]; - vec[1] /= vec[2]; + uv[0] = vec[0] / vec[2]; + uv[1] = vec[1] / vec[2]; - copy_v2_v2(uv, vec); -} - -BLI_INLINE void resolveUVAndDxDy(const float x, const float y, float matrix[3][3], - float *u_r, float *v_r, float *dx_r, float *dy_r) -{ - float inputUV[2]; - float uv_a[2], uv_b[2]; - - float dx, dy; - float uv_l, uv_r; - float uv_u, uv_d; - - warpCoord(x, y, matrix, inputUV); - - /* adaptive sampling, red (U) channel */ - warpCoord(x - 1, y, matrix, uv_a); - warpCoord(x + 1, y, matrix, uv_b); - uv_l = fabsf(inputUV[0] - uv_a[0]); - uv_r = fabsf(inputUV[0] - uv_b[0]); - - dx = 0.5f * (uv_l + uv_r); - - /* adaptive sampling, green (V) channel */ - warpCoord(x, y - 1, matrix, uv_a); - warpCoord(x, y + 1, matrix, uv_b); - uv_u = fabsf(inputUV[1] - uv_a[1]); - uv_d = fabsf(inputUV[1] - uv_b[1]); - - dy = 0.5f * (uv_u + uv_d); - - *dx_r = dx; - *dy_r = dy; - - *u_r = inputUV[0]; - *v_r = inputUV[1]; + deriv[0][0] = (matrix[0][0] - matrix[0][2] * uv[0]) / vec[2]; + deriv[1][0] = (matrix[0][1] - matrix[0][2] * uv[1]) / vec[2]; + deriv[0][1] = (matrix[1][0] - matrix[1][2] * uv[0]) / vec[2]; + deriv[1][1] = (matrix[1][1] - matrix[1][2] * uv[1]) / vec[2]; } PlaneTrackWarpImageOperation::PlaneTrackWarpImageOperation() : PlaneTrackCommonOperation() @@ -98,9 +55,6 @@ PlaneTrackWarpImageOperation::PlaneTrackWarpImageOperation() : PlaneTrackCommonO this->addOutputSocket(COM_DT_COLOR); this->m_pixelReader = NULL; this->setComplex(true); - - /* Currently hardcoded to 8 samples. */ - this->m_osa = 8; } void PlaneTrackWarpImageOperation::initExecution() @@ -109,8 +63,6 @@ void PlaneTrackWarpImageOperation::initExecution() this->m_pixelReader = this->getInputSocketReader(0); - BLI_jitter_init(this->m_jitter[0], this->m_osa); - const int width = this->m_pixelReader->getWidth(); const int height = this->m_pixelReader->getHeight(); float frame_corners[4][2] = {{0.0f, 0.0f}, @@ -127,41 +79,32 @@ void PlaneTrackWarpImageOperation::deinitExecution() this->m_pixelReader = NULL; } -void PlaneTrackWarpImageOperation::executePixelSampled(float output[4], float x_, float y_, PixelSampler sampler) +void PlaneTrackWarpImageOperation::executePixelSampled(float output[4], float x, float y, PixelSampler sampler) { - float color_accum[4]; + float xy[2] = {x, y}; + float uv[2]; + float deriv[2][2]; - zero_v4(color_accum); - for (int sample = 0; sample < this->m_osa; sample++) { - float current_x = x_ + this->m_jitter[sample][0], - current_y = y_ + this->m_jitter[sample][1]; - if (isPointInsideQuad(current_x, current_y, this->m_frameSpaceCorners)) { - float current_color[4]; - float u, v, dx, dy; + pixelTransform(xy, uv, deriv); - resolveUVAndDxDy(current_x, current_y, m_perspectiveMatrix, &u, &v, &dx, &dy); - - /* derivatives are to be in normalized space.. */ - dx /= this->m_pixelReader->getWidth(); - dy /= this->m_pixelReader->getHeight(); - - this->m_pixelReader->readFiltered(current_color, u, v, dx, dy, COM_PS_NEAREST); - add_v4_v4(color_accum, current_color); - } - } + m_pixelReader->readFiltered(output, uv[0], uv[1], deriv[0], deriv[1], COM_PS_BILINEAR); +} - mul_v4_v4fl(output, color_accum, 1.0f / this->m_osa); +void PlaneTrackWarpImageOperation::pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2]) +{ + warpCoord(xy[0], xy[1], m_perspectiveMatrix, r_uv, r_deriv); } bool PlaneTrackWarpImageOperation::determineDependingAreaOfInterest(rcti *input, ReadBufferOperation *readOperation, rcti *output) { float UVs[4][2]; + float deriv[2][2]; /* TODO(sergey): figure out proper way to do this. */ - warpCoord(input->xmin - 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[0]); - warpCoord(input->xmax + 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[1]); - warpCoord(input->xmax + 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[2]); - warpCoord(input->xmin - 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[3]); + warpCoord(input->xmin - 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[0], deriv); + warpCoord(input->xmax + 2, input->ymin - 2, this->m_perspectiveMatrix, UVs[1], deriv); + warpCoord(input->xmax + 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[2], deriv); + warpCoord(input->xmin - 2, input->ymax + 2, this->m_perspectiveMatrix, UVs[3], deriv); float min[2], max[2]; INIT_MINMAX2(min, max); diff --git a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h index ed9cc2fe5aa..ceb002c8039 100644 --- a/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h +++ b/source/blender/compositor/operations/COM_PlaneTrackWarpImageOperation.h @@ -36,8 +36,6 @@ class PlaneTrackWarpImageOperation : public PlaneTrackCommonOperation { protected: SocketReader *m_pixelReader; - int m_osa; - float m_jitter[32][2]; float m_perspectiveMatrix[3][3]; public: @@ -47,6 +45,7 @@ public: void deinitExecution(); void executePixelSampled(float output[4], float x, float y, PixelSampler sampler); + void pixelTransform(const float xy[2], float r_uv[2], float r_deriv[2][2]); bool determineDependingAreaOfInterest(rcti *input, ReadBufferOperation *readOperation, rcti *output); }; diff --git a/source/blender/compositor/operations/COM_ReadBufferOperation.cpp b/source/blender/compositor/operations/COM_ReadBufferOperation.cpp index b93e338ad6c..47d5fc6bcca 100644 --- a/source/blender/compositor/operations/COM_ReadBufferOperation.cpp +++ b/source/blender/compositor/operations/COM_ReadBufferOperation.cpp @@ -81,14 +81,16 @@ void ReadBufferOperation::executePixelExtend(float output[4], float x, float y, } } -void ReadBufferOperation::executePixelFiltered(float output[4], float x, float y, float dx, float dy, PixelSampler sampler) +void ReadBufferOperation::executePixelFiltered(float output[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler) { if (m_single_value) { /* write buffer has a single value stored at (0,0) */ m_buffer->read(output, 0, 0); } else { - m_buffer->readEWA(output, x, y, dx, dy, sampler); + const float uv[2] = { x, y }; + const float deriv[2][2] = { {dx[0], dx[1]}, {dy[0], dy[1]} }; + m_buffer->readEWA(output, uv, deriv, sampler); } } diff --git a/source/blender/compositor/operations/COM_ReadBufferOperation.h b/source/blender/compositor/operations/COM_ReadBufferOperation.h index cdc6e50f6d4..cce519c6eb3 100644 --- a/source/blender/compositor/operations/COM_ReadBufferOperation.h +++ b/source/blender/compositor/operations/COM_ReadBufferOperation.h @@ -43,7 +43,7 @@ public: void executePixelSampled(float output[4], float x, float y, PixelSampler sampler); void executePixelExtend(float output[4], float x, float y, PixelSampler sampler, MemoryBufferExtend extend_x, MemoryBufferExtend extend_y); - void executePixelFiltered(float output[4], float x, float y, float dx, float dy, PixelSampler sampler); + void executePixelFiltered(float output[4], float x, float y, float dx[2], float dy[2], PixelSampler sampler); const bool isReadBufferOperation() const { return true; } void setOffset(unsigned int offset) { this->m_offset = offset; } unsigned int getOffset() const { return this->m_offset; } |