Tài liệu Giáo trình Xử lý ảnh - Chương 13: Nén dữ liệu ảnh: Chương
13
Nén dữ liệu ảnh
13.1 Chỉ dẫn
Nén ảnh là một kỹ thuật mã hoá hiệu suất cao ảnh số nhằm làm giảm số bit cần cho biểu diễn ảnh. Chức năng của kỹ thuật này là giảm độ lớn dữ liệu phải lưu trữ cùng với thời gian truyền trong khi vẫn giữ nguyên chất lượng của ảnh. Để đánh giá sự cần thiết của nén ảnh, chúng ta xem xét về yêu cầu bộ nhớ và thời gian truyền khi dùng một modem 9600 baud (bit/s) cho các ảnh sau đây:
■ Một ảnh 512 ´ 512 điểm, 8 bit cho một điểm, ảnh mức xám yêu cầu 2,097,152 bit cho lưu giữ và mất 3.64 phút để truyền.
■ Một ảnh màu RGB có cùng các bước xử lý như trường hợp trên yêu cầu xấp xỉ 6 triệu bít cho lưu trữ và mất gần 11 phút để truyền.
■ Một phim âm bản có kích thước 24 ´ 36 mm (35 mm) chia bằng các khoảng cách nhau 12 àm, vào khoảng 3000 ´ 2000 điểm, 8 bit cho một điểm, yêu cầu 48 triệu bit cho lưu giữ ảnh và 83 phút để truyền. Một phim âm bản màu sẽ yêu cầu một số lớn gấp ba lần cho lưu giữ và truyền.
Rõ ràng, việc truyền và lưu giữ các ảnh sẽ có...
85 trang |
Chia sẻ: hunglv | Lượt xem: 1754 | Lượt tải: 1
Bạn đang xem trước 20 trang mẫu tài liệu Giáo trình Xử lý ảnh - Chương 13: Nén dữ liệu ảnh, để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên
Chương
13
Nén dữ liệu ảnh
13.1 Chỉ dẫn
Nén ảnh là một kỹ thuật mã hoá hiệu suất cao ảnh số nhằm làm giảm số bit cần cho biểu diễn ảnh. Chức năng của kỹ thuật này là giảm độ lớn dữ liệu phải lưu trữ cùng với thời gian truyền trong khi vẫn giữ nguyên chất lượng của ảnh. Để đánh giá sự cần thiết của nén ảnh, chúng ta xem xét về yêu cầu bộ nhớ và thời gian truyền khi dùng một modem 9600 baud (bit/s) cho các ảnh sau đây:
■ Một ảnh 512 ´ 512 điểm, 8 bit cho một điểm, ảnh mức xám yêu cầu 2,097,152 bit cho lưu giữ và mất 3.64 phút để truyền.
■ Một ảnh màu RGB có cùng các bước xử lý như trường hợp trên yêu cầu xấp xỉ 6 triệu bít cho lưu trữ và mất gần 11 phút để truyền.
■ Một phim âm bản có kích thước 24 ´ 36 mm (35 mm) chia bằng các khoảng cách nhau 12 àm, vào khoảng 3000 ´ 2000 điểm, 8 bit cho một điểm, yêu cầu 48 triệu bit cho lưu giữ ảnh và 83 phút để truyền. Một phim âm bản màu sẽ yêu cầu một số lớn gấp ba lần cho lưu giữ và truyền.
Rõ ràng, việc truyền và lưu giữ các ảnh sẽ có nhiều vấn đề. Có rất nhiều ví dụ khác mà sẽ dễ dàng làm sáng tỏ vai trò của nén ảnh, và rất có nhiều nghiên cứu tập trung vào lĩnh vực này. Fax, một tài liệu đồ hoạ được truyền qua đường dây điện thoại, nén dữ liệu ảnh y học, truyền hình là một vài trong số nhiều ứng dụng tiềm tàng của nén ảnh. Sự phát triển của kỹ thuật vi điện tử và sự phát triển của rất nhiều ứng dụng thương mại dẫn dắt sự phát triển cho các tiêu chuẩn và phần cứng của bộ nén ảnh theo thời gian thực.
Nén ảnh là đạt được bâừng cách loại bỏ các phần thừa trong ảnh. Các phần thừa này có thể ở trong miền không gian, miền phổ, hoặc là phần thừa trong thời gian.
■ Phần thừa không gian là kết quả do mối quan hệ tương quan giữa các điểm gần nhau.
■ Phần thừa phổ là kết quả do mối tương quan giữa các mặt phẳng màu khác nhau.
■ Phần thừa thời gian là kết quả mối tương quan của các khung khác nhau một dãy các ảnh như trong truyền hình .
Trong chương này tôi sẽ trình bày với các bạn một số thuật toán nén ảnh có kết quả tốt được thừa nhận rộng rãi. Chúng ta sẽ phát triển thuật toán thành mã chương trình C, sau đó kiểm tra kết quả của các kỹ thuật này qua các ví dụ chạy thử. Bạn sẽ có nhiều kinh nghiệm bằng cách tự chạy các chương trình này.
13.2 Mã thống kê
Những ảnh mà chúng ta thu nhận được mã hoá và lưu giữ dưới dạng "mã tự nhiên". Một mức xám của giá trị được mã hoá bằng 8 bit nhị phân bằng nhau. Ví dụ một mức xám giá trị 6 được mã hoá là 0000 0110. Một sự sắp xếp mã hoá luân phiên nhau được dùng trong mã mức xám. Loại mã này có đặc tính là bất kỳ hai từ mã liền nhau nào cũng chỉ khác nhau một vị trí. Bảng 13.1 trình bày hai kiểu mã khác nhau cho một tín hiệu mẫu có giá trị vào khoảng từ 0 đến 7. Một kiểu cho ta thấy rằng tín hiệu không nhất thiết phải có giá trị thực sự từ 0 đến 7, nhưng phải có 8 mức riêng biệt.
Bảng 13.1 Các mã khoảng cách bằng nhau
Vào
Mã tự nhiên
Mã mức xám
0
1
2
3
4
5
6
7
000
001
010
011
100
101
110
111
111
110
100
101
001
000
010
011
Những loại mã này thường gọi là mã khoảng cách bằng nhau. Mã khoảng cách bằng nhau không được dùng trong trong thống kê dữ liệu. Sự thừa nhận này được tạo ra khi ta coi rằng tất cả các mức xám (hoặc giá trị tín hiệu chói) có cùng số lần xuất hiện trong ảnh. Nếu điều này không đúng, dạng mã này không phải tốt nhất. Nếu chúng ta phát triển một mã mà một số ít bít hơn được kí hiệu cho các từ mã biểu diễn các mức xám có khả năng xuất hiện cao hơn, thì trung bình độ dài từ mã sẽ nhỏ nhất và loại mã mà chúng ta vừa phát triển là cơ bản cho mã phần thừa tối thiểu. Tất cả các loại mã này được biết với tên mã có độ dài thay đổi hoặc đôi khi gọi là mã entropy. Câu hỏi đặt ra cho chúng ta lúc này là :
■ Chiều dài từ mã trung bình tối thiểu mà có thể dùng giải mã để sửa lại mã chính xác là gì?
■ Làm cách nào chúng ta tạo ra mã này?
Câu trả lời cho câu hỏi đầu tiên có thể tìm thấy trong lý thuyết thông tin. Nếu ta cho rằng một mức xám g của xác suất p(g) được cho bằng từ mã dài L(g) bit. Chiều dài từ mã trung bình, trong một ảnh mức xám 8 bit, được cho bởi
bit/ pixel (13.1)
Một thừa nhận hợp lý nữa có thể suy ra là sự kiện có số lần xuất hiện ít, thì sẽ cung cấp nhiều thông tin hơn sự kiện số lần xuất hiện nhiều hơn. Sự thừa nhận này dẫn chúng ta đến mối quan hệ
(13.2)
Cơ số 2 dùng khi L(g) được cho dưới dạng đơn vị nhị phân hoặc bit. Chiều dài từ nhỏ nhất mà có thể được dùng cho bởi
bít/pixel (13.3)
Biểu thức này gọi là entropy của tín hiệu. Entropy thì không bao giờ âm vì p(g) nằm trong khoảng [0,1]. Đạo hàm của biểu thức entropy có thể tìm thấy trong các sách nói về tin học hoặc thông tin. Chú ý rằng cho một ảnh 256 mức xám mà tất cả các mức có khả năng xuất hiện bằng nhau khi dùng biểu thức (13.3) chúng ta có:
bit/pixel
Điều này có nghĩa là một mã có độ dài bằng nhau có thể dùng trên một ảnh mà có hàm phân bố cường độ sáng đồng đều.
Câu trả lời cho câu hỏi thứ hai đề cập đến mã có phần thừa nhỏ nhất (mã tối ưu) được Huffman tìm ra. Loại mã này gọi là mã Huffman và được áp dụng rộng rãi trong các kỹ thuật mã hoá bằng phần cứng cũng như bằng phần mềm trong các ứng dụng thương mại. Bây giờ chúng ta sẽ xem xét sơ đồ mã hoá Huffman.
Thuật toán mã hoá Huffman tuân theo các giới hạn sau:
Không có hai thông báo nào có sự sắp xếp của từ mã giống nhau.
Từ mã của thông báo được mã hóa theo cách mà không cần một sự chỉ dẫn nào thêm để chỉ ra đâu là nơi bắt đầu và đâu là nơi kết thúc của từ mã.
Hạn chế thứ hai chỉ ra rằng không có thông báo nào được mã hoá theo cách mà khi từ mã xuất hiện, bit nối bit, như là một phần của từ mã lớn hơn. Ví dụ, 01, 102, và 202 là các từ mã hợp lệ. Một dãy của các từ mã xuất hiện có dạng 1111022020101111102 có thể tách ra thành 111-102-202-01-01-111-102. Tất cả các vấn đề mà chúng ta cần quan tâm khi giải mã là bộ mã gốc. Nếu như một bộ mã bao gồm 11, 111, 102, 02 thì khi một thông báo bắt đầu vói 11, ta sẽ không biết liệu đây là thông báo 11 hay đây là phần bắt đầu của thông báo 111. Nếu một thông báo 11102 xuất hiện thì ta sẽ không biết liệu đây là 111-02 hoặc là 11-102 được truyền đi.
Mã Huffman được mã hoá theo hai hạn chế trên đây và gọi là mã có độ dư thừa tối thiểu hay gọi là mã tối ưu. Phương pháp mã hoá này theo hai bước: bước thu gọn và bước mở rộng. Để xem xét phương pháp mã hoá này ta coi rằng các thông báo để xây dựng từ mã được sắp xếp theo thứ tự xác suất xuất hiện giảm dần.
p(0) ³ p(1) ³ p(2) ³ ... ³ p(N - 1)
ở đây N là số của các thông báo. Như tôi đã chỉ ra ban đầu, cho bộ mã hoá tối ưu thì độ dài của từ mã được xắp xếp theo thứ tự
L(0) Ê L(1) Ê L(2) Ê ... Ê L(N - 1)
Các bước dưới đây trình bày giải thuật mã hoá Huffman. Giải thuật này, cũng như phần lớn các giải thuật khác trong cuốn sách này, được phát triển bởi chính tác giả.
13.2.1 Giải thuật thu gọn
Các bước của giải thuật thu gọn được trình bày tốt nhất theo các bước sau đây:
Đặt M = N và coi đây là một mảng tuyến tính có kích thước N - 2.
Cho i = 0 đến N - 3 lặp lại các bước sau:
{
Cộng p(M) và p(M - 1). Thay p(M - 1) bằng kết quả thu được.
Xác định vị trí, loc, trong M - 1 vị trí đầu tiên của mảng p trong đó p(M - 1) > p(loc).
Đặt temp = p(M - 1).
Chuyển các giá trị từ p(loc) đến p(M - 2) xuống một vị trí.
Đặt giá trị trung gian vào loc.
Lưu giá trị trung gian trong mảng tuyến tính v theo
v(N - 2 - i) = loc
Giảm M đi một giá trị.
}
Để hiểu giải thuật thu gọn ta xem xét ví dụ sau. Giả sử rằng khả năng xuất hiện của các thông tin là:
p = {0.25, 0.25, 0.125, 0.125, 0.0625, 0.0625, 0.0625, 0.0625}
Giải thuật thu gọn được trình bày ở trên, cho trường hợp này, trong hình 13.1, từ đó chúng ta có thể viết:
0.25
0.25
0.25
0.25
0.125
0.125
0.125
0.125
0.0625
0.0625
0.125
0.125
0.125
0.125
0.25
0.125
0.125
0.125
0.5
0.25
0.25
0.25
0.25
0.25
8
7
6
5
4
3
2
0.5
0.5
1
0.25
0.25
0.25
0.25
0.25
0.0625
0.0625
0.0625
0.0625
Hình 13.1 Giải thuật thu gọn.
Mảng tuyến tính v được dùng để xây dựng mã hoá Huffman theo bước mở rộng được trình bày ở dưới đây.
13.3.2 Bước mở rộng
Bước giải thuật này có thể trình bày theo các bước sau:
1) Gán các giá trị ban đầu cho một mảng H theo:
Chú ý là phần tử thứ hai bằng 1 còn các phần tử khác bằng 0.
Véc tơ tuyến tính
v = [v[1] v[2] v[3] ... v[N - 2]]T
được tính trong bước thu gọn. Chú ý là vị trí của các phần tử được cho bằng 1,2,3 ...
2) Đặt M = 2 và i = 1.
3) Lưu phần tử ở vị trí v[i] trong mảng H vào một biến trung gian temp.
temp = H[v(i)]
4) Dịch chuyển các phần tử các phần tử ở vị trí dưới H[v(i)] một vị trí lên phía trên. Chèn giá trị 0 vào vị trí cuối cùng.
5) Sao chép temp tới vị trí thứ M trong H.
6) Dịch trái H(M) một bít.
1
0
0
1
0
1
0
1
0
1
0
1
0
2
1
1
0
1
1
1
1
1
1
1
1
1
1
3
1
1
0
0
0
0
0
1
0
0
1
0
0
1
0
4
0
1
0
1
0
0
1
1
0
1
1
0
1
1
5
0
1
1
0
0
0
0
0
0
0
0
1
0
6
0
0
1
0
0
1
0
0
0
1
1
7
0
0
1
1
0
0
0
0
8
0
0
0
1
Hình 13.2 Giải thuật mở rộng.
7) H[M + 1] = H[M] + 1.
8) Tăng M và i lên một.
9) Nếu i Ê (N - 2) quay lại bước thứ ba.
Nếu không thì đã hoàn thành, và mã nằm trong bảng H.
Các bước trên giải thích qua hình 13.2 dùng ví dụ hình 13.1.
Các bước trên được lập ra bởi Huffman. Lu và Chen đã nhận ra rằng thuật toán của Huffman không phải lúc nào cũng tạo ra một mã có độ dài đơn diệu tăng dần qua ví dụ của họ dưới đây:
Xem xét khả năng xuất hiện thông tin dưới đây:
Theo giải thuật thu gọn ở trên chúng ta có
và tạo bộ mã
Bộ mã trên không thoả mãn điều kiện về chiều dài của từ mã đơn điệu tăng dần; tuy nhiên, bộ mã này có thể sử dụng nếu nó thoả mãn điều kiện có khả năng giải mã được. Lỗi này có thể sửa được bằng một thay đổi nhỏ theo Lu và Chen trong bước thứ ba từ
p(M - 1) > p(loc) thành p(M - 1) ³ p(loc)
Theo phương pháp này cho ta kết quả
và tạo ra một bộ mã
Chú ý là thuật toán Lu và Chen tạo ra bộ mã thì hoàn thiện hơn thuật toán mô tả ở trên.
Chương trình dưới đây tạo ra bộ mã dùng các bước thu gọn và mở rộng ở trên. Chương trình này sử dụng thuật toán Lu và Chen mô tả ở trên. Chương trình cũng tính trung bình độ dài từ. Chương trình không dùng trên ảnh, mà chỉ minh hoạ các bước thực hiện việc sinh mã Huffman. Chú ý điều kiện dùng trong bước 3 của quá trình thu nhỏ được thay bằng điều kiện của Lu và Chen.
Chương trình 13.1 Chương trình ví dụ sinh mã Huffman.
/*Program 13.1 "HUFFMAN.C". Example program for generating the Huffman code.*/
/* Example program to demonstrate the
algorithm for the generation procedure
of the Huffman code. */
#include
#include
/*float p[]={0.2, 0.18, 0.1, 0.1, 0.1, 0.06, 0.06,
0.04, 0.04, 0.04, 0.04, 0.03, 0.01};
int N=13;*/
float p[]={0.5,0.1,0.1,0.1,0.1,0.1};
int N=6;
void main()
{
int i,M,j,loc;
unsigned char v[11],L[13],code[13];
unsigned char ctemp,Ltemp;
float temp,sum,pt[13];
for(j=0;j<N;j++)
pt[i]=p[j];
clrscr();
for(i=0;i<(N-2);i++)
v[i]=0;
/* Contraction. */
M=N;
for(i=0;i<(N-2);i++)
{
p[M-2]=p[M-1]+p[M-2];
loc=M-2;
for(j=0;j<(M-1);j++)
{
if(p[M-2]>=p[j])
{
loc=j;
break;
}
}
temp=p[M-2];
for(j=M-2;j>=loc;j --)
p[j]=p[j-1];
p[loc]=temp;
M--;
v[(N-3)-i]=loc;
}
printf("\n The v vector is given by:\n");
for(j=0; j<(N-2); j++)
printf ("%d\n", v[j]+1);
/*Expansion.*/
for(j=0; j<N; j++)
{
code[j]=0 ;
L[j]=1;
}
printf("\n\n ");
code[0]=0;
code[1]=1;
M=1;
for(i=0; i<(N-2);i++)
{
if(v[i]==M)
{
code[M+1]=(code[M]<<1)+1;
code[M]=(code[M])<<1;
L[M]++;
L[M+1]=L[M];
}
else
{
ctemp=code[v[i]];
Ltemp=L[v[i]];
for(j=v[i];j<M;j++)
{
code[j]=code[j+1];
L[j]=L[j+1];
}
code[M]=ctemp;
L[M]=Ltemp;
code[M+1]=(code[M]<<1)+1;
code[M]=code[M]<<1;
L[M]++;
L[M+1]=L[M];
}
M++;
}
printf("Code words Length\n");
for(j=0;j<N;j++)
printf(" %d %d\n",code[j],L[j]);
sum=0.0;
for(j=0;j<N;j++)
sum+=pt[j]*(float)L[j];
printf("Lav.=%f",sum);
getch();
}
Độ dài từ trung bình được tính bởi công thức (13.1). Chú ý rằng chương trình cũng tính độ dài các bit của mỗi mã và độ dài từ trung bình. Độ dài sẽ mã cần thiết cho việc mã hoá và giải mã các tín hiệu.
Bài tâp 13.1
1. Chạy chương trình 13.1 dùng ví dụ cho trong tài liệu của Huffman:
p = {0.2, 0.18, 0.1, 0.1, 0.1, 0.06, 0.06, 0.04, 0.04, 0.04, 0.04, 0.03, 0.01}
2. Đổi giả thiết:
if (p[M - 1] ³ p[j]
thành
if (p[M - 1] > p[j])
trong chương trình. Chạy chương trình và so sánh với bộ mã cho bởi Lu và Chen.
Dùng ví dụ cho bởi Lu và Chen, chú ý sự khác nhau trong chiều dài từ mã vói hai điều kiện khác nhau.
Nhiệm vụ tiếp theo của chúng ta thật rõ ràng. Chúng ta cần viết một chương trình C để mã hoá và giải mã ảnh số dùng lược đồ mã Huffman. Chương trình sẽ bao gồm: tính toán khả năng xảy ra của mỗi mức xám và xắp xếp lại theo thứ tự giảm dần của khả năng có thể. Tất nhiên, chương trình cũng bao gồm một vector để lưu giữ các thông tin cần thiết về các mức xám theo khả năng của các sự kiện, ví dụ, một bảng tra cứu. Chương trình cũng sẽ loại trừ các sự kiện không có khả năng xảy ra. Khi các bước ở trên được thực hiện, mã Huffman được sinh ra theo giải thuật được mô tả trên. Trong việc mã hoá ảnh số bạn cần nhớ rằng chương trình sẽ cần nhiều byte trên đĩa. Khi viết một thủ tục giải mã ta cần chú ý đến giới hạn này.
13.2.3 Mã hoá ảnh số
Trước khi mã hoá ta cần tạo ra hai bảng tra cứu:
1. Một bảng tra cứu (LUT) chứa quan hệ các giá trị mức xám trên ảnh vói bộ mã hoá Huffman. Bảng LUT này có thể phát triển dùng quan hệ:
table_code[gray[i]] = code[i]
ở đây code[i] là mã Huffman.
2. Một bảng tra cứu xác lập quan hệ giữa các mức xám trên ảnh vói chiều dài từ mã Huffman. LUT có thể tạo ra theo mối quan hệ sau:
table_length [ gray[i]] = L[i]
Mã hoá theo các bước sau:
a. Mở một file để chứa ảnh đã được mã hoá.
b. Đặt Len = 0 và aux = 0.
aux là một thanh ghi bốn byte để truyền mã.
Len chứa số các bít còn lại trong aux.
c. Đọc giá trị điểm ảnh.
d. Xác định mã Huffman của nó, code[i], và độ dài của từ mã L[i], dùng hai bảng LUT để tạo ra giá trị ban đầu.
e. Dịch aux sang trái L[i] bit.
f. Thực hiện phép logic OR code[i] với aux.
g. Len = Len + L[i].
h. Nếu Len ³ 8 thì chuyển 8 bit cuối đến file đầu ra và giảm Len đi 8.
i. Nếu chưa hết file thì chuyển đến bước c.
k. Chuyển các bít còn lại trong thanh ghi aux ra file đầu ra.
Để giải mã được ảnh chúng ta cần phải tạo thêm phần header của file. Header của file bao gồm cả chiều dài thực sự của file tính theo bit. Chú ý là chiều dài thực sự có thể lớn hơn hoặc bằng chiều dài thực sự của file. Điều này bởi vì chúng ta chỉ có thể chứa dưới dạng đơn vị byte, còn file ảnh đã mã hoá phải có chiều dài không chia hết cho 8. Phần này chứa đủ thông tin để thiết lập các bảng LUT, dùng cho việc giải mã ảnh. Nó bao gồm các mức xám trên ảnh, dạng mã Huffman và chiều dài tương đương của chúng. Chú ý rằng các từ mã được sắp xếp theo thứ tự giảm dần theo xác suất xuất hiện của chúng.
Chương trình C sau đây sẽ trình bày các bước trên.
Chương trình 13.2 "HUCODING" . Mã hoá Huffman ảnh số.
/*Program 13.2 "HUCODING.C". Huffman coding of digital images.*/
/* This program is for coding a binary file
using Huffman codes. */
#include
#include
#include
#include
#include
#include
void main()
{
int i,j,N,M,loc,k,ind,xt,yt;
unsigned char *v,*L,*gray, *table_length;
unsigned long int *code, *table_code, ctemp2;
unsigned long int aux1,aux2,Lmask,act_len,flength;
unsigned char mask,Len;
int ch;
unsigned char ctemp,Ltemp;
float temp, sum, *pt, *p, big;
unsigned long int *histo;
double nsq;
char file_name1[14], file_name2[14];
FILE *fptr,*fptro;
clrscr();
printf("Enter file name of image -->");
scanf("%s",file_name1);
fptr=fopen(file_name1,"rb");
if(fptr==NULL)
{
printf("%s does not exist. ",file_name1 );
exit(1);
}
printf("\nEnter file name for compressed image -->");
scanf("%s",file_name2);
k=stricmp(file_name1 ,file_name2);
if (k==0)
{
printf("\nInput and output files cannot share the same name.");
fcloseall();
exit(1);
}
ind=access(file_name2,0);
while(!ind)
{
k=stricmp(file_name1, file_name2);
if(k==0)
{
gotoxy(1,8);
printf( "Input and output files cannot share the same name.");
exit(1);
}
gotoxy(1,4);
printf("File exists. Wish to overwrite? (y or n)-->");
while(((ch=getch())!='y')&&(ch!='n'));
putch(ch);
switch(ch)
{
case 'y':
ind=1;
break;
case 'n':
gotoxy (1,4);
printf( " ");
gotoxy (1,3);
printf ( " ");
gotoxy (1,3);
printf ( " Enter file name -->");
scanf ( "%s " , file_name2 ) ;
ind=access ( file_name2, 0);
}
}
fptro=fopen(file_name2, "wb");
xt=wherex();
yt=wherey();
gotoxy(70, 25);
textattr (RED+(LIGHTGRAY<<4) +BLINK);
cputs ( "WAIT" ) ;
/* Generate histogram. */
histo=( unsigned long int *) malloc(256*sizeof(long int));
for(i=0; i<256; i++)
histo[i]=0;
while((ch=getc(fptr))!=EOF)
histo[ch]++;
p=(float *)malloc(256*sizeof(float));
gray=(unsigned char *)malloc(256*sizeof(char));
k=0;
for(i=0;i<256;i++)
{
if(histo[i]!=0)
{
p[k]=(float)histo[j];
gray[k]=i ;
k++;
}
}
free(histo);
N=k;
/* Normalize.*/
sum=0.0;
for(i=0;i<N;i++)
sum+=p[i];
for(i=0;i<N;i++);
p[i]/=sum;
/* Rearranging the probability values in ascending order. */
for(i=0;i<(N-1);i++)
{
big=p[i];
loc=i;
for(j=i+1;j<N;j++)
{
if(p[j]>big)
{
big=p[i];
loc=j;
}
}
if(loc==i) continue;
else
{
temp=p[i];
ctemp=gray[j];
p[i]=p[loc];
gray[i]=gray[loc];
p[loc]=temp;
gray[loc]=ctemp;
}
}
pt=(float *)malloc(N*sizeof(float));
for(j=0;j<N;j++)
pt[j]=p[j];
v=(unsigned char *)malloc((N-2)*sizeof(char));
code=(unsigned long int *)malloc(N*sizeof(unsigned long int));
L=(unsigned char *)malloc(N*sizeof(char));
for(i=0; i<(N-2); i++)
v[i]=0;
/* Contraction steps in the generation of the Huffman's codes. */
M=N ;
for(i=0; i<(N-2) ;i++)
{
p[M-2]=p[M-1]+p[M-2];
loc=M-2;
for(j=0;j<(M-1);j++)
{
if(p[M-2]>=p[j])
{
loc=j;
break;
}
}
temp=p[M-2];
for(j=M-2;j>=loc;j--)
p[j]=p[j-1];
p[loc]=temp;
M--;
v[(N-3)-i]=loc;
}
/*Expansion steps in the generation of the Huffman's codes.*/
for(j=0;j<N;j++)
{
code[j]=0;
L[j]=1;
}
code[0]=0;
code[1]=1 ;
M=1;
for(i=0; i<(N-2);i++)
{
if(v[i]==M)
{
code[M+1]=(code[M]<<1)+1;
code[M]=(code[M])<<1;
L[M]++;
L[M+1]=L[M];
}
else
{
ctemp2=code[v[i]];
Ltemp=L[v[i]];
for(j=v[i];j<M;j++)
{
code[j]=code[j+1];
L[j]=L[j+1];
}
code[M]=ctemp2;
L[M]=Ltemp;
code[M+1]=(code[M]<<1)+1;
code[M]=code[M]<<1;
L[M]++;
L[M+1]=L[M];
}
M++;
}
/*printf("Code words Length\n");
for(j=0;j<N;j++)
printf(" %lu %u\n",code[j],L[j]);*/
sum=0.0;
for(j=0;j<N;j++)
sum+=pt[j]*(float)L[j];
gotoxy(xt,yt);
printf("\nAverage number of bits/pixel.=%f",sum);
free(v);
free(p) ;
free(pt);
/* Coding */
/* Writing the header in the output file.
The first 4 bytes stores the true length in
bits of the stored image with the MSB stored
first. The 5th byte is the number of Huffman
codes=N. The following N bytes contain the N
natural codes for the gray levels, followed by N
bytes containing the Huffman code lengths. These
are then followed by the actual Huffman codes
stored in packed form. */
for(i=0;i<4;i++)
putc((int)0,fptro);
/* reserve the first 4 bytes for the true length of
the file in bits.*/
putc(N,fptro);
for(i=0;i<N;i++)
putc(gray[i],fptro);
for(i=0;i<N;i++)
putc(L[i],fptro);
aux1=0;
Len=0;
for(i=0;i<N;i++)
{
Len+=L[i];
aux1=(aux1<<L[i])|code[i];
if(Len<8) continue;
else
while(Len>=8)
{
aux2=aux1;
Lmask=255;
Lmask<<=(Len-8);
aux2=(aux2&Lmask);
aux2>>=(Len-8);
ch=(char)aux2;
putc(ch,fptro);
Lmask=~Lmask;
aux1=aux1&Lmask;
Len-=8;
}
}
if(aux1!=0)
{
ch=(char)(aux1<<(8-Len));
putc(ch,fptro);
}
table_code=(unsigned long int *)malloc(256*sizeof(long int));
table_length=(unsigned char *)malloc(256*sizeof(char));
for(i=0; i<N; i++)
{
table_code[gray[i]]=code[i];
table_length[gray[i]]=L[i];
}
rewind(fptr);
Len=0 ;
act_len=0; /*variable to save true length of file in bits. */
aux1=aux2=0;
k=0;
while((ch=getc(fptr))!=EOF)
{
k++;
Len+=table_length[ch];
act_len+=table_length[ch];
aux1=(aux1<<table_length[ch])|table_code[ch];
if(Len<8) continue;
else
while(Len>=8)
{
aux2=aux1;
Lmask=255;
Lmask<<=(Len-8);
aux2=aux2&Lmask;
Lmask=~Lmask;
aux1=aux1&Lmask;
aux2>>=(Len-8);
ch=(char)aux2;
Len-=8 ;
putc(ch,fptro);
}
}
if(aux1!=0) /* Transmit remaining bits.*/
{
ch=(char)(aux1<<(8-Len));
putc(ch,fptro);
}
rewind(fptro); /* Write true length of file. */
M=8*sizeof(long int)-8;
for(i=0;i<4;i++)
{
Lmask=255;
Lmask<<=M;
aux1=(act_len&Lmask);
aux1>>=M;
ch=(int)aux1;
putc(ch,fptro);
M-=8;
}
fclose(fptro);
fclose(fptr);
gotoxy(70,25);
textattr(WHITE+(BLACK<<4));
cputs(" ");
gotoxy(1,yt+2);
getch( );
}
Bài tập 13.2
áp dụng chương trình 13.2 cho IKRAM.IMG. Tên file ra là IKRHUFF.IMG. Bạn sẽ cần đến file này để kiểm tra chương trình giải mã, chương trình này sẽ được trình bày sau.
Tính độ dài từ trung bình từ kích thước file ra. So sánh với độ dài từ đã được tính bằng chương trình.
13.2.4 Giải mã
Giải mã của một ảnh mã hoá bằng mã Huffman thực hiện qua các bước sau:
1. Đặt Len = 0, flength = 0; aux = 0.
Len là số đếm của các bít được giải mã.
Flength là một số đếm của các để so sánh với chiều dài thực sự của file.
Aux là một thanh ghi bốn byte chứa các từ mã sẽ được giải mã.
Lặp lại các bước sau đây cho tới khi chiều dài flength > true_length
{
2. Chuyển một byte từ file đến thanh ghi aux.
Lặp lại các bước sau 8 lần:
{
3. Dịch trái thanh ghi aux đi một bít.
4. Dịch chuyển bit dấu lớn nhất của ch đến vị trí bít dấu nhỏ nhất của ch.
5. Dịch trái 1 bit ch.
6. Tăng Len lên 1.
7. Tăng flength lên 1.
8. Kiểm tra lại tất cả các mã Huffman có chiều dài length = Len.
Nếu một mã được tìm ra:
a. Sao chép mức xám tương ứng ra file xuất ra.
b. Đặt aux bằng không.
c. Đặt Len bằng không.
} ( cụ thể chuyển tói bước 3).
} ( cụ thể chuyển tới bước 2).
Các thủ tục trên tạo ra các bước giả mã. Bạn cần nhớ rằng header của file chứa các thông tin cho việc giải mã. Một chương trình C cho việc giải mã được trình bày ở phần dưới đây.
Chương trình 13.3 "HUDECDNG". Chương trình giải mã ảnh được mã hoá Huffman.
/*Program 13.3 "HUDECDNG.C". Program for decoding a Huffman-coded image.*/
/*This program is for decoding binary files
coded in Huffman codes. */
#include
#include
#include
#include
#include
#include
#include
void main()
{
unsigned int i,j,N,M,k,xt,yt;
unsigned char *L,*gray;
unsigned long int *code;
unsigned long int aux1,aux2,Lmask,act_len,flength;
unsigned char mask,Len;
int ch,ind;
unsigned char ctemp,Ltemp;
char file_name1 [14],file_name2[14];
FILE *fptr,*fptro;
clrscr();
printf("Enter input file name for compressed image -->");
scanf("%s", file_name1);
fptr=fopen(file_name1,"rb");
if(fptr==NULL)
{
printf("\n No such file exists.");
exit(1);
}
printf("\nEnter file name for uncompressed image-->");
scanf("%s",file_name2);
k=stricmp(file_name1,file_name2);
if(k==0)
{
printf
("\nInput and output files cannot share the same name.");
exit(1);
}
ind=access(file_name2,0);
while(!ind)
{
k=stricmp(file_name1 ,file_name2);
if(k==0)
{
gotoxy(1,8);
printf("Input and output files cannot share the same name.");
exit(1);
}
gotoxy(1,4);
printf("File exists. Wish to overwrite? (y or n)-->");
while(((ch=getch())!='y')&&(ch!='n'));
putch(ch);
switch(ch)
{
case 'y':
ind=1;
break;
case 'n':
gotoxy(1,4);
printf(" ");
gotoxy(1,3);
printf(" ");
gotoxy(1,3);
printf(" Enter file name -->");
scanf ( "%s", file_name2);
ind=access(file_name2, 0);
}
}
fptro=fopen(file_name2, "wb");
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(RED+(LIGHTGRAY<<4)+BLINK);
cputs("WAIT");
/*Reading the header.*/
act_len=0;
M=8*sizeof(long int)-8;
for(j=0;i<4;i++)
{
ch=getc(fptr);
aux1=(unsigned long int)ch;
aux1<<=M;
act_len|=aux1;
aux1=(long int)0;
M-=8;
}
N=getc(fptr);
gray=( unsigned char * )malloc(N*sizeof(char));
L=(unsigned char *)malloc(N*sizeof(char));
code=( unsigned long int *)malloc ( N*sizeof(long int));
for(i=0; i<N; i++)
gray[i]=getc(fptr);
for(i=0; i<N ; i++)
L[i]=getc(fptr);
mask=128;
Len=0;
aux2=0;
aux1=0;
ind=1;
i=0;
while(ind)
{
ch=getc(fptr);
for(k=0;k<8;k++)
{
aux1=ch&mask;
ch<<=1;
aux2<<=1;
if(aux1!=0)
aux2|=1;
Len++;
if(Len==L[i])
{
code[i]=aux2;
aux2=0;
Len=0;
i++;
if(i==N){
ind=0;
break;
}
}
}
}
Len=0;
aux1=aux2=0;
flength=0; /*a parameter to measure against actual file length.*/
mask=128;
ind=1;
while(ind)
{
ch=getc(fptr);
for(k=0;k<8;k++)
{
/* test if the actual end of file has been reached.*/
if(flength>act_len){ind=0; break;}
aux1=(ch&mask);
ch<<=1;
aux2<<=1;
if(aux1!=0)
aux2|=1;
Len++;
for(j=0;j<N;j++)
{
if(L[j]>Len) break;
if(L[j]==Len)
{
if(code[j]==aux2)
{
aux2=0;
putc((int)gray[j],fptro);
Len=0 ;
break;
}
}
}
flength++;
}
}
gotoxy(70, 25);
textattr(WHITE+(BLACK<<4));
cputs( " ");
gotoxy (xt, yt);
printf(" \ n Done decoding. ");
fcloseall();
getch();
}
0
1
1
1
1
2
3
4
5
6
7
8
0
0
0
0
0
0
1
1
Nút gốc
Nút lá
Các nút gốc
Hình 13.3 Cây nhị phân giải mã Huffman hình 13.2.
Bài tập 13.3
áp dụng chương trình 13.3 cho việc giải mã ảnh đã mã hoá ở bài tập 13.2.
Mã hoá Huffman có thể đạt được kết quả hơn nhờ sử dụng cây nhị phân. Cây nhị phân trong hình 13.3 biểu diễn mã Huffman ở hình 13.2.
Viết chương trình C sử dụng mã Huffman đặt trong phần header của file ảnh đã mã hoá để tạo một cây nhị phân.
Mở rộng cho chương trình giải mã dùng cây nhị phân. Chương trình phải đạt được một vài yêu cầu quan trọng nhanh hơn phương pháp đã mô tả trong phần này. giải thích tại sao.
13.3 Mã chiều dài thay đổi
Mã chiều dài thay đổi (RLC) là một phương pháp nén ảnh dựa trên sự cắt bớt các dư thừa không gian. Cho mã hoá chiều dài thay đổi một chiều, một mã chiều dài thay đổi được định nghĩa là một số các phần tử điểm ảnh liên tục có chung một giá trị. Một ảnh có thể mã hoá dùng một cặp (mã chiều dài thay đổi, mã mức xám). Một chương trình như vậy sẽ không thể làm giảm kích thước của ảnh nếu ảnh không chứa các điểm có cùng các giá trị mức xám. Điều kiện này xuất hiện trong một ảnh nhiều chi tiết. Dù có thế đi chăng nữa thì định nghĩa của RLC có thể là một phương pháp tốt để mã hoá mà có thể khắc phục các vấn đề xuất hiện dựa theo các điều kiện sau:
Một mã chiều dài thay đổi được xác định bằng ba bít cuối cùng có ý nghĩa của nó được xác lập bằng 1. Còn 5 bít thấp của nó cung cấp một bộ đếm từ 1 đến 31 cho byte đi theo nó.
Nếu giá trị một điểm có mã chiều dài thay đổi bằng không, thì nó được mã hoá như sau:
Nếu 3 bit cuối của nó đều xác lập lên 1, châửng hạn, ³ 224, thì nó được mã hoá thành (11100000, giá trị điểm), cụ thể, mã chiều dài thay đổi bằng không theo sau bằng giá trị điểm.
Cho các trường hợp còn lại, nó được mã hoá như giá trị điểm.
Các bước trên giả thiết rằng trong một ảnh bình thường, mã chiều dài thay đổi lớn hơn 31 ít xuất hiện, và các điểm có giá trị lớn hơn 224 cũng ít xuất hiện. Chương trình C sau sẽ thực hiện các bước trên.
Chương trình 13.4 "RLC.C". Chương trình cho giải thuật RLC.
/*Program 13.4 "RLC.C". Program for RLC.*/
/* Run length code. */
/* This program can be used for either
coding images in RLC or decoding them. */
#include
#include
#include
#include
#define MASK 224 /* MASK=11100000 */
void main()
{
int buff[256],p;
int N,i,k;
FILE *fptr,*fptro;
char file_name[14],ch;
int NMASK=(~MASK);
clrscr();
printf("This program can be used for both coding");
printf("\nand decoding images using 1-D RLC.");
printf("\n\n Enter choice:");
printf("\n 1.Coding.");
printf("\n 2.Decoding --- >" );
while(((ch=getch())!='1')&&(ch!='2'));
putch(ch);
printf("\n\nEnter input file name-->");
scanf("%s",file_name);
fptr=fopen(file_name,"rb");
if(fptr==NULL)
{
printf("\n file %s does not exist.\n");
exit(1);
}
switch(ch)
{
case '1':
printf("Enter output file name-->");
scanf("%s",file_name);
fptro=fopen(file_name,"wb");
/* Read first line.*/
N=0;
while((buff[N]=getc(fptr))!=EOF)
{
N++;
if(N==256) break;
}
while(N!=0)
{
i=0;
while(i<N)
{
p=buff[i]; /*Read reference value.*/
k=0;
if(i==(N-1)) /*Is reference value last value in line? */
{
if((p&MASK)==MASK)
{
k=k|MASK;
putc(k,fptro);
putc(p,fptro);
}
else
putc(p,fptro);
break;
}
i++;
while(p==buff[i])
{
i++;
k++;
if((k==31)||(i==N)) break;
}
if(k==0)
{
if((p&MASK)==MASK)
{
k=MASK;
putc(k,fptro);
putc(p,fptro);
}
else
putc(p,fptro);
}
else
{
k|=MASK;
putc(k,fptro);
putc(p,fptro);
}
}
N=0;
while((buff[N]=getc(fptr))!=EOF)
{
N++;
if(N==256) break;
}
}
fcloseall();
printf("\nDone coding");
break;
case '2': /* Decoding procedure.*/
printf("\nEnter file name for storing decoded image-->");
scanf("%s",file_name);
fptro=fopen(file_name,"wb");
while(1)
{
k=getc(fptr);
if(k==EOF) break;
if((k&MASK)==MASK)
{
p=getc(fptr);
if(p==EOF)
{
putc(k,fptro);
break;
}
k=k&(NMASK);
for(i=0;i<=k;i++)
putc(p,fptro);
}
else
putc(k,fptro);
}
fcloseall();
break;
}
}
Bài tập 13.4
1. Dùng chương trình 13.4 trên ảnh IKRAM.IMG.
2. Dùng mã Hufman để nén ảnh và so sánh với các phương pháp mã hoá khác.
3. áp dụng các bước giải mã Huffman và RLC để giải nén ảnh.
4. Viết một chương trình C dùng các khái niệm cơ bản (mã chiều dài thay đổi, mã) mà không sử dụng các bước mô tả cho chương trình 13.4. áp dụng chương trình của bạn lên ảnh IKRAM.IMG. Khi đó kích thước của ảnh sau khi nén so với ảnh gốc sẽ thế nào ?
5. áp dụng mã Huffman với:
a. File rút ra từ phần 4.
b. Chiều dài thay đổi và giá trị điểm chia bằng cách xem xét chúng như hai file.
Có thể chuyển mã 1-D RLC sang mã 2-D RLC bằng cách kiểm tra các dòng trước, hoặc kiểm tra bốn hướng khác nhau (trên, dưới, trái, phải). Các 2-D RLC này có thể nén ảnh ở mức độ cao hơn.
13.4 Mã chuyển đổi
Nhắc lại là biến đổi Fourier cho một ảnh thì có phần lớn các giá trị lớn nhất nằm ở miền tần số thấp. Mật độ các giá trị này giảm xuống nhanh chóng khi tần số tăng lên. Tính chất này, tính chất mà chúng ta áp dụng để lọc ảnh, cũng được áp dụng trong khi nén ảnh. Có một số phép biến đổi thuận tiện hơn phép biến đổi Fourier. Phép biến đổi tối ưu nhất là phép biến đổi được đề xuất bởi Karhunen-Loeve (KL). Tuy nhiên, phép biến đổi này tự nó không thể đưa ra các bước tính toán nhanh hoặc hiệu quả hơn. Một phép biến đổi xem có vẻ giống như biến đổi KL nhưng có thể tính toán như biến đổi Fourier rời rạc là phép biến đổi cosin rời rạc. Biến đổi cosin có một sự thay đổi nhỏ tối ưu hoá trong miền tập trung năng lượng so với biến đổi KL; nhưng do ưu điểm của kỹ thuật tính toán nên nó được áp dụng như một tiêu chuẩn trong kỹ thuật nén ảnh.
Phép biến đổi được áp dụng trên toàn bộ ảnh nhưng thông thường người ta hay áp dụng trên các khối nhỏ hơn có kích thước 8 ´ 8 hoặc 16 ´ 16. Lý do là:
Biến đổi của các khối nhỏ thì dễ tính hơn là biến đổi cho toàn bộ ảnh.
Quan hệ giữa các điểm ảnh ít thay đổi giữa các điểm ảnh gần nhau.
Chúng tôi sẽ trình bày dưới đây phép biến đổi cosin và các kỹ thuật tính toán hoàn thiện của nó. Một số phép biến đổi khác như Hadamard, Walsh, ..., không được nghiên cứu trong cuốn sách này, bởi vì chúng không tối ưu bằng phép biến đổi Fourier và chúng có nhiều giới hạn trong lĩnh vực này.
13.4.1 Biến đổi cosin
Biến đổi một chiều cosin rời rạc (DCT-Discrete Cosin Transform) cho bởi
(13.4)
ở đây
với k = 0
= 1 với các trường hợp còn lại.
và k = 0,1,2, ..., N.
Một số phương pháp biến đổi nhanh cosin (FCT-Fast Cosine Transform) đã được phát triển. Một trong số đó là biến đổi của Chan và Ho. Thuật toán này tương tự như thuật toán FFT, ngoại trừ một số biến đổi ở đầu vào dữ liệu. Phép biến đổi này coi rằng N là bội số của 2.
Dữ liệu đầu vào x(n) được sắp xếp lại theo thứ tự:
n = 0,1,2,...,N/2 - 1 (13.5)
Thay vào biểu thức (13.4) chúng ta được
Sau khi nhân ta được:
(13.6)
Cũng giống như phương pháp hệ thập phân trong miền tần số cùng trong FFT, chúng ta cũng chia X(k) thành các giá trị chẵn và lẻ theo các bước sau:
Các mục chỉ số chẵn:
có thể biểu diễn thành:
(13.7)
Các mục chỉ số lẻ:
biểu thức này sau khi nhân một vài phép nhân, có thể biểu diễn thành:
(13.8)
Nhưng cos[(2k + 1)f] = 2cosfcos(2kf) - cos[(2k - 1) f]. Vì vậy,
(13.9)
Biểu thức thứ hai trong (13.9) là X(2k - 1). Bởi vậy, nên (13.9) rút gọn thành
(13.10)
Đặt (13.11)
và (13.12)
Vì vậy:
(13.13)
Và
(13.14)
Biểu thức (3.13) và (3.14) dẫn chúng ta đến sơ đồ hình 13.4 cho DCT 8 điểm. Chú ý rằng X(1) = X(-1). Biểu thức (13.11) và (13.12) là biểu thức cho các phép toán bướm:
Đặt
(13.15)
-1
-1
-1
-1
DCT
4 điểm.
x(0)
x(2)
x(4)
x(6)
DCT
4 điểm.
2x(1)
x(3)+x(1)
x(5)+x(3)
x(7)+x(5)
Hình 13.4 Lưu đồ cho biểu thức (13.13) và (13.14).
Vì thế: (13.16)
(13.17)
Chia Y00(k) và Y01(k) thành hai dãy chỉ số chẵn và lẻ như trên chúng ta được:
(13.18)
s (13.19)
(13.20)
(13.21)
Đặt (13.22)
(13.23)
(13.24)
(13.25)
ở đây (13.26)
-1
-1
-1
-1
2C58X11(1)
2C8X11(0)
DCT
DCT
X00(0)
X00(1)
X00(2)
X00(3)
Y00(0)
Y00(2)
2Y00(1)
Y00(3)+Y00(1)
X10(0)
X10(1)
2C58X13(1)
2C8X13(0)
DCT
DCT
X01(0)
X01(1)
X01(2)
X01(3)
Y01(0)
Y01(2)
2Y01(1)
Y01(3)+Y01(1)
X12(0)
X12(1)
Hình 13.5 Bước thứ hai của thuật toán biến đổi cosin.
Biểu thức (13.22) đến (13.25) biểu diễn cho toán tử bướm. Thay các biểu thức này vào các biểu thức từ (13.18) đến (13.21) chúng ta có:
(13.27)
(13.28)
(13.29)
(13.30)
Các biểu thức trên cho kết quả của bước thứ hai trên sơ đồ hình (13.5) trong trường hợp N = 8.
Đặt (13.31)
(13.32)
Vì vậy (13.33)
(13.34)
(13.35)
(13.36)
(13.37)
(13.38)
Nếu N = 8, các biểu thức trên có dạng:
(13.39)
(13.40)
(13.41)
(13.42)
ở đây k = 0,1.
Các biểu thức cuối có thể biểu diễn thành
Y10(0) = x10 (0) + x10(1)
Y10 = (x10(0) - x10(1))cos(p/4)
...
v..v...
Các biểu thức này dẫn chúng ta đến lưu đồ bướm cuối cùng trình bày ở hình 13.6.
Để rút ra tín hiệu đầu ra của lưu đồ FCT chúng ta cần quay trở lại. Cho ví dụ, từ biểu thức (13.21) và (13.15) chúng ta có thể viết:
(13.43)
x11(0)
x11(1)
Y11(0)
Y11(1)
1
C4
-1
x12(0)
x12(1)
Y12(0)
Y12(1)
1
C4
-1
x14(0)
x14(1)
Y14(0)
Y14(1)
1
C4
-1
x13(0)
x13(1)
Y13(0)
Y13(1)
1
C4
-1
x10(0)
x10(1)
Y10(0)
Y10(1)
1
C4
-1
Hình 13.6 Bướm cuối cùng trong FCT.
Bước cuối cùng của các thao tác bướm.
Đầu ra
Bitdịch chuyển
000
001
010
011
100
101
110
111
0
1
2
3
4
5
6
7
X(0)
X(4)
2X(2)
X(6)+X(2)
2X(1)
X(5)+X(3)
2(X(3)+X(1))
X(7)+X(5)+X(3)+X(1)
Dịch chuyển bit
000
100
010
110
001
101
011
111
0
4
2
6
1
5
3
7
Hình 13.7 Tín hiệu ra cuối cùng của thuật toán bướm tính dưới dạng hệ số FCT.
X(0)
2X(1)
2X(2)
2(X(3)+X(1))
X(4)
X(5)+X(3)
X(6)+X(2)
X(7)+X(5)+X(3)+X(1)
0
1
2
3
4
5
6
7
Dịch chuyển bit
Hình 13.8 Tín hiệu ra sau dịch chuyển bít.
-1
-1
-1
-1
-1
X(0)
X(1)
X(2)
X(3)
X(4)
X(5)
X(6)
X(7)
Dịch chuyển bit
Hình 13.9 Bước cộng truy hồi.
Tương tự, chúng ta có thể rút ra từ các biểu thức (13.32) đến (13.34) và (13.15) các biểu thức dưới đây:
(13.44)
(13.45)
(13.46)
Vì vậy, tín hiệu ra từ các bước cuối cùng của các thao tác bướm có thể tính dưới dạng các hệ số của FCT như trong hình 13.7. Nếu chúng ta sắp xếp lại vị trí của các giá trị bít dùng dịch chuyển bít, chúng ta rút ra tín hiệu ra như hình (13.8). Sau đó,các tín hiệu ra này có thể được dùng để rut ra FCT như hình (13.9). Bước cuối cùng này gọi là bước cộng truy hồi.
Có rất nhiều bước (các thao tác bướm, dịch chuyển bít, và cộng truy hồi) bây giờ có thể nằm trong một bước trong hình 13.10. Từ sơ đồ này, chương trình FCT có thể phát triển. Chương trình này dùng các thuật toán phát triển cho FFT. Chương trình dùng một bảng tra cứu để chứa các giá trị cosin, một bảng cho dịch chuyển bit. Chi tiết của chương trình này để lại cho người dùng như một bài tập.
Bài tập13.5
Phát triển lưu đồ cho DCT 16 và 32 điểm.
Sử dụng logic tương tự đã được dùng cho việc phát triển chương trình FFT, để viết một thuật toán cho FCT.
2-D FCT của một dãy 2-D thực được cho bởi
(13.47)
Thuật toán cho 2-D FCT có thể phát triển dùng phương pháp hàng-cột bình thường như thuật toán 2-D FFT. Chương trình 13.5 rút ra biến đổi FCT của các khối người dùng tự xác định kích thước, thông thường là 8 ´ 8 hoặc là 16 ´ 16, bằng các chia nhỏ các khối của ảnh.
ảnh được giả sử là có chiều dài bằng chiều rộng với cac chiều là bội của 2. Chương trình sử dụng thuật toán 1-D FCT và tận dụng các bảng tra cứu (LUT) như các trạng thái trước.
Chương trình 13.5 “FCT2D.C” 2-D FCT các khối ảnh (người dùng tự xác định kích thước).
/*Program 13.5 "FCT2D.C". 2-D FCT of image blocks (user-pecified).*/
/* This program is for carrying out
2-D Fast Cosine Transform of blocks
subdividing a given image. Block size is
user specified. */
#define PI 3.14159
#include
#include
#include
#include
#include
#include
void FCT(float *, unsigned int *, float *, int , int );
void bit_reversal(unsigned int *, int , int );
void WTS(float *,int , int );
void main()
{
int m,N,i,j,k1,k2,NB,NS,NT,k;
float *x,*C,**w;
unsigned int *L;
double nsq;
FILE *fptri,*fptro;
char file_name[14];
float *buffi;
unsigned char *buffo;
clrscr();
printf("Enter name of input file --- >");
scanf("%s",file_name);
fptri=fopen(file_name,"rb");
if(fptri==NULL)
{
printf("No such file exists.\n");
exit(1);
}
nsq=(double)filelength(fileno(fptri));
/*----------------------
Assume image is square.
----------------------*/
N=(int)sqrt(nsq);
m=(int)(log10((double)N)/log10((double)2.0));
k=1 ;
for(i=0;i<m;i++)
k<<=1;
if(k!=N)
{
printf
("Image must have dimensions which are multiples of 2.\n");
exit(1);
}
printf("Enter name for output file -->");
scanf("%s",file_name);
fptro=fopen(file_name,"wb");
printf ("\nEnter block size (e.g. 8x8, 16x 16, etc.)--> ");
scanf ( "%dx%d", &NB,&NB);
m=(int)(log10((double)NB)/log10((double)2.0));
NT=N*NB;
/*------------------------------
Blocks of NBxNB will be considered.
-------------------------------*/
/* ------------------
Assigning memory.
----------------- */
buffi=(float *)malloc(NT*sizeof(float));
buffo=(unsigned char *)malloc(NT*sizeof(char));
x=(float *)malloc(NB*sizeof(float));
L=(unsigned int *)malloc(NB*sizeof(int));
C=(float *)malloc((NB-1)*sizeof(float));
w=(float **)malloc(NB*sizeof(float *));
for(i=0;i<NB;i++)
*(w+i)=(float *)malloc(NB*sizeof(float));
bit_reversal(L,m,NB);
WTS(C,m,NB);
NS=(N>>m);
/* ------------------------------------
2-D FCT using the row-column approach.
----------------------------------- */
for(i=0;i<NS;i++)
{
fread(buffi,NT,1,fptri);
for(j=0;j<NS;j++)
{
for(k1=0;k1<NB;k1++)
for(k2=0;k2<NB;k2++)
w[k1][k2]=buffi[k1*N+k2+j*NB];
for(k1=0;k1<NB;k1++)
{
for(k2=0;k2<NB;k2++)
x[k2]=w[k1][k2];
FCT(x,L,C,m,NB);
for(k2=0;k2<NB;k2++)
w[k1][k2]=x[k2];
}
for(k2=0;k2<NB;k2++)
{
for(k1=0;k1<NB;k1++)
x[k1]=w[k1][k2];
FCT(x,L,C,m,NB);
for(k1=0;k1<NB;k1++)
w[k1][k2]=x[k1];
}
for(k1=0;k1<NB;k1++)
for(k2=0;k2<NB;k2++)
buffo[k1*N+k2+j*NB]=w[k1][k2];
}
fwrite(buffo,NT,sizeof(float),fptro);
}
fcloseall();
}
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
C4
C4
C4
C4
Dịch chuyển bit.
X(0)
X(1)
X(2)
X(7)
X(3)
X(4)
X(5)
X(6)
Hình 13.10 Biểu đồ chuyển đổi cosin nhanh.
void FCT(float *x, unsigned int *L,
float *C, int m, int N)
{
int NK1,NK2,i,j,k,kk,ip,incr,L1,k1,k2,iter;
float T;
/*Rearranging the order of the input sequence.*/
NK1=N>>1;
T=x[2];
x[2]=x[1];
x[1]=T;
k1=2;
k2=4;
for(i=1;i<(NK1-1);i++)
{
T=x[k2];
for(j=0;j<=(k2-k1);j++)
{
x[k2-j]=x[k2-j-1];
}
x[k1]=T;
k1++;
k2+=2;
}
NK2=N-1;
for(i=0; i>1);i++)
{
T=x[NK2-i];
x[NK2-i]=x[NK1+i];
x[NK1+i]=T;
}
/*Forward operation. */
ip=NK1;
kk=0;
incr=N;
for(iter=0;iter<m;iter++)
{
for(k=0;k<ip;k++)
{
for(j=k;j<N;j+=incr)
{
i=j+ip;
T=x[j];
x[j]=T+x[i];
x[i]=T-x[i];
x[i]*=C[kk];
}
kk++;
}
ip>>=1;
incr>>=1;
}
/*Bit reversal. */
for(i=0; i<(N-1); i++)
{
if(i<=L[i]) continue;
else
{
T=x[i];
x[i]=x[L[i]];
x[L[i]]=T;
}
}
/*Recursive subtraction. */
for(i=1;i<NK1;i++)
x[i]*=0.5;
NK1=(N>>2);
NK2=(N>>1);
kk=1;
for(iter=1;iter<m;iter++)
{
kk<<=1;
L1=kk-1; /*L1=2^iter-1*/
for(k=1;k<=L1;k++)
for(i=0;i<NK1;i++)
x[NK1+k*NK2+i]=x[NK1+k*NK2+i]-x[NK1+(k-1)*NK2+i];
NK1>>=1;
NK2>>=1;
}
x[0]*=1.414/(float)N;
for(i=1;i<N;i++)
x[i]*=2.0/(float)N;
}
void bit_reversal(unsigned int *L, int m, int N)
{
unsigned int MASK,C,A,j,k,i;
for(k=0;k<N;k++)
{
MASK=1;
C=0;
for(i=0,j=m-1;i<m;i++,j--)
{
A=(k&MASK)>>i;
A<<=j;
C|=A;
MASK<<=1;
}
L[k]=C;
}
}
void WTS(float *C,int m, int N)
{
int NK1,NK2,i,k,iter;
NK1=N>>1;
NK2=N<<1;
k=0;
for(iter=0;iter<m;iter++)
{
for(i=0;i<NK1;i++)
{
C[k]=
(float)cos((double)((PI*(float)(4*i+1))/(float)NK2));
k++;
}
NK1>>=1;
NK2>>=1;
}
for(i=0;i<(k-1);i++)
C[i]*=2.0;
}
Biến đổi ngược DCT được cho bỏi
(13.48)
ở đây
Thay vì lập lại cả một chương trình để tính biến đổi ngược FCT, chúng ta sẽ dùng lưu đồ của FCT tiến (forward). Để rút ra FCT ngược, tất cả các việc mà chúng ta cần làm là đảo ngược biểu đồ ở hình 13.10. Hình 13.11 giới thiệu phép biến đổi ngược của một bướm. Kết quả của thuật toán biến đổi ngược của lưu đồ FCT cho trong hình 13.12. Từ lưu đồ của hình 13.11 biến đổi ngược của FCT có thể phát triển bình thường từ FCT tiến. Chương trình tính 2-D FCT ngược cho trong chương trình 13.6. Giá trị của khối ảnh gốc phải được cho trước bởi người dùng.
0.5
A
B
C
D
CX
-1
A=C+D
B=(C-D)CX
C
D
A
B
-1
C=0.5A+B/(2CX)
D=0.5A-B/(2CX)
1/2CX
Hình 13.11 Phép đổi ngược của một bướm.
0.5
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
+
+
+
+
+
Dịch chuyển bit.
X(0)
X(1)
X(2)
X(3)
X(4)
X(6)
X(7)
X(5)
Hình 13.12 Biểu đồ đảo ngược giải thuật FCT.
Chương trình 13.6 "IFCT2D.C". Đảo ngược FCT. Kích thước khối sử dụng trên ảnh gốc.
/*Program 13.6 "IFCT2D.C". Inverse FCT. Block size used on the original image should be known to the user.*/
/* This program is for carrying out
the inverse of the 2-D Fast Cosine Transform. */
#define PI 3.14159
#include
#include
#include
#include
#include
#include
void IFCT(float *, unsigned int *, float *, int , int );
void bit_reversal(unsigned int *, int , int );
void WTSINV(float *,int, int);
void main()
{
int m,N,i,j,k1,k2,NB,NT,NS,k;
float *x,*C,**w;
unsigned int *L;
double nsq;
FILE *fptri, *fptro;
char file_name[14] ;
float *buffi ;
unsigned char *buffo;
clrscr();
printf ( "This program is for the inverse 2-D FCT . \n" ) ;
printf ( "Enter name of input file --->") ;
scanf ("%s " , file_name) ;
fptri=fopen(file_name, " rb");
if( fptri ==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
nsq=(double)filelength(fileno(fptri));
/*--------------------
Assume image is square,
----------------------*/
N=(int)sqrt(nsq/sizeof(float));
m=(int)(log10((double)N)/log10((double)2.0));
k=1;
for(i=0; i<m;i++)
k<<=1;
if(k!=N)
{
printf("\nTransform image must have dimensions");
printf(" which are multiples of 2. \n ");
exit(1);
}
printf("Enter name for Output file - - -> " ) ;
scanf("%s", file_name);
fptro=fopen(file_name,"wb");
printf("\nEnter block size used (e.g. 8x8,16xI6,etc.) -->");
scanf("%dx%d",&NB,&NB);
m=(int)(log10((double)NB)/log10((double)2.0));
NT=N*NB;
/*-----------------------------
Blocks of NBxNB were considered.
--------------------------------*/
/* ------------------
Assigning memory.
----------------- */
buffi=(float *)malloc(NT*sizeof(float));
buffo=(unsigned char *)malloc(NT*sizeof(char));
x=(float *)malloc(NB*sizeof(float));
L=(unsigned int *)malloc(NB*sizeof(int));
C=(float *)malloc((NB-1)*sizeof(float));
w=(float **)malloc(NB*sizeof(float *));
for(i=0;i<NB;i++)
*(w+i)=(float *)malloc(NB*sizeof(float));
bit_reversal(L,m,NB);
WTSINV(C,m,NB);
NS=(N>>m) ;
/*------------------------------------
2-D inverse FCT using the row-column approach.
----------------------------------- */
for(i=0;i<NS;i++)
{
fread(buffi,NT,sizeof(float),fptri);
for(j=0;j<NS;j++)
{
for(k1=0;k1<NB;k1++)
for(k2=0;k2<NB;k2++)
w[k1][k2]=buffi[k1*N+k2+j*NB];
for(k1=0;k1<NB;k1++)
{
for(k2=0;k2<NB;k2++)
x[k2]=w[k1][k2];
IFCT(x,L,C,m,NB);
for(k2=0;k2<NB;k2++)
w[k1][k2]=x[k2];
}
for(k2=0;k2<NB;k2++)
{
for(k1=0;k1<NB;k1++)
x[k1]=w[k1][k2];
IFCT(x,L,C,m,NB);
for(k1=0;k1<NB;k1++)
w[k1][k2]=x[k1];
}
for(k1=0;k1<NB;k1++)
for(k2=0;k2<NB;k2++)
buffo[k1*N+k2+j*NB]=(char)fabs((double)w[k1][k2]);
}
fwrite(buffo,NT,1,fptro);
}
fcloseall();
}
void bit_reversal(unsigned int *L, int m, int N)
{
unsigned int MASK,C,A,j,k,i;
for(k=0;k<N;k++)
{
MASK=1;
C=0;
for(i=0,j=m-1;i<m;i++,j--)
{
A=(k&MASK)>>i;
A<<=j;
C|=A;
MASK<<=1;
}
L[k]=C;
}
}
void WTSINV(float *C,int m,int N)
{
int NK,k,i,iter,LK;
NK=4;
k=0;
LK=1;
for(iter=1;iter<=m;iter++)
{
for(i=0;i<LK;i++)
{
C[k]=1.0/(float)cos(PI*(float)(4*i+1)/(float)NK);
k++;
}
NK<<=1;
LK<<=1;
}
C[0]/=2.0;
for(i=1;i<(N-1);i++)
C[i]/=4.0;
}
void IFCT(float *x, unsigned int *L, float *C, int m, int N)
{
int NK1,NK2,i,j,k,kk,ip,incr,L1,k1,k2,iter;
float T;
x[0]*=(float)N/1.414;
for(i=1;i<N;i++)
x[i]*=(float)N/2.0;
/*Recursive addition. */
NK1=1 ;
NK2=2 ;
kk=1;
for(i=0;i<m;i++)
kk*=2 ;
for(iter=1;iter<m;iter++)
{
kk>>=1;
L1=kk- 1 ;
for(k=L1;k>=1;k--)
for(i=0;i<NK1;i++)
x[NK1+k*NK2+i]=x[NK1+k*NK2+i]+x[NK1+(k-1)*NK2+i];
NK1<<=1;
NK2<<=1;
}
NK1=N>>1;
for(i=1;i<NK1;i++)
x[i]*=2.0;
/*Bit reversal. */
for(i=0;i<(N-1);i++)
{
if(i<=L[i]) continue;
else
{
T=x[i];
x[i]=x[L[i]];
x[L[i]]=T;
}
}
/*Forward operation. */
ip=1;
kk=0;
incr=2;
for(iter=0;iter<m;iter++)
{
for(k=0;k<ip;k++)
{
for(j=k;j<N;j+=incr)
{
i=j+ip;
T=x[j]*0.5;
x[i]=x[i]*C[kk];
x[j]=T+x[i];
x[i]=T-x[i];
}
kk++;
}
ip<<=1;
incr<<=1;
}
/*Rearranging the order of the input sequence.*/
kk=1;
for(i=1;i>1);i++)
{
T=x[N-1];
k=1;
for(j=kk;j<N;j++)
{
x[N-k]=x[N-k-1];
k++;
}
x[kk]=T;
kk+=2;
}
}
Bài tập 13.6
1. Chạy chương trình 2-D FCT (Chương trình 13.5) trên ảnh “IKRAM.IMG” và dùng một khối có kích thước 8 ´ 8. Tên của file xuất ra đặt là IKRAMFCT.IMG.
2. Lấy 2-D FCT của IKRAMFCT.IMG bằng chương trình 13.6. Chứa kết quả trên ảnh IKRAMR.IMG.
3. Hiển thị IKRAMR.IMG. So sánh ảnh này với ảnh IKRAM.IMG.
Tiếp theo tôi sẽ giới thiệu một chương trình mà sẽ hiển thị các khối 8 ´ 8 của kết quả đã chuyển đổi trên màn hình văn bản. Chú ý rằng giá trị sẽ giảm xuống một cách nhanh chóng kể từ góc trái của màn hình, hay là điểm tần số (0,0). Để thoát khỏi chương trình này, bấm ESC.
Chương trình 13.7 "DISPFCT.C". Chương trình hiển thị khối 2-D FCT.
/*Program 13.7 "DISPFCT.C". Program to display 2-D FCT blocks.*/
/* This program displays the result of the 2-D FCT. */
#include
#include
#include
#include
#include
#include
void main( )
{
int i,j,k1,k2,N,N1,NB,NS;
float *buff,nsq;
FILE *fptr;
char file_name[14];
int ch;
clrscr();
printf("Enter input file name containing 2-D FCT result.-->");
scanf("%s",file_name);
fptr=fopen(file_name,"rb");
if(fptr==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
nsq=(float)filelength(fileno(fptr));
N=sqrt((double)(nsq/sizeof(float)));
printf("\nEnter block size used (8x8, 16x16, etc.)-->");
scanf("%dx%d",&NB,&NB);
N1=N/NB;
NS=NB*N;
buff=(float *)malloc(NS*sizeof(float));
printf("\nPress ESC to exit.");
clrscr();
for(i=0;i<N1;i++)
{
fread(buff,sizeof(float),NS,fptr);
{
for(j=0;j<N1;j++)
{
for(k1=0;k1<NB;k1++)
{
for(k2=0;k2<NB;k2++)
printf(" %5.2f",buff[k1*N+k2+j*NB]);
printf("\n");
}
ch=getch();
if(ch==27) {
clrscr();
fclose(fptr);
exit(1);
}
clrscr();
gotoxy(1,1);
}
}
fclose(fptr) ;
}
}
Một mẫu của 2-D FCT từ một khối rút ra từ ảnh IKRAM.IMG giới thiệu trong hình 13.13. Các thành phần một chiều (các khối có tần số 0,0) có giá trị lớn nhất. Thành phần này biểu diễn giá trị trung bình của khối. Các thành phần khác (được hiểu là các thành phần xoay chiều) phải có các giá trị nhỏ hơn.
13.4.2 Sự phóng to và sự thu nhỏ của các ảnh
Một trong các ứng dụng của FCT là thay đổi kích thước của ảnh, bao gồm sự phóng to và sự thu nhỏ ảnh. Chúng ta sẽ chấp nhận rằng một ảnh phóng to dùng FCT thì sẽ không có gì khác so với một ảnh phóng to dùng FFT, trong đó ta dùng các thuộc tính vốn có của FFT của ảnh. Toàn bộ ý tưởng là ta sẽ chia tần số lấy mẫu bên trong với 2 và như vậy sẽ tăng gấp đôi khoảng cách tần số. Khi FCT của một khối giảm rất nhanh kể từ điểm tần số (0,0), chúng ta có thể nhân đôi khoảng cách tần số bằng cách thêm vào các điểm 0. Điều này làm cho chúng ta giảm được một nửa tần số lấy mẫu bên trong và tăng gấp đôi kích thước của ảnh. Vì vậy trong tương lai, kỹ thuật truyền hình phải sử dụng kỹ thuật ảnh số, thì tính chất này có thể giúp cho ta tạo ra được các tivi có màn ảnh rộng. Sự thu nhỏ thực hiện bằng cách bớt đi các giá trị của tần số cao.
Chương trình sau đây sẽ phóng to của một ảnh. Chương trình này yêu cầu FCT của một ảnh và kích thước của các khối trên ảnh gốc. Kết quả chạy chương trình này là một ảnh có kích thước gấp đôi.
Chương trình 13.8 "IFCT2DX.C". Chương trình phóng to ảnh gấp đôi dùng 2-D FCT.
/*Program 13.8 "IFCT2DX.C". Program for image
doubling using the 2-D FCT.
The 2-D FCT of the image is the input file.*/
/* This program is for carrying out the inverse
of the 2-D Fast Cosine Transform.
It doubles the size of the image at the same time. */
#define PI 3.14159
#include
#include
#include
#include
#include
#include
void IFCT(float *, unsigned int *, float *, int , int );
void bit_reversal(unsigned int *, int , int );
void WTSINV(float *,int, int);
void main()
{
int m,N,i,j,k1,k2,NB,NS,N1,NB1,m1,NT,NT1,k;
float *x,*C,**w;
unsigned int *L;
double nsq;
FILE *fptri,*fptro;
char file_name[14];
float *buffi;
unsigned char *buffo;
clrscr();
printf ("This program is for the inverse 2-D FCT. \n" ) ;
printf("It doubles the size of the image at the same time.");
printf("\nEnter name of input file --- > " ) ;
scanf("%s",file_name);
fptri=fopen(file_name,"rb");
if(fptri==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
nsq=(double)filelength(fileno(fptri));
/*---------------------
Assume image is square.
-----------------------*/
N=(int)sqrt(nsq/sizeof(float));
m=(int)(log10((double)N)/log10((double)2.0));
k=1;
for(i=0;i<m;i++)
k<<=1;
if(k!=N)
{
printf ( "\nTransform Image must have dimensions" ) ;
printf ( " which are multiples of 2. \n" );
exit(1);
}
printf("Enter name for output file --->");
scanf( "%s" , file_name);
fptro=fopen(file_name,"wb");
printf("\nEnter block size used (e.g. 8x8,16xl6) --- >");
scanf ( "%dx%d" , &NB , &NB ) ;
m=(int)(log10((double)NB)/log10((double)2.0));
NT=N*NB;
/*-----------------------------
Blocks of NBxNB were considered by the
FCT. Output blocks are (2xNB)x(2xNB).
-------------------------------*/
m1=m+1; NB1=NB<<1;
N1=(N<<1);
NT1=N1*NB1;
/*------------------
Assigning memory.
-----------------*/
buffi=(float *)malloc(NT*sizeof(float));
buffo=(unsigned char *)malloc(NT1*sizeof(char));
x=(float *)malloc(NB1*sizeof(float));
L=(unsigned int *)malloc(NB1*sizeof(int));
C=(float *)malloc((NB1-1)*sizeof(float));
w=(float **)malloc(NB1*sizeof(float *));
for(i=0;i<NB1;i++)
*(w+i)=(float *)malloc(NB1*sizeof(float));
bit_reversal(L,m1,NB1);
WTSINV(C,m1,NB1);
NS=(N>>m);
/* ------------------------------------
2-D IFCT using the row-column approach.
-----------------------------------*/
for(i=0;i<NS;i++)
{
fread(buffi,NT,sizeof(float),fptri);
for(j=0;j<NS;j++)
{
for(k1=0;k1<NB1;k1++)
for(k2=0;k2<NB1;k2++)
w[k1][k2]=0.0;
for(k1=0;k1<NB;k1++)
for(k2=0;k2<NB;k2++)
w[k1][k2]=buffi[k1*N+k2+j*NB];
for(k1=0;k1<NB1;k1++)
{
for(k2=0;k2<NB1;k1++)
x[k2]=w[k1][k2];
IFCT(x,L,C,m1,NB1);
for(k2=0;k2<NB1;k2++)
w[k1][k2]=x[k2];
}
for(k2=0;k2<NB1;k2++)
{
for(k1=0;k1<NB1;k1++)
x[k1]=w[k1][k2];
IFCT(x,L,C,m1,NB1);
for(k1=0;k1<NB1;k1++)
w[k1][k2]=x[k1];
}
for(k1=0;k1<NB1;k1++)
for(k2=0;k2<NB1;k2++)
buffo[k1*N1+k2+j*NB1]=(unsigned char)fabs(w[k1][k2]);
}
fwrite(buffo,NT1,1,fptro);
}
fcloseall();
}
Bài tập 13.7
Sử dụng chương trình 13.8 cho IKRAMFCT.IMG (ảnh mà bạn đã thu được trước đây). Hiển thị ảnh đã được phóng to.
Lặp lại công việc, sử dụng khối 16 ´ 16. Chú ý rằng bạn sẽ phải chạy lại chương trình 2-D FCT.
Thay đổi chương trình trên cho ảnh thu nhỏ.
13.5 Lượng tử hoá
Trong toàn bộ cuốn sách này chúng ta làm việc với các ảnh đã lượng tử hoá. Chúng ta coi rằng các ảnh được chia làm các mức xám dùng 8 bit với các giá trị từ 0 đến 255. Tuy vậy, chúng ta không bao giờ để ý đến các giá trị này. Lý do là các giá trị này bản thân nó không có ý nghĩa gì cả. Cái mà chúng ta thật sự quan tâm là biểu diễn dưới dạng nhị phân của giá trị cường độ sáng tín hiệu lấy mẫu. Tất cả các bít biểu diễn khoảng của tín hiệu chói. Nếu tín hiệu chói có giá trị trong khoảng từ giá trị nhỏ nhất yL cho đến giá trị lớn nhất yu, thì hàm lượng tử hoá sẽ chia tín hiệu ra thành N miền và gán cho mỗi miền một giá trị nhị phân như trong hình 13.14. Sự biểu diễn này không có ý nghĩa gì về mặt vật lý, và chức năng của sự biểu diễn này, như chúng ta muốn, chỉ dùng trong lĩnh vực xử lý tín hiệu số. Tất cả các tín hiệu số này gọi là điều mã xung (Pulse Code Modulated - PCM). Để có thể thực sự thấy giá trị các mức xám chúng ta cần lượng tử hoá ngược. Trong bước này, các giá trị nhị phân biểu diễn một độ chói cụ thể. Các bước thực sự của quá trình này biểu diễn trong hình 13.15. Trong lĩnh vực tương tự và lĩnh vực số quá trình này gọi là chuyển đổi từ tương tự sang số (A/D) và chuyển đổi từ số sang tương tự (D/A).
Trong các ứng dụng như trường hợp biến đổi cosin 2-D thì có một chút khác biệt. Cái mà chúng ta cần làm trong trường hợp này là biến đổi từ một tập hợp các dấu phẩy động sang một tập hợp các bít nhị phân và ngược lại. Biểu vì biến đổi ngược của lượng tử hoá là biến đổi từ nhiều vào một, nên quá trình này không thể tiến hành một cách thông thường được. May mắn thay, có một số phương pháp để lượng tử hoá và lượng tử hoá ngược. Chúng ta sẽ nghiên cứu các phương pháp này ở phần dưới đây.
dN
dN-1
d5
d4
d3
d2
d1
d0
N-1
N-2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4
3
2
1
0
11111111
11111110
.
.
.
.
.
.
.
.
.
.
.
.
.
.
00000100
00000011
00000010
00000001
00000000
YU
YL
Khoảng tín hiệu.
Dạng thập phân.
Dạng nhị phân.
di { i= 0,...,1} là các mức chia.
Hình 13.14 Lượng tử hoá.
13.5.1 Lượng tử hoá đồng đều
Đây là dạng đơn giản nhất của lượng tử hoá. Trong dạng lượng tử hoá này, khoảng (yu - yL) được chia thành N khoảng cách đều nhau (xem trong hình 13.14). Các giá trị từ d0 đến dN được gọi là các mức chia. Các mức lượng tử biểu diễn giá trị thực của các mức chia trong khoảng từ di đến di+1 dưới dạng số nhị phân bằng i.
Vì vậy, nếu như di+1 < y Ê di thì giá trị của lượng tử đầu ra = i.
Mức lượng tử đầu ra có thể biểu diễn theo công thức:
ở đây dấu có nghĩa là làm tròn thành số nguyên gần nhất.
Lượng tử hoá ngược dùng một bảng biến đổi ngược giữa giá trị lượng tử i và biến đổi ngược của nó ri
Tuy nhiên lượng tử hoá đồng đều không quan tâm đến khả năng xảy ra của sự kiện với các giá trị được đưa ra. Tổng quát, quá trình lượng tử hoá này áp dụng cho trường hợp tất cả các mức có khả năng xuất hiện bằng nhau. Điều này, trong hầu hết các trường hợp là không đúng. Dễ nhận thấy là các mức lượng tử hoá tập trung nhiều nhất vào miền mà khả năng xuất hiện của các mức xám nhiều nhất. Điều này dẫn chúng ta đến phương pháp thiết kế lượng tử hoá dưới đây.
Tín hiệu ánh sáng.
Tín hiệu tương tự.
t
t
Lượng tử hoá ngược hay biến đổi số sang tương tự.
Bộ lọc thông thấp.
Tín hiệu dạng nhị phân.
Mức được tái thiết.
Hình 13.15 Lượng tử hoá ngược.
13.5.2 Lượng tử hoá không đồng đều
Trong phần này chúng ta sẽ xem xét phương pháp lượng tử hoá không đồng đều tối ưu nhất trong hệ thống PCM. Các nghiên cứu cho phương pháp này đã được Panter và Dite đưa ra trong một cuốn sách xuất bản vào năm 1949. Trong cuốn sách này họ đã đưa giải thuật cho lượng tử hoá không đồng đều. Họ đưa ra một phương pháp xấp xỉ tối ưu cho lượng tử hoá không đồng đều. Giải thuật này sẽ không đúng cho các trường hợp quá trình lượng tử hoá có quá ít mức chia. Tuy nhiên các giải thuật này được phát triển một các trọn vẹn trong một báo cáo chưa được xuất bản của Lloyd vào năm 1957 và được Max kiểm nghiệm vào năm 1960. Một phương pháp lượng tử hoá kết hợp cả hai phương pháp của Lloyd và Max thường được gọi là phương pháp lượng tử hoá Lloyd-Max. Trong phần báo cáo xuất bản sau đó của Lloyd xuất bản vào năm 1982 đã cho thấy có rất nhiều ứng dụng rất thú vị của phương pháp này. Bản báo cáo này có hai phương pháp thiết kế, một phương pháp trong đó giống phương pháp của Max. Phương pháp này gọi là phương pháp II. Phương pháp I tỏ ra có nhiều ứng dụng và dễ tính toán hơn phương pháp II. Cả hai phương pháp thiết kế này đều được trình bày ở phần dưới đây.
Nếu chúng ta coi rằng các mức lượng tử hoá được cho bởi
di , i = 0 ... N
(xem hình 13.16) và các mức khôi phục cho bởi
ri , i = 0 ... N
và giá trị đo của tất cả các mức này cho bởi:
(13.49)
ở đây y là tín hiệu đầu vào còn p(y) là khả năng xuất hiện của y.
Hình 13.16 Các mức lấy mẫu và khôi phục.
dN
dN-1
d5
d4
d3
d2
d1
d0
Các mức chia.
yL
yu
y
mẫu i
Lượng tử hoá.
Lượng tử hoá ngược.
y(t)
ri
N-1
3
2
1
0
Mẫu
i ri
r0
. .
. .
. .
N-1 rN-1
LUT cho lượng tử hoá ngược.
Hình 13.17 Các xử lý lấy mẫu và khôi phục.
Hình 13.17 cung cấp sơ đồ khối của quá trình lượng tử hoá và lượng tử hoá ngược. Tín hiệu vào y(t) phải được coi là đã biết khả năng xuất hiện. Vấn đề đặt ra là phải xác định các mức lấy mẫu và các mức khôi phục sao cho méo tín hiệu là nhỏ nhất.
Viết lại biểu thức (13.49 ):
(13.50)
Vi phân (13.50) theo dI và cho biểu thức này bằng không chúng ta được
(13.51)
i = 1,2,3, ... ,N
Lấy vi phân (13.50) theo ri chúng ta được.
(13.52)
i = 0, 1, ... , N - 1.
Hình 13.18 Phương pháp Newton-Raphson cho tính các biểu thức trong ngoặc.
Biểu thức (13.50) và biểu thức (13.51) đưa ra phương pháp xác định các mức lấy mẫu và các mức khôi phục dùng cho cả phương pháp của Lloyd-Max hoặc là phương pháp Lloyd.
Phương pháp Lloyd-Max Lloyd và Max đã phát triển độc lập thuật toán để giải quyết biểu thức (13.50) và (13.51). Các chi tiết của thuật toán này vẫn chưa được cung cấp. Trong phần này tôi sẽ cung cấp cho bạn một thuật toán dựa trên thuật toán Lloyd-Max nhưng có các chi tiết cụ thể hơn. Tôi cũng sẽ cung cấp cho bạn phần mềm thiết kế N mức lấy mẫu và khôi phục.
Thuật toán này gồm các bước sau:
1. Chọn một giá trị cho r0. d0 và dN được coi là đã biết.
2. Cho i = 1,2,...,N - 1.
a. Tính di từ
b. Tính ri từ
3. Tính
4. Nếu rN-1 ạ r/, thay đổi lại r0 và lặp lại các bước từ bước 2 cho đến bước 4.
Bây giờ tôi sẽ cung cấp cho bạn các chi tiết cần thiết để tạo ra thuật toán trên.
Chi tiết cho việc tính di trong bước 2a của thuật toán Lloyd-Max. di có thể tính theo hàm sau đây:
(13.53)
Có thể rút ra biểu thức gốc theo công thức lặp Newton-Raphon được cho bởi:
(13.54)
ở đây l là số lần lặp và f/(di) là đạo hàm của f(di) theo di cho theo công thức:
(13.50)
d0i là giá trị ban đầu. Phép lặp diễn ra cho đến khi
Giá trị ban đầu cho d1 là d0 + D, cho d2 là d1 + D, ..., ở đây D là một giá trị nhỏ. Giá trị gốc của f(di) có thể tính theo dùng phương pháp nửa lặp (bisection). Ưu điểm của phương pháp Newton-Raphson là khả năng hội tụ nhanh. Nhược điểm là đạo hàm của một hàm thường có giá trị rất nhỏ và dễ dẫn đến giá trị zero gây nên sự không ổn định của các số.
Các chi tiết cho việc thay đổi r0. Giá trị r0 có thể thay đổi lại nếu chúng ta nhận thấy rằng giá trị gốc của hàm:
(13.56)
Giá trị gốc này có thể rút ra dùng các giả thiết của Newton-Raphson theo:
(13.57)
Đạo hàm của g(r0) có thể rất khó khăn cho việc phân tích. Trong trường hợp này cần có một công cụ tính toán khác. Chúng ta có thể thay thế đạo hàm bằng một giá trị hằng số có cùng dấu như biểu thức trong dấu ngoặc (xem hình 13.8). Trong trường hợp này thì tích phân lâu hội tụ hơn.
Tiếp theo là một chương trình cho tính các mức lượng tử hoá theo phương pháp Lloyd-Max. Tích phân được đưa ra dùng phương pháp tích phân Romberg bởi vì nó chính xác hơn phương pháp tích phân Simpson. Chương trình cho phép bạn thiết kế lượng tử hoá đồng đều, Gauss hoặc là Laplace.
Gauss
Laplace
Đồng dạng:
Chương trình này có thể thay kiểm tra trên cả ba loại phân bố (với s = 1) với 6 bit cụ thể là 64 mức lượng tử hoá. Bạn cần chú ý khi bạn chạy chương trình này bạn cần cho số của lựa chọn. Lựa chọn đầu tiên là số bit bạn muốn xây dựng cho mức lượng tử, lựa chọn thứ hai là lựa chọn về các mức phân bố (phân bố Gauss, phân bố đồng đều, hoặc phân bố Laplace), và thứ ba cũng là lựa chọn cuối cùng là lựa chon trong khi tính biểu thức trong ngoặc của g(r0). Lựa chọn của g(r0) có thể là lựa chọn cho một số, một giá trị hoặc là một giá trị tăng dần. Lý do của sự lựa chọn này là kết quả của phương pháp tính toán có thể là một số không thể tính được, ví dụ khi ta chia cho một số quá nhỏ. Qua kinh nghiệm tôi nhận thấy rằng phân bố Gauss và phân bố đồng đều thường cho kết quả tính toán rất tốt. Còn cho phân bố Laplace, khi dùng tăng giá trị của biểu thức trong ngoặc sẽ không gặp một khó khăn nào việc tạo mã 6 bit.
Chương trình 13.9 “MAXQ1.C” Chương trình cho tính các mức lượng tử Lloyd-Max.
/*Program 13.9 "MAXQI.C". Program for the Lloyd-Max quantizer.*/
/* Program for designing the Lloyd-Max quantizer
for Gauss, uniform or Laplace distribution.*/
#include
#include
#include
#include
#include
double decision_level(double , double);
double f1(double);
double f2(double);
double p(double);
double Romberg(double, double, double (*)(double));
char ch;
void main()
{
double *r,*d,r1,delta,alpha;
int i,m,N,k,xt,yt;
double der,rt,deltal, delta2;
char ch1,file_name[16];
FILE *fptr;
clrscr();
printf("Enter number of bits --- >");
scanf("%d",&m);
N=1<<m;
r=(double *)malloc((N+1)*sizeof(double));
d=(double *)malloc((N+1)*sizeof(double));
printf("Enter choice of distribution:");
printf("\n 1. Gauss.");
printf("\n 2. Uniform.");
printf("\n 3. Laplace. -->(1,2 or 3): ");
while(((ch=getch())!='1')&&(ch!='2')&&(ch!='3'));
putch(ch);
switch(ch)
{
case '1':
r[0]=-5.0; d[0]=-20.0;
break;
case '2':
r[0]=-0.98; d[0]=-1.0;
break;
case '3':
r[0]=-6.0; d[0]=-100.0;
if(m==6) r[0]=-8.0;
}
printf("\n\n The program uses Newton-Raphson's method");
printf("\n to compute the proper value of r[01. You have");
printf("\n a choice between using calculated values, a");
printf("\n fixed value or decreasing values for the ");
printf("\n derivatives at every update or calculated ");
printf("\n values. Always pick the first choice unless ");
printf("\n you encounter numerical problems.");
xt=wherex();
yt=wherey();
gotoxy(1,21);
printf("Recommendations: Guassian or Uniform select I or 2.\n");
printf(" Laplace select 3.");
gotoxy(xt,yt);
printf("\n Enter choice:");
printf("\n 1. Calculated derivative of the error function.");
printf("\n 2. Fixed value.");
printf("\n 3. Decreasing derivative.-->");
while(((ch1=getch())!='1')&&(ch1!='2')&&(ch1!='3'));
putch(ch1);
gotoxy(1,21);
delline();
delline();
delta=20.0;
k=0;
xt=wherex(); yt=wherey();
gotoxy(70,25);
textattr(WHITE+(GREEN<<4)+BLINK);
cputs("WAIT");
while(fabs(delta)>=1.e-7)
{
k++;
if(k>=1000) break;
/* Computing the derivative of the function:
f(r[0])=r[N-1]-r1, numerically. */
switch(ch1)
{
case '1':
rt=r[0]; r[0]=r[0]+0.001;
for(i=1;i<N;i++)
{
d[i]=decision_level(d[i-1],r[i-1]);
r[i]=2.0*d[i]-r[i-1];
}
r1=Romberg(d[N-1],-d[0],f1)/Romberg(d[N-1],-d[0],f2);
deltal=r[N-1]-r1;
r[0]=rt; r[0]=r[0]-0.001;
for(i=1;i<N;i++)
{
d[i]=decision_level(d[i-1],r[i-1]);
r[i]=2.0*d[i]-r[i-1];
}
r1=Romberg(d[N-1],-d[0],f1)/Romberg(d[N-1],-d[0],f2);
delta2=r[N-1]-r1;
r[0]=rt;
der=(deltal-delta2)/0.002;
break;
case '2':
der=5.0;
if(m==6) der=10.0;
break;
case '3':
if(k==1) der=5.0;
else if((k%5)==0) der/=2.0;
break;
}
/*Computing the values of dEil and r[i]. */
for(i=1;i<N;i++)
{
d[i]=decision_level(d[i-1],r[i-1]);
r[i]=2.0*d[i]-r[i-1];
}
r1=Romberg(d[N-1],-d[0],f1)/Romberg(d[N-1],-d[0],f2);
delta=(r[N-1]-r1);
gotoxy(1,18);
printf("\n %d derivative=%f error=%f ", k,der,delta);
if(fabs(delta)>=500.0)
{
printf("\n A numerical problem was encountered.");
printf("\n Restart problem with a different choice.");
exit(1);
}
r[0]-delta/der;
}
gotoxy(70,25);
textattr(WHITE+(BLACK<<4));
cputs(" ");
gotoxy(xt,yt+2);
switch(ch)
{
case '1':
case '3':
printf("\n\n -%c %9.6f ",(char)236,r[0]);
break;
case '2':
printf("\n\n %9.6f %9.6f ",d[0],r[0]);
}
for(i=1;i<N;i++)
printf("\n %9.6f %9.6f ",d[i],r[i]);
switch(ch)
{
case '1':
case '3':
printf("\n +%c",(char)236);
break;
case '2':
printf("\n %9.6f ",-d[0]);
}
printf("\nEnter file name for storing quantizer-->");
scanf("%s",file_name);
fptr=fopen(file_name,"w");
for(i=0;i<N;i++)
fprintf(fptr,"%f %f \n",d[i],r[i]);
fclose(fptr);
}
double p(double x)
{
double prob;
switch(ch)
{
case '1':
prob=exp(-x*x*0.5);
break;
case '2':
prob=1.0;
break;
case '3':
prob=exp(-1.4142136*fabs(x));
break;
}
return prob;
}
#define EPS 1.0e-6
/*----------------------------------
Calculating the next decision level,
given values for d[i-11 and r[i-1].
--------------------------------- */
double decision_level(double di_1, double ri_1)
{
double delta,di,fun,dfun,ff,area1,area2;
int i;
if(fabs(di_1)>10.0) delta=2.0;
else delta=0.1;
di=ri_1+delta;
i=0;
/* Using Newton-Raphson's method for root finding. */
while(fabs(delta)>=EPS)
{
i++;
if(i>=1000)
{
printf("\n no convergence.\n");
return di;
}
area1=Romberg(di-1,di,f1);
area2=Romberg(di-1,di,f2);
ff=area1/area2;
fun=(ri_1)-ff;
if(fabs(area2)<=1.0e-10)
{
printf("\n A numerical problem was encountered.");
printf("\n Restart problem with a different choice.");
exit(1);
}
dfun=-p(di)*(di-ff)/area2;
if(fabs(dfun)<=1.0e-10)
{
gotoxy(1,20);
printf("\n derivative=%f",dfun);
printf("\n A numerical problem was encountered.");
printf("\n Restart problem with a different choice.");
exit(1);
}
delta=-fun/dfun;
di+=delta;
}
return di;
}
/*-------------------------
functions used by numerical integration routine.
-------------------------------- */
double f1(double y)
{
return (y*(p(y)));
}
double f2(double y)
{
return p(y);
}
/* -----------------------------
Numerical integration using
Romberg method.
---------------------------- */
double Romberg(double di_1,double di,double (*f)(double))
{
double T[10][10],h,x,Area;
int N,k,kk,i,j;
N=9;
k=1 ;
h=di- di_1;
T[0][0]=h*((*f)(di_1)+(*f)(di))/2.0;
for(i=1;i<=N;i++)
{
k<<=1;
h/=2.0;
x=di_1 - h;
Area=0.0;
T[0][i]=(T[0][i-1])/2.0;
for(j=0;j<(k-1);j+=2)
{
x+=2.0*h ;
Area+=(*f)(x);
}
T[0][i]+=Area*h;
kk=1;
for(j=1;j<=i;j++)
{
kk<<=2;
T[j][i]=(kk*T[j-1][i]-T[j-1][i-1])/(float)(kk-1);
}
}
return T[N][N];
}
Bài tập 13.8
1. Chạy chương trình 13.9 dùng các lựa chọn sau:
3 bit, 4 bit, 5 bit, 6 bit.
2. Với mỗi lựa chọn số bit ta chọn:
Gauss.
Đồng đều.
Laplace.
3. Các mức lượng tử có thể dùng cho thiết kế có trên đĩa, ví dụ dùng Q3G.DAT cho lượng tử 3 bit dùng lượng tử Gauss, Q4L.DAT cho 4 bit lượng tử Laplace.
Phương pháp Lloyd Lloyd cũng đưa ra một phương pháp thứ hai cho phép xác định các mức lượng tử mà ông gọi là phương pháp Lloyd I. Phương pháp này có nhiều ưu điểm hơn phương pháp II (giải thuật Lloyd-Max), vì nó dễ dàng cho tính toán và các vector lượng tử hoá có thể mở rộng. Chú ý là vấn đề mà chúng ta quan tâm ở đây là khoảng cách lượng tử hoá, lượng tử hoá của hàm một biến đã biết được phân tán. Vector lượng tử hoá là một vector của nhiều biến mà với các biến này ta đã biết được phân tán.
Thuật toán Lloyd theo các bước sau:
1. Rút ra ước lượng cho phạm vi của các biến
di {i = 0, 1, 2, ..., N}
(Một ước lượng có thể rút ra bằng cách dùng các giá trị từ lượng tử hoá đồng đều hoặc từ các mức lượng tử trước mà ta cần một kết quả tốt hơn).
2. Đặt một biến D1= 0. D1 dùng để lưu lại tình trạng không chính xác lúc trước.
3. Tính
i = 0, 1, ..., N - 1.
4. Tính
i = 0, 1, ..., N - 1
5. Tính tình trạng không chính xác
Có thể dễ dàng mở rộng biểu thức trạng thái không chính xác theo
(13.58)
6. Nếu
thì một giải pháp đã được tìm ra. Lưu lại kết quả và thoát khỏi chương trình.
7. Đặt D1 = D2
8. Quay lại bước 3.
Một chương trình C cho giải thuật trên được đề cập đến ở dưới đây.
Chương trình 13.10 “LLOYDQ.C” Thuật toán Lloyd cho việc thiết kế các mức lượng tử.
/*Program13.10 "LLOYDO.C".Lloyd algorithm for quantizer design.*/
/* Program for designing the Lloyd-quantizer
for a Gauss, uniform or Laplace distribution.*/
#include
#include
#include
#include
#include
double decision_level(double, double);
double f1(double);
double f2(double);
double f3(double);
double p(double);
double Romberg(double, double, double (*)(double));
char ch;
void main( )
{
double *r,*d,step,sum,I1,I2,I3,D1,D2;
int i,j,m,N,xt,yt,niter,ind;
char file_name[16],ch1;
FILE *fptr;
clrscr( ) ;
printf("Enter number of bits--->");
scanf("%d",&m);
printf("Enter number of iterations-->");
scanf("%d",&niter);
N=1<<m;
r=(double *)malloc((N+1)*sizeof(double));
d=(double *)malloc((N+1)*sizeof(double));
printf("Enter choice of distribution:");
printf("\n 1. Gauss.");
printf("\n 2. Uniform.");
printf("\n 3. Laplace. -->(1,2 or 3): ");
while(((ch=getch())!='1')&&(ch!='2')&&(ch!='3'));
putch(ch);
printf("\n Do you wish to start from a previous design?
(y or n)-->");
while(((ch1=getch())!='y')&&(ch1!='n'));
putch(ch1);
ind=0;
switch(ch1)
{
case 'y':
printf("\n Enter name Of file containing quantizer data--Y");
scanf("%s",file_name);
fptr=fopen(file_name,"r");
if(fptr==NULL)
{
printf("\n No such file.");
ind=1;
}
if(ind!=1)
{
for(i=0;i<N;i++)
fscanf(fptr,"%f %f ",&d[i],&r[i]);
d[N]=20.0;
fclose(fptr);
break;
}
case 'n':
d[N]=20.0;
d[0]=-20.0;
if(ch=='2')
{d[N]=1.0; d[0]=-1.0;}
step=(d[N]-d[0])/(float)N;
for(i=1;i<N;i++)
d[i]=d[i-1]+step;
}
D1=0.0;
for(j=0;j<niter;j++)
{
gotoxy(1,15);
printf("iter=%d",j);
D2=0.0;
for(i=0;i<N;i++)
{
I1=Romberg(d[i],d[i+1],f3),
I2=Romberg(d[i],d[i+1],f1);
I3=Romberg(d[i],d[i+1],f2);
r[i]=I2/I3;
D2+=I1-2.0*r[i]*I2+r[i]*r[i]*13;
}
for(i=1;i<N;i++)
d[i]=(r[i]+r[i-1])/2.0;
printf("\nDistortion=%lf",D2);
if(fabs((D2-D1)/D2)<=1.e-6)
{
printf("\n\nConvergence was reached.\n");
break;
}
D1=D2;
}
switch(ch)
{
case '1':
case '3':
printf("\n\n -%c %9.6f ",(char)236,r[0]);
break;
case '2':
printf("\n\n %9.6f %9.6f ", d[0], r[0]);
}
for(i=1; i<N; i++)
printf("\n %9.6f %9.6f ",d[i],r[i]);
switch(ch)
{
case '1':
case '3':
printf("\n +%c",(char)236);
break;
case '2':
printf("\n %9.6f ",-d[0]);
}
printf("\nEnter file name for storing quantizer-->");
scanf("%s",file_name);
fptr=fopen(file_name,"W");
for(i=0; i<N; i++)
fprintf(fptr,"%f %f \n",d[i], r[i]);
fclose(fptr);
}
double p(double x)
{
double prob;
switch(ch)
{
case '1':
prob=exp(-x*x*0.5);
break;
case '2':
prob=1.0;
break;
case '3':
prob=exp(-1.4142136*fabs(x));
break;
}
return prob;
}
/* -------------------------
functions used by numerical integration routine.
-------------------------------- */
double f1(double y)
{
return (y*(p(y)));
}
double f2(double y)
{
return p(y);
}
double f3(double y)
{
return y*y*p(y);
}
Bài tập 13.10
Làm lại bài tập 13.8 nhưng lần này dùng chương trình 13.10 cho giải thuật Lloyd.
2. So sánh thời gian tính toán khi dùng giải thuật Lloyd-Max và khi dùng giải thuật Lloyd.
Từ biểu thức (13.52) và (13.58) chúng ta có thể phát triển một chương trình cho tình trạng méo tối thiểu:
(13.59)
13.6 Lượng tử hoá các hệ số của FCT
Trong phần 13.4 chúng ta đã bắt đầu vấn đề của biến đổi cho mã hoá. Phương pháp chúng ta áp dụng là chia ảnh thành các khối hình vuông; mỗi khối có kích thước 8 ´ 8 và 16 ´ 16. Biến đổi cosin nhanh cho mỗi khối này đã được rút ra. Chúng ta nhận thấy rằng hầu hết các hệ số này có biên độ rất nhỏ so với các giá trị xung quanh khối (một chiều) DC. Câu hỏi đặt ra lúc này là các hệ số nào chúng ta cần lưu giữ và bằng phương pháp nào chúng ta có thể lưu giữ tốt nhất các giá trị này? Câu trả lời cho vấn đề này có thể tìm thấy trong phần lượng tử hoá mà chúng ta đã nghiên cứu ở trên.
Chú ý là các hệ số của FCT xác định một dạng biến dạng. Cho ví dụ, một ảnh có 256 ´ 256 điểm và kích thước của các khối là 8 ´ 8 điểm, có tất cả 64 hệ số cho mỗi khối và 32 ´ 32 khối. Mỗi hệ số có 1024 giá trị khi chúng ta xem xét tất cả các khối, và tạo nên một biến dạng riêng. Đánh giá biến dạng cho hệ số thứ j có thể cho bởi
(13.60)
j = 0, 1, 2, ..., L - 1.
ở đây L là số các hệ số cho một khối và Nj số các mức lượng tử cho hệ số j. Tổng số các biến dạng sẽ là
(13.61)
Làm theo các bước trong phần 13.5 chúng ta được
(13.62)
và (13.63)
Nếu chúng ta coi rằng bất kỳ hệ số nào có thể xác định bằng cùng một hàm khả năng xuất hiện độ sáng, thì thay thế giá trị các hệ số này (mà được biểu diễn trong biểu thức trên là y) bằng
(13.64)
Chúng ta sẽ cho tất cả các hệ số với các phân bố xuất hiện giống nhau, với giá trị trung bình và chuẩn của độ lệch cho bởi m = 0 và s = 1. Kết quả sau khi tính toán cho ta các mức chia và các mức khôi phục cho tất cả hệ số “chia”. Điều này tất nhiên chỉ áp dụng với điều kiện là các hệ số có cùng một số các bit. Trước khi đưa ra các mức lượng tử chúng ta có thể bỏ bớt một số hệ số. Nếu hệ số (0, 0) hay còn gọi là thành phần một chiều DC biểu diễn cho giá trị trung bình của độ sáng của một khối, chúng ta không thể bỏ điểm này đi được. Các hệ số khác trong một khối (còn gọi là các hệ số xoay chiều AC) mang các thông tin về các chi tiết của ảnh. Có thể nhận thấy là các chi tiết có độ lệch lớn hơn độ lệch chuẩn thì mang nhiều tin tức hơn các chi tiết có độ lệch ít hơn độ lệch chuẩn. Vì vậy mà chúng ta bắt đầu lược bỏ các hệ số bắt đầu từ vùng có trải rộng ít nhất. Vậy bao nhiêu hệ số sẽ được chúng ta giữ lại? Điều này phụ thuộc vào mức độ mà chúng ta muốn nén ảnh và phụ thuộc vào bao nhiêu các chi tiết bị mất trên ảnh mà chúng ta có thể chấp nhận được.
Dựa trên các giả thiết trên chúng ta có thể phát triển một thuật toán cho nén ảnh và lượng tử hoá. Các bước sau mô tả cho cả việc lượng tử hoá các hệ số FCT.
1. Tính m và s cho tất cả các hệ số FCT. (Chú ý là độ lệch chuẩn và trung bình có thể tính trong một dải thông của ảnh dùng biểu thức sau cho s:
ở đây xi biểu diễn các giá trị cho một trong các hệ số). m được tính từ tổng của xi.
2. áp dụng các hệ số cho các chi tiết được giữ lại cụ thể là 0.25 , 0.5.
3. Giữ lại các hệ số đã nhân thêm phân số chia có sai lệch cao hơn sai lệch chuẩn.
4. Định dạng một ma trận T có dạng
5. Chia khoảng cách các hệ số cij (cụ thể một cho các giá trị mà Tij = 1) trong tất cả các khối, ngoại trừ các giá trị một chiều cho mỗi khối, như sau:
(13.65)
6. Tính phân bố của các giá trị AC “chia” cho mảng FCT.
7. Tính ss của phân bố rút ra từ các bước trước.
8. Dùng lượng tử hoá Lloyd-Max mức N, và sửa lại các mức chia và khôi phục các mức theo:
(13.66)
i = 0, 1, 2, ..., N - 1.
chú ý dN = -d0
Hàm phân bố Laplace cung cấp một xấp xỉ tốt hơn cho phân bố của các hệ số chia như chúng ta thấy ở phần dưới đây. Sự lựa chọn của N cũng như các hệ số chia của các hệ số phụ thuộc mức độ nén.
9. Lượng tử các hệ số AC “ chia” dùng lượng tử hoá của bước 8.
10. Chia mỗi giá trị một chiều với 2. Điều này đảm bảo rằng các giá trị một chiều không vượt quá 255 (một biểu diễn 8 bit).
11. Định dạng một phần đầu file chứa đầy đủ thông tin để khôi phục lại ảnh bị nén. Phần này chứa thông tin về các mức chia và các mức một chiều bị cắt bớt và tập hợp các giá trị AC cho các hệ số giữ lại.
12. áp dụng mà mã hoá Huffman cho file chứa các giá trị AC.
Ma trận T thường gọi là ma trận khu vực, và cần cung cấp các hệ số cho chức năng khôi phục. Phần đầu file dùng bốn thông tin theo thứ tự sau.
Chiều rộng của khối (ta coi khối là một hình vuông): 1byte.
Số các mức lượng tử: 1 byte.
Chiều rộng của ảnh: 2 byte.
Ma trận T: 1 bit cho một phần tử.
Độ lệch chuẩn, trung bình: 4 byte trong biểu diễn số thực ´ số các hệ số.
Các mức khôi phục: 4 byte trong biểu diễn số thực/ mức.
Chương trình 13.11 đưa ra lượng tử của các hệ số biến đổi cosin dùng các bước trên. Kết quả của chương trình là 3 file có cùng tên do người dùng xác định, nhưng có phần mở rộng khác nhau. Nếu người dùng cho tên ‘image1q’ khi trả lời cho thông báo “Enter file name to store quantized image ...” xuất hiện khi chạy chương trình, ba file sau đây được tạo ra:
‘image1q.hdr’, là file chứa thông tin về header.
‘image1q.dc’ là file chứa thông tin về các hệ số một chiều bị cắt bớt của một loạt các khối.
‘image1q.ac’ file chứa một loạt các lượng tử cho các hệ số AC của một loạt các khối. Phần này không chứa thông tin về các giá trị bị loại bỏ, nếu chúng ta coi rằng nó có giá trị 0 và ma trận vùng chỉ ra vị trí của nó trong khối.
Chương trình nhập các tên sau đây do người dùng đặt:
Tên file chứa biến đổi cosin của một ảnh bị nén, và kích thước của khối dùng trong FCT.
Tên file chứa dữ liệu lượng tử Lloyd-Max. Nếu bạn chạy MAXQ1.EXE (xem bài tập 13.8), và cho tên file chứa các mức lượng tử mà bạn muốn dùng. Lượng tử 5 bit Laplace cho kết quả tuyệt vời. Lượng tử 4 bit Laplace cho bạn một kết quả có thể chấp nhận được, như là một kinh nghiệm cho bạn.
Phần trăm của các hệ số AC được giữ lại. Bạn có thể lựa chọn bất kỳ giá trị nào: thông thường là 0.25 hoặc là 0.5. Các phân số nhỏ hơn hay mức thấp hơn các mức lượng tử gây nên một sai số cao hơn khi khôi phục lại ảnh.
Tên file chứa các lượng tử hoá của ảnh. Để trả lời, bạn cho tên file không có phần mở rộng như đã nói trước đây.
Chương trình 13.11 “QUANTIZE.C” Chương trình cho lượng tử hoá và nén kết quả của biến đổi FCT.
/*---------------------------------
This program does quantization on the block
DCT of an image . The input file should be the
block DCT transform o the image. The program
also designs a quantizer based on the statistics
of the DCT coefficients.
You should first run either the program for
the Lloyd-Max or the one for Lloyd quantizer
design if you haven't already done so.
-----------------------------------------*/
#include
#include
#include
#include
#include
#include
#include
void main()
{
int N,i,j,k,M,xt,yt,N1,NS,k1,k2,NB,Nt;
int m,kk,kk1,kk2,kt;
unsigned long int NC,NQ,loc,loc1;
float nsq,*sigma,*d,sum,*histo,*x,*mu;
float denom,mut,sigmat,*sumd,*sumdsq;
float Max,Min,h,*r,dt,rt,fdct,smax;
int b,*T;
FILE *fptri,*fptro,*fptro1;
char file_name[14],*imaget,temp[14];
float *buffi;
unsigned char ch;
clrscr();
printf("This program carries out quantization on the block");
printf("\n DCT of an image. The 2-D FCT should be run on the");
printf("\n image prior to this program. If you haven't, just");
printf("\n enter NULL to the first prompt this will return you");
printf("\n to DOS.\n");
printf("\n You will need also to run the program MAXQ for");
printf("\n designing the Max quantizer. The program will");
printf("\n modify the decision and reconstruction levels of");
printf("\n the quantizer to fit the results to the statistics");
printf("\n of the DCT data.\n\n");
printf("Enter name of input file containing 2-D FCT results-->");
scanf("s",file_name);
fptri=fopen(file_name,"rb");
if(fptri==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
nsq=filelength(fileno(fptri));
/*---------------------
Assume image is square.
-----------------------*/
N=(int)sqrt((double)(nsq/sizeof(float)));
m=(int)(log10((double)N)/log10((double)2.0));
k=1;
for( i=0; i<m; i++)
k<<=1;
if(k!=N)
{
printf ("\nTransformed Image should have dimensions" ) ;
printf ( "\nthat are multipies of 2. Check size of \n" ) ;
printf ( "original image used in the transform.");
exit(1);
}
printf("\nEnter block size used (e.g. 8x8, 16x16, etc.)-->");
scanf("%dx%d", &NB,&NB);
Nt=NB*NB ; /* size of block.*/
N1=N/NB; /* # of blocks in the horizontal or
vertical directions*/
NS=NB*N ; /* Size of input buffer. */
NC=N1*N1; /* Total # of blocks in image.
/*Allocating memory.*/
buffi=(float *)malloc(NS*sizeof(float));
sumd=(float *)malloc(Nt*sizeof(float));
sumdsq=(float *)malloc(Nt*sizeof(float));
sigma=(float *)malloc(Nt*sizeof(float));
mu=(float *)malloc(Nt*sizeof(float));
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(WHITE+(GREEN<<4)+BLINK);
cputs( "WAIT" );
gotoxy(xt,yt);
for(j=0;j<Nt ;j++)
{
sumd[j]=0.0;
sumdsq[j]=0.0;
}
for(i=0; i<N1; i++)
{
fread( buffi, sizeof(float),NS, fptri);
for(j=0;j<N1 ;j++)
{
kk=j*NB ;
for(k1=0;k1<NB;k1++)
{
loc1=kk+k1*N;
kk2=k1*NB;
for(k2=0;k2<NB;k2++)
{
loc=loc1+k2;
k=kk2+k2 ,
sumd[k]+=buffi[loc];
sumdsq[k]+=buffi[loc]*buffi[loc];
} /* k2 loop. */
} /* k1 loop. */
} /*j-loop.*/
} /* i-100p*/
/* Compute s and g for each coefficient.*/
printf("\n Computing the Standard deviation for ");
printf("each coefficient.\n");
denom=(float)NC*(float)(NC-1.0);
for(i=0;i<Nt;i++)
{
sigma[i]=fabs((((float)NC)*sumdsq[i]-(sumd[i]*sumd[i]))/denom);
sigma[i]=(float)sqrt((double)sigma[i]);
mu[i]=sumd[i]/(float)NC;
}
free(sumd);
free(sumdsq);
rewind(fptri); /*Rewind input file.*/
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(WHITE+(BLACK<<4));
cputs(" ");
gotoxy(xt,yt);
printf("\nDo you wish to save histogram data");
printf("\n for plotting (y or n)-->");
while(((ch-getch())!='y')&&(ch!='n'));
if(ch=='y')
{
printf("\nEnter name for output file to store histogram");
printf("\n of a.c components of the DCT of each block-->");
scanf("%s",file_name);
fptro=fopen(file_name,"w");
printf("\nEnter number of data points you wish generated on the");
printf("\n histogram distribution-->");
scanf("%d",&M);
}
else M=64;
/*Allocating memory.*/
d=(float *)malloc((M+1)*sizeof(float));
histo=(float *)malloc(M*sizeof(float));
x=(float *)malloc(M*sizeof(float));
clrscr ();
gotoxy(70, 25);
textattr(WHITE+(GREEN<<4)+BLINK);
cputs("WAIT");
gotoxy (1,1);
printf ("\nGenerating the distribution of the a.c. " );
printf ( "\n coefficients of the Cosine transform. ");
/*scaling and generating histogram of
a.c. coeff i ci ents . */
Max=1.e-20;
Min=1.e20 ;
for (i=0; i<M; i++)
histo[i]=0.0;
for(i=0; i<N1; i++)
{
fread(buffi,sizeof(float),NS,fptri);
for(j=0;j<N1;j++)
{
kk=j*NB ;
for(k1=0;k1<NB;k1++)
{
loc1=kk+k1*N;
kk2=k1*NB;
for(k2=0;k2<NB;k2++)
{
loc=loc1+k2;
k=kk2+k2 ;
buffi[loc]-(buffi[loc]-mu[k])/sigma[k];
if(buffi[loc]<Min) Min=buffi[loc];
if(buffi[loc]>Max) Max=buffi[loc];
}/* k2 loop. */
}/* k1 loop.*/
}/* j-loop.*/
}/* i-loop.*/
rewind(fptri); /*Rewind input file.*/
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(WHITE+(BLACK<<4));
cputs ( " ");
gotoxy(xt,yt+1);
printf("\n\nEnter fraction of DCT coefficients to be");
printf("\n retained per block (e.g. 0.25, 0.5, etc.) --->");
scanf ( "%f", &fdct);
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(WHITE+(GREEN<<4)+BLINK);
cputs("WAIT");
gotoxy(xt,yt);
k=fdct*Nt;
T=(int *)malloc(Nt*sizeof(int));
sum=0.0;
for(i=0;i<Nt;i++)
T[i]=0;
T[0]=1; /* dc component.*/
while(k>1)
{
smax=-1.0;
for(i=0;i<Nt;i++)
{
if(T[i]==1) continue;
if(sigma[i]>smax)
{
smax=sigma[i];
loc=i;
}
}
T[loc]=1;
k--;
}
h=(Max-Min)/(float)m;
d[0]=Min;
for(i=0;i<(m-1);i++)
d[i+1]=d[i]+h;
d[M]=Max;
for(i=0;i<N1;i++)
{
fread(buffi,sizeof(float),NS,fptri);
for(j=0;j<N1;j++)
{
kk=j*NB ;
for(k1=0;k1<NB;k1++)
{
loc1=kk+k1*N;
kk2=k1*NB ;
for(k2=0;k2<NB;k2++)
{
loc=loc1+k2;
k=kk2+k2 ;
if(T[k]==0) continue;
buffi[loc]=(buffi[loc]-mu[k])/sigma[k];
for(kt=0;kt<M;kt++)
{
if((buffi[loc]=d[kt]))
{
histo[kt]+=1.0;
break;
}
}
} /* k2 loop. */
} /* k1 loop.*/
}/* j-loop.*/
} /* i-loop.*/
rewind(fptri ) ; /*Rewind input file. */
for( i=0; i<M; i++)
x[i]=(d[i]+d[i+1])/2.0;
NQ=(long int)(Nt-1)*(long int)(N1*N1);
if (ch=='y' )
{
printf ( "\n\nPrepari ng data for output plotting.");
/* Preparing data for X-Y plotting.
Data is prepared in a format suitable for plotting by GRAFTOOL,
a graphical analysis software tool by 3-D Visions ,
(412 S.Pacific Coast Hwy. Suite 201, Redondo Beach, CA. 90277).*/
fprintf(fptro,"%d %d\n",M,2);
for(i=0;i<M;i++)
{
fprintf(fptro,"%e %e\n",x[i],histo[i]);
}
fclose(fptro);
}
xt=wherex();
yt=wherey();
gotoxy(70, 25);
textattr(WHITE+(BLACK<<4));
cputs ( " " );
gotoxy (xt, yt+1);
free (d);
printf(" \nDesigning the quantizer.");
mut=0.0;
sigmat=0.0;
for(j=0; j<M; j++)
mut+=histo[j]*x[j];
sum=0.0;
for(j=0; j<M ; j++)
sum+=x[j]*x[j]*histo[j];
sigmat= NQ*sum-mut*mut;
sigmat/=(float)NQ*(float)(NQ-1);
sigmat=(float)sqrt((double) sigmat);
printf("\n Enter file name containing the");
printf("\n Lloyd-Max quantizer--> ");
scanf("%s",file_name);
fptro=fopen(file_name,"r");
if(fptro==NULL)
{
printf("\n That file does not exist. Run the ");
printf("\n program for designing the Lloyd-Max quantizer.\n");
exit(1);
}
i=0;
while(1)
{
k=fscanf(fptro," %f %f",&dt,&rt);
i++;
if(k==EOF) break;
}
NQ=i-1;
d=(float *)malloc(NQ*sizeof(float));
r=(float *)malloc(NQ*sizeof(float));
rewind(fptro);
i=0;
while(1)
{
k=fscanf(fptro," %f %f",&d[i],&r[i]);
if(k==EOF) break;
d[i]=d[i]*sigmat;
r[i]=r[i]*sigmat;
i++;
}
fclose(fptro); /*Close file containing distribution.*/
/* Header information:
Block size, image size, number of quantization levels,
zonal sampling matrix, standard deviations & means,
reconstruction levels. */
printf("\nEnter file name to store quantized image,\n ");
textattr(BLUE+(YELLOW<<4)+BLINK);
cputs(" Do NOT give the name an extension --> ");
scanf("%s",file_name);
textattr(WHITE+(BLACK<<4));
imaget=strchr(file_name,'.');
if(imaget!=NULL)
{
j=strcspn(file_name,".");
file_name[j]='\0';
}
imaget=strcpy(temp,file_name);
imaget=strcat(temp,".hdr");
fptro=fopen(temp,"wb");
/*save size of block,
size of image and # of quant. levels .*/
putc(NB,fptro);
putc(NQ,fptro);
k=(255)&N;
putc(k,fptro);
k=(255)&(N>>8);
putc(k,fptro);
/* Zonal Sampling matrix. */
j=0; k=0; loc=0;
for(i=0;i<Nt;i++)
{
j=(j<<1)|T[i];
k++;
if(k==8)
{
putc(j,fptro);
k=0;
}
}
/* standard deviations and means.*/
for(i=0; i<Nt; i++)
{
k=i<<1;
buffi[k]=sigma[i];
buffi[k+1]=mu[i];
}
/* reconstruction levels.*/
k=Nt<<1;
for(i=0; i<NQ; i++)
{
buffi[k+i]=r[i];
}
fwrite(buffi , sizeof(float), (Nt<<1)+NQ,fptro);
fclose(fptro);
clrscr();
printf ( " \n Quantizing.");
imaget=strcpy(temp, file_name);
imaget=strcat(temp, ".DC");
fptro1=fopen(imaget,"wb");
imaget=strcpy(temp, file_name);
imaget=strcat(temp,".AC");
fptro=fopen(imaget,"wb");
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(RED+(LIGHTGRAY<<4)+BLINK);
cputs("WAIT");
gotoxy(xt,yt);
for(i=0;i<N1;i++)
{
fread(buffi,sizeof(float),NS,fptri);
for(j=0;j<N1;j++)
{
kk=j*NB;
for(k1=0;k1<NB;k1++)
{
kk1=k1*N;
kk2=k1*NB;
for(k2=0;k2<NB;k2++)
{
if((k1+k2)==(int)0)
{
ch=(unsigned char)(buffi[kk]/2.0);
putc(ch,fptro1);
continue;
}
loc=kk1+k2+kk;
loc1=kk2+k2;
if(T[loc1]==0) continue;
else
buffi[loc]=(buffi[loc]-mu[loc1])/sigma[loc1];
if(buffi[loc]<d[1]) ch=(char)0;
else if(buffi[loc]>=d[NQ-1]) ch=(char)(NQ-1);
else
{
for(k=1;k<(NQ-1);k++)
{
if((buffi[loc]=d[k]))
{
ch=(char)k;
break;
} /* if */
} /* k-loop. */
} /* else */
putc(ch,fptro);
} /* k2 loop. */
}/* k1 loop. */
} /* j-loop. */
} /* i-loop. */
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(WHITE+(BLACK<<4));
cputs ( " ");
gotoxy(xt,yt);
printf(" \nDone !!");
fcloseall();
}
Trước khi bạn dùng chương trình 13.11 bạn cần thiết kế lượng tử hoá Lloyd-Max như các bước trong bài tập 13.8. Bạn cũng sẽ cần chạy FCT2D.EXE trên ảnh đã bị nén.
Bài tập 13.11
Bạn sẽ tìm thấy trên đĩa đi kèm theo cuốn sách này ảnh “KAREN.IMG”. Đưa ra các bước sau:
Chạy 2-D FCT trên ảnh “KAREN.IMG” dùng một khối 8 ´ 8. File FCT này có tên là “KARENFCT.IMG”.
Chạy chương trình “QUANTIZE.EXE” trên “KARENFCT.IMG” dùng các bước sau:
Tỷ lệ các hệ số giữ lại là 0.25.
Lượng tử hoá 32 mức (5 bit). Chứa ảnh đã lượng tử hoá và nén trong “KarenQ”.
3. Chạy chương trình mã hoá Huffman trên “KarenQ.AC”.
Bài tập 13.11 cho kết quả sau đây:
Kích thước ảnh gốc: 65536 byte.
Kích thước file chứa các điểm AC: 15360 byte.
Kích thước file chứa các điểm DC: 1024 byte.
Kích thước của header: 652 byte.
Kích thước của file mã hoá Huffman: 8973 byte.
Tỷ lệ nén thường được tính theo bit trên điểm ảnh (bpp). Giá trị này có thể tính từ các biểu thức trên theo:
ằ 1.3
Giá trị này chứng tỏ rằng 83.75 phần trăm kích thước của ảnh gốc đã được nén lại. Bạn có thể giảm kích thước của file chứa các điểm DC một chút dùng RLC. Nếu bạn làm như vậy kích thước của file này sẽ là 1017 bit. Mã RLC có thể cho kết quả tốt trên ảnh có kích thước lớn.
Nếu bạn dùng lượng tử hoá 4 bit thì tỷ lệ nén là 1.1 bpp.
Nếu bạn trả lời “y” để tính dữ liệu lược đồ mức xám cho các khối AC (chỉ các dứ liệu còn giữ lại), bạn sẽ có được file dữ liệu dùng để in ra. Kết quả dữ liệu lược đồ mức xám của ảnh “KARENFCT.IMG” cho trong hình 13.19. Lược đồ mức xám này biểu diễn phân bố Laplace như tôi đã đề cập ở phần đầu cuốn sách này.
Bài tập 13.12
Làm lại bài tập 13.10 trên ảnh “IKRAM.IMG”.
Tính bpp.
Vẽ lược đồ mức xám của khoảng cách các hệ số AC.
Để nghiên cứu ảnh lượng tử hoá bạn cần hiện tất cả các khối một lúc. Chương trình sau thực hiện chính xác điều đó. Chương trình này sẽ in các dữ liệu không giữ lại là các điểm zero.
Chương trình 13.12 “QUANCHK.C” Chương trình hiện thị các khối lượng tử hoá.
/*Program 13.12 "QUANCHK.C". Program to display quantized blocks.*/
/* This program displays the quantized file block by block.
Values that were not retained by the quantizer
program are typed out as zeros. */
#include
#include
#include
#include
#include
#include
Hình 13.19 Lược đồ mức xám của các hệ số AC đã được chia.
void main()
{
int i,j,k1,k2,N,N1,NB,NS,Nt,*T,k,kk;
unsigned char *buff;
FILE *fptr,*fptr2;
char file_name[14],ch,*imaget,temp[14];
clrscr();
printf("\n");
printf("Enter file name in which you stored quantized image,\n");
textattr( BLUE+(YELLOW<<4)+BLINK);
cputs( " Do NOT give any extension - -> " ) ;
scanf( "%s" , file_name);
textattr(WHITE+(BLACK<<4));
imaget=strchr(file_name,'.');
if( imaget!=NULL)
{
j=strcspn(file_name,".");
file_name[j]='\0';
}
imaget=strcpy(temp,file_name);
imaget=strcat(temp,".hdr");
fptr=fopen(temp,"rb");
if(fptr==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
NB=getc(fptr);
N=getc(fptr);
N=getc(fptr)+(getc(fptr)<<8);
Nt=NB*NB ;
NS=NB*N ;
N1=N/NB ;
buff=(unsigned char *)malloc(NS*sizeof(char));
T=(int *)malloc(Nt*sizeof(int));
kk=Nt/8 ;
fread(buff,sizeof(char),kk,fptr);
/* Zonal Sampling matrix. */
k=0;
for(i=0;i<kk;i++)
{
for(j=0;j<8;j++)
{
T[k]=buff[i]&(128);
T[k] >>=7;
k++;
buff[i]<<=1;
}
}
fclose(fptr);
imaget=strcpy(temp,file_name);
imaget=strcat(temp,".AC");
fptr=fopen(temp,"rb");
if(fptr==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
imaget=strcpy(temp,file_name);
imaget=strcat(temp,".DC");
fptr2=fopen(temp,"rb");
if(fptr2==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
clrscr();
gotoxy(1,20);
printf("Press any key to continue, or ESC to exit.");
gotoxy(1,1);
for(i=0;i<N1;i++)
{
for(j=0;j<N1;j++)
{
for(k1=0;k1<NB;k1++)
{
for(k2=0;k2<NB;k2++)
{
if(T[k1*NB+k2]==0) printf(" %3d",0);
else
{
if((k1+k2)==0) printf(" %3d",getc(fptr2)*2);
else printf(" %3d",getc(fptr));
}
}
printf("\n");
}
ch=getch();
if(ch==27) exit(1);
clrscr();
gotoxy(1,20);
printf("Press any key to continue, or ESC to exit.");
gotoxy(1,1);
}
}
}
Thuật toán khôi phục lượng tử hoá được mô tả bởi chương trình C sau đây. Chú ý là biểu thức khôi phục được mô tả bằng
(13.67)
ở đây là các hệ số lượng tử.
Chương trình 13.13 "DEQUANTZ.C" Chương trình đảo ngược lượng tử hoá của tệp đã được lượng tử hoá trước.
/*Program 13.13 "DEQUANTZ.C". Program to dequantize a
previously quantized file.*/
/* This program restores a file quantized
by program 'QUANTIZE.EXE' . The output file
can be used as input to IFCT2D.EXE to obtain
the original image.*/
#include
#include
#include
#include
#include
#include
#include
void main()
{
int i,j,NB,N,N1,NS,NQ,xt,yt,k;
int kj,k1,k2,kk1,kkb,loc,ch,kk2,Nt,kk,NQ1;
float *buffo,*sigma,*r,*mu;
char file_name[14],*imaget,temp[14];
unsigned char *buffi;
int *T;
FILE *fptri,*fptro,*fptri2;
clrscr();
printf("\nEnter file name in which you stored quantized image,\n ");
textattr(BLUE+(YELLOW<<4)+BLINK);
cputs(" Do NOT give any extension --> ");
scanf("%s",file_name);
textattr(WHITE+(BLACK<<4));
imaget=strchr(file_name,'.');
if(imaget!=NULL)
{
j=strcspn(file_name,".");
file_name[j]='\0';
}
imaget=strcpy(temp,file_name);
imaget=strcat(temp,".hdr");
fptri=fopen(temp,"rb");
if(fptri==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
NB=getc(fptri);
NQ=getc(fptri);
N=getc(fptri)+(getc(fptri)<<8);
Nt=NB*NB;
NS=NB*N ;
N1=N/NB;
buffi=(unsigned char *)malloc(NS*sizeof(char));
buffo=(float *)malloc(NS*sizeof(float));
T=(int *)malloc(Nt*sizeof(int));
sigma=(float *)malloc(Nt*sizeof(float));
mu=(float *)malloc(Nt*sizeof(float));
r=(float *)malloc(NQ*sizeof(float));
kk=Nt/8;
fread(buffi,sizeof(char),kk,fptri);
/* Zonal Sampling matrix. */
k=0 ;
for(i=0;i<kk;i++)
{
for(j=0;j<8;j++)
{
T[k]=buffi[i]&(128);
T[k]>>=7;
k++;
buffi[i]<<=1;
}
}
/* standard deviations and means.*/
fread(buffo,sizeof(float),(Nt<<1)+NQ,fptri ) ;
for(i=0; i<Nt; i++)
{
k=i<<1;
sigma[i]=buffo[k];
mu[i]=buffo[k+1];
}
k=Nt<<1;
/* reconstruction levels.*/
for(i=0; i<NQ ; i++)
{
r[i]=buffo[k+i];
}
fclose(fptri);
imaget=strcpy(temp, file_name);
imaget=strcat(temp, ".AC");
fptri=fopen(temp,"rb");
if(fptri==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
imaget=strcpy(temp,file_name);
imaget=strcat(temp,".DC");
fptri2=fopen(temp,"rb");
if(fptri2==NULL)
{
printf("\nNo such file exists.\n");
exit(1);
}
printf("\nEnter file name to store dequantized file-->");
scanf("%s",file_name);
fptro=fopen(file_name,"wb");
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(RED+(LIGHTGRAY<<4)+BLINK);
cputs("WAIT");
gotoxy(xt,yt);
NQ1=NQ>>1;
for(i=0;i<N1;j++)
{
for(j=0;j<N1;j++)
{
kj=j*NB;
for(k1=0;k1<NB;k1++)
{
kk1=k1*N;
kkb=k1*NB;
for(k2=0;k2<NB;k2++)
{
if((k1+k2)==0)
{
buffo[kj]=(float)getc(fptri2)*2.0;
continue;
}
else
{
loc=kk1+k2+kj;
kk2=kkb+k2;
if(T[kk2]==0) ch=NQ1;
else ch=getc(fptri);
buffo[loc]=r[ch]*sigma[kk2]+mu[kk2];
}
}
}
}
fwrite(buffo,sizeof(float),NS,fptro);
}
fcloseall();
xt=wherex();
yt=wherey();
gotoxy(70,25);
textattr(WHITE+(BLACK<<4));
cputs(" ");
gotoxy(xt,yt);
}
Bài tập 13.13
Khôi phục ảnh KARENQ rút ra từ bài tập 13.12.
áp dụng IFCT.EXE trên file đã khôi phục.
Hiển thị kết quả và so sánh với ảnh gốc.
Làm lại bài tập 13.11 và làm lại các bước từ 1 đến 3 của bài tập này, nhưng lần này dùng:
Tỷ lệ các hệ số giữ lại là 0.25. 3 bit lượng tử hoá.
Tỷ lệ các hệ số giữ lại là 0.15. 3 bit lượng tử hoá.
Hệ số các chi tiết giữ lại là 0.1. 5 bit lượng tử hoá.
So sánh kết quả rút ra từ bước 4.
Bài tập 13.14
Làm lại bài tập 13.13 trên ảnh IKRAM.IMG.
Một phương pháp nén ảnh khác bằng mã hoá ảnh là phương pháp định vị bit. Với phương pháp này bạn không cần mã hoá Huffman. Phương pháp này dựa trên sự định vị số khác nhau của các bit cho hệ số dựa trên giá trị đóng góp của các hệ số này khi khôi phục lại ảnh từ các giá trị FCT của nó. Cũng như các chú ý từ phần trên, các hệ số với sai lệch cao hơn sai lệch chuẩn thì sẽ mang nhiều tin tức hơn các hệ số có sai lệch ít hơn sai lệch chuẩn. Vì vậy số bit cho các hệ số phụ thuộc vào sai lệch chuẩn của nó. Nếu bpp được cho, tổng số các bit xác định cho khối được cho bởi
B = bpp ´ NB ´ NB (13.68)
ở đây NB ´ NB là kích thước của khối ảnh; và bit quy đổi được cho là
(13.69)
ở đây bij là số các bit quy đổi cho hệ số (i,j) và q là một biến được thay đổi cho đến khi biểu thức sau thoả mãn:
(13.70)
ảnh nén.
ảnh nén.
ảnh.
2-D FCT
Lượng tử hoá
Mã hoá Huffman cho hệ số AC.
RLC của hệ số DC.
Giải mã Huffman hệ số AC.
Giải mã RLC của hệ số DC.
Lượng tử hoá ngược.
Đảo 2-D FCT
ảnh giải nén.
Hình 13.20 Trình bày các bước thực hiện mã hoá nén ảnh và giải nén ảnh.
Một thuật toán dựa trên giả thiết một nửa quãng cho thiết kế một ma trận xác định bit, mà ma trận
Các file đính kèm theo tài liệu này:
- CHUONG13-Nén dữ liệu ảnh.DOC