In [4]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.feature import register_translation
from skimage.feature.register_translation import _upsampled_dft
from scipy.ndimage import fourier_shift
%matplotlib inline
In [7]:
MAPI_01Path='./data/MAPI_01.txt'
MAPI_00Path='./data/MAPI_00.txt'
#load in txt files corresponding to height
MAPI01Height=np.loadtxt(MAPI_01Path, delimiter="\t", dtype='float')
MAPI00Height=np.loadtxt(MAPI_00Path, delimiter="\t", dtype='float')
ImageProduct=np.fft.fft2(MAPI00Height)*np.fft.fft2(MAPI01Height).conj()
CrossCorrelation=np.fft.fftshift(np.fft.ifft2(ImageProduct))
fig=plt.figure(figsize=(12,10))
MAPI01_plt=fig.add_subplot(131)
MAPI00_plt=fig.add_subplot(132)
CrossCorrelation_plt=fig.add_subplot(133)
MAPI01_plt.set_title('MAPI01', fontsize=14)
MAPI00_plt.set_title('MAPI00', fontsize=14)
#what does this cross correlation show? what does the dark spot correspond to?
CrossCorrelation_plt.set_title('Cross Correlation', fontsize=14)
MAPI01_plt.imshow(MAPI01Height,cmap='viridis')
MAPI00_plt.imshow(MAPI00Height,cmap='viridis')
CrossCorrelation_plt.imshow(CrossCorrelation.real)
plt.show()
print("Detected pixel offset (y, x): {}".format(shift))
#5x5 um images of roughly the same area
#256x256 image
In [6]:
#01 is offset from 00
shift,error,diffphase=register_translation(MAPI00Height, MAPI01Height)
#what about rotations?
#what is the difference between the pixel precision and subpixel precision and why would you use one over the other?
#how do you practically shift one image over and crop both to a common area so they are the same dimension?
In [8]:
#offset the image to match 00:
MAPI01HeightRolled=np.roll(MAPI01Height,(1,7), axis=(1,0))
#also try slicing, but not sure if i am slicing the correct region.
MAPI01HeightSlice=MAPI01Height[7:,1:]
#display with offset and sliced image
fig=plt.figure(figsize=(15,10))
MAPI01_plt=fig.add_subplot(142)
MAPI00_plt=fig.add_subplot(141)
MAPI01HeightRolled_plt=fig.add_subplot(143)
MAPI01HeightSlice_plt=fig.add_subplot(144)
MAPI01_plt.set_title('MAPI01', fontsize=14)
MAPI00_plt.set_title('MAPI00', fontsize=14)
#what does this cross correlation show? what does the dark spot correspond to?
MAPI01HeightRolled_plt.set_title('MAPI01 Rolled', fontsize=14)
MAPI01HeightSlice_plt.set_title('MAPI01 Slice', fontsize=14)
MAPI01_plt.imshow(MAPI01Height,cmap='viridis')
MAPI00_plt.imshow(MAPI00Height,cmap='viridis')
MAPI01HeightRolled_plt.imshow(MAPI01HeightRolled, cmap='viridis')
MAPI01HeightSlice_plt.imshow(MAPI01HeightSlice, cmap='viridis')
plt.show()
MAPI01HeightSlice.shape
Out[8]:
In [9]:
ImageProduct2=np.fft.fft2(MAPI00Height)*np.fft.fft2(MAPI01HeightRolled).conj()
CrossCorrelation2=np.fft.fftshift(np.fft.ifft2(ImageProduct))
fig=plt.figure(figsize=(12,10))
MAPI01Rolled_plt=fig.add_subplot(131)
MAPI00_plt=fig.add_subplot(132)
CrossCorrelation_plt2=fig.add_subplot(133)
MAPI01Rolled_plt.set_title('MAPI01 Rolled', fontsize=14)
MAPI00_plt.set_title('MAPI00', fontsize=14)
CrossCorrelation_plt2.set_title('Cross Correlation', fontsize=14)
MAPI01Rolled_plt.imshow(MAPI01HeightRolled,cmap='viridis')
MAPI00_plt.imshow(MAPI00Height,cmap='viridis')
CrossCorrelation_plt2.imshow(CrossCorrelation2.real)
plt.show()
shift,error,diffphase=register_translation(MAPI00Height, MAPI01HeightRolled)
print("Detected pixel offset (y, x): {}".format(shift))
#hmm, expect to be 0 , but it's not cropped so maybe that's why.
In [10]:
# conda install dipy -c conda-forge
import numpy as np
from dipy.viz import regtools
from dipy.data import fetch_stanford_hardi, read_stanford_hardi
from dipy.data.fetcher import fetch_syn_data, read_syn_data
from dipy.align.imaffine import (transform_centers_of_mass,
AffineMap,
MutualInformationMetric,
AffineRegistration)
from dipy.align.transforms import (TranslationTransform3D,
RigidTransform3D,
AffineTransform3D)
from dipy.align.transforms import (TranslationTransform2D,
RigidTransform2D,
AffineTransform2D)
In [14]:
nbins = 32
sampling_prop = None
metric = MutualInformationMetric(nbins, sampling_prop)
level_iters = [10000, 1000, 100]
sigmas = [3.0, 1.0, 0.0]
factors = [4, 2, 1]
affreg = AffineRegistration(metric=metric,
level_iters=level_iters,
sigmas=sigmas,
factors=factors)
In [58]:
transform = AffineTransform2D()
params0 = None
static=MAPI00Height
moving=MAPI01Height
affine = affreg.optimize(static, moving, transform, params0)
In [59]:
transformed=affine.transform(moving)
In [60]:
plt.matshow(transformed, cmap='viridis')
plt.matshow(static,cmap='viridis')
plt.matshow(transformed-static, cmap='viridis')
Out[60]:
In [61]:
from dipy.align.imwarp import SymmetricDiffeomorphicRegistration
from dipy.align.metrics import SSDMetric, CCMetric, EMMetric
from dipy.viz import regtools
In [127]:
#take center 128x128
static=MAPI00Height128=np.copy(MAPI00Height[63:191,63:191])
moving=MAPI01Height128=np.copy(MAPI01Height[63:191,63:191])
#static=MAPI00Height
#moving=MAPI01Height
In [128]:
dim = static.ndim
metric = CCMetric(dim)
In [129]:
level_iters = [200, 100, 50, 25]
sdr = SymmetricDiffeomorphicRegistration(metric, level_iters, inv_iter = 50)
In [130]:
mapping = sdr.optimize(static, moving)
In [131]:
regtools.plot_2d_diffeomorphic_map(mapping, 10, 'diffeomorphic_map.png')
Out[131]:
In [132]:
warped_moving = mapping.transform(moving, 'linear')
regtools.overlay_images(static, warped_moving, 'Static','Overlay','Warped moving',
'direct_warp_result.png')
Out[132]:
In [133]:
plt.matshow(warped_moving-static, cmap='viridis')
Out[133]:
In [ ]:
mapping.forward.shape
In [ ]:
plt.matshow(mapping.forward[:, :, 0], cmap='viridis')
plt.matshow(mapping.forward[:, :, 1], cmap='viridis')
In [ ]: