Guest User

Untitled

a guest
Jun 25th, 2018
71
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.78 KB | None | 0 0
  1. from matplotlib import cm,pyplot as plt
  2. from scipy.misc import *
  3. from sys import exit,argv
  4. from numpy import *
  5. from scipy.optimize import fmin
  6. from skimage.io import imread
  7.  
  8.  
  9. ###################################################
  10. #ROUTINES
  11. ###################################################
  12. def imageCentroid(imagen,threshold=50):
  13. height,width=imagen.shape
  14. cond=imagen>=threshold
  15. X,Y=meshgrid(arange(width),arange(height))
  16. xm=X[cond].mean()
  17. ym=Y[cond].mean()
  18. return xm,ym
  19.  
  20. def detectBorder(imagen,threshold=50,R=900):
  21. cond=imagen>=threshold
  22. xm,ym=imageCentroid(imagen,threshold=threshold)
  23. plt.close("all")
  24. contour=plt.contour(cond,levels=[0.2])
  25. border=[]
  26. for path in contour.collections[0].get_paths():
  27. for points in path.vertices:
  28. if abs(points[0]-xm)>R or abs(points[1]-ym)>R:
  29. continue
  30. border+=[points.tolist()]
  31. border=array(border)
  32. return xm,ym,border
  33.  
  34. ###################################################
  35. #PROCESSING IMAGES
  36. ###################################################
  37. filename=("moon.JPG")
  38. basename=filename.split(".")[0]
  39. print ("Processing image %s..."%filename)
  40.  
  41. #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  42. #LOAD IMAGE
  43. #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  44. imagen=imread(filename)[::,0]
  45. height,width=imagen.shape
  46. print ("Image size:"),width,height
  47.  
  48. #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  49. #DETECT CENTROID AND BORDER
  50. #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  51. def errorR(params):
  52. thresh=params[0]
  53. xm,ym,border=detectBorder(imagen,threshold=thresh)
  54. xborder=border[:,0]
  55. yborder=border[:,1]
  56. xm=xborder.mean()
  57. ym=yborder.mean()
  58. rborder=((xborder-xm)**2+(yborder-ym)**2)**0.5
  59. R=rborder.mean()
  60. dR=rborder.std()
  61. print ("tThreshold: %.1f, Measured Radius : %.4f +/- %.2f")%(thresh,R,dR)
  62. return dR
  63.  
  64. print ("Minimizing uncertainties:")
  65. solution=fmin(errorR,[80.0],xtol=0.5,maxiter=10)
  66. thrmin=solution[0]
  67. print ("Threshold:"),thrmin
  68.  
  69. xmc,ymc,borderc=detectBorder(imagen,threshold=thrmin)
  70. dRmin=errorR([thrmin])
  71. rborder=((borderc[:,0]-xmc)**2+(borderc[:,1]-ymc)**2)**0.5
  72. Rmin=rborder.mean()
  73. print ("Radius at minimum uncertainty: %.4f +/- %.2f")%(Rmin,dRmin)
  74.  
  75. ###################################################
  76. #PLOT IMAGE
  77. ###################################################
  78. xm,ym,borderm=detectBorder(imagen,threshold=20)
  79. xm,ym,borderM=detectBorder(imagen,threshold=120)
  80.  
  81. fig=plt.figure(figsize=(8,8))
  82. ax=fig.add_axes([0.0,0.0,1.0,1.0])
  83. plt.imshow(imagen,cmap=cm.gray)
  84. border=borderc
  85. plt.plot(border[::,0],border[::,1],'c.',markersize=0.5)
  86. border=borderm
  87. plt.plot(border[::,0],border[::,1],'y.',markersize=0.5)
  88. border=borderM
  89. plt.plot(border[::,0],border[::,1],'r.',markersize=0.5)
  90. plt.xlim((xmc-1.2*Rmin,xmc+1.2*Rmin));plt.ylim((ymc-1.2*Rmin,ymc+1.2*Rmin))
  91. plt.savefig(basename+"-border.png")
Add Comment
Please, Sign In to add comment