/*
 * 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 "stdosl.h"
#include "node_texture.h"

/* 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(point p, float H, float lacunarity, float octaves)
{
	float rmd;
	float value = 0.0;
	float pwr = 1.0;
	float pwHL = pow(lacunarity, -H);
	int i;

	for (i = 0; i < (int)octaves; i++) {
		value += safe_noise(p, "signed") * pwr;
		pwr *= pwHL;
		p *= lacunarity;
	}

	rmd = octaves - floor(octaves);
	if (rmd != 0.0)
		value += rmd * safe_noise(p, "signed") * pwr;

	return value;
}

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

float noise_musgrave_multi_fractal(point p, float H, float lacunarity, float octaves)
{
	float rmd;
	float value = 1.0;
	float pwr = 1.0;
	float pwHL = pow(lacunarity, -H);
	int i;

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

	rmd = octaves - floor(octaves);
	if (rmd != 0.0)
		value *= (rmd * pwr * safe_noise(p, "signed") + 1.0); /* correct? */

	return value;
}

/* 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(point p, float H, float lacunarity, float octaves, float offset)
{
	float value, increment, rmd;
	float pwHL = pow(lacunarity, -H);
	float pwr = pwHL;
	int i;

	/* first unscaled octave of function; later octaves are scaled */
	value = offset + safe_noise(p, "signed");
	p *= lacunarity;

	for (i = 1; i < (int)octaves; i++) {
		increment = (safe_noise(p, "signed") + offset) * pwr * value;
		value += increment;
		pwr *= pwHL;
		p *= lacunarity;
	}

	rmd = octaves - floor(octaves);
	if (rmd != 0.0) {
		increment = (safe_noise(p, "signed") + offset) * pwr * value;
		value += rmd * increment;
	}

	return value;
}

/* 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(point p, float H, float lacunarity,
                                          float octaves, float offset, float gain)
{
	float result, signal, weight, rmd;
	float pwHL = pow(lacunarity, -H);
	float pwr = pwHL;
	int i;

	result = safe_noise(p, "signed") + offset;
	weight = gain * result;
	p *= lacunarity;

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

		signal = (safe_noise(p, "signed") + offset) * pwr;
		pwr *= pwHL;
		result += weight * signal;
		weight *= gain * signal;
		p *= lacunarity;
	}

	rmd = octaves - floor(octaves);
	if (rmd != 0.0)
		result += rmd * ((safe_noise(p, "signed") + offset) * pwr);

	return result;
}

/* 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(point p, float H, float lacunarity,
                                          float octaves, float offset, float gain)
{
	float result, signal, weight;
	float pwHL = pow(lacunarity, -H);
	float pwr = pwHL;
	int i;

	signal = offset - fabs(safe_noise(p, "signed"));
	signal *= signal;
	result = signal;
	weight = 1.0;

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

	return result;
}

/* Shader */

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 Type = "fBM",
	float Dimension = 2.0,
	float Lacunarity = 1.0,
	float Detail = 2.0,
	float Offset = 0.0,
	float Gain = 1.0,
	float Scale = 5.0,
	point Vector = P,
	output float Fac = 0.0,
	output color Color = 0.0)
{
	float dimension = max(Dimension, 1e-5);
	float octaves = clamp(Detail, 0.0, 16.0);
	float lacunarity = max(Lacunarity, 1e-5);
	float intensity = 1.0;

	point p = Vector;

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

	p = p * Scale;

	if (Type == "Multifractal")
		Fac = intensity * noise_musgrave_multi_fractal(p, dimension, lacunarity, octaves);
	else if (Type == "fBM")
		Fac = intensity * noise_musgrave_fBm(p, dimension, lacunarity, octaves);
	else if (Type == "Hybrid Multifractal")
		Fac = intensity * noise_musgrave_hybrid_multi_fractal(p, dimension, lacunarity, octaves, Offset, Gain);
	else if (Type == "Ridged Multifractal")
		Fac = intensity * noise_musgrave_ridged_multi_fractal(p, dimension, lacunarity, octaves, Offset, Gain);
	else if (Type == "Hetero Terrain")
		Fac = intensity * noise_musgrave_hetero_terrain(p, dimension, lacunarity, octaves, Offset);
	
	Color = color(Fac, Fac, Fac);
}