It should be not too difficult, following the next steps.

1. Find the 11 regions (connected components) R_1,...,R_11
and fix a point P_k in each of them
(each region is determined by 2, 3 or 4 inequalities of the form lines[i]<0 or lines[i]>0)
2.  OK:={}
3. For each S subset {P_1,...,P_11} having 4 elements do
    If S is admissible (i.e. in each halfspace there are exactly 2 points of S), then OK:=OK union {S}.

The set OK will contain all the desired configurations.

Here is a DirectSearch approch for finding a configuration. It is possible to find several configurations if we repeat the call.

lines := [-(1/2)*x+y, (1/2)*x+y, 1-x+(1/2)*y, 3-x-(1/2)*y]:
pl:=(i,j) -> eval(lines[i],[x=x[j],y=y[j]]):
cond1:=seq(seq(pl(i,j)^2>eps,i=1..4),j=1..4):
cond:= seq( op( [
      pl(i,1)*pl(i,2)*pl(i,3)*pl(i,4)>eps,
      min(pl(i,1),pl(i,2),pl(i,3),pl(i,4))<-eps,
      max(pl(i,1),pl(i,2),pl(i,3),pl(i,4))>eps  ])  ,
      i=1..4):
bounds:=seq(x[i]=-10..10,i=1..4), seq(y[i]=-10..10,i=1..4),  eps=0..10:
vars:=[seq(x[i],i=1..4),seq(y[i],i=1..4),eps]:
with(DirectSearch):
ds := Search(-eps, [cond, cond1,bounds], variables=vars):
pts:= [seq( [ds[2](i),ds[2](i+4)],i=1..4 )];

plot1:=plots[pointplot](pts,color=red,symbol=solidcircle,symbolsize=8):
plot2:=plots[implicitplot](lines,x=-10..10,y=-10..10):
plots[display](plot1,plot2,axes=none);

You will have to use IntegrationTools (or op()) for this.

J:=int( int( a(x)*b(y),x=-infinity..infinity), y=-infinity..infinity):

JJ:=expand(J);

with(IntegrationTools):

J1:=op(1,JJ); J2:=op(2,JJ):
int(int(GetIntegrand(J1)*GetIntegrand(J2),GetVariable(J1)=GetRange(J1)),GetVariable(J2)=GetRange(J2));

You should show how to reproduce this.

Otherwise, see my worksheet:

12x.mw

I get correct results.

Why should  int( fD(x)*fA(x), x=a..b ) be between 0 and 1? If you integrate the product of two densities, the result could be even infinity!

This is the first problem of the "Hundred-dollar, Hundred-digit Challenge problems":

https://en.wikipedia.org/wiki/Hundred-dollar,_Hundred-digit_Challenge_problems

The Maple solution was obtained by Robert Israel:

https://www.math.ubc.ca/~israel/challenge/challenge1.html

Edit: Mathematica seems to have better algorithms for computing oscillatory integrals.

Too many errors ...

neville:=proc(X::list,Y::list,x) # polynomial interpolation
local T,i,j,n;
n:=nops(Y)-1;
T:=Array(0..n, Y);
for j to n do
  for i from 0 to n-j do
  T[i] :=  normal( (T[i]*(X[j+i+1]-x) + T[i+1]*(x-X[i+1]) )/(X[j+i+1]-X[i+1]) )
od; od;
end:

Test:

f:=x->sin(3*x);
r:=neville([$1..10],[f(n)$n=1..10],x):
plot([f(x),r],x=1..10);

You were told that f is not a pdf unless d=0.

For d=0  f is a Negative binomial distribution which is implemented in Maple, so you can use it.

But you want d>0.

In this case, you probably want to mimic the moments defined for a pdf i.e.

M[k] = Sum( n^k * f(n), n=0..infinity);

But it can be shown that f(n) is equivalent (as n --> infinity) to K/n^2 where K is a nonzero constant.
Hence all the moments are infinite.

Actually, your with your procedure a sequence with both elements is produced. In Maple 2015, Eqplot(2,3) returns 2x+3 and the plot in a reduced size. 

If you use r := Eqplot(2,3), then r[1] equals 2*x+3 and r[2] equals a PLOT structure such that r[2] will display the plot when evaluated.

It would be better to use:

Eqplot:=proc(a,b)
local y,x;
y:=a*x+b;
print(plot(y,x));
y;
end proc:

This way, Eqplot(2,3) will return 2*x+3 ant the plot itself is produced as a side effect.

f:=2*x-2:
an:=2*int( f*cos(Pi*x*n),x=0..1) assuming n::integer;

a0:=2*int(f,x=0..1);

To compare graphically f with a partial Fourier sum:

s3:=a0/2 + add( eval(an,n=k) * cos(k*Pi*x), k=1..3):

plot([f, s3], x=0..1, view=[-0.2..1.2,-2.2..0.5]);

Install the package in a directory where you have write permission. After the .hdb conversion you can move it into Program Files if you want.

For d=0, f is a PDF  (actually a Negative binomial distribution).

It seems to me (after some investigations) that this is the only case when f is a PDF.
According to Maple, for [r=3, b=2, d=0.01411053479]  f seemed to be a distribution
just because  Sum_x f(x)   is very close to 1 for a small d.

Why don't you use 1D math? Later, you may convert to 2D if you want. All these problems will disappear.

I think that 1D for input and 2D for output is the best choice for now. 

se:=proc(f::procedure,a::realcons,b::realcons)
local x,s;
s:=(f(b)-f(a))/(b-a)*(x-a)+f(a);
print(plot([f(x),s], x=a..b));
s;
end proc:

f:=x -> x^2:

se(f,-1,2);

rad:=u -> `*`(op(map( z->op(1,z),factors(u)[2])));

g:=expand((x*sqrt(2)+y*sqrt(3)+z*sqrt(6)+x*y*sqrt(6)+x)^6*(x+2*y)^2):
rad(g);

If _C1 is constant just use _C1 instead of _C1(t).

If _C1(t) etc are already in an expression, use

subs( [_C1(t)=_C1, _C2(t)=_C2, _C3(t)=C3], expr);