fit_ellipse.py

# -*- coding: utf-8 -*-
import numpy
def fitEllipse(x,y):

    #x=cont[:,0]
    #y=cont[:,1]
    x=x[:,None]
    y=y[:,None]
    D=numpy.hstack([x*x,x*y,y*y,x,y,numpy.ones(x.shape)])
    S=numpy.matmul(D.T,D)
    C=numpy.zeros([6,6])
    C[0,2]=C[2,0]=2
    C[1,1]=-1
    invS=numpy.linalg.inv(S)
    invSdotC=numpy.matmul(invS,C)
    E,V=numpy.linalg.eig(invSdotC)
    n=numpy.argmax(numpy.abs(E))
    a=V[:,n]
    #-------------------Fit ellipse-------------------
    b,c,d,f,g,a=a[1]/2., a[2], a[3]/2., a[4]/2., a[5], a[0]
    num=b*b-a*c
    cx=(c*d-b*f)/num
    cy=(a*f-b*d)/num
    angle=0.5*numpy.arctan(2*b/(a-c))*180/numpy.pi
    up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
    down1=(b*b-a*c)*( (c-a)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    down2=(b*b-a*c)*( (a-c)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
    a=numpy.sqrt(abs(up/down1))
    b=numpy.sqrt(abs(up/down2))
    params=[cx,cy,a,b,angle]
    return params

x=[1229,1230,1231,1231,1232,1232,1233,1233,1234,1234,1235,1235,1236,1236,1237,1237,1238,1238,1239,1239,1240,1240,1241,1241,1242,1242,1243,1244,1244,1245,1245,1246,1246,1247,1247,1248,1248,1249,1249,1250,1250,1251,1252,1252,1253,1253,1254,1254,1255,1255,1256,1256,1257,1258,1258,1259,1259,1260,1260,1261,1261,1262,1263,1263,1264,1264,1265,1265,1266,1266,1267,1268,1268,1269,1269,1270,1270,1271,1272,1272,1273,1273,1274,1275,1275,1276,1276,1277,1277,1278,1279,1279,1280,1280,1282,1282,1282,1283,1283,1284,1285,1285,1286,1286,1287,1288,1288,1289,1289,1290,1291,1291,1292,1292,1293,1294,1294,1295,1296,1296,1297,1297,1298,1299,1299,1300,1301,1301,1302,1302,1303,1304,1304,1305,1306,1306,1307,1308,1308,1309,1309,1310,1311,1311,1312,1313,1313,1314,1315,1315,1316,1317,1318,1319,1319,1320,1321,1321,1322,1323,1323,1324,1325,1326,1327,1328,1328,1329,1330,1330,1331,1332,1333,1334,1335,1336,1337,1337,1338,1339,1340,1340,1341,1342,1342,1343,1344,1345,1346,1347,1348,1349,1350,1351,1352,1353,1354,1355,1356,1357,1358,1359,1360,1361,1362,1363,1364,1365,1366,1367,1368,1369,1370,1371,1372,1373,1374,1375,1376,1377,1378,1379,1380,1381,1382,1383,1384,1384,1385,1385,1386,1387,1388,1389,1390,1391,1392,1393,1394,1395,1396,1397,1398,1399,1400,1401,1402,1403,1404,1405,1406,1407,1408,1409,1410,1411,1412,1413,1414,1419,1420,1421,1422,1423,1424,1425,1426,1427,1428,1429,1430,1431,1432,1433,1434,1435,1436,1437,1438,1439,1440,1441,1442,1443,1444,1445,1446,1447,1448,1449,1450,1451,1451,1452,1453,1453,1454,1455,1455,1456,1457,1457,1458,1458,1459,1460,1460,1461,1461,1462,1463,1463,1464,1464,1465,1465,1466,1465,1465,1466,1466,1467,1468,1468,1469,1469,1470,1470,1471,1471,1472,1472,1473,1473,1474,1474,1474,1475,1475,1476,1476,1477,1477,1478,1478,1478,1479,1479,1480,1480,1480,1481,1481,1482,1482,1482,1483,1484,1485,1485,1485,1486,1486,1487,1487,1487,1487,1488,1488,1489,1489,1489,1490,1490,1490,1491,1491,1491,1491,1492,1492,1492,1492,1492,1493,1493,1493,1493,1493,1494,1494,1494,1494,1494,1494,1494,1495,1495,1491,1491,1491,1491,1492,1492,1492,1492,1493,1493,1493,1493,1493,1493,1494,1494,1494,1494,1494,1494,1494,1494,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1495,1484,1485,1485,1485,1486,1486,1486,1487,1487,1487,1488,1488,1488,1488,1489,1489,1489,1490,1490,1490,1490,1491,1461,1461,1462,1462,1463,1463,1463,1464,1464,1465,1465,1465,1466,1466,1466,1467,1467,1467,1468,1468,1468,1469,1469,1469,1470,1470,1471,1471,1472,1472,1473,1473,1474,1474,1474,1475,1475,1475,1475,1476,1476,1476,1477,1477,1478,1478,1478,1479,1479,1480,1480,1480,1481,1481,1481,1482,1482,1483,]
y=[432,435,434,435,433,434,432,433,431,432,430,431,429,430,428,429,427,428,426,427,425,426,424,425,423,424,423,422,423,421,422,420,421,419,420,418,419,417,418,416,417,416,415,416,414,415,413,414,412,413,411,412,411,410,411,409,410,408,409,407,408,407,406,407,405,406,404,405,403,404,403,402,403,401,402,400,401,400,399,400,398,399,398,397,398,396,397,395,396,395,394,395,393,394,389,392,393,391,392,391,390,391,389,390,389,388,389,387,388,387,386,387,385,386,385,384,385,384,383,384,382,383,382,381,382,381,380,381,379,380,379,378,379,378,377,378,377,376,377,375,376,375,374,375,374,373,374,373,372,373,372,371,371,370,371,370,369,370,369,368,369,368,368,367,367,366,367,366,365,366,365,364,364,364,363,363,362,363,362,362,361,361,361,360,361,360,360,359,359,359,358,358,358,357,357,357,356,356,356,355,355,355,355,354,354,354,353,353,353,353,352,352,352,352,351,351,351,351,351,350,350,350,350,350,349,349,346,349,346,349,349,349,349,349,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,348,349,349,349,349,349,349,349,350,350,350,350,351,351,351,351,352,352,352,352,353,353,353,354,354,355,355,356,356,356,357,357,358,358,358,359,359,359,360,360,360,361,361,362,362,362,363,363,364,364,364,365,365,366,366,367,367,366,367,367,368,368,368,369,369,370,370,371,371,372,372,373,373,374,374,375,376,376,377,377,378,378,379,379,380,381,381,382,382,383,384,384,385,385,386,387,387,389,389,390,391,391,392,387,393,394,395,396,397,398,399,400,401,402,403,404,405,405,406,407,408,409,410,411,412,413,414,415,416,417,418,419,420,421,422,423,424,425,468,469,470,471,464,465,466,467,458,459,460,461,462,463,450,451,452,453,454,455,456,457,427,428,429,430,431,432,433,434,435,436,437,438,439,440,441,442,443,444,445,446,447,448,449,493,490,491,492,487,488,489,484,485,486,480,481,482,483,477,478,479,473,474,475,476,472,539,540,538,539,536,537,538,535,536,533,534,535,531,532,533,529,530,531,527,528,529,525,526,527,524,525,521,523,519,520,517,518,515,516,517,512,513,514,515,510,511,512,509,510,506,507,508,504,505,501,502,503,498,499,500,496,497,495,]

fitEllipse(numpy.array(x),numpy.array(y) )