Calculate actual Landsat image corner coordinates to derive azimuth (heading)

This is the continuation of an unsolved problem I had (see this thread if you like). I have a Landsat image and I would like to estimate the satellite azimuth for my location. As the satellite follows a near-polar orbit (see Landsat Handboook), the satellite must have an azimuth angle different than 0°, getting bigger as it moves poleward. If I have the coordinates of the corners of my image, I can estimate how the image is tilted against North direction (estimate the azimuth). So, the idea is:

1. 
   

   to choose one of the two sides of my image (wether the left or right) to get the angle of the normalized direction (see also How to calculate orientation of line segments using open source GIS?)

   
   


2. 
   
   

   to use the coordinates of the chosen corners (e.g. lower and upper left, as a and b variables, see below) in the following MATLAB function (inspired gist.github.com/604912)

   
   

   function d = line_dir( ax, ay, bx, by )
   
   dx= bx - ax;
   dy= by - ay;
   m= sqrt((dx^2)+(dy^2));
   r= atan2(dx/m, dy/m);
   d= (-r*180)/(pi);
   
   if d<0
   
       d=-d
   end
   end
   

However, the coordinates of the corners in the metadata file refers to the background image, not the real image (see the figure below to understand), so can I find now the real corners of the image, to then find the azimuth angle?

Answer

Eventually, I came up with a good result. This is the procedure I followed to estimate Landsat azimuth at my location.

1. I drew two segments in a GIS, one for each side of my scene (left and right, see figure 1), and added four ("real") corner points on the end of them (green points). This is done in the Reference System of the specific scene (in my case is WGS84-UTM32N).


2. I estimated the x and y coordinates of the four points with the field calculator (I use QGIS).


3. Finally, I use the coordinates in my function, using the left or the right corners, and choosing the lower corners as ax,y variable and the upper as bx,y.

My results are 13.3695 for the left side, and 13.7344 for the right part.

The proof the Landsat satellite is inclinated relies in Figure 2, showing a scene of the equator and my scene, and it underlines the distortion being higher the closer the satellites get to the poles (in this case, Northward).

My Landsat scene is from Italy (path/row= 193/028). I hope this can be of any help. Thanks to user30184 for the help!