| Re: Polynomial Fitting question [message #68139] |
Sat, 03 October 2009 09:15  |
Andrew Rodger
Messages: 5
Registered: October 2009
|
Junior Member |
|
|
OK so I am a nerd and worked on this at home, over the weekend :(
But I have my solution now and it does not require an inner FOR loop
and could be extended to encompass Gaussian curves as well.
I went back to basics and used the tried and true least squares
fitting of the following form
coefficients=INVERSE(TRANSPOSE(X)X)TRANSPOSE(X)Y I didn't really
know how to write this without using some pseudo syntax.
Anyhoo, the little code snippet below is what I used to test it. I am
sure that I have used the # and ## incorrectly and I have also
discovered that I could have used matrix_multiply to do the same thing
in a lot of cases (I may still do that). The thing I like about this
is that the values BIG_X and TRANSPOSE(x) can be calculated outside of
any loop and I only need to calculate XtY within the original outer
loop (BIG_Y is simply a one liner that is supplied via an ASSOC
variable). So I should only need two lines (the last two) within the
outer loop to get my coefficients on a line by line and pixel by pixel
basis.
PRO test,coefficients
;Some X array: These represent my wavelengths which will never change
in a given image
;(12,1) array
x_sing=[890,900,910,920,930,940,950,960,970,980,990,1000]/10 00.
;(3,12) array
x=double(transpose([[fltarr(12)+1],[x_sing],[x_sing^2]]))
;(3,3) array
temp_x=x#transpose(x)
;(3,3) array
BIG_X=REAL_PART(lu_complex(temp_x,/inverse))
;make some Y arrays. These represent the surface reflectance for a
given pixel
;in this case 5 pixels. In actuality they will be read in 512 pixels
at a time from my image data
;(12,1) array
y1=0.5+0.3*x_sing+0.77*x_sing^2
;(12,1) array
y2=2.0+0.5*x_sing+1.2*x_sing^2
;(12,1) array
y3=4.0+1.5*x_sing+0.2*x_sing^2
;(12,1) array
y4=10.0+0.6*x_sing+3.0*x_sing^2
;(12,1) array
y5=1.0+2.0*x_sing+3.0*x_sing^2
;(5,12) array
BIG_Y=TRANSPOSE([[y1],[y2],[y3],[y4],[y5]])
;(5,3) array
XtY=TRANSPOSE(x)##BIG_Y
;These are the coefficients that I seek
coefficients=BIG_X##xty
END
Cheers
Andrew
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
