デコンボリューション
2016/04/25 作成
2017/04/27 プログラム修正 配列メモリーの使用量を少なくしました。
2020/05/09 変換値の平均値の計算方法を正規の方法に修正しました。（スペクトル表示が無い為）
デコンボリューションとは、ピントのあっている部分に対して、ピントのあってない部分のボケ画像がかぶって、ボケが発生しているのを、光学系の特性値により、逆計算をしてかぶったボケ画像を取り除く手法のようですが、一般の写真画像は８ビットなので、演算上精度が不足するようです。
精度として１６ビット位必要らしいのですが、８ビットのデーターでも、ある程度効果がありそうなので、プログラムを作成してみました。
此処でのプログラムは、レンズのボケ特性を、ガウス分布で近似して、デコンボリューションをおこないます。
画像を、ガウス分布で、ボケ画像を作成したデーターにし、そのデーターをそのまま、同じガウス分布で、デコンボリューションを行うと、殆ど画質を低下させずに、元へ戻せますが、一旦８ビット画像にしてから、デコンボリューションを行うと、完全には元に戻せません。
実際にカメラで撮影したものは、デジタル化する時に、Jpeg圧縮によりデーターの欠落が生じるのと、レンズのボケ関数が不明なので、実際にボケが生じた画像のボケをとることは、かなり困難です。
しかし、普通に撮影された、ピントのあっている画像であれば、画像のシャープさを改善することが出来ます。
現在のデジタルカメラの解像度は、かなり高く、パソコンで写真を見るにしても、プリントアウトするにしても、一般に使用するなら必要以上に解像度が高いので、デコンボリューションの必要性は無いでしょう。
デコンボリューションの計算はフーリエ変換、逆フーリエ変換を使用しますが、画像のサイズが、２N
になっていないと、高速フーリエ変換(FFT)が使用できません。
そこで、多少無駄が生じても、画像サイズを２N にして高速フーリエ変換(FFT)が使える様に、FFT計算のプログラムを修正してあります。
画像サイズが大きい場合はFFTデコンボリューションを使用してください。
デコンボリューションの手順
１．画像を実数部と虚数部にフーリエ変換します。
２．デコンボリューション用フィルター値(マスク)を実数部と虚数部にフーリエ変換。
３．画像部のフーリエ値の実数部と虚数部をフィルターのフーリエ値実数部と虚数部で除算。
４．逆フーリエを実行して画像データーに戻す。
上記 ３．の除算をする時には、０で除算をしないようにするのは当然ですが、あまり小さい値で除算をしないようにします。
閾値を設定し、状況に応じて、値を変更します。
除算の方法ですが、通常のフィルターの計算のまま、フィルターの実数部と、虚数部を逆数にして処理をしても、フィルターの実数部と虚数部をそれぞれ自乗し加算して、逆数にして処理をしても、ほとんど同じ結果がえられます。
フィルターの虚数部は、非常に小さな値で、閾値の値より小さくなるので、結果にほとんど影響しません。
上画像、左側が、元の画像で、右側がデコンボリューションを行った画像です。
よくわからないので、胸元の刺繍を拡大してみました。
刺繍が鮮明になっています。
次にガウスフィルター(マスク)でボカしたものを、同じ値で、デコンボリューションした結果です。
結構元の状態に戻っています。
試しに、もっと ぼかした場合もテストしてみましたが、上画像の例ぐらいが限界で、元の状態まで戻すのは無理なようです。
デコンボリューション 1
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
rdr := conv_real[x][y] * conv_real[x][y] + conv_image[x][y] * conv_image[x][y];
if rdr > threshold then rdr := 1 / rdr
else rdr := 0;
// デコンボリューションの生成
r := (arBlu[x][y] - aiBlu[x][y]) * rdr;
i := (arBlu[x][y] + aiBlu[x][y]) * rdr;
arBlu[x][y] := r;
aiBlu[x][y] := i;
end;
デコンボリューション 2
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
rdr := sqrt(conv_real[x][y] * conv_real[x][y] + conv_image[x][y] * conv_image[x][y]);
if rdr > threshold then rdr := 1 / rdr
else rdr := 0;
// デコンボリューションの生成
r := (arBlu[x][y] - aiBlu[x][y]) * rdr;
i := (arBlu[x][y] + aiBlu[x][y]) * rdr;
arBlu[x][y] := r;
aiBlu[x][y] := i;
end;
デコンボリューション 3
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
rdr := abs(conv_real[x][y]);
rdi := abs(conv_image[x][y]);
if rdr > threshold then rdr := 1 / conv_real[x][y]
else rdr := 0;
if rdi > threshold then rdi := 1 / conv_image[x][y]
else rdi := 0;
// デコンボリューションの生成
r := arBlu[x][y] * rdr - aiBlu[x][y] * rdi;
i := arBlu[x][y] * rdi + aiBlu[x][y] * rdr;
arBlu[x][y] := r;
aiBlu[x][y] := i;
end;
デコンボリューション 4
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
rdr := absconv_real[x][y]);
if rdr > threshold then rdr := 1 / conv_real[x][y]
else rdr := 0;
// デコンボリューションの生成
arBlu[x][y] := arBlu[x][y] * rdr;
aiBlu[x][y] := aiBlu[x][y] * rdr;
end;
デコンボリューション 5
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
rdr := conv_real[x][y] * conv_real[x][y] + conv_image[x][y] * conv_image[x][y];
if rdr > threshold then rdr := 1 / rdr
else rdr := 0;
// デコンボリューションの生成
r := (arBlu[x][y] * conv_real[x][y] + aiBlu[x][y] * conv_image[x][y]) * rdr;
i := (aiBlu[x][y] * conv_real[x][y] - arBlu[x][y] * conv_image[x][y]) * rdr;
arBlu[x][y] := r;
aiBlu[x][y] := i;
end
conv_real conv_image がフィルター(マスク)のフーリエ変換値の実数部と虚数部です。
conve_image(虚数部)は、値が非常に小さいので計算から省略しても結果には影響しません。
デコンボリューション 1 は、実数部と虚数部をそれぞれ自乗してパワーレベルとして処理をしています。
２ は、１の値の平方根、リニアレベルとしたもので、３ は、完全にフィルター処理の逆の計算となってますが、どの方法で計算しても、デコンボリューション時のガウス関数のシグマ値と閾値を調整すると、殆ど同じ結果となります。
4 は虚数部の計算を省略した例です。
プログラム
プログラムはラプラシアン、ガウス フィルター(マスク)のプログラムに、デコンボリューションを追加しました。
ガウス計算のシグマ値を大きくする場合、デフォルトのガウス用配列の大きさ ９ では、有効範囲が不足するので、ガウス用配列を大きくする必要があります。
フィルターの拡張
フィルターの配列サイズを、画像のサイズと同じにして、フーリエ変換をする必要があります。
次の様に、四分割して、四隅にフィルターデーターを配置します。
unit DftDconvMain;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.ExtDlgs, Vcl.StdCtrls, Vcl.ExtCtrls, system.Math,
system.UITypes, Vcl.Grids;
type
D2array = array of array of Double;
B2array = array of array of Byte;
TDarray = array[0..0] of Double;
PDarray = ^TDarray;
TForm1 = class(TForm)
FileOpenBtn: TButton;
OpenPictureDialog1: TOpenPictureDialog;
MaskingBtn: TButton;
RadioGroup1: TRadioGroup;
FileSaveBtn: TButton;
SavePictureDialog1: TSavePictureDialog;
SigumaGrid: TStringGrid;
SigumaEdit: TLabeledEdit;
sourceBtn: TButton;
UnsharpEdit: TLabeledEdit;
RangeEdit: TLabeledEdit;
Timer1: TTimer;
RepaintBtn: TButton;
ScrollBox1: TScrollBox;
Image1: TImage;
magnificationEdit: TLabeledEdit;
GaussEdit: TLabeledEdit;
Memo1: TMemo;
FulBtn: TButton;
procedure FileOpenBtnClick(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure FormDestroy(Sender: TObject);
procedure MaskingBtnClick(Sender: TObject);
procedure FileSaveBtnClick(Sender: TObject);
procedure sourceBtnClick(Sender: TObject);
procedure Timer1Timer(Sender: TObject);
procedure RadioGroup1Click(Sender: TObject);
procedure SigumaEditChange(Sender: TObject);
procedure RangeEditChange(Sender: TObject);
procedure RepaintBtnClick(Sender: TObject);
procedure UnsharpEditChange(Sender: TObject);
procedure FulBtnClick(Sender: TObject);
private
{ Private 宣言 }
function SigumaMatSet: boolean;
procedure DFT(var real, image: D2array; inv, Y_SIZE, X_SIZE: integer);
procedure LaplacianMask;
procedure Imageout(Image: TBitmap);
procedure GaussianRoutine;
procedure trimming_imageout(Image: TBitmap);
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
const
OpenFileFilter =
'画像ファイル|*.png;*.jpg;*.gif;*.bmp;*.tif;*.ico;*.wdp'+
'|*.png|*.png' +
'|*.jpg|*.jpg' +
'|*.gif|*.gif' +
'|*.bmp|*.bmp' +
'|*.tif|*.tif' +
'|*.ico|*.ico' +
'|*.wdp|*.wdp';
SaveFileFilter =
'画像ファイル|*.png;*.jpg;*.gif;*.bmp;*.tif;*.wdp' +
'|*.png|*.png' +
'|*.jpg|*.jpg' +
'|*.gif|*.gif' +
'|*.bmp|*.bmp' +
'|*.tif|*.tif' +
'|*.wdp|*.wdp';
ImageHWC = 440; // 表示枠サイズ
gaus = 9; // ガウシアン配列サイズ
Laplacian = 0;
Gaussians = 1;
Deconvolution = 2;
Unsharp = 3;
var
InputDBitmap : TBitmap; // 入力データー表示用ビットマップ
OutputBitmap : TBitmap; // 回転画像表示用ビットマップ
GHeight : Integer; // 入力データー画像高さ
GWidth : Integer; // 入力データー画像幅
VRect : Trect; // 表示サイズ設定用
InFilename : string; // 入力ファイル名
threshold : Double; // 範囲指定
Gaussian : array[0..gaus - 1, 0..gaus - 1] of double;
const
kr = 0.299; // R 輝度変換係数
kg = 0.587; // G 輝度変換係数
kb = 0.114; // B 輝度変換係数
//------------------------
// 変倍出力
//------------------------
procedure TForm1.Imageout(Image: TBitmap);
var
Rect0 : Trect;
magnification : Double;
Check : Integer;
MW, MH : Integer;
begin
Val(magnificationEdit.Text,magnification,Check);
if Check <> 0 then begin
application.MessageBox('表示倍率入力に間違いがあります。','注意',0);
exit;
end;
MW := Round(GWidth * magnification);
MH := Round(GHeight * magnification);
Rect0 := Rect(0, 0, MW, MH);
Image1.Width := MW;
Image1.Height := MH;
Image1.Picture.Bitmap.SetSize(MW, MH);
Image1.Canvas.StretchDraw(Rect0, Image); // 出力枠に変倍出力
end;
//--------------------------------------------------
// 開いたファイルの画像表示
//--------------------------------------------------
procedure TForm1.trimming_imageout(Image: TBitmap);
begin
if Image.Width >= Image.Height then begin
Image1.Width := ImageHWC;
image1.Height := round((ImageHWC) * image.Height / Image.Width);
end
else begin
Image1.Height := ImageHWC;
image1.Width := round((ImageHWC) * image.Width / Image.Height);
end;
Image1.Picture.Bitmap.SetSize(Image1.Width, image1.Height);
Image1.Canvas.StretchDraw(Rect(0, 0, Image1.Width, image1.Height), Image); // 出力枠に変倍出力
magnificationEdit.Text := floatTostrF(Image1.Width / Image.Width, ffGeneral, 6, 4);
end;
//--------------------------------------------------------------------------------
// このルーチンは順変換（画像→フーリエ）、逆変換（フーリエ→画像）両方に使えます
// realには実数（順変換では元画像の値）、imageには虚数（はじめは0）が入っています
// invは順変換では1、逆変換では-1 マスク設定は2
// 変換後の値はreal、imageを書き換えて保存します
// 順変換時はシグマ値をデーター数で除算
// 逆変換時シグマ値をデーター数で除算しません。
// マスクはデーター数で除算しません
//--------------------------------------------------------------------------------
procedure TForm1.DFT(var real, image: D2array; inv, Y_SIZE, X_SIZE: integer);
var
M2 : Double;
Wnorm, Hnorm : Double;
RW, RH : Double;
Tre, Tim : D2array;
WcTable, WsTable : D2array;
HcTable, HsTable : D2array;
i, j, y, x, xx, yy: Integer;
re, im : double;
begin
// inv = 2 の設定を初期値にします
M2 := - 2 * PI;
Wnorm := 1;
Hnorm := 1;
// 逆変換設定
if inv = -1 then begin
M2 := 2 * PI;
end;
// 正変換設定
if inv = 1 then begin
Wnorm := 1 / X_SIZE;
Hnorm := 1 / Y_SIZE;
end;
// 作業用配列確保
SetLength(Tre, X_SIZE, Y_SIZE);
SetLength(Tim, X_SIZE, Y_SIZE);
SetLength(WcTable, X_SIZE, X_SIZE);
SetLength(WsTable, X_SIZE, X_SIZE);
SetLength(HcTable, Y_SIZE, Y_SIZE);
SetLength(HsTable, Y_SIZE, Y_SIZE);
// 変換テーブル生成 直流、低い周波数から 高いほうへ設定
RW := M2 / X_SIZE;
RH := M2 / Y_SIZE;
for i := 0 to X_SIZE - 1 do
for j := 0 to X_SIZE - 1 do begin
Wctable[i][j] := cos(i * j * RW);
Wstable[i][j] := sin(i * j * RW);
end;
for i := 0 to Y_SIZE - 1 do
for j := 0 to Y_SIZE - 1 do begin
Hctable[i][j] := cos(i * j * RH);
Hstable[i][j] := sin(i * j * RH);
end;
// X軸方向DFT
// y軸毎に変換
for y := 0 to Y_SIZE - 1 do
// 一次元DFT
for x := 0 to X_SIZE - 1 do begin
Tre[x][y] := 0.0;
Tim[x][y] := 0.0;
if inv = -1 then begin
for xx := 0 to X_SIZE - 1 do begin
re := real[xx][y];
im := image[xx][y];
Tre[x][y] := Tre[x][y] + re * Wctable[x][xx] - im * Wstable[x][xx];
Tim[x][y] := Tim[x][y] + re * Wstable[x][xx] + im * Wctable[x][xx];
end;
end
else begin
for xx := 0 to X_SIZE - 1 do begin
re := real[xx][y];
Tre[x][y] := Tre[x][y] + re * Wctable[x][xx];
Tim[x][y] := Tim[x][y] + re * Wstable[x][xx];
end;
end;
if inv = 1 then begin
Tre[x][y] := Tre[x][y] * wnorm;
Tim[x][y] := Tim[x][y] * wnorm;
end;
end;
// Y軸方向DFT
// X軸毎に変換
for x := 0 to X_SIZE - 1 do
// 一次元DFT
for y := 0 to Y_SIZE - 1 do begin
real[x][y] := 0.0;
if inv = -1 then begin
for yy := 0 to Y_SIZE - 1 do
real[x][y] := real[x][y] + Tre[x][yy] * Hctable[y][yy] - Tim[x][yy] * Hstable[y][yy];
end
else begin
image[x][y] := 0.0;
for yy := 0 to Y_SIZE - 1 do begin
re := Tre[x][yy];
im := Tim[x][yy];
real[x][y] := real[x][y] + re * Hctable[y][yy] - im * Hstable[y][yy];
image[x][y] := image[x][y] + re * Hstable[y][yy] + im * Hctable[y][yy];
end;
if inv = 1 then
image[x][y] := image[x][y] * Hnorm;
end;
if inv = 1 then
real[x][y] := real[x][y] * Hnorm;
end;
end;
//-----------------------------
// ガウシアン用データーの計算
//-----------------------------
function TForm1.SigumaMatSet: boolean;
const
MatSize = gaus;
var
X, Y : integer;
KS, Z : Double;
CH : integer;
Total : Double;
MDIV2 : integer;
caution : Pchar;
// シグマ値からガウス分布値の計算
function Maskdatacalc(X, Y : Integer; Z: Double): double;
begin
Result := 1 / 2 / pi / Z / Z * Exp(-(X * X + Y * Y) / 2 / Z / Z); // マスクデーターの計算
end;
begin
Result := False;
Z := 1;
CH := 0;
caution := '注意';
if RadioGroup1.ItemIndex = Gaussians then begin
val(GaussEdit.Text, Z, CH);
caution := 'Gaussian';
end;
if RadioGroup1.ItemIndex = Deconvolution then begin
val(SigumaEdit.Text, Z, CH);
caution := 'Deconvolution';
end;
if RadioGroup1.ItemIndex = Unsharp then begin
val(UnsharpEdit.Text, Z, CH);
caution := 'Unsharp';
end;
if CH <> 0 then begin
application.MessageBox('σ 入力に間違いがあります。', caution, 0);
exit;
end;
if Z <= 0 then begin
application.MessageBox('σ ゼロ 及び ゼロ以下ではいけません。', caution, 0);
exit;
end;
if RadioGroup1.ItemIndex = Deconvolution then begin
val(RangeEdit.Text, threshold, CH);
caution := 'Deconvolution';
if CH <> 0 then begin
application.MessageBox('Deconv Range 入力に間違いがあります。', caution, 0);
exit;
end;
if threshold <= 0 then begin
application.MessageBox('Deconv Range ゼロ及び負数ではいけません。', caution, 0);
exit;
end;
if threshold > 1 then begin
application.MessageBox('Deconv Range 1を超えた場合画像はなくなります。', caution, 0);
exit;
end;
end;
// ストリンググリッドの設定
SigumaGrid.ClientWidth := (SigumaGrid.DefaultColWidth + 1) * (MatSize + 1) - 1;
SigumaGrid.ColCount := MatSize + 1;
SigumaGrid.ClientHeight := (SigumaGrid.DefaultRowHeight + 1) * (MatSize + 1) - 1;
SigumaGrid.RowCount := MatSize + 1;
// 固定セルにナンバーリング
MDIV2 := MatSize div 2;
for Y := - MDIV2 to MDIV2 do SigumaGrid.Cells[0 , Y + MDIV2 + 1] := intTostr(Y);
for X := - MDIV2 to MDIV2 do SigumaGrid.Cells[X + MDIV2 + 1, 0] := intTostr(X);
// σ値によるアンシャープマスクデーター計算と表示
Total := 0;
for Y := - MDIV2 to MDIV2 do
for X := - MDIV2 to MDIV2 do begin
KS := Maskdatacalc( X, Y, Z);
Gaussian[X + MDIV2, Y + MDIV2] := KS;
// 合計計算
Total := Total + KS;
end;
// Gaussianフィルターの値変換 合計値が１になるように変換します
for Y := - MDIV2 to MDIV2 do
for X := - MDIV2 to MDIV2 do begin
KS := Gaussian[X + MDIV2, Y + MDIV2] / Total;
Gaussian[X + MDIV2, Y + MDIV2] := KS;
// グリッドに表示
SigumaGrid.Cells[X + MDIV2 + 1, Y + MDIV2 + 1] := FloatTostrF(KS, ffFixed, 6,5);
end;
application.ProcessMessages;
Result := True;
end;
type
TPrgbarry = array[0..0] of Trgbtriple; // 24ビットカラーレコード 32ビット用はTRGBQuadArray
Prgbarray = ^TPrgbarry; // ポインター
// 配列のポインターが必要なだけなので、長さは1 [0..0]で問題ありません。
var
conv_real : D2array; // ラプラシアン、ガウスフィルター用実数部
conv_image : D2array; // ラプラシアン、ガウスフィルター用虚数部
half_cw : integer;
half_ch : integer;
arBlu : D2array; // データ実数部（入出力兼用)
aiBlu : D2array; // データ虚数部（入出力兼用)
r, i : Double;
PBA : Prgbarray;
DI : Integer;
rdr : Double;
rdi : Double;
//------------------------------
// フィルター計算選択
//------------------------------
procedure TForm1.MaskingBtnClick(Sender: TObject);
begin
MaskingBtn.Enabled := False;
FileOpenBtn.Enabled := False;
FileSaveBtn.Enabled := False;
RepaintBtn.Enabled := False;
sourceBtn.Enabled := False;
Memo1.Clear;
Memo1.Lines.Add('計算時間');
Memo1.Lines.Add(' Start ' + DateTimeToStr(Now));
Timer1.Enabled := True;
end;
//----------------------------
// 再計算許可
//----------------------------
procedure TForm1.SigumaEditChange(Sender: TObject);
begin
MaskingBtn.Enabled := True;
end;
procedure TForm1.RadioGroup1Click(Sender: TObject);
begin
MaskingBtn.Enabled := True;
end;
procedure TForm1.RangeEditChange(Sender: TObject);
begin
MaskingBtn.Enabled := True;
end;
procedure TForm1.UnsharpEditChange(Sender: TObject);
begin
MaskingBtn.Enabled := True;
end;
//--------------------------------
// ボタン表示更新待ち合わせ
//--------------------------------
procedure TForm1.Timer1Timer(Sender: TObject);
begin
Timer1.Enabled := False;
if RadioGroup1.ItemIndex = Laplacian then LaplacianMask;
if (RadioGroup1.ItemIndex = Gaussians) or
(RadioGroup1.ItemIndex = Deconvolution) or
(RadioGroup1.ItemIndex = Unsharp) then GaussianRoutine;
// MaskingBtn.Enabled := True;
FileSaveBtn.Enabled := True;
FileOpenBtn.Enabled := True;
RepaintBtn.Enabled := True;
sourceBtn.Enabled := True;
fulbtn.Enabled := True;
Memo1.Lines.Add(' End ' + DateTimeToStr(Now));
end;
//------------------------------------------
// ラプラシアン計算
//------------------------------------------
procedure TForm1.LaplacianMask;
const
// ラプラシアンフィルター
cvnLap = 3;
convLap : array[0..cvnLap - 1, 0..cvnLap - 1] of double = ((1, 1, 1),
(1, -8, 1),
(1, 1, 1));
MatSize = cvnLap;
var
nx, ny, x, y: integer;
MDIV2 : integer;
KS : double;
begin
// ストリンググリッドの設定
SigumaGrid.ColCount := MatSize + 1;
SigumaGrid.RowCount := MatSize + 1;
SigumaGrid.ClientWidth := (SigumaGrid.DefaultColWidth + 1) * (MatSize + 1) - 1;
SigumaGrid.ClientHeight := (SigumaGrid.DefaultRowHeight + 1) * (MatSize + 1) - 1;
// 固定セルにナンバーリング
MDIV2 := MatSize div 2;
for Y := - MDIV2 to MDIV2 do SigumaGrid.Cells[0 , Y + MDIV2 + 1] := intTostr(Y);
for X := - MDIV2 to MDIV2 do SigumaGrid.Cells[X + MDIV2 + 1, 0] := intTostr(X);
// フィルター値表示
for Y := - MDIV2 to MDIV2 do
for X := - MDIV2 to MDIV2 do begin
KS := convLap[X + MDIV2, Y + MDIV2];
// グリッドに表示
SigumaGrid.Cells[X + MDIV2 + 1, Y + MDIV2 + 1] := FloatTostrF(KS, ffFixed, 2,1);
end;
application.ProcessMessages;
// 配列の確保
setlength(conv_real, GWidth, GHeight);
setlength(conv_image, GWidth, GHeight);
setlength(arBlu, GWidth, GHeight);
setlength(aiBlu, GWidth, GHeight);
// 配列のクリア
for y := 0 to GHeight - 1 do
for x := 0 to GWidth - 1 do begin
conv_real[x, y] := 0;
conv_image[x, y] := 0;
end;
// ラプラシアンフィルターの設定
half_cw := cvnLap div 2;
half_ch := cvnLap div 2;
for y := 0 to cvnLap - 1 do
for x := 0 to cvnLap -1 do begin
nx := GWidth - half_cw + x;
ny := Gheight - half_ch + y;
if nx >= GWidth then nx := nx - GWidth;
if ny >= Gheight then ny := ny - GHeight;
conv_real[nx][ny] := convLap[x][y];
end;
// 原画像を読み込みグレー変換
// ２次元DFTの入力データに変換します
for y := 0 to Gheight - 1 do begin
PBA := InputDBitmap.ScanLine[Y];
for x := 0 to GWidth - 1 do begin
// ラプラシアン用グレイ画像生成 輝度変換 実数部
arBlu[x, y] := Round(PBA[x].rgbtBlue * kb
+ PBA[x].rgbtGreen * kg
+ PBA[x].rgbtRed * kr);
aiBlu[x, y] := 0.0; // 虚数部
end;
end;
// 拡張したフィルターをフーリエ変換します
// 注 フィルターの場合、変換時、データー数による除算処理をしないようにします
DFT(conv_real, conv_image, 2, GHeight, GWidth);
// ２次元DFTを実行し 原画像を周波数成分に変換します
DFT(arBlu, aiBlu, 1, GHeight, GWidth);
// ガウシアン ラプラシアン フィルター処理
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
// フィルターの実施
r := arBlu[x][y] * conv_real[x][y] - aiBlu[x][y] * conv_image[x][y];
i := arBlu[x][y] * conv_image[x][y] + aiBlu[x][y] * conv_real[x][y];
arBlu[x][y] := r;
aiBlu[x][y] := i;
end;
// 逆DFTを実行して，フィルタ処理された周波数成分を画像に戻します
DFT(arBlu, aiBlu, -1, GHeight, GWidth);
// 結果を画像データに変換
for Y := 0 to Gheight - 1 do begin
PBA := OutputBitmap.ScanLine[Y];
for X := 0 to GWidth - 1 do begin
DI := Round(arBlu[x][y]);
if DI > 255 then DI := 255;
if DI < 0 then DI := 0;
PBA[X].rgbtBlue := DI;
PBA[X].rgbtGreen := DI;
PBA[X].rgbtRed := DI;
end;
end;
Imageout(OutputBitmap); // 出力枠に変倍出力
end;
//------------------------------------------
// ガウスフィルターによるボケ計算 Gaussian
// デコンボリューション Deconvolution
// アンシャープ Unsharp
//------------------------------------------
procedure TForm1.GaussianRoutine;
var
nx, ny, x, y: integer;
Step : integer;
begin
// ガウシアンフィルターデーターの生成
if not SigumaMatSet then exit;
// 配列の確保
setlength(conv_real, GWidth, GHeight);
setlength(conv_image, GWidth, GHeight);
setlength(arBlu, GWidth, GHeight);
setlength(aiBlu, GWidth, GHeight);
// 配列のクリア
for y := 0 to GHeight - 1 do
for x := 0 to GWidth - 1 do begin
conv_real[x, y] := 0;
conv_image[x, y] := 0;
end;
// 拡張フィルターの設定
half_cw := gaus div 2;
half_ch := gaus div 2;
for y := 0 to gaus - 1 do
for x := 0 to gaus -1 do begin
nx := GWidth - half_cw + x;
ny := Gheight - half_ch + y;
if nx >= GWidth then nx := nx - GWidth;
if ny >= Gheight then ny := ny - GHeight;
conv_real[nx][ny] := Gaussian[x][y];
end;
// 拡張したフィルターをフーリエ変換します
// 注 フィルターの場合、変換時、データー数による除算処理をしないようにします
DFT(conv_real, conv_image, 2, GHeight, GWidth);
// 原画像を読み込み，２次元DFTの入力データに変換します
for Step := 0 to 2 do begin
// 原画像を読み込み，２次元DFTの入力データに変換します
for y := 0 to Gheight - 1 do begin
PBA := InputDBitmap.ScanLine[Y];
for x := 0 to GWidth - 1 do begin
// カラー用三色分解 実数部
case Step of
0 : arBlu[x, y] := PBA[x].rgbtBlue;
1 : arBlu[x, y] := PBA[x].rgbtGreen;
2 : arBlu[x, y] := PBA[x].rgbtRed;
end;
aiBlu[x, y] := 0.0; // 虚数部
end;
end;
// ２次元DFTを実行し 原画像を周波数成分に変換します
DFT(arBlu, aiBlu, 1, GHeight, GWidth);
// ガウシアン ボカシ処理 ---------------------------
if RadioGroup1.ItemIndex = Gaussians then begin
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
// フィルターの実施
r := arBlu[x][y] * conv_real[x][y] - aiBlu[x][y] * conv_image[x][y];
i := arBlu[x][y] * conv_image[x][y] + aiBlu[x][y] * conv_real[x][y];
arBlu[x][y] := r;
aiBlu[x][y] := i;
end;
end;
//---------------------------------------------------
// デコンボリューション処理 -------------------------
if RadioGroup1.ItemIndex = Deconvolution then begin
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
// デコンボリューション係数
rdr := abs(conv_real[x][y]);
rdi := abs(conv_image[x][y]);
// ボカシ処理とは逆なので、逆数にします
// thresholdは値が小さすぎた場合逆数の値が大きくなりすぎるので制限する値です
// thresholdはゼロでの除算の防止もします
if rdr > threshold then rdr := 1 / conv_real[x][y]
else rdr := 0;
if rdi > threshold then rdi := 1 / conv_image[x][y]
else rdi := 0;
// デコンボリューションの生成
r := arBlu[x][y] * rdr - aiBlu[x][y] * rdi;
i := arBlu[x][y] * rdi + aiBlu[x][y] * rdr;
arBlu[x][y] := r;
aiBlu[x][y] := i;
end;
end;
//---------------------------------------------------
// アンシャープ処理 ---------------------------------
if RadioGroup1.ItemIndex = Unsharp then begin
for x := 0 to GWidth - 1 do
for y := 0 to GHeight - 1 do begin
// ボケ画像の生成
r := arBlu[x][y] * conv_real[x][y] - aiBlu[x][y] * conv_image[x][y];
i := arBlu[x][y] * conv_image[x][y] + aiBlu[x][y] * conv_real[x][y];
// 差分計算
r := arBlu[x][y] - r;
i := aiBlu[x][y] - i;
// 元データーに差分加算
arBlu[x][y] := arBlu[x][y] + r;
aiBlu[x][y] := aiBlu[x][y] + i;
end;
end;
//---------------------------------------------------
// 逆DFTを実行して，フィルタ処理された周波数成分を画像に戻します
DFT(arBlu, aiBlu, -1, GHeight, GWidth);
// 結果を画像データに変換
for Y := 0 to Gheight - 1 do begin
PBA := OutputBitmap.ScanLine[Y];
for X := 0 to GWidth - 1 do begin
DI := Round(arBlu[x][y]);
if DI > 255 then DI := 255;
if DI < 0 then DI := 0;
case Step of
0: PBA[X].rgbtBlue := DI;
1: PBA[X].rgbtGreen := DI;
2: PBA[X].rgbtRed := DI;
end;
end;
end;
end;
Imageout(OutputBitmap); // 出力枠に変倍出力
end;
//--------------
// 元画像の表示
//--------------
procedure TForm1.sourceBtnClick(Sender: TObject);
begin
Imageout(InputDBitmap); // 出力枠に変倍出力
end;
//----------------------------
// 計算結果の再表示
//----------------------------
procedure TForm1.RepaintBtnClick(Sender: TObject);
begin
Imageout(OutputBitmap); // 出力枠に変倍出力
end;
//------------------------------------------------
// ファイルのオープンとグレイデーターの生成
//------------------------------------------------
procedure TForm1.FileOpenBtnClick(Sender: TObject);
var
WIC : TWICImage;
begin
VRect := Rect(0, 0, Image1.Width, Image1.Height);
Image1.Canvas.Brush.Style := bsSolid;
Image1.Canvas.Brush.Color := clBtnface;
Image1.Canvas.FillRect(VRect); // Canvas 画像消去
OpenPictureDialog1.Filter := OpenFileFilter; // ファイルオープンフィルターの設定
if OpenPictureDialog1.Execute = true then // ファイルが指定されたら
begin
WIC := TWICImage.Create; // TWICImageの生成
try
InFilename := OpenPictureDialog1.FileName; // ファイル名の取得
WIC.LoadFromFile(InFilename); // 画像の読み込み
GHeight := WIC.Height; // 画像高さ取得
GWidth := WIC.Width; // 画像幅
Image1.Width := GWidth;
Image1.Height:= GHeight;
Image1.Picture.Bitmap.SetSize(GWidth, GHeight);
InputDBitmap.Width := GWidth;
InputDBitmap.Height := GHeight;
InputDBitmap.Canvas.Draw(0, 0, WIC); // ビットマップに描画フォーマット変換
finally
WIC.Free; // TWICImage 解放
end;
end
else exit;
OutputBitmap.Width := GWidth;
OutputBitmap.Height := GHeight;
trimming_imageout(InputDBitmap); // 出力枠合わせてに変倍出力
// Imageout(InputDBitmap); // 出力枠に変倍出力
MaskingBtn.Enabled := True;
sourceBtn.Enabled := True;
radiogroup1.Enabled := True;
fulbtn.Enabled := False;
end;
//------------------------------
// 画像のファイルへの保存
//------------------------------
procedure TForm1.FileSaveBtnClick(Sender: TObject);
var
WIC : TWicImage;
WICF : TWicImageFormat;
Fname : String;
ExeStr : String;
FnameTop: String;
Findex : integer;
function WFormatSet: Boolean; // 拡張子によるファイルフォーマット設定
begin
Result := false;
ExeStr := LowerCase(ExeStr);
if ExeStr = '.jpg' then begin WICF := Wifjpeg; Result := True; end;
if ExeStr = '.jpeg' then begin WICF := Wifjpeg; Result := True; end;
if ExeStr = '.tif' then begin WICF := Wiftiff; Result := True; end;
if ExeStr = '.tiff' then begin WICF := Wiftiff; Result := True; end;
if ExeStr = '.png' then begin WICF := Wifpng; Result := True; end;
if ExeStr = '.gif' then begin WICF := Wifgif; Result := True; end;
if ExeStr = '.bmp' then begin WICF := Wifbmp; Result := True; end;
if ExeStr = '.wdp' then begin WICF := WifWMPhoto; Result := True; end;
if ExeStr = '.hdp' then begin WICF := WifWMPhoto; Result := True; end;
end;
begin
SavePictureDialog1.Filter := SaveFileFilter;
// SavePictureDialog1.DefaultExt := GraphicExtension(TWicImage);
if not SavePictureDialog1.Execute then exit;
ExeStr := ExtractFileExt(SavePictureDialog1.FileName);
if ExeStr = '' then begin // 拡張子がなかったら
Findex := SavePictureDialog1.FilterIndex; // FilterIndexによる拡張子の設定
case Findex of
1, 3 : Fname := ChangeFileExt(SavePictureDialog1.FileName,'.jpg'); // 拡張子の設定
2 : Fname := ChangeFileExt(SavePictureDialog1.FileName,'.png'); // 拡張子の設定
4 : Fname := ChangeFileExt(SavePictureDialog1.FileName,'.gif'); // 拡張子の設定
5 : Fname := ChangeFileExt(SavePictureDialog1.FileName,'.bmp'); // 拡張子の設定
6 : Fname := ChangeFileExt(SavePictureDialog1.FileName,'.tif'); // 拡張子の設定
7 : Fname := ChangeFileExt(SavePictureDialog1.FileName,'.wdp'); // 拡張子の設定
end;
end
else
Fname := SavePictureDialog1.FileName;
ExeStr := ExtractFileExt(Fname); // 拡張子だけ取り出し
if not WFormatSet then begin // 拡張子によるファイルフォーマット設定と確認
application.MessageBox('ファイルの拡張子が間違っています。','注意', 0);
exit;
end;
FnameTop := ExtractFileName(Fname); // ファイル名だけ取り出し
if Length(FnameTop) = Length(ExeStr) then begin // ファイル名の長さ確認
application.MessageBox('ファイル名がありません。','注意', 0);
exit;
end;
if FileExists(Fname) then // ファイル名によるファイル検索
if MessageDlg('既に同じ名前のファイルがあります上書きしますか ' + ExtractFileName(Fname) + '?',
mtConfirmation, [mbYes, mbNo], 0, mbNo) = IDNo then exit;
WIC := TWicImage.Create; // TWicImage生成
try
WIC.Assign(OutputBitmap); // TWicImageにビットマップデーター割り付け
WIC.ImageFormat := WICF; // 保存フォーマットセット
WIC.SaveTofile(Fname); // ファイルの書き出し
finally
WIC.Free; // TWicImage解放
end;
end;
// 初期設定
procedure TForm1.FormCreate(Sender: TObject);
begin
Timer1.Enabled := False;
Image1.Width := ImageHWC div 2;
Image1.Height := ImageHWC div 2;
ScrollBox1.Height := ImageHWC;
ScrollBox1.Width := ImageHWC;
ScrollBox1.Height := ImageHWC + ScrollBox1.BevelWidth * 4;
ScrollBox1.Width := ImageHWC + ScrollBox1.BevelWidth * 4;
Image1.Top := 0;
Image1.Left := 0;
InputDBitmap := TBitmap.Create;
OutputBitmap := TBitmap.Create;
InputDBitmap.PixelFormat := pf24bit;
OutputBitmap.PixelFormat := pf24bit;
radiogroup1.Items.Add('Laplacian');
radiogroup1.Items.Add('Gaussian');
radiogroup1.Items.Add('Deconvolution');
radiogroup1.Items.Add('Unsharp');
radiogroup1.ItemIndex := 2;
radiogroup1.Enabled := False;
MaskingBtn.Enabled := False;
FileSaveBtn.Enabled := False;
sourceBtn.Enabled := False;
RepaintBtn.Enabled := False;
fulbtn.Enabled := False;
end;
// 終了処理
procedure TForm1.FormDestroy(Sender: TObject);
begin
InputDBitmap.Free;
OutputBitmap.Free;
end;
procedure TForm1.FulBtnClick(Sender: TObject);
begin
trimming_imageout(OutputBitmap); // 出力枠合わせてに変倍出力
end;
end.
Deconvolution.zip