Use a computer algebra program to compute the projections onto the isotypic components.
The subgroup $G\colon=T_d$ of $\OO(3)$ is gernerated by $C_2,C_3,S_4$ and $\s$ (cf. subsection). In sage the projections $P_j$ can be calculated as follows - actually we calculate $12P_j$:
import numpy as np
C3=matrix(3,[0,1,0,0,0,1,1,0,0])
C2=matrix(3,[1,0,0,0,-1,0,0,0,-1])
S4=matrix(3,[-1,0,0,0,0,1,0,-1,0])
S=matrix(3,[0,0,-1,0,1,0,-1,0,0])
P=matrix(3)
G=MatrixGroup(C2,C3,S4,S)
Chars=G.irreducible_characters()
for chi in Chars:
PP=matrix(np.kron(P,P))
for g in G:
mg=g.matrix()
mgg=matrix(np.kron(mg,mg))
PP=PP+chi(g)*mgg
print(PP)
$T_d\simeq S(4)$ and $S(4)$ as any symmetric group has two generators; however we took $C_2,C_3,S_4$ and $\s$, which are representatives of the conjugacy classes. In particular we have $P_3e_2=P_3e_3=P_3e_4=P_3e_6=P_3e_7=P_3e_8=0$ and
$$
P_3e_1=\frac{2e_1-e_5-e_9}3,\quad
P_3e_5=\frac{-e_1+2e_5-e_9}3,\quad
P_3e_9=\frac{-e_1-e_5+2e_9}3
$$
i.e. $P_3e_1$ and $P_3e_5$ form a basis for $\im P_3$ (here $e_1,\ldots,e_9$ denotes the standard basis of $\R^9$); it's easily seen that e.g.
$$
b_1\colon=\frac{2e_1-e_5-e_9}{\sqrt6}
\quad\mbox{and}\quad
b_2=\frac{e_5-e_9}{\sqrt2}
$$
form an orthonormal basis of $\im P_3$. In tensor notation the orthonormal basis is given by (here $e_1,e_2,e_9$ denotes the standard basis of $\R^3$):
$$
b_1=\frac{2e_1\otimes e_1-e_2\otimes e_2-e_3\otimes e_3}{\sqrt6},\quad
b_2=\frac{e_2\otimes e_2-e_3\otimes e_3}{\sqrt2}~.
$$
Finally the matrices of $C_3\otimes C_3$, $C_2\otimes C_2$, $S_4\otimes S_4$ and $\s\otimes\s$ with respect to the basis $b_1,b_2$:
$$
\left(\begin{array}{cc}
-1/2&\sqrt3/2\\
-\sqrt3/2&-1/2
\end{array}\right),
\left(\begin{array}{cc}
1&0\\
0&1
\end{array}\right),
\left(\begin{array}{cc}
1&0\\
0&-1
\end{array}\right),
\left(\begin{array}{cc}
-1/2&-\sqrt3/2\\
-\sqrt3/2&1/2
\end{array}\right)~.
$$
Of course this can also be done in sage: Find the image of 'PP' and determine an orthonormal basis of the image:
ImP=PP.image().matrix().transpose()
Q,R=np.linalg.qr(ImP)
print matrix(Q).str(rep_mapping=lambda x: str(round(x,4))
The columns of the matrix $Q$ are the components of an orthonormal basis with respect to $e_1\otimes e_1,\ldots,e_3\otimes e_3$ (precision: 4 digits). You may download the list of sage commands and run in sage:
load('181.sage')