Devlog #83 – Fala 62: Sieć neuronów, wzrost komórek, giroid, generatywne miasto i fotografia obliczeniowa

Nowe symulacje

🧠 Sieć neuronów — dynamika Leaky Integrate-and-Fire

Model LIF

Neuron Leaky Integrate-and-Fire (LIF) to koń roboczy neuronauki obliczeniowej. Mimo swojej prostoty odtwarza kluczowe cechy prawdziwych neuronów: integrację wejść synaptycznych, spike po przekroczeniu progu oraz podprogowy zanik napięcia z powrotem do spoczynku. Potencjał błonowy V każdego neuronu podlega równaniu:

τ_m · dV/dt = -(V - V_rest) + R · I(t)

gdzie:
  τ_m   = stała czasowa błony (typowo 10–20 ms)
  V_rest = potencjał spoczynkowy (−65 mV)
  R     = opór błony
  I(t)  = całkowity prąd synaptyczny w chwili t

Reguła spike'a:
  jeśli V ≥ V_thresh:
      zarejestruj spike; V ← V_reset; wejdź w refrakcję na t_ref ms

Symulacja dyskretyzuje to metodą Eulera w przód przy dt = 1 ms — wystarczająco grubo, by działać w czasie rzeczywistym, ale wystarczająco dokładnie, by wiernie uchwycić odstępy międzyimpulsowe. Przy τ_m = 10 ms swobodny zanik na krok wynosi exp(−dt/τ_m) ≈ 0,905, więc neuron nieotrzymujący żadnego wejścia traci około 9,5% swojego nadmiarowego napięcia co milisekundę.

Model synaptyczny i równowaga E/I

Każdy spike z neuronu j dostarcza impuls prądu do neuronu i po ustalonym opóźnieniu synaptycznym d_ij (1–5 ms). Prąd zanika wykładniczo ze stałą czasową synapsy τ_syn:

I_syn(t) = w_ij · exp(-(t - t_spike - d_ij) / τ_syn)   [dla t > t_spike + d_ij]

Synapsa pobudzająca: w_ij > 0  (podobna do AMPA, τ_syn ≈ 5 ms)
Synapsa hamująca: w_ij < 0  (podobna do GABA, τ_syn ≈ 10 ms)

Stosunek E/I: 80% neuronów pobudzających, 20% hamujących (prawo Dale'a)

Równowaga E/I ma kluczowe znaczenie. Gdy dominuje pobudzenie, sieć wchodzi w stan niekontrolowanego odpalania (wyładowania podobne do padaczkowych). Gdy dominuje hamowanie, aktywność zanika. Preset zbalansowany (mniej więcej równy całkowity napęd pobudzający i hamujący) daje asynchroniczny nieregularny stan (AI), przypominający korę mózgową podczas aktywnej czuwania.

Mapa cieplna napięcia i raster spike'ów

Lewy panel renderuje mapę cieplną 2D, w której każdy piksel odpowiada neuronowi (ułożonemu w siatce), zabarwioną wg napięcia błonowego: ciemnoniebieski przy V_rest, przechodząc przez cyjan i żółty do jasnobiałego przy progu. Prawy panel pokazuje przewijalny raster spike'ów — każda kropka to jeden spike jednego neuronu, naniesiony w chwili wystąpienia. Raster pozwala na natychmiastową wizualną diagnostykę stanu sieci: rozproszone losowe kropki wskazują stan AI, poziome pasma synchronicznych kropek wskazują aktywność oscylacyjną lub napadową.

Statystyki synchronii i częstotliwości

Dwie statystyki na żywo są obliczane w każdej klatce. Częstotliwość odpalania populacji r(t) to liczba spike'ów w ostatnich 50 ms podzielona przez N × 0,050 s, dając częstotliwość w Hz. Synchronia mierzona jest jako współczynnik zmienności odstępu międzyimpulsowego (CV_ISI), uśredniony po neuronach: CV = σ_ISI / μ_ISI. CV ≈ 1 wskazuje na odpalanie podobne do Poissona (AI); CV ≪ 1 to regularne (zegarowe); CV ≫ 1 wskazuje na aktywność napadową. Preset oscylacyjny wpędza sieć w rytm gamma ~40 Hz, widoczny jako rytmiczne kolumny na rastrze i CV ≈ 0,3.

🦠 Wzrost komórek i morfogeneza — podział, pakowanie i formowanie wzoru

Pakowanie komórek oparte na siłach

Komórki modelowane są jako miękkie dyski o pozycji x_i i promieniu r_i. W każdym kroku czasu na każdą komórkę działają trzy siły:

F_push(i,j)     = k_push  · max(0, r_i + r_j - |x_i - x_j|) · n̂_ij   // odpychanie nakładania
F_adhesion(i,j) = k_adh   · max(0, |x_i - x_j| - (r_i+r_j+d_adh)) · n̂_ij  // przyciąganie krótkiego zasięgu
F_boundary(i)   = k_wall  · max(0, r_i - dist_to_boundary) · n̂_wall   // miękka ściana kontenera

Wypadkowa: x_i += (ΣF) · dt / ζ    (nadtłumiona; bez bezwładności)

Siatka przestrzenna przyspiesza wyszukiwanie sąsiadów do O(n) poprzez haszowanie komórek do koszyków siatki o rozmiarze 2 r_max. Tylko komórki w tym samym lub sąsiednim koszyku wymagają obliczenia siły — utrzymuje to symulację na poziomie 60 fps nawet blisko limitu 600 komórek.

Reakcja-dyfuzja Graya-Scotta

Równolegle, ciągła siatka 80×60 ewoluuje zgodnie z równaniami Graya-Scotta dla dwóch związków chemicznych U i V:

∂U/∂t = D_U · ∇²U  -  U·V²  +  f·(1 - U)
∂V/∂t = D_V · ∇²V  +  U·V²  -  (f + k)·V

Typowe parametry dla wzoru plamek:
  D_U = 0.16,  D_V = 0.08,  f = 0.035,  k = 0.065

Laplasjan (5-punktowy szablon, granice periodyczne):
  ∇²U[i,j] = U[i-1,j] + U[i+1,j] + U[i,j-1] + U[i,j+1] - 4·U[i,j]

Współczynnik podaży f i współczynnik zabijania k determinują, jaki wzór Turinga powstaje: plamki, paski, labiryntowe struktury lub chaos czasoprzestrzenny. Symulacja mapuje pozycje komórek na ciągłą siatkę, aby próbkować lokalne stężenie V, które steruje losem komórki: wysokie V → wyspecjalizowany (zróżnicowany) los; niskie V → macierzysty.

Podział komórki i tryby koloru

Komórka rośnie ze stałą szybkością dr/dt = growth_rate. Gdy r > r_divide, dzieli się na dwie komórki potomne pod kątem ±45° od losowej orientacji, każda o r = r_divide / √2 (zachowując pole powierzchni). Zdarzenie podziału generuje krótkozasięgowy impuls ciśnienia, który propaguje się przez kolonię poprzez system sił, tworząc charakterystyczną dynamikę pakowania komórek obserwowaną w rosnących nabłonkach.

• Tryb losu — kolor koduje stężenie V Graya-Scotta próbkowane w pozycji komórki na siatce. Gradienty morfogenu dają przestrzennie spójne domeny kolorystyczne.
• Tryb wieku — kolor koduje numer pokolenia komórki (czas od ostatniego podziału), od młodego (cyjan) do starego (magenta).
• Tryb rozmiaru — kolor koduje aktualny promień względem progu podziału, pokazując, które komórki zaraz się podzielą (jasnobiały).

💠 Giroid — potrójnie okresowe powierzchnie minimalne

Czym jest TPMS?

Potrójnie okresowa powierzchnia minimalna (TPMS) to powierzchnia okresowa we wszystkich trzech kierunkach przestrzennych, mająca zerową krzywiznę średnią w każdym punkcie (H = 0, definiująca właściwość powierzchni minimalnej). Natura wytwarza te struktury samorzutnie: pojawiają się w łuskach skrzydeł motyli, membranach kopolimerów blokowych oraz w systemach dwuwarstw lipidowych. Giroid w szczególności jest strukturą biciągłą — obie strony powierzchni tworzą połączone labirynty — co czyni go wysoko cenionym w inżynierii materiałowej do filtracji, magazynowania energii i zastosowań fotonicznych.

Równania powierzchni uwikłanej

Każda powierzchnia jest zdefiniowana równaniem uwikłanym F(x, y, z) = t, gdzie t jest izowartością (domyślnie t = 0 daje dokładną powierzchnię minimalną). We współrzędnych przeskalowanych (u = 2πx/a, v = 2πy/a, w = 2πz/a):

Giroid:   sin(u)cos(v) + sin(v)cos(w) + sin(w)cos(u) = t
Schwarz P: cos(u) + cos(v) + cos(w) = t
Diament:  sin(u)sin(v)sin(w) + sin(u)cos(v)cos(w)
        + cos(u)sin(v)cos(w) + cos(u)cos(v)sin(w) = t
Neovius:  3(cos(u)+cos(v)+cos(w)) + 4cos(u)cos(v)cos(w) = t

Giroid został odkryty przez Alana Schoena w 1970 roku podczas badania NASA nad lekkimi strukturami. W przeciwieństwie do powierzchni Schwarza P i D, nie ma prostych linii ani płaszczyzn symetrii lustrzanej, co nadaje mu chiralny, labiryntowy charakter.

Algorytm maszerujących sześcianów

Symulacja oblicza F na regularnej siatce N³ wokseli (domyślnie N = 48) i wyodrębnia izopowierzchnię F = t, korzystając z 256-przypadkowej tabeli wyszukiwania maszerujących sześcianów Bourke'a. Każdemu sześcianowi woksela przypisywany jest 8-bitowy indeks (jeden bit na wierzchołek, ustawiony, gdy F > t) — indeks wskazuje jedną z 256 konfiguracji trójkątów, dając do 5 trójkątów na sześcian. Wierzchołki są umieszczane na krawędziach sześcianu przez interpolację liniową:

// Dla krawędzi od wierzchołka A (wartość fA) do wierzchołka B (wartość fB):
t_lerp = (isovalue - fA) / (fB - fA)
wierzchołek  = A + t_lerp · (B - A)

Łączna liczba trójkątów: O(N² · pole powierzchni w wokselach)
Dla N=48, kafelkowanie=1: typowo ~80 000–120 000 trójkątów

Normalne wierzchołków obliczane są analitycznie z ∇F, a nie przez uśrednianie normalnych ścian — daje to idealnie gładkie cieniowanie przy dowolnej rozdzielczości. Siatka jest wgrywana do Three.js BufferGeometry i renderowana z MeshStandardMaterial pod mapą otoczenia, dając metaliczny/opalizujący wygląd.

Sterowanie kafelkowaniem i izowartością

Suwak kafelkowania powtarza komórkę elementarną 1–4 razy w każdej osi, demonstrując okresowy charakter powierzchni. Zwiększanie izowartości t przesuwa powierzchnię w stronę jednej strony biciągłego podziału: przy t = 0 oba labirynty mają równą objętość; przy t > 0 jeden labirynt kurczy się i ostatecznie zamyka, przechodząc od TPMS do topologii pianki o zamkniętych komórkach.

🌃 Generatywne miasto — proceduralny krajobraz miejski

Rekurencyjny podział kwartałów

Układ miasta generowany jest przez rekurencyjny podział prostokątnego płótna. Zaczynając od pełnego prostokąta płótna, algorytm dzieli go wzdłuż osi poziomej lub pionowej, dodając niewielki losowy jitter (±10% długości boku) do pozycji podziału, aby uniknąć wrażenia siatki. Podział trwa, aż szerokość lub wysokość kwartału spadnie poniżej minimalnego rozmiaru działki. Oś podziału zmienia się deterministycznie na podstawie proporcji bieżącego prostokąta: zawsze dzieli się dłuższy bok.

function subdivide(rect, depth):
    jeśli rect.w < MIN_LOT lub rect.h < MIN_LOT: zwróć [rect]
    axis = rect.w > rect.h ? 'x' : 'y'
    split = lerp(0.35, 0.65, lcg_random()) * (axis='x' ? rect.w : rect.h)
    [left, right] = split_rect(rect, axis, split + jitter())
    zwróć subdivide(left, depth+1) + subdivide(right, depth+1)

Daje to nieregularny, ale realistyczny wzorzec kwartałów, przypominający organiczny wzrost miasta, a nie sztywną siatkę kartezjańską.

Strefowanie i gradient gęstości

Każdy kwartał liściowy otrzymuje strefę na podstawie odległości od centrum miasta (znormalizowanej do [0, 1]):

• Handlowa (d < 0,25) — gęste wysokie budynki; prostokąty wypełniają większość kwartału, pozostawiając jedynie wąską szczelinę uliczną. Kolory: ciepłe szarości i błękity łupkowe.
• Mieszkalna (0,25 ≤ d < 0,65) — mniejsze budynki z ogrodami. Kwartały zawierają 2–8 obrysów domów ze wspólną szczeliną podwórka. Kolory: ciepłe brązy.
• Park (losowe 12% kwartałów środkowego pierścienia) — wypełniony kropkowanym zielonym wzorem. Zapewnia wizualną przestrzeń i realistyczny teren otwarty.
• Woda (ciągłe plamy w zewnętrznym pierścieniu, ~8% powierzchni) — renderowana z animowanymi liniami migotania. Rzeki i jeziora powstają ze złączenia sąsiadujących kwartałów wodnych.

Wysokość budynku kodowana jest jasnością koloru: wyższe budynki otrzymują jaśniejsze wypełnienie i białą smugę-„czapkę”, dając widokowi z góry niejawną wskazówkę 3D.

Losowość LCG z ziarnem dla odtwarzalności

Cała losowość przepływa przez liniowy generator kongruentny zasiany wartością ziarna użytkownika:

// parametry LCG (Numerical Recipes)
a = 1664525;  c = 1013904223;  m = 2^32
seed = (a * seed + c) % m
random() { seed = ...; zwróć seed / m; }

To samo ziarno zawsze daje to samo miasto. Cztery motywy (noc, dzień, cyberpunk, zachód słońca) zmieniają jedynie paletę kolorów, nie układ — można więc porównać, jak ta sama struktura miejska wygląda o różnych porach dnia.

🖼️ Fotografia obliczeniowa — cztery algorytmy przetwarzania obrazu

Sztuka ASCII poprzez mapowanie luminancji

Obraz jest próbkowany w dół do siatki znaków (domyślnie 120 kolumn). Średnia luminancja każdej komórki L = 0,299R + 0,587G + 0,114B (ITU-R BT.601) mapowana jest na jeden z 70 znaków posortowanych według gęstości wizualnej od najciemniejszego do najjaśniejszego:

charset = "$@B%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/\\|()1{}[]?-_+~<>i!lI;:,\"^`'. "
char_index = floor(L / 255 * (charset.length - 1))
// Renderowanie czcionką o stałej szerokości; kolor z oryginalnego piksela dla trybu kolorowego ASCII

Fotomozaika przez uśrednianie koloru kafelka

Efekt fotomozaiki dzieli obraz na prostokątne kafelki (rozmiar regulowany od 4×4 do 32×32 pikseli) i zastępuje każdy kafelek jednolitym wypełnieniem jego średnim kolorem. Jest to cyfrowy odpowiednik klasycznej mozaiki fotograficznej, a przy dużych rozmiarach kafelków funkcjonuje jako filtr pikselizacji/prywatności. Uśrednianie wykorzystuje pełne skanowanie wszystkich pikseli w każdym kafelku — a nie przybliżenie próbkowane — zapewniając dokładność koloru.

Kropkowanie metodą dysku Poissona

Kropkowanie zamienia obraz w zbiór kropek, których gęstość koduje jasność. Symulacja wykorzystuje przybliżenie metodą próbkowania odrzuceniowego: pozycje kandydackich kropek losowane są jednostajnie; kropka jest umieszczana tylko wtedy, gdy lokalna luminancja (odwrócona, więc ciemne obszary otrzymują więcej kropek) przekracza próg. Minimalna odległość między kropkami wymusza rozkład dysku Poissona, zapobiegając zlepianiu:

dla każdego kandydata (x, y):
    L = luminance(x, y)
    density = 1 - L/255           // ciemny piksel → wysoka gęstość
    jeśli random() < density i no_dot_within_r(x, y, r_min):
        umieść kropkę w (x, y)

Dithering z dyfuzją błędu Floyda-Steinberga

Klasyczny dithering Floyda-Steinberga daje reprezentację 1-bitową (czarno-białą), która percepcyjnie odpowiada ciągłemu oryginałowi. Algorytm skanuje piksele od lewej do prawej, od góry do dołu. Każdy piksel jest kwantyzowany do czerni (0) lub bieli (255); błąd kwantyzacji rozprowadzany jest do czterech sąsiadów:

error = original_pixel - quantised_pixel

pixel[x+1][y  ] += error × 7/16
pixel[x-1][y+1] += error × 3/16
pixel[x  ][y+1] += error × 5/16
pixel[x+1][y+1] += error × 1/16

Błąd propaguje się do przodu i w dół, więc ciemne obszary kumulują błąd dodatni (sprawiając, że późniejsze piksele są bardziej skłonne do zaokrąglenia do bieli) i odwrotnie. Wynikiem jest wzór podobny do rastra półtonowego z wysokoczęstotliwościowym szumem ditheringu, który z normalnej odległości oglądania wygląda jak odcień szarości odpowiadający oryginalnej luminancji. Dithering Floyda-Steinberga został opublikowany w 1976 roku i pozostaje standardowym algorytmem w edytorach obrazów i drukarkach na całym świecie.

Co dalej

Fala 63 przechodzi do dynamiki płynów, matematyki kombinatorycznej i optyki. Nadchodzi: pełzający przepływ Stokesa z demonstracją odwracalności, współczynniki upakowania kul w 2D i 3D, wypełniające przestrzeń krzywe Hilberta i Peana, spirala Ulama liczb pierwszych, maksymalny przepływ Edmondsa-Karpa z animowanymi ścieżkami powiększającymi oraz raytracer halo lodowego metodą Monte Carlo odtwarzający halo 22° i słońca poboczne.