FracGen.cpp 9.7 KB
Newer Older
Jehferson Mello's avatar
Jehferson Mello committed
1 2 3 4 5 6 7 8 9 10 11
#include "FracGen.hpp"

#include <iostream>
#include <cfloat>
#include <functional>
#include <atomic>
#include <vector>

#include <tbb/parallel_for.h>
#include <tbb/blocked_range2d.h>

12 13
//static Vec3<tbFPType> boxMaxes;
//static Vec3<tbFPType> boxMins;
Jehferson Mello's avatar
Jehferson Mello committed
14 15


16
using  estimatorFunction = std::function<tbFPType(Vec3<tbFPType>, const Parameters<tbFPType>&)>;
Jehferson Mello's avatar
Jehferson Mello committed
17 18 19 20 21 22 23 24 25

/******************************************************************************
 *
 * Distance estimator functions and helpers
 *
 ******************************************************************************/



26
void sphereFold(Vec3<tbFPType>& z, tbFPType& dz, const Parameters<tbFPType>& params)
Jehferson Mello's avatar
Jehferson Mello committed
27
{
28 29
    tbFPType r2 = z.sqMod();
    if ( r2 < params.MinRadiusSq())
Jehferson Mello's avatar
Jehferson Mello committed
30 31
    {
        // linear inner scaling
32
        tbFPType temp = (params.FixedRadiusSq()/params.MinRadiusSq());
Jehferson Mello's avatar
Jehferson Mello committed
33 34 35
        z *= temp;
        dz *= temp;
    }
36
    else if(r2<params.FixedRadiusSq())
Jehferson Mello's avatar
Jehferson Mello committed
37 38
    {
        // this is the actual sphere inversion
39
        tbFPType temp =(params.FixedRadiusSq()/r2);
Jehferson Mello's avatar
Jehferson Mello committed
40 41 42 43 44
        z *= temp;
        dz*= temp;
    }
}

45
void boxFold(Vec3<tbFPType>& z, const Parameters<tbFPType>& params)
Jehferson Mello's avatar
Jehferson Mello committed
46
{
47
    z = z.clamp(-params.FoldingLimit(), params.FoldingLimit())* static_cast<tbFPType>(2.0f) - z;
Jehferson Mello's avatar
Jehferson Mello committed
48 49
}

50
tbFPType boxDist(const Vec3<tbFPType>& p, const Parameters<tbFPType>& params)
Jehferson Mello's avatar
Jehferson Mello committed
51 52 53 54 55 56 57
{
    /**
     * Distance estimator for a mandelbox
     *
     * Distance estimator adapted from
     * https://http://blog.hvidtfeldts.net/index.php/2011/11/distance-estimated-3d-fractals-vi-the-mandelbox/
     */
58 59 60 61
    const Vec3<tbFPType>& offset = p;
    tbFPType dr = params.BoxScale();
    Vec3<tbFPType> z{p};
    for (size_t n = 0; n < params.BoxIterations(); n++)
Jehferson Mello's avatar
Jehferson Mello committed
62
    {
63 64
        boxFold(z, params);       // Reflect
        sphereFold(z,dr, params);    // Sphere Inversion
Jehferson Mello's avatar
Jehferson Mello committed
65

66 67
        z = z * params.BoxScale() + offset;  // Scale & Translate
        dr = dr * std::abs(params.BoxScale()) + 1.0f;
Jehferson Mello's avatar
Jehferson Mello committed
68 69 70


    }
71
    tbFPType r = z.mod();
Jehferson Mello's avatar
Jehferson Mello committed
72 73 74 75
    return r/std::abs(dr);
}


76
tbFPType bulbDist(const Vec3<tbFPType>& p)
Jehferson Mello's avatar
Jehferson Mello committed
77 78 79 80 81 82 83 84 85 86
{

    /**
     * Distance estimator for a mandelbulb
     *
     * Distance estimator adapted from
     * https://www.iquilezles.org/www/articles/mandelbulb/mandelbulb.htm
     * https://www.shadertoy.com/view/ltfSWn
     */

87 88
    Vec3<tbFPType> w = p;
    tbFPType m = w.sqMod();
Jehferson Mello's avatar
Jehferson Mello committed
89 90

    //vec4 trap = vec4(abs(w),m);
91
    tbFPType dz = 3.0f;
Jehferson Mello's avatar
Jehferson Mello committed
92 93 94 95 96


    for( int i=0; i<4; i++ )
    {
#if 1
97 98
        tbFPType m2 = m*m;
        tbFPType m4 = m2*m2;
Jehferson Mello's avatar
Jehferson Mello committed
99 100
        dz = 8.0f*sqrt(m4*m2*m)*dz + 1.0f;

101 102 103
        tbFPType x = w.X(); tbFPType x2 = x*x; tbFPType x4 = x2*x2;
        tbFPType y = w.Y(); tbFPType y2 = y*y; tbFPType y4 = y2*y2;
        tbFPType z = w.Z(); tbFPType z2 = z*z; tbFPType z4 = z2*z2;
Jehferson Mello's avatar
Jehferson Mello committed
104

105 106 107 108
        tbFPType k3 = x2 + z2;
        tbFPType k2 = 1/sqrt( k3*k3*k3*k3*k3*k3*k3 );
        tbFPType k1 = x4 + y4 + z4 - 6.0f*y2*z2 - 6.0f*x2*y2 + 2.0f*z2*x2;
        tbFPType k4 = x2 - y2 + z2;
Jehferson Mello's avatar
Jehferson Mello committed
109 110 111 112 113 114 115 116

        w.setX(p.X() +  64.0f*x*y*z*(x2-z2)*k4*(x4-6.0f*x2*z2+z4)*k1*k2);
        w.setY(p.Y() + -16.0f*y2*k3*k4*k4 + k1*k1);
        w.setZ(p.Z() +  -8.0f*y*k4*(x4*x4 - 28.0f*x4*x2*z2 + 70.0f*x4*z4 - 28.0f*x2*z2*z4 + z4*z4)*k1*k2);
#else
        dz = 8.0*pow(sqrt(m),7.0)*dz + 1.0;
        //dz = 8.0*pow(m,3.5)*dz + 1.0;

117 118 119 120
        tbFPType r = w.mod();
        tbFPType b = 8.0*acos( w.Y()/r);
        tbFPType a = 8.0*atan2( w.X(), w.Z() );
        w = p + Vec3<tbFPType>( sin(b)*sin(a), cos(b), sin(b)*cos(a) ) * pow(r,8.0);
Jehferson Mello's avatar
Jehferson Mello committed
121 122 123 124 125 126 127 128 129 130 131 132
#endif

       // trap = min( trap, vec4(abs(w),m) );

        m = w.sqMod();
        if( m > 256.0f )
            break;
    }

    return 0.25f*log(m)*sqrt(m)/dz;
}

133
tbFPType sphereDist(Vec3<tbFPType> p)
Jehferson Mello's avatar
Jehferson Mello committed
134
{
135
    tbFPType radius = 2.f;
Jehferson Mello's avatar
Jehferson Mello committed
136 137 138 139 140 141 142 143 144
    return p.mod() - radius;
}

/******************************************************************************
 *
 * Coulouring functions and helpers
 *
 ******************************************************************************/

145
RGBA HSVtoRGB(int H, tbFPType S, tbFPType V)
Jehferson Mello's avatar
Jehferson Mello committed
146 147 148 149 150 151 152 153
{

    /**
     * adapted from
     * https://gist.github.com/kuathadianto/200148f53616cbd226d993b400214a7f
     */

    RGBA output;
154 155 156 157
    tbFPType C = S * V;
    tbFPType X = C * (1 - std::abs(std::fmod(H / 60.0, 2) - 1));
    tbFPType m = V - C;
    tbFPType Rs, Gs, Bs;
Jehferson Mello's avatar
Jehferson Mello committed
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202

    if(H >= 0 && H < 60)
    {
        Rs = C;
        Gs = X;
        Bs = 0;
    }
    else if(H >= 60 && H < 120)
    {
        Rs = X;
        Gs = C;
        Bs = 0;
    }
    else if(H >= 120 && H < 180)
    {
        Rs = 0;
        Gs = C;
        Bs = X;
    }
    else if(H >= 180 && H < 240)
    {
        Rs = 0;
        Gs = X;
        Bs = C;
    }
    else if(H >= 240 && H < 300)
    {
        Rs = X;
        Gs = 0;
        Bs = C;
    }
    else {
        Rs = C;
        Gs = 0;
        Bs = X;
    }

    output.setR(Rs + m);
    output.setG(Gs + m);
    output.setB(Bs + m);
    output.setA( 1.0f );

    return output;
}

203
RGBA getColour(const Vec4<tbFPType>& steps, const Parameters<tbFPType>& params )
Jehferson Mello's avatar
Jehferson Mello committed
204 205 206
{
    RGBA colour;

207
    Vec3<tbFPType> position(steps.X(),steps.Y(),steps.Z());
Jehferson Mello's avatar
Jehferson Mello committed
208 209


210
    if(steps.W() >= params.MaxRaySteps())
Jehferson Mello's avatar
Jehferson Mello committed
211
    {
212
        return RGBA(params.BgRed(),params.BgGreen(),params.BgBlue(),params.BgAlpha());
Jehferson Mello's avatar
Jehferson Mello committed
213 214
    }

215
//        This is a good place to tbFPType check your bounds if you need to
Jehferson Mello's avatar
Jehferson Mello committed
216 217 218 219 220 221 222 223 224

//        boxMins.setX(std::min(boxMins.X(),position.X()));
//        boxMins.setY(std::min(boxMins.Y(),position.Y()));
//        boxMins.setZ(std::min(boxMins.Z(),position.Z()));

//        boxMaxes.setX(std::max(boxMaxes.X(),position.X()));
//        boxMaxes.setY(std::max(boxMaxes.Y(),position.Y()));
//        boxMaxes.setZ(std::max(boxMaxes.Z(),position.Z()));

225 226
    tbFPType saturation = params.SatValue();
    tbFPType hueVal = (position.Z() * params.HueFactor()) + params.HueOffset();
227 228
    int hue = static_cast<int>( trunc(fmod(hueVal, 360.0f) ) );
    hue = hue < 0 ? 360 + hue: hue;
229
    tbFPType value = params.ValueRange()*(1.0f - std::min(steps.W()*params.ValueFactor()/params.MaxRaySteps(), params.ValueClamp()));
Jehferson Mello's avatar
Jehferson Mello committed
230 231 232 233

    colour = HSVtoRGB(hue, saturation, value);

//    Simplest colouring, based only on steps (roughly distance from camera)
234
    //colour = RGBA(value,value,value,1.0f);
Jehferson Mello's avatar
Jehferson Mello committed
235 236 237 238 239 240 241 242 243 244

    return colour;
}

/******************************************************************************
 *
 * Ray marching functions and helpers
 *
 ******************************************************************************/

245
Vec4<tbFPType> trace(const Camera<tbFPType>& cam, const Parameters<tbFPType>& params, size_t x, size_t y, const estimatorFunction& f)
Jehferson Mello's avatar
Jehferson Mello committed
246 247 248 249 250 251
{
    /**
     * This function taken from
     * http://blog.hvidtfeldts.net/index.php/2011/06/distance-estimated-3d-fractals-part-i/
     */

252
    tbFPType totalDistance = 0.0f;
Jehferson Mello's avatar
Jehferson Mello committed
253 254
    unsigned int steps;

255
    Vec3<tbFPType> pixelPosition = cam.ScreenTopLeft() + (cam.ScreenRight() * static_cast<tbFPType>(x)) + (cam.ScreenUp() * static_cast<tbFPType>(y));
Jehferson Mello's avatar
Jehferson Mello committed
256

257
    Vec3<tbFPType> rayDir = pixelPosition - cam.Pos();
Jehferson Mello's avatar
Jehferson Mello committed
258 259
    rayDir.normalise();

260 261
    Vec3<tbFPType> p;
    for (steps=0; steps < params.MaxRaySteps(); steps++)
Jehferson Mello's avatar
Jehferson Mello committed
262
    {
263
        p = cam.Pos() + (rayDir * totalDistance);
264
        tbFPType distance = f(p, params);
Jehferson Mello's avatar
Jehferson Mello committed
265
        totalDistance += distance;
266
        if (distance < params.CollisionMinDist()) break;
Jehferson Mello's avatar
Jehferson Mello committed
267 268
    }
    //return both the steps and the actual position in space for colouring purposes
269
    return Vec4<tbFPType>{p,static_cast<tbFPType>(steps)};
Jehferson Mello's avatar
Jehferson Mello committed
270 271 272
}

void traceRegion(colourVec& data,
273 274
                 const Camera<tbFPType>& cam,
                 const Parameters<tbFPType>& params,
Jehferson Mello's avatar
Jehferson Mello committed
275 276 277 278 279 280 281
                 const estimatorFunction& f,
                 const tbb::blocked_range2d<size_t,size_t>& range)
{
    for(size_t h = range.rows().begin(); h < range.rows().end(); h++)
    {
        for(size_t w = range.cols().begin(); w < range.cols().end(); w++ )
        {
282
            data[(h*cam.ScreenWidth())+w] = getColour(trace(cam, params, w, h, f), params);
Jehferson Mello's avatar
Jehferson Mello committed
283 284 285 286 287 288 289 290 291 292
        }
    }
}

/******************************************************************************
 *
 * Thread spawning section
 *
 ******************************************************************************/

293
bool FracGen::Generate()
Jehferson Mello's avatar
Jehferson Mello committed
294
{
295 296 297 298
    if(outBuffer->size() != cam->ScreenWidth()*cam->ScreenHeight())
    {
        outBuffer->resize(cam->ScreenWidth()*cam->ScreenHeight());
    }
Jehferson Mello's avatar
Jehferson Mello committed
299 300

    bool finishedGeneration = false;
301
    int heightStep = bench ? cam->ScreenHeight() : 10;
Jehferson Mello's avatar
Jehferson Mello committed
302

303
    //estimatorFunction bulb(bulbDist);
Jehferson Mello's avatar
Jehferson Mello committed
304
    estimatorFunction box(boxDist);
305
    //estimatorFunction sphere(sphereDist);
Jehferson Mello's avatar
Jehferson Mello committed
306

307
    auto tbbTrace = [this, f = box](const tbb::blocked_range2d<size_t,size_t>& range)
308
            {traceRegion(*(this->outBuffer), *(this->cam), *(this->parameters), f, range); };
Jehferson Mello's avatar
Jehferson Mello committed
309 310 311 312


    // Minor optimization in TBB: Since we can define the range here, we don't need to check bounds in
    // traceRegion. This is very similar to how SYCL and OpenCL do things
313 314
    tbb::blocked_range2d<size_t,size_t> range( lastHeight, std::min(lastHeight+heightStep, static_cast<size_t>(cam->ScreenHeight()))
                                             , 0, cam->ScreenWidth());
Jehferson Mello's avatar
Jehferson Mello committed
315 316
    tbb::parallel_for(range, tbbTrace);

317
    lastHeight+= heightStep;
Jehferson Mello's avatar
Jehferson Mello committed
318 319 320 321 322 323
    /*
     * TBB doesn't seem to like the way of doing things where I know how
     * many tasks there are and am collecting the return from them.
     * At least not in the initial tutorials and beginner documentation
     */

324
    if (lastHeight >=static_cast<size_t>(cam->ScreenHeight()) )
Jehferson Mello's avatar
Jehferson Mello committed
325
    {
326
        lastHeight = 0;
Jehferson Mello's avatar
Jehferson Mello committed
327 328 329 330 331 332
        finishedGeneration = true;
    }

    return finishedGeneration;
}

333
FracGen::FracGen(bool benching, CameraPtr c, ParamPtr p)
Jehferson Mello's avatar
Jehferson Mello committed
334
    : bench{benching}
335 336
    , cam{c}
    , parameters{p}
337
    , lastHeight{0}
Jehferson Mello's avatar
Jehferson Mello committed
338
{
339
    outBuffer = std::make_shared< colourVec >(cam->ScreenWidth()*cam->ScreenHeight());
Jehferson Mello's avatar
Jehferson Mello committed
340 341 342 343 344 345 346 347 348 349 350 351 352 353

    static bool once = false;
    if(!bench || !once)
    {
      std::cout << "Running on TBB "<< TBB_VERSION_MAJOR <<"."<< TBB_VERSION_MINOR << std::endl;
      once = true;
    }
}

FracGen::~FracGen()
{}