Saturday, 10 March 2018

Grant me the serenity to accept the things I cannot change, courage to change the things I can, and wisdom to know the difference.

"Briefly stated, the Gell-Mann Amnesia effect is as follows. You open the newspaper to an article on some subject you know well. In Murray’s case, physics. In mine, show business. You read the article and see the journalist has absolutely no understanding of either the facts or the issues. Often, the article is so wrong it actually presents the story backward—reversing cause and effect. I call these the 'wet streets cause rain' stories. Paper’s full of them. In any case, you read with exasperation or amusement the multiple errors in a story, and then turn the page to national or international affairs, and read as if the rest of the newspaper was somehow more accurate about Palestine than the baloney you just read. You turn the page, and forget what you know." - Michael Crichton

Single resistor VGA output from an FPGA. VGA pins 13 and 14, vertical and horizontal sync are connected directly to the FPGA (Altera Cyclone II), pins 1-3 are RGB and connected all parallel and through a 100 ohm resistor to the FPGA. Pins 6-8 are ground.

PLL is running at 25 MHz giving monochrome VGA output. Adding color would be relatively trivial.

magick convert -fuzz 40% -colors 8 -layers Optimize -delay 5 output.gif out.gif

module test(clock, out0i, out1i, out0q, out1q);

input clock;
output out0i;
output out1i;
output out0q;
output out1q;

reg vga_HS, vga_VS;
reg [9:0] CounterX;
reg [8:0] CounterY;
reg [9:0] temporary;
reg [9:0] cntx;
reg [9:0] cnty;
reg pixvalue;

reg [7:0] fnt[6399:0];
initial $readmemh("fnt.hex", fnt);

always @(posedge pll_clock)
CounterX <= CounterX + 1;
CounterX <= 0;
CounterY <= CounterY + 1;

CounterY <= 0;
temporary <= temporary + 1;
vga_HS <= (CounterX>640+16 & CounterX<800-48);
vga_VS <= (CounterY>400+12 & CounterY<449-35);
cntx <= CounterX + temporary;
cnty <= CounterY + (temporary >> 2);
pixvalue <= ((CounterX<640 & CounterY<400) & (CounterX[2] & CounterY[2])) | ((CounterX>128 & CounterX<512 & CounterY>64 & CounterY<300) & (cntx[4] & cnty[4]));
if(CounterX<640) pixvalue <= (fnt[CounterX+CounterY*760]>1);

assign out0i = ~vga_HS;
assign out1i = ~vga_VS;
assign out0q = pixvalue;

altpll0 altpll0 (

Output from a 1k resistor ladder DAC (3 resistors) with Cyclone II.
Driven by LVDS serializers at 500 MHz.


I think the concept of "culture" as used today should be clearly separated from the concept of religion. Religions are for the most part dogmas, quite different from culture which consists of temporary trends and behavioral models. Pretty much everyone knows the difference, but for some reason it is often ignored. We can judge the historical mistakes without holding anyone responsible for the sins of their fathers, but it is only natural and practical to judge the ideas and behaviors we now clearly see as inferior.

Simulation of quadrature signals from parametrically amplified quantum noise in a superconducting tunnel junction.

Stability through slavery.

Monday, 12 February 2018

LG OLED557V, an excellent panel, but a bit fucked up scaler

I recently purchased LG OLED55B7V. It's a great 3840x2160 panel with 12 bpc (36 bits per pixel) color depth, perfect blacks, HDR and 120 Hz refresh rate. However, its software (WebOS with Smart Share) for playing videos seems to be not performing as well as one would hope. It looks like the scaler has some problem with chrominance (when the video is not natively 4k). The same problem exists for all setting, and all videos regardless of their bitrate, but it's easiest to see with low resolution videos and high color saturation. This is an expensive OLED device so I would have expected the scalers to match the panel.

WebOS doesn't seem to be doing a good job at rendering videos. Pure colors tend to appear very blocky.
The flaw is a little bit odd because it doesn't quite correlate with the chrominance resolution, but is instead more blocky than one would expect from simply incorrect rendering of chrominance.
The clip is from a movie called The Black Hole (1979).
The same video rendered by VLC in Windows 10 (connected with HDMI).
No such annoying blocking is visible.

The source image is a single red line (2 pixels in diameter). Clearly the TV applies some kind of nasty filter on it. This effect is not visible when the video is natively 4k.
ffmpeg -r 24 -f image2 -i rb.png -filter_complex loop=399:1 -pix_fmt yuv420p -vcodec libx264 -crf 1 rb.mp4 
VLC in Windows 10 (connected with HDMI).
The original image used in these tests.
ffmpeg -r 24 -f image2 -i image.png -filter_complex loop=799:1 -an vcodec mpeg4 -qscale 1 video.mp4
It would also be very nice if there were connections which allowed one to actually plug into the panel directly in such a manner that one could use the full potential of the panel (4k/12bpcRGB/120Hz). Apparently the version of the HDMI standard used allows only 4k/8bpcRGB/60Hz, 4k/12bpcYUV420/60Hz or lesser resolutions with higher refresh rate. Never the less, it's the first time in a couple of decades I've been happy with the black levels. Luckily rendering problems can always be avoided by simply using a PC.


Nothing to do with my new TV, but an association developed to the stupid gamma error in most picture scaling software which is nothing much more than due to the fact that pixel brightness (power) is proportional to its value (voltage) squared and most scaling is done on pixel values rather than on pixel brightness. The result is this.

The left image is the original. Middle one is scaled by averaging pixels values.
The right one is scaled by correctly averaging pixel powers.

Sometimes the result of incorrect scaling is significantly distorted colors.
The left image is the original. Middle one scaled by averaging pixels values.
The right one is scaled by averaging pixel powers.

RGBA = double(imread('gamma_dalai_lama_gray.jpg')); figure; image(RGBA/255); truesize;
RGBB = double(imread('gamma_colors.jpg')); figure; image(RGBB/255); truesize;

RGB1 = imresize(RGBA, 0.5, 'bicubic')/255;
RGB2 = abs(sqrt(imresize(RGBA.^2, 0.5, 'bicubic')))/255;

RGB3 = imresize(RGBB, 0.5, 'bicubic')/255;
RGB4 = abs(sqrt(imresize(RGBB.^2, 0.5, 'bicubic')))/255;

figure; image(imresize(RGB1, 2, 'nearest')); truesize; imwrite(RGB1, '001.jpg');
figure; image(imresize(RGB2, 2, 'nearest')); truesize; imwrite(RGB2, '002.jpg');
figure; image(imresize(RGB3, 2, 'nearest')); truesize; imwrite(RGB3, '003.jpg');
figure; image(imresize(RGB4, 2, 'nearest')); truesize; imwrite(RGB4, '004.jpg');


While we're at it, let's talk about YUV420. This is a signal format used in video compression as human vision is much less sensitive to changes in color than it is to changes in brightness. Instead of recording the primary colors red, green and blue, it records luminance (which is overall brightness of the pixel) with full resolution and chrominance (which records color separate from luminance) with half the original resolution. So for example 4k material recorded in YUV420 as usual has a resolution of 3840x2160 for luminance and a resolution of 1920x1080 for chrominance.

As usual these images can be rendered correctly and incorrectly. Since most display devices are natively RGB these images are eventually always converted back to RGB, but this conversion can be done poorly. A typical mistake is that while rendering the RGB data the chrominance is not interpolated. This results in images where colors are blocky. This is especially apparent in situations where the luminance is the same for both colors, but the colors are different. Below there is an example of how this effect looks like.

Improperly rendered YUV420 image data. Luminance is rendered properly, but the chrominance data is without interpolation and makes the image look much blockier than with proper interpolation. This effect is visible only in regions with high color saturation and especially in regions where there is no difference in the luminance of the signal, but large differences in the chrominance.

The original RGB image vs. properly rendered YUV420 where the chrominance signal is interpolated to match the luminance resolution.

RGB = imread('chroma3.png');
YCBCR = rgb2ycbcr(RGB);
N = 2;
YCBCR1(:, :, 2) = imresize(imresize(YCBCR(:, :, 2), 1/N, 'bicubic'), N, 'nearest');
YCBCR1(:, :, 3) = imresize(imresize(YCBCR(:, :, 3), 1/N, 'bicubic'), N, 'nearest');
YCBCR2(:, :, 2) = imresize(imresize(YCBCR(:, :, 2), 1/N, 'bicubic'), N, 'bicubic');
YCBCR2(:, :, 3) = imresize(imresize(YCBCR(:, :, 3), 1/N, 'bicubic'), N, 'bicubic');
RGB2 = ycbcr2rgb(YCBCR1);
RGB3 = ycbcr2rgb(YCBCR2);

So let's have a quick look at 30 bit HDR material as well. The compression codec these videos typically use is H.265 also known as HEVC. Pixel format is YUV420p10le where le refers to little endian (signifying order of bit significance). While standard dynamic range is 256 different values given by 8 bpc (8 bits for each of the three colors), the high dynamic range is typically in these videos 1024 different values (10 bpc). The main advantage typically being that more detail can exist in the extremely dark and extremely bright parts of the image.

If you consider the middle image as the standard range of light levels your camera would see, then you can observe that it is missing details both extremely bright and extremely dark or in other words the image is simultaneously overexposed and underexposed. Having additional bits allows one to capture the details in these missing regions. The middle image shows middle 8 bits of the total of 10 and top and bottom show the brightest and darkest 8.

Suppose your standard dynamic range captures brightness values between the blue bars. Additional bits allow capture of the information outside of your standard range.

The increased number of bits allows more brightness values which can either be used to increase the dynamic range (more darks and brights) or decrease the minimum distance between different brightness values or some compromise between these two.

Top most picture "simulates" (by exaggerated banding and limited maximum brightness) standard dynamic range with standard number of bits. Limited number of bits results in banding of the brightness values in image which would ideally be continuous change from total darkness to maximum brightness. Middle picture "simulates" a situation with standard dynamic range with two extra bits used to decrease the minimum spacing between brightness values. The third picture "simulates" high dynamic range with two extra bits. One of the bits is used to decrease minimum spacing between brightness values and one to increase the dynamic range now allowing colors brighter than standard dynamic range. This type of high dynamic range still shows more banding than standard dynamic range with two more bits, but less banding than standard dynamic range with standard number of bits. These cases are analogous to SDR8, SDR10 and HDR10.

I'm personally quite capable of seeing the banding in regular SDR8 stripes. Though, I will admit that under normal circumstances when the image contains certain amount of noise, it results in dithering which pretty much masks any banding that might be otherwise visible.

Both of the images above have the same amount of colors.
The lower one is just "dithered" in such a way that spatial pixel noise reflects continuous gradient.

One could also dither in time by making the pixels flicker at suitable duty cycles.
Of course combining some noise improves the result significantly.

X = 480;
img = zeros(64, X);
for x = 1:64
    img(x, :) = linspace(0, 1, X);
%     img(x, :) = linspace(0.3, 0.4, X);

M = 32;
% M = 256;

im = round((M-1)*img);

pix = zeros(64, X);
for y = 1:64
    for x = 1:X
        p = (M-1)*img(y, x);
        q = round((M-1)*img(y, x));
        err = p-q;
        if rand-0.5>err
            pix(y, x) = q;
            pix(y, x) = q + 1;



Sunday, 7 January 2018

One of my greatest fears is that one day I'll wake up in a gridlock I predicted ages ago

Einstein field equations (EFE)

describe the four dimensional geometry g associated with the points of spacetime and corresponding stress-energy tensor T. The differential terms are a function of the metric as described below.

Variable g with upper indices (contravariant) is the inverse of the (covariant) metric.

T describes flux of four momentum between different dimensions, including time. The first term (with upper indices tt) is simply the invariant mass-energy. For a single isolated particle the following (in natural units) applies (t-direction component of 4-dimensional v is c, v with lower index 3 represent the regular 3-dimensional velocity vector).

Contravariant stress-energy tensor can be converted into covariant form by the following.

EFE are a kind of generalization of Newtons gravitational potential (Poisson equation).

EFE describe a kind of 4-dimensional array of matrices (in naive numerical sense) corresponding to the structure of such spacetime. This sort of block universe can be seen as static. Theoretically it's possible that in this sort of block there could exist "flows" and geometries that form a closed loop in what we would normally call time (allowing a kind of time travel into the past). The naive interpretation of EFE would suggest that if such loops exist, they ought to form consistent spacetime. That is to say in order for them to be a solution to the equations, no grandfather paradoxes can form. You will (for one reason or another) only do such things in the "past" that are consistent with the block as a whole. Quantum mechanics may significantly complicate the situation if something along the lines of these loops are possible, but instead lead to alternative realities. The block of all universes must still be consistent as a whole, but as there are now several alternative realities and bridges between them, the consistency requirement can still apply and allow some particular futures to be consistent with time traveller "causing" them.

How to actually compute a solution to some practical problem with EFE is a topic for some other time.
The end justifies the means as long as the means are factored into the end.

If you could actually rationalize your opinion, it wouldn't be an opinion.

I don't really care about subjective changes all that much, it's the objective ones that I find interesting.

Every January I remember my mortality, but then I forget it for a while again.

Sometimes I consider will a flaw.

What to do when you reach a conclusion that is highly likely to be sound, yet highly undesirable?'s just not in my nature.

A true explanation must come to an end where no further explanations are required and none can exist, but can it be?
BTW. The Great Philosophers by Stephen Law is pretty nice and compact book about philosophy.
double g[4][4][N*N*N*N]; /* spacetime metric tensors */
double T[4][4][N*N*N*N]; /* spacetime stress-energy tensors */

int npow(int i) {
  int val = 1;
  for(int j=0; j<i; j++)
    val = val*N;
  return val;

/* computes values of individual Christoffel symbols */
double christoffel(int i, int j, int l, int loc) {
  double val = 0.0, d;
  for(int k=0; k<4; k++) {
    d = \
      (g[k][i][loc + npow(j)] - g[k][i][loc - npow(j)]) +  \
      (g[k][j][loc + npow(i)] - g[k][j][loc - npow(i)]) -  \
      (g[i][j][loc + npow(k)] - g[i][j][loc - npow(k)]);
    val = val + g[l][k][loc]*d;
  return 0.5*val;

/* computes values of individual components of Ricci curvature tensor */
double ricci(int i, int j, int loc) {
  double val = 0.0;
  double a = 0.0, b = 0.0;
  for(int l=0; l<4; l++) val += \
    christoffel(i, j, l, loc + npow(l)) - christoffel(i, j, l, loc - npow(l));
  for(int l=0; l<4; l++) val -= \
    christoffel(i, l, l, loc + npow(j)) - christoffel(i, j, l, loc - npow(j));
  for(int m=0; m<4; m++)
    for(int l=0; l<4; l++) val += \
      christoffel(i, j, m, loc)*christoffel(m, l, l, loc);
  for(int m=0; m<4; m++)
    for(int l=0; l<4; l++) val -= \
      christoffel(i, l, m, loc)*christoffel(m, j, l, loc);
  return val;

/* computes scalar curvature */
double scalar_curvature(int i, int j, int loc) {
  double gi[4][4], A[4][4], val = 0.0;
  double A2[4][4], A3[4][4], trA, trA2, trA3, det = 0.0;

  /* find metric tensor and initialize gi to zero */
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++) {
      A[i][j] = g[i][j][loc];
      gi[i][j] = 0.0;

  /* compute g^2 and g^3 */
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      for(int k=0; k<4; k++)
A2[i][j] += A[i][k]*A[k][j]; 
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      for(int k=0; k<4; k++)
A3[i][j] += A2[i][k]*A[k][j]; 

  /* use Cayley-Hamilton method to compute inverse of g */
  trA  =  A[0][0] +  A[1][1] +  A[2][2] +  A[3][3];
  trA2 = A2[0][0] + A2[1][1] + A2[2][2] + A2[3][3];
  trA3 = A3[0][0] + A3[1][1] + A3[2][2] + A3[3][3];
  for(int i=0; i<4; i++) {
    gi[i][i] = (trA*trA*trA - 3*trA*trA2 + 2*trA3)/6.0;
    for(int j=0; j<4; j++)
      gi[i][j] += A2[i][j]*trA - A3[i][j] - 0.5*A[i][j]*(trA*trA-trA2);

  /* division of gi such that (g*gi)_{00} = 1 */
  for(int k=0; k<4; k++)
    det += A[0][k]*gi[k][0];
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      gi[i][j] /= det;

  /* compute Ricci scalar */
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      val += gi[i][j]*ricci(i, j, loc);

  return val;

Sunday, 19 November 2017

Spacetime tells matter how to move, matter tells spacetime how to curve.

"Beyond the corridor of our spacetime there are infinite number of universes, each of them governed by their own set of laws."
Not general relativity, maybe next time if my computer has enough memory and time for that.

That's all Folks!
% ffmpeg -r 30 -f image2 -i %04d.png -ss 00:00:5 -vcodec libx264 -crf 2 output.mp4
clear all; close all; map = parula(256);
N    = 256;
phi  = zeros(N, N);
dphi = zeros(N, N);
ddphi= zeros(N, N);
rho  = zeros(N, N);
[x y] = meshgrid(linspace(-1, 1, N));
f = 0;
for t = 1:3*N
    A = 0.05;
    if t>N
        f = 0.01;
    x1 = cos(2*pi*f*t)*A;
    y1 = sin(2*pi*f*t)*A;
    x2 = -cos(2*pi*f*t)*A;
    y2 = -sin(2*pi*f*t)*A;
    rho = exp(-2000*((x+x1).^2+(y+y1).^2)) + exp(-2000*((x+x2).^2+(y+y2).^2));

    ddphi = rho + del2(phi);
    dphi = dphi + ddphi;
     phi(1, :)   =   phi(2, :);
    dphi(1, :)   = -dphi(2, :);
     phi(end, :) =   phi(end-1, :);
    dphi(end, :) = -dphi(end-1, :);
     phi(:, 1)   =   phi(:, 2);
    dphi(:, 1)   = -dphi(:, 2);
     phi(:, end) =   phi(:, end-1);
    dphi(:, end) = -dphi(:, end-1);

    phi = phi + dphi;
    c = phi/max(max(phi));
    surface(c, 'edgecolor', 'none');
    caxis([0 1]);
    axis([1 N 1 N]);
    axis square;
    axis off;
    shading interp;
    imwrite(c*255, map, ['png/' num2str(t, '%04.f') '.png']);

Saturday, 11 November 2017

50 Things I'd like to do if I had all the time in the world

 1. Solve quantum gravity.
 2. Crack RSA (factoring) or prove it impossible (polynomial time).
 3. Build a working quantum computer (capable of factoring RSA) or prove it impossible (a separate question from the possibility of classical factoring).
 4. Learn to play piano (I can a bit, but learn it well).
 5. Build my own CPU from transistor level up.
 6. Build my own compiler.
 7. Build my own operating system.
 8. Build my own AI that rivals humans in general creative tasks.
 9. Build my own house.
10. Build my own car.
11. Acquire a few airplanes.
12. Build my own spaceship.
13. Build my own research satellite.
14. Visit the moon.
15. P vs. NP.
16. Climb Aconcagua.
17. Climb Mt. Denali.
18. Climb Mt. Everest.
19. Learn to fly a plane.
20. Visit Mt. Thor (Baffin island).
21. Sail to Ball's Pyramid from Sydney.
22. Visit Svalbard.
23. Directly observe the event horizon of the Milky Way central black hole.
24. Solve the mystery of baryon asymmetry.
25. Solve the mystery of dark energy and dark matter.
26. Crack FTL/time travel or prove them impossible.
27. Find extraterrestrial life.
28. Learn creative drawing.
29. Walk across the Canadian wilderness (~5000 km).
30. If there is a mystery related to human consciousness, solve it, if not, make it obvious for everone there is none.
31. See the day when religions basically are no more.
32. See the day when aging and all disease are cured.
33. See the Milky Way-Andromeda collision.
34. Visit a planet inhabited by estraterrestrial life.
35. Develop more meaningful social relationships.
36. Acquire a nice pacific island and build my own winter retreat there.
37. Acquire an underground missile base from the US and make it nice.
38. Acquire some weapons and explosives for fun and distribute them to my retreats all around the world.
39. Stockpile the retreats with supplies and build bunkers.
40. Visit Antarctica and Cape Horn area, Greenland, Siberia, India, China, New Zealand, Jordan, Africa...
41. Build an internationally successful company that manufactures useful hardware, solves a number of practical problems and employs a number of good people.
42. Make a game that tells an interesting story, teaches a few things and presents practical problems requiring decent problem solving skills.
43. Work for a while as national park ranger, chemist, philosopher, politician and as a neuroscientist.
44. Work for a while on genetic research and related biological issues.
45. Work for a while for NASA, nvidia, AI research, weapons research, some startups...
46. Live a few years in the wilderness and a monestary.
47. Try living in different cities around the world.
48. Get a cabin from Norway from some nice remote mountain slope next to a small lake with a small airstrip where I can fly from Finland with my own small plane.
49. Get a nice place from somewhere along the Californian cost, not too far from San Francisco, but in a not too densely populated location either.
50. Write a book that will stir up people and teach them some critical thinking.

* The word "prove" used in colloquial sense here.

Though I didn't like this particular turn of events in that particular episode as in my opinion while it is not possible to avoid all mistakes given the initial conditions, I never the less believe that losing is due to mistakes that under perfect knowledge and logic should not exist. Choose your mistakes wisely and you'll be fine.
Unexamined life is not worth living is often taken to imply that rather than not being worth living some lives are simply missing something important. I would, however, argue that all lives are missing something important and on the other than that no single thing in any life is something so important that one should consider it essential. Most life has intrinsic value that is independent of its fullness.

It is often also said that the saying refers to one being able to say why things are important, right or wrong rather than simply insisting they are. However, this too is a bit arbitrary as in the end there never seems to be any why. We are what we are and want what we want, because that's the way we have come to be. Some goals are better served by some actions, but no action is right or wrong without reference to motivations which in end are arbitrary. Without such motivations we would not be humans, but rather simple machines simply waiting to do what we are told.

It is of course worthwhile to know how to influence other people by justifying actions and inactions for the betterment of humanity in general.

Wednesday, 1 November 2017

Origin of the power of quantum computing

Predicting the behavior of many particle quantum system requires solving the Schrödinger equation for a wavefunction that has complex dimensionality proportinal to the number of particles. Such a wavefunctions contains a value for all possible combinations of particle positions, i.e. it describes a kind of "multiverse" where everything conceivable happens simultaneously.

The system evolves as a whole and thus in order to predict phenomena such as entanglement, it is necessary to solve the time evolution of the "multiverse" rather than time evolution of separate wavefunctions corresponding to single particles. Single separate wavefunctions cannot yield predictions corresponding to observations in every circumstance.
The power of quantum computing is based on the complexity of the time evolution of this kind of multiparticle wavefunction. The classical computational resources required to do full simulations scale extremely poorly, but it is possible to demonstrate the idea with a simple code (below).

Full simulation of the time evolution of three electrons starting at rest relative to each other. Variable containing the wavefunction is a 6D grid of 26 complex points (3D double complex 26x26) and requires about 2.3 GBs of memory. Adding one more particle would make it 8D and require about 1556 GBs of memory. Simulating 30 particles in this configuration would require more bytes of memory than there are atoms in the observable universe. Many problems can be approximated effectively without conducting such simulations, but there are a few cases like quantum computers which cannot. Qubit are typically in a simple static grid rather than free space such as electrons in this simulation so larger number of qubits can still be simulated. The point at which quantum computers are expected to become more powerful than classical supercomputers (so called quantum supremacy) is a the moment 56 qubits.
n = 0.05;
N = 26;
x10 = -0.4; y10 =  0.4; vx1 = 0; vy1 = 0;
x20 =  0.4; y20 =  0.4; vx2 = 0; vy2 = 0;
x30 =  0.0; y30 = -0.4; vx3 = 0; vy3 = 0;
frequency_spread = 40;
psi = zeros(N, N, N, N, N, N, 'single');
V = zeros(N, N, N, N, N, N, 'single');
S = linspace(-1.0, 1.0, N);
for y1 = 1:N
    for x1 = 1:N
        for y2 = 1:N
            for x2 = 1:N
                for y3 = 1:N
                    for x3 = 1:N
                        envelope = exp(-frequency_spread*((S(x1)+x10)^2 + (S(y1)+y10)^2 + ...
                                                          (S(x2)+x20)^2 + (S(y2)+y20)^2 + ...
                                                          (S(x3)+x30)^2 + (S(y3)+y30)^2));
                        velocity = exp(-j*(vx1*S(x1) + vy1*S(y1) + ...
                                           vx2*S(x2) + vy2*S(y2) + ...
                                           vx3*S(x3) + vy3*S(y3)));
                        psi(x1, y1, x2, y2, x3, y3) =  envelope.*velocity;
                        V(x1, y1, x2, y2, x3, y3) = 1/sqrt(n + (S(x1)-S(x2))^2 + ...
                                                               (S(y1)-S(y2))^2)^2 + ...
                                                    1/sqrt(n + (S(x1)-S(x3))^2 + ...
                                                               (S(y1)-S(y3))^2)^2 + ...
                                                    1/sqrt(n + (S(x2)-S(x3))^2 + ...
for t = 1:500
    PSI1 = squeeze(mean(mean(mean(mean(psi, 1), 2), 3), 4));
    PSI2 = squeeze(mean(mean(mean(mean(psi, 1), 2), 5), 6));
    PSI3 = squeeze(mean(mean(mean(mean(psi, 3), 4), 5), 6));
    a = abs(PSI1).^2; a = a/max(max(a)); b = abs(PSI2).^2; b = b/max(max(b));
    c = abs(PSI3).^2; c = c/max(max(c)); val = a + b + c;
    imwrite(val/max(max(val)), ['png/' num2str(t, '%04.f') '.png']);
    z = zeros(N, N, N, N, N, N, 'single');
    z(2:N-1, :, :, :, :, :) = z(2:N-1, :, :, :, :, :) + diff(psi, 2, 1);
    z(:, 2:N-1, :, :, :, :) = z(:, 2:N-1, :, :, :, :) + diff(psi, 2, 2);
    z(:, :, 2:N-1, :, :, :) = z(:, :, 2:N-1, :, :, :) + diff(psi, 2, 3);
    z(:, :, :, 2:N-1, :, :) = z(:, :, :, 2:N-1, :, :) + diff(psi, 2, 4);
    z(:, :, :, :, 2:N-1, :) = z(:, :, :, :, 2:N-1, :) + diff(psi, 2, 5);
    z(:, :, :, :, :, 2:N-1) = z(:, :, :, :, :, 2:N-1) + diff(psi, 2, 6);
    % that's just del2(psi)

    psi = psi + 0.005*j*(z - V.*psi);


mergesort(int a[], int n) { // O(n log n)
  int *b, *c;
  int x=0, y=0, z;
  if(n>3) {
    b = (int*)malloc(4*n/2);
    c = (int*)malloc(4*n/2);
    for(z=0; z<n/2; z++) {
      b[z] = a[z];
      c[z] = a[z+n/2];
    mergesort(b, n/2);
    mergesort(c, n/2);
    for(z=0; z<n; z++)
      if(y<n/2 && x<n/2)
        if(b[x]>c[y]) a[z] = c[y++];
        else a[z] = b[x++];
        if(y>x) a[z] = b[x++];
        else a[z] = c[y++];
  } else
    if(a[0]>a[1]) { // swap a[0] and a[1] (xor swap)                            
      a[0] = a[0]^a[1];
      a[1] = a[1]^a[0];
      a[0] = a[0]^a[1];

Standard bicubic upscaler (4x).
A = imread('smb.png');
A = double(A)/255;

s = size(A);
S = 2;
for y = 1:S*s(2)-1
    for x = 1:S*s(1)-1
        r(x, y) = A(1+floor(x/S), 1+floor(y/S), 1);
        g(x, y) = A(1+floor(x/S), 1+floor(y/S), 2);
        b(x, y) = A(1+floor(x/S), 1+floor(y/S), 3);

r = r + del2(r);
g = g + del2(g);
b = b + del2(b);


[x y] = gradient(r); r = r + abs(x) + abs(y);
[x y] = gradient(g); g = g + abs(x) + abs(y);
[x y] = gradient(b); b = b + abs(x) + abs(y);
[x y] = gradient(r); r = r - abs(x) - abs(y);
[x y] = gradient(g); g = g - abs(x) - abs(y);
[x y] = gradient(b); b = b - abs(x) - abs(y);
Nabla squared filter applied in two stages to nearest neighbourhood doubler and sharpen in post prosessing. Pixel edges in 45 degree angle are totally eliminated and regular raster becomes solid color.