LevelSetUtils-Impl.h

Go to the documentation of this file.

/*************************************************************************
> File Name: LevelSetUtils-Impl.h
> Project Name: CubbyFlow
> Author: Chan-Ho Chris Ohk
> Purpose: Level set util functions for CubbyFlow.
> Created Time: 2017/07/13
> Copyright (c) 2018, Chan-Ho Chris Ohk
*************************************************************************/

// Function fractionInside is originally from Christopher Batty's code:
//
// http://www.cs.ubc.ca/labs/imager/tr/2007/Batty_VariationalFluids/
//
// and
//
// https://github.com/christopherbatty/Fluid3D

#ifndef CUBBYFLOW_LEVEL_SET_UTILS_IMPL_H
#define CUBBYFLOW_LEVEL_SET_UTILS_IMPL_H

#include <Core/Utils/Constants.h>

#include <algorithm>
#include <cmath>

namespace CubbyFlow
{
    template <typename T>
    bool IsInsideSDF(T phi)
    {
        return phi < 0;
    }

    template <typename T>
    inline T SmearedHeavisideSDF(T phi)
    {
        if (phi > 1.5)
        {
            return 1;
        }
        
        if (phi < -1.5)
        {
            return 0;
        }

        return 0.5f + phi / 3.0 + 0.5f * (1 / PI<T>()) * std::sin(PI<T>() * phi / 1.5);
    }

    template <typename T>
    inline T SmearedDeltaSDF(T phi)
    {
        if (std::fabs(phi) > 1.5)
        {
            return 0;
        }

        return 1.0 / 3.0 + 1.0 / 3.0 * std::cos(PI<T>() * phi / 1.5);
    }

    template <typename T>
    T FractionInsideSDF(T phi0, T phi1)
    {
        if (IsInsideSDF(phi0) && IsInsideSDF(phi1))
        {
            return 1;
        }

        if (IsInsideSDF(phi0) && !IsInsideSDF(phi1))
        {
            return phi0 / (phi0 - phi1);
        }

        if (!IsInsideSDF(phi0) && IsInsideSDF(phi1))
        {
            return phi1 / (phi1 - phi0);
        }
        
        return 0;
    }

    template <typename T>
    void CycleArray(T* arr, int size)
    {
        T t = arr[0];
        
        for (int i = 0; i < size - 1; ++i)
        {
            arr[i] = arr[i + 1];
        }
        
        arr[size - 1] = t;
    }

    template <typename T>
    T FractionInside(T phiBottomLeft, T phiBottomRight, T phiTopLeft, T phiTopRight)
    {
        int insideCount = (phiBottomLeft < 0 ? 1 : 0) + (phiTopLeft < 0 ? 1 : 0) + (phiBottomRight < 0 ? 1 : 0) + (phiTopRight < 0 ? 1 : 0);
        T list[] = { phiBottomLeft, phiBottomRight, phiTopRight, phiTopLeft };

        if (insideCount == 4)
        {
            return 1;
        }
        
        if (insideCount == 3)
        {
            // Rotate until the positive value is in the first position
            while (list[0] < 0)
            {
                CycleArray(list, 4);
            }

            // Work out the area of the exterior triangle
            T side0 = 1 - FractionInsideSDF(list[0], list[3]);
            T side1 = 1 - FractionInsideSDF(list[0], list[1]);

            return 1 - 0.5f * side0 * side1;
        }
        
        if (insideCount == 2)
        {
            // Rotate until a negative value is in the first position, and the next
            // negative is in either slot 1 or 2.
            while (list[0] >= 0 || !(list[1] < 0 || list[2] < 0))
            {
                CycleArray(list, 4);
            }

            // The matching signs are adjacent
            if (list[1] < 0)
            {
                T side_left = FractionInsideSDF(list[0], list[3]);
                T side_right = FractionInsideSDF(list[1], list[2]);

                return 0.5f * (side_left + side_right);
            }

            // Matching signs are diagonally opposite
            // determine the center point's sign to disambiguate this case
            T middle_point = 0.25f * (list[0] + list[1] + list[2] + list[3]);
            if (middle_point < 0)
            {
                T area = 0;

                // First triangle (top left)
                T side1 = 1 - FractionInsideSDF(list[0], list[3]);
                T side3 = 1 - FractionInsideSDF(list[2], list[3]);

                area += 0.5f * side1 * side3;

                // Second triangle (top right)
                T side2 = 1 - FractionInsideSDF(list[2], list[1]);
                T side0 = 1 - FractionInsideSDF(list[0], list[1]);

                area += 0.5f * side0 * side2;

                return 1 - area;
            }
            
            T area = 0;

            // First triangle (bottom left)
            T side0 = FractionInsideSDF(list[0], list[1]);
            T side1 = FractionInsideSDF(list[0], list[3]);
            area += 0.5f * side0 * side1;

            // Second triangle (top right)
            T side2 = FractionInsideSDF(list[2], list[1]);
            T side3 = FractionInsideSDF(list[2], list[3]);
            area += 0.5f * side2 * side3;
            return area;
        }

        if (insideCount == 1)
        {
            // Rotate until the negative value is in the first position
            while (list[0] >= 0)
            {
                CycleArray(list, 4);
            }

            // Work out the area of the interior triangle, and subtract from 1.
            T side0 = FractionInsideSDF(list[0], list[3]);
            T side1 = FractionInsideSDF(list[0], list[1]);
            
            return 0.5f * side0 * side1;
        }
        
        return 0;
    }

    template <typename T>
    T DistanceToZeroLevelSet(T phi0, T phi1)
    {
        if (std::fabs(phi0) + std::fabs(phi1) > std::numeric_limits<double>::epsilon())
        {
            return std::fabs(phi0) / (std::fabs(phi0) + std::fabs(phi1));
        }
        
        return static_cast<T>(0.5);
    }
}

#endif