Adding to what @Philipp said (especially the "caching blended result"). You can reverse calculate the second image from the
blend(grey,????) = colored equation for the multiply,overlay and screen functions. The problem is that the range of the
????? image may be out of the range of
0-255. So you might need to save them as a different format or encode something like a signed int into the 4 bytes of the RGBA values in one pixel for each channel.
As others said you can also use other color-models for help. HSV is the one that appears to be the obvious choice and it will look like this:
But as you can see the separated saturation and hue image still contains a good amount of the shading information, so a color model like YUV might seem more appropriate. It would look like this:
The problem here is that the inter-relationship that the RGB values in your blending source have with the final image might not always be so apparent.
Using just the raw components of these color models though has the issue of not being able to make image darker whatsoever. So I'd say use the reverse blending modes or combine the other color-spaces with the blending modes to get it done.
Edit: to visualize what I meant with the multiplication look at these examples:
The top 2 use multiplication with the regular color spectrum (0-255 mapped to 0-1), you can see that the image can only be made darker, making the red channel bright enough for a natural skin tone is basically impossible within the default spectrum. However if you look at the lower example of brightening it, then you can see that if you multiply it with values over 1 you can once again recreate the original result
Btw. if you're worried about the underwear being colored too, you need to mask it out (just use the grey base image where you don't want any coloration):
To calculate the scale value lets take for example the gray scale has the color [120 , 120, 120] at one point and the colored image has the color [152 107 74] at the same spot then the scale factors would be s_red = 152/120 = 1.27 , s_green = 107/120 = 0.892 and s_blue = 74/120 = 0.617, so [1.27 0.892 0.617].
You could reasonably map the re-color RGB value from [0-255] to [0-2] or something like that to easier use existing structures as well as limit the brightening to reasonable values. So your setup would look something like this (with the masking applied of course):
Don't forget to clamp your resulting image to [0-255] though.