PCV.cpp

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// -                        CloudViewer: www.cloudViewer.org                  -
// ----------------------------------------------------------------------------
// Copyright (c) 2018-2024 www.cloudViewer.org
// SPDX-License-Identifier: MIT
// ----------------------------------------------------------------------------
 
#include "PCV.h"
 
#include "PCVContext.h"
 
// Qt
#include <QString>
 
#ifdef USE_VLD
// VLD
#include <vld.h>
#endif
 
// System
#include <assert.h>
#include <string.h>
 
#include <algorithm>
#include <vector>
 
using namespace CVLib;
 
static int gcd(int num1, int num2) {
    int remainder = (num2 % num1);
 
    return (remainder != 0 ? gcd(remainder, num1) : num1);
}
 
 
static bool SampleSphere(unsigned N, std::vector<CCVector3>& dirs) {
    static const double c_eps = 2.2204e-16;
    static const double c_twist = 4.0;
 
    if (N == 0) {
        assert(false);
        return false;
    }
 
    try {
        dirs.resize(N, CCVector3(0, 0, 1));
    } catch (const std::bad_alloc&) {
        // not enough memory
        return false;
    }
 
    if (N == 1) {
        return true;
    }
 
    try {
        double area = (4 * M_PI) / N;
        double beta = acos(1.0 - 2.0 / N);  // return in [0,pi/2] as '1-2/N'
                                            // goes from 0 to 1 when N --> inf
        double gamma = M_PI - 2 * beta;     // in [0,pi]
        double fuzz = c_eps * 2 * N;
 
        int Ltemp = static_cast<int>(ceil(gamma / sqrt(area) - fuzz));
        int L = 2 + std::max(Ltemp, 1);
 
        // init mbar
        std::vector<double> mbar;
        mbar.resize(L, 0);
        assert(L >= 3);
        {
            mbar[0] = 1.0;
            double theta = gamma / (L - 2);
            for (int i = 1; i < L - 1; ++i) {
                mbar[i] =
                        N *
                        (cos(theta * (i - 1) + beta) - cos(theta * i + beta)) /
                        2;
            }
            mbar[L - 1] = 1.0;
        }
 
        // init m
        std::vector<int> m;
        m.resize(L, 0);
        {
            m[0] = 1;
 
            double alpha = 0.0;
            for (int i = 1; i < L; ++i) {
                double f = floor(mbar[i] + alpha + fuzz);
                if ((mbar[i] - f) >= 0.5) f = ceil(mbar[i] + alpha - fuzz);
                m[i] = static_cast<int>(f);
 
                alpha += mbar[i] - m[i];
            }
        }
 
        // now we can build the rays
        {
            std::vector<double> offset;
            offset.resize(L - 1, 0);
 
            double z = 1.0 - static_cast<double>(2 + m[1]) / N;
 
            unsigned int rayIndex = 1;  // the first one is already (0,0,1)
            for (int i = 1; i < L - 1; ++i) {
                if (m[i - 1] != 0 && m[i] != 0) {
                    offset[i] = offset[i - 1] +
                                static_cast<double>(gcd(m[i], m[i - 1])) /
                                        (2 * m[i] * m[i - 1]) +
                                std::min<double>(c_twist,
                                                 floor(m[i - 1] / c_twist)) /
                                        m[i - 1];
                } else {
                    offset[i] = 0.0;
                }
 
                double temp = static_cast<double>(m[i]) / N;
 
                double h = cos((acos(z + temp) + acos(z - temp)) / 2);
                double r = sqrt(1.0 - h * h);
 
                for (int j = 0; j < m[i]; ++j) {
                    double theta = 2.0 * M_PI *
                                   (offset[i] + static_cast<double>(j) / m[i]);
 
                    dirs[rayIndex++] = CCVector3(
                            static_cast<PointCoordinateType>(r * cos(theta)),
                            static_cast<PointCoordinateType>(r * sin(theta)),
                            static_cast<PointCoordinateType>(h));
                }
 
                z -= static_cast<double>(m[i] + m[i + 1]) / N;
            }
 
            assert(rayIndex + 1 == N);
        }
    } catch (const std::bad_alloc&) {
        // not enough memory
        return false;
    }
 
    dirs[N - 1] = CCVector3(0, 0, -1);
 
    return true;
}
 
bool PCV::GenerateRays(unsigned numberOfRays,
                       std::vector<CCVector3>& rays,
                       bool mode360 /*=true*/) {
    // generates light directions
    unsigned rayCount = numberOfRays * (mode360 ? 1 : 2);
    if (!SampleSphere(rayCount, rays)) {
        return false;
    }
 
    // we keep only the light directions that meets input parameters (non
    // predictible if not in 360?mode!)
    if (!mode360) {
        unsigned lastIndex = 0;
        for (size_t i = 0; i < rays.size(); ++i) {
            if (rays[i].z < 0) {
                if (lastIndex != i) {
                    rays[lastIndex] = rays[i];
                }
                ++lastIndex;
            }
        }
        rayCount = lastIndex;
        rays.resize(rayCount);
    }
 
    return true;
}
 
int PCV::Launch(unsigned numberOfRays,
                GenericCloud* vertices,
                GenericMesh* mesh /*=0*/,
                bool meshIsClosed /*=false*/,
                bool mode360 /*=true*/,
                unsigned width /*=1024*/,
                unsigned height /*=1024*/,
                CVLib::GenericProgressCallback* progressCb /*=0*/,
                QString entityName /*=QString()*/) {
    // generates light directions
    std::vector<CCVector3> rays;
    if (!GenerateRays(numberOfRays, rays, mode360)) {
        return -2;
    }
 
    if (!Launch(rays, vertices, mesh, meshIsClosed, width, height, progressCb,
                entityName)) {
        return -1;
    }
 
    return static_cast<int>(rays.size());
}
 
bool PCV::Launch(std::vector<CCVector3>& rays,
                 CVLib::GenericCloud* vertices,
                 CVLib::GenericMesh* mesh /*=0*/,
                 bool meshIsClosed /*=false*/,
                 unsigned width /*=1024*/,
                 unsigned height /*=1024*/,
                 CVLib::GenericProgressCallback* progressCb /*=0*/,
                 QString entityName /*=QString()*/) {
    if (rays.empty()) return false;
 
    if (!vertices || !vertices->enableScalarField()) return false;
 
    // vertices/points
    unsigned numberOfPoints = vertices->size();
    // rays
    unsigned numberOfRays = static_cast<unsigned>(rays.size());
 
    // for each vertex we keep count of the number of light directions for which
    // it is "illuminated"
    std::vector<int> visibilityCount;
    try {
        visibilityCount.resize(numberOfPoints, 0);
    } catch (const std::bad_alloc&) {
        // not enough memory?
        return false;
    }
 
    /*** Main illumination loop ***/
 
    CVLib::NormalizedProgress nProgress(progressCb, numberOfRays);
    if (progressCb) {
        if (progressCb->textCanBeEdited()) {
            progressCb->setMethodTitle("ShadeVis");
            QString infoStr;
            if (!entityName.isEmpty()) infoStr = entityName + "\n";
            infoStr.append(QString("Rays: %1").arg(numberOfRays));
            if (mesh)
                infoStr.append(QString("\nFaces: %1").arg(mesh->size()));
            else
                infoStr.append(QString("\nVertices: %1").arg(numberOfPoints));
            progressCb->setInfo(qPrintable(infoStr));
        }
        progressCb->update(0);
        progressCb->start();
    }
 
    bool success = true;
 
    // must be done after progress dialog display!
    PCVContext win;
    if (win.init(width, height, vertices, mesh, meshIsClosed)) {
        for (unsigned i = 0; i < numberOfRays; ++i) {
            // set current 'light' direction
            win.setViewDirection(rays[i]);
 
            // flag viewed vertices
            win.GLAccumPixel(visibilityCount);
 
            if (progressCb && !nProgress.oneStep()) {
                success = false;
                break;
            }
        }
 
        if (success) {
            // we convert per-vertex accumulators to an 'intensity' scalar field
            for (unsigned j = 0; j < numberOfPoints; ++j) {
                ScalarType visValue =
                        static_cast<ScalarType>(visibilityCount[j]) /
                        numberOfRays;
                vertices->setPointScalarValue(j, visValue);
            }
        }
    } else {
        success = false;
    }
 
    return success;
}