Skip to content

Commit

Permalink
Merge pull request #355 from jonwright/master
Browse files Browse the repository at this point in the history 
fix geometry to match module and cleanup and fit a CeO2
  • Loading branch information
@jonwright
jonwright authored Nov 20, 2024
2 parents e6e5209 + b5f8a20 commit bd380f4
Show file tree
Hide file tree
Showing 2 changed files with 555 additions and 37 deletions.
101 changes: 64 additions & 37 deletions ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,18 +115,9 @@
"outputs": [],
"source": [
"# powder patterns\n",
"hsI, dsbinedges = np.histogram( cf_4d.ds, weights=cf_4d.sum_intensity, bins=np.linspace(0, cf_4d.ds.max(), 5000))\n",
"hs1, dsbinedges = np.histogram( cf_4d.ds, bins=np.linspace(0, cf_4d.ds.max(), 5000))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"pbins = 5000\n",
"hsI, dsbinedges = np.histogram( cf_4d.ds, weights=cf_4d.sum_intensity, bins=np.linspace(0, cf_4d.ds.max(), pbins))\n",
"hs1, dsbinedges = np.histogram( cf_4d.ds, bins=np.linspace(0, cf_4d.ds.max(), pbins))\n",
"dsbincens = 0.5*(dsbinedges[1:]+dsbinedges[:-1])"
]
},
Expand Down Expand Up @@ -166,7 +157,9 @@
"outputs": [],
"source": [
"# These peaks are already indexed by the main phase\n",
"ringmask = ImageD11.peakselect.rings_mask( cf_4d, 0.005, 1.0, cell = ucell )"
"dstar_tol = 0.005 # range in dstar to mask\n",
"dsmax = 1.0 # max dstar (highest angle)\n",
"ringmask = ImageD11.peakselect.rings_mask( cf_4d, dstar_tol, 1.0, cell = ucell )"
]
},
{
Expand All @@ -178,7 +171,7 @@
"outputs": [],
"source": [
"# cring = columnfile of things not indexed by the main rings\n",
"cring = cf_4d.copyrows( (cf_4d.ds < 1.)&( ~ringmask ))"
"cring = cf_4d.copyrows( (cf_4d.ds < dsmax)&( ~ringmask ))"
]
},
{
Expand All @@ -193,7 +186,7 @@
"ax.plot( dsbincens, hsI/hsI.max() , '-')\n",
"ax.plot( ucell.ringds, np.full( len(ucell.ringds), 1e-5) , \"|\", ms=20, lw =1 )\n",
"ax.plot( cring.ds, cring.sum_intensity/cring.sum_intensity.max(), \".\", ms=1 )\n",
"ax.set(yscale='log', xlim=(0.1,1), ylim=(1e-5,1))"
"ax.set(yscale='log', xlim=(0.1,1), ylim=(1e-5,1));"
]
},
{
Expand Down Expand Up @@ -240,22 +233,36 @@
},
"outputs": [],
"source": [
"# Position where the beam goes through the rotation axis\n",
"y0 = 0.0"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# Interactive plot to select a peak\n",
"\n",
"fig, ax = plt.subplots(1,2)\n",
"ax[0].scatter( cpk.omega, cpk.dty, \n",
" c=np.log(cpk.sum_intensity) * np.sign( cpk.yl ), \n",
" s=4, cmap='RdBu', picker=True, pickradius = 5)\n",
"\n",
"selected, = ax[0].plot( [], [], \"o\", mfc = 'none')\n",
"\n",
"\n",
"ax[1].scatter( cring.omega, cring.dty,\n",
" c=np.log(cring.sum_intensity) * np.sign( cring.yl ), \n",
" s=4, cmap='RdBu')\n",
" \n",
"\n",
"def find_pair( idx, cf ):\n",
" \"\"\"\n",
" Given a peak idx, locate the Friedel pair in cf\n",
" Given a peak idx, locate the Friedel pair in colfile\n",
" Matches on gvector only (not intensity, etc)\n",
" \"\"\"\n",
" dg2 = (cf.gx[idx] + cf.gx)**2 + (cf.gy[idx] + cf.gy)**2 + (cf.gz[idx] + cf.gz)**2 \n",
" pair = np.argmin( dg2 )\n",
Expand All @@ -264,32 +271,42 @@
"om = np.linspace( cring.omega.min(), cring.omega.max(), 90 )\n",
"fitline, = ax[1].plot( om, np.zeros_like(om), '-')\n",
"\n",
"### GLOBALS\n",
"p1 = None\n",
"p2 = None\n",
"xy = None\n",
"peak_ycalc = None\n",
"### GLOBALS, also y0\n",
"p1 = None # index of first peak in cring\n",
"p2 = None # index of second peak in cring\n",
"xy = None # fitted position\n",
"peak_ycalc = None # dty computed\n",
"\n",
"def fit_pair( cf, i, j ):\n",
" # fits a sine wave\n",
" global xy, p1, p2, peak_ycalc\n",
" global xy, p1, p2, peak_ycalc, y0\n",
" p1 = i\n",
" p2 = j\n",
" # sin/cos omega\n",
" ci = np.cos( np.radians( cf.omega[i] ) )\n",
" cj = np.cos( np.radians( cf.omega[j] ) )\n",
" si = np.sin( np.radians( cf.omega[i] ) )\n",
" sj = np.sin( np.radians( cf.omega[j] ) )\n",
" yi = cf.dty[i]\n",
" yj = cf.dty[j]\n",
" # yi = xp.ci + yp.si\n",
" # yj = xp.cj + yp.sj\n",
" #\n",
" # ImageD11.sinograms.geometry:\n",
" # def dtycalc_sincos(sinomega, cosomega, x, y, y0):\n",
" # ....\n",
" # return y0 - x * sinomega - y * cosomega\n",
" #\n",
" # yi = y0 - xp.si - yp.ci\n",
" # yj = y0 - xp.sj - yp.cj\n",
" #\n",
" yi = cf.dty[i] - y0\n",
" yj = cf.dty[j] - y0\n",
" #\n",
" # Only two peaks, cannot fit x/y/y0\n",
" # Y = R.x\n",
" R = np.array( [ [ ci, si ], [cj, sj] ] )\n",
" xy = np.linalg.inv(R).dot( [ yi, yj ] )\n",
" ycalc = np.cos(np.radians(om))*xy[0]+np.sin(np.radians(om))*xy[1]\n",
" R = np.array( [ [ -si, -ci ], \n",
" [ -sj, -cj] ] )\n",
" xy = x, y = np.linalg.inv(R).dot( [ yi, yj ] )\n",
" ycalc = ImageD11.sinograms.geometry.x_y_y0_omega_to_dty(om, x, y, y0)\n",
" fitline.set_ydata( ycalc )\n",
" peak_ycalc = np.cos(np.radians(cring.omega))*xy[0]+np.sin(np.radians(cring.omega))*xy[1]\n",
" peak_ycalc = ImageD11.sinograms.geometry.x_y_y0_omega_to_dty(cring.omega, x, y, y0)\n",
" \n",
"def onpick( evt ):\n",
" art = evt.artist\n",
Expand All @@ -300,6 +317,7 @@
" selected.set_ydata( [cpk.dty[ind[-1]],cpk.dty[pair]] )\n",
" ax[0].set( title = f'{ind} {pair}' )\n",
" fig.canvas.draw_idle()\n",
"\n",
"ax[1].set(title=\"Pick a peak on the other plot\")\n",
"fig.canvas.mpl_connect( 'pick_event', onpick );"
]
Expand All @@ -312,12 +330,13 @@
},
"outputs": [],
"source": [
"dtyrange = 5.\n",
"print(\"Selected points are\",p1,p2)\n",
"print(\"Omega, dty\")\n",
"print(cpk.omega[p1],cpk.dty[p1])\n",
"print(cpk.omega[p2],cpk.dty[p2])\n",
"f,ax = plt.subplots()\n",
"ax.hist( cring.dty - peak_ycalc, bins = np.arange(-2,2,ds.ystep) );\n",
"ax.hist( cring.dty - peak_ycalc, bins = np.arange(-dtyrange,dtyrange,ds.ystep) );\n",
"ax.set(title='dty error histogram', xlabel='dty', ylabel='frequency');"
]
},
Expand All @@ -330,8 +349,9 @@
"outputs": [],
"source": [
"# now we select the peaks from the sinogram based on that position in space\n",
"dty_cut = 0.4 # from the last plot\n",
"\n",
"ymask = abs(peak_ycalc - cring.dty)<0.4\n",
"ymask = abs(peak_ycalc - cring.dty)<dty_cut\n",
"cgrain = cring.copyrows( ymask )\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot( projection='3d', proj_type='ortho')\n",
Expand Down Expand Up @@ -448,7 +468,8 @@
},
"outputs": [],
"source": [
"pks, dists = find_overlaps_3d( cgrain, (cpk.spot3d_id[p1], cpk.spot3d_id[p2] ), gtol=0.002 )\n",
"gtol = 0.002 # from the last plot\n",
"pks, dists = find_overlaps_3d( cgrain, (cpk.spot3d_id[p1], cpk.spot3d_id[p2] ), gtol=gtol )\n",
"\n",
"fig = plt.figure(constrained_layout=True)\n",
"a = fig.add_subplot(1,1,1, projection='3d', proj_type='ortho')\n",
Expand Down Expand Up @@ -486,7 +507,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Position selectiion on the sinogram"
"## Position selection on the sinogram"
]
},
{
Expand All @@ -498,10 +519,12 @@
"outputs": [],
"source": [
"# make a sinogram and recon\n",
"nbins_angle = min( len( ds.obincens), int(len(ds.ybincens)*1.5) )\n",
"sino, obinedge, ybinedge = np.histogram2d( cring.omega, cring.dty, weights = np.log(cring.sum_intensity),\n",
" bins = (np.linspace(-90, 90, 180), ds.ybinedges) )\n",
" bins = (np.linspace(ds.obinedges.min(), ds.obinedges.max(), nbins_angle), ds.ybinedges) )\n",
"obincen = 0.5*(obinedge[:-1] + obinedge[1:])\n",
"recon = ImageD11.sinograms.roi_iradon.iradon( sino.T, obincen, filter_name='shepp-logan', workers=10)"
"recon = ImageD11.sinograms.roi_iradon.iradon( sino.T, obincen, filter_name='shepp-logan', \n",
" workers=ImageD11.cImageD11.cores_available())"
]
},
{
Expand All @@ -518,10 +541,12 @@
" ybinedge,\n",
" sino.T);\n",
"ax[1].pcolormesh( -ds.ybinedges,ds.ybinedges, recon, vmin=0, vmax = recon.max()/2 )\n",
"location, = ax[1].plot( [],[], 'ow', mfc='none', ms=20)\n",
"ax[1].set(xlim=(ds.ybinedges[-1],ds.ybinedges[0]),\n",
" xlabel='laboratory Y', ylabel='laboratory X')\n",
"ax[0].set(title='sinogram')\n",
"ax[1].set(title='click on the place you want')\n",
"\n",
"om = np.linspace( cring.omega.min(), cring.omega.max(), 90 )\n",
"fitline1, = ax[0].plot( om, np.zeros_like(om), 'w-')\n",
"fitline2, = ax[0].plot( om, np.zeros_like(om), 'w-')\n",
Expand All @@ -538,6 +563,8 @@
" ycalcall = ImageD11.sinograms.geometry.dtycalc( cring.omega, x, y, 0 )\n",
" fitline1.set_ydata( ycalc + 1 )\n",
" fitline2.set_ydata( ycalc - 1 )\n",
" location.set_xdata( [y,] )\n",
" location.set_ydata( [x,] )\n",
" fig.canvas.draw_idle()\n",
" \n",
"fig.canvas.mpl_connect( 'button_press_event', onclick );"
Expand Down
Loading

0 comments on commit bd380f4

Please sign in to comment.