アフィン変換とその問題
アフィン変換とは
アフィン変換は、回転、移動、変形、拡大、縮小、移動を一発で変換するアルゴリズムである。
詳しい説明はWIkipedeiaでも読めば良い。
JpegのIDCTに続いて行列が出てくる。要するに画像処理の世界は離散数学と行列処理ばかりの世界なのだよ。この演算は、三次元に拡張することも可能である。ちなみに三次元に拡張した場合、四次元行列を使うか、四元数を使うかはケースバイケースらしいぞ（3Dは専門外だ）
ところで2Dでは虚数を使う選択肢はあまりない。コンピュータが遅かりし時代から存在した画像処理の世界は虚数どころか小数点演算すら避けたい世界だからだろう。
ところでアフィン変換のrustインプリである（長いので最後に置いた）。アフィン変換自体は非常に簡単なのである。
アフィン変換の問題
しかし、拡大や回転を行うと以下の問題がおきる
そう変換前の座標に無い点が抜け落ちるのである。これは果たしてどうすれば良いのか？そこで出てくるのが、拡大縮小でおなじみの補完アルゴリズムのである。
拡大縮小アルゴリズム
良く知られて居るアルゴリズムとしては、
ニアレストネイバー法
バイリニア法
バイキュービック法
Lanczos3
などである。ニアレストネイバーは正直汚いので整数倍以外はバイリニアかバイキュービックの二択になろうかと。Lanczos3は少し重い。昔DCTを応用した方法があると聞いた記憶があるが特許が成立しているようだ。
そう、アフィン変換より、こいつのインプリが面倒なのだ。
今までアフィン変換のテストをするためにJPEGの実装していたのだ。残る相手はWebPとdeflateだな（何と戦っているのだ）。GIFとTIFFはもう勘弁してください。IDCTだけではなくbit readerも最適化したので、ややでかい画像でも2s切る様になったかな。デバッグ用コードを削るともう少し速くなるかも、それ以外に後はunsafeなコードを入れる方法があるけど辞めておく。fetchの途中のデータを奪えればもう少し速くなる気がする。やはり、IDCTがネックなのでGPUやSMIDがフルに使えないとネイティブ実装には勝てない。
アフィン変換コード
/*
* affine.rs Mith@mmk (C) 2022
* create 2022/03/13
*/
use core::f32::consts::PI;
use super::canvas::*;
pub struct Affine {
affine: [[f32;3];3], // 3*3
}
impl Affine {
pub fn new() -> Self {
let affine = [
[1.0,0.0,0.0],
[0.0,1.0,0.0],
[0.0,0.0,1.0]];
Self {
affine,
}
}
fn matrix(self:&mut Self,f: &[[f32;3];3]) {
let affine = self.affine;
let mut result:[[f32;3];3] = [[0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]];
for i in 0..3 {
for j in 0..3 {
result[i][j] = affine[i][0] * f[0][j]
+ affine[i][1] * f[1][j]
+ affine[i][2] * f[2][j];
}
}
self.affine = result;
}
pub fn translation(self:&mut Self,x:f32,y:f32) {
self.matrix(&[[1.0 ,0.0 , x],
[0.0 ,1.0 , y],
[0.0 ,0.0 ,1.0]]);
}
pub fn invert_x(self:&mut Self) {
self.matrix(&[[-1.0 , 0.0 , 0.0],
[ 0.0 , 1.0 , 0.0],
[ 0.0 , 0.0 , 1.0]]);
}
pub fn invert_y(self:&mut Self) {
self.matrix(&[[ 1.0 , 0.0 , 0.0],
[ 0.0 , -1.0 , 0.0],
[ 0.0 , 0.0 , 1.0]]);
}
pub fn invert_xy(self:&mut Self) {
self.matrix(&[[-1.0 , 0.0 , 0.0],
[ 0.0 ,-1.0 , 0.0],
[ 0.0 , 0.0 , 1.0]]);
}
pub fn scale(self:&mut Self,x:f32,y:f32) {
self.matrix(&[[ x ,0.0 , 0.0],
[ 0.0 , y , 0.0],
[ 0.0 ,0.0 , 1.0]]);
}
pub fn rotate_by_dgree(self:&mut Self,theta:f32) {
let theta = PI * theta / 180.0;
self.rotate(theta);
}
pub fn rotate(self:&mut Self,theta:f32) {
let c = theta.cos();
let s = theta.sin();
self.matrix(&[[ c , -s , 0.0],
[ s , c , 0.0],
[ 0.0 ,0.0 , 1.0]]);
}
pub fn shear(self:&mut Self,x:f32,y:f32) {
self.matrix(&[[ 1.0 , y , 0.0],
[ x ,1.0 , 0.0],
[ 0.0 ,0.0 , 1.0]]);
}
/* not implement Scaling Routine */
pub fn conversion(self:&mut Self,input_canvas: &Canvas,output_canvas:&mut Canvas) {
let min_x = 0;
let max_x = output_canvas.width() as i32;
let min_y = 0;
let max_y = output_canvas.height() as i32;
let ox = max_x / 2;
let oy = max_y / 2;
// |X| |a00 a01 a02||x| X = a00 * x + a01 * y + a02
// |Y| = |a10 a11 a12||y| Y = a10 * x + a11 * y + a12
// |Z| |a20 a21 a22||1| _ do not use
let x0 = self.affine[0][0];
let x1 = self.affine[0][1];
let x2 = self.affine[0][2];
let y0 = self.affine[1][0];
let y1 = self.affine[1][1];
let y2 = self.affine[1][2];
for y in 0..input_canvas.height() as usize {
let input_base_line = input_canvas.width() as usize * 4 * y;
for x in 0..input_canvas.width() as usize {
let offset = input_base_line + x * 4;
let x_ = x as i32 - ox;
let y_ = y as i32 - oy;
let xx_ = x_ as f32 * x0 + y_ as f32 * x1 + x2;
let xx = (xx_ + ox as f32).round() as i32;
if xx < min_x || xx >= max_x {continue} // Out of bound
let yy_ = x_ as f32 * y0 + y_ as f32 * y1 + y2;
let yy = (yy_ + oy as f32).round() as i32;
if yy < min_y || yy >= max_y {continue} // Out of bound
let output_offset = ((yy * max_x + xx) *4) as usize;
output_canvas.buffer[output_offset ] = input_canvas.buffer[offset ];
output_canvas.buffer[output_offset + 1] = input_canvas.buffer[offset + 1];
output_canvas.buffer[output_offset + 2] = input_canvas.buffer[offset + 2];
output_canvas.buffer[output_offset + 3] = input_canvas.buffer[offset + 3];
}
}
}
}