Thanks C0ryMcG,
I must admit I like the final result I got from Adam Majewski's code. The main drawbacks seem to be the need to calculate derivatives and the need to store all the derivatives for post processing.
I should have put superheel's slope functions in the code segment as well so you can see how they work.
double transformResultToHeight(DWORD iteration, long size)
{
return (double)iteration;
// if (size > 0L)
// return (double)iteration / (double)size;
// else return -1;
}
double getGradientX(DWORD *wpixels, int index, int size)
{
// size = ydots
int x = index % size;
double it = /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index), threshold);
if (x == 0) {
return (/*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index + 1), threshold) - it) * 2;
}
else if (x == size - 1) {
return (it - /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index - 1), threshold)) * 2;
}
else {
double diffL = it - /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index - 1), threshold);
double diffR = it - /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index + 1), threshold);
return diffL * diffR >= 0 ? 0 : diffL - diffR;
}
}
double getGradientY(DWORD *wpixels, int index, int size) {
int y = index / size;
double it = /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index), threshold);
if (y == 0) {
return (it - /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index + size), threshold)) * 2;
}
else if (y == size - 1) {
return (/*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index - size), threshold) - it) * 2;
}
else {
double diffU = it - /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index - size), threshold);
double diffD = it - /*ColorAlgorithm.*/transformResultToHeight(*(wpixels + index + size), threshold);
return diffD * diffU >= 0 ? 0 : diffD - diffU;
}
}
Here is the code written by Adam Majewski that I used to get the early results:
/*
c program: console, 1-file
normal.c
normal = Normal map effect Mandelbrot set
https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set
thx for help:
* Claude Heiland-Allen http://mathr.co.uk/blog/
--------------------------------
1. draws Mandelbrot set for Fc(z)=z*z +c
using Mandelbrot algorithm ( boolean escape time )
-------------------------------
2. technique of creating ppm file is based on the code of Claudio Rocchini
http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
create 24 bit color graphic file , portable pixmap file = PPM
see http://en.wikipedia.org/wiki/Portable_pixmap
to see the file use external application ( graphic viewer)
-----
it is example for :
https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set
-------------
compile :
gcc normal.c -lm -Wall
./a.out
*/
#include < stdio.h >
#include < math.h >
#include < complex.h >
#define M_PI 3.14159265358979323846 /* pi */
/* screen ( integer) coordinate */
int iX, iY;
const int iXmax = 10000;
const int iYmax = 10001; // for main antenna
/* world ( double) coordinate = parameter plane*/
double Cx, Cy;
const double CxMin = -2.2;
const double CxMax = 0.8;
const double CyMin = -1.5;
const double CyMax = 1.5;
/* */
double PixelWidth; //=(CxMax-CxMin)/iXmax;
double PixelHeight; // =(CyMax-CyMin)/iYmax;
/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
const int MaxColorComponentValue = 255;
FILE * fp;
char * filename = "n10000.ppm"; // https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/File:M-bump-1.png
char * comment = "# "; /* comment should start with # */
static unsigned char color[3]; // 24-bit rgb color
/* */
const int IterationMax = 2000; // N in wiki
/* bail-out value for the bailout test for escaping points
radius of circle centered ad the origin, exterior of such circle is a target set */
const double EscapeRadius = 1000; // big !!!!
double complex give_c(int iX, int iY) {
double Cx, Cy;
Cy = CyMax - iY * PixelHeight;
Cx = CxMin + iX * PixelWidth;
return Cx + Cy * I;
}
/*
The dot product of two vectors a = [a1, a2, ..., an] and b = [b1, b2, ..., bn] is defined as:[1]
d = a1b1 + a2b2
*/
double cdot(double complex a, double complex b) {
return creal(a) * creal(b) + cimag(a) * cimag(b);
}
//
// output
//
double GiveReflection(double complex C, int iMax) {
int i = 0; // iteration
double complex Z = 0.0; // initial value for iteration Z0
double complex dC = 0.0; // derivative with respect to c
double reflection = FP_ZERO; // inside
double h2 = 1.5; // height factor of the incoming light
double angle = 45.0 / 360.0; // incoming direction of light in turns
double complex v = cexp(2.0 * angle * M_PI * I); // = exp(1j*angle*2*pi/360) // unit 2D vector in this direction
// incoming light 3D vector = (v.re,v.im,h2)
double complex u;
for (i = 0; i < iMax; i++) {
dC = 2.0 * dC * Z + 1.0;
Z = Z * Z + C;
if (cabs(Z) > EscapeRadius) { // exterior of M set
u = Z / dC;
u = u / cabs(u);
reflection = cdot(u, v) + h2;
reflection = reflection / (1.0 + h2); // rescale so that t does not get bigger than 1
if (reflection < 0.0) reflection = 0.0;
break;
}
}
return reflection;
}
int compute_color(complex double c, unsigned char color[3]) {
double reflection;
unsigned char b;
// compute
reflection = GiveReflection(c, IterationMax);
if (reflection == FP_ZERO) { /* interior of Mandelbrot set = inside_color = blue */
color[0] = 0; // M_waves
color[1] = 0;
color[2] = 127;
} else // exterior of Mandelbrot set = normal
{
b = 255 * reflection;
//b = (unsigned char) (255 - 255*potential );
color[0] = b; /* Red*/
color[1] = b; /* Green */
color[2] = b; /* Blue */
};
return 0;
}
void setup() {
PixelWidth = (CxMax - CxMin) / iXmax;
PixelHeight = (CyMax - CyMin) / iYmax;
/*create new file,give it a name and open it in binary mode */
fp = fopen(filename, "wb"); /* b - binary mode */
/*write ASCII header to the file*/
fprintf(fp, "P6\n %s\n %d\n %d\n %d\n", comment, iXmax, iYmax, MaxColorComponentValue);
}
void info() {
double distortion;
// widt/height
double PixelsAspectRatio = (double) iXmax / iYmax; // https://en.wikipedia.org/wiki/Aspect_ratio_(image)
double WorldAspectRatio = (CxMax - CxMin) / (CyMax - CyMin);
printf("PixelsAspectRatio = %.16f \n", PixelsAspectRatio);
printf("WorldAspectRatio = %.16f \n", WorldAspectRatio);
distortion = PixelsAspectRatio - WorldAspectRatio;
printf("distortion = %.16f ( it should be zero !)\n", distortion);
// file
printf("file %s saved. It is called M-bump-1.png in A Cheritat wiki\n", filename);
}
void close() {
fclose(fp);
info();
}
// ************************************* main *************************
int main() {
complex double c;
setup();
printf(" render = compute and write image data bytes to the file \n");
for (iY = 0; iY < iYmax; iY++)
for (iX = 0; iX < iXmax; iX++) { // compute pixel coordinate
c = give_c(iX, iY);
/* compute pixel color (24 bit = 3 bytes) */
compute_color(c, color);
/*write color to the file*/
fwrite(color, 1, 3, fp);
}
close();
return 0;
}
However, it relies on derivative of the original fractal and is calculated in real time. In order to use this algorithm in post processing, I don't know how to transform the outputs from superheel's slope calculations to be the equivalent to the "reflection" calculations the above code.
I have tried to merge the two approaches without success.
Any suggestions? Thanks.