A few of the stuff that appears in this response were initially packed together in a compact comment. Inevitably, in this answer, I would have to provide code, which I wouldn't like in this case. You have to get what is going on. So, I am going to try and talk a little bit more "around" the code and hopefully by the end of this, the code will have written itself.
In just a few words: $H$ is your frequency response in the frequency domain, it is based on a bandpass pattern but it is supposed to be applied in a directional way. So, build a bandpass and rotate it.
Now then...Bandpass "pattern". You need something that goes up and then down and when it goes up it stays up as flat as possible. Take a plot of some $y_1 = \frac{x}{a}$. Well, obviously, that's a line with a slope of $\frac{1}{a}$ and it crosses zero at zero. Nothing fancy going on there.
Now, offset that zero a bit to the right, so $y_2 = \frac{b-x}{a}$. What did you just do? You brought zero bang on on $b$ but also you introduced some NEGATIVE values for $x>b$.
To see why this will come handy, do a $y_3=(\frac{b-x}{a})^2$. This is a parabola. It goes down and then up. It is also curved, let's remove this curve by undoing the square with a root, so $y_4=\sqrt{(\frac{b-x}{a})^2}$. What's going on now? We have two lines, one that arrives down to $b$ and one that departs upwards from $b$. Let's now reverse these trends.
$y_5=\frac{1}{\sqrt{(\frac{b-x}{a})^2}}$.
What is any $u=\frac{1}{v}$? A hyperbola. Remember those two lines that meet at $b$? If you were to substitute $v$ for $y_5$, then, when $y_5$ goes down, towards $b$, it tends to smaller values, therefore, $u$ goes upwards!. And when $y_5$ departs from $b$ towards increasing values then $y_5$ goes downwards.
What happens bang on on $b$? That's when $y_5$ hits zero and if $v$ hits zero then $u$ goes to infinity. We don't like that. So, let's make that $v$ a sum, to control the denominator and ensure it never hits zero. $y_6 = \frac{1}{1+\sqrt{(\frac{b-x}{a})^2}}$. So now, when that root term goes to zero, $y_6$ stays at $1$.
Let's deal with that stay flat business. Remember that $\sqrt{q^2}$ trick we did to get rid of possible negative values of $q$? Well, we could have done this with $|q|$, which means taking the absolute value of some $q$. But this power raising has an added benefit. So, in general, if you raise to some power, you increase the slope of a function. That is, how sharply the curve evolves. Take some $q^2$ and $q^4$, for any two $q, q+\delta_q$ where $\delta_q > 0$, then, $q^4$ grows faster than $q^2$.
BUT!, raising to a power has another property as well, raising large numbers to a power makes them larger and raising small numbers ($ < 1$) to a power makes them...even smaller!. Why is this important for us? Because as our $y_2$ approaches its $b$ limit, it grows smaller and smaller and smaller and we can accelerate this decrease by raising to some power. Therefore, we now go $\sqrt{q^n}$ and use $n$ to control how fast that function will approach small values and stay small leading to that flatness characteristic.
Now, here is the thing with multiplication, if you raise some $q$ to an odd power, you retain the original sign of $q$ and if that was negative then $\sqrt(q)$ will become complex. We don't like that in this case. So, what do we do? We ensure that some $q$ is raised to an even power by multiplying by two. All together now:
$$y_7 = \frac{1}{1+\sqrt{(\frac{b-x}{a})^{2n}}}$$
OK, fine, how do I apply $y_7$ as a filter to a signal? Pretend that your $x$ is now frequency, mirror $y_7$ around some $x_{\frac{N}{2}}$ (where $N$ is the full length of your impulse response), take its inverse Fourier Transform, throw away the complex part (which should have been zero anyway) and treat what you get back as your impulse response. That is your $h$.
Hold on, if $x$ now represents frequency then what are the $b,a,n$? Consider how $y_7$ goes up towards something and then comes down from that something. So, $b$, is your Center Frequency. The centre frequency of the bandpass filter. $a$ is now your bandwidth, that is how fat that lobe of values that go up and down is. And finally, $n$ is now your filter's degree, that is, how flat it tends to be. In fact, if you raise $n$ to some absurd value (which, remember, it will be doubled), then it will look very much like a square pulse which would be an ideal filter. But, the more your crank up $n$ the more Fourier coefficients you need to synthesize it accurately, and then the impulse response gets longer (that $N$ from earlier) and...anyway, just keep that in mind, it's not important for now.
WHAT THE?!?!? HOW DOES THIS APPLY TO MY ORIGINAL QUESTION?!?!!?!?!?!
We are very close now, all that we have to do is bring in another dimension, so that our $y_7$ now gets a $y$ as well. (I hope using $y$ and $y_7$ is not confusing, they are just symbols)
I am not going to go all over this again because we need to focus on what is really important here. When you bring in another dimension you can now define direction. I am not going to go through the full two dimensional Fourier spectrum explanation but I will just say this. When you look at the frequency spectrum of an image, the low frequencies are at the centre and the high frequencies are at the periphery (by convention).
Therefore, filters that treat all directions equally look like cakes or donuts. It's not hard to imagine why. Take our bandpass profile $y_7$ and mentally rotate it about the end where its zero lies. This will give you a ring (or doughnut) about the centre of the image. Pick any direction, what you see is the bandpass filter and its symmetric image, remember earlier that I said, "...mirror it around some $x_{\frac{N}{2}}$..."?.
But, hey, we now have the luxury to filter along a specific direction and not treat all of them equally. So, let's take our filter profile and apply it to just a specific direction.
How are we going to do this?
Remember our $y_7$? This is effectively an $f(x)$, a function of some $x$ variable which we can assign any meaning we like. It could be some generic $x$, it could be a frequency axis segment. But, at the end of the day, what we did there was to scan $x$ and assign a $y$ value, so the height of $y$ depends on the value of $x$.
What we are going to do now, is scan the frequency plane and assign a height to it, so imagine your $H$ as being a height map. But! We want to be able to control for the direction which translates to rotation.
So, while I am scanning some $u,v$ point on the frequency plane, I am assigning its height on some $u^*, v^*$ which are the rotated coordinates, around the centre $t_x, t_y$.
At this point, I have to say, I can see why the authors want to call $D_v, D_h$ their cutoff frequencies and $tx,ty$ their "sub-band centre frequency" but from the viewpoint of the above explanation, $D_v, D_h$ would correspond to the bandwidth and $tx,ty$ would correspond to the cutoff frequency.
Now, you work in C# and you probably have very good reasons for this but my quick test was performed in Octave (which is very similar to MATLAB). It is a piece of FOSS and that's where you could have just done a imagesc(H) and looking at your results we could more or less understand where the problem is.
I am not going to give you the full code (yet) because that's not going to do any good. But I will write here the line where I suspect you are running into trouble (straightforward translation to C#):
H(i+N2+1,j+N2+1) = 1./(1+(0.414*(sqrt(((su./Dh)+(sv./Dv)).^(2.*n)))));
I will also suggest an alternative way that you can achieve the same thing, much more straightforwardly, just re-using the stuff I have been talking about so far. Less lines of code, less bugs to worry about.
Now, you might notice that the authors in their paper, state that $H$ should be made symmetric around the origin along the orientation you are trying to filter. So, at the end of the code that is building $H$ you have an extra step (which I suspect you are not doing at the moment) to mirror it. If you try to mirror this function AFTER it has been captured in some $H$ matrix, it is going to be hard work.
For the following to work, I am assuming that in C# you are using some elementary image processing library.
So, write a function that builds you a $y_7$ from above. Nice, clean, simple, one dimensional bandpass filter in the frequency domain. Once you have that, MIRROR IT by something like $y_8 = [y_7(N_{y_7}:-1:1) y_7]$. The indexing convention is $StartValue:Step:StopValue$, so effectively, if your $y_7$ was "Y7VALUES" your $y_8$ is now "SEULAV7Y7VALUES", notice that we skipped the value of 0. Otherwise, the "Y" would appear twice and that would mess the spectrum.
So, your $y_8$ is now one row of your $H$. To turn it into an "image mask", replicate $N$ copies of $y_8$, one below the other. Once you are done, you have your $H$ which if you were to look at as if it was an image, it looks like two parallel lines, two fat bands, one somewhere to the left of the middle of the image and one to the right. This is your $H$ at zero degrees.
To generate the same $H$ at a different orientation. Simply rotate that initial $H$.
DON'T ROTATE A ROTATED VERSION!!!! (so if you did 45 degrees, don't rotate that another 45 to go to 90), because that brings in distortions due to interpolation. Your image processing library should have elementary rotation capabilities. If it doesn't....good luck, you have the Hough transform coming up in the same paper...
I hope this helps.
imagesc(H)and post it? $\endgroup$