(This may take a while. I will try to keep it short, in which case I may have to go back to my answer to answer the questions.) RGB color space interpolation often uses trilinear interpolation, which can be built on top of a pair bilinear interpolations. But there is no need to use trilinear interpolation. In fact, other interpolates are often better, for example, a simplicial (or tetrahedral) interpolator is usually preferable for a number of reasons due to trilinearity. There are several such tetrahedral lattice openings that can be used. One of them is pretty standard. (I will not go into details, at least not yet.) In addition, there is no reason why one MUST interpolate to RGB instead of any other space, although it can be argued that RGB has its own special problems, usually around the interpolation of neutrals and close neutrals.
A characteristic related to RGB and interpolation is that the neutral is defined as a point such that R = G = B. The trilinear interpolator will have a maximum error along this neutral axis and usually has a characteristic (notched) shape for errors along the neutral path through the color space.
So how do we interpolate in 3rd? I assume that it is interpolated in the regular lattice of points in the color space. In this case, you can define a cube containing any one point. If you interpolate inside a scattered set of points, then the simplest solution is usually to construct a triangulation of these points, and then make a simplicial (linear) interpolation inside any given tetrahedron. In any case, higher order interpolators are problematic here, as in some cases they can cause color problems. For example, I would not want to see turns along the gradients. This can happen because ringing is a serious problem for spline interpolations in areas with relatively high curvature. And if gamma mapping is involved, then such transitions will undoubtedly be a problem. Even if there is no need to display gamma, there are still problems with gamma.
There are several ways to construct domain triangulation from scattered data. Alpha forms are based on Delaunay triangulation and are a smart choice. But assuming you have a regular lattice and you want to do trilinear interpolation, the problem boils down to interpolation inside a simple cube in the 3rd.
Note that trilinear interpolation is not truly linear interpolation than bilinear interpolation. These schemes are linear ONLY along the axes of the lattice, but along any other path through the color space they are polynomial in nature. Thus, the trilinear interpolator will display cubic polynomial behavior along the main diagonal or along most common paths through the cube. We can convince ourselves that trilinear interpolation is NOT truly linear, since we interpolate between 8 points. at 3, 4 points define truly linear interpolation as a function of these independent variables, but we have 8 points that define the cube. That is, we will consider the mapping from one RGB space to another as really 3 independent mappings, so RGB → UVW (I chose UVW here to represent some common other color space, which may or may not be RGB in the character. )
The trick is that we are building a trilinear interpolator, interpolating between two bilinear interpolations. We build these bilinear interpolations, interpolating linearly between two points along one edge, and then do a third interpolation between them. Thus, we can consider a trilinear interpolator consisting of 7 simple linear interpolations. It is interesting that it can be shown that it does not matter which axes we do primarily with interpolations. Thus, we can first interpolate along R, then B, then the G axis or choose any other order - the trilinear interpolator will be unique and identical for any selected order. (The same goes for a bilinear interpolator.)
So the trick is how do we do linear interpolation between two triads of points? First, we need to determine where on the line segment between these points we lie. For example, consider two points in our color space that lie along the red (R) edge of a cube. I will use the same values that you specified for these points, this way:
Q1 = [66, 51, 77] Q2 = [55, 66, 77]
These are the values that we will interpolate between, essentially, the result of our mapping, but we also need to know where these points lie in the RGB input space. Assume that these coordinates, based on the coordinates of the cube from which they came, are as follows:
P1 = [0, 0, 0] P2 = [1, 0, 0]
This is a single cube in the 3rd, as I wrote it, so the rest of the points lie in
P3 = [0, 1, 0] P4 = [1, 1, 0] P5 = [0, 0, 1] P6 = [1, 0, 1] P7 = [0, 1, 1] P8 = [1, 1, 1]
Of course, any common cube also works, and there is no reason for it to be a true cube. Any trihedral rectangular 4-sided prism will also work here. You can always convert things into a single cube.
Now suppose we want to interpolate along this edge of the cube between P1 and P2 into the region defined by Q1 and Q2? Select a point along this edge. You can see that only R changes along this edge between these points, therefore we only care about the value of R at the point at which we interpolate. Think of it as a percentage of the distance along the edge. Interpolation is simply the weighted average of two endpoints, a linear combination. Thus, for a point with a red value of r along the edge from 0 to 1 in the red channel, our interpolation will be
Q(r) = Q1*(1-r) + Q2*r
As you can see, when r is 1/2, so halfway along the edge, our interpolation will decrease to
Q(1/2,0,0) = (Q1 + Q2)/2
Logically, the average will be the average of the two endpoints. You interpolate for EVERY output channel independently.
Q(1/2,0,0) = ([66, 51, 77] + [55, 66, 77])/2 = [60.5, 58.5, 77]
Does this help restore endpoints? Of course. When r = 0 or r = 1, you can see that it returns exactly the corresponding Q1 or Q2.
Again, you are doing this interpolation along each of the four red edges for a trilinear interpolator. Then you do two more interpolations, perhaps along the green edges of the four results that we got above. Finally, you do another interpolation along the blue edge to get trilinear interpolation. Again, it doesn't matter in which order you choose the interpolation axis. The result will be mathematically the same.
You stopped at bilinear interpolation, then there are three such linear interpolations. Yes, it is true that bilinear interpolation or trilinear interpolation can also be performed as a weighted combination of all 4 (or 8) corners of a rectangle (or cube). This can be left in the future.