From 06b0562154c130ed56533adf5dc51115e6f2594e Mon Sep 17 00:00:00 2001 From: Wenzel Jakob Date: Mon, 22 Oct 2012 19:30:47 -0400 Subject: [PATCH] new branch to implement displacements based on work by Milos. Step 1: heightfields! --- src/shapes/SConscript | 1 + src/shapes/heightfield.cpp | 371 +++++++++++++++++++++++++++++++++++++ 2 files changed, 372 insertions(+) create mode 100644 src/shapes/heightfield.cpp diff --git a/src/shapes/SConscript b/src/shapes/SConscript index 6abc53ce..4c48f497 100644 --- a/src/shapes/SConscript +++ b/src/shapes/SConscript @@ -12,5 +12,6 @@ plugins += env.SharedLibrary('hair', ['hair.cpp']) plugins += env.SharedLibrary('shapegroup', ['shapegroup.cpp']) plugins += env.SharedLibrary('instance', ['instance.cpp']) plugins += env.SharedLibrary('animatedinstance', ['animatedinstance.cpp']) +plugins += env.SharedLibrary('heightfield', ['heightfield.cpp']) Export('plugins') diff --git a/src/shapes/heightfield.cpp b/src/shapes/heightfield.cpp new file mode 100644 index 00000000..5d91aa92 --- /dev/null +++ b/src/shapes/heightfield.cpp @@ -0,0 +1,371 @@ +/* + This file is part of Mitsuba, a physically based rendering system. + + Copyright (c) 2007-2011 by Wenzel Jakob and others. + + Mitsuba is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License Version 3 + as published by the Free Software Foundation. + + Mitsuba 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 this program. If not, see . +*/ + +#include +#include +#include +#include +#include +#include + +MTS_NAMESPACE_BEGIN + +/*!\plugin{heightfield}{Height field} + * \order{12} + * + * Developed by Milo\^{s} Ha\^{s}an + */ +class Heightfield : public Shape { +private: + inline Float getIJ(int i, int j) const { + i = std::max(0, std::min(m_rows-1, i)); + j = std::max(0, std::min(m_cols-1, j)); + return m_data[i*m_cols + j]; + } + + inline Normal getNormalIJ(int i, int j) const { + i = std::max(0, std::min(m_rows-1, i)); + j = std::max(0, std::min(m_cols-1, j)); + return m_normals[i*m_cols + j]; + } + + inline Float getXY(int x, int y) const { + return getIJ(m_rows - y - 1, x); + } + + inline Normal getNormalXY(int x, int y) const { + return getNormalIJ(m_rows - y - 1, x); + } + + inline Float getBilinear(Float x, Float y) const { + // cout << "getBilinear: " << x << ' ' << y << ' '; + int xi = (int) std::floor(x); + int yi = (int) std::floor(y); + Float u = x - xi, v = y - yi; + Float z00 = getXY(xi, yi); + Float z01 = getXY(xi, yi+1); + Float z10 = getXY(xi+1, yi); + Float z11 = getXY(xi+1, yi+1); + Float result = (1-u)*(1-v)*z00 + (1-u)*v*z01 + u*(1-v)*z10 + u*v*z11; + // cout << result << '\n'; + return result; + } + + inline Vector2 getBilinearGradient(Float x, Float y) const { + int xi = (int) std::floor(x); + int yi = (int) std::floor(y); + Float u = x - xi, v = y - yi; + Float z00 = getXY(xi, yi); + Float z01 = getXY(xi, yi+1); + Float z10 = getXY(xi+1, yi); + Float z11 = getXY(xi+1, yi+1); + Float du = (1-v) * (z10 - z00) + v * (z11 - z01); + Float dv = (1-u) * (z01 - z00) + u * (z11 - z10); + return Vector2(du, dv); + } + + inline Normal getInterpolatedNormal(Float x, Float y) const { + int xi = (int) std::floor(x); + int yi = (int) std::floor(y); + Float u = x - xi, v = y - yi; + Normal z00 = getNormalXY(xi, yi); + Normal z01 = getNormalXY(xi, yi+1); + Normal z10 = getNormalXY(xi+1, yi); + Normal z11 = getNormalXY(xi+1, yi+1); + Normal result = (1-u)*(1-v)*z00 + (1-u)*v*z01 + u*(1-v)*z10 + u*v*z11; + return normalize(result); + } + +public: + Heightfield(const Properties &props) : Shape(props) { + m_pixelSize = props.getFloat("pixelSize", 0.01f); + Float multiplier = props.getFloat("multiplier", 1); + m_flipNormal = props.getBoolean("flipNormal", false); + + // read bitmap; use only red channel + fs::path hfPath = props.getString("hf"); + ref stream = new FileStream(hfPath); + ref hf = new Bitmap(Bitmap::EOpenEXR, stream); + m_rows = hf->getHeight(); + m_cols = hf->getWidth(); + Vector4* hfData = (Vector4*) hf->getFloatData(); + + m_data = new float[m_rows * m_cols]; + for (int i = 0; i < m_rows * m_cols; i++) + m_data[i] = hfData[i].x * multiplier; + + // compute bounding box + Float xMax = (m_cols-1) * m_pixelSize / 2; + Float yMax = (m_rows-1) * m_pixelSize / 2; + Float zMin = std::numeric_limits::infinity(), zMax = -zMin; + for (int i = 0; i < m_rows * m_cols; i++) { + Float v = m_data[i]; + zMin = std::min(zMin, v); + zMax = std::max(zMax, v); + } + zMin -= Epsilon; zMax += Epsilon; + m_aabb = AABB(Point(-xMax, -yMax, zMin), Point(xMax, yMax, zMax)); + // cout << m_aabb.toString() << '\n'; + + // compute normals + m_normals = 0; + if (props.getBoolean("interpolateNormals", false)) { + m_normals = new Normal[m_rows * m_cols]; + int index = 0; + + for (int i = 0; i < m_rows; i++) + for (int j = 0; j < m_cols; j++) { + Float dx = getIJ(i, j+1) - getIJ(i, j-1); + Float dy = getIJ(i-1, j) - getIJ(i+1, j); // reversal between i and y + Normal n(-dx / m_pixelSize, -dy / m_pixelSize, 2); + m_normals[index++] = normalize(n); + } + } + } + + Heightfield(Stream *stream, InstanceManager *manager) + : Shape(stream, manager) { + Assert(false); + } + + ~Heightfield() { + delete[] m_data; + if (m_normals) delete[] m_normals; + } + + void serialize(Stream *stream, InstanceManager *manager) const { + Shape::serialize(stream, manager); + Assert(false); + } + + AABB getAABB() const { + return m_aabb; + } + + Float getSurfaceArea() const { + return 0; + } + + /// find smallest t >= 0 such that a*t + b is integer + inline static Float nextInt(Float a, Float b) { + if (a == 0) return std::numeric_limits::infinity(); + else if (a > 0) return (std::ceil(b) - b) / a; + else return (std::floor(b) - b) / a; + } + + /// clip ray against bounding box + bool clipRay(const Ray& ray, Float& mint, Float& maxt) const { + Float tmp1, tmp2; + if (!m_aabb.rayIntersect(ray, tmp1, tmp2)) return false; + mint = std::max(mint, tmp1); + maxt = std::min(maxt, tmp2); + if (maxt <= mint) return false; + return true; + } + + /// transform from world to grid space + Point world2grid(Point p) const { + p.x /= m_pixelSize; + p.y /= m_pixelSize; + p.x += (m_cols - 1) / 2.0f; + p.y += (m_rows - 1) / 2.0f; + return p; + } + + /// Subtract ray, fit quadratic function and solve. + /// x0, xMid and x1 assumed to be 0, 0.5 and 1. + /// Return -1 if no solution in [0,1]. + inline Float fitAndSolve(Float f0, Float fMid, Float f1, Float z0, Float z1) const { + // subtract the "z" line of the ray + f0 -= z0; f1 -= z1; fMid -= (z0 + z1) / 2; + + // fit quadratic function to values + Float a = 2*f0 - 4*fMid + 2*f1; + Float b = -3*f0 + 4*fMid - f1; + Float c = f0; + + // solve and return the smaller root within [0,1] (if any) + Float t0, t1; + if (!solveQuadratic(a, b, c, t0, t1)) return -1; + if (t0 >= 0 && t0 <= 1) return t0; + if (t1 >= 0 && t1 <= 1) return t1; + return -1; + } + + /// find nearest intersection of ray with bilinear interpolant + Float intersectPatch(Float x0, Float y0, Float z0, + Float x1, Float y1, Float z1) const { + Float f0 = getBilinear(x0, y0); + Float fMid = getBilinear((x0 + x1) / 2, (y0 + y1) / 2); + Float f1 = getBilinear(x1, y1); + return fitAndSolve(f0, fMid, f1, z0, z1); + } + + bool rayIntersect(const Ray &ray, Float mint, Float maxt, Float &tResult, void *tmp) const { + if (!clipRay(ray, mint, maxt)) return false; + + // ray segment endpoints in grid space + Point start = world2grid(ray(mint)); + Point end = world2grid(ray(maxt)); + + // the 2D vector from start to end of ray + Vector d = end - start; + if (d.lengthSquared() == 0) return false; + + // handle trivial case of perpendicular ray + if (d.x == 0 && d.y == 0) { + Assert(d.z != 0); + Float z = getBilinear(start.x, start.y); + + if ((z - start.z) * (z - end.z) > 0) return false; + + Float t = (z - start.z) / (end.z - start.z); + tResult = (1-t) * mint + t * maxt; + return true; + } + + // the 2D ray parameter and its deltas (can be infinite) + Float t = 0; + Float tMax = std::sqrt(d.x*d.x + d.y*d.y); + Float tDeltaX = std::abs(tMax / d.x); + Float tDeltaY = std::abs(tMax / d.y); + d /= tMax; + Float tMaxX = nextInt(d.x, start.x); + Float tMaxY = nextInt(d.y, start.y); + Float tNext; + + while (t < tMax) { + if (tMaxX < tMaxY) { + // advance in x + tNext = std::min(tMaxX, tMax); + tMaxX += tDeltaX; + } else { + // advance in y + tNext = std::min(tMaxY, tMax); + tMaxY += tDeltaY; + } + + // get intersection point (shifted into [0,1]) + Point p = start + t * d; + Point q = start + tNext * d; + Float tPatch = intersectPatch(p.x, p.y, p.z, q.x, q.y, q.z); + + if (tPatch >= 0) { + // found intersection, shift back from [0,1] + tPatch = t + tPatch * (tNext - t); + tResult = mint + tPatch * (maxt - mint) / tMax; + return true; + } + + t = tNext; + } + + return false; + } + + bool rayIntersect(const Ray &r, Float mint, Float maxt) const { + Float t; + return rayIntersect(r, mint, maxt, t, 0); + } + + void fillIntersectionRecord(const Ray &ray, + const void *temp, Intersection &its) const { + Float t = its.t; + memset(&its, 0, sizeof(Intersection)); + + // position + its.t = t; + its.p = ray(t); + + // normals + Point pGrid = world2grid(its.p); + Vector2 grad = getBilinearGradient(pGrid.x, pGrid.y); + Normal gn = normalize(Vector(-grad.x, -grad.y, m_pixelSize)); + Normal sn = m_normals ? getInterpolatedNormal(pGrid.x, pGrid.y) : gn; + + // primitive uv-s + its.uv = Point2(pGrid.x, pGrid.y); + its.dpdu = Vector(m_pixelSize, 0, grad.x); + its.dpdv = Vector(0, m_pixelSize, grad.y); + + if (m_flipNormal) { gn *= -1; sn *= -1; } + its.geoFrame = Frame(gn); + its.shFrame = Frame(sn); + + its.wi = its.toLocal(-ray.d); + its.shape = this; + } + + void getNormalDerivative(const Intersection &its, + Vector &dndu, Vector &dndv, bool shadingFrame) const { + if (!m_normals) { + dndu = Vector(0.0f); + dndv = Vector(0.0f); + return; + } + + Float x = its.uv.x, y = its.uv.y; + int xi = (int) std::floor(x); + int yi = (int) std::floor(y); + Float u = x - xi, v = y - yi; + Normal n00 = getNormalXY(xi, yi); + Normal n01 = getNormalXY(xi, yi+1); + Normal n10 = getNormalXY(xi+1, yi); + Normal n11 = getNormalXY(xi+1, yi+1); + + Normal N((1-u)*(1-v)*n00 + (1-u)*v*n01 + u*(1-v)*n10 + u*v*n11); + Float il = 1.0f / N.length(); N *= il; + + dndu = (1-v) * (n10 - n00) + v * (n11 - n01); + dndu *= il; dndu -= N * dot(N, dndu); + dndv = (1-u) * (n01 - n00) + u * (n11 - n10); + dndv *= il; dndv -= N * dot(N, dndv); + } + + size_t getPrimitiveCount() const { + return 1; + } + + size_t getEffectivePrimitiveCount() const { + return 1; + } + + std::string toString() const { + std::ostringstream oss; + oss << "Heightfield[" << endl + << " pixelSize = " << m_pixelSize << ", " << endl + << " rows = " << m_rows << ", " << endl + << " cols = " << m_cols << ", " << endl + << "]"; + return oss.str(); + } + + MTS_DECLARE_CLASS() +private: + AABB m_aabb; + Float m_pixelSize; + int m_rows, m_cols; + float* m_data; + Normal* m_normals; + bool m_flipNormal; +}; + +MTS_IMPLEMENT_CLASS_S(Heightfield, false, Shape) +MTS_EXPORT_PLUGIN(Heightfield, "Height field intersection primitive"); +MTS_NAMESPACE_END +