renderer_detail_normal_mrf.cpp

#include "renderer_detail_normal_mrf.h"
#include "renderer_detail_normal_mrf_helpers.h"

#include <limits>
#include <cmath>




namespace Renderer {


namespace detail {




/*==============================================================================
    NormalMRF methods
==============================================================================*/


void NormalMRF::CalculateNormals(
    Vector< float, 3 >* const normals,
    Vector< double, 2 > const fieldOfView,
    Vector< double, 3 > const xx,
    Vector< double, 3 > const yy,
    Vector< double, 3 > const zz
)
{
    assert( m_coarseness == 0 );

    Vector< double, 2 > const spread = 2 * fieldOfView / m_dimension;
    Vector< double, 3 > const spreadX = spread[ 0 ] * xx;
    Vector< double, 3 > const spreadY = spread[ 1 ] * yy;

    #pragma omp parallel for
    for ( int ii = 0; ii < static_cast< int >( m_dimension[ 1 ] ); ++ii ) {

        Vector< double, 3 > ray = (
            (      0.5 - 0.5 * m_dimension[ 0 ] ) * spreadX -
            ( ii + 0.5 - 0.5 * m_dimension[ 1 ] ) * spreadY -
            zz
        );

        unsigned int const offset = ii * m_dimension[ 0 ];
        for ( int jj = 0; jj < static_cast< int >( m_dimension[ 0 ] ); ++jj ) {

            Vector< double, 3 > normal = { { 0, 0, 0 } };

            double const reciprocalDepth = m_states[ offset + jj ].first;
            if ( reciprocalDepth > 0 ) {

                double const depth = 1 / reciprocalDepth;
                Vector< double, 2 > const gradient = -Square( depth ) * m_states[ offset + jj ].second;

                Vector< double, 3 > horizontalDifference = gradient[ 0 ] * ray + ( gradient[ 0 ] + depth ) * spreadX;
                Vector< double, 3 > verticalDifference   = gradient[ 1 ] * ray - ( gradient[ 1 ] + depth ) * spreadY;

                normal = Cross( horizontalDifference, -verticalDifference );
                normal /= normal.Norm();
            }
            normals[ offset + jj ] = normal;

            ray += spreadX;
        }
    }
}


void NormalMRF::SmoothHorizontal( bool const odd, bool const decreasing ) {

    unsigned int const size = ( 1u << m_coarseness );

    int jjBegin = 0;
    int jjEnd   = ( ( m_dimension[ 0 ] - 1 ) & ~( size - 1 ) );
    int jjDelta = size;
    if ( decreasing ) {

        std::swap( jjBegin, jjEnd );
        jjDelta = -jjDelta;
    }
    jjEnd += jjDelta;

    #pragma omp parallel for
    for ( int ii = ( odd ? size : 0 ); ii < static_cast< int >( m_dimension[ 1 ] ); ii += size * 2 ) {

        for ( int jj = jjBegin; jj != jjEnd; jj += jjDelta ) {

            assert( ( jj >= 0 ) && ( jj < m_dimension[ 0 ] ) );
            Vector< unsigned int, 2 > const coordinates = { { jj, ii } };
            SolveCardinal(
                m_depths.get(),
                m_states.get(),
                m_dimension,
                coordinates,
                size,
                size,
                m_spread,
                m_voxelSize,
                Square( size ) * m_observationPrecision,
                m_gradientPrecision,
                m_minimumEdgeDelta,
                m_maximumEdgeDelta
            );
        }
    }
}


void NormalMRF::SmoothVertical( bool const odd, bool const decreasing ) {

    unsigned int const size = ( 1u << m_coarseness );

    int jjBegin = 0;
    int jjEnd   = ( ( m_dimension[ 1 ] - 1 ) & ~( size - 1 ) );
    int jjDelta = size;
    if ( decreasing ) {

        std::swap( jjBegin, jjEnd );
        jjDelta = -jjDelta;
    }
    jjEnd += jjDelta;

    #pragma omp parallel for
    for ( int ii = ( odd ? size : 0 ); ii < static_cast< int >( m_dimension[ 0 ] ); ii += size * 2 ) {

        for ( int jj = jjBegin; jj != jjEnd; jj += jjDelta ) {

            assert( ( jj >= 0 ) && ( jj < m_dimension[ 1 ] ) );
            Vector< unsigned int, 2 > const coordinates = { { ii, jj } };
            SolveCardinal(
                m_depths.get(),
                m_states.get(),
                m_dimension,
                coordinates,
                size,
                size,
                m_spread,
                m_voxelSize,
                Square( size ) * m_observationPrecision,
                m_gradientPrecision,
                m_minimumEdgeDelta,
                m_maximumEdgeDelta
            );
        }
    }
}


void NormalMRF::SmoothCardinalCheckerboard( bool const odd ) {

    unsigned int const size = ( 1u << m_coarseness );

    #pragma omp parallel for
    for ( int ii = 0; ii < static_cast< int >( m_dimension[ 1 ] ); ii += size ) {

        for ( int jj = ( ( ( ii & size ) != 0 ) == odd ) ? 0 : size; jj < static_cast< int >( m_dimension[ 0 ] ); jj += size * 2 ) {

            Vector< unsigned int, 2 > const coordinates = { { jj, ii } };
            SolveCardinal(
                m_depths.get(),
                m_states.get(),
                m_dimension,
                coordinates,
                size,
                size,
                m_spread,
                m_voxelSize,
                Square( size ) * m_observationPrecision,
                m_gradientPrecision,
                m_minimumEdgeDelta,
                m_maximumEdgeDelta
            );
        }
    }
}


void NormalMRF::SmoothDiagonalCheckerboard( bool const odd ) {

    unsigned int const size = ( 1u << m_coarseness );
    double const scale = size * std::sqrt( 2 );

    #pragma omp parallel for
    for ( int ii = odd ? size : 0; ii < static_cast< int >( m_dimension[ 1 ] ); ii += size * 2 ) {

        for ( int jj = odd ? size : 0; jj < static_cast< int >( m_dimension[ 0 ] ); jj += size * 2 ) {

            Vector< unsigned int, 2 > const coordinates = { { jj, ii } };
            SolveDiagonal(
                m_depths.get(),
                m_states.get(),
                m_dimension,
                coordinates,
                size,
                scale,
                m_spread,
                m_voxelSize,
                Square( scale ) * m_observationPrecision,
                m_gradientPrecision,
                m_minimumEdgeDelta,
                m_maximumEdgeDelta
            );
        }
    }
}


void NormalMRF::InitializeStates() {

    unsigned int const size = ( 1u << m_coarseness );

    Vector< double, 2 > const gradient0 = { { 0, 0 } };

    #pragma omp parallel for
    for ( int ii = 0; ii < static_cast< int >( m_dimension[ 1 ] ); ii += size ) {

        unsigned int const offset = ii * m_dimension[ 0 ];
        for ( int jj = 0; jj < static_cast< int >( m_dimension[ 0 ] ); jj += size ) {

            double const reciprocalDepth = 1 / m_depths[ offset + jj ];

            m_states[ offset + jj ].first = reciprocalDepth;
            m_states[ offset + jj ].second = gradient0;
        }
    }
}




}    // namespace detail


}    // namespace Renderer