Krdytkyu

does anyone know how to "manually" implement the built-in Sharpen function?

I'll show you my problem below:

Say I have a certain dataset, for example the one generated below:

data = Table[

        (Exp[-2 (x - 1)^2] + Exp[-2 (x + 1)^2])*Exp[-(y)^2], 

       {x, -5, 5, 0.2}, {y, -5, 5, 0.2}];

data = data/Total[data, 2];

And that I apply some Blur-like noise in the following way (I'm aware of the built-in function Blur, but in this case I prefer to have full control on what is happening):

BlurMatrix[mat_List, hBlur_List, vBlur_List] := 

 Block[{dim, row, col, h1, h2, v1, v2},

  dim = Dimensions[mat];

  {col, row} = ConstantArray[0, #] & /@ dim;

  v1 = Total[

    Drop[Join[ArrayReshape[ConstantArray[col, #], {dim[[1]], #}], mat,

          2], None, -#]*vBlur[[#]] & /@ Range[Length[vBlur]]];

  v2 = Total[

    Drop[Join[mat, ArrayReshape[ConstantArray[col, #], {dim[[1]], #}],

          2], None, #]*vBlur[[#]] & /@ Range[Length[vBlur]]];

  h1 = Total[

    Drop[Join[ConstantArray[row, #], mat], -#]*hBlur[[#]] & /@ 

     Range[Length[hBlur]]];

  h2 = Total[

    Drop[Join[mat, ConstantArray[row, #]], #]*hBlur[[#]] & /@ 

     Range[Length[hBlur]]];

  Total[{mat, v1, v2, h1, h2}]/Total[{mat, v1, v2, h1, h2}, 3]

  ]

This function takes the initial dataset and add noise to each pixel by mixing it with the adjacent pixels with different (tunable) weights. For example:

noiseH = {0.9, 0.7, 0.5, 0.1};

noiseV = {0.9, 0.7, 0.5, 0.1};

blurMatrix = BlurMatrix[data, noiseH, noiseV];

blurMatrix is clearly a blurred version of the original data:

GraphicsRow[

 ListDensityPlot[#, ImageSize -> Small, PlotRange -> All, 

    InterpolationOrder -> 0] & /@ {data, blurMatrix}]

And the difference between the original data and the blurred one is:

difference = Abs[blurMatrix - data];

ListDensityPlot[difference, ImageSize -> Small, PlotRange -> All, 

 InterpolationOrder -> 0]

My question is, how do I implement a function that makes the blurred image as close as possible to the original data?

I have tried with the built-in function Sharpen, but I'd like to have a "manual" version of it, where I know what the algorithm is doing and I can modify it accordingly to each specific noise scenario, i.e. knowing what kind of noise is acting on the "true" data (in my case I know the form of noiseH and noiseV) how do I apply the function to get rid of the noise as much as possible.

Anyway, for the benefit of others, see below my implementation with the Sharpen function:

distance = Table[

   sharpen = 

    ImageData[Sharpen[Image[blurMatrix/Max[blurMatrix]], i]];

   sharpen = sharpen/Total[sharpen, 2];

   {i, Total[Abs[sharpen - data], 2]}

   , {i, .1, 4, .05}]; 



First[SortBy[distance, Last]][[{1, 2}]]

ListLinePlot[distance[[;; , {1, 2}]], Frame -> True, 

 FrameLabel -> {"Sharpen argument", "Distance"}]

ListDensityPlot[First[SortBy[distance, Last]][[3]], 

 ImageSize -> Small, PlotRange -> All, InterpolationOrder -> 0]

{2.45, 0.0310388}

In this case, the best argument for the sharpen function is 2.45, but I found it completely empirically not knowing what the sharpen function does (the documentation doesn't help much in this sense).

In my view, you have two unrelated questions here

does anyone know how to "manually" implement the built-in Sharpen function?

This is what Mathematica does under the hood

mySharpen[img_Image, r_] := 

 Image[3.*ImageData[img] - 2.*ImageData[GaussianFilter[img, {r}, Padding -> "Fixed"]]]

In this case, the best argument for the sharpen function is 2.45, but I found it completely empirically not knowing what the sharpen function does (the documentation doesn't help much in this sense).

Here, it gets a bit more involved. What you are looking for is an image deconvolution. What you are referring to when you say "the best argument" is the point spread function (PSF) which degenerated (or blurred) your image. The PSF can be thought of what happens when you convolve a single white pixel with it. This is called an "impulse response". Now, what you basically want is blind deconvolution which tries to improve the image without knowing the specific PSF.

Without going into further detail, I suggest you first read all documentation of ImageDeconvolve and look up the wiki pages I referenced. This is a large topic because image deconvolution is an ill-posed problem.

It looks like it's a convolution with a fixed filter. You can find the convolution kernel by applying Sharpen to an impulse signal (i.e. a single white pixel):

impulse = Image[SparseArray[{{8, 8} -> 1}, {16, 16}]];

impulseResponse = Sharpen[impulse];

The impulse response looks something like this:

MatrixPlot[ImageData[impulseResponse]]

Then Sharpen is just a convolution with that impulse response:

rnd = RandomImage[1];

ImageDistance[Sharpen[rnd], ImageConvolve[rnd, impulseResponse]]

0.

ADD: I just noticed that your question has a second part; You want to invert your custom blurring function.

I won't pretend that I fully understood your BlurMatrix function. I'll just assume that it performs a linear shift-invariant filter (i.e. a convolution) and will use the same technique as above to find the kernel:

impulse = SparseArray[{{8, 8} -> 1}, {15, 15}];    

impulseResponse = BlurMatrix[Normal[impulse], noiseH, noiseV];    

MatrixPlot[impulseResponse]

Meaning: Your BlurMatrix function should be equivalent to ListConvolve[impulseResponse, data, {8,8}].

Now we can use this kernel matrix to deconvolve the blurred matrix to reconstruct the original:

deconvolution = ListDeconvolve[impulseResponse, blurMatrix];



GraphicsRow[

 ListDensityPlot[#, ImageSize -> Small, PlotRange -> All, 

    InterpolationOrder -> 0] & /@ {data, blurMatrix, deconvolution}]

The plot on the left is the original data, the middle plot shows the blurred version and the plot on the right shows the result of the deconvolution, which is close to the original:

MatrixPlot[deconvolution - data, PlotLegends -> Automatic]