/**************************************************************************
 *{@C
 *      Copyright:      1988-2024 Paul Obermeier (obermeier@poSoft.de)
 *
 *      Module:         ImageProcessing
 *      Filename:       IP_Fractal.c
 *
 *      Author:         Paul Obermeier
 *
 *      Description:    Functions to generate fractal noise, Julia and
 *                      Mandelbrot sets.
 *
 *      Additional documentation:
 *                      None.
 *
 *      Exported functions:
 *                      IP_GenNoise
 *                      IP_GenJuliaSet
 *                      IP_GenMandelbrot
 *
 **************************************************************************/

#include <stdio.h>

#include "UT_Compat.h"
#include "UT_Macros.h"
#include "UT_Error.h"
#include "UT_Math.h"

#define  USE_IMG_MACROS
#include "FF_Image.h"

#include "IP_Image.h"
#include "IP_ImagePrivate.h"


/* Set channel "chn" of pixel (x, y) in image "img" to "v". */

static void set_pixel (IP_ImageId img, Int32 chn, Int32 x, Int32 y, Float32 v)
{
    v = UT_MAX (0.0, UT_MIN (v, 1.0));
    switch (img->pimg.channel[chn]) {
        case FF_ImgFmtTypeNone:
            break;
        case FF_ImgFmtTypeUByte:
            *FF_ImgUBytePixel (&img->pimg, chn, x, y) = (UInt8) (v * 255.0);
            break;
        case FF_ImgFmtTypeFloat:
            *FF_ImgFloatPixel (&img->pimg, chn, x, y) = v;
            break;
    }
    return;
}

static void juliaset (IP_ImageId img,
                      Float32 cx, Float32 cy, Float32 cenx, Float32 ceny, Float32 rlen, 
                      Int32 numit, UInt8 rmap[], UInt8 gmap[], UInt8 bmap[])
{
    Float64 rmin, imin, dr, di, ilen;
    Float64 zr2, zi2;
    Float64 sx, sy, zx, zy, znewx, znewy;
    Int32 k, x, y;

    /* Do some precalculation. */
    ilen = rlen * (Float64)img->pimg.desc.height/(Float64)img->pimg.desc.width;
    rmin = cenx - rlen * 0.5;
    imin = ceny - ilen * 0.5;
    dr = rlen / img->pimg.desc.width;
    di = ilen / img->pimg.desc.height;
    
    sy = imin;
    for (y = 0; y < img->pimg.desc.height; y++) {
        sx = rmin;
        for (x=img->pimg.desc.width-1; x>=0; x--) {
            zx = sx; 
            zy = sy;
            k = 0;
            do {
                zr2 = zx * zx;
                zi2 = zy * zy;
                znewx = zr2 - zi2 + cx;
                znewy = 2.0*zx*zy + cy;
                zx = znewx;
                zy = znewy;
            } while ((zx*zx + zy*zy)<4.0 && ++k<numit);

            if (k >= numit) {
                /* Point belongs to the set */
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeRed,   x, y) = 0;
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeGreen, x, y) = 0;
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeBlue,  x, y) = 0;

            } else {
                /* Point doesn't belong to the set */
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeRed,   x, y) = rmap[k % 255];
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeGreen, x, y) = gmap[k % 255];
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeBlue,  x, y) = bmap[k % 255];
            }
            sx += dr;
        }
        sy += di;
    }
    return;
}

static void mandelbrot (IP_ImageId img, Float32 cenx, Float32 ceny, Float32 rlen, 
                        Int32 numit, UInt8 rmap[], UInt8 gmap[], UInt8 bmap[])
{
    Float64 rmin, imin, dr, di, ilen;
    Float64 sx, sy, zx, zy, znewx, znewy;
    Int32 k, x, y;
    Float64 zr2, zi2;

    ilen = rlen * (Float64)img->pimg.desc.height/(Float64)img->pimg.desc.width;
    rmin = cenx - rlen * 0.5;
    imin = ceny - ilen * 0.5;
    dr = rlen / img->pimg.desc.width;
    di = ilen / img->pimg.desc.height;
    
    sy = imin;
    for (y=0; y<img->pimg.desc.height; y++) {
        sx = rmin;
        for (x=img->pimg.desc.width-1; x>=0; x--) {
            zx = 0.0; 
            zy = 0.0;
            k = 0;
            do {
                zr2 = zx*zx;
                zi2 = zy*zy;
                znewx = zr2 - zi2 + sx;
                znewy = 2.0*zx*zy + sy;
                zx = znewx;
                zy = znewy;
            } while ((zr2 + zi2)<100.0 && ++k<numit);

            if (k >= numit) {
                /* Point belongs to the set */
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeRed,   x, y) = 0;
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeGreen, x, y) = 0;
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeBlue,  x, y) = 0;
            } else {
                /* Point doesn't belong to the set */
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeRed,   x, y) = rmap[k % 255];
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeGreen, x, y) = gmap[k % 255];
                *FF_ImgUBytePixel (&img->pimg, FF_ImgChanTypeBlue,  x, y) = bmap[k % 255];
            }
            sx += dr;
        }
        sy += di;
    }
    return;
}

/* Fill channel "chn" of image "img" with fractal noise. */

static void fractal_noise
        (IP_ImageId img, FF_ImgChanType chn,
         Int32 seed, Int32 p, Int32 coh, Float32 z)
{
    UT_ImgNoiseGen gen;
    Int32       imin, i, x, y;
    Float32     ff[UT_ImgNoiseBits + 1], aa[UT_ImgNoiseBits + 1], f, a, v;

    UT_ImgMakeNoiseGen (seed, &gen);
    a = 1.0 / (1.0 - 1.0 / (1 << (p - coh)));
    imin = coh + 1;
    for (i = imin; i <= p; ++i) {
        ff[i] = 1.0 / (1 << (p - i));
        aa[i] = a / (1 << (i - coh));
    }
    for (y = 0; y < img->pimg.desc.height; ++y) {
        for (x = 0; x < img->pimg.desc.width; ++x) {
            v = 0.0;
            for (i = imin; i <= p; ++i) {
                f = ff[i];
                v += (UT_ImgNoise (&gen, i, x * f, y * f, z * f) - 0.5) * aa[i];
            }
            set_pixel (img, chn, x, y, v + 0.5);
        }
    }
    return;
}

/***************************************************************************
 *[@e
 *      Name:           IP_GenNoise
 *
 *      Usage:          Fill one channel of an image with fractal noise.
 *
 *      Synopsis:       UT_Bool IP_GenNoise
 *                              (IP_ImageId img, FF_ImgChanType chn,
 *                              Int32 seed, Int32 p, Int32 coh, Float32 z)
 *
 *      Description:    Channel "chn" of image "img" is filled with fractal
 *                      noise: A pseudo-random value in the range from 0.0 to
 *                      1.0 is assigned to every pixel so that pixels, which
 *                      are close to each other, are likely to receive similar
 *                      values.
 *                      The fractal noise pattern is generated from a seed
 *                      value, "seed", for a pseudo-random number generator.
 *                      The pattern is three-dimensional, and the values
 *                      written into "img" represent a two-dimensional slice
 *                      of it. The position of the slice is selected by the
 *                      value of "z", so that sequences of slices from the
 *                      same noise pattern can be generated by stepping
 *                      through different "z" values.
 *                      The noise pattern is periodic in x, y, (i.e.
 *                      horizontal and vertical) and z direction with a
 *                      period length of (1 << p) pixels.
 *                      The "coh" parameter determines how coherent the noise
 *                      pattern is. The pattern becomes more and more
 *                      incoherent when "coh" increases. "Coh" should be in
 *                      the range from 0 to p-1.
 *
 *      Return value:   UT_True if successful, else UT_False.
 *
 *      See also:
 *
 ***************************************************************************/

UT_Bool IP_GenNoise
        (IP_ImageId img, FF_ImgChanType chn,
         Int32 seed, Int32 p, Int32 coh, Float32 z)
{
    if (chn < 0 || chn >= FF_NumImgChanTypes) {
         UT_ErrSetNum (UT_ErrParamInvalid, str_imgchan, 
                       FF_ImgGetChanTypeName (chn));
        return UT_False;
    }
    if (p < 1) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_inotless, str_periodlength, 1);
        return UT_False;
    }
    p = UT_MIN (p, UT_ImgNoiseBits);
    coh = UT_MAX (0, UT_MIN (coh, p - 1));
    fractal_noise (img, chn, seed, p, coh, z);
    return UT_True;
}

/***************************************************************************
 *[@e
 *      Name:           IP_GenJuliaSet
 *
 *      Usage:          Generate an image of a Julia set.       
 *
 *      Synopsis:       UT_Bool IP_GenJuliaSet
 *                              (IP_ImageId img, Float32 cr, Float32 ci,
 *                               Float32 rcen, Float32 icen, Float32 rlen,
 *                               Int32 numit, UInt8 rmap[], UInt8 gmap[], 
 *                               UInt8 bmap[])
 *
 *      Description:    img:            Image with red, green and blue 
 *                                      channels in FF_ImgFmtTypeUByte format.
 *                      cr, ci:         Real and imaginary part of the additive
 *                                      constant.
 *                      rcen, icen:     Real and imaginary part of the center
 *                                      of the screen.
 *                      rlen:           Length of the displayed interval in
 *                                      real direction.
 *                      numit:          Number of iterations per pixel.
 *                      rmap, gmap, bmap:
 *                                      Array of 256 color values for the red,
 *                                      green and blue channels.
 *
 *      Return value:   UT_True if successful, else UT_False.
 *
 *      See also:       IP_GenMandelbrot
 *
 ***************************************************************************/

UT_Bool IP_GenJuliaSet
        (IP_ImageId img,
         Float32 cr, Float32 ci, Float32 rcen, Float32 icen, Float32 rlen, 
         Int32 numit, UInt8 rmap[], UInt8 gmap[], UInt8 bmap[])
{
    Float32 cx, cy, cenx, ceny;

    if (img->pimg.channel[FF_ImgChanTypeRed] != FF_ImgFmtTypeUByte ||
        img->pimg.channel[FF_ImgChanTypeGreen] != FF_ImgFmtTypeUByte ||
        img->pimg.channel[FF_ImgChanTypeBlue] != FF_ImgFmtTypeUByte) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_rgb_ubyte);
        return UT_False;
    }
    cx = cr;
    cy = ci;
    cenx = rcen;
    ceny = icen;
    juliaset (img, cx, cy, cenx, ceny, rlen, numit, rmap, gmap, bmap);
    return UT_True;
}

/***************************************************************************
 *[@e
 *      Name:           IP_GenMandelbrot
 *
 *      Usage:          Generate an image of a Mandelbrot set.  
 *
 *      Synopsis:       UT_Bool IP_GenMandelbrot
 *                              (IP_ImageId img, Float32 rcen, Float32 icen, 
 *                               Float32 rlen, Int32 numit, 
 *                               UInt8 rmap[], UInt8 gmap[], UInt8 bmap[])
 *
 *      Description:    img:            Image with red, green and blue 
 *                                      channels in FF_ImgFmtTypeUByte format.
 *                      rcen, icen:     Real and imaginary part of the center
 *                                      of the screen.
 *                      rlen:           Length of the displayed interval in
 *                                      real direction.
 *                      numit:          Number of iterations per pixel.
 *                      rmap, gmap, bmap:
 *                                      Array of 256 color values for the red,
 *                                      green and blue channels.
 *
 *      Return value:   UT_True if successful, else UT_False.
 *
 *      See also:       IP_GenJuliaSet
 *
 ***************************************************************************/

UT_Bool IP_GenMandelbrot
        (IP_ImageId img, Float32 rcen, Float32 icen, Float32 rlen, 
         Int32 numit, UInt8 rmap[], UInt8 gmap[], UInt8 bmap[])
{
    Float32 cenx, ceny;

    if (img->pimg.channel[FF_ImgChanTypeRed] != FF_ImgFmtTypeUByte ||
        img->pimg.channel[FF_ImgChanTypeGreen] != FF_ImgFmtTypeUByte ||
        img->pimg.channel[FF_ImgChanTypeBlue] != FF_ImgFmtTypeUByte) {
        UT_ErrSetNum (UT_ErrParamInvalid, str_rgb_ubyte);
        return UT_False;
    }
    cenx = rcen;
    ceny = icen;
    mandelbrot (img, cenx, ceny, rlen, numit, rmap, gmap, bmap);
    return UT_True;
}
