/*
 * Copyright 2011-2013 Blender Foundation
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "node_noise.h"
#include "stdcycles.h"
#include "vector2.h"
#include "vector4.h"

#define vector3 point

/* 1D Musgrave fBm
 *
 * H: fractal increment parameter
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 *
 * from "Texturing and Modelling: A procedural approach"
 */

float noise_musgrave_fBm_1d(float co, float H, float lacunarity, float octaves)
{
  float p = co;
  float value = 0.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value += safe_snoise(p) * pwr;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * safe_snoise(p) * pwr;
  }

  return value;
}

/* 1D Musgrave Multifractal
 *
 * H: highest fractal dimension
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 */

float noise_musgrave_multi_fractal_1d(float co, float H, float lacunarity, float octaves)
{
  float p = co;
  float value = 1.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value *= (pwr * safe_snoise(p) + 1.0);
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value *= (rmd * pwr * safe_snoise(p) + 1.0); /* correct? */
  }

  return value;
}

/* 1D Musgrave Heterogeneous Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hetero_terrain_1d(
    float co, float H, float lacunarity, float octaves, float offset)
{
  float p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  /* first unscaled octave of function; later octaves are scaled */
  float value = offset + safe_snoise(p);
  p *= lacunarity;

  for (int i = 1; i < (int)octaves; i++) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += increment;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += rmd * increment;
  }

  return value;
}

/* 1D Hybrid Additive/Multiplicative Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hybrid_multi_fractal_1d(
    float co, float H, float lacunarity, float octaves, float offset, float gain)
{
  float p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float value = safe_snoise(p) + offset;
  float weight = gain * value;
  p *= lacunarity;

  for (int i = 1; (weight > 0.001) && (i < (int)octaves); i++) {
    if (weight > 1.0) {
      weight = 1.0;
    }

    float signal = (safe_snoise(p) + offset) * pwr;
    pwr *= pwHL;
    value += weight * signal;
    weight *= gain * signal;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * ((safe_snoise(p) + offset) * pwr);
  }

  return value;
}

/* 1D Ridged Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_ridged_multi_fractal_1d(
    float co, float H, float lacunarity, float octaves, float offset, float gain)
{
  float p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float signal = offset - fabs(safe_snoise(p));
  signal *= signal;
  float value = signal;
  float weight = 1.0;

  for (int i = 1; i < (int)octaves; i++) {
    p *= lacunarity;
    weight = clamp(signal * gain, 0.0, 1.0);
    signal = offset - fabs(safe_snoise(p));
    signal *= signal;
    signal *= weight;
    value += signal * pwr;
    pwr *= pwHL;
  }

  return value;
}

/* 2D Musgrave fBm
 *
 * H: fractal increment parameter
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 *
 * from "Texturing and Modelling: A procedural approach"
 */

float noise_musgrave_fBm_2d(vector2 co, float H, float lacunarity, float octaves)
{
  vector2 p = co;
  float value = 0.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value += safe_snoise(p) * pwr;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * safe_snoise(p) * pwr;
  }

  return value;
}

/* 2D Musgrave Multifractal
 *
 * H: highest fractal dimension
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 */

float noise_musgrave_multi_fractal_2d(vector2 co, float H, float lacunarity, float octaves)
{
  vector2 p = co;
  float value = 1.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value *= (pwr * safe_snoise(p) + 1.0);
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value *= (rmd * pwr * safe_snoise(p) + 1.0); /* correct? */
  }

  return value;
}

/* 2D Musgrave Heterogeneous Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hetero_terrain_2d(
    vector2 co, float H, float lacunarity, float octaves, float offset)
{
  vector2 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  /* first unscaled octave of function; later octaves are scaled */
  float value = offset + safe_snoise(p);
  p *= lacunarity;

  for (int i = 1; i < (int)octaves; i++) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += increment;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += rmd * increment;
  }

  return value;
}

/* 2D Hybrid Additive/Multiplicative Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hybrid_multi_fractal_2d(
    vector2 co, float H, float lacunarity, float octaves, float offset, float gain)
{
  vector2 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float value = safe_snoise(p) + offset;
  float weight = gain * value;
  p *= lacunarity;

  for (int i = 1; (weight > 0.001) && (i < (int)octaves); i++) {
    if (weight > 1.0) {
      weight = 1.0;
    }

    float signal = (safe_snoise(p) + offset) * pwr;
    pwr *= pwHL;
    value += weight * signal;
    weight *= gain * signal;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * ((safe_snoise(p) + offset) * pwr);
  }

  return value;
}

/* 2D Ridged Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_ridged_multi_fractal_2d(
    vector2 co, float H, float lacunarity, float octaves, float offset, float gain)
{
  vector2 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float signal = offset - fabs(safe_snoise(p));
  signal *= signal;
  float value = signal;
  float weight = 1.0;

  for (int i = 1; i < (int)octaves; i++) {
    p *= lacunarity;
    weight = clamp(signal * gain, 0.0, 1.0);
    signal = offset - fabs(safe_snoise(p));
    signal *= signal;
    signal *= weight;
    value += signal * pwr;
    pwr *= pwHL;
  }

  return value;
}

/* 3D Musgrave fBm
 *
 * H: fractal increment parameter
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 *
 * from "Texturing and Modelling: A procedural approach"
 */

float noise_musgrave_fBm_3d(vector3 co, float H, float lacunarity, float octaves)
{
  vector3 p = co;
  float value = 0.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value += safe_snoise(p) * pwr;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * safe_snoise(p) * pwr;
  }

  return value;
}

/* 3D Musgrave Multifractal
 *
 * H: highest fractal dimension
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 */

float noise_musgrave_multi_fractal_3d(vector3 co, float H, float lacunarity, float octaves)
{
  vector3 p = co;
  float value = 1.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value *= (pwr * safe_snoise(p) + 1.0);
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value *= (rmd * pwr * safe_snoise(p) + 1.0); /* correct? */
  }

  return value;
}

/* 3D Musgrave Heterogeneous Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hetero_terrain_3d(
    vector3 co, float H, float lacunarity, float octaves, float offset)
{
  vector3 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  /* first unscaled octave of function; later octaves are scaled */
  float value = offset + safe_snoise(p);
  p *= lacunarity;

  for (int i = 1; i < (int)octaves; i++) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += increment;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += rmd * increment;
  }

  return value;
}

/* 3D Hybrid Additive/Multiplicative Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hybrid_multi_fractal_3d(
    vector3 co, float H, float lacunarity, float octaves, float offset, float gain)
{
  vector3 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float value = safe_snoise(p) + offset;
  float weight = gain * value;
  p *= lacunarity;

  for (int i = 1; (weight > 0.001) && (i < (int)octaves); i++) {
    if (weight > 1.0) {
      weight = 1.0;
    }

    float signal = (safe_snoise(p) + offset) * pwr;
    pwr *= pwHL;
    value += weight * signal;
    weight *= gain * signal;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * ((safe_snoise(p) + offset) * pwr);
  }

  return value;
}

/* 3D Ridged Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_ridged_multi_fractal_3d(
    vector3 co, float H, float lacunarity, float octaves, float offset, float gain)
{
  vector3 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float signal = offset - fabs(safe_snoise(p));
  signal *= signal;
  float value = signal;
  float weight = 1.0;

  for (int i = 1; i < (int)octaves; i++) {
    p *= lacunarity;
    weight = clamp(signal * gain, 0.0, 1.0);
    signal = offset - fabs(safe_snoise(p));
    signal *= signal;
    signal *= weight;
    value += signal * pwr;
    pwr *= pwHL;
  }

  return value;
}

/* 4D Musgrave fBm
 *
 * H: fractal increment parameter
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 *
 * from "Texturing and Modelling: A procedural approach"
 */

float noise_musgrave_fBm_4d(vector4 co, float H, float lacunarity, float octaves)
{
  vector4 p = co;
  float value = 0.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value += safe_snoise(p) * pwr;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * safe_snoise(p) * pwr;
  }

  return value;
}

/* 4D Musgrave Multifractal
 *
 * H: highest fractal dimension
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 */

float noise_musgrave_multi_fractal_4d(vector4 co, float H, float lacunarity, float octaves)
{
  vector4 p = co;
  float value = 1.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);

  for (int i = 0; i < (int)octaves; i++) {
    value *= (pwr * safe_snoise(p) + 1.0);
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value *= (rmd * pwr * safe_snoise(p) + 1.0); /* correct? */
  }

  return value;
}

/* 4D Musgrave Heterogeneous Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hetero_terrain_4d(
    vector4 co, float H, float lacunarity, float octaves, float offset)
{
  vector4 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  /* first unscaled octave of function; later octaves are scaled */
  float value = offset + safe_snoise(p);
  p *= lacunarity;

  for (int i = 1; i < (int)octaves; i++) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += increment;
    pwr *= pwHL;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    float increment = (safe_snoise(p) + offset) * pwr * value;
    value += rmd * increment;
  }

  return value;
}

/* 4D Hybrid Additive/Multiplicative Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hybrid_multi_fractal_4d(
    vector4 co, float H, float lacunarity, float octaves, float offset, float gain)
{
  vector4 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float value = safe_snoise(p) + offset;
  float weight = gain * value;
  p *= lacunarity;

  for (int i = 1; (weight > 0.001) && (i < (int)octaves); i++) {
    if (weight > 1.0) {
      weight = 1.0;
    }

    float signal = (safe_snoise(p) + offset) * pwr;
    pwr *= pwHL;
    value += weight * signal;
    weight *= gain * signal;
    p *= lacunarity;
  }

  float rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    value += rmd * ((safe_snoise(p) + offset) * pwr);
  }

  return value;
}

/* 4D Ridged Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_ridged_multi_fractal_4d(
    vector4 co, float H, float lacunarity, float octaves, float offset, float gain)
{
  vector4 p = co;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;

  float signal = offset - fabs(safe_snoise(p));
  signal *= signal;
  float value = signal;
  float weight = 1.0;

  for (int i = 1; i < (int)octaves; i++) {
    p *= lacunarity;
    weight = clamp(signal * gain, 0.0, 1.0);
    signal = offset - fabs(safe_snoise(p));
    signal *= signal;
    signal *= weight;
    value += signal * pwr;
    pwr *= pwHL;
  }

  return value;
}

shader node_musgrave_texture(
    int use_mapping = 0,
    matrix mapping = matrix(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
    string musgrave_type = "fBM",
    string dimensions = "3D",
    point Vector = P,
    float W = 0.0,
    float Dimension = 2.0,
    float Scale = 5.0,
    float Detail = 2.0,
    float Lacunarity = 2.0,
    float Offset = 0.0,
    float Gain = 1.0,
    output float Fac = 0.0)
{
  float dimension = max(Dimension, 1e-5);
  float octaves = clamp(Detail, 0.0, 16.0);
  float lacunarity = max(Lacunarity, 1e-5);

  vector3 s = Vector;

  if (use_mapping)
    s = transform(mapping, s);

  if (dimensions == "1D") {
    float p = W * Scale;
    if (musgrave_type == "multifractal") {
      Fac = noise_musgrave_multi_fractal_1d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "fBM") {
      Fac = noise_musgrave_fBm_1d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "hybrid_multifractal") {
      Fac = noise_musgrave_hybrid_multi_fractal_1d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "ridged_multifractal") {
      Fac = noise_musgrave_ridged_multi_fractal_1d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "hetero_terrain") {
      Fac = noise_musgrave_hetero_terrain_1d(p, dimension, lacunarity, octaves, Offset);
    }
    else {
      Fac = 0.0;
    }
  }
  else if (dimensions == "2D") {
    vector2 p = vector2(s[0], s[1]) * Scale;
    if (musgrave_type == "multifractal") {
      Fac = noise_musgrave_multi_fractal_2d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "fBM") {
      Fac = noise_musgrave_fBm_2d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "hybrid_multifractal") {
      Fac = noise_musgrave_hybrid_multi_fractal_2d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "ridged_multifractal") {
      Fac = noise_musgrave_ridged_multi_fractal_2d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "hetero_terrain") {
      Fac = noise_musgrave_hetero_terrain_2d(p, dimension, lacunarity, octaves, Offset);
    }
    else {
      Fac = 0.0;
    }
  }
  else if (dimensions == "3D") {
    vector3 p = s * Scale;
    if (musgrave_type == "multifractal") {
      Fac = noise_musgrave_multi_fractal_3d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "fBM") {
      Fac = noise_musgrave_fBm_3d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "hybrid_multifractal") {
      Fac = noise_musgrave_hybrid_multi_fractal_3d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "ridged_multifractal") {
      Fac = noise_musgrave_ridged_multi_fractal_3d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "hetero_terrain") {
      Fac = noise_musgrave_hetero_terrain_3d(p, dimension, lacunarity, octaves, Offset);
    }
    else {
      Fac = 0.0;
    }
  }
  else if (dimensions == "4D") {
    vector4 p = vector4(s[0], s[1], s[2], W) * Scale;
    if (musgrave_type == "multifractal") {
      Fac = noise_musgrave_multi_fractal_4d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "fBM") {
      Fac = noise_musgrave_fBm_4d(p, dimension, lacunarity, octaves);
    }
    else if (musgrave_type == "hybrid_multifractal") {
      Fac = noise_musgrave_hybrid_multi_fractal_4d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "ridged_multifractal") {
      Fac = noise_musgrave_ridged_multi_fractal_4d(
          p, dimension, lacunarity, octaves, Offset, Gain);
    }
    else if (musgrave_type == "hetero_terrain") {
      Fac = noise_musgrave_hetero_terrain_4d(p, dimension, lacunarity, octaves, Offset);
    }
    else {
      Fac = 0.0;
    }
  }
  else {
    Fac = 0.0;
  }
}