17
$\begingroup$

Consider the contour plot below, of a function f(x,y) - which is calculated numerically by solving some equations. If one looks carefully, one can observe two types of waves, to the left and right of the depression at (0,0). To the right we see waves of long wavelength (only two crests are visible actually). The waves to the left, on the other hand, are much more closely spaced and much less intense, but one can still distinguish them. Let us ignore the vertical stripes, which are of no interest.

real-space image

Now if we take the 2D Fourier transform of the data, we get a picture like this:

fourier

The very small circle in the center corresponds to the waves of large separation, and the other 2 circles correspond to the washed out, closely-spaced waves.

My question is, what type of manipulation can I do (some filtering/convolution comes to mind) in order to give more weight to the washed-out waves to the left (make them more visible)?

Edit:

Just to clarify the above. I am not interested in any of the parallel vertical lines in the real-space plot, either to the right or to the left. What I am interested in are the very small parabolic waves very close to the center of the image, on the left side. I try to point them out in the zoom below, because they might be hard to spot in the original plot.

zoom

Edit 2:

I would prefer a way of selecting just the waves which correspond to a certain ring in the Fourier transformed image, something like a band pass filter which selects one ring. That is because I suspect that the close-spaced waves on the left of the center correspond to one particular ring only (the one pointed to by the blue arrow in the fourier plot).

The main motivation for doing things like this is that then I want to apply the same method to the problem below:

data = Import["http://pastebin.com/raw.php?i=ieiLy8AD"];
img = Image[Rescale[data], ImageSize -> 400]

space-mat

In this image, the washed out waves to the left of the dark central point are even weaker. My gut feeling is that they can be best brought out by filtering in Fourier space.

However, I am willing to accept any answer that can transform this last image to one where the small ripples are seen clearly and then make a cut of that 2D plot to calculate their wavelength (along the y=0 axis for instance).

$\endgroup$
3
  • $\begingroup$ In case of confusion: I am in fact only interested in the last image (the one corresponding to the raw data). $\endgroup$
    – Andrei
    Jun 14, 2013 at 15:48
  • $\begingroup$ So why not just do what you suggest yourself? Multiply your fourier domain picture by something which highlights the modes you are interested in, since you know where they are in fourier space. $\endgroup$ Jun 14, 2013 at 15:49
  • $\begingroup$ @AndrewJaffe simply because I don't know how to do it in practice :)) $\endgroup$
    – Andrei
    Jun 14, 2013 at 16:54

4 Answers 4

13
$\begingroup$

You could simply remove the vertical waves (e.g. by subtracting the median of each column) and histogram modification. Using @bill s's cropped image:

img = Image[
  ImageData[
    ColorConvert[Import["https://i.stack.imgur.com/EvjuW.png"], 
     "Grayscale"]][[;; , ;; , 1]]]; (* remove alpha channel *)
columnMedian = Median[ImageData[img]];
medianRemoved = # - columnMedian & /@ ImageData[img];
ImageAdjust[Image[medianRemoved]]

This effectively removes the vertical waves:

enter image description here

To enhance the contrast of the waves left and right of the center, I'd simply apply a function like this to each pixel:

contrastEnhancement = #/Sqrt[Abs[#] + .05] &;
Plot[contrastEnhancement[x], {x, -1, 1}]

enter image description here

ImageAdjust[Image[contrastEnhancement[medianRemoved]]]

enter image description here

Or, using ListPlot3D:

ListPlot3D[GaussianFilter[-contrastEnhancement[medianRemoved], 5], PlotRange -> All, ImageSize -> 600]

enter image description here

$\endgroup$
7
  • $\begingroup$ great work @nikie. one thing remains - isn't there some way of selecting just the waves which correspond to a certain ring in the Fourier transformed image? say something like a band pass filter which selects one ring? because i know that the close-spaced waves on the left of the center correspond to one particular ring only (namely the bright ring on the left of the 2D fft plot). $\endgroup$
    – Andrei
    Jun 13, 2013 at 19:18
  • $\begingroup$ @Andrei: Could you upload the raw images (without ticks, legend)? In the cropped image, the Fourier transform contains all kinds of artifacts. $\endgroup$ Jun 13, 2013 at 20:47
  • $\begingroup$ in fact, i would like to upload the raw data, but its a large matrix and i'm not sure how to do that $\endgroup$
    – Andrei
    Jun 13, 2013 at 22:25
  • $\begingroup$ @Andrei, you could probably put the raw data on Pastebin. $\endgroup$ Jun 14, 2013 at 0:30
  • $\begingroup$ @nikie, i uploaded a better image and the raw data corresponding to it is now on Pastebin. $\endgroup$
    – Andrei
    Jun 14, 2013 at 14:30
8
$\begingroup$

It actually works pretty well to just subtract off a Gaussian-blurred copy of the image. This is effectively a cheap high-pass filter, but one that does not introduce any ringing (i.e. extra waves that were not present in the original image).

i1 = Image[ImageData[ColorConvert[Import["https://i.stack.imgur.com/EvjuW.png"], 
      "Grayscale"]][[;; , ;; , 1]]]; (* image loading code borrowed from @nikie's answer *)
j1 = ImageAdjust@ImageSubtract[i1, GaussianFilter[i1, 5]]

enter image description here

i2 = Import["https://i.stack.imgur.com/aIFO0.png"];
j2 = ImageAdjust@ImageSubtract[i2, GaussianFilter[i2, 5]]

enter image description here

$\endgroup$
2
  • $\begingroup$ nice! is there any way to do the reverse, i.e. a kind of low-pass filter that gets rid of these small waves? $\endgroup$
    – Andrei
    Jun 14, 2013 at 15:44
  • 1
    $\begingroup$ Just use the Gaussian-filtered image instead of subtracting it from the original? $\endgroup$
    – user484
    Jun 14, 2013 at 16:08
7
$\begingroup$

Some deconvolve algorithm can not recover fine details very well, which is of course a shortcoming. However, by taking advantage of this kind of shortcoming, we can separate fine detailed structures from large scale structures effectively.

We borrow the same image loader from nikie's answer:

img = Image[
        ImageData[ColorConvert[
          Import["https://i.stack.imgur.com/EvjuW.png"], 
          "Grayscale"]][[;; , ;; , 1]]];

First we blur img with a kernel whose "width" is close to the characteristic scale of the fine structure we want, which will suppress it maximally:

imgblur   = ImageConvolve[img, GaussianMatrix[7]];

Then we choose a proper (in the sense of relatively incapable of obtaining sharp images) deconvolve algorithm, say the Wiener filter:

imgdeblur = ImageDeconvolve[imgblur, GaussianMatrix[7], Method -> "Wiener"];

At last we subtract imgdeblur from img and adjust the contrast, brightness and gamma correction. Hopfully we have filtered out most of the large scale structures:

ImageSubtract[img, imgdeblur] // ImageCrop[#, 120] & // ImageAdjust[#, {0, 7, .5}] &

result image

For the remained vertical bars, in my opinion they look like well-fitted for some sine wave, so it should not be hard to filter out by 1-D Fourier transformation along the horizontal axis.

Or we can simply take the average along vertical axis to separate the sine-like background:

sinebg = Image[ConstantArray[
                    Mean[ImageData[imgsubed]],
                    {ImageDimensions[imgsubed][[2]]}]]

ImageSubtract[imgsubed, sinebg] // ImageAdjust

result 2

$\endgroup$
4
$\begingroup$

Another approach would be to use a wavelet method: to multiply the wavelet coefficients at a specified level by some given factor, and then invert. In this implementation, start with the removal of the vertical features as suggested by nikie.

img = Image[ImageData[ColorConvert[Import["https://i.stack.imgur.com/EvjuW.png"], 
        "Grayscale"]][[;; , ;; , 1]]];
columnMedian = Median[ImageData[img]];
medianRemoved = # - columnMedian & /@ ImageData[img];
imgProc = ImageAdjust[Image[medianRemoved]];

The technique is to pick a given level/scale, multiply the wavelet coefficients at that level/scale by a constant, and then take the inverse wavelet transform. Because it's hard to predict what levels and scales might be the best for a given image, it's implemented as a manipulate so it's easy to scan through and pick the most informative images.

Manipulate[swt = StationaryWaveletTransform[imgProc];
   mapWave = WaveletMapIndexed[ImageMultiply[#, m] &, swt, level];
   ImageAdjust[InverseWaveletTransform[mapWave]],
   {{level, {0,0,0,2}, "level"}, Dynamic[swt[All][[All, 1]]]}, {{m, 20, ""}, 0.1, 40}]

enter image description here

Here are a couple of other snapshots with different settings:

enter image description hereenter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.