Трассиврока лучей и мягкие тени, алгоритм генерации лучей внутри конуса (ray tracing + soft shadows)

Решил разобраться в трассировке лучей. Для начала думаю все это дело реализовать чисто на процессоре, то есть при помощи обычного C++ кода, а затем уже буду пробовать это все перенести на видео-карту (OpenCL или Compute Shaders), превратив рекурсивный алгоритм в итеративный (это еще предстоит). Все по началу шло хорошо, сцена с плоскостью (двумя треугольниками) и несколькими сферами рисовалась довольно быстро (даже с преломлениями и отражением), и тут я решил добавить мягкие тени...

Поскольку у меня источники света пока-что исключительно сферические, то алгоритм по сути заключается в том чтобы от каждой точки пересечения пустить несколько лучей, которые бы "покрывали" весь диск сферы источника света, а затем смотреть - сколько из запущенных лучей достигли источника и сколько было перекрыто преградой. Отношение достигших источника лучей ко всем - и будет интенсивностью освещенности источником. Как-то так:

Но как этот самый набор лучей получить?

Я рассуждал так - все лучи должны откланяться от "основного луча" не более чем на угол, тангенс которого это отношение радиуса сферы и длины вектора до ее центра. Но даже если я этот угол выясню, что мне это даст? Как мне это поможет в построении векторов (пусть даже случайных) в пределах этого конуса?

В итоге я решил действовать чуть иначе и оставить этот угол в покое. Я решил, что во первых, нужно построить 4 вектора перпендикулярных основному вектору (направленному в центр источника). Достичь этого я решил следующим образом:

1. Получаем вектор №1 векторно перемножая основной (центральный) вектор и чуть отклоненный в сторону по всем осям. Нормализуем получившийся вектор.
2. Получаем вектор №2 векторно перемножая вектор №1 и основной, нормализуем результат.
3. Получаем вектор №3 путем инверсии вектора №1 (добавив минус)
4. Получаем вектор №4 путем инверсии вектора #2 (добавив минус)

Таким образом у нас получается "крест" из 4 векторов, все из которых равномерно расходятся во все сторону от центрального вектора:

Далее я планировал умножить эти векторы на отношение радиуса источника к длине центрального (не нормированного) вектора, получив тем самым такие вектора, которые нужно прибавить к центральному (нормированному), чтобы отклонить его таким образом, чтобы он стал указывать на границы диска сферы. Ну а если этот самый коэффициент еще поделить - можно еще добавить промежуточные лучи! Но как бы сильно я не делил коэффициент все они будут распространяться исключительно по этому "кресту". Что же делать?

Я решил что нужно добавить промежуточных векторов, перед всем этим. И достичь этого можно путем суммы соседних векторов (и последующей нормализацией результатов). В итоге окончательный алгоритм получился таким:

1. Строим 4 вектора ортогональных центральному
2. При помощи суммирования векторов сегментируем все покуда не будет получено достаточно этих ортогональных векторов идущих во все стороны от центрального.
3. Получаем коэффициент (отношение радиуса к длине центрального не нормированного)
4. В цикле "дробим коэффициент" и получаем столько итоговых векторов, сколько нужно. Создаем луч для каждого вектора

Итоговая функция которая все это делает у меня получилась такой:

std::vector<Ray> RayTracer::generateSphericalRaySet(
        const math::Vec3<float> &origin,
        const math::Vec3<float> &sphereCenter,
        const float &sphereRadius,
        size_t segmentationLevel,
        size_t perSegment)
{
    // Центральный вектор направленный в сторону сферы
    auto toSphereCenter = sphereCenter - origin;
    // Нормализованный центральный вектор
    auto central = math::Normalize(toSphereCenter);
    // Расстрояние до цента источника от начала
    auto distance = math::Length(toSphereCenter);
    // Результат
    std::vector<Ray> result = {Ray(origin,central)};
    // Получаем первые 4 основные вектора ортогональных основному вектору (направленному к центру сферы)
    auto v1 = math::Normalize(math::Cross(central,central + math::Vec3<float>(1e-3,1e-3,1e-3)));
    auto v2 = math::Normalize(math::Cross(central,v1));
    auto v3 = -v1;
    auto v4 = -v2;
    // Добавляем вектора в связанный список
    std::forward_list<math::Vec3<float>> vectors;
    vectors.push_front(v1);vectors.push_front(v2);
    vectors.push_front(v3);vectors.push_front(v4);
    // Если нужно сегментировать (получить вектора находящиеся между)
    if(segmentationLevel > 0)
    {
        // Каждая итерация удваивает кол-во векторов, таким образом окружность сегменитируется
        // и появляется возможность для большей точности
        for(size_t l = 0; l < segmentationLevel; l++)
        {
            for(auto current = vectors.begin(); current != vectors.end(); ++current)
            {
                auto next = std::next(current) != vectors.end() ? std::next(current) : vectors.begin();
                vectors.insert_after(current,math::Normalize(*current + *next));
                ++current;
            }
        }
    }
    // Коэффициент пропорциональности, отношение радиуса диска сферы и длинны до центра
    // Он дает понять какое нужно отклонение (какой должна быть длина перпедикулярного вектора) добавить к
    // базовому чтобы новый вектор указывал в сторону границ диска сферы
    auto fullBias = (sphereRadius-0.01f) / distance;
    // Создаем лучи
    for(const auto& segmentVector : vectors)
    {
        for(size_t i = 0; i < perSegment; i++)
        {
            auto cof = (static_cast<float>(i + 1) / static_cast<float>(perSegment)) * fullBias;
            result.emplace_back(Ray(origin,(segmentVector * cof) + central));
        }
    }
    return result;
}

И когда я применил этот алгоритм, подогнав кол-во лучей, я получил вот такую картинку:

НО!

Рисование этой картинки с мягкими тенями заняло около 20 секунд! В то время как картинка БЕЗ мягких теней рисовалась меньше секунды. При этом сцена здесь крайне проста. А что будет если объектов будет больше? Я решил проверить что будет если буду просто генерировать набор, но не использовать его - оказалось что даже если не делать проверку лучей, сама их генерация занимает тоже порядочно. При этом, я понимаю, что в дальнейшем, когда я буду превращать это все в итеративный алгоритм, для работы на видео-карте, все равно придется ради мягких теней делать проверку НАБОРА лучей для каждого луча, и в отдельных потоках это никак сделать не выйдет...

Как же тогда правильно? Есть какой-то более оптимальный алгоритм? Наверняка ведь есть уже что-то общепринятое? Сколько я не искал, так и не нашел ничего внятного. Заранее спасибо.