Fit a line to a matrix in python
Question: I have a matrix of shape 256x256 to which I'm trying to find a line of best
fit. This is an image by the way, so these are simply intensity values. Let's
assume I want to find the line of best fit through all the intensities, how
would I go about doing that?
[This](http://stackoverflow.com/questions/2298390/fitting-a-line-in-3d) link
describes how to do so on a 3d dataset, using svd. However I'm a bit confused
as to how that would be applied to my numpy array?

EDIT: Here's an example with random values of doubles that I've profiled with
%timeit:

    
    
    ran = [25,50,75,100,125]
    for i in ran:                                                                                 
        J = np.random.random((i,i))                                                     
        y,x=np.indices(J.shape)                                                                   
        x = x.ravel()                                                                             
        y = y.ravel()                                                                             
        J=J.ravel()                                                                               
        data = np.concatenate((x[:, np.newaxis],
                              y[:, np.newaxis],
                              J[:, np.newaxis]),axis=1)        
        datamean = data.mean(axis=0)                                                              
        print "Doing %d now"  %i                                                                    
        %timeit U, S, V = np.linalg.svd(data - datamean) 
    

I get the following output:

    
    
    Doing 25 now
    100 loops, best of 3: 10.4 ms per loop
    Doing 50 now
    1 loops, best of 3: 285 ms per loop
    Doing 75 now
    1 loops, best of 3: 3 s per loop
    Doing 100 now
    1 loops, best of 3: 5.83 s per loop
    Doing 125 now
    1 loops, best of 3: 15.1 s per loop
    

EDIT2: [Here's my actual
array](https://www.dropbox.com/s/8a8zj686j4k0dfn/arr.npy?dl=0). I just saved
it in numpy's npy format


Answer: [The answer you pointed out](http://stackoverflow.com/a/2333251/832621) is
directly applicable to your problem by doing:

    
    
    import numpy as np
    
    z = your_matrix_256_x_256
    y, x = np.indices(z.shape)
    x = x.ravel()
    y = y.ravel()
    z = z.ravel()
    

Note that the intervals for `x` and `y` can be reajusted multiplying these
arrays by proper scalars.

* * *

## EDIT:

having a look to your data it seems more that **your problem is a 2-D curve
fit** , which can be done with `np.polyfit()`, as shown in the example below.

    
    
    z = np.load('arr.npy').astype(np.int32)
    y, x = np.indices(z.shape)
    valid_z = (y.ravel()>0) & (z.ravel()>0)
    x_valid = x.ravel()[valid_z]
    y_valid = y.ravel()[valid_z]
    z_valid = z.ravel()[valid_z]
    # fitting best curve
    fig = plt.figure(figsize=(5,3))
    ax = fig.add_subplot(111)
    ax.scatter(x_valid, y_valid, c=z_valid, alpha=0.2, s=20, edgecolor='none',
            cmap=plt.cm.jet)
    # finding best-fit curve
    z = np.polyfit(x_valid, y_valid, w=z_valid**0.5, deg=1)
    p = np.poly1d(z)
    # plotting
    x_plot = np.linspace(x_valid.min(), x_valid.max(), 100)
    y_plot = p(x_plot)
    ax.plot(x_plot, y_plot, '-r', lw=2)
    ax.set_xlim(0, x.shape[1])
    ax.set_ylim(0, y.shape[0])
    
    ax.legend(loc='lower left', frameon=False, fontsize=8)
    fig.savefig('test.png', bbox_inches='tight')
    

which gives:

![enter image description here](http://i.stack.imgur.com/3Yu6w.png)

