WebGPUとRustでネコチャンを点描した話

69370 ワード

テキストリンク

はじめに

昨年１２月にこんなツイートを見かけました。

かわいいですね。幸いにして論文と実装が公開されていたので、自分でもやってみようと思って、その結果を書いたのが前回の記事です。

読んでいただければわかるとおり、前回の記事の中ではGPUを使わずにアルゴリズムやデータ構造を工夫して近似的に計算しています。結果の画像についてはそんなに悪くないと思っていますが、限界もありました。パーティクルの数も少ないし、一部の画像ではうまく行きませんでした。

やっぱり、もっとちゃんとネコチャンを点描してみたいですよね。

なので、今回の記事ではGPUを使ってアルゴリズムを再現し、よりヴィヴィッドなネコチャンの点描を作成しようと思います。GPUを使って計算するために、WebGPUのRust実装であるwgpuを使用します。ネコチャンの画像を点描にしたい人と、WebGPUに入門してcompute shaderで何か計算していみたいという人には参考になると思います。

なおこの記事ではWebの話はしません。WASMの話もしません。

WebGPU

WebGPUってなんだ？　という人もいるでしょう。簡単に言うと、W3Cによって作成されているブラウザからGPUを操作するためのWeb APIです。こうしたAPIとしてはWebGL2.0が有名で普及もしていますが、WebGPUはWebGL2.0の後継となることを目指して作成されているものです。歴史は比較的浅く現在はまだWorking Draftのステータスですが、次世代の標準になる（かもしれない）APIといったところでしょうか。

WebGPUにはwgpuというRustの実装があります。これはFirefoxのバックエンドで使用されているものです（ちなみに前回の記事で利用したnannouもバックエンドでwgpuを使っています）。今回はこのwgpuというクレートを使って実装します。

点描化のアルゴリズム

点描化のアルゴリズムについては前回の記事にも書いたので、ここでは簡単な説明にとどめます。

電荷を持った粒子（パーティクル）の動きを以下の条件のもとにシミュレーションします。

パーティクルは負の電荷を持ち、互いに反発する
画像はグレースケールに変換する。各ピクセルはグレースケールの数値に比例する正の電荷を持ち、パーティクルを引きつける

この二つの力を同時に計算することで、色の濃い領域ほどたくさんのパーティクルが集まり、一方でパーティクル同士が近づくほど反発する力が強くなるので、程よい密度になるよう距離が保たれるという効果を生み出すことができます。

前回の記事では計算量を減らすためにパーティクルの近くのピクセルだけを計算に使っていました。パーティクルの相互作用も近傍だけ計算できるようにデータ構造を工夫しています。

今回はもちろんすべてのピクセルとパーティクルを使って計算します。

実装について

wgpuはmasterブランチのものを使います。crates.ioに登録されている最新版はv0.12ですが、このバージョンではwgsl（WebGPU Shading Language）のいくつかの記法がコンパイルできなかったためです。

# Cargo.toml
[dependencies]
wgpu = { git = "https://github.com/gfx-rs/wgpu", branch = "master" }
...

実装の方針

詳しい実装を見る前に、大まかな方針について説明します。処理の流れを素直に書くと以下の様になります

1. 画像を読み込む
2. GPUとのインターフェースを用意する
3. 画像をグレースケールに変換する
4. loop {
5.   パーティクルが画像から受ける力を計算する
6.   パーティクル相互の力を計算する
7.   パーティクルの位置を更新する
8.   描画する
9. }

このまま実装してもよいのですが、工夫の余地があります。

「5. パーティクルが画像から受ける力を計算する」の処理は結構重い処理です。全ピクセル分の力を計算するとなると、500x500ピクセルの画像なら250000回計算する必要があります。全パーティクルについて都度この計算を行うと、パーティクルの位置の更新に結構時間がかかります。

ですが実はこの部分、繰り返しの中で計算する必要はないのです。というのも画像から受ける力はパーティクルの位置と画像にのみ依存して決まり、他のパーティクルの位置を参照する必要がないからです。だから、あらかじめすべての位置に対して画像から受ける力を計算し保存しておくことで（キャッシュしておくことで）、パーティクルの位置更新時には保存した力を参照するだけで良くなります。

このことをふまえ、次のような流れで実装します

1. 画像を読み込む
2. GPUとのインターフェースを用意する
3. 画像をグレースケールに変換する
4. 画像から受ける力を計算し、保存する
5. loop {
6.   パーティクルが画像から受ける力を保存した結果から読み出す
7.   パーティクル相互の力を計算する
8.   パーティクルの位置を更新する
9.   描画する
10. }

以下、具体的な実装について主要な部分に絞って説明しようと思います。

deviceとqueue

GPUとやり取りするインターフェースとなるのがdeviceです。queueは処理を実行したりする時に使います。

let instance = wgpu::Instance::new(wgpu::Backends::all());
let surface = unsafe { instance.create_surface(window) };
let adapter = instance
    .request_adapter(&wgpu::RequestAdapterOptions {
        power_preference: wgpu::PowerPreference::default(),
        compatible_surface: Some(&surface),
        force_fallback_adapter: false,
    })
    .await
    .unwrap();
let (device, queue) = adapter
    .request_device(
        &wgpu::DeviceDescriptor {
            label: None,
            features: wgpu::Features::empty(),
            limits: wgpu::Limits::default(),
        },
        None,
    )
    .await
    .unwrap();

deviceとqueueと通じて、データ領域を確保したり、シェーダーを実行したりします。

グレースケールへの変換

Texture

画像データはimageクレートを使って読み込みます。

let img_bytes = include_bytes!("../assets/cat.png");
let input_image = image::load_from_memory(img_bytes).unwrap();

ただし、このままではwgpuでは利用できません。画像データはTextureというフォーマットでGPUと受け渡しする必要があるので、Textureを作成します。

let input_texture = device.create_texture(&wgpu::TextureDescriptor {
    label: Some("input"),
    size: texture_size,
    mip_level_count: 1,
    sample_count: 1,
    dimension: wgpu::TextureDimension::D2,
    format: wgpu::TextureFormat::Rgba8UnormSrgb,
    usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST,
});

queueのメソッドを使って先程の画像データをTextureに書き込みます。

let (width, height) = input_image.dimensions();
queue.write_texture(
    input_texture.as_image_copy(),
    &input_image.to_rgba8(),
    wgpu::ImageDataLayout {
        offset: 0,
        bytes_per_row: std::num::NonZeroU32::new(4 * width),
        rows_per_image: std::num::NonZeroU32::new(height),
    },
    texture_size,
);

グレースケールに変換した結果を格納するTextureも用意します。

let gray_texture = device.create_texture(&wgpu::TextureDescriptor {
    label: Some("gray texture"),
    size: texture_size,
    mip_level_count: 1,
    sample_count: 1,
    dimension: wgpu::TextureDimension::D2,
    format: wgpu::TextureFormat::Rgba32Float,
    usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::STORAGE_BINDING,
});

usageにSTORAGE_BINDINGを指定しておくことで、結果を格納するTextureとして利用できます。

pipeline

pipelineは読み込んだシェーダーやデータをGPUに渡すためのものです。

let pipeline = device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
    label: Some("grayscale pipeline"),
    layout: None,
    module: &grayscale_shader,
    entry_point: "main",
});

シェーダー

ここで、シェーダーの中身も見てみましょう。シェーダーはwgsl（WebGPU Shading Language）というDSLで記述します。

@group(0) @binding(0) var input_texture : texture_2d<f32>;
@group(0) @binding(1) var output_texture : texture_storage_2d<rgba32float, write>;

@stage(compute)
@workgroup_size(16, 16)
fn main(@builtin(global_invocation_id) global_invocation_id : vec3<u32>) {
    let dimensions = textureDimensions(input_texture);
    let coords = vec2<i32>(global_invocation_id.xy);
    if (coords.x >= dimensions.x || coords.y >= dimensions.y) {
        return;
    }

    let color = textureLoad(input_texture, coords.xy, 0);
    let gray = dot(vec3<f32>(0.299, 0.578, 0.114), color.rgb);

    textureStore(output_texture, coords.xy, vec4<f32>(gray, gray, gray, color.a));
}

残念ながらzennではシンタックスハイライトが効かないようなのですが、それはともかく、wgslを全く知らないという人でもどんな処理が実装されているか、おおよその想像はつくのではないでしょうか。input_textureから色情報を取り出し、内積をとってグレースケールに変換し、output_textureに格納するというだけです。

わかりにくいのは@group、@binding、@workgroup_sizeあたりの記述ではないかと思います。これらを理解するために、再びRustの世界に戻りましょう。

BindGroup

BindGroupはTextureなどのデータを格納します。BindGroupはpieplineを通じてGPUに渡されます。

let bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
    label: Some("gray scale bind group"),
    layout: &pipeline.get_bind_group_layout(0),
    entries: &[
        wgpu::BindGroupEntry {
            binding: 0,
            resource: wgpu::BindingResource::TextureView(
                &input_texture.create_view(&wgpu::TextureViewDescriptor::default()),
            ),
        },
        wgpu::BindGroupEntry {
            binding: 1,
            resource: wgpu::BindingResource::TextureView(
                &output_texture.create_view(&wgpu::TextureViewDescriptor::default()),
            ),
        },
    ],
});

entries にBindGroupEntryの配列を渡していますね。BindGroupEntryの中にbindingという要素が見えるかと思いますが、これが先程見たシェーダーの中の@bindingの数値に対応しています。binding: 0にinpute_textureを渡しておけば、シェーダーの中で@binding(0)でinput_textureを受け取ることができるというわけです。

では@groupは何なのかと思われるかもしれません。実はBindGroupはpipelineに対して複数セットできます。@groupの数値はどのBindGroupかを識別するためのものです。

グレースケールへの変換

GPUで画像を変換しましょう。

let mut encoder =
    device.create_command_encoder(&wgpu::CommandEncoderDescriptor { label: None });
{
    let dispatch_width = (texture_size.width + 15) / 16;
    let dispatch_height = (texture_size.height + 15) / 16;
    let mut compute_pass =
        encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { label: None });
    compute_pass.set_pipeline(&pipeline);
    compute_pass.set_bind_group(0, &bind_group, &[]);
    compute_pass.dispatch(dispatch_width, dispatch_height, 1);
}
queue.submit(Some(encoder.finish()));

dispatch_width、dispatch_heightというものを計算しています。これは先程のシェーダーに出てきた@workgroup_sizeに関連するものです。@workgroup_size(16, 16)とは要するに16x16ピクセルずつ並列で処理するということを表しています。例えば1000x1000ピクセルの画像を処理するのであれば、縦横63(=ceil(1000 / 16))個のworkgroupが作成されることになります。compute_pass.dispatchにはworkgroupの数を渡す必要があるので、それを計算しているというわけです（compute_pass.dispatchの最後の引数はz軸のworkgroup数です。今回は使いません）

画像のキャッシュを計算する

つづいて、画像のキャッシュを計算する部分です。

まず、シミュレーションのパラメータを格納するBufferを用意しています。BufferはTextureと同様GPUに渡すためのデータを格納するものです。Textureは画像データ専用ですが、Bufferは汎用的に用いられます。

let param = vec![BLANK_LEVEL, NUM_PARTICLES as f32 * Q_CHARGE];
let param_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
    label: Some("cache param buffer"),
    contents: bytemuck::cast_slice(&param),
    usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST,
});

シェーダーは少し長いですが、次のようになります。

struct Params {
    blank_level: f32,
    dot_total_charge: f32
};

@group(0) @binding(0) var<uniform> params : Params;
@group(0) @binding(1) var input_texture : texture_2d<f32>;
@group(0) @binding(2) var output_texture : texture_storage_2d<rgba32float, write>;

fn rotate(v: f32) -> f32 {
    if (v > 1.0) {
        return v - 2.0;
    } else if (v < -1.0) {
        return v + 2.0;
    } else {
        return v;
    }
}

fn coord_to_pos(x: i32, y: i32, width: i32, height: i32) -> vec2<f32> {
    let pos = vec2<f32>(
        2.0 * f32(x) / f32(width) - 1.0,
        1.0 - 2.0 * f32(y) / f32(height)
    );
    return pos;
}

@stage(compute)
@workgroup_size(16, 16)
fn main(
    @builtin(global_invocation_id) global_invocation_id: vec3<u32>
) {
    let dimensions = textureDimensions(input_texture);
    let width = dimensions.x;
    let height = dimensions.y;
    let coords = vec2<i32>(global_invocation_id.xy);
    let Epsilon = 0.0000001;

    if (coords.x >= width || coords.y >= height) {
        return;
    }

    let pos = coord_to_pos(coords.x, coords.y, width, height);

    var acc : vec2<f32> = vec2<f32>(0.0, 0.0);
    var total_q : f32 = 0.0;
    var ix: i32 = 0;
    loop {
        if (ix >= width) {
            break;
        }
        var iy: i32 = 0;
        loop {
            if (iy >= height) {
                break;
            }
            let texel = textureLoad(input_texture, vec2<i32>(ix, iy), 0);
            let gray = texel[0];
            let bmpQ = params.blank_level - gray;
            total_q += bmpQ;

            let texpos = coord_to_pos(ix, iy, width, height);
            var dp : vec2<f32> = texpos - pos;
            dp.x = rotate(dp.x);
            dp.y = rotate(dp.y);

            // ゼロ除算を防ぐため適当に小さい数を足す
            let d2 = dot(dp, dp) + 0.00003;
            let d = sqrt(d2);

            if (abs(d) > Epsilon) {
                // 距離の二乗に反比例する力を足している
                let q = bmpQ / d2;
                acc += q * dp / d;
            }
            continuing {
                iy = iy + 1;
            }
        }
        continuing {
            ix = ix + 1;
        }
    }
    acc *= params.dot_total_charge / total_q;
    textureStore(output_texture, coords.xy, vec4<f32>(acc, f32(coords.x), f32(coords.y)));
}

「距離の二乗に反比例する力を足している」とコメントをつけた箇所がアルゴリズムに直接関連する部分で、ほかはデータの準備や後続の処理につなぐための後処理です。acc *= params.dot_total_charge / total_q;の部分で画像の総電荷とパーティクルの総電荷が釣り合うように調整しています。計算の結果はもはや画像として解釈できるものではありませんが、textureStoreでTextureに格納しています。

loopの記述が冗長だな、と思った人もいるかもしれません。wgslのドキュメントにはfor文についての記載があり、そのうち利用できるようになるとは思うのですが、今のところ最新のwgpuでもコンパイルできないようなので、やむを得ずloopを使用しています。

シミュレーション

では、いよいよ本命のシミュレーションを実装していきましょう。

必要なのは

上で計算した結果を格納したTexture
シミュレーションのパラメータを格納するバッファ
パーティクルのデータ（位置、速度、加速度）を格納するバッファ
計算用のパイプライン

です。バッファやパイプラインの作り方はこれまでと同じですが、パーティクルのバッファだけやや込み入ったことをする必要があります。以下で詳しく説明します。

パーティクルのバッファ

今から説明する仕組みはwgpuのexampleを参考にしています。

パーティクルは位置、速度、加速度を状態として持ちます。それぞれx軸、y軸の量があるので6次元の配列で表現できますね。例えばパーティクルが1024個なら6x1024の長さの配列を用意して初期化すれば良いわけです。

では、それだけのデータを格納するBufferを用意すればよいかというと、それでは不十分です。問題はパーティクルの位置を更新する際に起こります。in-placeに位置を更新すると、後から計算するパーティクルでは計算がずれてしまうのです。だから「計算元のデータを格納するBuffer」と「計算結果を格納するBuffer」の二つを用意しなくてはなりません。

単純な話ですが、実現するための仕組みは少々複雑です。実際のコードを見てみましょう。

パーティクルの初期化とBufferへの格納は以下のとおりです。位置、速度、加速度の順で格納します。位置のみ-1~1の範囲の乱数で初期化しています。

let mut initial_particle_data: Vec<f32> = vec![0.0; (6 * NUM_PARTICLES) as usize];
let mut rng = rand::rngs::StdRng::seed_from_u64(111);
let unif = Uniform::new_inclusive(-1.0, 1.0);
for particle_instance_chunk in initial_particle_data.chunks_mut(6) {
    particle_instance_chunk[0] = unif.sample(&mut rng); // 位置x座標
    particle_instance_chunk[1] = unif.sample(&mut rng); // 位置y座標
    particle_instance_chunk[2] = 0.0;                   // 速度x成分
    particle_instance_chunk[3] = 0.0;                   // 速度y成分
    particle_instance_chunk[4] = 0.0;                   // 加速度x成分
    particle_instance_chunk[5] = 0.0;                   // 加速度y成分
}

let mut particle_buffers = Vec::<wgpu::Buffer>::new();
for i in 0..2 {
    particle_buffers.push(
        device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
            label: Some(&format!("Particle Buffer {}", i)),
            contents: bytemuck::cast_slice(&initial_particle_data),
            usage: wgpu::BufferUsages::VERTEX
                | wgpu::BufferUsages::STORAGE
                | wgpu::BufferUsages::COPY_DST
                | wgpu::BufferUsages::MAP_READ,
        }),
    );
}

particle_buffersというBufferを格納する配列を作成し、まったく同じデータを格納したBufferを二つpushしています。二つのBufferは説明のためBufferA、BufferBと呼ぶことにしましょう。particle_buffersには[BufferA, BufferB]の順で格納されています。

つづいてBindGroupです。

let mut particle_bind_groups = Vec::<wgpu::BindGroup>::new();
for i in 0..2 {
    particle_bind_groups.push(device.create_bind_group(&wgpu::BindGroupDescriptor {
        layout: &compute_bind_group_layout,
        entries: &[
            wgpu::BindGroupEntry {
                binding: 0,
                resource: particle_buffers[i].as_entire_binding(),
            },
            wgpu::BindGroupEntry {
                binding: 1,
                resource: particle_buffers[(i + 1) % 2].as_entire_binding(),
            },
        ],
        label: None,
    }));
}
return particle_bind_groups;

やはりparticle_bind_groupsというBindGroupを格納する配列を作成し、そのなかに二つBindGroupを格納しています。ただし今度は全く同じものを格納するわけではありません。entriesの中身のインデックスに注目すればわかるとおり、１つ目のBindGroupにはBufferA、BufferB、２つ目のBindGroupにはBufferB、BufferAという順序でBufferが渡されています。

ここで渡されるのはいうならばBufferの参照です。だから例えばBufferAが更新されてBufferA'になると、１つ目のBindGroupの中身はBufferA'、BufferB、２つ目のBindGroupはBufferB、BufferA'というふうに更新されます。

この二つのBindGroupを交互に利用することで、BufferAとBufferBを交互に「計算元のバッファ」と「計算結果を格納するバッファ」として機能させることができるのです。

どういうことか。位置の更新と描画を実際に行う部分をコードで見てみましょう。まず、位置の更新をするのが以下の部分

let mut command_encoder = self
    .device
    .create_command_encoder(&wgpu::CommandEncoderDescriptor { label: None });
{
    let mut compute_pass =
        command_encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { label: None });
    compute_pass.set_pipeline(&self.compute_pipeline);
    compute_pass.set_bind_group(0, &parameter_bind_group, &[]);
    compute_pass.set_bind_group(1, &self.particle_bind_groups[self.frame_num % 2], &[]); // particle_bind_groupを渡している
    compute_pass.set_bind_group(2, &self.texture_bind_group, &[]);
    compute_pass.dispatch(self.work_group_count, 1, 1);
}
// 以下render用の処理

この処理はloopの中で呼ばれます。self.frame_numはloopの中でインクリメントするカウンタで、二つあるparticle_bind_groupがframe_numに応じて交互に渡されることになるのがわかると思います。

描画は以下の部分

let frame = self.surface.get_current_texture()?;
let view = frame
    .texture
    .create_view(&wgpu::TextureViewDescriptor::default());
let color_attachments = [wgpu::RenderPassColorAttachment {
    view: &view,
    resolve_target: None,
    ops: wgpu::Operations {
        load: wgpu::LoadOp::Clear(wgpu::Color {
            r: 1.0,
            g: 1.0,
            b: 1.0,
            a: 1.0,
        }),
        store: true,
    },
}];
let render_pass_descriptor = wgpu::RenderPassDescriptor {
    label: None,
    color_attachments: &color_attachments,
    depth_stencil_attachment: None,
};
{
    let mut render_pass = command_encoder.begin_render_pass(&render_pass_descriptor);
    render_pass.set_pipeline(&self.render_pipeline);
    render_pass.set_vertex_buffer(0, self.particle_buffers[(self.frame_num + 1) % 2].slice(..)); // particle_bufferを渡している
    render_pass.set_vertex_buffer(1, self.vertices_buffer.slice(..));
    render_pass.draw(0..24, 0..NUM_PARTICLES);
}
self.queue.submit(Some(command_encoder.finish()));

いろいろ書いていますが、ポイントはrender_pass.set_vertex_buffer(0, self.particle_buffers[(self.frame_num + 1) % 2].slice(..));の行です。self.frame_numに応じて、二つあるparticle_bufferが交互に描画されるという処理になっています。

この実装により「計算元のバッファ」と「計算結果を格納するバッファ」を実現されるのですが、イメージしづらいかもしれません。文章だけでは説明しづらいので図を用意しました。

まずBindGroupはそれぞれ以下のようになっています。初期化したタイミングではいずれのバッファにも同じデータが入っています。

はじめにBindGroup①を使って更新と描画を行います。BufferAのパーティクルの位置を元に１フレーム分だけ進んだ位置を計算し、BufferBに格納し描画します。

この時、Bindgroup②はこのようになっています。BufferBが更新されてBufferB'になっていますね。

次のフレームではBindGroup②を使って計算します。BufferB'を元に計算した結果がBufferAに格納されBufferA'となり、それが描画されます。

あとは、これを繰り返していくことで、順次パーティクルが更新されるというわけです。

シェーダーは引用するには少々長いのですが、以下のとおりです。

struct Particle {
  pos : vec2<f32>,
  vel : vec2<f32>,
  acc : vec2<f32>,
};

struct SimParams {
    q_charge: f32,
    blank_level: f32,
    time_delta: f32,
    d_max: f32,
    sustain: f32,
    dt: f32,
    dt_2: f32,
    dtt: f32,
    v_max: f32,
    fc: f32,
    pi: f32
};

@group(0) @binding(0) var<uniform> params : SimParams;

@group(1) @binding(0) var<storage, read> particlesSrc : array<Particle>;
@group(1) @binding(1) var<storage, read_write> particlesDst : array<Particle>;

@group(2) @binding(0) var t: texture_2d<f32>;

fn rotate(v: f32) -> f32 {
    if (v > 1.0) {
        return v - 2.0;
    } else if (v < -1.0) {
        return v + 2.0;
    } else {
        return v;
    }
}

fn get_texel(texture: texture_2d<f32>, x: f32, y: f32) -> vec4<f32> {
    let texture_size = textureDimensions(texture);
    let width = texture_size[0];
    let height = texture_size[1];
    let ix = i32(floor((x + 1.0) * f32(width) / 2.0));
    let iy = i32(floor((1.0 - y) * f32(height) / 2.0));
    return textureLoad(texture, vec2<i32>(ix, iy), 0);
}

@stage(compute)
@workgroup_size(64)
fn main(@builtin(global_invocation_id) global_invocation_id: vec3<u32>) {
  let total = arrayLength(&particlesSrc);
  let index = global_invocation_id.x;
  if (index >= total) {
    return;
  }

  var vPos : vec2<f32> = particlesSrc[index].pos;
  var vVel : vec2<f32> = particlesSrc[index].vel;
  var vAcc : vec2<f32> = particlesSrc[index].acc;

  let Epsilon = 0.0000001;
  let texel = get_texel(t, vPos.x, vPos.y); // 画像から受ける力を取り出す
  var acc : vec2<f32> = vec2<f32>(texel[0], texel[1]);

  // パーティクル同士の相互作用を計算する
  var i: u32 = 0u;
  loop {
      if (i >= total) {
          break;
      }
      if (i == index) {
          continue;
      }
      let dotPos = particlesSrc[i].pos;
      var dp : vec2<f32> = dotPos - vPos;
      dp.x = rotate(dp.x);
      dp.y = rotate(dp.y);

      let d2 = dot(dp, dp) + 0.00003;
      let d = sqrt(d2);
      if (abs(d) > Epsilon) {
          let q = params.q_charge / d2;
          acc -= q * dp / d;
      }
      continuing {
          i = i + 1u;
      }
  }
  acc *= params.q_charge;

  var vel : vec2<f32> = (params.sustain * vVel) + params.dt_2 * (vAcc + acc);
  let vlen = length(vel);
  if (vlen > params.v_max) {
      vel *= params.v_max / vlen;
  }
  var dp : vec2<f32> = (vel * params.dt) + (0.5 * acc * params.dtt);
  let dlen = length(dp);
  if (dlen > params.d_max) {
      dp *= params.d_max / dlen;
  }
  var pos : vec2<f32> = vPos + dp;

  pos.x = rotate(pos.x);
  pos.y = rotate(pos.y);

  particlesDst[index] = Particle(pos, vel, acc);
}