Thanks again. It all makes more sense now. I’m just playing with making a small texture paint utility, no research project, I’m afraid 
Before I even saw your solution, I made a quick hack with a vertex program that handled the projection:
// TOOL ray info
//
// The incoming coordinates are expressed in the tools own coordinate system.
// They will have to be changed to world coords first, thus we begin with
// constructing the tool basis axes:
// I = JxK
// J = up
// K = - direction (shoots paint along negative Z)
vec3 tool_I = - normalize (cross (tool_up, -tool_direction));
vec3 tool_J = normalize (cross (tool_I , -tool_direction));
// And we also find the direction-vector of
// the RAY: tool_position ==> tool_position + x*tool_I + y*tool_J + z*tool_K
//
vec3 ray_dir = gl_Vertex.x * tool_I + gl_Vertex.y * tool_J + tool_direction;
// PLANE of the face
//
// Next we want to build info for the PLANE of the FACE that we are painting at,
// i.e. the normal of the plane: AxB
vec3 plane_OA = plane_A - plane_O;
vec3 plane_OB = plane_B - plane_O;
vec3 plane_N = normalize (cross (plane_OB , plane_OA));
// RAY-PLANE INTERSECTION
//
// Then we find the intersection point in global coordinates of the RAY and the PLANE
// Plane: (V-plane_O) dot plane_N = 0
// Ray : t*ray_dir + tool_position
//
// This yields: t = (plane_O - tool_position) dot plane_N
// -------------------------------------
// tool_dir dot plane_N
float RdotN = dot (ray_dir, plane_N);
if (RdotN < 0.001f)
{ // Problem!
gl_Position = vec4 (0.0f, 0.0f, 1.0f, 1.0f);
no_hit = 100.0f;
}
else
{
no_hit = 0.0f;
float t = dot((plane_O - tool_position) , plane_N) / RdotN;
vec3 hit = tool_position + t * ray_dir;
// Construct a basis I,J,K based on O,OA,OB
// P is the point we want to project into O,OA,OB coordinates
vec3 P = hit - plane_O;
vec3 I = normalize (plane_OA);
vec3 J = - normalize(cross (I,plane_N));
// Doing some math gives the stuff below as a way to find a matrix
// M = [m o]
// [n p] that satisfies:
// I = [a,0]
// J = [b,c]
// M*OA = [1,0] this operation maps the vector OA onto (1,0) or E1 in the O,OA,OB coordinate space
// M*OB = [0,1] this operation maps the vector OB onto (0,1) or E2 in the O,OA,OB coordinate space
//
// Then we project P onto I,J to get (i,j)
// Then we do a M*(i,j) to get the coordinates in O,OA,OB jargon: (s,t)
float a = length(plane_OA); // [a 0] is OA projected onto I,J
float b = dot(plane_OB, I); // [b c] is OB projected onto I,K
float c = dot(plane_OB, J);
float m = 1.0f/a;
float n = 0.0f;
float o = - m*b / c;
float p = 1.0f/c;
mat2 M = mat2(m,n,o,p);
float PonI = dot (P,I);
float PonJ = dot (P,J);
vec2 PinAB= M * vec2(PonI,PonJ);
// Going from P in O,OA,OB coordinates to the UV texture space is easy:
// PinUV = To + Toa*s + Tob*t
vec2 PinUV= tex_O + PinAB.x*(tex_A-tex_O) + PinAB.y*(tex_B-tex_O);
gl_Position = vec4(PinUV.x, PinUV.y, 1.0f, 1.0f) * 2 - 1;
}
gl_FrontColor = gl_Color;
gl_BackColor = gl_Color;
gl_TexCoord[0] = gl_MultiTexCoord0;
Though it required some hard debugging, it is running now, but the quality is not good as I suspected.
