Q:

How to detect the intensity gradient direction

Having a Mat that is square area of grayscale pixels. How to create a straight line whose direction is created as a perpendicular to most pixel values change direction (average gradient, aerage over the whole Mat, the result would be just one direction (which can be then drawn as a line))?
For example having

it would look like

How can one do such thing in OpenCV (in python or C++)?

A:

An OpenCV implementation would look something like the following. It solves the problem in a similar fashion as explained in the answer by Mark Setchell, except that normalising the image does not have any effect on the resulting direction.
Mat img = imread("img.png", IMREAD_GRAYSCALE);

// compute the image derivatives for both the x and y direction
Mat dx, dy;
Sobel(img, dx, CV_32F, 1, 0);
Sobel(img, dy, CV_32F, 0, 1);

Scalar average_dx = mean(dx);
Scalar average_dy = mean(dy);

double average_gradient = atan2(-average_dy[0], average_dx[0]);
cout << "average_gradient = " << average_gradient << endl;

And to display the resulting direction
Point center = Point(img.cols/2, img.rows/2);
Point direction = Point(cos(average_gradient) * 100, -sin(average_gradient) * 100);

Mat img_rgb = imread("img.png"); // read the image in colour
line(img_rgb, center, center + direction, Scalar(0,0,255));
imshow("image", img_rgb);
waitKey();

A:

I can't easily tell you how to do it with OpenCV, but I can tell you the method and demonstrate using ImageMagick just at the command-line.
First, I think you need to convert the image to grayscale and normalise it to the full range of black to white - like this:
convert gradient.png -colorspace gray -normalize stage1.png

Then you need to calculate the X-gradient and the Y-gradient of the image using a Sobel filter and then take the inverse tan of the Y-gradient over the X-gradient:
convert stage1.png -define convolve:scale='50%!' -bias 50% \
  \( -clone 0 -morphology Convolve Sobel:0 \) \
  \( -clone 0 -morphology Convolve Sobel:90 \) \
  -fx '0.5+atan2(v-0.5,0.5-u)/pi/2' result.jpg

Then the mean value of the pixels in result.jpg is the direction of your line.
You can see the coefficients used in the convolution for X- and Y-gradient like this:
convert xc: -define morphology:showkernel=1 -morphology Convolve Sobel:0 null:

Kernel "Sobel" of size 3x3+1+1 with values from -2 to 2
Forming a output range from -4 to 4 (Zero-Summing)
 0:         1         0        -1
 1:         2         0        -2
 2:         1         0        -1

convert xc: -define morphology:showkernel=1 -morphology Convolve Sobel:90 null:
Kernel "Sobel@90" of size 3x3+1+1 with values from -2 to 2
Forming a output range from -4 to 4 (Zero-Summing)
 0:         1         2         1
 1:         0         0         0
 2:        -1        -2        -1

See Wikipedia here - specifically this line:

A:

Convert the image to grayscale and classify its pixels based on the gray level. For classification, you can use something like Otsu method or kmeans with 2 clusters. Then take the morphological gradient to detect the doundary.
Here are the classified pixels and the boundary using Otsu method.

Now find the non-zero pixels of the boundary image and fit a 2D line to those pixels using the fitLine function that finds a weighted least squares line or use this RANSAC implementation. fitLine gives a normalized vector collinear to the line. Using this vector, you can find an orthogonal vector to it.
I get [0.983035, -0.183421] for the collinear vector using the code below. So, [0.183421 0.983035] is orthogonal to this vector.
Here, in the left image, the red line is the least squares line and the blue line is a perpendicular line to the red one. In the right image, red line is the least squares line and the green one is the line fitted using the RANSAC library mentioned above.

Mat im = imread("LP24W.png", 0);

Mat bw, gr;

threshold(im, bw, 0, 255, CV_THRESH_BINARY|CV_THRESH_OTSU);

Mat kernel = getStructuringElement(MORPH_ELLIPSE, Size(3, 3));
morphologyEx(bw, gr, CV_MOP_GRADIENT, kernel);

vector<vector<Point>> contours;
findContours(gr, contours, CV_RETR_LIST, CV_CHAIN_APPROX_SIMPLE);
vector<Point> points;
for (vector<Point>& cont: contours)
{
    points.insert(points.end(), cont.begin(), cont.end());
}
Vec4f line;
fitLine(points, line, CV_DIST_L2, 0, 0.01, 0.01);
cout << line << endl;

