Saturday, December 31, 2016

Bab 3. Pemrosesan Citra Digital Menggunakan MATLAB


3
Akuisisi Citra



3.1 Introduksi


Citra digital diakuisisi oleh kamera menggunakan array sensor foto. Sensor-sensor tersebut terbuat dari semikonduktor. Elemen-elemen detektor foto di dalam array didesain dalam ukuran tertentu untuk menentukan resolusi citra pada kamera. Untuk menangkap suatu citra berwarna, kamera digital harus memiliki suatu array filter warna (color filter array, CFA) untuk memisahkan cahaya yang masuk menjadi komponen spektral merah, biru dan merah. Setiap komponen spektral tersebut memicu array sensor foto untuk menangkap citra. Setiap komponen citra yang ditangkap berukuran sama.



Selanjutnya akan dijelaskan persamaan-persamaan matematik yang mendeskrisikan proses pencuplikan suatu citra kontinyu dan bagaimana mengkonstuksinya dari hasil pencuplikan.



































































































































































































































Gambar 3.3 Illustrasi domain Fourier atas distorsi aliasing akibat pencuplikan suatu citra kontinyu: (a) Transformasi Fourier atas suatu citra lowpass kontinyu; (b) Transformasi Fourier atas suatu citra tercuplik dengan frekuensi-frekuensi pencuplikan sama dengan frekuensi-frekuensi Nyquist; (c) Transformasi Fourier atas suatu cita tercuplik dengan undersampling; (d) Transformasi Fourier atas suatu cita tercuplik dengan oversampling


Solusi: Berikut diberikan skript MATLAB yang dibutuhkan.

Contoh 3-1
Efek undersampling

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
% Contoh 3_1.m
% Mencuplik suatu citra asli untuk
% mempelajari efek undersampling
%
A = imread('lena.png');
[Height,Width,Depth] = size(A);
if Depth == 1
   f = A;
else
   f = A(:,:,1);
end

% downsample dengan suatu faktor M, tidak ada pra-filter
M = 4;
f1 = f(1:M:Height,1:M:Width);
figure,imshow(f1),
title(['Didownsample dengan ' num2str(M) ':Tidak ada pra-filter'])
f1 = imresize(f1,[Height,Width],'bicubic');% rekonstruksi ukuran penuh
figure,imshow(f1)
title('Direkonstruksi dari citra terdownsample: Tidak ada pra-filter')

% downsample dengan suatu faktor M, setelah pra-filter
% menggunakan suatu filter gaussian
f = imfilter(f,fspecial('gaussian',7,2.5),'symmetric','same');
f1 = f(1:M:Height,1:M:Width);
figure,imshow(f1), title(['Didownsample dengan ' num2str(M) ' setelah pra-filter'])
f1 = imresize(f1,[Height,Width],'bicubic');
figure,imshow(f1)
title('Direkonstruksi dari citra terdownsample: dengan pra-filter')











































Gambar 3.4 Illustrasi distorsi aliasing pada suatu citra riil: (a) Citra asli lena.png; (b) Citra asli yang didownsample dengan faktor 4 pada dua dimensi spasial tanpa pra-filter; (c) Citra asli yang didownsample dengan faktor 4 pada dua dimensi spasial dengan pra-filter menggunakan filter lowpass Gaussian 7 x 7; (d) Citra terekonstruksi dari (b); Citra terekonstruksi dari (c).


Pencuplikan Ideal
Divais-divais pencuplikan pada kenyataannya menggunakan pulsa-pulsa persegi-panjang dengan lebar terbatas. Yaitu, fungsi pencuplikan adalah suatu array pulsa-pulsa persegi-panjang, bukan array impuls. Oleh karena itu, fungsi pencuplikan dapat ditulis menjadi

















Contoh 3-2: Perhatikan citra pada Gambar 3.5a, yang merupakan array persegi-panjang berukuran 256 x 256. Citra tersebut dicuplik dengan suatu pulsa persegi-panjang berukuran 32 x 32, seperti yang ditunjukkan pada Gambar 3.5b. Citra tercuplik ditampilkan pada Gambar 3.5c. Terlihat nyata efek pengaburan yang terjadi akibat pencuplikan non-ideal ini. Transformasi Fourier 2D atas citra masukan, pulsa pencuplik, dan citra tercuplik ditampilkan pada Gambar 3.5d-f. Karena efek pengaburan ekivalen dengan efek pemilteran lowpass, maka transformasi Fourier atas citra tercuplik sedikit lebih sempit bila dibandingkan dengan transformasi Fourier atas citra masukan. Transformasi Fourier atas suatu impuls ideal dan transformasi Fourier atas citra yang tercuplik ideal berturut-turut ditampilkan pada Gambar 3.5g,h. Dapat diperhatikan bahwa transformasi Fourier atas citra yang tercuplik ideal sama dengan transformasi Fourier atas citra masukan. Citra terekonstruksi menggunakan suatu filter lowpass ideal ditunjukkan pada Gambar 3.5i. Citra terekonstruksi tersebut identik dengan citra tercuplik karena digunakan pencuplikan dengan suatu pulsa berlebar-terbatas.


Solusi: Berikut skript MATLAB yang diperlukan.

Contoh 3-2
Efek pencuplikan citra dengan pulsa berlebar terbatas

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
% Contoh 3_2.m
% Menunjukkan efek pencuplikan citra dengan
% suatu pulsa berlebar terbatas
% Dua kasus yang mungkin: ’delta’ and ’pulsa’
% delta berkaitan dengan pencuplikan impuls dan pulsa
% berkaitan dengan pencuplikan dengan pulsa berlebar terbatas.
% ukuran array citra adalah NxN piksel
% area citra aktual yang dicuplik adalah MxM
% TxT adalah area pulsa pencuplikan

samplingType = 'pulsa'; % 'delta' atau 'pulsa’
N = 256; % ukuran array is NxN
f = zeros(N,N); % citra yang akan dicuplik
M = 32; % lebar pulsa citra
M1 = N/2 - M + 1; % titik mulai citra
M2 = N/2+M; % titik akhir
f(M1:M2,M1:M2) = 1; % citra uji dengan nilai 1
p = zeros(N,N); % pulsa pencuplik persegi-panjang
if strcmpi(samplingType,'delta')
  p(128,128) = 1;
end

if strcmpi(samplingType,'pulsa')
   T = M/2; % lebar pulsa citra
   T1 = N/2 - T + 1; % titik mulai pulsa pencuplik
   T2 = N/2+T; % titik akhir pulsa pencuplik
   p(T1:T2,T1:T2) = 1; % pulsa pencuplik; lebar = 1/2 dari citra
end

fs = conv2(f,p,'same'); % konvolusi citra dengan pulsa
figure,imshow(f,[]),title('Citra Masukan')
figure,imshow(p,[]), title('Pulsa Pencuplik')
figure,imshow(fs,[]),title('Citra Tercuplik')

figure,imshow(log10(1+5.5*abs(fftshift(fft2(f)))),[])
title('Transformasi Fourier atas citra masukan')

if strcmpi(samplingType,'delta')
   figure,imshow(log10(1+1.05*abs(fftshift(fft2(p)))),[])
   title('Transformasi Fourier atas fungsi Delta')
end

if strcmpi(samplingType,'pulsa')
   figure,imshow(log10(1+5.5*abs(fftshift(fft2(p)))),[])
   title('Transformasi Fourier atas pulsa persegi-panjang')
end

Fs = fftshift(fft2(fs));
figure,imshow(log10(1+5.5*abs(Fs)),[])
title('Transformasi Fourier atas atas citra tercuplik')

% h adalah filter rekonstruksi 2D
%h = (sinc(0.1*pi*(-128:127)))’ * sinc(0.1*pi*(-128:127));
h = (sinc(0.2*pi*(-128:127)))' * sinc(0.2*pi*(-128:127));
H = fft2(h);
figure,imshow(ifftshift(ifft2(fft2(fs).*H,'symmetric')),[])
title('Citra Terekonstruksi')

















































Gambar 3.5 Efek pencuplikan ideal


Pencuplikan Non-Ideal

Sejauh ini, diskusi yang telah dilakukan adalah pencuplikan citra pada suatu grid persegi-panjang, karena merupakan struktur yang paling banyak digunakan pada sistem-sistem penampil seperti TV. Begitu juga dengan kebanyakan kamera digital yang memiliki sensor-sensor yang dibuat dalam array grid persegi-panjang. Akan tetapi, adalah hal yang memingkinkan untuk menggunakan grid non-persegi-panjang, seperti grid pencuplikan heksagonal, untuk mengakuisisi suatu citra digital.

Keuntungan menggunakan grid pencuplikan heksagonal adalah bahwa citra yang diakuisisi memiliki 13.4% data lebih rendah daripada citra yang diakuisisi dengan grid pencuplikan persegi-panjang. Diketahui juga bahwa deteksi tepi lebih effisien pada grid pencuplikan heksagonal. Pencuplikan heksagonal secara luas digunakan dalam visi mesin dan pencitraan biomedik. Meskipun metode-metode kompressi standar, seperti Motion Picture Experts Group (MPEG), menggunakan struktur grid persegi-panjang untuk mengkode citra diam dan citra bergerak, khususnya dalam kompensasi dan estimasi pergerakan (motion estimation and compensation), tetapi akan lebih effisien menggunakan grid heksagonal dalam mendapatkan rasio kompressi yang lebih tinggi.
































Contoh 3-3
Mengkonversi grid persegi-panjang menjadi heksagonal


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
% Contoh3_3.m
% mengkonversi grid persegi-panjang menjadi heksagonal
clear
close all
M = 256; N = 256;
f = uint8(zeros(N,N));
g = uint8(zeros(M+N,M+N));
% menciptakan suatu citra, yang memuat titik-titik hitam-putih
% saling bergantian berukuran R/2xR/2

R = 2; % faktor replikasi piksel
for k = 1:R:M
for l = 1:R:N
for k1 = 1:R/2
for l1 = 1:R/2
f(k+k1-1,l+l1-1) = 255;
end
end
end
end

figure,imshow(f), title('Grid Persegi-panjang')
Fig1 = figure;
imshow(f), title('Grid Persegi-panjang')
zoom(Fig1,4)

% transformasi persegi-panjang menjadi heksagonal
V = [1 1; 1 -1];
for n1 = 1:M
for n2 = 1:N
t1 = V(1,1)*n1 + V(1,2)*n2;
t2 = V(2,1)*n1 + V(2,2)*n2;
g(513-t1,256+t2) = f(n1,n2);
end
end

figure,imshow(g), title('Grid Heksagonal')
g1 = imcrop(g,[130 130 230 230]);
Fig2 = figure;
imshow(g1), title('Grid Heksagonal')
zoom(Fig2,4)








































Gambar 3.6 Mengkonversi grid persegi-panjang menjadi grid heksagonal; (a) Piksel-piksel pada suatu grid persegi-panjang dengan spasi seragam pada dua arah vertikal dan horisontal; (b) Piksel-piksel pada grid heksagonal setelah transformasi; (c) Versi diperbesar dari (a); (d) Versi diperbesar dari (b). Kedua pembesaran menggunakan faktor 4.


Contoh 3-4: Pada contoh ini, citra lena.png akan dikonversi dari grid pencuplikan persegi-panjang menjadi grid pencuplikan heksagonal. Pada Gambar 3.7a ditampilkan citra lena.png yang dipotong menggunakan grid persegi-panjang. Pada Gambar 3.7b ditampilkan citra lena.png yang dipotong menggunakan grid heksagonal. Jika  M dan N adalah baris dan kolom suatu citra, maka ukuran citra tercuplik secara heksagonal adalah (M+N-1) x (M+N)  . Berikut skript MATLAB yang dibutuhkan.


























Gambar 3.7 Mengkonversi suatu citra riil dari grid persegi-panjang menjadi grid heksagonal

Contoh 3-4
Mengkonversi grid persegi-panjang menjadi grid heksagonal


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
% Contoh 3_4.m
% mengkonversi persegi-panjang menjadi grid heksagonal
% menggunakan transformasi linier
% spasi cuplik adalah DX dalam dimensi baris
% dan DY dalam dimensi kolom
% transformasinya adalah [t1 t2]’ = V*[n1 n2]’;
% dimana V = [DX DY; DX -DY];
% Pada contoh ini diasumsikan DX = DY = 1;

clear
close all
%A = imread('kameraman.tif');
A = imread('lena.png');
%{
[M,N,D] = size(A);
N2 = 2*N;
M2 = 2*M;
%}
D = 3;
%M = 216; N = 284; M2 = 2*M; N2 = 2*N;
M = 164; N = 205; M2 = 2*M; N2 = 2*N;

if D > 1
%f = A(:,:,1);
f = A(60:275,302:512,1);
%f = A(182:345,76:280,1);
else
f = A;
end
figure,imshow(f), title('Grid Persegi-Panjang')
g = uint8(zeros(M+N,M+N));% citra ditransformasi dua kali lebih besar dari asli

% mengkonversi persegi-panjang menjadi heksagonal
V = [1 1; 1 -1]; % transformasi persegi-panjang menjadi heksagonal
for n1 = 1:M
for n2 = 1:N
t1 = V(1,1)*n1 + V(1,2)*n2;
t2 = V(2,1)*n1 + V(2,2)*n2;
g(t1,N+t2) = f(n1,n2);
end
end
figure,imshow(g'),title('Grid Heksagonal')
% kembali ke grid persegi-panjang
VI = inv(V);
f1 = uint8(zeros(M+N,M+N));
for t1 = 2:M+N
%for t2 = -N+1:N-1
for t2 = -N+1:M-1
n1 = floor(VI(1,1)*t1 + VI(1,2)*t2);
n2 = floor(VI(2,1)*t1 + VI(2,2)*t2);
%f1(n1+M/2,n2+N/2) = g(N2+1-t1,M+t2);
%f1(n1+round((N-1)/2),n2+round((M-1)/2)) = g(N2+1- t1,M+t2);
f1(n1+round((N-1)/2),n2+round((M-1)/2)) = g(M+N+1- t1,N+t2);
end
end
f1 = f1';
%f1(:,1:N2) = f1(:,N2:-1:1);
f1(:,1:M+N) = f1(:,M+N:-1:1);
f1(1:M+N,:) = f1(M+N:-1:1,:);
figure,imshow(f1(round(N/2)+2:round(N/2)+M+1,round(M/2)+2:round(M/2)+N+1))
title('Grid Persegi-Panjang')


3.3 Kuantisasi Citra

Seperti yang diindikasikan pada Gambar 3.1, cuplik-cuplik citra memiliki nilai-nilai kontinyu. Yaitu, cuplik-cuplik analog dan harus direpresentasi dalam format digital untuk penyimpanan dan transmisi. Karena jumlah bit (dijit biner) untuk merepresentasikan setiap cuplik citra dibatasi---biasanya 8 atau 10 bit, maka cuplik-cuplik analog harus dikuantisasi menjadi sejumlah level yang kemudian dikodekan dalam bilangan biner. Karena proses kuantisasi mengkompressi nilai-nilai analog menjadi nilai-nilai diskrit, maka akan terjadi distorsi. Distorsi tersebut dikenal dengan derau kuantisasi.























































































































































































































Gambar 3.8 Suatu contoh kuantisator seragam: (a) Citra asli; (b) Citra asli yang dikuantisasi menjadi 3 bit per piksel; (c) Citra asli yang dikuantisasi menjadi 5 bit per piksel; (d) SNR akibat kuantisasi vs bpp; (e) Histogram citra masukan; (f) Error akibat kuantisasi: atas: error kuantisasi L=32; bawah: histogram error; (g) Grafik batas-batas keputusan versus level-level rekonstruksi

















Contoh 3-5
Perancangan kuantisator seragam


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
% Contoh 3_5.m
% Perancangan kuantisator seragam
% MenghitComung interval-interval keputusan dan level-level keluaran
% untuk bit kuantisasi ntuktization dari 1 sampai B
% Mengkuantisasi citra masukan menjadi B bit dan
% Menghitung SNR dalam dB
% untuk 1 sampai B bit dan menggrafikkan hasil-hasilnya
%
clear
close all
%
A = imread('lena.png');
%A = imread(’cameraman.tif’);
[Height,Width,Depth] = size(A);

if Depth == 1
f = double(A);
else
f = double(A(:,:,1));
end

if Height>512 || Width > 512
figure,imshow(uint8(f(60:275,302:585)))
[Counts,Bins] = imhist(uint8(f(60:275,302:585)),256);
else
figure,imshow(uint8(f))
[Counts,Bins] = imhist(uint8(f),256);
end

title('Citra dikuantisasi menjadi 8 bit')
% Menghitung histogram citra masukan
figure,plot(Bins,Counts,'k', 'LineWidth', 2)
title('Histogram citra masukan')
xlabel('Nilai piksel')
ylabel('Jumlah')
%
f1 = zeros(size(f));
E = zeros(size(f));
sigmaf = std2(f);
B = 7;
snr = zeros(B,1);

% metode kuantisasi 1 adalah proses menggunakan
% interval-interval keputusan
% dan level-level rekonstruksinya
% metode 2 sama dengan aturan kuantisasi JPEG
MetodeQuantizer = 1;% opsi 1 atau 2
switch MetodeQuantizer

case 1
% Merancang quatizer seragam untuk bit 1 sampai 7
for b = 1:B
L = 2^b;
q = 255/L;
q2 = q/2;
d = linspace(0,255,L+1);
r = zeros(L,1);
for k = 1:L
r(k) = d(k)+q2;
end
for k = 1:L
[x,y] = find(f>d(k) & f<=d(k+1));
for j = 1:length(x)
f1(x(j),y(j)) = round(r(k));
E(x(j),y(j)) = double(f(x(j),y(j)))-r(k);
% termasuk surf plot
end
if b == 5 && k == 32
figure,subplot(2,1,1),plot(E(150,:),'k',
'LineWidth',2)
title(['Error kuantisasi: L = ' num2str(k)])
xlabel('piksel # pada baris 150')
ylabel('Error kuantisasi')
[Counts,Bins]= hist(E(150,:));
subplot(2,1,2),plot(Bins,Counts,'k',
'LineWidth',2)
title('Histogram error kuantisasi')
xlabel('Error bins'), ylabel('Jumlah')
end
clear x y
end
clear E
if b == 3 || b == 5
if Height>512 || Width > 512
figure,imshow(uint8(f1(60:275,302:585)))
else
figure,imshow(uint8(f1))
end
title(['Citra dikuantisasi menjadi ' num2str(b) ' bit'])
end
snr(b) = 20*log10(sigmaf/std2(f-f1));
end

figure,plot(1:B,snr,'k','LineWidth',2)
title('SNR Vs Bit/pixel atas suatu quantizer seragam')
xlabel('# bit')
ylabel('SNR(dB)')
figure,stairs(d(1:L),r,'k','LineWidth',2)
title('Quantizer seragam: interval-interval keputusan & level-level keluaran')
xlabel('Region keputusan')
ylabel('Level-level keluaran')

case 2
for b = 1:B
L = 2^b;
q = 255/L;
f1 = round(floor(f/q + 0.5)*q);
if b == 5
E = f - f1;
figure,subplot(2,1,1),plot(E(150,:),'k','LineWidth',2)
title(['Error kuantisasi: L = ' num2str(L)])
xlabel('piksel # pada baris 150')
ylabel('Error kuantisasi')
[Counts,Bins]= hist(E(150,:));
subplot(2,1,2),plot(Bins,Counts,'k','LineWidth',2)
title('Histogram error kuantisasi')
xlabel('Error bins'), ylabel('Jumlah')
end
if b == 3 || b == 5
if Height>512 || Width > 512
figure,imshow(uint8(f1(60:275,302:585)))
else
figure,imshow(uint8(f1))
end
title(['Citra dikuantisasi menjadi ' num2str(b) ' bit'])
end
snr(b) = 20*log10(sigmaf/std2(f-f1));
end
figure,plot(1:B,snr,'k','LineWidth',2)
title('SNR Vs Bit/piksel atas suatu quantizer seragam')
xlabel('# bit')
ylabel('SNR(dB)')
end































Contoh 3-6
Perancangan kuantisator seragam

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
% Contoh 3_6.m
% Merancang suatu kuantisator tak-seragam untuk suatu citra masukan
% Dimulai dengan suatu himpunan level masukan inisial dan
% secara iteratif membangun himpunan hasil akhir
% dalam sudut pandang MSE antara citra masukan dan
% citra terekonstruksi
clear
close all
A = imread('kameraman.tif');
%A = imread(’barbara.tif’);
[Height,Width,Depth] = size(A);

if Depth == 1
f = double(A);
else
f = double(A(:,:,1));
end

%figure,imshow(A(60:275,302:585,1));
figure,imshow(A(:,:,1));
title('Citra asli 8bpp')
%
sigmaf = std2(f); % deviasi standar citra masukan
B = 7; % bit-bit kuantisasi
snr = zeros(B,1); % derau sinyal-thd-quantisasi dalam dB
err = 0.1e-01; % akurasi untuk memberhentikan iterasi

for b = 5:5
L = 2^b; % jumlah level-level keluaran
%C = linspace(0.5,255,L); % level-level keluaran inisial
C = linspace(2,253,L);
PartitionRegion = zeros(L,Height,Width);
% partisi-partisi yang memuat piksel-piksel masukan
RegionCount = zeros(L,1); % jumlah elemen dalam setiap partisi
RegionSum = zeros(L,1); % jumlah piksel dalam region partisi
Distortion = zeros(L,1); % distorsi minimum dalam suatu region partisi
TotalDistortion = 1.0e+5; % distorsi minimum total
%TD = abs(TotalDistortion-sum(Distortion)); % distorsi total inisial
% ukuran konvergensi
TD = (TotalDistortion-sum(Distortion))/TotalDistortion;
% iterasi mulai
It = 1; % jumlah iterasi
while TD >= err
for r = 1:Height
for c = 1:Width
d1 = abs(f(r,c)-C(1));
J = 1;
for k = 2:L
d2 = abs(f(r,c)-C(k));
if d2<d1
d1 = d2;
J = k;
end
end
RegionCount(J) = RegionCount(J) + 1;
RegionSum(J) = RegionSum(J) + f(r,c);
PartitionRegion(J,r,c) = f(r,c);
Distortion(J) = Distortion(J)+(f(r,c)-
C(J))*(f(r,c)-C(J));
end
end
Index = logical(RegionCount == 0);% periksa region-region kosong
RegionCount(Index) = 1;%menetapkan region-region kosong jadi 1
Distortion = Distortion ./RegionCount;%distorsi minimum rata-rata
%TD = abs(TotalDistortion-sum(Distortion));
% ukuran konvergensi
TD = (TotalDistortion-sum(Distortion))/TotalDistortion;
C = RegionSum ./RegionCount; % level-level keluaran yang telah diperbaiki
TotalDistortion = sum(Distortion);% distorsi minimum total
It = It + 1;
end
clear Index
sprintf('Jumlah iterasi = %d; \t bpp = %d',It-1,b)
%sprintf(’SNR = %4.2f dB’,10*log10(sigmaf*sigmaf/TotalDistortion))
DR = zeros(L+1,1);
DR(L+1) = 255;
f1 = zeros(size(f));
for k = 2:L
DR(k) = (C(k-1)+C(k))/2;
end

% Grafik level-level keluaran terhadap region-region keputusan
if b == 5
figure,stairs(DR(1:L),C,'k','LineWidth',2)
title(['Kuantisator tak-seragam: ' num2str(b) ' bpp'])
xlabel('Region-region keputusan')
ylabel('Level level keluaran')
figure,subplot(2,1,1),stem(RegionCount/(It-1),'k')
title(['Perbandinganran & histogram masukan pada
' num2str(b) ' bpp'])

[Counts,Bins] = imhist(uint8(f),32);
subplot(2,1,2),stem(Counts,'k')
xlabel('Jumlah level'), ylabel('Jumlah')
end

% dekode menggunakan buku-kode yang dibangkitkan
for r = 1:Height
for c = 1:Width
d1 = abs(f(r,c)-C(1));
for k = 2:L
d2 = abs(f(r,c)-C(k));
if d2<d1
d1 = d2;
J = k;
end
end
f1(r,c) = C(J);
end
end
% dekode menggunakan region-region keputusan dan buku-kode
%{
for k = 1:L
[x,y] = find(f>DR(k) & f<=DR(k+1));
if ?isempty(x)
for j = 1:length(x)
f1(x(j),y(j)) = C(k);
end
end
clear x y
end
%}
if b == 3 || b == 5
%figure,imshow(uint8(f1(60:275,302:585figure,imshow(uint8(f1))
title(['Citra terkuantisasi-ulang pada ' num2str(b) ' bpp'])
end
mse = sum(sum((f-f1).*(f-f1)))/(Height*Width);
snr(b) = 10*log10(sigmaf*sigmaf/mse);
end
figure,plot(1:B,snr,'k','LineWidth',2)
title('SNR Vs Bit/piksel atas suatu kuantisator tak-seragam')
xlabel('# bit')
ylabel('SNR(dB)')










































Gambar 3.9 Suatu contoh kuantisator tak-seragam: (a) Citra asli; (b) Citra terkuantisasi 5 bpp; (c) Atas: Jumlah piksel terkuantisasi pada tiap level-level keluaran yang berbeda; Bawah: Histogram citra masukan; (d) Grafik region-region keputusan masukan versus level-level keluaran suatu kuantisator tak-seragam