Shock filtering using structure tensor
I think that the question here is "how is the dominant direction of the local slope actually estimated?" with the outlook of using it in a shock filter.

A Shock Filter is applied iteratively and each time it propagates grayscale values along the direction of an edge but not across the edge, therefore, progressively, the edge is preserved but the area towards the broad direction of the structure is "filled". This is similar to applying dilation or erosion but taking into account the local contrast gradient. The edge can be defined at various scale levels, which results to the average direction of an edge defined over a bigger area surrounding a pixel of interest $(x,y)$.

The code seems alright, in that it sets up the normalised Gaussian g and then derives its first order derivative in gd. It uses that to obtain the gradients Ix,Iy but sx, sxy should be pointwise, rather than matrix, operations. Otherwise, Ix, Iy, "filter" each other. By controlling for the variance parameter of g (and by extension, gd too), we control how big the area around the pixel of interest $(x,y)$ is (or, in terms of dilation and erosion, how big the disk -or structuring element- is).

At this point, Ix, Iy could be used to estimate the local contrast gradient but if we come across areas with roughly opposite slopes, then averaging their components would result to 0 and the (false) perception that the area is flat. Instead, the Structure Tensor is used and more specifically the eigenvector that corresponds to the largest eigenvalue (and there is one of those for each $(x,y)$ pair). For more information please see this link and specifically sections 2, 3 and 4.

So, with that vector indicating the direction of the gradient of an edge, we find the perpendicular vector (with that sign flipping) and propagate values in that direction.

Hope this helps.