cancel
Showing results for 
Search instead for 
Did you mean: 
cancel
Showing results for 
Search instead for 
Did you mean: 

Community Tip - You can subscribe to a forum, label or individual post and receive email notifications when someone posts a new topic or reply. Learn more! X

resample function

Jbryant61
4-Participant

resample function

Hi.

Eden kindly gave me this fantastic function to resample a matrix at higher resolution. It is only now that I am using it and can't quite work out whats actually happening - although it does the job beautifully. Would any experts be kind enough to give me a quick explanation.

I can see it starts with taking the FFT and then checking for odd number of pixels in the image. But from there I am lost.

thanks
Jason
83 REPLIES 83

Before any explanation, check if it does something on real images. The function seems to spread the spot, like zooming the spot ... looks very complicated just for spreading the spot. I may be wrong immediately and wish so.

Jean

... I have added PIX so you can see what the function does on real images. Interesting but maybe not so much as it distorts the image quite badly, though it really wants to "oversample" ! Comparing with zoom the next would be to find a "spreading function" that would associate with zoom and produce less distortion ... RemToDo .

Jean
PhilipOakley
5-Regular Member
(To:ptc-1368288)

It takes the fourier transform which will create the frequency spectrum and then extends (spreads out) that to fill the bigger space. finally it inverts the spectrum back to the image.

It is a fairly standard approach.

Philip Oakley

Thanks.
Exactly how doe sit fill in (Im trying to understand the mathcad approach), Im guessing its done by the submatrix part - how does it know which frequencies to fill in?

Fascinating!
Jason
PhilipOakley
5-Regular Member
(To:Jbryant61)

On 8/6/2009 12:16:09 PM, Jbryant61 wrote:
>Thanks.
>Exactly how doe sit fill in
>(Im trying to understand the
>mathcad approach), Im guessing
>its done by the submatrix part
>- how does it know which
>frequencies to fill in?
>
>Fascinating!
>Jason

Basically it just doubles everything (that is the spectrum's indexes) and then scales the amplitudes to match the spreading out.

result - same central amplitude, picture twice as big, reasonably smooth infill in the gaps.

To have 'fun' try using say the 'lena' image and look for ripples near the edges, which is where the fft assumes the picture could be used like a repetitive wallpaper pattern, and smooths across that 'gap' as well. Big steps in the butted pictures produce bigger artefacts.

Philip Oakley

Is it possible to see the spectra of the raw image and then compare to the spectra of the high res image so I can see ij frequency space what it has done?

Thanks
RichardJ
19-Tanzanite
(To:Jbryant61)

On 8/6/2009 12:31:33 PM, Jbryant61 wrote:
>Is it possible to see the
>spectra of the raw image and
>then compare to the spectra of
>the high res image so I can
>see ij frequency space what it
>has done?

Easy enough.

It adds zeros to the array. It's called zero-padding or zero-filling. There's a good explanation here:

http://www.dspguru.com/howto/tech/zeropad.htm

Richard
Jbryant61
4-Participant
(To:RichardJ)

Very clear - thanks Richard.

On 8/6/2009 2:09:24 PM, Jbryant61 wrote:
>Very clear - thanks Richard.
_____________________________

Both functions [Eden, Richard] need either form of recentering. Proof not possible from your asymmetric A but quite obvious from my very near symmetric 'I'... RemToDo.
You can probably make 'I' totally symmetric to complete the proof.

Jean

On 8/6/2009 4:39:37 PM, jmG wrote:
>On 8/6/2009 2:09:24 PM, Jbryant61 wrote:
>>Very clear - thanks Richard.
>_____________________________
>
>Both functions [Eden, Richard] need
>either form of recentering. Proof not
>possible from your asymmetric A but
>quite obvious from my very near
>symmetric 'I'... RemToDo.
>You can probably make 'I' totally
>symmetric to complete the proof.
>
>Jean
>_________________________________

Proof confirmed: yellow 'II' symmetric & centered.
resample2D_freq(II,30,30) = incorrect.
Just to help complete "oversampling".

Jean



Jbryant61
4-Participant
(To:RichardJ)

Thanks for the link - its a great explanation and I have learned a new technique for interpolating - its really fascinating.

I have tried to work thru and example in 1D and I am almost there. The interpolated space domain curve doesn't quite line up with the raw data.

Thanks
Jason

On 8/7/2009 4:47:37 AM, Jbryant61 wrote:
>Thanks for the link - its a
>great explanation and I have
>learned a new technique for
>interpolating - its really
>fascinating.
>
>I have tried to work thru and
>example in 1D and I am almost
>there. The interpolated space
>domain curve doesn't quite
>line up with the raw data.
>
>Thanks
>Jason
______________________________

You have probably more chance winning the lottery than finding Fourier on the web. Learning Fourier is a personal step by step process. If you want to stay with the "Spot", better use my clean one rather than your crappy one. Some plots in your original don't plot ? By padding 0's, we are back to months ago with my 'U' [below Marlett], and Lou. And we are back to the pixelation by "zoom", which "zoom" we don't have and desperately need to conclude pixelation, though the Zoom does perfectly demonstrate the right process of pixelation, as exemplified last night. Pixelation goes by spreading (you call it interpolation). In lieu of Fourier and all the proposals, you can use the "Diagonal Spline interpolator" and CreateMesh for the spreading. It may not be as perfect but we have yet to wait some guru make resample2D work as it should ... RemToDo and there lot of them in this collab.

Read again the 3 Fourier conjectures about pixelation. Work with real images, apply your techniques and compare (my last work sheet yesterday). You have the perfect pixelation in hand, what is the problem not to use as passed ?


Jean
Jbryant61
4-Participant
(To:ptc-1368288)

Hi Jean.

Are you saying that using the creatmesh does the interplolation?

why does the resample2D not work?

thanks
Jason

On 8/7/2009 11:48:07 AM, Jbryant61 wrote:
>Hi Jean.
>
>Are you saying that using the
>creatmesh does the
>interplolation?
>
>why does the resample2D not
>work?
>
>thanks
>Jason
_____________________________

Oh ! sure it does interpolate/pixelate.
Is there a "Diagonal Fourier interpolator" ?
Never heard of it ... RemToDo.
Ordinary data and image data aren't the same.

Don't know why resample2D does not work as it would apparently ? You said from Eden, then from Lou ... does not matter. Once compacted it's more difficult to debug. My interest in all those things is limited to educative. Maybe the concept should be on CDFT [Cosine Discrete Fourier Transform].

Jean



...(suite),

... your last double cumulative integral is all correct ! It sums the pixel values in both directions under the supplementary constrains.
>"D" is the result from method 1 - fitting in the space domain to a functional form and then integrating it.
Is the discrepancy expected ?< Answer = theoretically YES.... really NO.
... because you are integrating spix4 from resample2D, which if not perfect is sufficiently good and as it seems, respects (�) the noisy data set you had originally. Based on that consideration, you can only expect an approximate radius, immediately from spix4 pixelated image. We don't have much points from spix4 that would describe a circle seating at some level in the spix4, but enough to detect that the coordinates of the center of spix4 is different than you had assumed. Here again we could fit the diagonal spline and CreateMesh for + points, but nothing will change really except for an estimated more accurate level of the � cumulative 2d integral. In the attached, what is done is an attempt to locate the closest radius in comparison to the theoretical one you had from pure fitting followed by the Given/Find the scalar integral. The good news is that the difference between the two cumulative integrals (the fit and the pixelated) just can't be materialised !. That's what you were looking for...is it ?

Jean

...errata:

What a gross mistake I made !
Corrections in blue, the two curves separate
with terminal values for you to appreciate.

Jean

BTW, Jason
Your Volume 2d cumulative integral is a lovely tool,
that I saved and added in my tool box, thanks.

Jean

On 8/8/2009 12:07:29 AM, jmG wrote:
>BTW, Jason
>Your Volume 2d cumulative
>integral is a lovely tool,
>that I saved and added in my
>tool box, thanks.
>
>Jean
____________________________

The attached is the 1d Fourier interpolator. This work results from a long lasting collaboration between Theodore & Jean [original work from Theodore]. This tool has done wonderful curve fitting and quite a few in this collab ! For unknown reasons, no collabs have commented ... lack of understanding ? from two old wolves ? from the wrong continent ? Whatever, you can appreciate this direct interpolator. I might be able to make it 2d and done with pixelation at will. Time is a low ingredient.

For unknown reasons your original converted 11 work sheet had to be made compatible, not so bad. As you can see, the original 19 rows spreads on so many rows as you put npts, which in 2d would mean spreading over a 76 x 76 image size ... (76=4*19 in this example).
It is a remarkable fact that most collabs consider their stuff first and don't adapt easily or not at all to others. You might have comments about the Fourier interpolator considering other avenues. See you sometimes.

Jean

Jao, back to your thread:

http://collab.mathsoft.com/read?126658,63

I have abandoned the 2d Fourier interpolator for two reasons:

1. Don't know how to do 2d
2. It does not work well at low pixel values.
It is in fact only the Fourier polynomial fit and we are really at the the other Fourier in pixelation. Read what I posted last night "Pixelation Oversampling", comparing the 3 methods. I want to study again "resampe2D" [Lou/Eden]. I had another question to myself: what about the DCT [Discrete Cosine Transform] where there is no sine. Even that one, I can't make it work ... DCT is fine but not iDCT. Any way, you had a good reply about the Volume integral, i.e: comparing the fit and the real.

A lot of numerical maths aren't popular on the web, if at all !

Jean
LouP
11-Garnet
(To:Jbryant61)

This looks like the 2D reasample I had in my "pixelation_sampling(5bQAQA2).mcd" sheet. The explanation there, repeated here, was obviously not very clear.

This is an interpolation using DFT's.Here's the explanation in 1-D. For a given sample spacing xs corrrsponding to a spatial frequency fs = 1/xs, a sampled system can represent only spatial frequencies from 0 to a spatial freq fNyquist = fs/2. The spectrum is periodic at fs, and the spectrum from fs/2 to fs is the conjugate mirror of the spectrum from 0 to fs/2.

Having a denser spatial sampling grid means a higher spatial sampling frequency. Using a numeric example, let fs = 10 correspond to 0.1 mm = 100 micron sample spacing. Anything represented by this sampling grid has unique frequenncy components from 0-5, with the spectrum from 5-10 being the conjugate mirror. If we want the same function with the same spectrum but sampled at 20 micron spacing, then we need fs=50.

Create a new spectrum as follows:
freq 0-5: use orig 0-5 spectrum
freq 5-45: add zero intermediate values
freq 45-50: use orig 5-10 spectrum

This new spectrum has the same freq components as the original, but the mirror is now withrespect to 50, not f=10. Take the inverse DFT, and you have the original bandlimited function, but now sampled at a higher finer spatial pitch (sampling freq is higher).

Now what I hope is a bit clearer:

Take a simple example of a, 5 point, 1D "image" (a vector) with a DFT = (2,4+2j,7-j,7+j,4-2j). Note the conjugate symmetry, assuming the original fct. is real. Assume these correspond to frequencies(spatial or temporal - take your pick) 0,1,2,3,4, with a sampling freq fs =5. The spectrum repeats for both positive and negative freq's. The complete spectrum is:

frq DFT value
...
-2 7+j
-1 4-2j
0 2
1 4+2j
2 7-j
3 7+j
4 4-2j
5 2
6 4+2j
...


Only the one period is need for the description. The resampling preserves the frequency components that are already present, but increases the sampling frequency to get more points covering the same span.The unique components here are at freq's 0,1,2. The components at freq=-1,-2 are the conjugates of the f=1,2 components. Since the DFT is periofdic at f=5, the f=-2,-1 components are the same as the f=+3,+4 components. Only freq's up to fn = fs/2, called the Nyquist freq., can be uniquely defined as result of the periodicity and conjugate symmetry constraints. In our example , this is up to 2.5, (=2 for the discrete DFT case).

Let's resample to fs=8, so that we will have 8 "pixels covering the same overall image span, rather than the original 5. We need to preserve the conjugae symmetry and periodicity, so if we need to keep the f=0,1,2 components, then for fs=8, we already know what the fs=8,7,6 components will be. fn is now 8/2 = 4, so we have additional components we can define. since we don't want to add new, nonzero freq componets, we make all the new ones zero.

The resampled period becomes:

frq DFT value
...
0 2
1 4+2j
2 7-j
3 0
4 0
5 0
6 7+j
7 4-2j
8 2 (repeat of f=0)
...

The inverse DFT of this will have more spatial sample points, but the identical frequency content, of the original.

If there is a component exactly at fs/2 (fct. of odd vs. even # points), then this component need to be split in half when inserting the new zero components, part going to the upper half and part to the lower half. Here are simple examples:


orig 5 pt spectrum(no pt at fn): 5,3,9,9,3
new 10 pt spectrum:5,3,9,0,0,0,0,0,9,3 (,5...)

orig 6 pt spectrum(pt at fn with val = 8): 5,3,9,8,9,3
new 10 pt spectrum: 5,3,9,4,0,0,0,4,9,3 (,5...)

The bulk of the routine is taking care of the odd/even logistics, and doing it for both rows and cols. Hope this is a bit clearer. Follow a simple example yourself.

Lou


Jbryant61
4-Participant
(To:LouP)

Lou - I apologise, it was in fact all of your work and not Eden's.

Thanks for the explanation.
Your worksheet is fantastic - still using it!

Not sure I will be able to improve Eden pixelation "resample2D(a,nr,nc)". It can't be that "simple" from literature, though a long program. I like my technique on two counts:
1. simplicity & modularity
2. superior quality [trustworthy]

Appreciate or discard. Work with real images.
The total zoom factor 19.75


Jean

This work sheet "ad absurdum QED" is my conclusive conjecture:
1. You can't pixelate w/o the circular convolution,
2. over the kernel 1
3. and not w/o the cyclic permutation "Center".

Note that here for the demo I have used my "Integer Zoom". For not integer zoom, the Mathcad built-in Improc zoom would be needed. For a little while, Eden resample2D seemed to be the one (Improc zoom), but it is not. For sure resample2D does pixelate, but crappy ... so does not technically pixelate.

Jean

... same work sheet but collapsed conveniently. Now you can click on rnd. Experiment one step further, select example 2 and see about same pixelation (visually). You can see my old point about real image rather than meaningless "spots". With all that material in hand, you could try convolving resample2D with the kernel and center. It may be the missing link that would make everybody happy. Bed time.

Jean
PhilipOakley
5-Regular Member
(To:ptc-1368288)

Jason,
Just to complement Jean's comment on re-centering:

The fourier transform comes in a few 'flavours'. In the school text style you have separate sine and cosine terms, and a DC term.

When you have numbers that are complex, you then start to have 'negative frequencies' and complex amplitutes.

What this means is that when you have N samples you have frequencies from 0 to N/2. Then, depending on which function you call, and in particular the ones that process 2d data, you still have values in the arry from N/2 to N-1. When the input is real, these upper values are just the complex complement of the values from 0 to N/2 in reverse order i.e. the value at 1 has its complex complement at -1 -> N-1.

Re-centering is wher we shift the 0,0 point of the array to the middkle so we can see the symmetry.
This reflection also explains why some of the zero insersions are not immediately obvious.

As an aside, when the numbers are complex, those sine and cosines begin to look like spirals along the t-axis when they are plotted as an argand plane square to the t-axis. The negative/positive frequencies are then simply clockwise/anticlockwise rotations.

Philip Oakley

I wanted to compare this frequency domain method to a previous method (that you have all helped with and I am very grateful) in which I fit a function to the data in the space domain. The desired metric at the end (encircled energy) shows a slight discrepancy between the two methods.

Is this expected?

Thanks.
Jason

From yesterday, it came clear that resample2D does not work properly for real images, why should it work for your asymmetric spot ? It does not work properly because it does not pixelate as it should. Maybe your spot is not asymmetric, only not centered ? The project starts by analyzing the spot at some point of its original data set, or as you seem to indicate : analyze past the fitting ... a complicated fitting that I respect but my first preference would go for Tom's method. If you don't have a circle, the circle will be an ellipse. But for all that the route is via the diagonal spline interpolator, CreateMesh, Given/Find for the contour, fit circle/ellipse, then integrate for the volume.

By reading, I couldn't follow your 2d integrator.
I think the project must be started from the very beginning, or at least from some confirmed steps that if you like those confirmed steps, just start at the resulting something. What I'm saying here is: plug an image in a work sheet and say find the overall volume, find the cutting circle/ellipse that will give so much proportion of the total volume.

Same work sheet just made more visible.

Jean
LouP
11-Garnet
(To:Jbryant61)

I can think of two structural sources of errors in addition to any numerical approximation errors:

1) The freq domain analysis integrates an interpolated version of the original image, while the spatial domain analysis integrates the fitted function. these are not identical, so variations in integral vs. radius should not be surprising. I plotted the interpolated image and the fit, and they are close, but certainly not identical.

2) The area integrals vs. radius do not converge rapidly to a final value within the domain of the image (non-negligible area in the outer "dark" region), so that the radius used to define the normalization matters. For the 2D DFT method, the function was simply offset by its global minimum before integration(via scale fct.), leaving all contributions in the "dark" region to be positive.

Lou
Top Tags