Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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).

share|improve this question
    
In case of confusion: I am in fact only interested in the last image (the one corresponding to the raw data). –  Andrei Jun 14 '13 at 15:48
    
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. –  Andrew Jaffe Jun 14 '13 at 15:49
    
@AndrewJaffe simply because I don't know how to do it in practice :)) –  Andrei Jun 14 '13 at 16:54
add comment

4 Answers

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["http://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

share|improve this answer
    
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). –  Andrei Jun 13 '13 at 19:18
    
@Andrei: Could you upload the raw images (without ticks, legend)? In the cropped image, the Fourier transform contains all kinds of artifacts. –  nikie Jun 13 '13 at 20:47
    
in fact, i would like to upload the raw data, but its a large matrix and i'm not sure how to do that –  Andrei Jun 13 '13 at 22:25
    
@Andrei, you could probably put the raw data on Pastebin. –  J. M. Jun 14 '13 at 0:30
    
@nikie, i uploaded a better image and the raw data corresponding to it is now on Pastebin. –  Andrei Jun 14 '13 at 14:30
show 2 more comments

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["http://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

share|improve this answer
add comment

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["http://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["http://i.stack.imgur.com/aIFO0.png"];
j2 = ImageAdjust@ImageSubtract[i2, GaussianFilter[i2, 5]]

enter image description here

share|improve this answer
    
nice! is there any way to do the reverse, i.e. a kind of low-pass filter that gets rid of these small waves? –  Andrei Jun 14 '13 at 15:44
1  
Just use the Gaussian-filtered image instead of subtracting it from the original? –  Rahul Narain Jun 14 '13 at 16:08
add comment

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["http://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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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