Merge pull request #13 from LaurentRDC/develop

Develop

RELEASE.rst

@@ -1,7 +1,7 @@
 Steps to Release scikit-ued
 ===========================
 
-These are the steps to take to create a release of the module ``skued``:
+These are the steps to take to create a release of the package ``skued``:
 
 1. Switch to the release branch
 

docs/source/api.rst

@@ -10,7 +10,9 @@ Parallel Utilities
 
 Plot Utilities
 ==============
-.. automodule:: skued.plot_utils
+.. autofunction:: skued.plot_utils.spectrum_colors
+
+.. autofunction:: skued.plot_utils.rgb_sweep
 
 Array Utilities
 ===============

docs/source/conf.py

@@ -16,6 +16,7 @@
 #
 import os
 import sys
+
 sys.path.insert(0, os.path.abspath('../..'))
 
 import skued

docs/source/tutorials/baseline.rst

@@ -25,7 +25,7 @@ example polycrystalline vanadium dioxide pattern:
 	import matplotlib.pyplot as plt
 	import numpy as np
 
-	s, intensity = np.load('data\\powder.npy')
+	s, intensity = np.load('powder.npy')
 	fig = plt.figure()
 	ax = fig.add_subplot(111)
 	ax.plot(s, intensity, 'k')
@@ -50,7 +50,7 @@ substrates, as well as inelastic scattering effects::
 	import numpy as np
 	from skued import gaussian
 
-	s, intensity = np.load('data\\powder.npy')
+	s, intensity = np.load('powder.npy')
 
 	# Double exponential inelastic background and substrate effects
 	background = 75 * np.exp(-7 * s) + 55 * np.exp(-2 * s)
@@ -78,7 +78,7 @@ Scikit-ued offers two ways of removing the background.
 Iterative Baseline Determination using the Discrete Wavelet Transform
 =====================================================================
 
-The prodecude and rational for the :code:`baseline_dwt` routine is described in detail in:
+The procedure and rational for the :code:`baseline_dwt` routine is described in detail in:
 
 Galloway et al. 'An Iterative Algorithm for Background Removal in Spectroscopy by Wavelet 
 Transforms', Applied Spectroscopy pp. 1370 - 1376, September 2009.
@@ -89,7 +89,7 @@ Here is a usage example for the data presented above::
 	from skued import gaussian
 	from skued.baseline import baseline_dwt
 
-	s, intensity = np.load('data\\powder.npy')
+	s, intensity = np.load('powder.npy')
 
 	# Double exponential inelastic background and substrate effects
 	diffuse = 75 * np.exp(-7 * s) + 55 * np.exp(-2 * s)
@@ -105,7 +105,7 @@ Here is a usage example for the data presented above::
 	from skued import gaussian, spectrum_colors
 	from skued.baseline import baseline_dwt
 
-	s, intensity = np.load('data\\powder.npy')
+	s, intensity = np.load('powder.npy')
 
 	# Double exponential inelastic background and substrate effects
 	diffuse = 75 * np.exp(-7 * s) + 55 * np.exp(-2 * s)
@@ -155,7 +155,7 @@ Here is a usage example for the data presented above::
 	from skued import gaussian
 	from skued.baseline import baseline_dt
 
-	s, intensity = np.load('data\\powder.npy')
+	s, intensity = np.load('powder.npy')
 
 	# Double exponential inelastic background and substrate effects
 	diffuse = 75 * np.exp(-7 * s) + 55 * np.exp(-2 * s)
@@ -171,7 +171,7 @@ Here is a usage example for the data presented above::
 	from skued import gaussian, spectrum_colors
 	from skued.baseline import baseline_dt
 
-	s, intensity = np.load('data\\powder.npy')
+	s, intensity = np.load('powder.npy')
 
 	# Double exponential inelastic background and substrate effects
 	diffuse = 75 * np.exp(-7 * s) + 55 * np.exp(-2 * s)

docs/source/tutorials/image.rst

@@ -13,6 +13,7 @@ Contents
 ========
 
 * :ref:`streaming`
+* :ref:`alignment`
 * :ref:`powder`
 
 .. _streaming:
@@ -105,43 +106,127 @@ Here is a recipe for it::
 	    errors = isem(stream2)
 	    yield from zip(averages, errors)
 
-.. _powder:
+.. _alignment:
 
-Image analysis on polycrystalline diffraction patterns
-======================================================
+Diffraction pattern alignment
+=============================
 
-Center-finding
---------------
-Polycrystalline diffraction patterns display concentric rings, and finding
-the center of those concentric rings is important.
+Diffraction patterns can drift over a period of a few minutes, and for reliable data synthesis
+it is important to align patterns to a reference.
+
+The procedure of detecting, or registering, the translation between two similar images is usually
+done by measuring the cross-correlation between images. When images are very similar, this procedure
+is fine; take a look at scikit-image's :code:`skimage.feature.register_translation` for example. 
+
+However, diffraction patterns all have a fixed feature: the position of the beam-block. Therefore, some pixels 
+in each diffraction pattern must be ignored in the computation of the cross-correlation. 
+
+Setting the 'invalid pixels' to 0 will not work, at those will correlate with the invalid pixels from the reference. One must use
+the **masked normalized cross-correlation** through scikit-ued's :code:`mnxc2`.
+
+All of this is taken care of in scikit-ued's :code:`diff_register` function. Let's look at some polycrystalline Chromium:
+
+.. plot::
 
-Let's load a test image::
-	
-	from skimage import img_as_uint
 	from skimage.io import imread
 	import matplotlib.pyplot as plt
 
-	path = '\\data\\vo2.tif'
-	im = img_as_uint(imread(path, plugin = 'tifffile'))
+	ref = imread('Cr_1.tif')
+	im = imread('Cr_2.tif')
 
-	mask = np.zeros_like(im, dtype = np.bool)
-	mask[0:1250, 700:1100] = True
-	im[mask] = 0
+	fig, (ax1, ax2, ax3) = plt.subplots(nrows = 1, ncols = 3, figsize = (9,3))
+	ax1.imshow(ref, vmin = 0, vmax = 200)
+	ax2.imshow(im, vmin = 0, vmax = 200)
+	ax3.imshow(ref - im, cmap = 'RdBu_r')
+
+	for ax in (ax1, ax2, ax3):
+		ax.get_xaxis().set_visible(False)
+		ax.get_yaxis().set_visible(False)
+
+	ax1.set_title('Reference')
+	ax2.set_title('Data')
+	ax3.set_title('Difference')
+
+	plt.tight_layout()
+	plt.show()
+
+From the difference pattern, we can see that the 'Data' pattern is shifted from 'Reference' quite a bit.
+To determine the exact shift, we need to use a mask that obscures the beam-block and main beam::
+
+	from skued.image import diff_register, shift_image
+	import numpy as np
+
+	ref = imread('Cr_1.tif')
+	im = imread('Cr_2.tif')
+
+	mask = np.zeros_like(ref, dtype = np.bool)
+	mask[0:1250, 950:1250] = True
+
+	shift = diff_register(im, reference = ref, mask = mask)
+	im = shift_image(im, shift)
 
-	plt.imshow(im, vmin = 1000, vmax = 1200)
+Let's look at the difference:
+
+.. plot::
+
+	from skimage.io import imread
+	import matplotlib.pyplot as plt
+	import numpy as np
+	from skued.image import diff_register, shift_image
+
+	ref = imread('Cr_1.tif')
+	im = imread('Cr_2.tif')
+
+	mask = np.zeros_like(ref, dtype = np.bool)
+	mask[0:1250, 950:1250] = True
+
+	shift = diff_register(im, ref, mask)
+	shifted = shift_image(im, -shift)
+
+	fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows = 2, ncols = 3, figsize = (9,6))
+	ax1.imshow(ref, vmin = 0, vmax = 200)
+	ax2.imshow(im, vmin = 0, vmax = 200)
+	ax3.imshow(ref - im, cmap = 'RdBu_r')
+	ax4.imshow(mask, vmin = 0, vmax = 1, cmap = 'binary')
+	ax5.imshow(shifted, vmin = 0, vmax = 200)
+	ax6.imshow(ref - shifted, cmap = 'RdBu_r')
+
+	for ax in (ax1, ax2, ax3, ax4, ax5, ax6):
+		ax.get_xaxis().set_visible(False)
+		ax.get_yaxis().set_visible(False)
+
+	ax1.set_title('Reference')
+	ax2.set_title('Data')
+	ax3.set_title('Difference')
+	ax4.set_title('Mask')
+	ax5.set_title('Aligned data')
+	ax6.set_title('Diff. after shift')
+
+	plt.tight_layout()
 	plt.show()
 
+.. _powder:
+
+Image analysis on polycrystalline diffraction patterns
+======================================================
+
+Center-finding
+--------------
+Polycrystalline diffraction patterns display concentric rings, and finding
+the center of those concentric rings is important. Let's load a test image:
+
 .. plot::
 
-	from skimage import img_as_uint
 	from skimage.io import imread
 	import matplotlib.pyplot as plt
-	path = 'data\\vo2.tif'
-	im = img_as_uint(imread(path, plugin = 'tifffile'))
+	path = 'Cr_1.tif'
+
+	im = imread(path, plugin = 'tifffile')
 	mask = np.zeros_like(im, dtype = np.bool)
-	mask[0:1250, 700:1100] = True
+	mask[0:1250, 950:1250] = True
+
 	im[mask] = 0
-	plt.imshow(im, vmin = 1000, vmax = 1200)
+	plt.imshow(im, vmin = 0, vmax = 200)
 	plt.show()
 
 This is a noisy diffraction pattern of polycrystalline vanadium dioxide. 
@@ -163,20 +248,19 @@ Finding the center of such a symmetry pattern can be done with the
 
 .. plot::
 
-	from skimage import img_as_uint
 	from skimage.io import imread
 	import numpy as np
 	import matplotlib.pyplot as plt
-	path = 'data\\vo2.tif'
-	im = img_as_uint(imread(path, plugin = 'tifffile'))
+	path = 'Cr_1.tif'
+	im = imread(path, plugin = 'tifffile')
 	from skued.image import powder_center
 	mask = np.zeros_like(im, dtype = np.bool)
-	mask[0:1250, 700:1100] = True
+	mask[0:1250, 950:1250] = True
 	ic, jc = powder_center(im, mask = mask)
 	ii, jj = np.meshgrid(np.arange(im.shape[0]), np.arange(im.shape[1]),indexing = 'ij')
 	rr = np.sqrt((ii - ic)**2 + (jj - jc)**2)
-	im[rr < 100] = 0
-	plt.imshow(im, vmin = 1000, vmax = 1200)
+	im[rr < 100] = 1e6
+	plt.imshow(im, vmin = 0, vmax = 200)
 	plt.show()
 
 Angular average

external/MaskedFFTRegistrationCode/MaskedTranslationRegistration.m

@@ -0,0 +1,47 @@
+function [transform,maxC,C,numberOfOverlapMaskedPixels] = MaskedTranslationRegistration(fixedImage,movingImage,fixedMask,movingMask,overlapRatio)
+
+% [transform,maxC,C,numberOfOverlapMaskedPixels] =
+% MaskedTranslationRegistration(fixedImage,movingImage,fixedMask,movingMask,overlapRatio) 
+%   Masked FFT normalized cross-correlation registration of movingImage and
+%   fixedImage under masks movingMask and fixedMask.
+%   movingMask and fixedMask should consist of only 1s and 0s, where 1
+%   indicates locations of useful information in the corresponding image,
+%   and 0 indicates locations that should be masked (ignored).
+%   fixedImage and movingImage need not be the same size, but fixedMask
+%   must be the same size as fixedImage, and movingMask must be the same
+%   size as movingImage.
+%   If a mask is not needed for either the fixedImage or the movingImage,
+%   the fixedMask and/or movingMask can be set to an image of all ones of
+%   the same size as the corresponding fixedImage and/or movingImage.
+%   The optional overlapRatio specifies the number of pixels needed in the
+%   overlap region for meaningful results.  It is specified as a ratio of the
+%   maximum number of overlap pixels.  Regions in the resulting correlation
+%   image that have fewer than this number of pixels will be set to 0.   
+%
+%   References: 
+%   D. Padfield. "Masked Object Registration in the Fourier Domain".
+%   Transactions on Image Processing. 
+%   D. Padfield. "Masked FFT registration". In Proc. Computer Vision and
+%   Pattern Recognition, 2010. 
+%
+%   Author: Dirk Padfield, GE Global Research, padfield@research.ge.com
+%
+
+if( nargin < 5 )
+    overlapRatio = 3/10;
+end
+
+[C,numberOfOverlapMaskedPixels] = normxcorr2_masked(fixedImage,movingImage,fixedMask,movingMask);
+
+imageSize = size(movingImage);
+
+% Mask the borders;
+numberOfPixelsThreshold = overlapRatio * max(numberOfOverlapMaskedPixels(:));
+C(numberOfOverlapMaskedPixels < numberOfPixelsThreshold) = 0;
+
+[maxC, imax] = max(C(:));
+[ypeak, xpeak] = ind2sub(size(C),imax(1));
+transform = [(xpeak-imageSize(2)) (ypeak-imageSize(1))];
+
+% Take the negative of the transform so that it has the correct sign.
+transform = -transform;

external/MaskedFFTRegistrationCode/MaskedTranslationRegistrationTest.m

@@ -0,0 +1,42 @@
+% MaskedTranslationRegistrationTest
+%   Test the MaskedTranslationRegistration code.
+%
+%   References: 
+%   D. Padfield. "Masked Object Registration in the Fourier Domain".
+%   Transactions on Image Processing. 
+%   D. Padfield. "Masked FFT registration". In Proc. Computer Vision and
+%   Pattern Recognition, 2010. 
+%
+%   Author: Dirk Padfield, GE Global Research, padfield@research.ge.com
+%
+
+
+close all;
+clear variables;
+clc;
+
+% Perform the registration on several sets of images.
+x = [75 -130 130];
+y = [75 130 130];
+overlapRatio = 1/10;
+
+for i = 1:3
+    fixedImage = imread(sprintf('OriginalX%2iY%2i.png',x(i),y(i)));
+    movingImage = imread(sprintf('TransformedX%2iY%2i.png',x(i),y(i)));
+    fixedMask = fixedImage~=0;
+    movingMask = movingImage~=0;
+
+    [translation,maxC] = MaskedTranslationRegistration(fixedImage,movingImage,fixedMask,movingMask,overlapRatio);
+    [transformedMovingImage] = TransformImageTranslation(movingImage,translation);
+
+    % Given the transform, transform the moving image.
+    [overlayImage] = OverlayRegistration(fixedImage,transformedMovingImage);
+    figure; imagesc(overlayImage); title(['Test ' num2str(i) ': Registered Overlay Image']);
+
+    disp(['Test ' num2str(i) ':']);
+    disp(['Computed translation: ' num2str([translation(1) -translation(2)])]);
+    disp(['Correlation score: ' num2str(maxC)]);
+    trueTranslation = [x(i),-y(i)];
+    disp(['Transformation error: ' num2str(translation - trueTranslation)]);
+    disp(' ');
+end