Nhập môn Xử lý ảnh - Chương 2: Các hệ thống hai chiều

chương 2 Các hệ thống hai chiều 2.1 Chỉ dẫn Các tín hiệu vốn là hai chiều thì phải được xử lý bằng kỹ thuật xử lý tín hiệu hai chiều. Các điểm ảnh được coi là một dạng tín hiệu hai chiều, vì vậy việc xem xét các kiến thức cơ bản của hệ thống hai chiều (2-D) là cần thiết. Trong chương này chúng ta sẽ xem xét biểu diễn tần số và biểu diễn không gian của tín hiệu hai chiều, và cung cấp kiến thức cần thiết cho việc phát triển các kỹ thuật xử lý ảnh hai chiều như là nổi ảnh, trơn ảnh và khôi phục ảnh. 2.2 Các tín hiệu hai chiều Một tín hiệu lấy mẫu hai chiều được biểu diễn dưới dạng một mảng hai chiều x(n1TV,n2TH) (V: Vertical - chiều dọc, H: Horizontal - chiều ngang). ở đây n1, n2 là số nguyên và TV, TH là khoảng cách lấy mẫu theo chiều ngang và chiều dọc. x(n1TV, n1TH) thường được ký hiệu tắt là x(n1,n2) và biểu diễn cho cường độ sáng hay là biên độ tín hiệu tại điểm (n1,n2) trong miền không gian. Hình 2.1 giới thiệu một phép biểu diễn trong miền không gian, ở đây cường độ sáng của tín hiệu được ký hiệu bởi chiều cao tại điểm x(n1TV, n2TH). Chú ý rằng TV là khoảng cách lấy mẫu giữa hai cột dọc kế tiếp nhau, TH là khoảng cách lấy mẫu giữa hai dòng ngang kế tiếp nhau của tín hiệu 2-D. Mặc dù một tín hiệu lấy mẫu 2-D có thể xử lý như các dãy của tín hiệu lấy mẫu một chiều bằng cách xử lý tất cả các hàng (cột) một cách tuần tự song cách tiếp cận này không cho một kết quả mong đợi như khi xử lý hai chiều. Ví dụ, nếu chúng ta dùng một bộ lọc làm nổi đường biên ảnh 1-D, cụ thể đó là một bộ lọc thông cao, trên một ảnh bằng cách xử lý từng hàng một, thì đường biên sẽ phần lớn được làm nổi bật dọc theo các đường thẳng đứng. Các đường biên ảnh nằm theo các đường nằm ngang sẽ không được làm nổi một chút nào và các đường biên nằm theo các hướng khác ngoài hai hướng này sẽ nhận được hiệu ứng làm nổi ảnh ít hơn các đường biên dọc. Để đạt được hiệu quả như nhau theo mọi hướng, tín hiệu được lấy mẫu hai chiều phải được xử lý qua một hệ thống 2-D (Hình 2.2). Trong hệ thống tuyến tính bất biến - TTBB (Linear Shift Invariant - LSI), đáp ứng đầu ra có thể tính theo công thức : (2.1) Dấu * được hiểu là tích chập và h(n1,n2) là đáp ứng xung của hệ thống 2-D. Biểu thức (2.1) có thể viết là: (2.2) n1Tv 2Tv Tv 2TH TH x(n1,Tv,n2,TH) n2TH Hình 2.1 Biểu diễn trong miền khoảng cách. 2.3 Một số dãy 2-D thông dụng Chúng bao gồm: 1. Dãy xung đơn vị : (2.3) 2. Dãy nhảy bậc đơn vị : (2.4) 3. Dãy hàm mũ: (2.5) 4. Dãy tín hiệu hình sin (phức): -Ơ <n1,n2< +Ơ (2.6) h(n1,n2) x(n1,n2) y(n1,n2) Hình 2.2 Xử lý tín hiệu 2-D. 2.4 Đáp ứng tần số của hệ thống 2-D -TTBB Đặt Đáp ứng ra có thể rút ra khi dùng biểu thức (2.2). (2.7) hoặc (2.8) Công thức này có thể viết lại thành Tín hiệu ra là tín hiệu hình sin phức (sinusoid) hoàn toàn có cùng tần số như tín hiệu vào, nhưng biên độ và góc pha thì bị thay đổi bởi hàm khuyếch đại phức H(w1,w2). Hàm khuếch đại này gọi là đáp ứng tần số và được cho bởi (2.9) Biểu thức được gọi là nhân. Nếu khoảng cách cách lấy mẫu TV,TH đã được biết thì biểu thức (2.9) có thể viết lại thành (2.10) w1, w2 có thứ nguyên là radian/đơn vị, còn u và v có thứ nguyên là vòng/đơn vị. Đơn vị ở đây có thể là đơn vị khoảng cách (như cm, inch) hoặc là đơn vị thời gian (như giây). Việc chọn đơn vị (thời gian hoặc khoảng cách) phụ thuộc nguồn gốc của ảnh, đó là một phép chiếu từ không gian ba chiều lên mặt phẳng hai chiều. Nếu ta xử lý với một ảnh lấy ra trực tiếp từ ma trận CCD camera thì TV và TH (và do đó là đơn vị) phải tính theo chiều không gian (xem hình 2.3). Mặt khác, với một ảnh truyền hình thì TV và TH phải theo chiều thời gian (xem hình 2.4). Từ (2.9) ta có thể viết (2.11) Và từ (2.10) ta có thể viết (2.12) TV TH Hình 2.3 TV và TH cho lấy mẫu ảnh trên một ma trận camera CCD. Hình 2.4 TV và TH cho một ảnh quét xen kẽ. Hàm H(w1,w2) xác định trên toàn bộ miền và là hàm tuần hoàn trong miền tần số với chu kì tuần hoàn là 2p đối với w1 và w2. H(u,v) xác định trên miền và là hàm tuần hoàn với chu kì 1/TV và 1/TH cho u và v. Có thể chiếu H(w1, w2) hoặc H(u, v) lên miền chuẩn hoá, ở đây w/1, w/2 bằng cách đặt w/1=w1/p; w/2=w2/p hoặc w/1=2uTV; w/2=2vTh . w/1 và w/2 gọi là tần số chuẩn hoá, hàm H(w/1,w /2 ) có thể viết lại (2.13) Nếu chúng ta hạn chế h(n1,n1) chỉ lấy các giá trị thực thì đáp ứng tần số thoả mãn: (2.14) H* = liên hợp phức của H. Điều này dẫn đến H(w1,w2) đối xứng (Hình 2.5). A B B* A* w1 w2 Hình 2.5 Đối xứng tâm. Chú ý rằng nếu x(n1,n2) = d(n1,n2), thì biểu thức (2.2) trở thành y(n1,n2) = h(n1,n2). Vì lý do này mà h(n1,n2) được gọi là đáp ứng xung, hoặc là đáp ứng biên độ, của hệ thống 2-D. Bài tập 2.1 Tính biểu thức đáp ứng tần số của một hệ thống với đáp ứng xung cho bởi Chứng minh rằng công thức tính đáp ứng tần số có thể tách được. 2.5 Tính đáp ứng xung từ đáp ứng tần số Đáp ứng tần số của h(n1,n2) được cho bởi : (2.15) Xét tích phân (2.16) Thay biểu thức (2.15) vào biểu thức (2.16) chúng ta được Và có thể viết thành Và biến đổi thành Điều này có nghĩa là đáp ứng xung có thể tính từ đáp ứng tần số qua mối quan hệ: h(n1,n2) = (2.17) Nếu đáp ứng tần số được cho dưới dạng hàm của u,v (vòng/đơn vị), thì biểu thức (2.17) có thể viết thành (2.18) Hoặc cho tần số chuẩn hoá: (2.19) Ví dụ 2.3 Cho đáp ứng tần số (xem hình 2.10), hãy tính đáp ứng xung. w1 a -a b b p p -p -p w2 Hình 2.10 Ví dụ 2.3. Giải Từ phương trình (2.17) chúng ta có thể viết : Bởi vì đáp ứng tần số là hàm tách được của hai biến và nên đáp ứng xung cũng là một hàm hai biến tách được. Khái niệm “tách được” ở đây nghĩa là có thể phân tích h(n1,n2) = f1(n1).f2(n2). Ví dụ 2.4 Tìm đáp ứng xung của một bộ lọc thông thấp đối xứng vòng tròn lý tưởng được mô tả như sau (xem hình 2.11 và 2.12): Giải Có thể dễ dàng thấy nếu là một hàm đối xứng vòng tròn lý tưởng, cụ thể là thì cũng là một hàm tuần hoàn đối xứng vòng tròn, tức là . Vì vậy cách dễ dàng nhất để tìm là tìm h(n1, 0) và hàm theo n1. Chúng ta rút ra từ: Ta có hoặc R -R p p -p -p w2 w1 Hình 2.11 Ví dụ 2.4. Hình 2.12 Ví dụ 2.4. Biểu thức này có thể biểu diễn thành ở đây J1(x) = hàm Bessel loại 1. Vì thế : Không phải lúc nào cũng dễ dàng rút ra được một biểu thức phân tích đáp ứng xung như trên. Ví dụ dưới đây minh hoạ việc sử dụng phép tích phân số để thu được . Ví dụ 2.5 Tính đáp ứng xung của bộ lọc thông thấp Butterworth đối xứng vòng tròn có đáp ứng tần số cho bởi ở đây và D0 = 3-dB. (Giả thiết rằng ). Dùng biểu thức trong công thức 2.17 để tính đáp ứng xung với Giải Vì là thực và đối xứng vòng tròn, công thức (2.17) có thể viết đơn giản lại là đơn giản hơn nữa ta có thể viết (2.20) Lời giải của bài toán này sẽ được đưa ra thông qua một chương trình viết bằng ngôn ngữ C. Ngôn ngữ C được sử dụng rộng rãi bởi tính linh hoạt của nó. C phù hợp cho các ứng dụng khác nhau. Trước khi bạn viết chương trình, bạn cần nghiên cứu phương pháp tính tích phân kép. Sau đây ta sẽ phát triển thêm qui tắc Simpson để tính tích phân kép. Qui tắc tính tích phân của Simpson được viết như sau: ở đây và m là một số lẻ. Với tích phân kép Sử dụng qui tắc Simpson ta có thể viết như sau: vậy có thể viết Lại áp dụng qui tắc Simpson vào công thức trên ta có (2.21) ở đây m và n là số lẻ. Nếu là thực và có tính đối xứng vòng tròn, thì đáp ứng xung có thể được tính từ tích phân kép Dùng công thức trên có thể viết một chương trình C để thu được đáp ứng xung từ đáp ứng tần số. Dưới đây là chương trình được viết bằng Turbo C 2.0. Chương trình 2.1 "FILTERD.C". Tính đáp ứng xung từ đáp ứng tần số /*Program to compute the impulse from the frequency response.*/ /* Chuong trinh nay dung de tinh dap ung xung cua bo loc 2-D tuan hoan doi xung. Dap ung tan so cua bo loc duoc cho boi nguoi dung trong chuong trinh con tinh ham H(float w1, float w2). */ #include #include #include #include #include #define pi 3.1415926 int n1,n2; float f(float,float); float H(float,float); float Simpson2(float(*)(float,float),float,float,float,float,int,int); void main(void) { float f(float,float),xmin,xmax,ymin,ymax; int M,N; int NS,xt,yt; float h[11][11],temp; char file_name[40],ch; FILE *fptr; clrscr(); xmin=ymin=0.0; xmax=ymax=pi; M=N=20; /* So cac diem chia theo phuong x va y. cac gia tri nay phai la chan.*/ gotoxy(1,2); printf (" This program calculates the impulse response from the"); printf("\n frequency response."); printf("\n The impulse response is calculated over a window"); printf("\n centered around the origin."); printf ("\n Enter the number of samples in the impulse response."); printf("\n e.g. 7*7 (max. 11*11) (number have to be odd.)-->"); scanf("%d*%d",&NS,&NS); NS=(NS-1)>>1; /* chieu rong va chieu dai cua mot cua so*/ xt=wherex(); yt=wherey(); gotoxy(70,25); textattr(WHITE+(GREEN<<4)+BLINK); cputs("WAIT"); for(n1=0;n1<=NS;n1++) for(n2=0;n2<=NS;n2++) { if(n2>n1) continue; h[n1][n2]=simpson2(f,xmin,xmax,ymin,ymax,M,N); h[n1][n2]/=pi*pi; } for(n1=0;n1<=NS;n1++) for(n2=0;n2<=NS;n2++) { if(n2<n1) h[n2][n1]=h[n1][n2]; } gotoxy(70,25); textattr(WHITE+(BLACK<<4)); cputs(" "); gotoxy(xt,yt); printf("\n Impulse Response one quadrant.\n\n"); for(n1=0;n1<=NS;n1++) { for(n2=0;n2<=NS;n2++) printf("%7.6f ",h[n1][n2]); printf("\n"); } /* Dinh dang dap ung xung tuan hoan. */ temp=h[NS][NS]; for(n1=NS;n1<((NS<<1)+1);n1++) for(n2=NS;n2<((NS<<1)+1);n2++) h[n1][n2]=h[n1-NS][n2-NS]; h[(NS<<1)][(NS<<1)]=temp; for(n1=NS;n1<((NS<<1)+1);n1++) for(n2=0;n2<NS;n2++) h[n1][n2]=h[n1][(NS<<1)-n2]; for(n1=0;n1<NS;n1++) for(n2=0;n2<((NS<<1)+1);n2++) h[n1][n2]=h[(NS<<1)-n1][n2]; printf ("\n\nEnter Name of file for storing impulse response-->"); scanf("%s",file_name); fptr=fopen(file_name,"w"); for(n1=0;n1<((NS<<1)+1);n1++) for(n2=0;n2<((NS<<1)+1);n2++) fprintf(fptr,"%f",h[n1][n2]); fclose(fptr); printf("\nCalculating and storing frequency response "); printf("\nfor 3-D plotting.\n"); printf("----------------------------------------------------"); /* Chuan bi du lieu cho 3-D Du lieu duoc dinh dang thich hop voi chuong trinh in cua GRAFTOOL mot cong cu phan mem phan tich do hoa cua 3-D vision (412 S.Pacific Coast Hwy. Suite 201, Redondo Beach,CA. 90277). */ { /* Declaring local variable*/ float w1[33],w2[33],z[33],dw,sum; int row,col; register i,j; printf("\n Enter file name for storing 3-D plotting-->"); scanf("%s",file_name); fptr=fopen(file_name,"w"); row=col=34; fprintf(fptr,"%d %d\n",row,col); fprintf(fptr,"%e",(float)0.0); w1[0]=w2[0]=-pi; dw=pi/16.0; for(i=1;i<33;i++) w1[i]=w2[i]=w1[i-1]+dw; for(i=0;i<33;i++) fprintf(fptr," %e",w2[i]); xt=wherex();yt=wherey(); gotoxy(70,25); textattr(RED+(LIGHTGRAY<<4)+BLINK); cputs("WAIT"); gotoxy(xt,yt); for(i=0;i<33;i++) { fprintf(fptr,"\n"); for(j=0;j<33;j++) { sum=0.0; for(n1=-NS;n1<=NS;n1++) for(n2=-NS;n2<=NS;n2++) sum+=h[n1+NS][n2+NS]*(float)cos((double)(w1[i]*(float)n1+ w2[j]*(float)n2)); z[i]=sum; } fprintf(fptr,"%e",w1[i]); for(j=0;j<33;j++) fprintf(fptr,"%e",z[j]); } fclose(fptr); } xt=wherex();yt=wherey(); gotoxy(70,25); textattr(WHITE+(BLACK<<4)); cputs(" "); gotoxy(xt,yt); printf("\n Press any key to exit."); getch(); } /* Dinh nghia ham cho tich phan */ float f(float x,float y) { float H(float,float),a; a=H(x,y)*(float)cos((double)(x*n1))*cos((double)(y*n2)); return(a); } /*********************************************/ /*Chuong trinh con Simpson tinh tich phan kep*/ /*********************************************/ float simpson2( float(*f)(float,float),float xmin,float xmax, float ymin, float ymax, int M, int N) /* f la mot ham hai bien dinh nghia boi nguoi dung. xmin, xmax, va ymin, ymax la gioi han cua hai tich phan. M,N la so khoang cach tren huong x va y va chi co gia tri chan*/ { register i,j; float sum1,sum2,dx,dy,x,y,I; float *A; A=(float *) malloc(M*sizeof(float)); dx=(xmax-xmin)/(float) M; dy=(ymax-ymin)/(float) N; x=xmin; for (i=0;i<=M;i++) { sum1=sum2=0.0; y=ymin+dy; for(j=1;j<N; j++) { if((j%2)==0) sum1+=(*f)(x,y); else sum2+=(*f)(x,y); y+=dy; } *(A+i)=(*f)(x,ymin)+2.0*sum1+4.0*sum2+(*f)(x,ymax); x+=dx; } sum1=sum2=0.0; for(i=1;i<M;i++) { if((i%2)==0) sum1+=*(A+i); else sum2+=*(A+i); } I=*A+2.0*sum1+4.0*sum2+*(A+M); return(I*dx*dy/9.0); } /*************************************** * Dinh nghia ham bo loc H(w1,w2). Ham * * xac dinh boi nguoi dung * ***************************************/ float H(float w1, float w2) { float DO,a,R2; DO=0.5*pi; // cut-off R2=w1*w1+w2*w2; a=0.414*DO*DO/(R2+0.414*DO*DO); return(a); } /*************************************** Vi du ve mot so ham khac 1.Laplacian Operator. a=R2; 2. Phase contrast filters if(R2<(DO*DO)) a=1.0; else a=-1.0; 3. homomorphic filter float gamah, gamal, k; gamah=1.0; gamal=0.5; if (R2<(DO*DO)) { k=(1.0)-(gamah-gamal))/(gamah-gamal); a=(R2/(R2*DO*DO*k))+gamal; } else a=gamah; 4.Low pass filter a=0.414*DO*DO/(R2+0.414*DO*DO); 5.High pass filter a=R2/(R2+0.414*DO*DO); *****************************************/ Trong chương trình 2.1 đáp ứng xung được tính trên cửa sổ kích thước 11 ´ 11 mẫu mà tâm là gốc toạ độ. Một ảnh ba chiều của đáp ứng tần số tính từ 11 ´ 11 đáp ứng xung được chỉ trong hình 2.13. Bảng 2.1 Một phần tư đáp ứng xung của bộ lọc thông thấp D0 = 0.3p 0.103591 0.047501 0.017310 0.008810 0.003712 0.002135 0.047501 0.031265 0.014898 0.007494 0.003578 0.001873 0.017310 0.014898 0.009506 0.005306 0.002844 0.001496 0.008810 0.007494 0.005306 0.003352 0.001956 0.001102 0.003712 0.003578 0.002844 0.001956 0.001242 0.000743 0.002135 0.001873 0.001496 0.001102 0.000743 0.000473 Hình 2.13 Hình ảnh ba chiều đáp ứng biên độ của bộ lọc thông thấp. Ví dụ 2.6 Lặp lại bài toán trước với bộ lọc tuần hoàn đối xứng Butterworth được cho bởi trong đó Giải Thay hàm H(w1,w2) trong chương trình 2.1, như đã giải thích trong phần cuối của chương trình, chúng ta thu được hệ số được liệt kê trong bảng 2.2. Đáp ứng tần số được chỉ trong hình 2.14, và hầu như đúng với đáp ứng tần số gốc được tính từ biểu thức của . Hình 2.14 Đáp ứng tần số được tính từ 11 ´ 11 đáp ứng xung. Bảng 2.2 Một phần tư đáp ứng xung của bộ lọc thông cao 0.896408 -0.047501 -0.017311 -0.008811 -0.003712 -0.002135 -0.047501 -0.031265 -0.014898 -0.007494 -0.003578 -0.001873 -0.017311 -0.014898 -0.009506 -0.005306 -0.002844 -0.001496 -0.008811 -0.007494 -0.005306 -0.003352 -0.001956 -0.001102 -0.003712 -0.003578 -0.002844 -0.001956 -0.001242 -0.000743 -0.002135 -0.001873 -0.001496 -0.001102 -0.000743 -0.000473 Bài tập 2.2 Thay hàm H(w1,w2) trong chương trình 2.1 để tính đáp ứng xung trên cửa sổ 11 ´ 11 có tâm tại điểm gốc, theo hàm truyền đạt được cho dưới đây: 1. Lọc thông thấp: 2. Lọc thông cao: Với giả định đường tròn 3-dB nằm tại mức 0.3p.