The answer is fairly simple, but it depends what you expect. The basic concept isn't a MATLAB problem at all. Solving the problem starts with knowing how to do it in concept -- which is something for which there are plenty of photoshop tutorials.
That said, asking how to do it in MATLAB is a fair question. It can certainly be done with base MATLAB and IPT tools, but it won't be convenient or apparent. I'm going to use MIMT tools here, since MATLAB/IPT don't have any blending or composition tools. The displacement mapping can likely be done with imwarp(). The specific usage of imlnc() isn't directly replicable with other tools, but imadjust() is likely sufficient. FG = imread('crumpledpaper.jpg');
BG = imread('peppers.png');
FG = imresize(FG,imsize(BG,2));
FG = imlnc(FG,'independent','g',2.2,'k',0.6);
BG = displace(BG,[1.5 1.5],'xmap',FG,'ymap',FG,'edgetype','replicate');
R = imblend(FG,BG,1,'contrast',0.5);
Why is the FG image shifted to be centered about 50% gray?
Blend modes like 'overlay', 'soft light', etc are members of a category of (often misleadingly-named) bidirectional modes. They can both lighten and darken the BG image depending on the content of the FG image. These modes have a neutral response (i.e. R = BG) when FG = 50%. Shifting the FG to have a centered distribution means it will roughly lighten and darken in equal amounts.
Why are there four ouputs images? Which is the right one?
The simple answer is that there are four images to show that there is no right one. Most PS tutorials will probably suggest to use 'overlay', because that's the hammer PS has at hand. A GIMP tutorial might suggest layering multiple 'overlay' layers (because GIMP 'overlay' is actually a permanent bug). Another GIMP tutorial might suggest to use 'grain merge' and then reduce the opacity.
MIMT imblend() leaves that decision up to you. You can use 'overlay' or 'grain merge' or any other bidirectional mode you want. Note that most of these modes are also adjustable, so workarounds like stacking duplicate layers or muting blown-out results by abusing opacity aren't necessary.
The last two modes are different. While the others are a simple function of FG, BG, and the scalar adjustment parameter, 'scaleadd' and 'contrast' depend on the mean color of the FG image. While the neutral response of 'grainmerge' is fixed at FG = 0.5, the neutral response of 'scaleadd' is at mean(FG). In these two cases, the whole task of centering the FG histogram is typically unnecessary.
The comparison here is apropos, since 'scaleadd' is essentially an adjustable version of 'grainmerge' with (default) mean-centering behavior. They both have roughly the same intended purpose.
The synopsis and this PDF constitute a fair description of the modes you have available.