Program 1
このページを見るときは固定等幅フォントを使用してください
/****************************************************************/
/* ファジールールに基づくメジアンフィルタ */
/* A Vol.J78-A No.2 pp.123-131 1995 */
/****************************************************************/
/* @001 Initial 1998/01/16 TK */
/* @002 Membership 関数の修正 1998/01/23 TK */
/* @003 雑音部，信号部各々の MSE の算出 1998/01/23 TK */
/****************************************************************/
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "gf_file.h" // BMP File 書き出し
#define VER "1.03a" // Version Number
#define DEBUG 0
#define DEBUGF 0
#define INNAME 0 // ファイル名を入力するようにする
#define WRITE 1 // 画像出力をファイルに記録する
#if WRITE
#define BMPFILE 0 // BMP でのファイルの書き出し
#endif
//-------- 画像の選択 ---------
#define IMG_LENA 1 // Lena
#define IMG_BOAT 0 // Boat
#define IMG_LIGHTHOUSE 0 // LIGHTHOUSE
/* 各種定義 */
#define xsize 256 // 画像の横のサイズ
#define ysize 256 // 画像の縦のサイズ
#define winsize 5 // 局所情報窓サイズ
//--- プロトタイプ宣言 ---//
int file_read( char * , double * ) ; // File read function
void fmedian( int, int ) ; // Fuzzy UM Process
double mse( void ) ; // MSE Caliculation
void mapmse( void ) ; // 雑音部と信号部の MSE を別々に算出
//--- 各種配列宣言 --//
double org[ysize][xsize] ; // 理想画像配列
double id[ysize][xsize] ; // 原画像配列and入力画像配列
double od[ysize][xsize] ; // 出力画像配列
/* FIle Name */
#if IMG_LENA
char ogname[] = "lena" ; // 原信号ファイル名
char iname[] = "lena_i5.fob" ; // 入力(劣化）信号ファイル名
#if BMPFILE
char oname[] = "lena_i5r.bmp" ; // 出力信号ファイル名
#else
char oname[] = "lena_i5r.fob" ; // 出力信号ファイル名
#endif
#endif
void main(void)
{
FILE *ofp , *fp ; // 画像書込用 , 関数最適結果用
char dummy[100]; // Dummy 変数
int x, y ; // Counter
double minbl , maxbl ; // 局所情報の分散値の最小と最大値
double dv , bv , msed ; // DV / BV 用変数 , MSE 用変数
char ver[] = VER ; // Version Number
#if WRITE
int wdat ; // ファイル書込用変数
#endif
/******************************************
* 画像処理プログラム *
******************************************/
printf("ファジールールに基づくメジアンフィルタ V%s [1998.01.15]\n" , &ver[0] );
printf(" Copyright MIT [Tomoaki Kimura] 1997\n");
//--- 入力配列 id[y][x] に輝度値が入る ---//
#if INNAME
printf("劣化画像") ;
#endif
if ( file_read( &iname[0] , &id[0][0]) ) exit(1) ; // 劣化画像読み込み
#if INNAME
printf("理想画像名") ;
#endif
if ( file_read( &ogname[0] , &org[0][0] ) ) exit(1) ; // 理想画像読み込み
//-- 配列の初期化 --//
printf("配列の初期化\n") ;
for(y=0 ; y＜ysize ; y++) for(x=0 ; x＜xsize ; x++) od[y][x] = 0.0 ;
//-- use FUZZY --//
printf("ファージメジアンフィルタの計算\n") ;
for( y=2 ; y＜ysize-2 ; y++ ){
for( x=2 ; x＜xsize-2 ; x++ ){
#if DEBUG
printf("od[%d][%d]", y, x ) ;
#endif
fmedian( y, x ) ; // FUZZY による Median 処理
#if DEBUG
printf("od[%d][%d](%lf) , id(%lf) : ", y, x, od[y][x] , id[y][x] ) ;
#endif
} // X
} // Y
//********************************//
// MSE の算出 //
//********************************//
printf("MSE の計算\n") ;
msed = mse() ; // MSE
mapmse() ; // 雑音部と信号部各々の MSE
//******************************************//
// 出力ファイル //
//******************************************//
#if WRITE
printf("算出後のデータをファイルに出力しています\n") ;
#if BMPFILE
Gf_write_bmp( &oname[0] , &od[0][0] ) ;
#else // (BMPFILE)
printf("算出後のデータをファイルに出力しています\n") ;
ofp=fopen(oname,"wb");
for( y=0 ; y＜ysize ; y++ ){
for( x=0 ; x＜xsize ; x++ ){
if ( od[y][x] ＞ 255 ) od[y][x] = 255 ;
if ( od[y][x] ＜ 0 ) od[y][x] = 0 ;
wdat = (int)od[y][x] ; // データを変換する
putc( wdat, ofp ) ;
putc( (wdat＞＞8), ofp ) ;
}
}
fclose(ofp);
#endif // (BMPFILE)
#endif // (WRITE)
return;
} //-- end of main() --//
//////////////////////////////////////////////////////////////////////////
// SUB Routine //
//////////////////////////////////////////////////////////////////////////
/****************************/
/* File Read Function */
/****************************/
int file_read( char *fname , double *inp )
{
int x , y ; // 単なるカウンタ
int dat ; // データ用一時変数
#if INNAME
int sw=0 ; // File name 用 SW
#endif
FILE *fp ;
#if INNAME
do {
printf ( "ファイル名を入力してください：" ) ;
gets(fname);
if((fp=fopen( fname,"rb"))==NULL) printf("入力ファイルが見つかりません。") ;
else sw = 1 ;
} while ( !sw ) ;
#else
if((fp=fopen( fname,"rb"))==NULL){
printf("入力ファイルが見つかりません。");
return(1) ; // 異常終了
}
#endif
//--- input memory from file data ---//
for( y=0 ; y＜ysize ; y++ ){
for( x=0 ; x＜xsize ; x++){
dat = (double)getc(fp);
dat += (int )(getc(fp)＜＜8);
*inp = (double)dat ;
inp ++ ;
}
}
fclose(fp);
return(0) ; // 正常終了
}
//======================================================================//
// Fuzzy による メジアンフィルタ 処理 //
// //
// 入力変数 //
// int y, x 画像の縦横 //
// //
// 外部変数として以下の配列を指定する //
// id[][] : 入力画像 //
// od[][] : 出力画像 //
//======================================================================//
void fmedian( int y , int x )
{
//-- use FUZZY UM --//
double f[25] , g[9] ; // 一時的に局所情報に使用する変数
double m , p ; // 各種情報用変数
double ssort , med ; // ソート及びメジアン用変数
double tkim , tkip ; // 各情報のメンバーシップ適合度
double u[4] , uw[4] ; // 各メンバーシップの適合度
double sumu , sumuw , yr ;
int i , j , k ; // カウンタ変数
//==============================//
// FUZZY による Median 処理 //
//==============================//
/* 局所情報の各値をとってくる */
for ( i=-2 , k=0 ; i＜=2 ; i++ ) {
for ( j=-2 ; j＜=2 ; j++ , k++ ) f[k] = id[y+i][x+j] ; // 5x5
}
for ( i=-1 , k=0 ; i＜=1 ; i++ ) {
for ( j=-1 ; j＜=1 ; j++ , k++ ) g[k] = id[y+i][x+j] ; // 3x3
}
/*** メジアン値を求める ***/
for( j=0 ; j＜(25-1) ; j++ ){
for( i=j+1 ; i＜25 ; i++ ){
if( f[j] ＞ f[i]){
ssort = f[j] ; // 小さい方を前に持ってくる
f[j] = f[i] ;
f[i] = ssort ;
}
}
}
med = f[12] ; // Median 値の保存
#if DEBUG
printf(" median:[%lf] " , med ) ;
#endif
/*** 入力とメジアンの差に関するルール m ***/
m = fabs( id[y][x] - med ) ;
#if DEBUG
printf(" m:[%lf]\n " , m ) ;
#endif
/*** 入力の孤立に関するルール p ***/
#if DEBUG & DEBUGF
printf("\n 3x3 の情報 : ") ;
for ( i=0 ; i＜9 ; i++ ) printf("[%d]" , (int)g[i] ) ;
#endif
for ( i=0 ; i＜9 ; i++ ) g[i] = fabs( g[i] - id[y][x] ) ;
for( j=0 ; j＜(9-1) ; j++ ){
for( i=j+1 ; i＜9 ; i++ ){
if( g[j] ＞ g[i]){
ssort = g[j] ; // 小さい方を前に持ってくる
g[j] = g[i] ;
g[i] = ssort ;
}
}
}
p = ( (g[1]+g[2]) / 2 ) ;
#if DEBUG
printf("\n Sort 後の 3x3 の情報 id[x][y][%d] p[%lf]: " , (int)id[y][x], p ) ;
#if DEBUGF
for ( i=0 ; i＜9 ; i++ ) printf("[%d]" , (int)g[i] ) ;
#endif
#endif
//----------------------//
// メンバーシップ関数 //
//----------------------//
//--- m における適合度の設定 ---//
if ( m ＜= 20.0 ) tkim = 1.0 ;
else if ( m ＜= 60.0 ) tkim = 0.9 ;
else if ( m ＜= 80.0 ) tkim = 0.8 ;
else if ( m ＜= 100.0 ) tkim = 0.2 ;
else tkim = 0.1 ;
//--- p における適合度の設定 ---//
if ( p ＜= 10.0 ) tkip = 1.0 ;
else if ( p ＜= 20.0 ) tkip = 0.7 ;
else if ( p ＜= 30.0 ) tkip = 0.3 ;
else if ( p ＜= 50.0 ) tkip = 0.2 ;
else tkip = 0.1 ;
#if DEBUG
printf("m()=%lf : %lf p()=%lf :%lf \n", m, tkim, p, tkips );
#endif
//------------------------//
// ファジイルールの設定 //
//------------------------//
yr = tkim * tkip ; // 処理係数の設定
#if DEBUG
printf("%lf %lf %lf %lf \n",u[0], u[1], u[2], u[3] ) ;
#endif
od[y][x] = ( med + yr*(id[y][x]-med) ) ; // 最終的な画像の値
#if DEBUG
printf("uw[%lf] u[%lf] yr[%lf] id[%d] med[%lf] od[%d]\n", sumuw, sumu, yr, (int)id[y][x], med, (int)od[y][x] ) ;
getch() ;
#endif
return ;
}
//======================================================================//
// MSE の算出 //
// //
// 外部変数として以下の配列を指定する //
// org[][] : オリジナル画像 //
// od[][] : DV/BV を算出する画像 //
//======================================================================//
double mse( void )
{
int x, y ; // Counter
double mse , tmse ; // MSE 計算用
tmse = 0.0 ;
for( y=2 ; y＜ysize-2 ; y++ ){
for( x=2 ; x＜xsize-2 ; x++ ){
mse = ( od[y][x] - org[y][x] ) ;
tmse += ( mse * mse ) ;
}
}
tmse /= ( (double)(ysize-2) * (double)(xsize-2) ) ;
printf(" MSE = %f \n" , tmse ) ;
return ( tmse ) ;
}
//======================================================================//
// 雑音部と信号部それぞれにおけるMSE の算出 //
// //
// 外部変数として以下の配列を指定する //
// org[][] : オリジナル画像 //
// od[][] : MSE を算出する画像 //
// map[][] : 雑音部と信号部を分ける MAP //
//======================================================================//
int map[ysize][xsize] ; // 雑音部と信号部を分ける Map 用配列
void mapmse( void )
{
int x, y ; // Counter
double mse , tmsen , tmses ; // MSE 計算用
double nsize ; // 画像におけるノイズの数
//--- MAP 作成 (雑音部と信号部を分ける) ---//
for( y=0 ; y＜ysize ; y++ ){
for( x=0 ; x＜xsize ; x++ ){
if ( id[y][x] == org[y][x] ) map[y][x] = 0 ;
else map[y][x] = 1 ;
}
}
tmses = tmsen = nsize = 0.0 ;
for( y=2 ; y＜ysize-2 ; y++ ){
for( x=2 ; x＜xsize-2 ; x++ ){
if ( map[y][x] ) { // 雑音部分の MSE
mse = ( od[y][x] - org[y][x] ) ;
tmsen += ( mse * mse ) ;
nsize ++ ; // 雑音部のカウント
}
else { // 信号部分の MSE
mse = ( od[y][x] - org[y][x] ) ;
tmses += ( mse * mse ) ;
}
}
}
tmsen /= nsize ;
tmses /= ( (double)(ysize-2) * (double)(xsize-2) - nsize ) ;
printf(" MSE(Noise)[%f] = %f , MSE(Signal) = %f \n" , nsize , tmsen , tmses ) ;
return ;
}
/************************************************************************/
/* ＜＜ ＢＭＰ画像ファイル書き込みルーチン ( 8 bit save only) ＞＞ */
/* 対象: 256x256x65536 は 256x256x256 のデータに変換する */
/* */
/* input *fn 出力ファイル名ポインタ */
/* *rdat 計算後のデータポインタ */
/* */
/* return none */
/* */
/************************************************************************/
void Gf_write_bmp( char *fn , double *rdat )
{
char out_fn[14] ;
char *inp_fn = Inp_fn ;
unsigned long i ;
int j ;
BYTE fdat[256] ; // 一時画像エリアの確保
FILE *fp ;
unsigned long cv ;
printf(" ＢＭＰファイルの書き込み\n") ;
//--- 書出ファイルを ８ビットに変換する ---//
gf_write_to8chg( rdat ) ;
//--- ファイルのＯＰＥＮ --//
strcpy( Out_fn , fn ) ; // 入力ファイル名の設定
if ( Out_fn[0] == '\0' ) { // ファイル名の確認
printf("\nPlase Input OutputFile Name :") ; //
scanf("%s", &Out_fn[0] ) ; // ファイル名の入力
} //
printf(" Output File Name : [%s]\n" , Out_fn ) ; //
fp = fopen( Out_fn,"wb" ) ; // ファイルの OPEN
//--- ファイルの書き込み状態バーの表示等 --//
cv = 256.0 / 20.0 ;
printf("0%% 50%% 100%%\n") ;
printf("|--------+----------|\n") ;
//--- ヘッダファイルの書き込み ---//
put_header_bmp( fp , 0x100 , 0x100 ) ;
put_palette_bmp( fp ) ;
//--- ファイルの書き込み --//
for ( i=0 ; i＜Memsize ; i++ ) rdat ++ ; // 画像右下を始値にする
for ( i=0 ; i＜256 ; i++ ) {
for ( j=255 ; j＞=0 ; j-- ) {
rdat -- ;
fdat[j] = ( (BYTE)*rdat & 0xff ) ;
}
fwrite( fdat, 256, 1, fp) ; // 下位８ビットを書き出す
//--- ファイル書込状態の表示 ---//
if ( cv ＜ i ) {
printf("o") ;
cv += ( 256.0 / 20.0 ) ;
}
}
printf("\n") ;
fclose(fp) ;
return ;
}
/************************************************************************/
/* ＜＜ 書き出し用画像ファイルを 8 bit に変換するルーチン ＞＞ */
/* */
/* input *rdat 取込用データエリア */
/* */
/************************************************************************/
void gf_write_to8chg( double *wdat )
{
double max , cn, tmp ;
double *rdat = wdat ;
long i ;
printf("書出画像を８ビットに変換\n") ;
//--- 最高値の検出 ---//
for ( i=0 ; i＜Memsize ; i++ ) {
if ( max ＜ *rdat ) max = *rdat ;
rdat ++ ;
}
printf("書き出しファイルの最高値:[%f]", max ) ;
#if DEBUG
getch() ; // 一時停止
#endif
rdat = wdat ; // ポインタを最初の位置へ
if ( max ＞ 255.0 ) {
printf("．．変換中\n") ;
cn = max / 255.0 ; // 正規化係数の算出
for ( i=0 ; i＜Memsize ; i++ , rdat++ ) {
tmp = ( *rdat / cn ) ;
*rdat = tmp ; // データの保存
#if DEBUG
printf("(%lf)", *rdat) ;
#endif
}
}
else printf("．．無変換．．\n") ;
#if DEBUG
getch() ; // 一時停止
#endif
return ;
}
/************************************************************************/
/* Bit Map File 変換用 関数 */
/************************************************************************/
/************************************************************************/
/* ＜＜ BMP ヘッダ書出関数 (8 bit only) ＞＞ */
/* 色ビット数[bit/dot]は (1,4,8,24) のうち 8 bit のみ */
/* */
/* input *sput ファイルポインタ */
/* rx 画像Ｘサイズ [dot](1〜) */
/* ry 画像Ｙサイズ [dot](1〜) */
/* */
/* return 0 nomalend */
/* 1 画面サイズエラー */
/* 2 ファイル書出エラー */
/* */
/************************************************************************/
int put_header_bmp( FILE *sput , int rx , int ry )
{
int i , // ループ変数
color = 256 ; // 色数 [ 8 bit なので色数は 256 ]
unsigned long headsize , // ヘッダサイズ [byte]
imagesize ; // 画面サイズ
if ( (rx＜=0) || (ry＜=0) ) return(1) ; // 画面サイズの確認
if ( (sput==NULL) || ferror(sput) ) return(2) ; // ファイル書出の確認
headsize = 14 + 40 + 4*color ; // BMP File headre + BMP 情報 header + 色数
imagesize = (unsigned long)rx * (unsigned long)ry ;
//-- ビットマップ・ファイル・ヘッダの書出 --//
fputs("BM", sput) ; // 識別文字列
fputc4lh( headsize+imagesize , sput) ; // ファイルサイズ[byte]
fputc2lh( 0 , sput ) ; // 予約領域[byte]
fputc2lh( 0 , sput ) ; // 予約領域[byte]
fputc4lh( headsize , sput ) ; // ヘッダサイズ[byte]
if ( ferror(sput) ) return(2) ; // ファイル書出の確認
//-- ビットマップ情報ヘッダの書出 --//
fputc4lh( 40 , sput ) ; // 情報サイズ[byte]
fputc4lh( rx , sput ) ; // 画像Ｘサイズ [dot]
fputc4lh( ry , sput ) ; // 画像Ｙサイズ [dot]
fputc2lh( 1 , sput ) ; // 面数
fputc2lh( 8 , sput ) ; // 色ビット数[bit/dot] (8 bit)
fputc4lh( 0 , sput ) ; // 圧縮方式
fputc4lh( imagesize , sput ) ; // 圧縮サイズ [byte]
fputc4lh( 0 , sput ) ; // 水平解像度
fputc4lh( 0 , sput ) ; // 垂直解像度
fputc4lh( color , sput ) ; // 色数
fputc4lh( color , sput ) ; // 重要色数
if ( ferror(sput) ) return(2) ; // ファイル書出の確認
return(0) ;
}
/************************************************************************/
/* ＜＜ BMP パレット書出関数 (白黒256階調) ＞＞ */
/* */
/* input *sput ファイルポインタ */
/* */
/* return 0 nomalend */
/* 2 ファイル書出エラー */
/* */
/************************************************************************/
int put_palette_bmp( FILE *sput )
{
int i ;
for ( i=0 ; i＜256 ; i++ ) {
putc( i , sput ) ; // 青
putc( i , sput ) ; // 緑
putc( i , sput ) ; // 赤
putc( 0 , sput ) ; // α
}
if ( ferror(sput) ) return(2) ; // ファイル書出の確認
return(0) ;
}
/************************************************************************/
/* ＜＜ ファイル２バイト低高書出関数 ＞＞ */
/* */
/* input rdat 取込用データエリア */
/* *sput ファイルポインタ */
/* */
/************************************************************************/
int fputc2lh( WORD d , FILE *sput )
{
putc( d & 0xff , sput ) ;
return( putc( d ＞＞ CHAR_BIT, sput ) ) ;
}
/************************************************************************/
/* ＜＜ ファイル４バイト低高書出関数 ＞＞ */
/* */
/* input rdat 取込用データエリア */
/* *sput ファイルポインタ */
/* */
/************************************************************************/
int fputc4lh( unsigned long d , FILE *sput )
{
putc( d & 0xff , sput ) ;
putc( (d ＞＞ CHAR_BIT) & 0xff , sput ) ;
putc( (d ＞＞ CHAR_BIT*2) & 0xff , sput ) ;
return
putc( (d ＞＞ CHAR_BIT*3) & 0xff , sput ) ;
}
/*****************************************************************************/
/* EEEEEEEE NN NN DDDDD */
/* EE NNN NN DD DD */
/* EEEEEE NN N NN DD DD */
/* EE NN NNN DD DD */
/* EEEEEEEE NN NN DDDDDDDD */
/*****************************************************************************/