วันอาทิตย์ที่ 5 กรกฎาคม พ.ศ. 2552

How to computing camera matrix P

1. ทำการแรนดอมดาต้าของจุดใน 3D มาก่อน โดยทำการแรนดอมโดยฟิกค่าเฉลี่ยของระยะบนแกน Z อยู่ที่ 5 แล้วทำการแรนดอมค่าให้ได้เมตริกซ์ขนาด 4x30 โดยแต่ละคอลัมน์ของเมตริกซ์แสดงถึงค่าของ [ x ; y ; z+5; 1] ซึ่งผลลัพธ์ของขั้นตอนแรกนี้จะได้เมตริกซ์ลักษณะนี้ Xw= [x1,x2,....,x30 ; y1,y2,....,y30; z1+5,z2+5,...,z30+5;1,1,......,1]4x30 (Xw เป็นชื่อเมตริกซ์นะ)

2.จากนั้นทำการกำหนดค่า intrinsic camera parameter ลงในเมตริกซ์ขนาด 3x3 ดังนี้ K = [ 320 0 320 ; 0 320 240 ; 0 0 1 ]3x3 ในกรณีนี้ principal point ถูกเซ็ตค่าอยู่ที่ px=320(K(1,3)) py=240(K(2,3))

3.หลังจากกำหนดค่า intrinsic camera parameter แล้วก็ต้องกำหนด extrinsic camera parameter ด้วย กรณีนี้ให้ Rotation matrix มีค่าดังนี้ [1 0 0 ; 0 1 0 ; 0 0 1]3x3 และยังไม่มีการ translate จึงได้เมตริกซ์ t = [ 0 ; 0 ; 0 ] เมื่อได้ทั้ง extrin และ intrin แล้วจะได้เมตริกซ์ของกล้อง P ดังนี้ P = K*[ R t ] ซึ่งได้ผลลัพธ์ดังนี้ P= [320 0 320 0; 0 320 240 0; 0 0 1 0]3x4 เมตริกซ์ P นี้จะเป็นเมตริกซ์ที่ใช้ในการฉายภาพจากวัตถุใน 3D ลงมายังฉาก 2D (กล่าวคือ เปลี่ยนพิกัด 3D เป็น 2D บนภาพนั่นเอง)

4. เมื่อได้ P แล้วก็นำจุดที่แรนดอม Xw มาโปรเจคลงบนฉากโดยเราจะเรียกว่าเป็นการโปรเจคแบบอุดมคติ (ideal projection) โดยแต่ละจุดบนฉากจะแทนด้วยเมตริกซ์ x โดยที่ x = P*Xw ซึ่งมีขนาด 3x30 ซึ่งแต่ละคอลัมน์ของเมตริกซ์นี้จะเป็นจุดที่มีค่าทั้งแกน x , y และ z จากหลักการของ homogeneous matrix การแปลงจุดจาก 3D (x,y,z) สามารถแปลงไปเป็นจุดบนฉากได้โดยการนำ z ไปหาร x และ y ซึ่งจะได้พิกัดบนฉาก 2D เราจึงได้จุด (x/z,y/z) มาถึงตอนนี้เราจะได้พิกัดบนฉาก 2D ละ โดยเราเก็บจุดทั้ง 30 จุดไว้ในเมตริกซ์ x เหมือนเดิม เพื่อง่ายต่อการเข้าใจเราใช้คำสั่งนี้เพื่อให้ได้เมตริกซ์ของจุด 30 จุดบนฉาก2D x=x./repmat(x(3,:), 3 , 1 ) (ผลลัพธ์จะได้พิกัดจุดบนฉาก 30 จุด ซึ่งแถวที่สามในทุกคอลัมน์ของเมตริกซ์จะเป็น 1 หมด)

5.จากนั้นเราจะทำการเพิ่ม gaussian noise และทำเมตริกซ์ให้เป็น discrete มากขึ้นโดยการแรนดอมค่าระหว่าง 0 ถึง 1 แล้วบวกเข้ากับแถว 1 และ 2 ในเมตริกซ์ x จากนั้นก็ round (ปัดทศนิยมทิ้ง) แล้วเก็บไว้ในเมตริกซ์ x อันเดิม

6.หาค่าเฉลี่ยระยะห่างของจุดทั้ง 30 จุดใน 3D (ใช้ค่าในเมตริกซ์ Xw) กับตำแหน่งของกล้อง(ในที่นี้ตำแหน่งกล้องอยู่ที่ 0,0,0) โดยใช้ avg_dist_to_cam = mean( sqrt( sum(( Xw(1:3,:) - repmat( C, 1, n )).^2))); คำสั่งนี้จะทำการหาค่าระยะห่างระหว่างตำแหน่งกล้องและจุด 3D 30 จุด โดยในแต่ละคอลัมน์จะคำนวณดังนี้ ระยะห่างจากกล้องถึงจุด3D = sqrt((x-0)^2 + (y-0)^2 + (z-0)^2) (สาเหตุที่ x, y ,z ลบด้วย 0 เพราะว่า camera อยู่ที่ 0,0,0) ณ ขั้นตอนนี้เมตริกซ์ 3x30 เราจะถูกยุบเหลือ 1x30 แล้วเนื่องจากคำสั่ง sum จะทำการบวกแถวในแต่ละคอลัมน์นั่นเอง เมื่อ sum และ sqrt แล้วสิ่งที่เราได้ตอนนี้คือค่าระยะห่างจากกล้องถึงจุดทั้ง 30 จุดซึ่งเป็นเมตริกซ์ 1x30 นั่นเอง จากนั้นเราใช้คำสั่ง mean เพื่อหาค่าเฉลี่ยระยะห่างจากกล้องถึงจุดทั้ง 30 จุด แล้วเก็บไว้ในตัวแปร avg_dist_to_cam

7.ทำการนอร์มอลไลซ์ homogeneous coordinate
Xw = Xw ./ repmat( Xw(4,:), 4, 1 ); <--- จุดใน 3D
x = x ./ repmat( x(3,:), 3, 1 ); <-- จุดบนฉาก 2D

8.หา isotropic scaling เพื่อ normalize x กับ Xw (isotropic scaling คือ sx = sy = sz)
สำหรับ 3D จะหาเมตริกซ์ transform เพื่อจะนำจุดcentroid(ของจุดทั้ง 30 จุด)กลับมายัง origin โดยจะคำนวณ centroid ของจุดใน 3D โดยใช้ค่าเฉลี่ยของแต่ละแกน (x1+x2+...+x30)/30 แกน y ก็หาโดย (y1+y2+...+y30)/30 แกนz ก็หาโดย (z1+z2+...+z30)/30 ก็จะได้ centroid ของ จุด 30 จุด แต่เราจะต้องใส่เครื่องหมายลบให้ค่าของ centroid เนื่องจาก เป็นการนำ centroid กลับมายัง origin โดยใช้คำสั่งใน octave ดังนี้ Ttrans = [ eye( k ), -mean( Xw(1:k,:)' )' ; zeros( 1, k ), 1 ]; Ttrans คือเมตริกซ์ขนาด 4x4 ที่มีลักษณะดังนี้ [ 1, 0 , 0 , ค่าที่ใช้ย้ายในแกนx ; 0 , 1 , 0 , ค่าที่ใช้ย้ายในแกนy ; 0 , 0 , 1 , ค่าที่ใช้ย้ายในแกนz ; 0 , 0 , 0 ,1 ]4x4 (หากใครเคยเรียนCG มาคงคุ้น ๆ กับเมตริกซ์นี้ ซึ่งมันก็คือเมตริกซ์ที่เราใช้ในการ translate นั่นเอง)
มาถึงตอนนี้แล้วเราได้เมตริกซ์ที่ translate centroid มายัง 0,0,0 แล้ว ก็ใช้ให้มันเป็นประโยชน์ โดยการนำเมตริกซ์นี้มาใช้หาตำแหน่งของจุดทั้ง 30 จุดโดยอ้างอิงจากจุด 0,0,0 (แต่เดิมจุด30จุดมันไม่ได้อ้างอิงจากจุด 0,0,0 การจะ rotate , scale ก็ทำยากมากมาย) traslated_Xw = Ttrans*Xw ; (คูณเมตริกซ์เหมือน CG เลย เอา translate เมตริกซ์ไว้ข้างหน้าจุดไว้ด้านหลัง ไม่ใช่ Xw*Ttrans นะ)
เมื่อได้ translated_Xw แล้ว เราต้องหา scaling factor ที่เหมาะสมให้กับจุดทั้ง30จุดใน translated_Xw โดยหาระยะห่างจากจุดทั้ง 30 จุดกับ origin แล้วหาค่าเฉลี่ยรวมของระยะห่างทั้ง 30 จุด ดังนี้ s=sqrt(k)/mean(sqrt(sum(translatedXw(1:k,:).^2))) โดย k = size(Xw,1)-1 เมื่อได้ s แล้วเราก็จะสามารถสร้าง Tscale เมตริกซ์ซึ่งใช้สำหรับการสเกลได้ (ซึ่งเมตริกซ์นี้จะทำการ scale ทั้งแกน x,y,z เท่ากันเพราะใช้ s ตัวเดียวกัน หรือเรียกอีกอย่างว่าเป็น isotropic scale ไงล่ะ) ใช้คำสั่ง octave ได้ดังนี้ Tscale = [ s * eye(k), zeros(k,1) ; zeros(1,k), 1 ]; ซึ่ง Tscale จะมีลักษณะดังนี้ [ s , 0 , 0 ,0 ; 0 , s ,0 ,0 ; 0 , 0 , s , 0 ; 0 , 0 , 0 , 1 ]4x4 (เหมือนเช่นเคยมันคือเมตริกซ์ scaling ใน CG อีกนั่นเอง)

มาถึงตอนนี้เราได้เมตริกซ์ Trans(late) และ Tscale แล้ว จัดการรวมร่างมันซะ ให้ได้ Transform matrix ทำได้ดังนี้ T_for_Xw = Tscale*Trans;

ขั้นตอนที่กล่าวมานี้อธิบายสำหรับเมตริกซ์ Xw (จุดใน 3D) ซึ่งยังไม่จบ เราต้องหา Transform matrix สำหรับ 2D ด้วย ซึ่งก็ทำเหมือนกันเพียงแต่ค่า k จะเปลี่ยนเป็น k=size(x,1)-1
Ttrans = [ eye( k ), -mean( x(1:k,:)' )' ; zeros( 1, k ), 1 ];
translated_x = Ttrans*x;
s=sqrt(k)/mean(sqrt(sum(translated_x(1:k,:).^2)));
Tscale = [ s * eye(k), zeros(k,1) ; zeros(1,k), 1 ];
T_for_x=Tscale*Ttrans;

9. เมื่อได้ transform matrix มาแล้วเราก็ทำการ transform ทั้งจุดใน 3D(Xw) และก็จุดในฉาก(x) กันเลย โดยใช้คำสั่งดังนี้
xtilde = T_for_x*x;
Xwtilde = T_for_Xw*Xw;

10. เมื่อได้จุดที่ถูก transform แล้วเราจะมาประมาณค่าแบบเชิงเส้นของเมตริกซ์ P (เมตริกซ์ของกล้อง) โดยใช้ DLT algorithm เพื่อประมาณค่า P (DLT algorithm ไว้จะมาอธิบายวันหลัง) แต่มีโค้ดดังนี้
%------------------------------------------------------------------------

function P = dlt( xtilde, Xwtilde )

n = size( Xwtilde, 2 );
if n < 6
error( 'DLT for P requires at least 6 points' );
end;

A = [];

for i = 1:n

xi = xtilde( 1, i );
yi = xtilde( 2, i );
wi = xtilde( 3, i );
Xwi = Xwtilde( :, i );

Ai = [ 0, 0, 0, 0, -wi * Xwi', yi * Xwi' ;
wi * Xwi', 0, 0, 0, 0, -xi * Xwi' ];

A = [ A ; Ai ];
end;

[U,D,V] = svd( A );

% The solution P minimizing |Ap| subject to |p|=1 is the last column of V

P = reshape( V(:,12), 4, 3 )';
P = P / norm(P(3,1:3));

%------------------------------------------------------------------------

ไม่มีความคิดเห็น:

แสดงความคิดเห็น