// To reproduce the object in the demo by Henrik Rydgård/@greggman/alteredq//
/* AlteredQualia's Marching Cubes, C++ version, based on AQ's code based on Henrik Rydgård and @greggman.
 * https://github.com/WebGLSamples/WebGLSamples.github.io/blob/master/blob/marching_cubes.js
 *
 * Based on alteredq's version  https://github.com/mrdoob/three.js/blob/master/examples/js/MarchingCubes.js
 *
 * Port of greggman's ThreeD version of marching cubes to Three.js
 * http://webglsamples.googlecode.com/hg/blob/blob.html
*/
class meta_ball_example {

    meta_balls(int num_blobs, REAL time, REAL scale) {
        int numblobs = num_blobs;  // default: 4
        for (int ball_i = 0; ball_i < numblobs; ball_i++) {
            REAL D = 1;
            REAL ballx = sin(ball_i + 1.26 * time * (1.03 + 0.5*cos(0.21 * ball_i))) * 0.27 * D + 0.5;
            REAL bally = std::abs(cos(ball_i + 1.12 * time * cos(1.22 + 0.1424 * ball_i))) * 0.77 * D;  // dip into the floor
            REAL ballz = cos(ball_i + 1.32 * time * 0.1*sin((0.92 + 0.53 * ball_i))) * 0.27 * D + 0.5;
            REAL subtract = 12;
            REAL strength = 1.2 / ((sqrt(numblobs)- 1) / 4 + 1);
            mc.addBall(ballx, bally, ballz, strength, subtract, scale);
        }
    }
    void addBall( REAL ballx, REAL bally, REAL ballz, REAL strength, REAL subtract, REAL scale);
    void addPlaneX( REAL strength, REAL subtract );
    void addPlaneZ( REAL strength, REAL subtract );
    void addPlaneY( REAL strength, REAL subtract );

};




inline
void MarchingCubes::addBall(
        REAL ballx, REAL bally, REAL ballz,
        REAL strength, REAL subtract, REAL scale) {
    // Solves this equation:
    // 1.0 / (0.000001 + radius^2) * strength - subtract = 0
    REAL radius = this->resolution * sqrt(strength / subtract);

    REAL
        zs = ballz * this->resolution / scale,
        ys = bally * this->resolution / scale,
        xs = ballx * this->resolution / scale;

    int min_zi = floor( zs - radius ); if ( min_zi < 1 ) min_zi = 1;
    int max_zi = floor( zs + radius ); if ( max_zi > this->resolution - 1 ) max_zi = this->resolution - 1;
    int min_yi = floor( ys - radius ); if ( min_yi < 1 ) min_yi = 1;
    int max_yi = floor( ys + radius ); if ( max_yi > this->resolution - 1 ) max_yi = this->resolution - 1;
    int min_xi = floor( xs - radius ); if ( min_xi < 1  ) min_xi = 1;
    int max_xi = floor( xs + radius ); if ( max_xi > this->resolution - 1 ) max_xi = this->resolution - 1;


    // Don't polygonize_cube in the outer layer because normals aren't
    // well-defined there.

    // var x, y, z, y_offset, z_offset, fx, fy, fz, fz2, fy2, val;
    int x, y, z;
    REAL fx, fy, fz, fz2, fy2, val;  //Does doing like this make it faster?
    int y_offset, z_offset;

    for ( z = min_zi; z < max_zi; z++ ) {

        z_offset = this->size2 * z,
        fz = z / (REAL)this->resolution - ballz,
        fz2 = fz * fz;

        for ( y = min_yi; y < max_yi; y++ ) {

            y_offset = z_offset + this->resolution * y;
            fy = y / (REAL)this->resolution - bally;
            fy2 = fy * fy;

            for ( x = min_xi; x < max_xi; x++ ) {

                fx = x / (REAL)this->resolution - ballx;
                val = strength / ( (REAL)0.000001 + fx * fx + fy2 + fz2 ) - subtract;
                if ( val > 0.0 ) this->field[ y_offset + x ] += val / 100;
            }
        }
    }
}

void MarchingCubes::addPlaneX(REAL strength, REAL subtract ) {
    int x, y, z;
    REAL val;
    REAL xx, xdiv;
    int cxy;

    // cache attribute lookups
    int yd = this->yd;
    int resolution = this->resolution;
    int zd = this->zd;
    array1d& field = this->field;
    REAL dist = resolution * sqrt(strength / (REAL)subtract);

    if ( dist > resolution ) dist = resolution;
    for ( x = 0; x < dist; x++ ) {
        xdiv = x / (REAL)resolution;
        xx = xdiv * xdiv;
        val = strength / (REAL)( 0.0001 + xx ) - subtract;
        if ( val > 0.0 ) {
            for ( y = 0; y < resolution; y++ ) {
                cxy = x + y * yd;
                for ( z = 0; z < resolution; z++ ) {
                    field[ zd * z + cxy ] += val;
                }
            }
        }
    }
}


void MarchingCubes::addPlaneY(REAL strength, REAL subtract ) {
    int x, y, z;
    REAL yy;
    REAL val;
    REAL ydiv;
    int cy;
    int cxy;

    // cache attribute lookups
    int resolution = this->resolution;
    int yd = this->yd;
    int zd = this->zd;
    array1d& field = this->field;
    REAL dist = resolution * sqrt(strength / subtract);

    if ( dist > resolution ) dist = resolution;

    for ( y = 0; y < dist; y++ ) {
        ydiv = y / (REAL)resolution;
        yy = ydiv * ydiv;
        val = strength / (REAL)( 0.0001 + yy ) - subtract;
        if ( val > 0.0 ) {
            cy = y * yd;
            for ( x = 0; x < resolution; x++ ) {
                cxy = cy + x;
                for ( z = 0; z < resolution; z++ )
                    field[ zd * z + cxy ] += val;
            }
        }
    }
}

void MarchingCubes::addPlaneZ( REAL strength, REAL subtract )
{
    int x, y, z;
    REAL zz, val, zdiv;
    int cz, cyz;

    // cache attribute lookups
    int resolution = this->resolution;
    int yd = this->yd;
    int zd = this->zd;
    array1d& field = this->field;
    REAL dist = resolution * sqrt( strength / subtract );

    if ( dist > resolution ) dist = resolution;
    for ( z = 0; z < dist; z++ ) {
        zdiv = z / (REAL)resolution;
        zz = zdiv * zdiv;
        val = strength / (REAL)( 0.0001 + zz ) - subtract;
        if ( val > 0.0 ) {
            cz = zd * z;
            for ( y = 0; y < resolution; y++ ) {
                cyz = cz + y * yd;
                for ( x = 0; x < resolution; x++ )
                    field[ cyz + x ] += val;
            }
        }
    }
}