numpy - How do I fit a quadratic surface to some points in Python? -


i have points in 3-space , i'd fit quadratic surface through them.

i tried code 

import itertools import numpy np import matplotlib.pyplot plt   def main():     points = [ [ 175697888, -411724928, 0.429621160030365 ], [ 175697888, -411725144, 0.6078286170959473 ], [ 175698072, -411724640, 0.060898926109075546 ], [ 175698008, -411725360, 0.6184252500534058 ], [ 175698248, -411725720, 0.0771455243229866 ], [ 175698448, -411724456, -0.5925689935684204 ], [ 175698432, -411725936, -0.17584866285324097 ], [ 175698608, -411726152, -0.24736160039901733 ], [ 175698840, -411724360, -1.27967369556427 ], [ 175698800, -411726440, -0.21100902557373047 ], [ 175699016, -411726744, -0.12785470485687256 ], [ 175699280, -411724208, -2.472576856613159 ], [ 175699536, -411726688, -0.19858847558498383 ], [ 175699760, -411724104, -3.5765910148620605 ], [ 175699976, -411726504, -0.7432857155799866 ], [ 175700224, -411723960, -4.770215034484863 ], [ 175700368, -411726304, -1.2959377765655518 ], [ 175700688, -411723760, -6.518451690673828 ], [ 175700848, -411726080, -3.02254056930542 ], [ 175701160, -411723744, -7.941056251525879 ], [ 175701112, -411725896, -3.884831428527832 ], [ 175701448, -411723824, -8.661275863647461 ], [ 175701384, -411725720, -5.21607780456543 ], [ 175701704, -411725496, -6.181706428527832 ], [ 175701800, -411724096, -9.490276336669922 ], [ 175702072, -411724344, -10.066594123840332 ], [ 175702216, -411724560, -10.098011016845703 ], [ 175702256, -411724864, -9.619892120361328 ], [ 175702032, -411725160, -6.936516284942627 ] ]      n = len(points)     x, y, z = map(np.array, zip(*points))      plt.figure()     plt.subplot(1, 1, 1)      # fit 3rd order, 2d polynomial     m = polyfit2d(x,y,z, order=2)      # evaluate on grid...     nx, ny = 100, 100     xx, yy = np.meshgrid(np.linspace(x.min(), x.max(), nx), np.linspace(y.min(), y.max(), ny))     zz = polyval2d(xx, yy, m)      plt.scatter(xx, yy, c=zz, marker=2)     plt.scatter(x, y, c=z)     plt.show()  def polyfit2d(x, y, z, order=2):     ncols = (order + 1)**2     g = np.zeros((x.size, ncols))     ij = itertools.product(range(order+1), range(order+1))     k, (i,j) in enumerate(ij):         g[:,k] = x**i * y**j     m, _, _, _ = np.linalg.lstsq(g, z)     return m  def polyval2d(x, y, m):     order = int(np.sqrt(len(m))) - 1     ij = itertools.product(range(order+1), range(order+1))     z = np.zeros_like(x)     a, (i,j) in zip(m, ij):         z += * x**i * y**j     return z  main()

based on answer: python 3d polynomial surface fit, order dependent

but gives opposite result:

enter image description here

look @ colour of points compared surface. idea i'm doing wrong?

edit: update code remove imshow showing isn't issue.

there seems problem floating point accuracy. played code bit , change range of x , y made least square solution work. doing

x, y = x - x[0], y - y[0]

solved accuracy issue. can try:

import itertools import numpy np import matplotlib.pyplot plt mpl_toolkits.mplot3d import axes3d # matplotlib import cbook matplotlib import cm matplotlib.colors import lightsource   def poly_matrix(x, y, order=2):     """ generate matrix use lstsq """     ncols = (order + 1)**2     g = np.zeros((x.size, ncols))     ij = itertools.product(range(order+1), range(order+1))     k, (i, j) in enumerate(ij):         g[:, k] = x**i * y**j     return g   points = np.array([[175697888, -411724928, 0.429621160030365],                    [175697888, -411725144, 0.6078286170959473],                    [175698072, -411724640, 0.060898926109075546],                    [175698008, -411725360, 0.6184252500534058],                    [175698248, -411725720, 0.0771455243229866],                    [175698448, -411724456, -0.5925689935684204],                    [175698432, -411725936, -0.17584866285324097],                    [175698608, -411726152, -0.24736160039901733],                    [175698840, -411724360, -1.27967369556427],                    [175698800, -411726440, -0.21100902557373047],                    [175699016, -411726744, -0.12785470485687256],                    [175699280, -411724208, -2.472576856613159],                    [175699536, -411726688, -0.19858847558498383],                    [175699760, -411724104, -3.5765910148620605],                    [175699976, -411726504, -0.7432857155799866],                    [175700224, -411723960, -4.770215034484863],                    [175700368, -411726304, -1.2959377765655518],                    [175700688, -411723760, -6.518451690673828],                    [175700848, -411726080, -3.02254056930542],                    [175701160, -411723744, -7.941056251525879],                    [175701112, -411725896, -3.884831428527832],                    [175701448, -411723824, -8.661275863647461],                    [175701384, -411725720, -5.21607780456543],                    [175701704, -411725496, -6.181706428527832],                    [175701800, -411724096, -9.490276336669922],                    [175702072, -411724344, -10.066594123840332],                    [175702216, -411724560, -10.098011016845703],                    [175702256, -411724864, -9.619892120361328],                    [175702032, -411725160, -6.936516284942627]])  ordr = 2  # order of polynomial x, y, z = points.t x, y = x - x[0], y - y[0]  # improves accuracy  # make matrix: g = poly_matrix(x, y, ordr) # solve np.dot(g, m) = z: m = np.linalg.lstsq(g, z)[0]   # evaluate on grid... nx, ny = 30, 30 xx, yy = np.meshgrid(np.linspace(x.min(), x.max(), nx),                      np.linspace(y.min(), y.max(), ny)) gg = poly_matrix(xx.ravel(), yy.ravel(), ordr) zz = np.reshape(np.dot(gg, m), xx.shape)  # plotting (see http://matplotlib.org/examples/mplot3d/custom_shaded_3d_surface.html): fg, ax = plt.subplots(subplot_kw=dict(projection='3d')) ls = lightsource(270, 45) rgb = ls.shade(zz, cmap=cm.gist_earth, vert_exag=0.1, blend_mode='soft') surf = ax.plot_surface(xx, yy, zz, rstride=1, cstride=1, facecolors=rgb,                        linewidth=0, antialiased=false, shade=false) ax.plot3d(x, y, z, "o")  fg.canvas.draw() plt.show()

which gives 3dresult plot]

to evaluate quality of fit read documentation np.linalg.lstsq(). rank should size of result vector , residual divided number of data points gives average error (distance between point , plane).


Comments

Post a Comment

Popular posts from this blog

authentication - Mongodb revoke acccess to connect test database -

i have ec2 instance mongodb 3.2.0 installed. have user admin access. can connect user , perform operation remotely. 1 can connect database test remote if knows ip mongo 11.11.11.11. although not able make operations, want restrict access user credentials can connect server. first need create users , grant roles users, example: use admin db.createuser( { user: "myuseradmin", pwd: "abc123", roles: [ { role: "useradminanydatabase", db: "admin" } ] } ) then need start server --auth parameter mongod --auth more information here
Read more

r - Update two sets of radiobuttons reactively - shiny -

i have read ( how make choices in radiobuttons reactive in shiny? ) shows me how update radiobutton s in reactive way. however, when try , update 2 sets of buttons same data, 1 set renders. example: server: # create example data wafer <- rep(c(1:3), each=3) length <- c(1,1,2,1,1,1,3,5,1) width <- c(3,1,6,1,1,1,1,1,6) dd <- data.frame(wafer, length, width) shinyserver(function(input, output, session){ # create reactive dataframe store data values <- reactivevalues() values$df <- data.frame() # lengths , widths of wafer user input <- eventreactive(input$do, { subset(dd, wafer %in% input$wafer, select = length:width) }) # update reactive data frame lengths , widths have been selected user input observe({ if(!is.null(a())) { values$df <- rbind(isolate(values$df), a()) } }) output$wl <- renderdatatable({ a() }) # update radio buttons unique length , widths stored in values$df # ever "observe" put first in co...
Read more

ios - Realm over CoreData should I use NSFetchedResultController or a Dictionary? -

i'm working on app using core data , nsfetchedresultcontroller many reasons switch use realm. saw there "realm version" of nsfetchedresultcontroller on github wouldn't compatible current code using core data. the view controller displaying list of people same school, address book. list sublist of people studied in same city. so thinking make 1 unique request database retrieve list of people , filter locally list within 1 dictionary per school [string: anyobject] string section name within tableview , anyobject array of people. first solution, make nsfetchedresultcontroller each school predicate filtering location. delete actions etc.. handled delegates -> not compatible realm create dictionaries , update them each actions... -> works realm it's annoying code. any better solution? edit: i need clarify request: i'd write class inherit uitableviewcontroller. this class has list of people sorted in alphabetical order the tableview has ...
Read more