Fitting a curve between two sets of points

You can use FindMinimum to fit a smooth spline between these point sets.

First of all, let's extract the red and black points in your image:

img = Import["https://i.stack.imgur.com/D0KpG.png"];
pts = ComponentMeasurements[
      Binarize[Erosion[ColorNegate[ColorDistance[img, #]], 1], .5], 
      "Centroid", #Circularity > .99 &][[All, 2]] & /@ {Red, Black};

I'm going to normalize all X/Y coordinates to 0..1 for clarity:

pts[[All, All, 1]] = Rescale[pts[[All, All, 1]]];
pts[[All, All, 2]] = Rescale[pts[[All, All, 2]]];    
{red, black} = pts;
ListPlot[pts, AspectRatio -> 1, ImageSize -> 400]

enter image description here

Next step: Let's define the spline that will fit between these two point sets. A spline is a sum of spline basis functions. I'll use 16 spline basis functions on the intervall [0..1]:

nVars = 16;
vars = c /@ Range[nVars];

degree = 3;
knots = Join[ConstantArray[0, degree], Subdivide[nVars - degree], 
   ConstantArray[1, degree]];

These basis functions look like this:

Plot[Evaluate[
  Table[BSplineBasis[{degree, knots}, i, x], {i, 0, nVars - 1}]], {x, 
  0, 1}, PlotRange -> All, ImageSize -> 400]

enter image description here

And the spline we're looking for is a weighted sum of these basis functions:

fn[x_] := 
 Table[BSplineBasis[{degree, knots}, i, x], {i, 0, nVars - 1}].vars

so e.g. fn[.3] evaluates to

0.000166667 c3 + 0.221167 c4 + 0.657167 c[6] + 0.1215 c[7]

The smoothness of the spline is simply:

smoothness = Total[Differences[vars, 2]^2];

(actually, I'm not sure this is 100% correct. I think the spacing of the knots should make a difference here. But since the knots are spaced equally between 0..1, it probably doesn't matter.)

And we have the constraint that the line should be above all black and below all red points:

constraints = 
 Flatten[{fn[#[[1]]] <= #[[2]] & /@ red, 
   fn[#[[1]]] >= #[[2]] & /@ black}];

Which leads to the simple optimization:

solution = FindMinimum[{smoothness, constraints}, vars]

This takes a few seconds to evaluate, and if there is no spline with these knots that fits the constraints, FindMinimum will return an error. But for these points, and nVars=16, I get this nice looking solution:

Show[ListPlot[pts], 
 Plot[fn[x] /. solution[[2]], {x, 0, 1}, PlotStyle -> Red], 
 AspectRatio -> 1, ImageSize -> 400]

enter image description here

You can play with degree and nVars to get different curves. For example, for degree=0 and nVars=128, you get this piecewise-constant function:

enter image description here

and for degree=1, you get a piecewise linear function:

enter image description here


"I am already very interested in any solution assuming monotonicity of f(x)."

When the two lists are separable by a monotone curve (as is the case in OP) you can you can use Internal`List`Min twice to get the lower and upper envelopes, respectively, of the two lists.

Using the red and black from Niki's answer

{iFred, iFblack} = Interpolation[#, InterpolationOrder -> 1] & /@
     {{-1, 1} # & /@ SortBy[First][Internal`ListMin[{-1, 1} # & /@ red]],
      {1, -1} # & /@ SortBy[First][Internal`ListMin[{1, -1} # & /@ black]]}; 

Show[ListPlot[{red, black}, AspectRatio -> 1, 
    PlotStyle -> {Red, Black}, ImageSize -> Large], 
 Quiet @ Plot[{iFred[t], iFblack[t], (iFred[t] + iFblack[t])/2}, {t, 0, 1}, 
    AspectRatio -> 1, PlotStyle -> {Orange, Gray, Blue}, 
    Filling -> {1 -> {{2}, Opacity[.5, Green]}}]]

enter image description here