There is no complex-analytic derivative for zed-bar. But at least it doesn't have a "corner" like absolute value does, which gets replicated everywhere when things are iterated until it's nowhere smooth. It should be possible to do something with the Jacobian matrix, treating it as a system of two real variables.

Hi pauldelbrot
Thanks for replying. Please elaborate on "Jacobian matrix, treating it as a system of two real variables". This is way beyond my understanding at this time.
Many thanks.

A real analysis textbook (perhaps there's one at Wikibooks) could probably explain it better than I can.

I have been able to successfully create a Mandelbar (tricorn) fractal and implement slope using forward differencing:



However, when I try to create a slope using derivatives, I don't get the correct image



Here is the code:

Code: [Select]

    case 13: //  z1 = conj(z)); z = z1*z1+c; [Tricorn or Mandelbar]
Complex temp;
temp = Z;
temp.y = -temp.y; // get conjugate
dC = dC * (temp + temp) + 1.0;
sqr.x = temp.x * temp.x;
sqr.y = temp.y * temp.y;
real_imag = temp.x * temp.y;
Z.x = Q.x + sqr.x - sqr.y;
Z.y = Q.y + real_imag + real_imag;
break;

Here is the code for standard Mandelbrot and it works:

Code: [Select]

    case 0: //  z = z*z+c; Mandelbrot
dC = dC * (Z + Z) + 1.0;
sqr.x = Z.x * Z.x;
sqr.y = Z.y * Z.y;
real_imag = Z.x * Z.y;
Z.x = Q.x + sqr.x - sqr.y;
Z.y = Q.y + real_imag + real_imag;
break;

Can anyone see what I'm doing wrong in calculating the derivative? Is there some special trick needed to get the derivative of a complex conjugate?

Any help greatly appreciated. Thanks.

At least the second image generates some very nice coloring patterns!

Thanks FractalAlex,

My mistakes always generate "interesting" images, especially when I play with palettes 

Thanks Claude,

There's a lot to take in here. My immediate task is to learn how to use the Jacobian Matrix. Then I need to apply it to derivatives and finally to incorporate it into complex conjugates. I will spend some time looking at it.

I find it interesting that you can apply it to Burning Ship as it contains abs() and these are discontinuities in the curves and therefore there is no continuous derivative. I'm not sure how distance estimation relates to slope calculations, but then I haven't studied your links all that carefully as yet. I'll see what I can come up with.
Many thanks.

Hi Adam,
Please clarify what you need. Is it the Mandelbrot code for Forward differencing?

Maybe this code is what you're looking for...

/*
    FWDDIFF.CPP a module for determining the slope using forward differencing calculations of fractals.
   
    Written in Microsoft Visual 'C++' by Paul de Leeuw.

     https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/core/ThreadDraw.java#L3640
     https://github.com/hrkalona/Fractal-Zoomer/blob/master/src/fractalzoomer/main/app_settings/BumpMapSettings.java
*/

#include "slope.h"

/**************************************************************************
    Calculate functions
**************************************************************************/

void CSlope::DoSlopeFwdDiffFn(Complex *z, Complex *q)
    {
    Complex sqr, temp;
    Complex t = { 1.0, 0.0 };
    double real_imag;
    int k;
    int degree = (int)param[5];
    int degree1, degree2;

//    PaletteStart = (int)param[4];

    switch (subtype)
{
case 0: // Mandelbrot
    sqr.x = z->x * z->x;
    sqr.y = z->y * z->y;
    real_imag = z->x * z->y;
    z->x = q->x + sqr.x - sqr.y;
    z->y = q->y + real_imag + real_imag;
    break;
// other cases removed
}
    }

/**************************************************************************
    Get the gradients in the x and y directions
**************************************************************************/

double CSlope::getGradientX(double *wpixels, int index, int width)
    {
    int x = index % width;
    double it = *(wpixels + index);

    if (x == 0) {
return (*(wpixels + index + 1) - it) * 2;
}
    else if (x == width - 1) {
return (it - *(wpixels + index - 1)) * 2;
}
    else {
double diffL = it - *(wpixels + index - 1);
double diffR = it - *(wpixels + index + 1);
return diffL * diffR >= 0 ? 0 : diffL - diffR;
}
    }

double CSlope::getGradientY(double *wpixels, int index, int width, int height)
    {
    int y = index / width;
    double it = *(wpixels + index);

    if (y == 0) {
return (it - *(wpixels + index + width)) * 2;
}
    else if (y == height - 1) {
return (*(wpixels + index - width) - it) * 2;
}
    else {
double diffU = it - *(wpixels + index - width);
double diffD = it - *(wpixels + index + width);
return diffD * diffU >= 0 ? 0 : diffD - diffU;
}
    }

/**************************************************************************
   Brightness Scaling
**************************************************************************/

int CSlope::changeBrightnessOfColorScaling(int rgb, double delta)
    {
    int     new_color = 0;
    double bump_transfer_factor = param[0]; // future add to dialogue box

//    double mul = getBumpCoef(delta);
    double  mul = (1.5 / (fabs(delta * bump_transfer_factor) + 1.5));

    if (delta > 0) {
rgb ^= 0xFFFFFF;
int r = rgb & 0xFF0000;
int g = rgb & 0x00FF00;
int b = rgb & 0x0000FF;
int ret = (int)(r * mul + 0.5) & 0xFF0000 | (int)(g * mul + 0.5) & 0x00FF00 | (int)(b * mul + 0.5);
new_color = 0xff000000 | (ret ^ 0xFFFFFF);
}
    else {
int r = rgb & 0xFF0000;
int g = rgb & 0x00FF00;
int b = rgb & 0x0000FF;
new_color = 0xff000000 | (int)(r * mul + 0.5) & 0xFF0000 | (int)(g * mul + 0.5) & 0x00FF00 | (int)(b * mul + 0.5);
}

    return new_color;
    }

/**************************************************************************
Slope Fractal
**************************************************************************/

int CSlope::RunSlopeFwdDiff(HWND hwndIn, void(*plot)(WORD, WORD, DWORD), int user_data(HWND hwnd), char* StatusBarInfo, bool *ThreadComplete, int subtypeIn, int NumThreadsIn, int threadIn, Complex j, double mandel_width, double hor, double vert, double xgap, double ygap,
    double rqlim, long threshold, double paramIn[], CTrueCol *TrueCol, CDib *Dib, double *SlopeDataPtr, BYTE juliaflag)
    {
    Complex z = 0.0; // initial value for iteration Z0
    Complex c, q;
    BigComplex cBig, zBig, qBig;
    BigDouble BigTemp, BigTemp1;

    double log_zn, nu;
    int lastChecked = -1;

    int x, y, i;
    DWORD index;
    double iterations;

    thread = threadIn;
    NumThreads = NumThreadsIn;
    subtype = subtypeIn;
    hwnd = hwndIn;
    for (i = 0; i < NUMSLOPEDERIVPARAM; i++)
param[i] = paramIn[i];

    if (NumThreads == 0)
StripWidth = xdots;
    else
StripWidth = xdots / NumThreads;
    StripStart = StripWidth * thread;
    *ThreadComplete = false;

    for (y = 0; y < ydots; y++)
{
if (BigNumFlag)
    cBig.y = BigVert + BigWidth - Big_ygap * y;
else
    c.y = vert + mandel_width - y * ygap;
double progress = (double)y / ydots;
if (int(progress * 100) != lastChecked)
    {
    lastChecked = int(progress * 10);
    sprintf(StatusBarInfo, "Progess (%d%%), %d Threads", int(progress * 100), NumThreads);
    }
for (x = StripStart; x < StripStart + StripWidth; x++)
    {
    if (user_data(hwnd) < 0) // trap user input
return -1;

    c.x = hor + x * xgap;
       {
        if (juliaflag)
            {
            q = j;
            z = c;
            }
        else
            {
            q = c;
            z = 0.0;
            }
        }
    iterations = 0.0;
    for (;;)
{
iterations++;
if (iterations >= threshold)
    break;
DoSlopeFwdDiffFn(&z, &q);
if ((z.x * z.x + z.y * z.y) >= rqlim)
    break;
}
    if (iterations < threshold)
{
if (BigNumFlag)
    {
    // sqrt of inner term removed using log simplification rules.
    BigTemp = zBig.x * zBig.x + zBig.y * zBig.y;
    BigTemp1 = BigTemp.BigLog();
    log_zn = mpfr_get_d(BigTemp1.x, MPFR_RNDN) / 2;
    nu = log(log_zn / log(2.0)) / log(2.0);
    }
else
    {
    // sqrt of inner term removed using log simplification rules.
    log_zn = log(z.x * z.x + z.y * z.y) / 2;
    nu = log(log_zn / log(2.0)) / log(2.0);
    }
iterations = iterations + 1 - nu;
}
    index = ((DWORD)y * (DWORD)width) + (DWORD)x;
    if (x >= 0 && x < xdots - 1 && y >= 0 && y < ydots - 1)
*(SlopeDataPtr + index) = iterations;
    plot(x, y, (long)iterations);
    }
}
    *ThreadComplete = true;
    return 0;
    }

/**************************************************************************
Slope Fractal
**************************************************************************/

int CSlope::RenderSlope(long threshold, double paramIn[], CTrueCol *TrueCol, CDib *Dib, double *SlopeDataPtr)

    {
    double dotp, gradAbs, gradCorr, cosAngle, sizeCorr, smoothGrad, lightAngleRadians, lightx, lighty;
   
    for (int i = 0; i < NUMPARAM; i++)
param[i] = paramIn[i];

    double bumpMappingDepth = param[3]; // future add to dialogue box
    double bumpMappingStrength = param[4]; // future add to dialogue box
    double lightDirectionDegrees = param[2]; // future add to dialogue box
    int PaletteStart = (int)param[1];

    double iterations;
    int lastChecked = -1;
    DWORD index;
    int x, y;
    double gradx, grady;
    unsigned char r, g, b;
    RGBTRIPLE colour;
    int modified;

    lastChecked = -1;

    sizeCorr = 0.0;
    lightx = 0.0;
    lighty = 0.0;

    gradCorr = pow(2, (bumpMappingStrength - DEFAULT_BUMP_MAPPING_STRENGTH) * 0.05);
    sizeCorr = ydots / pow(2, (MAX_BUMP_MAPPING_DEPTH - bumpMappingDepth) * 0.16);
    lightAngleRadians = lightDirectionDegrees * PI / 180.0;
    lightx = cos(lightAngleRadians) * gradCorr;
    lighty = sin(lightAngleRadians) * gradCorr;

    for (y = 0; y < ydots; y++)
{
for (x = 0; x < xdots; x++)
    {
    index = ((DWORD)y * (DWORD)width) + (DWORD)x;
    if (SlopeDataPtr)
iterations = *(SlopeDataPtr + index);
    else
return 0; // oops, we don't actually have forward differencing

    if (iterations >= threshold)
{ //  interior of Mandelbrot set = inside_color
colour.rgbtRed = (BYTE)TrueCol->InsideRed; // M_waves
colour.rgbtGreen = (BYTE)TrueCol->InsideGreen;
colour.rgbtBlue = (BYTE)TrueCol->InsideBlue;
}
    else
{
// modified = rgbs[index];
if (iterations < PaletteStart)
    modified = 0x00FFFFFF;
else
    modified = 0xFF000000 | ((DWORD)*(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 0) << 16) | ((DWORD)*(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 1) << 8) | *(TrueCol->PalettePtr + (BYTE)(((long)iterations) % 256) * 3 + 2);
gradx = getGradientX(SlopeDataPtr, index, xdots);
grady = getGradientY(SlopeDataPtr, index, xdots, ydots);
dotp = gradx * lightx + grady * lighty;
// int original_color = modified; // not sure what this is for
if (dotp != 0)
    {
    gradAbs = sqrt(gradx * gradx + grady * grady);
    cosAngle = dotp / gradAbs;
    smoothGrad = -2.3562 / (gradAbs * sizeCorr + 1.5) + 1.57;
    //smoothGrad = Math.atan(gradAbs * sizeCorr);
    modified = changeBrightnessOfColorScaling(modified, cosAngle * smoothGrad);
    }
r = (modified >> 16) & 0xFF;
g = (modified >> 8) & 0xFF;
b = modified & 0xFF;
colour.rgbtRed = r;
colour.rgbtGreen = g;
colour.rgbtBlue = b;
}
    outRGBpoint(Dib, x, y, colour);
    }
}
    return 0;
    }


Let me know if you need more info.

slope relates to directional DE because both are based on spatial derivative of escape time
https://code.mathr.co.uk/kalles-fraktaler-2/blob/ec44abb522811ea93e03c8a27933b19472d54bc1:/fraktal_sft/fraktal_sft.cpp#l1053 is the implementation in KF with analytic DE (simple) vs forward differences (ugly, inaccurate, noisy when jitter is enabled)

re jacobians,
https://code.mathr.co.uk/kalles-fraktaler-2/blob/ec44abb522811ea93e03c8a27933b19472d54bc1:/formula/formula.xml#l175 is the implementation in KF including derivatives for Burning Ship
https://code.mathr.co.uk/kalles-fraktaler-2/blob/ec44abb522811ea93e03c8a27933b19472d54bc1:/formula/formula.xml#l402 is Mandelbar
Hopefully those examples are clear enough.  (daa, dab; dba, dba) is the identity matrix for normal views, zooming/rotating/skewing changes these.