rfilter.h

Go to the documentation of this file.

/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2014 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 <http://www.gnu.org/licenses/>.
*/

#pragma once
#if !defined(__MITSUBA_CORE_RFILTER_H_)
#define __MITSUBA_CORE_RFILTER_H_

#include <mitsuba/core/cobject.h>

MTS_NAMESPACE_BEGIN

/// Reconstruction filters will be tabulated at this resolution
#define MTS_FILTER_RESOLUTION 31

/**
 * \brief Generic interface to separable image reconstruction filters
 *
 * When resampling bitmaps or adding radiance-valued samples to a rendering in
 * progress, Mitsuba first convolves them with a so-called image reconstruction
 * filter. Various kinds are implemented as subclasses of this interface.
 *
 * Because image filters are generally too expensive to evaluate for
 * each sample, the implementation of this class internally precomputes
 * an discrete representation (resolution given by \ref MTS_FILTER_RESOLUTION)
 *
 * \ingroup librender
 * \ingroup libpython
 */
class MTS_EXPORT_CORE ReconstructionFilter : public ConfigurableObject {
public:
        /**
         * \brief When resampling data to a different resolution using
         * \ref Resampler::resample(), this enumeration specifies how lookups
         * <em>outside</em> of the input domain are handled.
         *
         * \see Resampler
         */
        enum EBoundaryCondition {
                /// Clamp to the outermost sample position
                EClamp = 0,
                /// Assume that the input repeats in a periodic fashion
                ERepeat,
                /// Assume that the input is mirrored along the boundary
                EMirror,
                /// Assume that the input function is zero outside of the defined domain
                EZero,
                /// Assume that the input function is equal to one outside of the defined domain
                EOne
        };

        /// Return the filter's width
        inline Float getRadius() const { return m_radius; }

        /// Return the block border size required when rendering with this filter
        inline int getBorderSize() const { return m_borderSize; }

        /// Evaluate the filter function
        virtual Float eval(Float x) const = 0;

        /// Perform a lookup into the discretized version
        inline Float evalDiscretized(Float x) const { return m_values[
                std::min((int) std::abs(x * m_scaleFactor), MTS_FILTER_RESOLUTION)]; }

        /// Serialize the filter to a binary data stream
        void serialize(Stream *stream, InstanceManager *manager) const;

        /** \brief Configure the object (called \a once after construction) */
        void configure();

        MTS_DECLARE_CLASS()
protected:
        /// Create a new reconstruction filter
        ReconstructionFilter(const Properties &props);

        /// Unserialize a filter
        ReconstructionFilter(Stream *stream, InstanceManager *manager);

        /// Virtual destructor
        virtual ~ReconstructionFilter();
protected:
        Float m_radius, m_scaleFactor;
        Float m_values[MTS_FILTER_RESOLUTION+1];
        int m_borderSize;
};

/**
 * \brief Utility class for efficiently resampling discrete datasets to different resolutions
 * \tparam Scalar
 *      Denotes the underlying floating point data type (i.e. <tt>half</tt>, <tt>float</tt>,
 *      or <tt>double</tt>)
 */
template <typename Scalar> struct Resampler {
        /**
         * \brief Create a new Resampler object that transforms between the specified resolutions
         *
         * This constructor precomputes all information needed to efficiently perform the
         * desired resampling operation. For that reason, it is most efficient if it can
         * be used over and over again (e.g. to resample the equal-sized rows of a bitmap)
         *
         * \param sourceRes
         *      Source resolution
         * \param targetRes
         *      Desired target resolution
         * \param bc
         *      Boundary conditions that should be observed when looking up samples
         *      outside of the defined input domain.
         */
        Resampler(const ReconstructionFilter *rfilter, ReconstructionFilter::EBoundaryCondition bc,
                        int sourceRes, int targetRes) : m_bc(bc), m_sourceRes(sourceRes), m_targetRes(targetRes),
                        m_start(NULL), m_weights(NULL) {
                SAssert(sourceRes > 0 && targetRes > 0);
                Float filterRadius = rfilter->getRadius(), scale = (Float)1.0, invScale = (Float)1.0;

                /* Low-pass filter: scale reconstruction filters when downsampling */
                if (targetRes < sourceRes) {
                        scale = (Float) sourceRes / (Float) targetRes;
                        invScale = 1 / scale;
                        filterRadius *= scale;
                }

                m_taps = math::ceilToInt(filterRadius * 2);
                if (sourceRes == targetRes && (m_taps % 2) != 1)
                        --m_taps;
                m_halfTaps = m_taps / 2;

                if (sourceRes != targetRes) { /* Resampling mode */
                        m_start = new int[targetRes];
                        m_weights = new Scalar[m_taps * targetRes];
                        m_fastStart = 0;
                        m_fastEnd = m_targetRes;

                        for (int i=0; i<targetRes; i++) {
                                /* Compute the fractional coordinates of the new sample i in the original coordinates */
                                Float center = (i + (Float) 0.5f) / targetRes * sourceRes;

                                /* Determine the index of the first original sample that might contribute */
                                m_start[i] = math::floorToInt(center - filterRadius + (Float) 0.5f);

                                /* Determine the size of center region, on which to run fast non condition-aware code */
                                if (m_start[i] < 0)
                                        m_fastStart = std::max(m_fastStart, i + 1);
                                else if (m_start[i] + m_taps - 1 >= m_sourceRes)
                                        m_fastEnd = std::min(m_fastEnd, i - 1);

                                Float sum = 0;
                                for (int j=0; j<m_taps; j++) {
                                        /* Compute the the position where the filter should be evaluated */
                                        Float pos = m_start[i] + j + (Float) 0.5f - center;

                                        /* Perform the evaluation and record the weight */
                                        Float weight = rfilter->eval(pos * invScale);
                                        m_weights[i * m_taps + j] = (Scalar) weight;
                                        sum += weight;
                                }

                                /* Normalize the contribution of each sample */
                                Float normalization = 1.0f / sum;
                                for (int j=0; j<m_taps; j++) {
                                        Scalar &value = m_weights[i * m_taps + j];
                                        value = (Scalar) ((Float) value * normalization);
                                }
                        }
                } else { /* Filtering mode */
                        m_weights = new Scalar[m_taps];
                        Float sum = 0;
                        for (int i=0; i<m_taps; i++) {
                                Scalar weight = (Scalar) rfilter->eval((Float) (i-m_halfTaps));
                                m_weights[i] = weight;
                                sum += (Float) weight;
                        }
                        Float normalization = 1.0f / sum;
                        for (int i=0; i<m_taps; i++) {
                                Scalar &value = m_weights[i];
                                value = (Scalar) ((Float) value * normalization);
                        }
                        m_fastStart = std::min(m_halfTaps, m_targetRes-1);
                        m_fastEnd = std::max(m_targetRes-m_halfTaps-1, 0);
                }

                /* Don't have overlapping fast start/end intervals when
                   the target image is very small compared to the source image */
                m_fastStart = std::min(m_fastStart, m_fastEnd);
        }

        /// Release all memory
        ~Resampler() {
                if (m_start)
                        delete[] m_start;
                if (m_weights)
                        delete[] m_weights;
        }

        /**
         * \brief Resample a multi-channel array and clamp the results
         * to a specified valid range
         *
         * This function is preferred if too large positive/negative values
         * due to ringing are unacceptable.
         *
         * \param source
         *     Source array of samples
         * \param target
         *     Target array of samples
         * \param sourceStride
         *     Stride of samples in the source array. A value
         *     of '1' implies that they are densely packed.
         * \param targetStride
         *     Stride of samples in the source array. A value
         *     of '1' implies that they are densely packed.
         * \param channels
         *     Number of channels to be resampled
         * \param min
         *     Minumum sample value after resampling
         * \param max
         *     Maximum sample value after resampling
         */
        void resampleAndClamp(const Scalar *source, size_t sourceStride,
                        Scalar *target, size_t targetStride, int channels,
                        Scalar min = (Scalar) 0, Scalar max = (Scalar) 1) {
                const int taps = m_taps, halfTaps = m_halfTaps;
                targetStride = channels * (targetStride - 1);
                sourceStride *= channels;

                if (m_start) {
                        /* Resample the left border region, while accounting for the boundary conditions */
                        for (int i=0; i<m_fastStart; ++i) {
                                int start = m_start[i];

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[i * taps + j];
                                    *target++ = std::min(max, std::max(min, result));
                                }

                                target += targetStride;
                        }

                        /* Use a faster, vectorizable loop for resampling the main portion */
                        for (int i=m_fastStart; i<m_fastEnd; ++i) {
                                int start = m_start[i];

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += source[sourceStride * (start + j) + ch] * m_weights[i * taps + j];
                                    *target++ = std::min(max, std::max(min, result));
                                }

                                target += targetStride;
                        }

                        /* Resample the right border region, while accounting for the boundary conditions */
                        for (int i=m_fastEnd; i<m_targetRes; ++i) {
                                int start = m_start[i];

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[i * taps + j];
                                    *target++ = std::min(max, std::max(min, result));
                                }

                                target += targetStride;
                        }
                } else {
                        /* Filter the left border region, while accounting for the boundary conditions */
                        for (int i=0; i<m_fastStart; ++i) {
                                int start = i - halfTaps;

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[j];
                                    *target++ = std::min(max, std::max(min, result));
                                }

                                target += targetStride;
                        }

                        /* Use a faster, vectorizable loop for filtering the main portion */
                        for (int i=m_fastStart; i<m_fastEnd; ++i) {
                                int start = i - halfTaps;

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += source[sourceStride * (start + j) + ch] * m_weights[j];
                                    *target++ = std::min(max, std::max(min, result));
                                }

                                target += targetStride;
                        }

                        /* Filter the right border region, while accounting for the boundary conditions */
                        for (int i=m_fastEnd; i<m_targetRes; ++i) {
                                int start = i - halfTaps;

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[j];
                                    *target++ = std::min(max, std::max(min, result));
                                }

                                target += targetStride;
                        }
                }
        }

        /**
         * \brief Resample a multi-channel array
         *
         * \param source
         *     Source array of samples
         * \param target
         *     Target array of samples
         * \param sourceStride
         *     Stride of samples in the source array. A value
         *     of '1' implies that they are densely packed.
         * \param targetStride
         *     Stride of samples in the source array. A value
         *     of '1' implies that they are densely packed.
         * \param channels
         *     Number of channels to be resampled
         */
        void resample(const Scalar *source, size_t sourceStride,
                        Scalar *target, size_t targetStride, int channels) {
                const int taps = m_taps, halfTaps = m_halfTaps;

                targetStride = channels * (targetStride - 1);
                sourceStride *= channels;

                if (m_start) {
                        /* Resample the left border region, while accounting for the boundary conditions */
                        for (int i=0; i<m_fastStart; ++i) {
                                int start = m_start[i];

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[i * taps + j];
                                        *target++ = result;
                                }

                                target += targetStride;
                        }

                        /* Use a faster, vectorizable loop for resampling the main portion */
                        for (int i=m_fastStart; i<m_fastEnd; ++i) {
                                int start = m_start[i];

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += source[sourceStride * (start + j) + ch] * m_weights[i * taps + j];
                                        *target++ = result;
                                }

                                target += targetStride;
                        }

                        /* Resample the right border region, while accounting for the boundary conditions */
                        for (int i=m_fastEnd; i<m_targetRes; ++i) {
                                int start = m_start[i];

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[i * taps + j];
                                        *target++ = result;
                                }

                                target += targetStride;
                        }
                } else {
                        /* Filter the left border region, while accounting for the boundary conditions */
                        for (int i=0; i<m_fastStart; ++i) {
                                int start = i - halfTaps;

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[j];
                                        *target++ = result;
                                }

                                target += targetStride;
                        }

                        /* Use a faster, vectorizable loop for filtering the main portion */
                        for (int i=m_fastStart; i<m_fastEnd; ++i) {
                                int start = i - halfTaps;

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += source[sourceStride * (start + j) + ch] * m_weights[j];
                                        *target++ = result;
                                }

                                target += targetStride;
                        }

                        /* Filter the right border region, while accounting for the boundary conditions */
                        for (int i=m_fastEnd; i<m_targetRes; ++i) {
                                int start = i - halfTaps;

                                for (int ch=0; ch<channels; ++ch) {
                                        Scalar result = 0;
                                        for (int j=0; j<taps; ++j)
                                                result += lookup(source, start + j, sourceStride, ch) * m_weights[j];
                                        *target++ = result;
                                }

                                target += targetStride;
                        }
                }
        }

private:
        FINLINE Scalar lookup(const Scalar *source, int pos, size_t stride, int offset) const {
                if (EXPECT_NOT_TAKEN(pos < 0 || pos >= m_sourceRes)) {
                        switch (m_bc) {
                                case ReconstructionFilter::EClamp:
                                        pos = math::clamp(pos, 0, m_sourceRes - 1);
                                        break;
                                case ReconstructionFilter::ERepeat:
                                        pos = math::modulo(pos, m_sourceRes);
                                        break;
                                case ReconstructionFilter::EMirror:
                                        pos = math::modulo(pos, 2*m_sourceRes);
                                        if (pos >= m_sourceRes)
                                                pos = 2*m_sourceRes - pos - 1;
                                        break;
                                case ReconstructionFilter::EZero:
                                        return (Scalar) 0;
                                case ReconstructionFilter::EOne:
                                        return (Scalar) 1;
                        }
                }
                return source[stride * pos + offset];
        }

private:
        ReconstructionFilter::EBoundaryCondition m_bc;
        int m_sourceRes;
        int m_targetRes;
        int *m_start;
        Scalar *m_weights;
        int m_fastStart, m_fastEnd;
        int m_taps, m_halfTaps;
};

extern MTS_EXPORT_CORE std::ostream &operator<<(std::ostream &os, const ReconstructionFilter::EBoundaryCondition &value);

MTS_NAMESPACE_END

#endif /* __MITSUBA_CORE_RFILTER_H_ */