In [ ]:
%matplotlib inline
In [ ]:
from __future__ import division
import os
import errno
import glob
try:
import configparser as ConfigParser
except ImportError:
import ConfigParser
import numpy as np
import pylab as plt
import skimage
import skimage.io
import h5py
import astra
import json
import cv2
In [ ]:
def mkdir_p(path):
try:
os.makedirs(path)
except OSError as exc: # Python >2.5
if exc.errno == errno.EEXIST and os.path.isdir(path):
pass
else:
raise
In [ ]:
def log_progress(sequence, every=None, size=None):
from ipywidgets import IntProgress, HTML, VBox
from IPython.display import display
is_iterator = False
if size is None:
try:
size = len(sequence)
except TypeError:
is_iterator = True
if size is not None:
if every is None:
if size <= 200:
every = 1
else:
every = size / 200 # every 0.5%
else:
assert every is not None, 'sequence is iterator, set every'
if is_iterator:
progress = IntProgress(min=0, max=1, value=1)
progress.bar_style = 'info'
else:
progress = IntProgress(min=0, max=size, value=0)
label = HTML()
box = VBox(children=[label, progress])
display(box)
index = 0
try:
for index, record in enumerate(sequence, 1):
if index == 1 or index % every == 0:
if is_iterator:
label.value = '{index} / ?'.format(index=index)
else:
progress.value = index
label.value = u'{index} / {size}'.format(
index=index,
size=size
)
yield record
except:
progress.bar_style = 'danger'
raise
else:
progress.bar_style = 'success'
progress.value = index
label.value = str(index or '?')
In [ ]:
def read_config(config_path):
def as_dict(config):
d = dict(config._sections)
for k in d:
d[k] = dict(config._defaults, **d[k])
d[k].pop('__name__', None)
return d
config = ConfigParser.RawConfigParser()
config.optionxform = str
config.read(config_path)
res = as_dict(config)
return res
In [ ]:
corrections = True
if corrections:
# data_root = '/home/makov/tmp/yaivan/thermal_only/'
data_root = '/home/makov/tmp/yaivan/all_corrections_bh92_rc16/'
else:
data_root = '/home/makov/tmp/yaivan/no_corrections_at_all/'
!ls {data_root}
In [ ]:
log_file = glob.glob(os.path.join(data_root,'*rec.log'))[0]
print(log_file)
In [ ]:
read_config(log_file)
In [ ]:
from scipy.interpolate import interp2d
def load_astra_sinogram(sinogram_number, thermal_correction=False):
# data_root = '/media/makov/buext4/yaivan/MMC_1'
data_root = '/diskmnt/storage0/yaivan/MMC_1'
raw_root = os.path.join(data_root,'Raw')
out_file = os.path.join(data_root, 'raw.h5')
# Without thermal correction
if thermal_correction == False:
with h5py.File(out_file,'r') as h5f:
sinogram = h5f['data'][:,-(sinogram_number+1),:]
else:
# # With thermal correction
# # Loading geometrical correction map generated by NRecon
# corrections_map = np.loadtxt(os.path.join(raw_root,'MMC1_2.82um__TS.crv'),skiprows=2)
# if thermal_correction == True:
# with h5py.File(out_file,'r') as h5f:
# sinogram = np.zeros_like(h5f['data'][:,-(sinogram_number+1),:])
# for ii in log_progress(range(sinogram.shape[0])):
# sinogram[ii] = np.roll(
# h5f['data'][ii,-(sinogram_number+1)-int(corrections_map[ii,2]),:],
# -int(corrections_map[ii,1])
# )
# With thermal correction
# Loading geometrical correction map generated by NRecon
corrections_map = np.loadtxt(os.path.join(raw_root,'MMC1_2.82um__TS.crv'),skiprows=2)
with h5py.File(out_file,'r') as h5f:
sinogram = np.zeros_like(h5f['data'][:,-(sinogram_number+1),:]) # fix sinogram numbering with nrecon
for ii in log_progress(range(sinogram.shape[0])):
c1 = int(np.floor(corrections_map[ii,1])) ## correction factor floor (horizontal)
c2 = int(np.floor(corrections_map[ii,2])) ## correction factor floor (vertical)
proj = h5f['data'][ii, -(sinogram_number+1)-c2:-(sinogram_number+1)-c2+2]
original_x, original_y = np.meshgrid(
np.arange(proj.shape[1]),
np.arange(proj.shape[0])
)
f = interp2d(proj, original_x, original_y)
new_x = np.arange(proj.shape[1])-c1
new_y = np.ones_like(len(new_x))*c2
proj_interp = f(new_x, new_y)
sinogram[ii] = proj_interp
# sinogram[ii] = np.roll(proj[0], -c1)
# select only central 180 deg for reconstruction
# sinogram = sinogram[int(sinogram.shape[0]/2-900):1800-int(sinogram.shape[0]/2-900)]
# sinogram = sinogram[:1800]
sinogram = np.array(sinogram/2.**16)
return sinogram
In [ ]:
sinogram_number = 960
ii=500
out_file = '/diskmnt/storage0/yaivan/MMC_1/raw.h5'
corrections_map = np.loadtxt('/diskmnt/storage0/yaivan/MMC_1/Raw/MMC1_2.82um__TS.crv',skiprows=2)
reference = skimage.io.imread(sinogram_file).astype('float32')[ii]
with h5py.File(out_file,'r') as h5f:
sinogram_no_fix = h5f['data'][-ii,-(sinogram_number+1),:]
c1 = int(np.floor(corrections_map[-ii,1])) ## correction factor floor (horizontal)
c2 = int(np.floor(corrections_map[-ii,2])) ## correction factor floor (vertical)
sinogram_simple_fix = np.roll(h5f['data'][-ii,-(sinogram_number+1)-c2,:],-c1)
frame = h5f['data'][-ii]
original_x, original_y = np.meshgrid(np.arange(frame.shape[1]), np.arange(frame.shape[0]))
f = interp2d(frame, original_x, original_y)
plt.figure(figsize=(15,5))
plt.subplot(121)
plt.plot(reference[1900:2100], label='reference')
plt.plot(sinogram_no_fix[1900:2100], label='no fix')
plt.plot(sinogram_simple_fix[1900:2100], label='simple fix, c1 = {}, c2 ={}'.format(c1,c2))
plt.legend(loc=0)
plt.grid(True)
plt.subplot(122)
plt.plot(reference[1900:2100]-sinogram_no_fix[1900:2100], label='diff no fix')
plt.plot(reference[1900:2100]-sinogram_simple_fix[1900:2100], label='diff simple fix')
plt.legend(loc=0)
plt.grid(True)
plt.show()
In [ ]:
def images_diff(im1, im2):
assert(im1.shape==im2.shape)
rec_diff = np.zeros(shape=(im1.shape[0],im1.shape[1],3), dtype='float32')
im1_t = im1.copy()
im1_t = (im1_t-im1_t.min())/(im1_t.max()-im1_t.min())
im2_t = im2.copy()
im2_t = (im2_t-im2_t.min())/(im2_t.max()-im2_t.min())
# nrecon_rec_t[nrecon_rec_t<0] = 0
diff_rec = im1_t-im2_t
rec_diff[...,0] = diff_rec*(diff_rec>0)
rec_diff[...,1] = -diff_rec*(diff_rec<0)
rec_diff[...,2] = rec_diff[...,1]
return rec_diff
In [ ]:
sinogram_file = glob.glob(os.path.join(data_root,'*__sinoraw_0960.tif'))[0]
sinogram_log_file = glob.glob(os.path.join(data_root,'*__sino0960.tif'))[0]
nrecon_rec_file = glob.glob(os.path.join(data_root,'*.png'))[0]
In [ ]:
sinogram = skimage.io.imread(sinogram_file).astype('float32')
sinogram_log = skimage.io.imread(sinogram_log_file).astype('float32')/2**16
nrecon_rec = skimage.io.imread(nrecon_rec_file).astype('float32')
nrecon_rec = nrecon_rec*(0.52+0.18)-0.18
In [ ]:
astra_sinogram = load_astra_sinogram(960, thermal_correction=corrections)
astra_sinogram = np.flipud(astra_sinogram)
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(sinogram, cmap=plt.cm.viridis)
plt.colorbar(orientation='horizontal');
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(sinogram_log, cmap=plt.cm.viridis)
plt.colorbar(orientation='horizontal');
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(astra_sinogram, cmap=plt.cm.viridis)
plt.colorbar(orientation='horizontal');
In [ ]:
plt.figure(figsize=(15,10))
plt.imshow(images_diff(astra_sinogram,sinogram)*10)
In [ ]:
def cogistogram(im1, im2):
def normalize(im):
t = im.copy()
t = (t-t.min())/(t.min()-t.max())
t = (t*256).astype('uint8')
return t
assert(im1.shape == im2.shape)
t1 = normalize(im1)
t2 = normalize(im2)
res = np.zeros((t1.max(),t2.max()))
for x in range(res.shape[0]):
for y in range(res.shape[1]):
res[t1[x,y],t2[x,y]] +=1
return res
In [ ]:
res = cogistogram(astra_sinogram,sinogram)[:80,:80]
plt.figure(figsize=(15,17))
plt.subplot(121)
plt.imshow(res, cmap=plt.cm.viridis, interpolation='nearest')
plt.subplot(122)
res_non_diag = res - np.diag(np.diag(res))
plt.imshow(res_non_diag, cmap=plt.cm.viridis, interpolation='nearest')
In [ ]:
def build_reconstruction_geomety(detector_size, angles):
# proj_geom = astra.create_proj_geom('parallel', 1.0, detector_size, angles)
#Object to Source (mm) = 56.135
#Camera to Source (mm) = 225.082
# All distances in [pixels]
# pixel_size = 2.82473e-3
pixel_size = 2.82e-3
os_distance = 56.135/pixel_size
ds_distance = 225.082/pixel_size
proj_geom = astra.create_proj_geom('fanflat', ds_distance/os_distance, detector_size, angles,
os_distance, (ds_distance-os_distance))
return proj_geom
In [ ]:
def astra_tomo2d_fanflat_fbp(sinogram, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = sinogram.shape[1]
rec_size = detector_size # size of reconstruction region
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = build_reconstruction_geomety(detector_size, angles)
sinogram_id = astra.data2d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data2d.create('-vol', vol_geom)
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('FBP_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
cfg['option']['ShortScan'] = False
# cfg['option']['MinConstraint'] = 0
# cfg['option']['MaxConstraint'] = 5
# Available algorithms:
# SIRT_CUDA, SART_CUDA, EM_CUDA, FBP_CUDA (see the FBP sample)
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
# Run 150 iterations of the algorithm
astra.algorithm.run(alg_id, 1)
# Get the result
rec = astra.data2d.get(rec_id)
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinogram_id)
astra.clear()
return rec, proj_geom, cfg
def astra_tomo2d_fanflat_sirt(sinogram, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = sinogram.shape[1]
rec_size = detector_size # size of reconstruction region
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = build_reconstruction_geomety(detector_size, angles)
sinogram_id = astra.data2d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data2d.create('-vol', vol_geom)
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('SIRT_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
# cfg['option']['MinConstraint'] = 0
# cfg['option']['MaxConstraint'] = 5
# Available algorithms:
# SIRT_CUDA, SART_CUDA, EM_CUDA, FBP_CUDA (see the FBP sample)
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
# Run 150 iterations of the algorithm
astra.algorithm.run(alg_id, 20)
# Get the result
rec = astra.data2d.get(rec_id)
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinogram_id)
astra.clear()
return rec, proj_geom, cfg
def astra_tomo2d_fanflat_sart(sinogram, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = sinogram.shape[1]
rec_size = detector_size # size of reconstruction region
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = build_reconstruction_geomety(detector_size, angles)
sinogram_id = astra.data2d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data2d.create('-vol', vol_geom)
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('SART_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
# cfg['option']['MinConstraint'] = 0
# cfg['option']['MaxConstraint'] = 5
# Available algorithms:
# SIRT_CUDA, SART_CUDA, EM_CUDA, FBP_CUDA (see the FBP sample)
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
# Run 150 iterations of the algorithm
astra.algorithm.run(alg_id, 2000)
# Get the result
rec = astra.data2d.get(rec_id)
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinogram_id)
astra.clear()
return rec, proj_geom, cfg
# Define the plugin class (has to subclass astra.plugin.base)
# Note that usually, these will be defined in a separate package/module
class SIRTPlugin(astra.plugin.base):
"""Example of an ASTRA plugin class, implementing a simple 2D SIRT algorithm.
Options:
'rel_factor': relaxation factor (optional)
"""
# The astra_name variable defines the name to use to
# call the plugin from ASTRA
astra_name = "SIRT-PLUGIN"
def initialize(self,cfg, rel_factor = 1):
self.W = astra.OpTomo(cfg['ProjectorId'])
self.vid = cfg['ReconstructionDataId']
self.sid = cfg['ProjectionDataId']
self.rel = rel_factor
def run(self, its):
v = astra.data2d.get_shared(self.vid)
s = astra.data2d.get_shared(self.sid)
print(s.shape)
W = self.W
for i in range(its):
v[:] += self.rel*(W.T*(s - (W*v).reshape(s.shape))).reshape(v.shape)/s.size
# from plugin import SIRTPlugin
def astra_tomo2d_fanflat_plugin(sinogram, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = sinogram.shape[1]
rec_size = detector_size # size of reconstruction region
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = build_reconstruction_geomety(detector_size, angles)
proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
sinogram_id = astra.data2d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data2d.create('-vol', vol_geom)
astra.plugin.register(SIRTPlugin)
print(astra.plugin.get_registered())
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('SIRT-PLUGIN')
cfg['ProjectorId'] = proj_id
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
cfg['option']['rel_factor'] = 1.5
# cfg['option']['MinConstraint'] = 0
# cfg['option']['MaxConstraint'] = 5
# Available algorithms:
# SIRT_CUDA, SART_CUDA, EM_CUDA, FBP_CUDA (see the FBP sample)
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
# Run 150 iterations of the algorithm
astra.algorithm.run(alg_id, 10)
# Get the result
rec = astra.data2d.get(rec_id)
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinogram_id)
astra.clear()
return rec, proj_geom, cfg
def create_sinogram(data, angles):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = data.shape[1]
rec_size = detector_size # size of reconstruction region
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = build_reconstruction_geomety(detector_size, angles)
proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
W = astra.OpTomo(proj_id)
P = data
sinogram = W * P
sinogram = sinogram.reshape([len(angles), detector_size])
return np.rot90(sinogram,3)
def get_reconstruction(sinogram, reconstruction_function, min_level=None):
angles = np.arange(sinogram.shape[0])*0.1#-11.493867*2
angles = angles.astype('float64')/180.*np.pi
if min_level is None:
astra_rec, proj_geom, cfg = reconstruction_function(np.flipud(sinogram), angles)
else:
astra_rec, proj_geom, cfg = reconstruction_function(np.flipud(sinogram), angles, min_level)
# print('Projection geometry: {}'.format(proj_geom))
# print('Reconstruction config: {}'.format(cfg))
astra_rec = np.flipud(astra_rec)
return astra_rec
def get_reconstruction_fbp(sinogram):
return get_reconstruction(sinogram, astra_tomo2d_fanflat_fbp)
def get_reconstruction_sirt(sinogram):
return get_reconstruction(sinogram, astra_tomo2d_fanflat_sirt)
def get_reconstruction_sart(sinogram):
return get_reconstruction(sinogram, astra_tomo2d_fanflat_sart)
def get_reconstruction_plugin(sinogram):
return get_reconstruction(sinogram, astra_tomo2d_fanflat_plugin)
In [ ]:
plt.figure(figsize=(10,10))
plt.plot(sinogram_log[0,:])
plt.plot(np.flipud(sinogram_log[1800,:]))
plt.grid()
plt.show()
In [ ]:
for i in log_progress([12]):
if i>=0:
rec = get_reconstruction_fbp(sinogram_log[-1800:,i:])
else:
rec = get_reconstruction_fbp(sinogram_log[-1800:,:i])
# rec[rec<0] = 0
plt.figure(figsize=(10,10))
plt.imshow(rec, cmap=plt.cm.viridis)
plt.title(i)
plt.show()
In [ ]:
plt.figure(figsize=(10,10))
tmp_rec = nrecon_rec.copy()
tmp_rec[tmp_rec<0]=0
plt.imshow(tmp_rec, cmap=plt.cm.viridis)
plt.show()
In [ ]:
rec_diff = np.zeros(shape=(rec.shape[0],rec.shape[1],3), dtype='float32')
rec_t = rec.copy()
rec_t = (rec_t-rec_t.min())/(rec_t.max()-rec_t.min())
# rec_t[rec<0] = 0
nrecon_rec_t = nrecon_rec.copy()
nrecon_rec_t = (nrecon_rec_t-nrecon_rec_t.min())/(nrecon_rec_t.max()-nrecon_rec_t.min())
# nrecon_rec_t[nrecon_rec_t<0] = 0
diff_rec = rec_t-nrecon_rec_t
rec_diff[...,0] = diff_rec*(diff_rec>0)
rec_diff[...,1] = -diff_rec*(diff_rec<0)
rec_diff[...,2] = rec_diff[...,1]
In [ ]:
plt.figure(figsize=(10,10))
plt.imshow(rec_diff*8)
# plt.colorbar(orientation='horizontal');
In [ ]:
# for nrecon_folder in log_progress(nrecon_folders):
# data_file = os.path.join(nrecon_folder, 'MMC1_2.82um__sino0960.tif')
# logging.info('Sinogram file: {}'.format(data_file))
# sinogram = plt.imread(data_file)
# logging.info('Sinogram angles, length: {}'.format(sinogram.shape))
# nrecon_config_file = os.path.join(nrecon_folder, 'MMC1_2.82um__rec.log')
# nrecon_config = read_config(nrecon_config_file)
# rec_sart = get_reconstruction_sart(sinogram)
# output_folder = os.path.join(astra_root_folder, nrecon_folder[len(nrecon_root_folder)+1:])
# mkdir_p(output_folder)
# astra_sart_file = os.path.join(output_folder, 'MMC1_2.82um__rec0960_astra_sart.png')
# logging.info('Output file: {}'.format(astra_sart_file))
# plt.imsave(astra_sart_file, rec_sart, cmap = plt.cm.gray)
# data_config = os.path.join(output_folder, 'MMC1_2.82um__rec.log')
# logging.info('Output config file: {}'.format(data_config))
# config = ConfigParser.RawConfigParser()
# config.optionxform = str
# config.add_section('Reconstruction')
# config.set('Reconstruction', 'Minimum for CS to Image Conversion', rec_sart.min())
# config.set('Reconstruction', 'Maximum for CS to Image Conversion', rec_sart.max())
# bh = nrecon_config['Reconstruction']['Beam Hardening Correction (%)']
# rc = nrecon_config['Reconstruction']['Ring Artifact Correction']
# config.set('Reconstruction', 'Beam Hardening Correction (%)', bh)
# config.set('Reconstruction', 'Ring Artifact Correction', rc)
# with open(data_config, 'wb') as configfile:
# config.write(configfile)
# # # sinogram = sinogram[-1800:]
# # nrecon_rec_file = os.path.join(nrecon_folder,'MMC1_2.82um__rec0960.png')
# # nrecon_rec = plt.imread(nrecon_rec_file)[...,0]
In [ ]:
nrecon_folder = [d for d in nrecon_folders if 'bh_92_rc_20' in d][0]
nrecon_rec = plt.imread(os.path.join(nrecon_folder, 'MMC1_2.82um__rec0960.png'))[...,0]
nrecon_rec = nrecon_rec*(0.52+0.18)-0.18
print nrecon_folder
In [ ]:
import matplotlib
font = {'size' : 18}
matplotlib.rc('font', **font)
plt.figure(figsize=(10,12))
plt.imshow(nrecon_rec, cmap=plt.cm.gray)
plt.colorbar()
# plt.colorbar(orientation='horizontal')
In [ ]:
# buzmakov
from numba import jit
import logging
from scipy import ndimage
import skimage.io
def calculate_background(data, zeros_mask):
labeled_mask, num_features = ndimage.measurements.label(zeros_mask)
logging.info('Found regions: {}'.format(num_features-1))
sigma = []
for nf in range(num_features):
if nf == 0 :
continue
data_constant = data[labeled_mask==nf]
s = np.std(data_constant)
sigma.append(s)
logging.info('STD for regions: {}'.format(sigma))
std = np.mean(sigma)
logging.info('Mean STD for regions: {}'.format(std))
mean_value = data.mean()
logging.info('Mean reconstructed value for all data: {}'.format(mean_value))
res = std/mean_value
logging.info('Normalized STD: {}'.format(res))
return res
#ingacheva
from scipy import misc
from scipy import ndimage
@jit
def base_value(distance_transform, original):
delta = distance_transform.max() * 0.9
threshold = distance_transform.max() - delta
xx = 1 * np.logical_and(distance_transform >= threshold, distance_transform > 0)
summ = xx.sum()
base_v = original[distance_transform >= threshold].sum() / summ
return base_v
@jit
def weighted_variance(mask, distance_transform, original):
base_v = base_value(distance_transform, original)
weight = np.zeros_like(mask)
weight[mask > 0.0] = np.sqrt(distance_transform[mask > 0.0])
threshold = distance_transform.max() * 0.8
xx = 1 * np.logical_and(distance_transform <= threshold, distance_transform > 0.0)
orig = np.zeros_like(mask)
orig = orig.astype('float64')
orig[xx > 0.0] = original[xx > 0.0] - base_v
res = weight[xx > 0.0] * np.power(orig[xx > 0.0], 2)
result = np.sqrt(res.sum() / weight[xx > 0.0].sum())
return result
@jit
def calck_square(mask, distance_transform, original):
base_v = base_value(distance_transform, original)
iter_max = int(distance_transform.max())
sq = 0
for i in range(1, iter_max):
value = original[(distance_transform>=i)*(distance_transform< i+1)].sum()
sq += np.abs(base_v - value)
sq = sq / (base_v * iter_max)
return sq
@jit
def calculate_cupping(original, mask):
mask[mask > 0.0] = 1.0
labeled, nr_objects = ndimage.label(mask > 0.0)
# ndimage.find_objects(labeled)
logging.info('Number of objects is {}'.format(nr_objects))
result = 0
square = 0
for i in range(1, nr_objects+1):
mask[mask > 0.0] = 0.0
mask[labeled == i] = 1.0
sx = mask.sum(axis=0)
sxx = np.argwhere(sx>0)
x_min = sxx.min()
x_max = sxx.max()
sy = mask.sum(axis=1)
syy = np.argwhere(sy>0)
y_min = syy.min()
y_max = syy.max()
dist = ndimage.distance_transform_edt(mask[y_min:y_max, x_min:x_max])
#res = weighted_variance(mask, dist, original)
#result += res
#data['weighted_variance'] = res
sq = calck_square(mask[y_min:y_max, x_min:x_max],
dist,
original[y_min:y_max, x_min:x_max])
square += sq
logging.info("square {} of the object {}".format(sq, i))
#result = result / nr_objects
#print 'mean weighted variance ', result
square = square / nr_objects
return square
mask_background = skimage.io.imread(
'/diskmnt/a/makov/yaivan/MMC_1/_tmp/binary_masks/MMC1_2.82um__rec0960_MASK_ZEROS_CONERS.png')[...,0]
mask_cup = skimage.io.imread(
'/diskmnt/a/makov/yaivan/MMC_1/_tmp/binary_masks/MMC1_2.82um__rec0960_Mask_objects.png')[...,0]
In [ ]:
plt.figure(figsize=(10,5))
plt.subplot(121)
plt.imshow(mask_background, cmap=plt.cm.viridis)
plt.subplot(122)
plt.imshow(mask_cup, cmap=plt.cm.viridis)
In [ ]:
data_file = os.path.join(nrecon_folder, 'MMC1_2.82um__sino0960.tif')
sinogram = plt.imread(data_file)
In [ ]:
plt.figure(figsize=(12,10))
plt.imshow(sinogram, cmap=plt.cm.viridis)
plt.colorbar(orientation='horizontal')
In [ ]:
rois = []
rois.append(np.ix_(np.r_[2100:2350],np.r_[1700:2000]))
rois.append(np.ix_(np.r_[2750:3200],np.r_[2350:2900]))
rois.append(np.ix_(np.r_[2300:2380],np.r_[1200:1300]))
rois.append(np.ix_(np.r_[2600:2750],np.r_[1100:1300]))
rois.append(np.ix_(np.r_[2400:2900],np.r_[600:1000]))
rois.append(np.ix_(np.r_[2700:3500],np.r_[1300:2100]))
rois.append(np.ix_(np.r_[1450:2000],np.r_[600:1150]))
rois.append(np.ix_(np.r_[1750:2050],np.r_[2500:2930]))
In [ ]:
def astra_my(sinogram, angles, min_level=0):
angles = angles.astype('float64') # hack for astra stability, may be removed in future releases
detector_size = sinogram.shape[1]
rec_size = detector_size # size of reconstruction region
vol_geom = astra.create_vol_geom(rec_size, rec_size)
proj_geom = build_reconstruction_geomety(detector_size, angles)
sinogram_id = astra.data2d.create('-sino', proj_geom, data=sinogram)
# Create a data object for the reconstruction
rec_id = astra.data2d.create('-vol', vol_geom)
# Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra.astra_dict('SART_CUDA')
cfg['ReconstructionDataId'] = rec_id
cfg['ProjectionDataId'] = sinogram_id
cfg['option'] = {}
cfg['option']['MinConstraint'] = min_level
# cfg['option']['MaxConstraint'] = 5
# Available algorithms:
# SIRT_CUDA, SART_CUDA, EM_CUDA, FBP_CUDA (see the FBP sample)
# Create the algorithm object from the configuration structure
alg_id = astra.algorithm.create(cfg)
# Run 150 iterations of the algorithm
astra.algorithm.run(alg_id, 10000)
# Get the result
rec = astra.data2d.get(rec_id)
# Clean up. Note that GPU memory is tied up in the algorithm object,
# and main RAM in the data objects.
astra.algorithm.delete(alg_id)
astra.data2d.delete(rec_id)
astra.data2d.delete(sinogram_id)
astra.clear()
return rec, proj_geom, cfg
def get_reconstruction_my(sinogram, min_level):
return get_reconstruction(sinogram, astra_my, min_level)
In [ ]:
# rec_my = rec_fbp
# plt.figure(figsize=(15,15))
# plt.imshow(rec_my/(rec_my.max()-rec_my.min())-nrecon_rec/(nrecon_rec.max()-nrecon_rec.min()), cmap=plt.cm.gray)
# plt.title('My')
# plt.show()
logging.getLogger('').setLevel(logging.ERROR)
rec_fbp = get_reconstruction_fbp(sinogram[-1800:])
for min_level in log_progress(np.arange(-3,-1,1)):
print 'Min_level = {}'.format(min_level)
rec_my = get_reconstruction_my(sinogram, min_level)
artifact_my_bg = calculate_background(rec_my, mask_background)
artifact_my_cup = calculate_cupping(rec_my, mask_cup)
print 'My: bg:{}, cup:{}'.format(artifact_my_bg, artifact_my_cup)
artifact_fbp_bg = calculate_background(rec_fbp, mask_background)
artifact_fbp_cup = calculate_cupping(rec_fbp, mask_cup)
print 'FBP: bg:{}, cup:{}'.format(artifact_fbp_bg, artifact_fbp_cup)
artifact_nrecon_bg = calculate_background(nrecon_rec, mask_background)
artifact_nrecon_cup = calculate_cupping(nrecon_rec, mask_cup)
print 'NRecon: bg:{}, cup:{}'.format(artifact_nrecon_bg, artifact_nrecon_cup)
for roi in rois:
d_my = rec_my[roi]
d_fbp = rec_fbp[roi]
d_nrecon = nrecon_rec[roi]
plt.figure(figsize=(15,15))
plt.subplot(221)
plt.imshow(d_my, cmap=plt.cm.gray, vmin=0)
plt.title('My')
plt.subplot(222)
plt.imshow(d_fbp, cmap=plt.cm.gray, vmin=0)
plt.title('FBP')
plt.subplot(223)
plt.imshow(d_nrecon, cmap=plt.cm.gray, vmin=0)
plt.title('NRecon')
plt.subplot(224)
pos = int(d_my.shape[0]/2)
d_my_1 = d_my[pos]
d_my_1 = d_my_1/ (d_my_1.max()-d_my_1.min())
d_fbp_1 = d_fbp[pos]
d_fbp_1 = d_fbp_1/ (d_fbp_1.max()-d_fbp_1.min())
d_nrecon_1 = d_nrecon[pos]
d_nrecon_1 = d_nrecon_1/ (d_nrecon_1.max()-d_nrecon_1.min())
plt.plot(d_my_1, label='my')
plt.plot(d_fbp_1, label='FBP')
plt.plot(d_nrecon_1, label='NRecon')
plt.grid()
plt.legend(loc=0)
# plt.colorbar(orientation='horizontal')
plt.show()
plt.figure(figsize=(15,15))
plt.imshow(rec_my, cmap=plt.cm.gray)
plt.title('My')
plt.show()
plt.figure(figsize=(15,15))
plt.imshow(nrecon_rec, cmap=plt.cm.gray)
plt.title('NRecon')
plt.show()
In [ ]:
np.eye?
In [ ]:
import scipy.interpolate
In [ ]:
scipy.interpolate.interp2d?
In [ ]:
np.trunc?
In [ ]:
a = np.array([-1.7, -1.5, -0.2, 0.2, 1.5, 1.7, 2.0])
np.floor(a)
In [ ]:
cv2.