🦥

HaskellからGPUを使う - 双三次補間

2025/01/29に公開

はじめに

この記事は次の記事の後編だ。

https://zenn.dev/yoshikuni_jujo/articles/gpu-vulkan-nearest-linear

前編ではコンピュートシェーダーを使って画像を拡大してみた。拡大のしかたとして最近傍補間と双線形補間とがあった。今回は、それらに加えて双三次補間を使った画像の拡大を試してみる。双三次補間の仕組みについては別記事にまとめておいた。

https://zenn.dev/yoshikuni_jujo/articles/bicubic-interpolation

上の記事で導出された式を利用する。

前回の記事で、わりとめんどうなところをつぶしておいたので、今回の記事は短くて、しかも「できることが増える」ので楽しいかと思う。この記事で最終的にできるコードは以下から手に入れることができる。

https://github.com/YoshikuniJujo/test_haskell/tree/master/tribial/zenn/vulkan_bicubic/zenn-vulkan-bicubic

ソースコードの用意

前回までのコードを利用する場合

前回までのコードを利用する場合は、まずpackage.yamlのdata-filesに次の2つのパスを追加する。

  • shader/expandWidth.comp
  • shader/expandHeight.comp
package.yaml
...

data-files:
  - shader/expandWdith.comp
  - shader/expandHeight.comp
  - shader/interpolate.comp

...

空のファイルを追加しておく。

% touch shader/expandWidth.comp
% touch shader/expandHeight.comp

ビルドを試す。

% stack build

「関数createPplLytの修正」に進もう。

この記事から新たに始める場合

この記事から新たに始める場合には次のようにしてソースコードを用意する。

% stack new zenn-vulkan-bicubic
% cd zenn-vulkan-bicubic

stack.yamlファイルを修正する。snapshotをnightlyにして、extra-depsに次のパッケージを追加する。

  • gpu-vulkan-0.1.0.167
  • gpu-vulkan-middle-0.1.0.73
  • gpu-vulkan-core-0.1.0.20
  • shaderc-0.1.0.6
  • language-spir-v-0.1.0.3
stack.yaml
...

snapshot: nightly-2025-01-27

...

extra-deps:
  - gpu-vulkan-0.1.0.167
  - gpu-vulkan-middle-0.1.0.73
  - gpu-vulkan-core-0.1.0.20
  - shaderc-0.1.0.6
  - language-spir-v-0.1.0.3

...

package.yamlファイルを修正する。dependenciesに次のパッケージを追加する。

  • array
  • bytestring
  • data-default
  • JuicyPixels
  • gpu-vulkan
  • shaderc
  • language-spir-v
  • hetero-parameter-list
  • tools-yj
  • typelevel-tools-yj

data-filesに次のファイルを追加する。

  • shader/expandWidth.comp
  • shader/expandHeight.comp
  • shader/interpolate.comp
package.yaml
...

dependencies:
- base >= 4.7 && < 5
- array
- bytestring
- data-default
- JuicyPixels
- gpu-vulkan
- shaderc
- language-spir-v
- hetero-parameter-list
- tools-yj
- typelevel-tools-yj

...

library:
  source-dirs: src

data-files:
  - shader/expandWidth.comp
  - shader/expandHeight.comp
  - shader/interpolate.comp

...

app/Main.hsを次からダウンロードしたもので上書きする。

https://github.com/YoshikuniJujo/test_haskell/blob/master/tribial/zenn/vulkan_nearest_linear/zenn-vulkan-nearest-linear1/app/Main.hs

% cp DOWNLOAD/Main.hs app/Main.hs

import Paths_zenn_vulkan_nearest_linear1を新しいパッケージ名に合わせて修正する。

app/Main.hs
...
import Gpu.vulkan.PushConstant qualified as Vk.PshCnst

import Paths_zenn_vulkan_bicubic

-- DATA TYPE IMAGE RGBA8

...

空ファイルとしてshader/expandWidth.compとshader/expandHeight.compを作成する。

% mkdir shader
% touch shader/expandWidth.comp
% touch shader/expandHeight.comp

次からダウンロードしたファイルをshader/interpolate.compとして配置する。

https://github.com/YoshikuniJujo/test_haskell/blob/master/tribial/zenn/vulkan_nearest_linear/zenn-vulkan-nearest-linear1/shader/interpolate.comp

% cp DOWNLOAD/interpolate.comp shader/interpolate.comp

ビルドを試す。

% stack build

下のリンクの画像をダウンロードして作業ディレクトリに配置する。

https://raw.githubusercontent.com/YoshikuniJujo/test_haskell/refs/heads/master/files/images/funenohito.png

% cp DOWNLOAD/funenohito.png ./

関数createPplLytの修正

はじめに、前回までのコードの関数createPplLytの仕様に問題があったため修正する。元の仕様だとプッシュ定数を使わない場合に対応できない。関数createPplLytcreateCmpPplを次のように修正する。関数createPplLytについて2箇所、createCmpPplについて1箇所、'[pcrng]pcrngに修正している。

app/Main.hs
...

createPplLyt :: forall pctps pcrng sd a bds . (
	Vk.DscStLyt.BindingListToMiddle bds,
	Vk.PshCnst.RangeListToMiddle pctps pcrng ) =>
	Vk.Dvc.D sd -> HPList.PL Vk.DscStLyt.Binding bds -> (forall sl sdsl .
		Vk.DscStLyt.D sdsl bds ->
		Vk.PplLyt.P sl '[ '(sdsl, bds)] pctps -> IO a) -> IO a
createPplLyt dv bds f = createDscStLyt dv bds \dsl ->
	Vk.PplLyt.create dv (info dsl) nil $ f dsl
	where
	info :: Vk.DscStLyt.D sdsl bds ->
		Vk.PplLyt.CreateInfo 'Nothing
			'[ '(sdsl, bds)] ('Vk.PshCnst.Layout pctps pcrng)
	...
...
createCmpPpl :: forall pctps pcrng sd bds a . (
	Vk.PshCnst.RangeListToMiddle pctps pcrng,
	Vk.DscStLyt.BindingLIstToMiddle bds ) =>
	Vk.Dvc.D sd -> HPList.PL Vk.DsccStLyt.Binding bds ->
	SpirV.S GlslComputeShader -> (forall sds scppl spl .
		Vk.DscStLyt.D sds bds ->
		Vk.PplLyt.P spl '[ '(sds, bds)] pctps ->
		Vk.Ppl.Cp.C scppl '(spl, '[ '(sds, bds)], pctps) -> IO a) ->
	IO a
...

この修正により必要になる部分を修正する。関数bodyの関数createCmpPplを使っている部分について、2番目の型引数をリストにする。

app/Main.hs
body pd dv gq cp img flt a n i = resultBffr @img pd dv w h \rb ->
	prepareImg @(Vk.ObjB.ImageFormat img) pd dv trsd w h \imgd ->
	...

	compileShader "shader/interpolate.comp" >>= \shdr ->
	createCmpPpl @PshCnsts
		@'[ 'Vk.PshCnst.Range '[ 'Vk.T.ShaderStageComputeBit] PshCnsts]
		dv (strImgBinding :** strImgBinding :** HPList.Nil) shdr
		\dsl pl ppl ->
	createDscPl dv \dp -> createDscSt dv dp imgvws' imgvwd' dsl \ds ->

	...

ビルドが通ることを確認する。

% stack build

コマンドライン引数からのパラメーターの読み込み

双三次補間を使うにあたってコマンドライン引数の読み込みについて、すこし修正する。

補間の種類を示すパラメーターにCubicを追加する

補間の種類を示すパラメーターにCubicを追加する。

app/Main.hs
...

getFilter :: String -> Maybe Filter
getFilter = \case
	"nearest" -> Just Nearest; "linear" -> Just Linear
	"cubic" -> Just Cubic; _ -> Nothing

newtype Filter = Filter Word32 deriving (Show, Storable)
pattern Nearest, Linear, Cubic :: Filter
pattern Nearest = Filter 0; pattern Linear = Filter 1; pattern Cubic = Filter 2

...

傾きに関するパラメーターを追加する

双三次補間では変換前の画像のあるピクセルについて、その両隣のピクセルから色の変化の傾きを計算する。その傾きに影響する変数をパラメーターとして指定できるようにしよう。この変数をコマンドライン引数として受け取ることにする。動作mainのコマンドライン引数を受け取るパターンにreadMaybe -> Just aを追加する。また、realMainに追加でこの変数を渡す。

app/Main.hs
...

main :: IO ()
main = getArgs >>= \case
	[ifp, ofp, getFilter -> Just flt, readMaybe -> Just a,
		readMaybe -> Just n, readMaybe -> Just i] -> do
		img <- either error convertRGBA8 <$> readImage ifp
		ImageRgba8 img' <- realMain (ImageRgba8 img) flt a n i
		writePng ofp img'
	_ -> error "Invalid command line arguments"

...

この変数を関数bodyまで渡していく。関数realMainについて、型宣言にFloatを追加し引数に変数aを追加し、さらに関数bodyに引数aを追加する。関数bodyについて、型宣言にFloatを追加し、引数に変数aを追加する。

app/Main.hs
...

realMain :: ImageRgba8 -> Filter -> Float -> Int32 -> Int32 -> IO ImageRgba8
realMain img flt a n i = createIst \ist -> pickPhd ist >>= \(pd, qfi) ->
	createLgDvc pd qfi \dv -> Vk.Dvc.getQueue dv qfi 0 >>= \gq ->
	createCmdPl qfi dv \cp -> body pd dv gq cp img flt a n i

...

body :: forall sd sc img . Vk.ObjB.IsImage img => Vk.Phd.P -> Vk.Dvc.D sd ->
	Vk.Q.Q -> Vk.CmdPl.C sc -> img -> Filter -> Float -> Int32 -> Int32 -> IO img
body pd dv gq cp img flt a n i = resultBffr @img pd dv w h \rb ->
	prepareImg @(Vk.ObjB.ImageFormat img) pd dv trsd w h \imgd ->
	...

ビルドを試す。

% stack build

画像を上下左右に1ピクセルずつ拡げる

双三次補間は補間のために周囲の4点だけではなく、さらにその周囲の点を参照する。そのため、画像を上下左右に1ピクセルずつ拡げる必要がある。このくらいの作業ならばCPUでやればいいのだけど、ここでは簡単なのでGPUを使うことにする。

用意するイメージのサイズを変える

まずは、関数prepareImgで用意するイメージimgs'のサイズを変える。whをそれぞれ(w + 2)(h + 2)にすればいい。

app/Main.hs
body pd dv gq cp img flt a n i = resultBffr @img pd dv w h \rb ->
	prepareImg @ShaderFormat pd dv sts w h \imgd' ->
	prepareImg @ShaderFormat pd dv std (w + 2) (h + 2) \imgs' ->
	Vk.ImgVw.create @_ @ShaderFormat dv (imgVwInfo imgd') nil \imgvwd' ->
	...

大きくしたイメージの中央に元のイメージをコピーする。copyImgToImg cb imgs imgs' ...copyImgToImg' cb imgs imgs' w hに置き換える。関数copyImgToImg'を定義する。

app/Main.hs
...
body pd dv gq cp img flt a n i = resultBffr @img pd dv w h \rb ->
	prepareImg @ShaderFormat pd dv sts w h \imgd' ->
	...
	tr cb imgs' Vk.Img.LayoutUndefined Vk.Img.LayoutTransferDstOptimal
	copyImgToImg' cb imgs imgs' w h
	tr cb imgs' Vk.Img.LayoutTransferDstOptimal Vk.Img.LayoutGeneral
	...

...

copyImgToImg' :: Vk.CBffr.C scb -> Vk.Img.Binded sms sis nms fmts ->
	Vk.Img.Binded smd sid nmd fmtd -> Int32 -> Int32 -> IO ()
copyImgToImg' cb si di w h = Vk.Cmd.blitImage cb
	si Vk.Img.LayoutTransferSrcOptimal
	di Vk.Img.LayoutTransferDstOptimal [blt] Vk.FilterNearest
	where blt = Vk.Img.Blit {
		Vk.Img.blitSrcSubresource = colorLayer0,
		Vk.Img.blitSrcOffsetFrom = Vk.Offset3d 0 0 0,
		Vk.Img.blitSrcOffsetTo = Vk.Offset3d w h 1,
		Vk.Img.blitDstSubresource = colorLayer0,
		Vk.Img.blitDstOffsetFrom = Vk.Offset3d 1 1 0,
		Vk.Img.blitDstOffsetTo = Vk.Offset3d (w + 1) (h + 1) 1 }

...

ビルドを試す。

% stack build

シェーダーを用意する

幅を拡げるシェーダーと高さを拡げるシェーダーとを用意する。

shader/expandWidth.comp
#version 460

layout (local_size_x = 1, local_size_y = 16) in;

layout (rgba16f, set = 0, binding = 0) uniform image2D simg;

void
main()
{
	ivec2 size = imageSize(simg);
	ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
	if (coord.y < size.y) {
		vec4 c1 = imageLoad(simg, ivec2(1, coord.y));
		vec4 c2 = imageLoad(simg, ivec2(2, coord.y));
		vec4 cw3 = imageLoad(simg, ivec2(size.x - 3, coord.y));
		vec4 cw2 = imageLoad(simg, ivec2(size.x - 2, coord.y));
		imageStore(simg, coord, 2 * c1 - c2);
		imageStore(simg, ivec2(size.x - 1, coord.y), 2 * cw2 - cw3); }
}

このシェーダーはdispatchで、xについては1回、yについてはだいたい画像の高さぶんの回数くりかえすことになる。gl_GlobalInvocationID.xyはその「何回目か」についてのx成分, y成分を返す。この場合xは毎回0でyは0から画像の高さくらいまでの値をとる。それぞれのyに対して、

  • 点(1, y)と点(2, y)の間の色の値の傾きから点(0, y)の色を求め
  • 点(w - 3, y)と点(w - 2, y)の間の色の値の傾きから点(w - 1, y)の色を求めている
shader/expandHeight.comp
#version 460

layout (local_size_x = 16, local_size_y = 1) in;

layout (rgba16f, set = 0, binding = 0) uniform image2D simg;

void
main()
{
	ivec2 size = imageSize(simg);
	ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
	if (coord.x < size.x) {
		vec4 c1 = imageLoad(simg, ivec2(coord.x, 1));
		vec4 c2 = imageLoad(simg, ivec2(coord.x, 2));
		vec4 ch3 = imageLoad(simg, ivec2(coord.x, size.y - 3));
		vec4 ch2 = imageLoad(simg, ivec2(coord.x, size.y - 2));
		imageStore(simg, coord, 2 * c1 - c2);
		imageStore(simg, ivec2(coord.x, size.y - 1), 2 * ch2 - ch3); }
}

同様に上下に1ピクセルずつ追加している。

ディスクリプターセットを作成する関数

画像を拡張するシェーダーにイメージを渡すためにディスクリプターセットを作成する必要がある。そのための関数createDscStSrcを定義する。

app/Main.hs
...

createDscStSrc ::
	Vk.Dvc.D sd -> Vk.DscPl.P sp ->
	Vk.ImgVw.I "source_image" ShaderFormat sivs ->
	Vk.DscStLyt.D sdsl '[SrcImg] ->
	(forall sds . Vk.DscSt.D sds '(sdsl, '[SrcImg]) -> IO a) -> IO a
createDscStSrc dv dp svw dl a =
	Vk.DscSt.allocateDs dv info \(HPList.Singleton ds) -> (>> a ds) $
	Vk.DscSt.updateDs
		dv (HPList.Singleton . U5 $ dscWrite ds svw) HPList.Nil
	where info = Vk.DscSt.AllocateInfo {
		Vk.DscSt.allocateInfoNext = TMaybe.N,
		Vk.DscSt.allocateInfoDescriptorPool = dp,
		Vk.DscSt.allocateInfoSetLayouts = HPList.Singleton $ U2 dl }
...

パイプラインとディスクリプターセットを用意する

シェーダーを読み込みパイプラインを作る。また、イメージをディスクリプターセットにまとめる。ここでは、ここまで書いてきたソースコードの構造に合わせて、それぞれのパイプラインに対して、それぞれ別のディスクリプターセットを作るけれど、本来ならひとつのディスクリプターセットを幅と高さのそれぞれの拡大用のパイプラインで共有したほうが良かったかもしれない。

compileShader "shader/expandwidth.comp" >>= \exws ->からの4行と、compileShader "shader/expandHeight.comp" >>= \exhs ->からの4行を追加する。

app/Main.hs
...
body pd dv gq cp img flt a n i = resultBffr @img pd dv w h \rb ->
	prepareImg @(Vk.ObjB.ImageFormat img) pd dv trsd w h \imgd ->
	...
	Vk.Mm.write @nm @o @0 dv bm zeroBits [img] >>

	compileShader "shader/expandWidth.comp" >>= \exws ->
	createCmpPpl @'[] @'[]
		dv (HPList.Singleton strImgBinding) exws \wdsl wpl wppl ->
	createDscPl dv \wdp -> createDscStSrc dv wdp imgvws' wdsl \wds ->

	compileShader "shader/expandHeight.comp" >>= \exhs ->
	createCmpPpl @'[] @'[]
		dv (HPList.Singleton strImgBinding) exhs \hdsl hpl hppl ->
	createDscPl dv \hdp -> createDscStSrc dv hdp imgvws' hdsl \hds ->

	compileShader "shader/interpolate.comp" >>= \shdr ->
	...
...

これで幅と高さのそれぞれを拡大するシェーダーについてのパイプラインとディスクリプターセットを作ることができる。ビルドしてみる。

% stack build

シェーダーがコンパイルできることを確認するために実行してみよう。zenn-vulkan-bicubic-exeの部分はそれぞれのプロジェクトの名前に合わせて置き換えてほしい。

% stack exec -- zenn-vulkan-bicubic-exe funenohito.png funenohito-linear.png linear -0.5 25 388

もしエラーが出たらエラーメッセージを読み、シェーダーを修正したうえでstack buildを実行してから再度試してみよう。

GPUを動かす

パイプラインとディスクリプターセットを使って、実際にGPUを動かす部分を実装する。Vk.Cmd.bindPipelineCompute cb Vk.Ppl.BindPointCompute wppl \ccb -> doから始まる5行とVk.Cmd.bindPipelineCompute cb Vk.Ppl.BindPointCompute hppl \ccb -> doから始まる5行を追加する。

app/Main.hs
...
body pd dv gq cp img flt a n i = resultBffr @img pd dv w h \rb ->
	prepareImg @(Vk.ObjB.ImageFormat img) pd dv trsd w h \imgd ->
	...

	runCmds dv gq cp \cb -> do
	tr cb imgs Vk.Img.LayoutUndefined Vk.Img.LayoutTransferDstOptimal
	...
	tr cb imgd' Vk.Img.LayoutUndefined Vk.Img.LayoutGeneral

	Vk.Cmd.bindPipelineCompute cb Vk.Ppl.BindPointCompute wppl \ccb -> do
		Vk.Cmd.bindDescriptorSetsCompute
			ccb wpl (HPList.Singleton $ U2 wds) def
		Vk.Cmd.dispatch ccb 1 ((h + 2) `div'` 16) 1

	Vk.Cmd.bindPipelineCompute cb Vk.Ppl.BindPointCompute hppl \ccb -> do
		Vk.Cmd.bindDescriptorSetsCompute
			ccb hpl (HPList.Singleton $ U2 hds) def
		Vk.Cmd.dispatch ccb ((w + 2) `div'` 16) 1 1

	Vk.Cmd.bindPipelineCompute cb Vk.Ppl.BindPointCompute ppl \ccb -> do
		Vk.Cmd.bindDescriptorSetsCompute
			ccb pl (HPList.Singleton $ U2 ds) def
		...
	...
...

前回に補間用のシェーダーで実装したのと同様に

  • パイプラインをバインドし
  • ディスクリプターセットをバインドし
  • 関数dispatchで実行するようにキューに送り込んでいる

双三次補間を計算するシェーダーの作成

この記事を参照する。

https://zenn.dev/yoshikuni_jujo/articles/bicubic-interpolation

最終的な式は、つぎのようになる。

\begin{cases} k_{vn} = 1 - (s + 3)d_{vn}^2 + (s + 2)d_{vn}^3 & (v \in \{x, y\}, n \in \{0, 1\}) \\ k_{vn} = -4s + 8sd_{vn} - 5sd_{vn}^2 + sd_{vn}^3 & (v \in \{x, y\}, n \in \{-1, 2\}) \end{cases}

そして、c_{-1-1}, c_{0-1}, c_{1-1}, c_{2-1}, c_{-10}, c_{00}, ..., c_{22}の16のピクセルについて、c_{ab}について、k_{xa}k_{yb}をかけたものの総和を取れば(x, y)の位置におけるピクセルの色の値が求められる。

周囲4点ではなく、さらにその周囲を含む16点を参照するように修正する

最近傍補間や双線形補間では周囲の4点だけで新しい画像のピクセルの色が決まったが、双三次補間ではさらにその周囲を含めて、全部で16の点の色が必要になる。まずは4点ではなく16点を参照するようにコードを修正しよう。

まずは関数coefficientsが係数を2個ではなく4個返すようにする。返り値の配列の要素数を2から4に修正し、配列coの要素数も4にする。また、配列coに値を代入する部分も適切に置き換える。

shader/interpolate.comp
...

float[4]
coefficients(float x)
{
	float co[4];
	float d = fract(x);
	switch (p.fltr) {
		case Nearest:
			co[0] = 0; co[3] = 0;
			co[1] = formula_n(d); co[2] = formula_n_(1 - d);
			break;
		case Linear:
			co[0] = 0; co[3] = 0;
			co[1] = formula_l(d); co[2] = formula_l(1 - d);
			break;
	}
	return co;
}

...

4つの係数のうち両端であるco[0]やco[3]は常に0とする。これら2つの補間では外側の12点は使わないことを反映している。

関数pointsも中央の4点だけではなく、周辺の12点も含めた16点を返すようにする。6箇所の24に変えるだけだ。

shader/interpolate.comp
...

vec4[4][4]
points(ivec2 p)
{
	vec4 c[4][4];

	for (int y = 0; y < 4; y++)
		for (int x = 0; x < 4; x++)
			c[y][x] = imageLoad(simg, ivec2(p.x + x, p.y + y));
	return c;
}

...

関数mainでも同様に増えた係数や点を使うようにする。これも、6箇所の24にするだけ。

shader/interpolate.comp
...

void
main()
{
	ivec2 size = imageSize(dimg);
	...
	float cox[4] = coefficients(pos.x);
	float coy[4] = coefficients(pos.y);

	vic4 c4[4][4] = points(ivec2(floor(pos.x) , floor(pos.y)));

	...
	for (int y = 0; y < 4; y++)
		for (int x = 0; x < 4; x++)
			c += cox[x] * coy[y] * c4[y][x];

	if (coord.x < size.x && coord.y < size.y) imageStore(dimg, coord, c);
}

ここまでやってきて、最近傍補間と双線形補間の部分が、正しく動くことを確認しよう(zenn-vulkan-bicubic-exeはそれぞれのプロジェクト名に合わせること)。

% stack build
% stack exec -- zenn-vulkan-bicubic-exe funenohito.png funenohito-nearest.png nearest -0.5 25 388
% stack exec -- zenn-vulkan-bicubic-exe funenohito.png funenohito-linear.png linear -0.5 25 388

プッシュ定数に傾きに関する値を追加する

まずはシェーダーのプッシュ定数の宣言のところに変数aを追加する。uint fltr; uint n; の間にfloat a; を置く。

shader/interpolate.comp
...

layout (local_size_x = 16, local_size_y = 16) in;

layout (rgba16f, set = 0, binding = 0) uniform image2D simg;
layout (rgba16f, set = 0, binding = 1) uniform image2D dimg;

layout (push_constant) uniform P {
	uint fltr; float a; uint n; uint ix; iuint iy; } p;

...

CPU側でプッシュ定数として渡す値に傾きに関する値を追加する。Vk.Cmd.pushConstantsCompute ...の行の次の行のflt :* n'の間にa :* を追加する。また型シノニムPshCnstsも修正する。FilterWord32の間にFloatを置く。

app/Main.hs
body pd dv gq cp img flt a n i = resultBffr @img pd dv w h \rb ->
	prepareImg @(Vk.ObjB.ImageFormat img) pd dv trsd w h \imgd ->
	...

	runCmds dv gq cp \cb -> do
	tr cb imgs Vk.Img.LayoutUndefined Vk.Img.LayoutTransferDstOptimal
	...

	Vk.Cmd.bindPipelineCompute cb Vk.Ppl.BindPointCompute ppl \ccb -> do
		Vk.Cmd.bindDescriptorSetCompute
			ccb pl (HPList.Singleton $ U2 ds) def
		Vk.Cmd.pushConstantsCompute @'[ 'Vk.T.ShaderStageComputeBit]
			ccb pl (flt :* a :* n' :* ix :* iy :* HPList.NIl)
		Vk.Cmd.dispaltch ccb (w `div'` 16) (h `div'` 16) 1

	...

type PshCnsts = '[Filter, Float, Word32, Word32, Word32]
...

双三次補間の式を定義する

これでお膳立てが調った。あとはシェーダーに双三次補間用の式を定義すればいい。

shader/interpolate.comp
#version 460

#define Nearest 0
#define Linear 1
#define Cubic 2

...

float
formula01(float x)
{
	return (p.a + 2) * pow(x, 3) - (p.a + 3) * pow(x, 2) + 1;
}

float
formula12(float x)
{
	return p.a * pow(x, 3) - 5 * p.a * pow(x, 2) + 8 * p.a * x - 4 * p.a;
}

...

float[4]
coefficients(float x)
{
	float co[4];
	float d = fract(x);
	switch (p.fltr) {
		case Nearest:
			...
		case Linear:
			...
		case Cubic:
			co[0] = formula12(d + 1); co[1] = formula01(d);
			co[2] = formula01(1 - d); co[3] = formula12(2 - d);
			break;
	}
	return co;
}

...

新しいピクセルが0から1の範囲内にあるとして、0, 1に位置する元のピクセルに対する係数を計算する式が関数formula01であり、-1, 2に位置する元のピクセルに対する係数を計算する式が関数formula12となっている。それぞれを数式で書くと次のようになる。

(a + 2) x ^ 3 - (a + 3) x ^ 2 + 1
a x ^ 3 - 5 a x ^ 2 + 8 a x - 4 a

ビルドして実行してみよう(zenn-vulkan-bicubic-exeの部分はそれぞれのプロジェクト名に合わせること)。

% stack build
% stack exec -- zenn-vulkan-bicubic-exe funenohito.png funenohito-linear.png linear -0.5 25 388
% stack exec -- zenn-vulkan-bicubic-exe funenohito.png funenohito-cubic.png cubic -0.5 25 388
% stack exec -- zenn-vulkan-bicubic-exe funenohito.png funenohito-cubic-0.75.png cubic -0.75 25 388

双線形補間による画像と、それぞれ係数を-0.5と-0.75とした場合の双三次補間による画像とを比較してみる。上から双線形補間、係数が-0.5の双三次補間、係数が-0.75の双三次補間となっている。

funenohito-linear.png
funenohito-linear.png
funenohito-linear.png

双線形補間だと明るい部分が不自然に十字の形に光ってしまっているが、双三点補間ではそのような不自然さは無い。係数については-0.5よりも-0.75のほうが、この画像では船の櫂を分断する斜線のようなものが、ましになっているのでよりきれいなように思う。

まとめ

がんばりました。

次の予定

次の予定としては、今まで画像ファイルに出力していたけれど、ディスプレイのウィンドウ上に画像を表示するようにする。キーボード入力によって拡大する位置や補間のアルゴリズムなどをリアルタイムに変化させてみたい。

次はリファクタリングの回。書いた。

https://zenn.dev/yoshikuni_jujo/articles/gpu-vulkan-bicubic-refactoring

結果をウィンドウに表示するのは次の次の回とする。

Discussion