The test image

image=ImageResize[ColorConvert[ExampleData[{"TestImage","Elaine"}],"Grayscale"],256];
f=ImageData[image];
Image[f]

Gaussian derivative based filters

Mathematica tip: You can make symbols by first pressing the "Esc" key and then type for example "Phi" and close off again by pressing "Esc". See also "Pallettes>Special Characters"
You can use the function With to locally define variables which make the final expression more readable.
The function GaussianFilter requires to specify a radius r of the filter, and a standard deviation σ of the Gaussian. We always set r = 3σ as for (Sqrt[(x^2)+(y^2)])>r the Gaussian has decayed to almost zero already.

(Subscript[Φ, 1]) and (Subscript[Φ, 2]) are first order derivatives along the x and y axis respectively (along the row-index and column index respectively)

(Subscript[Φ, 1])[f_,s_]:=With[{σ=(Sqrt[2s])},
GaussianFilter[f,{{3σ,3σ},{σ,σ}},{1,0}]
](*First order derivative in x-direction (row index direction)*)
(Subscript[Φ, 2])[f_,s_]:=With[{σ=(Sqrt[2s])},
GaussianFilter[f,{{3σ,3σ},{σ,σ}},{0,1}]
](*First order derivative in y-direction (column index direction)*)

(Subscript[Φ, 3]) is essentially the squre of the gradient magnitude

(Subscript[Φ, 3])[f_,s_]:=With[{σ=(Sqrt[2s])},
((GaussianFilter[f,{{3σ,3σ},{σ,σ}},{1,0}])^2)+((GaussianFilter[f,{{3σ,3σ},{σ,σ}},{0,1}])^2)](*Gradient magnitude squared*)

(∂_w) and (∂_v) are Gauge derivatives

partialw[f_,s_]:=With[{σ=(Sqrt[2s])},
With[{grad=Transpose[{GaussianFilter[f,{{3σ,3σ},{σ,σ}},{1,0}],GaussianFilter[f,{{3σ,3σ},{σ,σ}},{0,1}]},{3,1,2}]},
Map[With[{w=((#)/(Norm[#]))},#.w]&,grad,{2}](*inner product of the gradient with vector w*)
]]
partialv[f_,s_]:=With[{σ=(Sqrt[2s])},
With[{grad=Transpose[{GaussianFilter[f,{{3σ,3σ},{σ,σ}},{1,0}],GaussianFilter[f,{{3σ,3σ},{σ,σ}},{0,1}]},{3,1,2}]},
Map[With[{v=((RotationMatrix[π/2].#)/(Norm[#]))},#.v]&,grad,{2}](*inner product of the gradient with vector v*)
]]

((∂_v)^2) (second order (∂_v) derivative)

partialv2[f_,s_]:=With[{σ=(Sqrt[2s])},
With[{
(*The gradient encodes first order variations in the data*)
grad=Transpose[{GaussianFilter[f,{{3σ,3σ},{σ,σ}},{1,0}],GaussianFilter[f,{{3σ,3σ},{σ,σ}},{0,1}]},{3,1,2}],
(*The Hessian (here stored as a 2x2 matrix at each pixel) stores the second order variations*)
hes=Transpose[{{GaussianFilter[f,{{3σ,3σ},{σ,σ}},{2,0}],GaussianFilter[f,{{3σ,3σ},{σ,σ}},{1,1}]},{GaussianFilter[f,{{3σ,3σ},{σ,σ}},{1,1}],GaussianFilter[f,{{3σ,3σ},{σ,σ}},{0,2}]}},{3,4,1,2}]},
MapThread[With[{v=((RotationMatrix[π/2].#1)/(Norm[#1]))},v.(#2.v)]&,{grad,hes},2](*Computes the second order variations in the direction v, the Hessian is "projected" onto v*)
]]

(Subscript[Φ, 4]) is a corner detector

(Subscript[Φ, 4])[f_,s_]:=-partialv2[f,2]*((partialw[f,s])^2)

Apply the filters

Set a scale

s=1;

(Subscript[Φ, 1])

GraphicsRow[{Image[f],resultΦ1=ImageAdjust[Image[(Subscript[Φ, 1])[f,s]]]},ImageSize->Large]

(Subscript[Φ, 2])

GraphicsRow[{Image[f],resultΦ2=ImageAdjust[Image[(Subscript[Φ, 2])[f,s]]]},ImageSize->Large]

(Subscript[Φ, 3])

GraphicsRow[{Image[f],resultΦ3=ImageAdjust[Image[(Subscript[Φ, 3])[f,s]]]},ImageSize->Large]
(*It is the same as ((((∂_w)f))^2)*)
GraphicsRow[{Image[f],ImageAdjust[Image[partialw[f,s]^2]]},ImageSize->Large]

(Subscript[Φ, 4])

GraphicsRow[{Image[f],resultΦ4=ImageAdjust[Image[(Subscript[Φ, 4])[f,s]]]},ImageSize->Large]

Rotate the image and apply the filters again

Image rotation

f90=ImageData[ImageRotate[Image[f],π/2]];
Image[f90]

(Subscript[Φ, 1])

GraphicsRow[{Image[f90],resultΦ190=ImageAdjust[Image[(Subscript[Φ, 1])[f90,s]]]},ImageSize->Large]

(Subscript[Φ, 2])

GraphicsRow[{Image[f90],resultΦ290=ImageAdjust[Image[(Subscript[Φ, 2])[f90,s]]]},ImageSize->Large]

(Subscript[Φ, 3])

GraphicsRow[{Image[f90],resultΦ390=ImageAdjust[Image[(Subscript[Φ, 3])[f90,s]]]},ImageSize->Large]
(*It is the same as ((((∂_w)f))^2)*)
GraphicsRow[{Image[f90],ImageAdjust[Image[partialw[f90,s]^2]]},ImageSize->Large]

(Subscript[Φ, 4])

GraphicsRow[{Image[f90],resultΦ490=ImageAdjust[Image[(Subscript[Φ, 4])[f90,s]]]},ImageSize->Large]

Compare results by rotating the results of the 90 degree rotated inputs back

GraphicsRow[Image/@{f,f,f,f},ImageSize->1000,PlotLabel->"Original Input"]
GraphicsRow[{resultΦ1,resultΦ2,resultΦ3,resultΦ4},ImageSize->1000,PlotLabel->"Original Output"]
GraphicsRow[ImageRotate[#,-π/2]&/@{resultΦ190,resultΦ290,resultΦ390,resultΦ490},ImageSize->1000,PlotLabel->"(Subscript[R, -90])[Φ[(Subscript[R, 90])[f]]]"]