# -*- coding: utf-8 -*- """ Created on Tue Jul 21 23:01:40 2015 @author: xhxiao """ import tomopy import numpy as np from tomopy.recon.rotation import write_center from tomopy.recon.algorithm import recon import scipy import glob, fnmatch import h5py import shutil, os, time ####----------------- Section 1: Input file -- Staragt ----------------- print time.asctime() user_top_dir = '/home/projects/APS/xctRecon/rawDataFiles/June2017' user_out_dir = '/home/projects/APS/xctRecon/2017_06/new' centering_dir = '/home/projects/APS/xctRecon/2017_06/reconScripts/center' # Get the number of processors assigned by scheduler (per node) numProc = os.environ['PBS_NUM_PPN'] print 'numProc = ' + str(numProc) #reconList = {25} reconList = range(151,201) for expNum in reconList: print 'ExpID = ' + str(expNum) 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(expNum)+'*'))[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 (three digits) output_dir = os.path.join(data_out_dir,'Exp' + '{:03}'.format(expNum)) data_files = np.sort(fnmatch.filter(os.listdir(data_dir),'*.hdf')) if len(data_files) == 0: data_files = np.sort(fnmatch.filter(os.listdir(data_dir),'*.h')) file_data = data_files[0] file_flat = data_files[0] file_dark = data_files[0] print file_data print file_flat print file_dark file_name = os.path.join(data_dir, file_data) flat_name = os.path.join(data_dir, file_flat) dark_name = os.path.join(data_dir, file_dark) output_file = os.path.join(output_dir,'8bit','Exp'+'{:03}'.format(expNum)) # New for missing frames recon: 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 ###----------------- 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() offset = 40 # offset from the top slice; for full volume, set to 0 #numRecSlices = 1000 # 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 ###----------------- Section 2: Parameter configuration -- End ----------------- ###----------------- Section 2: Center Selectiong -- Start ----------------- # Centering: see centering_2017_06.py and centerList.txt. Copy the center list that results here. center = [1276.0, # Exp001 1289.0, # Exp002 - low contrast 1283.0, # Exp003 - low contrast 1283.5, # Exp004 - " 1283.0, # Exp005 - " 1285.0, # Exp006 - " 1285.0, # Exp007 - " 1282.0, # Exp008 - " 1283.5, # Exp009 - " 1283.5, # Exp010 - " 1278.0, # Exp011 1279.0, # Exp012 1278.5, # Exp013 1278.5, # Exp014 - low contrast 1279.0, # Exp015 1280.5, # Exp016 1278.5, # Exp017 1256.5, # Exp018 1285.5, # Exp019 1283.0, # Exp020 1284.5, # Exp021 1285.0, # Exp022 1283.5, # Exp023 1284.5, # Exp024 1284.0, # Exp025 1280.0, # Exp026 1279.5, # Exp027 1278.5, # Exp028 1279.0, # Exp029 1280.5, # Exp030 1281.0, # Exp031 1285.0, # Exp032 1283.5, # Exp033 1285.5, # Exp034 1286.0, # Exp035 1285.0, # Exp036 1285.5, # Exp037 1283.0, # Exp038 1282.0, # Exp039 1281.0, # Exp040 1287.0, # Exp041 1288.5, # Exp042 1287.0, # Exp043 1284.0, # Exp044 1282.0, # Exp045 1281.5, # Exp046 1282.5, # Exp047 1281.5, # Exp048 1288.5, # Exp049 1278.0, # Exp050 1285.0, # Exp101 1285.0, # Exp102 1281.0, # Exp103 1285.0, # Exp104 1275.0, # Exp105 1278.5, # Exp106 1279.5, # Exp107 1281.0, # Exp108 1281.5, # Exp109 1285.0, # Exp110 1282.5, # Exp111 1286.0, # Exp112 1283.0, # Exp113 1285.0, # Exp114 1275.0, # Exp115 1282.5, # Exp116 1286.5, # Exp117 1283.5, # Exp118 1286.0, # Exp119 1285.0, # Exp120 1285.5, # Exp121 1284.5, # Exp122 1280.0, # Exp123 1285.0, # Exp124 1278.5, # Exp125 1277.0, # Exp126 1281.5, # Exp127 1283.5, # Exp128 1281.5, # Exp129 1277.5, # Exp130 1281.5, # Exp131 1279.0, # Exp132 1279.5, # Exp133 1279.5, # Exp134 1286.0, # Exp135 1275.0, # Exp136 1282.0, # Exp137 1285.5, # Exp138 1284.5, # Exp139 1287.5, # Exp140 1285.0, # Exp141 1285.0, # Exp142 1285.0, # Exp143 1284.5, # Exp144 1281.5, # Exp145 1279.0, # Exp146 1280.5, # Exp147 1286.0, # Exp148 1281.0, # Exp149 1281.5, # Exp150 1284.0, # Exp151 1283.0, # Exp152 1283.5, # Exp153 1286.0, # Exp154 1278.5, # Exp155 1283.5, # Exp156 1283.0, # Exp157 1282.5, # Exp158 1280.0, # Exp159 1281.0, # Exp160 1273.5, # Exp161 1288.0, # Exp162 1283.0, # Exp163 1283.0, # Exp164 1284.0, # Exp165 1282.5, # Exp166 1280.0, # Exp167 1281.5, # Exp168 1279.0, # Exp169 1282.0, # Exp170 1282.5, # Exp171 1282.5, # Exp172 1277.5, # Exp173 1277.0, # Exp174 1281.5, # Exp175 1286.0, # Exp176 1285.5, # Exp177 1282.5, # Exp178 1279.0, # Exp179 1280.0, # Exp180 1274.5, # Exp181 1286.0, # Exp182 1285.0, # Exp183 1285.5, # Exp184 1283.5, # Exp185 1284.5, # Exp186 1285.5, # Exp187 1283.5, # Exp188 1285.5, # Exp189 1286.0, # Exp190 1286.0, # Exp191 1285.0, # Exp192 1286.0, # Exp193 1279.5, # Exp194 1285.0, # Exp195 1278.5, # Exp196 1283.0, # Exp197 1285.0, # Exp198 1286.0, # Exp199 1285.0, # Exp200 1287.0, # Exp201 1285.5, # Exp202 1283.5, # Exp203 1284.0, # Exp204 1285.0, # Exp205 1284.0, # Exp206 1279.0, # Exp207 1283.0, # Exp208 1281.5, # Exp209 1279.5, # Exp210 1280.0, # Exp211 1274.0, # Exp212 1273.5, # Exp213 1285.0] # Exp214 # Make sure the correct center is used: curCenter=center[expNum-1] print 'center = ' + str(curCenter) ###----------------- Section 2: Center Selectiong -- End ----------------- ###----------------- 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 = 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' print 'chunk ' + str(ii) roll_dist = int(2.0*(curCenter - 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) # phase retrieval data = tomopy.prep.phase.retrieve_phase(data,pixel_size=pxl,dist=z,energy=eng,alpha=rat,pad=True,ncore=numProc) # remove stripes in sinograms data = tomopy.prep.stripe.remove_stripe_ti(data,alpha=5,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) 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=curCenter,algorithm='gridrec',ncore=numProc) # Attempt to remove rings from recon data_recon = tomopy.misc.corr.adjust_range(data_recon, dmin=np.min(data_recon)/0.2, dmax=np.max(data_recon)*0.2) 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=500.0, thresh_max=300.0, thresh_min=-100.0, theta_min=20, rwidth=45, ncore=numProc, out=data_recon) # data_recon = tomopy.misc.corr.remove_outlier(data_recon, dif=100, size=5, axis=0, ncore=numProc, out=data_recon) # Subtract a gaussian blur data_blur = scipy.ndimage.filters.gaussian_filter(data_recon,3,1) data_recon = data_recon - data_blur # data_recon = tomopy.misc.corr.remove_outlier(data_recon, dif=100, size=3, axis=1, ncore=numProc, out=data_recon) data_recon = tomopy.misc.corr.remove_outlier(data_recon, dif=100, size=3, axis=2, ncore=numProc, out=data_recon) # data_recon = tomopy.misc.corr.gaussian_filter(data_recon, sigma=3, order=0, axis=0, ncore=numProc) # data_recon = tomopy.misc.corr.gaussian_filter(data_recon, sigma=3, order=0, axis=1, ncore=numProc) # data_recon = tomopy.misc.corr.gaussian_filter(data_recon, sigma=3, order=0, axis=2, ncore=numProc) # data_recon = tomopy.misc.corr.median_filter(data_recon, size=5, axis=0, ncore=numProc) # data_recon = tomopy.misc.corr.median_filter(data_recon, size=5, axis=1, ncore=numProc) # data_recon = tomopy.misc.corr.median_filter(data_recon, size=5, axis=2, ncore=numProc) # save reconstructions 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.001, vmin = -0.0011, fname = output_file, start = sliceStart+np.int(margin_slices/2), overwrite = True) #########----------------- Section 4: Full reconstruction -- End ----------------- # print 'reconstruction is saved' print time.asctime() print 'Reconstruction finishes!'