Author Topic: Domain Coloring (Read 11141 times)

_vince · « **on:** August 15, 2020, 02:13:58 pm »

A kind of method for plotting complex functions:
https://en.wikipedia.org/wiki/Domain_coloring

So pretty I had to share this

Code: [Select]

defdbl a-z
dim shared pi as double
pi = 4*atn(1)

const sw = 800
const sh = 600

declare function hrgb(h as double, m as double)

screen _newimage(sw, sh, 32)


'dim as double x, y, u, v, m, a, uu, vv, x0, y0

dim p as string

for j=0 to 7
for yy=0 to sh
for xx=0 to sw
        select case j
        case 0
                u = (xx - sw/2)*0.015
                v = (sh/2 - yy)*0.015

                x = u
                y = v

                p = "z"
        case 1
                u = (xx - sw/2)*0.01 + 2
                v = (sh/2 - yy)*0.01

                x = exp(u)*cos(v)
                y = exp(u)*sin(v)

                p = "exp(z)"
        case 2
                u = (xx - sw/2)*0.015
                v = (sh/2 - yy)*0.015

                x = sin(u)*_cosh(v)
                y = cos(u)*_sinh(v)

                p = "sin(z)"
        case 3
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 0
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu 
                        v = vv 
                next
                x = u
                y = v

                p = "f(z) = z^2 + c"
        case 4
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 1
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu 
                        v = vv 
                next
                x = u
                y = v

                p = "f(f(z))"
        case 5
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 3
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu 
                        v = vv 
                next
                x = u
                y = v

                p = "f(f(f(f(z))))"
        case 6
                u = (xx - sw/2)*0.005 - 0.5
                v = (sh/2 - yy)*0.005
                x0 = u
                y0 = v
                for i=0 to 9
                        uu = u*u - v*v + x0
                        vv = 2*u*v + y0

                        u = uu 
                        v = vv 
                next
                x = u
                y = v

                p = "f(f(f(f(f(f(f(f(f(f(z))))))))))"
	case 7
		u = (xx - sw/2)*0.0015
		v = (sh/2 - yy)*0.0015

		x = sin(u/(u*u + v*v))*_cosh(-v/(u*u + v*v))
		y = cos(u/(u*u + v*v))*_sinh(-v/(u*u + v*v))

		p = "sin(1/z)"
        end select

        m = sqr(x*x + y*y)
        a = (pi + _atan2(y, x))/(2*pi)

        pset (xx, yy), hrgb(a, m)

        locate 1,1: ? p
next
next

sleep
next

system

function hrgb(h as double, m as double)
	'dim as double r, g, b, mm

	r = 0.5 - 0.5*sin(2*pi*h - pi/2)
    g = (0.5 + 0.5*sin(2*pi*h*1.5 - pi/2)) * -(h < 0.66)
    b = (0.5 + 0.5*sin(2*pi*h*1.5 + pi/2)) * -(h > 0.33)

	mm = (m*100 mod 100)*0.01

	if ((m*2*pi*100 mod 2*pi*100) < 50) or ((h*2*pi*100 mod 4*2*pi) < 6) then
		'hrgb = _rgb(0,0,0)
		hrgb = _rgb(15*r, 155*g, 155*b)
	else
		hrgb = _rgb(255*r, 255*g, 255*b)
	end if
end function

_vince · « **Reply #1 on:** August 15, 2020, 03:23:42 pm »

Pretty enough to take pictures

bplus · « **Reply #2 on:** August 15, 2020, 03:33:36 pm »

@_vince Nice! How did you come across this?

_vince · « **Reply #3 on:** August 15, 2020, 03:42:55 pm »

I found the greatest book ever https://www.springer.com/gp/book/9783034801799

Lots of pretty pictures, I can email you a pdf copy if you want

euklides · « **Reply #4 on:** August 16, 2020, 04:18:04 am »

Nice to see !
And you can change the cases, like...

= (xx - sw / 2) * 0.035 + SIN(RND * 6)
v = (sh / 2 - yy) * 0.015
x = SIN(u + 1) * _COSH(v)
y = COS(u - 1) * _SINH(v)

or

u = (xx - sw / 3) * 0.035 + SIN(RND * 6)
v = (sh / 2 - yy + u * 15) * 0.025
x = SIN(u + 1) * 1 / _COSH(v)
y= 1 / (COS(u - 1) * _SINH(v)) + 2

and so on...

Ashish · « **Reply #5 on:** August 16, 2020, 10:24:26 am »

Beautilful! I like it.

bplus · « **Reply #6 on:** August 16, 2020, 11:59:04 am »

:) Almost has Mandelbrot looking like exotic butterfly.

and sounds just like "f(f(f(f(f(f(f(f(f(f(z))))))))))"

_vince · « **Reply #7 on:** August 17, 2020, 03:01:42 pm »

Some more picture, Taylor expansion of sine

You can clearly see all the complex zeros and the number of them is the order of the polynomial, not something you'd see in a real plot

_vince · « **Reply #8 on:** August 17, 2020, 03:12:50 pm »

Moving average filter

Again, a fun image of how zeros and poles look, you see how magnitude contours get closer and closer together around poles as the magnitude quickly increases to infinity

bplus · « **Reply #9 on:** August 17, 2020, 03:50:11 pm »

Hi @_vince

Looking over code today. Why is red discriminated? Could this be a typo? There seems to be no lack of color?

From the hrgb sub in original post:

Code: QB64: [Select]

        if ((m*2*pi*100 mod 2*pi*100) < 50) or ((h*2*pi*100 mod 4*2*pi) < 6) then
                'hrgb = _rgb(0,0,0)
                hrgb = _rgb(15*r, 155*g, 155*b)  '  <<< what's with red, typo ?
        else
                hrgb = _rgb(255*r, 255*g, 255*b)
        end if
 

_vince · « **Reply #10 on:** August 17, 2020, 03:54:50 pm »

Yes, just a typo. That's why the contours don't have the red tint around red.

This colouring scheme isn't particularly good, At some point I'm intending to make a really nice one like in pictures like this, https://en.wikipedia.org/wiki/File:Domain_coloring_x2-1_x-2-i_x-2-i_d_x2%2B2%2B2i.xcf

bplus · « **Reply #11 on:** August 17, 2020, 04:10:29 pm »

I discovered an effect while messing with the colors, it's a moving experience ;-))

Code: QB64: [Select]

DEFDBL A-Z
DIM SHARED pi AS DOUBLE
pi = 4 * ATN(1)
 
CONST sw = 800
CONST sh = 600
_TITLE "b+ fiddles... press spacebar "
declare function hrgb(h as double, m as double)
 
SCREEN _NEWIMAGE(sw, sh, 32)
 
 
'dim as double x, y, u, v, m, a, uu, vv, x0, y0
 
DIM p AS STRING
 
FOR j = 6 TO 6
    FOR yy = 0 TO sh
        FOR xx = 0 TO sw
            SELECT CASE j
                CASE 0
                    u = (xx - sw / 2) * 0.015
                    v = (sh / 2 - yy) * 0.015
 
                    x = u
                    y = v
 
                    p = "z"
                CASE 1
                    u = (xx - sw / 2) * 0.01 + 2
                    v = (sh / 2 - yy) * 0.01
 
                    x = EXP(u) * COS(v)
                    y = EXP(u) * SIN(v)
 
                    p = "exp(z)"
                CASE 2
                    u = (xx - sw / 2) * 0.015
                    v = (sh / 2 - yy) * 0.015
 
                    x = SIN(u) * _COSH(v)
                    y = COS(u) * _SINH(v)
 
                    p = "sin(z)"
                CASE 3
                    u = (xx - sw / 2) * 0.005 - 0.5
                    v = (sh / 2 - yy) * 0.005
                    x0 = u
                    y0 = v
                    FOR i = 0 TO 0
                        uu = u * u - v * v + x0
                        vv = 2 * u * v + y0
 
                        u = uu
                        v = vv
                    NEXT
                    x = u
                    y = v
 
                    p = "f(z) = z^2 + c"
                CASE 4
                    u = (xx - sw / 2) * 0.005 - 0.5
                    v = (sh / 2 - yy) * 0.005
                    x0 = u
                    y0 = v
                    FOR i = 0 TO 1
                        uu = u * u - v * v + x0
                        vv = 2 * u * v + y0
 
                        u = uu
                        v = vv
                    NEXT
                    x = u
                    y = v
 
                    p = "f(f(z))"
                CASE 5
                    u = (xx - sw / 2) * 0.005 - 0.5
                    v = (sh / 2 - yy) * 0.005
                    x0 = u
                    y0 = v
                    FOR i = 0 TO 3
                        uu = u * u - v * v + x0
                        vv = 2 * u * v + y0
 
                        u = uu
                        v = vv
                    NEXT
                    x = u
                    y = v
 
                    p = "f(f(f(f(z))))"
                CASE 6
                    u = (xx - sw / 2) * 0.005 - 0.5
                    v = (sh / 2 - yy) * 0.005
                    x0 = u
                    y0 = v
                    FOR i = 0 TO 9
                        uu = u * u - v * v + x0
                        vv = 2 * u * v + y0
 
                        u = uu
                        v = vv
                    NEXT
                    x = u
                    y = v
 
                    p = "f(f(f(f(f(f(f(f(f(f(z))))))))))"
                CASE 7
                    u = (xx - sw / 2) * 0.0015
                    v = (sh / 2 - yy) * 0.0015
 
                    x = SIN(u / (u * u + v * v)) * _COSH(-v / (u * u + v * v))
                    y = COS(u / (u * u + v * v)) * _SINH(-v / (u * u + v * v))
 
                    p = "sin(1/z)"
            END SELECT
 
            m = SQR(x * x + y * y)
            a = (pi + _ATAN2(y, x)) / (2 * pi)
 
            PSET (xx, yy), hrgb(a, m)
 
            LOCATE 1, 1: PRINT p
        NEXT
    NEXT
    DO
        _SCREENMOVE 200 + RND * 4 - 2, 80 + RND * 3 - 1.5
        _LIMIT 30
    LOOP UNTIL _KEYDOWN(32)
    SLEEP
NEXT
 
SYSTEM
 
FUNCTION hrgb (h AS DOUBLE, m AS DOUBLE)
    'dim as double r, g, b, mm
 
    r = 0.5 - 0.5 * SIN(2 * pi * h - pi / 2)
    g = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 - pi / 2)) * -(h < 0.66)
    b = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 + pi / 2)) * -(h > 0.33)
 
    mm = (m * 100 MOD 100) * 0.01
 
    IF ((m * 2 * pi * 100 MOD 2 * pi * 100) < 50) OR ((h * 2 * pi * 100 MOD 4 * 2 * pi) < 6) THEN
        'hrgb = _rgb(0,0,0)
        hrgb = _RGB(115 * r, 115 * g, 115 * b)
    ELSE
        hrgb = _RGB(255 * r, 255 * g, 255 * b)
    END IF
END FUNCTION
 
 

_vince · « **Reply #12 on:** August 17, 2020, 07:24:11 pm »

edit:

_vince · « **Reply #13 on:** August 19, 2020, 12:54:08 pm »

https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula

Cauchy's integral formula is a powerful result in complex analysis that says you can evaluate a function f(z0) inside a contour just by taking the contour integral of f(z)/(z - z0). It essentially says that the function inside the contour can be defined entirely by the contour.

To begin, in the first image I plot the function

f(z) = z - z^3/3! + z^5/5! - z^7/7!

and draw the contour:

x(t) = pi*cos(t) + 3*pi/2
y(t) = pi*sin(t)

Then I numerically evaluate the following line integral using trapezoidal rule: https://en.wikipedia.org/wiki/Trapezoidal_rule

f(z0) = (1/(2*pi*i)) \cint_C f(z)/(z - z0) dz

Firstly for N = 30 and you see all the distortion around the contour and then for N = 100 and N = 1000 which is dead on. And you see that it converges to the original function as you'd expect. The pretty artifacts around the contour are the errors due to numerical integration, it is like a complex visual of https://upload.wikimedia.org/wikipedia/commons/d/d1/Integration_num_trapezes_notation.svg

Code: QB64: [Select]

'#lang "fblite"
DEFLNG A-Z
 
CONST sw = 800
CONST sh = 600
 
DIM SHARED pi AS DOUBLE
pi = 4 * ATN(1)
 
declare function hrgb(h as double, m as double)
'declare function cosh(x as double) as double
'declare function sinh(x as double) as double
declare sub cexp(u as double, v as double, x as double, y as double, a as double, b as double)
declare sub cmul(u as double, v as double, x as double, y as double, a as double, b as double)
declare sub cdiv(u as double, v as double, x as double, y as double, a as double, b as double)
 
'screenres sw, sh, 32
SCREEN _NEWIMAGE(sw, sh, 32)
 
DIM m AS DOUBLE, a AS DOUBLE
DIM x AS DOUBLE, y AS DOUBLE, u AS DOUBLE, v AS DOUBLE, uu AS DOUBLE, vv AS DOUBLE
 
n = 30
DIM sum AS DOUBLE, delta AS DOUBLE, p AS DOUBLE, q AS DOUBLE, pp AS DOUBLE, qq AS DOUBLE
 
 
FOR i = 0 TO 2
 
    FOR yy = 0 TO sh
        FOR xx = 0 TO sw
    
            SELECT CASE i
                CASE 0
                    u = (xx - sw / 2) * 0.015
                    v = (sh / 2 - yy) * 0.015
                    x = u
                    y = v
                CASE 1
                    u = (xx - sw / 2) * 0.01 + 3 * pi / 2
                    v = (sh / 2 - yy) * 0.01
 
                    x = u
                    y = v
        
                    cexp uu, vv, u, v, 3, 0
                    x = x - uu / 6
                    y = y - vv / 6
        
                    cexp uu, vv, u, v, 5, 0
                    x = x + uu / 120
                    y = y + vv / 120
        
                    cexp uu, vv, u, v, 7, 0
                    x = x - uu / 5040
                    y = y - vv / 5040
                CASE 2
                    u = (xx - sw / 2) * 0.01 + 3 * pi / 2
                    v = (sh / 2 - yy) * 0.01
 
                    x = u
                    y = v
        
                    cexp uu, vv, u, v, 3, 0
                    x = x - uu / 6
                    y = y - vv / 6
        
                    cexp uu, vv, u, v, 5, 0
                    x = x + uu / 120
                    y = y + vv / 120
 
                    cexp uu, vv, u, v, 7, 0
                    x = x - uu / 5040
                    y = y - vv / 5040
        
                    sum = 0
                    FOR j = 0 TO n - 1
                        p = pi * COS(j * 2 * pi / n) + 3 * pi / 2
                        q = pi * SIN(j * 2 * pi / n)
            
                        pp = p
                        qq = q
                        cexp uu, vv, p, q, 3, 0
                        pp = pp - uu / 6
                        qq = qq - vv / 6
                        cexp uu, vv, p, q, 5, 0
                        pp = pp + uu / 120
                        qq = qq + vv / 120
                        cexp uu, vv, p, q, 7, 0
                        pp = pp - uu / 5040
                        qq = qq - vv / 5040
            
                        cdiv uu, vv, pp, qq, p - u, q - v
            
                        IF j = 0 OR j = n - 1 THEN
                            sum = sum + 0.5 * (-uu * pi * SIN(j * 2 * pi / n) - vv * pi * COS(j * 2 * pi / n))
                        ELSE
                            sum = sum + (-uu * pi * SIN(j * 2 * pi / n) - vv * pi * COS(j * 2 * pi / n))
                        END IF
                    NEXT
                    x = sum * 2 * pi / n
        
                    sum = 0
                    FOR j = 0 TO n - 1
                        p = pi * COS(j * 2 * pi / n) + 3 * pi / 2
                        q = pi * SIN(j * 2 * pi / n)
            
                        pp = p
                        qq = q
                        cexp uu, vv, p, q, 3, 0
                        pp = pp - uu / 6
                        qq = qq - vv / 6
                        cexp uu, vv, p, q, 5, 0
                        pp = pp + uu / 120
                        qq = qq + vv / 120
                        cexp uu, vv, p, q, 7, 0
                        pp = pp - uu / 5040
                        qq = qq - vv / 5040
            
                        cdiv uu, vv, pp, qq, p - u, q - v
            
                        IF j = 0 OR j = n - 1 THEN
                            sum = sum + 0.5 * (-vv * pi * SIN(j * 2 * pi / n) + uu * pi * COS(j * 2 * pi / n))
                        ELSE
                            sum = sum + (-vv * pi * SIN(j * 2 * pi / n) + uu * pi * COS(j * 2 * pi / n))
                        END IF
                    NEXT
                    y = sum * 2 * pi / n
        
                    cmul uu, vv, x, y, 0, -1 / (2 * pi)
        
                    x = uu
                    y = vv
 
 
            END SELECT
 
            m = SQR(x * x + y * y)
            a = (pi + _ATAN2(y, x)) / (2 * pi)
 
            PSET (xx, yy), hrgb(a, m)
        NEXT
    NEXT
 
    CIRCLE (sw / 2, sh / 2), pi / 0.01 - 2, _RGB(0, 0, 0)
    CIRCLE (sw / 2, sh / 2), pi / 0.01 - 1, _RGB(0, 0, 0)
    CIRCLE (sw / 2, sh / 2), pi / 0.01, _RGB(255, 255, 0)
    CIRCLE (sw / 2, sh / 2), pi / 0.01 + 1, _RGB(0, 0, 0)
    CIRCLE (sw / 2, sh / 2), pi / 0.01 + 2, _RGB(0, 0, 0)
 
    'for a = 0 to 2*pi step 0.01
    '    x = pi*cos(a) + 3*pi/2
    '    y = pi*sin(a)
 
    '    pset ((x/0.015 + sw/2), (sh/2 - y/0.015)), rgb(255,255,255)
    'next
 
    SLEEP
NEXT
 
SYSTEM
 
FUNCTION hrgb (h AS DOUBLE, m AS DOUBLE)
    DIM r AS DOUBLE, g AS DOUBLE, b AS DOUBLE, t AS DOUBLE
 
    r = 0.5 - 0.5 * SIN(2 * pi * h - pi / 2)
    g = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 - pi / 2)) * -(h < 0.66)
    b = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 + pi / 2)) * -(h > 0.33)
 
    n = 16
    't = 0.3/m
    t = 0.2
    h = h + (0.5 * t / n)
 
    IF ((m * 500 MOD 500) < 70) OR (ABS((h * n) - INT(h * n)) < t) THEN
        'hrgb = rgb(0,0,0)
        hrgb = _RGB(125 * r, 125 * g, 125 * b)
    ELSE
        hrgb = _RGB(255 * r, 255 * g, 255 * b)
    END IF
END FUNCTION
 
SUB cexp (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
    DIM lnz AS DOUBLE, argz AS DOUBLE, mag AS DOUBLE, ang AS DOUBLE
 
    lnz = x * x + y * y
 
    IF lnz = 0 THEN
        u = 0
        v = 0
    ELSE
        lnz = 0.5 * LOG(lnz)
        argz = _ATAN2(y, x)
 
        mag = EXP(a * lnz - b * argz)
        ang = a * argz + b * lnz
 
        u = mag * COS(ang)
        v = mag * SIN(ang)
    END IF
END SUB
 
SUB cmul (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
    u = x * a - y * b
    v = x * b + y * a
END SUB
 
SUB cdiv (u AS DOUBLE, v AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE, b AS DOUBLE)
    DIM d AS DOUBLE
    d = a * a + b * b
    u = (x * a + y * b) / d
    v = (y * a - x * b) / d
END SUB
 

bplus · « **Reply #14 on:** August 19, 2020, 01:10:47 pm »

Hey _vince your graphs remind me of my study of Inversion Points of 2D Circle:
https://www.qb64.org/forum/index.php?topic=2933.msg121929#msg121929

News:

Author Topic: Domain Coloring (Read 11141 times)

_vince

Domain Coloring

_vince

Re: Domain Coloring

bplus

Re: Domain Coloring

_vince

Re: Domain Coloring

euklides

Re: Domain Coloring

Ashish

Re: Domain Coloring

bplus

Re: Domain Coloring

_vince

Re: Domain Coloring

_vince

Re: Domain Coloring

bplus

Re: Domain Coloring

_vince

Re: Domain Coloring

bplus

Re: Domain Coloring

_vince

Re: Domain Coloring

_vince

Re: Domain Coloring

bplus

Re: Domain Coloring