%{
This code is to render a Mesh given a 3x4 camera matrix with an image resolution widthxheight. The rendering result is an ID map for facets, edges and vertices. This can usually used for occlusion testing in texture mapping a model from an image, such as the texture mapping in the following two papers.
--Jianxiong Xiao http://mit.edu/jxiao/
Citation:
[1] J. Xiao, T. Fang, P. Zhao, M. Lhuillier, and L. Quan
Image-based Street-side City Modeling
ACM Transaction on Graphics (TOG), Volume 28, Number 5
Proceedings of ACM SIGGRAPH Asia 2009
[2] J. Xiao, T. Fang, P. Tan, P. Zhao, E. Ofek, and L. Quan
Image-based Facade Modeling
ACM Transaction on Graphics (TOG), Volume 27, Number 5
Proceedings of ACM SIGGRAPH Asia 2008
%}
clear
clc
syms P00 P01 P02 P03 P10 P11 P12 P13 P20 P21 P22 P23
syms m_far m_near m_width m_height scale
projection = [P00 P01 P02 P03; P10 P11 P12 P13; P20 P21 P22 P23; 0 0 0 1];
% handle the near and far clip plane in OpenGL
protr = [1.0, 0.0, 0.0, 0.0;
0.0, 1.0, 0.0, 0.0;
0.0, 0.0, m_far / (m_far - m_near), - m_near * m_far /(m_far - m_near);
0.0, 0.0, 1.0, 0.0];
% handle half pixel inconsistency
offset = [1.0, 0.0, 0.5, 0.0;
0.0, 1.0, 0.5, 0.0;
0.0, 0.0, 1.0, 0.0;
0.0, 0.0, 0.0, 1.0];
% undo image aspect ratio and size
m0 = [m_width / 2, 0.0, 0.0, 0 + m_width / 2.0;
0.0, m_height / 2, 0.0, 0 + m_height / 2.0;
0.0, 0.0, 0.5, 0.5;
0.0, 0.0, 0.0, 1.0];
% handle scaling
m1=[1 / scale, 0.0, 0.0, 0.0;
0.0, 1 / scale, 0.0, 0.0;
0.0, 0.0, 1.0, 0.0;
0.0, 0.0, 0.0, 1.0];
% handle upside down in vertical direction in image
m2=[1, 0.0, 0.0, 0.0;
0.0, -1.0, m_height, 0.0;
0.0, 0.0, 1.0, 0.0;
0.0, 0.0, 0.0, 1.0];
M = (m0 \ protr * m2 * m1 * offset * projection);
M =
[ P20*(1/(m_width*scale) - 1) + (2*P00)/(m_width*scale), P21*(1/(m_width*scale) - 1) + (2*P01)/(m_width*scale), P22*(1/(m_width*scale) - 1) + (2*P02)/(m_width*scale), P23*(1/(m_width*scale) - 1) + (2*P03)/(m_width*scale)]
[ - P20*(1/(m_height*scale) - 1) - (2*P10)/(m_height*scale), - P21*(1/(m_height*scale) - 1) - (2*P11)/(m_height*scale), - P22*(1/(m_height*scale) - 1) - (2*P12)/(m_height*scale), - P23*(1/(m_height*scale) - 1) - (2*P13)/(m_height*scale)]
[ (P20*(m_far + m_near))/(m_far - m_near), (P21*(m_far + m_near))/(m_far - m_near), (P22*(m_far + m_near))/(m_far - m_near), (P23*(m_far + m_near))/(m_far - m_near) - (2*m_far*m_near)/(m_far - m_near)]
[ P20, P21, P22, P23]
inv_width_scale = 1/(m_width*scale);
inv_height_scale = 1/(m_height*scale);
inv_width_scale_1 =inv_width_scale - 1;
inv_height_scale_1_s = -(inv_height_scale - 1);
inv_width_scale_2 = inv_width_scale*2;
inv_height_scale_2_s = -inv_height_scale*2;
m_far_a_m_near = m_far + m_near;
m_far_s_m_near = m_far - m_near;
m_far_d_m_near = m_far_a_m_near/m_far_s_m_near;
M =
[ P20*(inv_width_scale - 1) + (2*P00)*inv_width_scale, P21*(inv_width_scale - 1) + (2*P01)*inv_width_scale , P22*(inv_width_scale - 1) + (2*P02)*inv_width_scale , P23*(inv_width_scale - 1) + (2*P03)*inv_width_scale]
[ - P20*(inv_height_scale - 1) - (2*P10)*inv_height_scale, - P21*(inv_height_scale - 1) - (2*P11)*inv_height_scale, - P22*(inv_height_scale - 1) - (2*P12)*inv_height_scale, - P23*(inv_height_scale - 1) - (2*P13)*inv_height_scale]
[ (P20*m_far_a_m_near)/m_far_s_m_near, (P21*m_far_a_m_near)/m_far_s_m_near, (P22*m_far_a_m_near)/m_far_s_m_near, (P23*m_far_a_m_near)/m_far_s_m_near - (2*m_far*m_near)/m_far_s_m_near]
[ P20, P21, P22, P23]
% final OpenGL matrix = M'
% column 1
P20*inv_width_scale_1 + P00*inv_width_scale_2,
P21*inv_width_scale_1 + P01*inv_width_scale_2,
P22*inv_width_scale_1 + P02*inv_width_scale_2,
P23*inv_width_scale_1 + P03*inv_width_scale_2,
% column 2
P20*inv_height_scale_1_s + P10*inv_height_scale_2_s,
P21*inv_height_scale_1_s + P11*inv_height_scale_2_s,
P22*inv_height_scale_1_s + P12*inv_height_scale_2_s,
P23*inv_height_scale_1_s + P13*inv_height_scale_2_s,
% column 3
P20*m_far_d_m_near,
P21*m_far_d_m_near,
P22*m_far_d_m_near,
P23*m_far_d_m_near - (2*m_far*m_near)/m_far_s_m_near,
% column 4
P20
P21
P22
P23
% final matrix to opengl
P20*inv_width_scale_1 + P00*inv_width_scale_2, P20*inv_height_scale_1_s + P10*inv_height_scale_2_s, P20*m_far_d_m_near, P20,
P21*inv_width_scale_1 + P01*inv_width_scale_2, P21*inv_height_scale_1_s + P11*inv_height_scale_2_s, P21*m_far_d_m_near, P21,
P22*inv_width_scale_1 + P02*inv_width_scale_2, P22*inv_height_scale_1_s + P12*inv_height_scale_2_s, P22*m_far_d_m_near, P22,
P23*inv_width_scale_1 + P03*inv_width_scale_2, P23*inv_height_scale_1_s + P13*inv_height_scale_2_s, P23*m_far_d_m_near - (2*m_far*m_near)/m_far_s_m_near, P23
inv_width_scale = 1/(m_width*scale);
inv_height_scale = 1/(m_height*scale);
inv_width_scale_1 =inv_width_scale - 1;
inv_height_scale_1_s = -(inv_height_scale - 1);
inv_width_scale_2 = inv_width_scale*2;
inv_height_scale_2_s = -inv_height_scale*2;
m_far_a_m_near = m_far + m_near;
m_far_s_m_near = m_far - m_near;
m_far_d_m_near = m_far_a_m_near/m_far_s_m_near;