# -*- coding: utf-8 -*- # Make sure to import tomopy 1.01 if it is not the default in your # environment import tomopy import os import numpy as np from tomopy.recon.rotation import write_center from tomopy.recon.algorithm import recon from scipy import misc import glob, fnmatch import h5py import shutil, os ####----------------- Section 1: Input file -- Start ----------------- # Point to the folder containing the "Exp###" hdf files for the given # specimen. Some adjustment of this specimen reading protocol may be # required depending on your environment. In this example, we have # HDF files for specimen numbers 18 to 25 in the folder "Nov_2017". user_top_dir = '/home/projects/APS/xctRecon/rawDataFile/Nov_2017' user_out_dir = '/home/projects/APS/xctRecon/2017_11' for expNum in range(18,25): print 'ExpID = ' + str(expNum) #sam_idx = 1 Exp_idx = expNum #data_top_dir = user_top_dir + 's' + '{:03}'.format(sam_idx) + '*/' data_top_dir = user_top_dir + '/' data_out_dir = user_out_dir + '/' data_dir = glob.glob(os.path.join(data_top_dir,'Exp'+'{:03}'.format(Exp_idx)+'*'))[0] if user_top_dir == user_out_dir: output_dir = data_dir # Output to a folder in the original directory else: # Make an output folder with JUST the experiment number in as the name output_dir = os.path.join(data_out_dir,'Exp' + '{:03}'.format(Exp_idx)) file_data = data_files[0] print file_data BlankAt = 90 if BlankAt == 0: first_missing_end = 87 second_missing_start = 1396 else: missing_start = 657 missing_end = 828 open_start = 739 open_end = 749 missing_width = missing_end - missing_start avg_width = 20 file_name = os.path.join(data_dir, file_data) output_file = output_dir+'/recon_'+file_data.split(".")[-2]+'/recon_'+file_data.split(".")[-2] +'_' ####----------------- Section 1: Input file -- End ----------------- ####----------------- Section 2: Parameter configuration -- Start ----------------- f = h5py.File(file_name,"r") try: arr = f["exchange/data"] except: try: arr = f["entry/data/data"] except: arr = f["MyDataDir/TomoShots"] numProj = (arr.shape)[0] numSlices = (arr.shape)[1] imgWidth = (arr.shape)[2] f.close() print numProj,numSlices,imgWidth offset = 0 # offset from the top slice; for full volume, set to 0 #numRecSlices = 600 # For full volume, comment this out numRecSlices = numSlices # For full volume, uncomment this line chunk_size = 200 if chunk_size > numRecSlices: chunk_size = numRecSlices margin_slices = 30 num_chunk = np.int(numRecSlices/(chunk_size-margin_slices)) + 1 if numRecSlices == chunk_size: num_chunk = 1 z = 20 eng = 60 pxl = 0.65e-4 rat = 1.0e-02 zinger_level = 300 center = [1289 , # proj_3 1291.5 , # proj_4 1290.5 , # proj_5 1293.5 , # proj_6 1293 , # proj_7 1291.5 , # proj_8 1297.5 , # proj_9 1297.5 , # proj_10 1290.5 , # proj_11 1285.5 , # proj_12 1285 , # proj_13 1286 , # proj_14 1287 , # proj_15 1284 , # proj_16 1295 , # proj_17 1286.5 , # proj_18 1287.5 , # proj_19 1286.5 , # proj_20 1286 , # proj_21 1286.5 , # proj_22 1288.5 , # proj_23 1288.5 , # proj_24 1285 , # proj_25 1289 , # proj_26 1290.5 , # proj_27 1289.5 , # proj_28 1286 , # proj_29 1283 , # proj_30 1286.5 , # proj_31 1287 , # proj_32 1296.5 , # proj_33 1296 , # proj_34 1286 , # proj_35 1286 , # proj_36 1285 , # proj_37 1290 , # proj_38 1293.5 , # proj_39 1290 , # proj_40 1290 , # proj_41 1284.5 , # proj_42 # proj_43 # empty 1290 , # proj_44 1286 , # proj_45 1286 , # proj_46 1286 , # proj_47 1287.5 , # proj_48 1290 , # proj_49 1290.5 , # proj_50 1290 , # proj_51 1297.5 , # proj_52 1289 , # proj_53 1285 , # proj_54 1285 , # proj_55 1286 , # proj_56 1288.5 , # proj_57 1290 , # proj_58 1296.5 , # proj_59 1295 , # proj_60 1294.5 , # proj_61 1300.5 , # proj_62 1287 , # proj_63 1290.5 , # proj_64 1299 , # proj_65 1287 , # proj_66 1288.5 , # proj_67 1299 , # proj_68 1286.5 , # proj_69 1286.5 , # proj_70 1286.5 , # proj_71 1294.5 , # proj_72 1287.5 , # proj_73 1298 , # proj_74 1285 , # proj_75 1290 , # proj_76 1299.5 , # proj_77 1289.5 , # proj_78 1289.5 , # proj_79 1300.5 , # proj_80 1291.5 , # proj_81 1292 , # proj_82 1292 , # proj_83 1289 , # proj_84 1290.5 , # proj_85 # proj_86 Mean of zero slice warning, black 1288.5 , # proj_87 1291 , # proj_88 1298 , # proj_89 1288.5 , # proj_90 1291 , # proj_91 1289.5 , # proj_92 1294.5 , # proj_93 1291.5 , # proj_94 1299 , # proj_95 1290.5 , # proj_96 1287.5 , # proj_97 1296.5 , # proj_98 1289 , # proj_99 1289 , # proj_100 1286 , # proj_101 1293 , # proj_102 1288 , # proj_103 1295.5 , # proj_104 1295 , # proj_105 1290.5 , # proj_106 1287.5 , # proj_107 1297.5 , # proj_108 1288 , # proj_109 1284 , # proj_110 1298.5 , # proj_111 1290 , # proj_112 1300 , # proj_113 1288.5 , # proj_114 1284.5 , # proj_115 1290.5 , # proj_116 1293.5 , # proj_117 1281.5 , # proj_118 1287.5 , # proj_119 1293.5 , # proj_120 1292 , # proj_121 1289 , # proj_122 1290 , # proj_123 1282.5 , # proj_124 1287 , # proj_125 1295.5 , # proj_126 1285.5 , # proj_127 1291 , # proj_128 1284 , # proj_129 1296.5 , # proj_130 1293.5 , # proj_131 1299 , # proj_132 1289.5 , # proj_133 1287.5 , # proj_134 1288 , # proj_135 1288 , # proj_136 1291 , # proj_137 1286 , # proj_138 1292 , # proj_139 1287 , # proj_140 1298.5 , # proj_141 1290.5 , # proj_142 1288 , # proj_143 1286] # proj_144 print 'center = ' + str(center[expNum]) ########----------------- Section 4: Full reconstruction -- Start ----------------- for ii in xrange(num_chunk): if ii == 0: sliceStart = offset + ii*chunk_size sliceEnd = offset + (ii+1)*chunk_size else: sliceStart = offset + ii*(chunk_size-margin_slices) sliceEnd = offset + sliceStart + chunk_size if sliceEnd > (offset+numRecSlices): sliceEnd = offset+numRecSlices if sliceEnd > numSlices: sliceEnd = numSlices for ii in range(num_chunk): if ii == 0: sliceStart = offset + ii*chunk_size sliceEnd = offset + (ii+1)*chunk_size else: sliceStart = offset + ii*(chunk_size-margin_slices) sliceEnd = sliceStart + chunk_size if sliceEnd > (offset+numRecSlices): sliceEnd = offset+numRecSlices if sliceEnd > numSlices: sliceEnd = numSlices if (sliceEnd - sliceStart) <= margin_slices: print 'Reconstruction finishes!' break f = h5py.File(file_name,"r") try: arr = f["exchange/data"] except: try: arr = f["MyDataDir/TomoShots"] except: arr = f["entry/data/data"] data = arr[0:numProj,sliceStart:sliceEnd,:] f.close() f = h5py.File(file_name,"r") try: arr = f["exchange/data_white"] except: try: arr = f["MyDataDir/TomoShots"] except: arr = f["entry/data/data"] white = arr[1:9,sliceStart:sliceEnd,:] f.close() f = h5py.File(file_name,"r") try: arr = f["exchange/data_dark"] except: try: arr = f["MyDataDir/TomoShots"] except: arr = f["entry/data/data"] dark = arr[1:9,sliceStart:sliceEnd,:] f.close() data_size = data.shape theta = np.linspace(0,np.pi,num=data_size[0]) print 'data is read' roll_dist = int(2.0*(center[expNum-1] - data_size[2]/2.0)+1) # remove zingers (pixels with abnormal counts) data = tomopy.misc.corr.remove_outlier(data,zinger_level,size=15,axis=0,ncore=numProc) white = tomopy.misc.corr.remove_outlier(white,zinger_level,size=15,axis=0,ncore=numProc) # normalize projection images; for now you need to do below two operations in sequence data = tomopy.prep.normalize.normalize(data,white,dark,ncore=numProc) data = tomopy.prep.normalize.normalize_bg(data,air=10,ncore=numProc) # remove stripes in sinograms data = tomopy.prep.stripe.remove_stripe_fw(data,level=6,wname='sym16',sigma=1,pad=True) # phase retrieval data = tomopy.prep.phase.retrieve_phase(data,pixel_size=pxl,dist=z,energy=eng,alpha=rat,pad=True,ncore=numProc) # more stripe removal in sinograms data = tomopy.prep.stripe.remove_stripe_ti(data,alpha=5,ncore=numProc) # data = tomopy.prep.stripe.remove_stripe_sf(data,size=63,ncore=numProc) data = tomopy.prep.stripe.remove_stripe_sf(data,size=31,ncore=numProc) # data = tomopy.prep.normalize.normalize_bg(data,air=10,ncore=numProc) # data = tomopy.prep.stripe.remove_stripe_fw(data,level=3,wname='sym16',sigma=1.5,pad=True) for ii in range(missing_start,open_start): data[ii,:,:] = np.roll(np.mean(np.roll(data[(ii - missing_width - avg_width):(ii - missing_width),:,:],roll_dist,axis=2),axis=0) + \ np.mean(np.roll(data[(ii + missing_width):(ii + missing_width + avg_width),:,:],roll_dist,axis=2),axis=0),-1*roll_dist)/2.0 for ii in range(open_end,missing_end): data[ii,:,:] = np.roll(np.mean(np.roll(data[(ii - missing_width - avg_width):(ii - missing_width),:,:],roll_dist,axis=2),axis=0) + \ np.mean(np.roll(data[(ii + missing_width):(ii + missing_width + avg_width),:,:],roll_dist,axis=2),axis=0),-1*roll_dist)/2.0 # tomo reconstruction data = tomopy.prep.normalize.minus_log(data,ncore=numProc) # data_recon = recon(data,theta,center=center[expNum-1],algorithm='mlem',num_iter=64,ncore=numProc) data_recon = recon(data,theta,center=center[expNum-1],algorithm='gridrec',ncore=numProc) # mask out a circle around the outside # data_recon = tomopy.misc.corr.circ_mask(data_recon,axis=0,ratio=0.999,val=0.0,ncore=numProc) # Attempt to remove rings from recon # data_recon = tomopy.misc.corr.adjust_range(data_recon, dmin=np.min(data_recon)/0.5, dmax=np.max(data_recon)*0.5) # data_recon = tomopy.misc.corr.remove_outlier(data_recon, dif=100, size=9, axis=0, ncore=numProc, out=data_recon) data_recon = tomopy.misc.corr.remove_ring(data_recon, thresh=300.0, thresh_max=1000.0, thresh_min=-100.0, theta_min=20, rwidth=45, ncore=numProc, out=data_recon) data_recon = tomopy.misc.corr.remove_ring(data_recon, thresh=100.0, thresh_max=2000.0, thresh_min=-200.0, theta_min=10, rwidth=90, ncore=numProc, out=data_recon) data_recon = tomopy.misc.corr.remove_ring(data_recon, thresh=10.0, thresh_max=100.0, thresh_min=-10.0, theta_min=30, rwidth=9, ncore=numProc, out=data_recon) # Subtract a gaussian blur data_blur = scipy.ndimage.filters.gaussian_filter(data_recon,3,1) # data_blur = scipy.ndimage.filters.gaussian_filter(data_blur,5,1) data_recon = data_recon - data_blur # save reconstructions in 8 bit tomopy.io.writer.write_tiff_stack(data_recon[np.int(margin_slices/2):(sliceEnd-sliceStart-np.int(margin_slices/2)),:,:], axis = 0, dtype = 'uint8', vmax = 0.0004, vmin = -0.0006, fname = output_file, start = sliceStart+np.int(margin_slices/2), overwrite = True) print 'reconstruction is saved' print time.asctime() print 'Reconstruction finishes!' # # save reconstructions in 16 bit # tomopy.io.writer.write_tiff_stack(data_recon[np.int(margin_slices/2):(sliceEnd-sliceStart-np.int(margin_slices/2)),:,:], # axis = 0, # fname = output_file, # start = sliceStart+np.int(margin_slices/2), # overwrite = True) ##########----------------- Section 4: Full reconstruction -- End ----------------- # #