You can consider the cross product of the vectors and take the arcsin of the magnitude of that cross product.  Specifically, for the case of 2D vectors $\vec{a}$ and $\vec{b}$, we have

$$\sin{\theta} = \frac{a_x\, b_y - a_y\,b_x}{\sqrt{a_x^2+a_y^2}\sqrt{b_x^2+b_y^2}}$$

How to get angle bewteen two vectors in range -1 to 1 without using arc cosine?
I need to visualize 3D data having 200*200 data points. 
ListDensityPlot is a good candidate, but it seems to have strange performance issues.
I created a small test that plots 40000 points. It has three cases: using 2D array as input, using array of 3D points with integer coordinates $(x,y)$ and using array of 3D points with float coordinates $(x,y)$.

    TestFunction[x_, y_] := Sin[Pi/20* Sqrt[x^2 + y^2]];
    testdata = Table[TestFunction[i, j], {i, -100, 100}, {j, -100, 100}];
    testdataPoints = 
      Flatten[Table[{i, j, TestFunction[i, j]}, {i, -100, 100}, {j, -100, 
         100}], 1];
    testdataPoints2 = 
      Flatten[Table[{i*0.01, j*0.01, TestFunction[i, j]}, {i, -100, 
         100}, {j, -100, 100}], 1];
Benchmark[d_, n_] := 
  Timing[ListDensityPlot[d, ColorFunction -&gt; "Rainbow", 
    PlotRange -&gt; Full, InterpolationOrder -&gt; n]];
TableForm[
 Table[Benchmark[data, 
   n], {data, {testdata, testdataPoints, testdataPoints2}}, {n, 0, 
   2}], TableHeadings -&gt; {{"Array", "Integer", "Float"}, {0, 1, 2}}]




It gives pretty strange result. When I simply scale axes and my coordinates are not integer anymore (the case is called "Float"). The performance drops almost 10 times (3 seconds vs 25 seconds).

Here's the output timing:

            0         1         2
Array   0.834276 2.96802    3.05685  
Integer 4.84968  3.42562    3.19835
Float   27.5574   26.1669   25.9262


Any explanation for such behavior?

Edit

As alternative to Michael's solution one can use ArrayPlot (as Silvia suggested) if interpolation is not important.
It can be easily scaled to look like "Float" case.
