This was never a complicated question, even in 2015. This is concise and works for logical and numeric masks. The background can be I/RGB, and so can the foreground color tuple.
inpict = imread('peppers.png');
mask = imread('redpepmask.png');
fgc = bsxfun(@times,mask,permute(fgcolor,[1 3 2]));
bgc = bsxfun(@times,(1-mask),im2double(inpict));
outpict = im2uint8(bsxfun(@plus,fgc,bgc));
If you really absolutely wanted code that works on nothing but logical masks, you can do that too! Believe it or not, this isn't significantly faster than doing all that multiplication.
inpict = imread('peppers.png');
mask = imread('redpepmask.png')>128;
fgchans = numel(fgcolor);
bgchans = size(inpict,3);
outpict = zeros(size(inpict));
for c = 1:max(fgchans,bgchans)
thischan = im2double(inpict(:,:,min(c,bgchans)));
thischan(mask) = fgcolor(min(c,fgchans));
outpict(:,:,c) = thischan;
outpict = im2uint8(outpict);
While FEX imoverlay() existed, that's basically a simple scalar alpha composition, so it wouldn't easily work for a masked task. The contemporary IPT imoverlay() would probably satisfy OP's needs, but it didn't exist then.
Circa 2015, MIMT did have tools to do this, though I don't remember if everything had been consolidated into replacepixels() yet. If nothing else, it should have at least supported logical masks. That said, even the current version of replacepixels() works at least back to R2009b. So if you wanted the "simplest" (easiest to use, fewest lines needed, fewest restrictions to accomodate), then MIMT replacepixels() would take the cake in 2015, same as it does today. inpict = imread('peppers.png');
mask = imread('redpepmask.png');
outpict = replacepixels(fgcolor,inpict,mask);
IPT imoverlay is limited in that it really can only do logical operations (no smooth masks, no transparency), and can only do so with (at most) a color tuple as the foreground. So even though IPT imoverlay() exists today, I can't help but see it as a tiny fraction of an answer to general masked composition.