// This work is licensed under the Creative Commons Attribution-ShareAlike 4.0
// International Public License.
// Visit https://creativecommons.org/licenses/by-sa/4.0/.
// SPDX-License-Identifier: CC-BY-SA-4.0
//     povr histogram.pov +w900 +h450 +a0.05 +am2 +r3 +p
//
// Change function just below to look at a histogram of the scalar values
// returned. Set up runs all in the parser and using sum() which is slow. Expect
// 15 - 45 seconds to parse and render with 256 buckets of 100k samples. On the
// whole probably not too different than running a set of external programs to
// get the histogram - and many of those are limited with respect to values
// considered.
#version unofficial 3.8; // povr
#if (file_exists("version.inc"))
    #include "version.inc"
#end
#if (!defined(Fork_povr))
    #error "This POV-Ray SDL code requires the povr fork."
#end

#declare Fn00 = function { pattern { wrinkles noise_generator 1 scale 1/pi } }

#declare Samples = int(1e5);
#declare Buckets = int(256)+1;   // Bucket value inside int()
#declare LoVal   = -0.0-(1e-5);  // Usually 0 or function_interval -1
#declare HiVal   = +1.0+(1e-5);
#declare Range   = HiVal - LoVal;
#declare Delta   = Range / (Buckets - 1);
#declare Xl = -1.0; #declare Xh = 1.0;
#declare Yl = -1.0; #declare Yh = 1.0;
#declare Normalize = yes;  // If yes, scales max bucket % to y = 1.

#include "functions.inc"
#declare Fn01 = function (_row,_col,_lo,_hi) {
    select ((Fn00(_row,_col,0)>=_lo & Fn00(_row,_col,0)<_hi),
        0,0,1)
}
#declare FnSum =
    function (x,_lo,_delta,_range,_samples,_samplesSqrt) {
        sum(I,0,_samples-1,
            Fn01(f_range2range(mod(I,_samplesSqrt),0.0,_samplesSqrt,Xl,Xh),
                 f_range2range(int(I/_samplesSqrt),0.0,_samplesSqrt,Yl,Yh),
                _lo+(x*_range),_lo+(x*_range)+_delta)
        ) / _samples
    }

#declare Pigment00 = pigment {
    image_map {
        function Buckets, 1 {
            FnSum(x,LoVal,Delta,Range,Samples,sqrt(Samples))
        }
        map_type 0
        once
    }
}
#declare FnHist = function { pigment { Pigment00 } }


//--------- Create the scene to render the histogram

#declare MaxBucketPct = 0.0;
#declare LocalDelta = (1.0/(Buckets-1))*(HiVal-LoVal);
#declare Union00 = union {
#for (II,0,Buckets-2,1)
    #declare JJ = II / (Buckets-1);
    #declare LocalX = f_range2range(JJ,0.0,1.0,LoVal,HiVal);
    #if (FnHist(JJ,0,0).red > 0.0)
        cylinder { <LocalX+LocalDelta/2.0,0.0,0.0>,
            <LocalX+LocalDelta/2.0,FnHist(JJ,0,0).red,0>,
            LocalDelta/2.0
        }
    #end
    #debug concat(str(LocalX,8,6)," to ",
        str(LocalX+LocalDelta,8,6)," --> ",
        str(FnHist(JJ,0,0).red,8,6),"\n")
    #declare MaxBucketPct = max(FnHist(JJ,0,0).red,MaxBucketPct);
#end
#if (Normalize)
    scale <1,1/MaxBucketPct,1>
#end
}
#debug concat("Max bucket % -> ",str(MaxBucketPct,8,6),"\n")

global_settings { assumed_gamma 1 }
#declare Grey05 = srgb <0.05,0.05,0.05>;
background { Grey05 }
#declare VarOrthoMult =
    3.0/max(image_width/image_height,image_height/image_width);
#declare Camera01z = camera {
    orthographic
    location <0,0.5,-2>
    direction z
    right VarOrthoMult*x*max(1,image_width/image_height)
    up VarOrthoMult*y*max(1,image_height/image_width)
}

object { Union00 pigment { rgb <0,1,0> } finish { emission 1 } }
camera { Camera01z }

// Mark origin and +-1
#declare CylY = cylinder {
    -1*y, 0*y, 0.005
    pigment { rgb <0,0.5,1> } finish { emission 1 }
}
object { CylY }
object { CylY translate -1*x }
object { CylY translate +1*x }