2004 2005 2 2 1G01P038-0

1 2 1.1.............................. 2 1.2......................... 2 1.3......................... 3 2 4 2.1............................ 4 2.2....................... 4 2.3....................... 5 2.3.1...................... 5 2.3.2 ( )....... 6 2.4............................. 9 3 10 3.1............................ 10 3.2.................... 10 3.3.......... 11 3.3.1......... 12 3.3.2...... 12 3.4 ( )............................ 12 3.5............................. 13 4 14 4.1............................ 14 4.2..................... 14 4.3............................. 15 1

5 16 5.1............................ 16 5.2............................ 16 5.3....................... 17 5.4............................. 17 5.5......................... 18 5.5.1 (V = 1, θ = 0)..................... 19 5.5.2 (T = 0.1, θ = 0).................... 20 5.6............................ 22 5.6.1......... 22 5.6.2...... 24 5.6.3 ( )....... 26 5.6.4................. 29 5.7............................ 31 5.7.1 (V = 1, θ = 0)..................... 31 5.7.2 (T = 0.1, θ = 0).................... 32 5.8 PSNR.......................... 33 5.8.1 PSNR(Peak Signal to Noise Ratio)......... 33 5.8.2...... 33 5.8.3 34 5.9............................. 35 6 36 6.1.............................. 36 6.2.............................. 36 6.3.............................. 37 38 39 2

A 40 A.1......................... 40 A.1.1 pgm.h......................... 40 A.1.2 pgm.c......................... 40 A.2......................... 43 A.2.1 FFT1.h........................ 43 A.2.2 FFT1.c........................ 43 A.2.3 FFT2.h........................ 46 A.2.4 FFT2.c........................ 46 A.3.......... 47 A.3.1 convertdata.h..................... 47 A.3.2 convertdata.c..................... 48 A.4......................... 49 A.4.1 blur.c......................... 49 A.5......................... 51 A.5.1 inverse.c........................ 51 A.5.2 inverse2........................ 52 A.5.3 wiener.c........................ 54 A.5.4 inverse2 2.c...................... 56 3

1 1.1 (PSF) ( ) PSF ( ) ( ) 1.2 1.1 PSF PSF 2 ( ) C 4

1.3 2 ( ) 3 PSF ( ) 4 5 6 5

2 2.1 ( ) 2.2 ( ) g(x, y) g(x, y) = h(x, y, α, β, γ)f(α, β)dαdβ (2.1) f(α, β) h(,,, ) (PSF) (2.1) 2 g(x, y) = h(x α, y β)f(α, β)dαdβ (2.2) g(x, y) = f(x, y) h(x, y) (2.3) (2.3) 2 G(µ, ν) = F{g(x, y) = F{f(x, y) h(x, y) = H(µ, ν)f (µ, ν) (2.4) 6

µ, ν x, y H(µ, ν) (OTF) n(x, y) N(µ, ν) (2.2) (2.4) g(x, y) = h(x α, y β)f(α, β)dαdβ + n(x, y) (2.5) G(µ, ν) = H(µ, ν)f (µ, ν) + N(µ, ν) (2.6) (2.2) f g g = [H]f + n (2.7) n 2.3 ˆf g [B] b(x, y) B(µ, ν) ˆf(x, y) ˆF (µ, ν) ˆf(x, y) = b(x α, y β)g(α, β)dαdβ (2.8) ˆF (µ, ν) = [H(µ, ν)f (µ, ν) + N(µ, ν)]b(µ, ν) (2.9) (2.8) (2.5) (deconvolution) 2.3.1 B(µ, ν) B(µ, ν) = 1 H(µ, ν) (2.10) 7

ˆF (µ, ν) ˆF (µ, ν) = F (µ, ν) + N(µ, ν) H(µ, ν) (2.11) (2.10) H(µ, ν) (inverse filter) F (µ, ν) H(µ, ν) (2.11) 2 (2.7) [H] [H] 1 ˆf ˆf = f + [H] 1 (2.12) [H] 1 (2.12) 2 (ill-conditioned problem) (ill-posed problem) (reguralization) 2.3.2 ( ) 2 E ɛ 2 = E[{f(x, y) ˆf(x, y) 2 ] (2.13) E[ ] 2 2 (least square filter) (Wiener filter) 8

2 B(µ, ν) = H (µ, ν) H(µ, ν) 2 + S n (µ, ν)/s f (µ, ν) (2.14) S f (µ, ν) S n (µ, ν) R f (x, y) R n (x, y) S f (µ, ν) = S n (µ, ν) = R f (x, y) exp{ j2π(µx + νy)dxdy (2.15) R n (x, y) exp{ j2π(µx + νy)dxdy (2.16) (2.14) B(µ, ν) H (µ, ν) H(µ, ν) 2 + Γ (2.17) Γ (2.7) 2 [B] ˆf ˆf = [B]g (2.18) ɛ 2 2 [B] min E[ɛɛ T ] = min E[(f ˆf)(f ˆf) T ] (2.19) [B] (2.7),(2.18) E[(f [B]g)(f [B]g) T ] = E[{ff T [B]([H]ff T + nf T ) (ff T [H] + fn T )[B] +[B]([H]ff T [H] T + nf[h] T + [H]fn T + nn T )[B] T ] (2.20) E[fn T ] = E[n T f] = 0 (2.21) 9

f n E[ff T ] = [R ff ], E[nn T ] = [R nn ] (2.22) (2.20) E[ɛɛ T ] = [R ff ] 2[B][H][R ff ] + [B][H][R ff ][H][B] T + [B][R nn ][B] T E[ɛɛ T ] [B] (2.23) = 0 (2.24) [B] = [R ff ][H] T ([H][R ff ][H] T + [R nn ]) 1 (2.25) 2 f 2 f [R ff ] = R 00 R 01... R 0(M 1) R 10 R 11.......... R (M 1)0...... R (M 1)(M 1) (2.26) (2.26) R ij N N R ij = E[f i f T j ] (2.27) f i f j i j R ij = R ji (2.26) (2.14) (2.25) 2 (2.13) (2.19) (2.14) H(µ, ν) = 0 B(µ, ν) 10

2.4 ( ) PSF ( ) 11

3 3.1 PSF ( ) 3.2 f(x, y) x y α(t) β(t) T g(x, y) g(x, y) = T 2 T 2 f(x α(t), y β(t))dt (3.1) = = G(µ, ν) ( T ) 2 f(x α(t), y β(t))dt exp{ j2π(µx + νy)dxdy T 2 T 2 T 2 dt f(x α(t), y β(t)) exp{ j2π(µx + νy)dxdy (3.2) u = x α(t), v = y β(t) (3.3) G(µ, ν) = T 2 T 2 dt f(u, v) exp{ j2π(µ(u + α(t)) + ν(v + β(t))dudv 12

= = T 2 T 2 T 2 T 2 dt f(u, v) exp{ j2π(µu + νv) exp{ j2π(µα(t) + νβ(t))dudv exp{ j2π(α(t)µ + β(t)ν)dt F (µ, ν) (3.4) H(µ, ν) = T 2 T 2 exp{ j2π(α(t)µ + β(t)ν)dt (3.5) G(µ, ν) = H(µ, ν)f (µ, ν) (3.6) (3.5) V θ H(µ, ν) = α(t) = V cos θ t, β(t) = V sin θ t (3.7) T 2 T 2 exp{ j2π(v cos θ tµ + V sin θ tνdt = sin πv (µ cos θ + ν sin θ)t πv (µ cos θ + ν sin θ) = T sin c{πv (µ cos θ + ν sin θ)t (3.8) 3.3 (2.10) (3.8) B(µ, ν) B(µ, ν) = = 1 H(µ, ν) 1 T sin c{πv (µ cos θ + ν sin θ)t ˆF (µ, ν) ˆF (µ, ν) = G(µ, ν)b(µ, ν) (3.9) 13

= G(µ, ν) T sin c{πv (µ cos θ + ν sin θ)t (3.10) ˆf (3.8) H(µ, ν) = 0 (2.11) 3.3.1 r B(µ, ν) = 1 H(µ,ν) µ 2 + ν 2 r 2 1 µ 2 + ν 2 > r 2 (3.11) 3.3.2 H(µ, ν) ɛ B(µ, ν) = 1 H(µ,ν) H(µ, ν) > ɛ 0 H(µ, ν) ɛ (3.12) 3.4 ( ) (2.14) (3.8) B(µ, ν) B(µ, ν) = = H (µ, ν) H(µ, ν) 2 + S n (µ, ν)/s f (µ, ν) T sin c{πv (µ cos θ + ν sin θ)t T sin c{πv (µ cos θ + ν sin θ)t 2 + S n (µ, ν)/s f (µ, ν) T sin c{πv (µ cos θ + ν sin θ)t (3.13) T sin c{πv (µ cos θ + ν sin θ)t 2 + Γ 14

ˆF (µ, ν) ˆF (µ, ν) = G(µ, ν)b(µ, ν) = G(µ, ν)t sin c{πv (µ cos θ + ν sin θ)t T sin c{πv (µ cos θ + ν sin θ)t 2 + S n (µ, ν)/s f (µ, ν) G(µ, ν)t sin c{πv (µ cos θ + ν sin θ)t T sin c{πv (µ cos θ + ν sin θ)t 2 + Γ (3.14) ˆf 3.5 PSF ( ) 15

4 4.1 (, ) (2.5) (2.11) 4.2 g(x, y) g(x, y) = F 1 {F{f(x, y)h(µ, ν) (4.1) g(x, y) g(x, y) = g(x, y) = g(x, y) + n(x, y) (4.2) n(x, y) or 0 g(x, y) g(x, y) g(x, y) + 1 ( ) F 1 {F{g(x, y)b(µ, ν) F 1 {F{ g(x, y)b(µ, ν) F 1 {F{g(x, y) + 1B(µ, ν) F 1 {F{g(x, y)b(µ, ν) F 1 {F{ g(x, y)b(µ, ν) F 1 {F{g(x, y) + 1B(µ, ν) F 1 {F{g(x, y)b(µ, ν) f(x, y) 16

F 1 {F{g(x, y) + 1B(µ, ν) or F 1 {F{g(x, y)b(µ, ν) f(x, y) F 1 {F{g(x, y) + 1B(µ, ν) (4.3) ˆf(x, y) = F 1 {F{g(x, y)b(µ, ν) + F 1 {F{g(x, y) + 1B(µ, ν) 2 (4.4) ˆf(x, y) 4.3 17

5 5.1 PSNR(Peak Signal to Noise Ratio) 5.2 CPU : Celeron 2.00GHz Memory : 768MB OS : Windows XP 18

5.3 pgm(portable graymap) ( ) 5.4 256 256 256 5.1: 19

5.5 (3.8) T, V, θ (4.1) PSF 5.2: T = 0.1, V = 1, θ = 0 5.3: T = 0.1, V = 1, θ = 0 5.4: T = 0.1, V = 1, θ = π/2 5.5: T = 0.1, V = 1, θ = π/2 20

5.5.1 (V = 1, θ = 0) 5.6: T = 0.01, V = 1, θ = 0 5.7: T = 0.01, V = 1, θ = 0 5.8: T = 0.1, V = 1, θ = 0 5.9: T = 0.1, V = 1, θ = 0 21

5.10: T = 1, V = 1, θ = 0 5.11: T = 1, V = 1, θ = 0 5.5.2 (T = 0.1, θ = 0) 5.12: T = 0.1, V = 0.1, θ = 0 5.13: T = 0.1, V = 0.1, θ = 0 22

5.14: T = 0.1, V = 1, θ = 0 5.15: T = 0.1, V = 1, θ = 0 5.16: T = 0.1, V = 10, θ = 0 5.17: T = 0.1, V = 10, θ = 0 23

5.6 (3.8) T = 0.1, V = 1, θ = 0 5.6.1 5.18: (r = 5) 5.19: PSF(r = 5) 5.20: (r = 9) 5.21: PSF(r = 9) 24

5.22: (r = 10) 5.23: PSF(r = 10) 5.24: ( ) 5.25: PSF( ) 25

5.6.2 ɛ = 0.004 5.26: (ɛ = 0.010) 5.27: PSF(ɛ = 0.010) 5.28: (ɛ = 0.004) 5.29: PSF(ɛ = 0.004) 26

5.30: (ɛ = 0.001) 5.31: PSF(ɛ = 0.001) 5.32: ( ) 5.33: PSF( ) 27

5.6.3 ( ) Γ = 0.00010 5.34: (Γ = 0.00100) 5.35: PSF(Γ = 0.00100) 5.36: (Γ = 0.00050) 5.37: PSF(Γ = 0.00050) 28

5.38: (Γ = 0.00010) 5.39: PSF(Γ = 0.00010) 5.40: (Γ = 0.00005) 5.41: PSF(Γ = 0.00005) 29

5.42: (Γ = 0.00001) 5.43: PSF(Γ = 0.00001) 30

5.6.4 5.44: (ɛ = 0.010) 5.45: PSF(ɛ = 0.010) 5.46: (ɛ = 0.005) 5.47: PSF(ɛ = 0.005) 31

5.48: (ɛ = 0.004) 5.49: PSF(ɛ = 0.004) 5.50: (ɛ = 0.001) 5.51: PSF(ɛ = 0.001) 32

5.7 5.7.1 (V = 1, θ = 0) 5.52: (T = 0.05) 5.53: (T = 0.05) 5.54: (T = 0.25) 5.55: (T = 0.25) 33

5.7.2 (T = 0.1, θ = 0) 5.56: (V = 0.5) 5.57: (V = 0.5) 5.58: (V = 2.5) 5.59: (V = 2.5) 34

5.8 PSNR 5.8.1 PSNR(Peak Signal to Noise Ratio) PSNR(Peak Signal to Noise Ratio) ( ( ) 2 ) P SNR = 10 log 10 sqrt({f(x, y) ˆf(x, (db) (5.1) y) 2 ) 5.8.2 PSNR inverse2 inverse2 2 ɛ inverse2 inverse2 2 0.001 26.852376 26.853394 0.002 33.251105 33.311163 0.003 33.320486 33.380404 0.004 33.340275 33.398870 0.005 33.339408 33.369434 0.010 33.280122 33.280385 5.1: PSNR(dB) 35

5.8.3 PSNR inverse2 inverse2 2 T inverse2 inverse2 2 0.01 0.05 33.414531 33.443846 0.10 33.340275 33.398870 0.15 33.338891 33.368337 0.20 33.263328 33.292726 0.25 33.144230 33.170711 1.00 5.2: PSNR(dB)(V = 1, θ = 0) V inverse2 inverse2 2 0.1 0.5 33.416530 33.447274 1.0 33.340275 33.398870 1.5 33.342936 33.372222 2.0 33.298598 33.326756 2.5 33.200220 33.227025 10.0 5.3: PSNR(dB)(T = 0.1, θ = 0) 36

5.9 PSNR 37

6 6.1 2 ( ) PSNR 2 ( ) 3 PSF ( ) 2 4 5 6.2 0 sin ( ) 38

PSNR(dB) 6.3 ( ) ( ) 39

1 2 1 4 40

1. N.Mastronardi, P.Lemmerling, and S.Van Huffel Fast structured total least squares algorithm for solving the basic deconvolution problem,siam Journal on Matrix Analysis and Applications,22(2):533-553,(2000) 2. 1998 3. 1978 4. 2002 41

A A.1 A.1.1 pgm.h /********************************************************* /* pgm.h /* pgm.c /********************************************************* /* #define MAX_IMAGESIZE #define MAX_BRIGHTNESS #define GRAYLEVEL #define MAX_FILENAME #define MAX_BUFFERSIZE 1024 /* 255 /* 256 /* (= +1) 256 /* 256 /* /* /* unsigned char image1[max_imagesize][max_imagesize], image2[max_imagesize][max_imagesize]; int x_size1, y_size1, /* image1 x_size2, y_size2; /* image2 /* void load_image_data( ); /* void save_image_data( ); /* A.1.2 pgm.c /************************************************************** /* pgm.c /* pgm /************************************************************** #include <stdio.h> #include <stdlib.h> 42

#include "pgm.h" /************************************************************** /* pgm /* image1[ ][ ] x_size1 y_size1 /************************************************************** void load_image_data( ) { char file_name[max_filename]; /* char buffer[max_buffersize]; /* FILE *fp; /* int max_gray; /* int x, y; /* /* printf("-----------------------------------------------------\n"); printf(" \n"); printf("-----------------------------------------------------\n"); printf(" pgm, \n"); printf(" (*.pgm) : "); scanf("%s",file_name); fp = fopen( file_name, "rb" ); if ( NULL == fp ){ printf(" \n"); exit(1); /* (=P5) fgets( buffer, MAX_BUFFERSIZE, fp ); if ( buffer[0]!= P buffer[1]!= 5 ){ printf(" P5 \n"); exit(1); /* x_size1, y_size1 # x_size1 = 0; y_size1 = 0; while ( x_size1 == 0 y_size1 == 0 ){ fgets( buffer, MAX_BUFFERSIZE, fp ); if ( buffer[0]!= # ){ sscanf( buffer, "%d %d", &x_size1, &y_size1 ); /* max_gray # max_gray = 0; while ( max_gray == 0 ){ fgets( buffer, MAX_BUFFERSIZE, fp ); if ( buffer[0]!= # ){ sscanf( buffer, "%d", &max_gray ); 43

/* printf(" = %d, = %d\n", x_size1, y_size1 ); printf(" = %d\n",max_gray); if ( x_size1 > MAX_IMAGESIZE y_size1 > MAX_IMAGESIZE ){ printf(" %d x %d \n", MAX_IMAGESIZE, MAX_IMAGESIZE); printf(" \n"); exit(1); if ( max_gray!= MAX_BRIGHTNESS ){ printf(" \n"); exit(1); /* for ( y = 0; y < y_size1; y ++ ){ for ( x = 0; x < x_size1; x ++ ){ image1[y][x] = (unsigned char)fgetc( fp ); printf(" \n"); printf("-----------------------------------------------------\n"); fclose(fp); /**************************************************************** /* image2[ ][ ], x_size2, y_size2 /* pgm /**************************************************************** void save_image_data( ) { char file_name[max_filename]; /* FILE *fp; /* int x, y; /* /* printf("-----------------------------------------------------\n"); printf(" pgm \n"); printf("-----------------------------------------------------\n"); printf(" (*.pgm) : "); scanf("%s",file_name); fp = fopen(file_name, "wb"); /* "P5" fputs( "P5\n", fp ); /* # fputs( "# Created by Image Processing\n", fp ); /* fprintf( fp, "%d %d\n", x_size2, y_size2 ); /* fprintf( fp, "%d\n", MAX_BRIGHTNESS ); 44

/* for ( y = 0; y < y_size2; y ++ ){ for ( x = 0; x < x_size2; x ++ ){ fputc( image2[y][x], fp ); printf(" \n"); printf("-----------------------------------------------------\n"); fclose(fp); A.2 A.2.1 FFT1.h /********************************* /* FFT1.h /* FFT1.c /********************************* #define FFT_MAX 1024 #define PI 3.141592653589 /* int calc_power_of_two( int number ); void make_initial_data( double *re_part, double *im_part, int num_of_data, int power ); void FFT1( double *re_part, double *im_part, int num_of_data, int flag ); A.2.2 FFT1.c /************************************** /* FFT1.c /* 1 FFT /************************************** #include "FFT1.h" /************************************** /* number /* 8 --> 3, 16 --> 4, 32 --> 5,... /************************************** int calc_power_of_two( int number ) { int power; power = 0; 45

while ( number!= 1 ){ power ++; number = number >> 1; return power; /* /******************************************************************** /* FFT /* double *re_part, /* *im_part; ( ) /* int num_of_data; /* int power; num_of_data power /******************************************************************** void make_initial_data( double *re_part, double *im_part, int num_of_data, int power ) { int i, j; /* int ptr, offset; /* int new_ptr; /* int dft; /* DFT ( i double re_buf[fft_max], im_buf[fft_max]; /* dft = num_of_data; for ( i = 1; i < power; i ++ ){ /* i new_ptr = 0; offset = 0; while( new_ptr < num_of_data ){ ptr = 0; while( ptr < dft ){ re_buf[new_ptr] = *( re_part + offset + ptr ); im_buf[new_ptr] = *( im_part + offset + ptr ); new_ptr ++; ptr = ptr + 2; if ( ptr == dft ) ptr = 1; offset = offset + dft; /* for ( j = 0; j < num_of_data; j ++ ){ *( re_part + j ) = re_buf[j]; *( im_part + j ) = im_buf[j]; dft = dft / 2; /******************************************************************** 46

/* FFT (flag = 1), FFT (flag = -1) /* double *re_part, /* *im_part; ( ) /* int num_of_data, flag; FFT /******************************************************************** void FFT1( double *re_part, double *im_part, int num_of_data, int flag ) { int i,j,k,power; /* double unit_angle,step_angle; /* int dft,half,num_of_dft; /* DFT DFT int num_out,num_in1,num_in2; /* DFT double re_buf,im_buf,angle; /* /* static double re_part_new[fft_max],im_part_new[fft_max]; /* FFT ( flg = -1 ) num_of_data /* if ( flag == -1){ for ( i = 0; i < num_of_data; i ++ ){ *( re_part + i ) = *( re_part + i ) / num_of_data; *( im_part + i ) = - *( im_part + i ) / num_of_data; /* num_of_data (power) power = calc_power_of_two( num_of_data ); /* make_initial_data( re_part, im_part, num_of_data, power ); /* 2 DFT, 4 DFT, 8 DFT,... unit_angle = 2.0 * PI / num_of_data; dft = 2; /* = DFT for ( i = 0; i < power; i ++ ){ /* DFT DFT /* i=0 -> 2, i=1 -> 4, i=2 -> 8,... num_of_dft = num_of_data / dft; /* DFT step_angle = unit_angle * num_of_dft; half = dft / 2; for ( j = 0; j < num_of_dft; j ++ ){ angle = 0.0; for ( k = 0; k < dft; k ++ ){ /* No.num_in1 No.num_in2 No.num_out /* No.num_in1 num_out = j * dft + k; if ( k < half ){ num_in1 = num_out; num_in2 = num_in1 + half; else { num_in2 = num_out; num_in1 = num_out - half; 47

/* (re_) (im_) re_buf = *( re_part + num_in2 ); im_buf = *( im_part + num_in2 ); re_part_new[num_out] = *( re_part + num_in1 ) + re_buf * cos(angle) + im_buf * sin(angle); im_part_new[num_out] = *( im_part + num_in1 ) + im_buf * cos(angle) - re_buf * sin(angle); /* angle = angle + step_angle; /* for ( j = 0; j < num_of_data; j ++ ){ *( re_part + j ) = re_part_new[j]; *( im_part + j ) = im_part_new[j]; dft = dft * 2; /* DFT /* FFT ( flg = -1 ) if ( flag == -1 ){ for ( j = 0; j < num_of_data; j ++ ){ *( im_part + j ) = - *( im_part + j ); A.2.3 FFT2.h /************************************* /* FFT2.h /* FFT2.c /************************************* /* double data[fft_max][fft_max]; /* double jdata[fft_max][fft_max]; /* int num_of_data; /* /* void FFT2( int flag ); A.2.4 FFT2.c /*************************** /* FFT2.c 48

/* 2 FFT /*************************** #include "FFT1.h" #include "FFT2.h" /********************************************************* /* data, jdata, num_of_data /* flag = 1 : FFT, flag = -1 : FFT /********************************************************* void FFT2( int flag ) { int i, j; /* static double re[fft_max], im[fft_max]; /* for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ re[j] = data[i][j]; im[j] = jdata[i][j]; FFT1( re, im, num_of_data, flag ); for ( j = 0; j < num_of_data; j ++ ){ data[i][j] = re[j]; jdata[i][j] = im[j]; for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ re[j] = data[j][i]; im[j] = jdata[j][i]; FFT1( re, im, num_of_data, flag ); for ( j = 0; j < num_of_data; j ++ ){ data[j][i] = re[j]; jdata[j][i] = im[j]; A.3 A.3.1 convertdata.h /******************************** /* convertdata.h /* convertdata.c 49

/******************************** /* void make_original_data( ); /* FFT void make_output_image( ); /* FFT A.3.2 convertdata.c /***************************************** /* convertdata.c /* <--> FFT /***************************************** #include <stdio.h> #include <stdlib.h> #include "pgm.h" #include "FFT1.h" #include "FFT2.h" /************************************************************* /* image1[y][x] data[y][x], jdata[y][x] /************************************************************* void make_original_data( ) { int i, j; /* if ( x_size1!= y_size1 ){ printf(" \n"); exit(-1); else { printf("\n \n"); num_of_data = x_size1; for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data[i][j] = (double)image1[i][j]; jdata[i][j] = 0.0; /**************************************************************** /* data[y][x], jdata[y][x] image2[y][x] /*, [0,255] /**************************************************************** void make_output_image( ) 50

{ int x, y; /* double max, min; /*, max = 0.0; min = MAX_BRIGHTNESS; printf("\n FFT \n"); x_size2 = num_of_data; y_size2 = num_of_data; for ( y = 0; y < y_size2; y ++ ){ for ( x = 0; x < x_size2; x ++ ){ if ( data[y][x] < 0 ) data[y][x] = 0; if ( data[y][x] < min ) min = data[y][x]; if ( data[y][x] > max ) max = data[y][x]; for ( y = 0; y < y_size2; y ++ ){ for ( x = 0; x < x_size2; x ++ ){ data[y][x] = MAX_BRIGHTNESS/(max-min)*(data[y][x]-min); image2[y][x] = (unsigned char)data[y][x]; A.4 A.4.1 blur.c /******************************************************************** /* blur.c /* /******************************************************************** #include <stdio.h> #include <stdlib.h> #include <math.h> #include "pgm.h" #include "FFT1.h" #include "FFT2.h" #include "convertdata.h" /*********************** /* sinc /*********************** double sinc(double x){ double n; if ( x == 0 ) n = 1; 51

else n = sin(x)/x; return n; /*************************************** /* /*************************************** void blur( ) { int i, j; /* double T; /* double V; /* double x; /* double h[num_of_data][num_of_data]; /* printf(" T :"); scanf("%lf",&t); printf(" V :"); scanf("%lf",&v); printf(" :"); scanf("%lf",&x); for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ /* h[i][j] = T*sinc(PI*V*T*(j*cos(x)+i*sin(x))); printf("\nfft \n"); for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data[i][j] *= h[i][j]; jdata[i][j] *= h[i][j]; /*********************** /* /*********************** main( ) { load_image_data( ); /* make_original_data( ); /* FFT2( 1 ); /* FFT blur( ); /* FFT2( -1 ); /* FFT make_output_image( ); /* 52

save_image_data( ); /* return 0; A.5 A.5.1 inverse.c /*********************************************************** /* inverse.c /* /* /*********************************************************** #include <stdio.h> #include <stdlib.h> #include <math.h> #include "pgm.h" #include "FFT1.h" #include "FFT2.h" #include "convertdata.h" /*********************** /* sinc /*********************** double sinc(double x){ double n; if ( x == 0 ) n = 1; else n = sin(x)/x; return n; /************************************************************** /* /* 1 /************************************************************** void inverse( ) { int i, j; /* double T; /* double V; /* double x; /* double r; /* double b[num_of_data][num_of_data]; /* 53

printf(" T :"); scanf("%lf",&t); printf(" V :"); scanf("%lf",&v); printf(" :"); scanf("%lf",&x); printf(" r :"); scanf("%lf",&r); /* for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ if ( i*i+j*j <= r*r ) b[i][j] = 1/(T*sinc(PI*V*T*(j*cos(x)+i*sin(x)))); else b[i][j] = 1; printf("\nfft \n"); for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data[i][j] *= b[i][j]; jdata[i][j] *= b[i][j]; /*********************** /* /*********************** main( ) { load_image_data( ); /* make_original_data( ); /* FFT2( 1 ); /* FFT inverse( ); /* FFT2( -1 ); /* FFT make_output_image( ); /* save_image_data( ); /* return 0; A.5.2 inverse2 /********************************************************* /* inverse2.c /* 54

/* /********************************************************* #include <stdio.h> #include <stdlib.h> #include <math.h> #include "pgm.h" #include "FFT1.h" #include "FFT2.h" #include "convertdata.h" /*********************** /* sinc /*********************** double sinc(double x){ double n; if ( x == 0 ) n = 1; else n = sin(x)/x; return n; /****************************************************************** /* /* 0 /****************************************************************** void inverse2( ) { int i, j; /* double T; /* double V; /* double x; /* double n; /* double b[num_of_data][num_of_data]; /* printf(" T :"); scanf("%lf",&t); printf(" V :"); scanf("%lf",&v); printf(" :"); scanf("%lf",&x); printf(" n :"); scanf("%lf",&n); /* for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ b[i][j] = T*sinc(PI*V*T*(j*cos(x)+i*sin(x))); if ( b[i][j] <= n && b[i][j] >= -n ) b[i][j] = 0; else b[i][j] = 1/b[i][j]; 55

printf("\nfft \n"); for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data[i][j] *= b[i][j]; jdata[i][j] *= b[i][j]; /*********************** /* /*********************** main( ) { load_image_data( ); /* make_original_data( ); /* FFT2( 1 ); /* FFT inverse2( ); /* FFT2( -1 ); /* FFT make_output_image( ); /* save_image_data( ); /* return 0; A.5.3 wiener.c /************************************************************* /* wiener.c /* /* /************************************************************* #include <stdio.h> #include <stdlib.h> #include <math.h> #include "pgm.h" #include "FFT1.h" #include "FFT2.h" #include "convertdata.h" /*********************** /* sinc /*********************** double sinc(double x){ 56

double n; if ( x == 0 ) n = 1; else n = sin(x)/x; return n; /******************************************* /* /******************************************* void wiener( ) { int i, j; /* double T; /* double V; /* double x; /* double n; /* SN double b[num_of_data][num_of_data]; /* printf(" T :"); scanf("%lf",&t); printf(" V :"); scanf("%lf",&v); printf(" :"); scanf("%lf",&x); printf("sn n :"); scanf("%lf",&n); /* for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ b[i][j] = T*sinc(PI*V*T*(j*cos(x)+i*sin(x))) /((T*sinc(PI*V*T*(j*cos(x)+i*sin(x))))*(T*sinc(PI*V*T*(j*cos(x)+i*sin(x))))+n); printf("\nfft \n"); for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data[i][j] *= b[i][j]; jdata[i][j] *= b[i][j]; /*********************** /* /*********************** main( ) { 57

load_image_data( ); /* make_original_data( ); /* FFT2( 1 ); /* FFT wiener( ); /* FFT2( -1 ); /* FFT make_output_image( ); /* save_image_data( ); /* return 0; A.5.4 inverse2 2.c /********************************************************* /* inverse2_2.c /* /* /* /********************************************************* #include <stdio.h> #include <stdlib.h> #include <math.h> #include "pgm.h" #include "FFT1.h" #include "FFT2.h" double data1[1024][1024]; double data2[1024][1024]; /*********************** /* sinc /*********************** double sinc(double x){ double n; if ( x == 0 ) n = 1; else n = sin(x)/x; return n; /************************************************************* /* image1[y][x] data[y][x], jdata[y][x] /************************************************************* void make_original_data2( ) { int i, j; /* if ( x_size1!= y_size1 ){ 58

printf(" \n"); exit(-1); else { printf("\n \n"); num_of_data = x_size1; for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data[i][j] = (double)image1[i][j]; jdata[i][j] = 0.0; /* image1[][] +1 /* if ( image1[i][j]!= MAX_BRIGHTNESS) data1[i][j] = (double)image1[i][j]+1; else data1[i][j] = (double)image1[i][j]; /****************************************************************** /* /* 0 /****************************************************************** void inverse2( ) { int i, j; /* double T; /* double V; /* double x; /* double n; /* double b[num_of_data][num_of_data]; /* printf(" T :"); scanf("%lf",&t); printf(" V :"); scanf("%lf",&v); printf(" :"); scanf("%lf",&x); printf(" n :"); scanf("%lf",&n); /* for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ b[i][j] = T*sinc(PI*V*T*(j*cos(x)+i*sin(x))); if ( b[i][j] <= n && b[i][j] >= -n ) b[i][j] = 0; else b[i][j] = 1/b[i][j]; 59

printf("\nfft \n"); for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data[i][j] *= b[i][j]; jdata[i][j] *= b[i][j]; /**************************************************************** /* data[y][x], jdata[y][x] image2[y][x] /*, [0,255] /**************************************************************** void make_output_image2( ) { int x, y; /* double max, min; /*, max = 0.0; min = MAX_BRIGHTNESS; printf("\n FFT \n"); x_size2 = num_of_data; y_size2 = num_of_data; for ( y = 0; y < y_size2; y ++ ){ for ( x = 0; x < x_size2; x ++ ){ /* data[y][x] = (data[y][x]+data2[y][x])/2; if ( data[y][x] < 0 ) data[y][x] = 0; if ( data[y][x] < min ) min = data[y][x]; if ( data[y][x] > max ) max = data[y][x]; for ( y = 0; y < y_size2; y ++ ){ for ( x = 0; x < x_size2; x ++ ){ data[y][x] = MAX_BRIGHTNESS/(max-min)*(data[y][x]-min); image2[y][x] = (unsigned char)data[y][x]+0.5; /*********************** /* /*********************** main( ) { int i, j; /* load_image_data( ); /* make_original_data2( ); /* 60

FFT2( 1 ); /* FFT inverse2( ); /* FFT2( -1 ); /* FFT /* for ( i = 0; i < num_of_data; i ++ ){ for ( j = 0; j < num_of_data; j ++ ){ data2[i][j] = data[i][j]; data[i][j] = data1[i][j]; jdata[i][j] = 0.0; FFT2( 1 ); /* FFT inverse2( ); /* FFT2( -1 ); /* FFT make_output_image2( ); /* save_image_data( ); /* return 0; 61