| Re: fft and least sqaues problem [message #81102] |
Fri, 17 August 2012 07:40  |
rogass
Messages: 200
Registered: April 2008
|
Senior Member |
|
|
Am Freitag, 17. August 2012 06:58:12 UTC+2 schrieb Craig Markwardt:
> On Thursday, August 16, 2012 5:13:49 AM UTC-4, (unknown) wrote:
>
>> Hi,
>
>>
>
>> the following was missing:
>
>
>
>
>
> Cool, thanks!
Hi Craig,
do you have a routine to perform iteratively reweighted least squares?
Cheers
CR
|
|
|
|
|
|
| Re: fft and least sqaues problem [message #81115 is a reply to message #81105] |
Thu, 16 August 2012 02:13  |
rogass
Messages: 200
Registered: April 2008
|
Senior Member |
|
|
Am Mittwoch, 15. August 2012 01:53:13 UTC+2 schrieb Craig Markwardt:
> On Tuesday, August 14, 2012 2:22:11 PM UTC-4, (unknown) wrote:
>
>> Finally, I found the mistake - sitting in front of this pc :)
>
>
>
> OK.
>
> So summarize for us: what was wrong, and what did you change to fix it?
>
>
>
> Craig
Hi,
the following was missing:
wx-=s[0]/2
wy-=s[1]/2
However, the solution is noise sensitive. Smoothing the data beforehand helps :)
Changing:
phase=atan(imaginary(corr)/real_part(corr))
to:
phase=atan(corr,/phase)
additionally reduces processing time. Thanks to Helder.
Cheers
CR
|
|
|
|
|
|
| Re: fft and least sqaues problem [message #81141 is a reply to message #81135] |
Tue, 14 August 2012 11:22 |
rogass
Messages: 200
Registered: April 2008
|
Senior Member |
|
|
Finally, I found the mistake - sitting in front of this pc :)
Here is the code:
s=size(im,/dim)*1.
im1 =im
dx=.25
dy=.3
im2=image_shift(im1,dx,dy);more precise then interpolate(im,findgen(s[0])+dx,findgen(s[1])+dy,/grid,/cubi c)
fim1=fft(im1,-1)
fim2=fft(im2,-1)
corr=fim1*conj(fim2)/abs(fim1*fim2)
corr=shift(corr,s/2)
phase=atan(imaginary(corr)/real_part(corr))
wx=(findgen(s) mod s[0])
wy=(rebin(findgen(1,s[1]),s))
wx-=s[0]/2
wy-=s[1]/2
wx*=2.*!pi/s[0]
wy*=2.*!pi/s[1]
r=5;fitting radius
phase2=phase[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
wx2=wx[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
wy2=wy[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
print,la_least_squares(transpose([[wx2[*]],[wy2[*]]]),phase2 [*])
Thanks Helder!
Cheers
CR
|
|
|
|
| Re: fft and least sqaues problem [message #81142 is a reply to message #81141] |
Tue, 14 August 2012 07:50 |
Helder Marchetto
Messages: 520
Registered: November 2011
|
Senior Member |
|
|
On Tuesday, August 14, 2012 4:06:53 PM UTC+2, (unknown) wrote:
> Hi Folks,
>
>
>
> I try to estimate the subpixelshifts if an image is compared with its shifted represenetation, but something is going wrong. Maybe somebody can help me.
>
>
>
> s=size(im,/dim)*1.
>
> im1 =im
>
> dx=.25
>
> dy=.3
>
> im2=image_shift(im1,dx,dy);more precise then interpolate(im,findgen(s[0])+dx,findgen(s[1])+dy,/grid,/cubi c)
>
> fim1=fft(im1,-1)
>
> fim2=fft(im2,-1)
>
> corr=fim1*conj(fim1)/abs(fim1*fim2)
>
> corr=shift(corr,s/2)
>
> phase=atan(imaginary(corr)/real_part(corr))
>
> wx=(findgen(s) mod s[0])*2.*!pi/s[0]
>
> wy=(rebin(findgen(1,s[1]),s))*2.*!pi/s[1]
>
> r=5;fitting radius
>
> phase2=phase[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
>
> wx2=wx[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
>
> wy2=wy[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
>
> print,la_least_squares(transpose([[wx2[*]],[wy2[*]]]),phase2 [*])
>
>
>
> The last line should give dx and dy but its erroneous. I don't really know why!
>
>
>
> Thanks in advance
>
>
>
> CR
Maybe it's a typo, but you defined the phase-correlation "corr" as:
corr=fim1*conj(fim1)/abs(fim1*fim2)
instead of:
corr=fim1*conj(fim2)/abs(fim1*fim2)
Try that...
Did you have a look at: http://en.wikipedia.org/wiki/Phase_correlation ?
Cheers,
Helder
|
|
|
|
| Re: fft and least sqaues problem [message #81143 is a reply to message #81142] |
Tue, 14 August 2012 07:13 |
rogass
Messages: 200
Registered: April 2008
|
Senior Member |
|
|
Am Dienstag, 14. August 2012 16:06:53 UTC+2 schrieb (unbekannt):
> Hi Folks,
>
>
>
> I try to estimate the subpixelshifts if an image is compared with its shifted represenetation, but something is going wrong. Maybe somebody can help me.
>
>
>
> s=size(im,/dim)*1.
>
> im1 =im
>
> dx=.25
>
> dy=.3
>
> im2=image_shift(im1,dx,dy);more precise then interpolate(im,findgen(s[0])+dx,findgen(s[1])+dy,/grid,/cubi c)
>
> fim1=fft(im1,-1)
>
> fim2=fft(im2,-1)
>
> corr=fim1*conj(fim1)/abs(fim1*fim2)
>
> corr=shift(corr,s/2)
>
> phase=atan(imaginary(corr)/real_part(corr))
>
> wx=(findgen(s) mod s[0])*2.*!pi/s[0]
>
> wy=(rebin(findgen(1,s[1]),s))*2.*!pi/s[1]
>
> r=5;fitting radius
>
> phase2=phase[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
>
> wx2=wx[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
>
> wy2=wy[s[0]/2 - r : s[0]/2 + r,s[1]/2 - r : s[1]/2 + r]
>
> print,la_least_squares(transpose([[wx2[*]],[wy2[*]]]),phase2 [*])
>
>
>
> The last line should give dx and dy but its erroneous. I don't really know why!
>
>
>
> Thanks in advance
>
>
>
> CR
Typo!
Please replace:
im2=image_shift(im1,dx,dy)
with:
im2=image_shift(im1,dx,dy,interp_type='F')
Cheers
CR
|
|
|
|
| Re: fft and least sqaues problem [message #81198 is a reply to message #81102] |
Fri, 17 August 2012 08:45 |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
On Friday, August 17, 2012 10:40:21 AM UTC-4, CR wrote:
...
> Hi Craig,
>
> do you have a routine to perform iteratively reweighted least squares?
You're kind of hijacking an existing conversation. The short answer is, no, nothing specific to IRLS, but there's nothing stopping you from doing the easy iterative part yourself.
Craig
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
