I have not actually tried out your code, but based on my reading of it, it appears to have a serious flaw. For the sake of discussion suppose that your polygon is a triangle whose plane is already parallel with the x and y axes. That will give a zero for the two theta's. Then in the line
"points = points - repmat([0,0,points(1,3)],size(points,1),1);"
the three vertices are all projected vertically so they lie within the x,y plane itself. Finally in the for-loop Rx and Ry are identity matrices and the triangle is simply moved without rotation in a direction so that vertex #1 would move to the origin, except that only vertices #2 and #3 are actually moved. Vertex #1 has remained behind where it was, and the triangle has now assumed a completely different shape. I am sure this is not what you intended.
I will make two suggestions. First, there is no need to do two rotations. All you need is to rotate your polygon about a line which is orthogonal to both the polygon's normal and the z-axis, which is to say, about their cross product, and the amount of rotation should be such that the polygon's normal rotates to become parallel to the z-axis.
The second suggestion has to do with determining the normal direction of a polygon. Your method depends entirely on the first three vertices. However, if you have a polygon whose vertices deviate to some degree from being coplanar, it would be best to use matlab's 'svd' (singular value decomposition) function to determine the best normal direction. First subtract the mean values of the x, y, and z coordinates for the vertices from each vertex and then apply the 'svd':
[U,S,V] = svd(bsxfun(@minus,points,mean(points,1)),'econ');
Your "best" normal vector will be the unit vector V(:,3) which has the smallest singular value.
