diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 8b2ff9ff..548f37c5 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -44,7 +44,8 @@ Incorporating CellProfiler's Sobel edge detector, build and bug fixes. - Emmanuelle Guillart - Total variation noise filtering + Total variation noise filtering, integration of CellProfiler's + mathematical morphology tools. - Maël Primet Total variation noise filtering @@ -60,3 +61,6 @@ - Kyle Mandli CSV to ReST code for feature comparison table. + +- The Scikit Learn team + From whom we borrowed the example generation tools. diff --git a/doc/examples/README.txt b/doc/examples/README.txt new file mode 100644 index 00000000..8a6c6bd5 --- /dev/null +++ b/doc/examples/README.txt @@ -0,0 +1,6 @@ + +General examples +------------------- + +General-purpose and introductory examples for the scikit. + diff --git a/doc/examples/plot_lena_tv_denoise.py b/doc/examples/plot_lena_tv_denoise.py new file mode 100644 index 00000000..55f167e1 --- /dev/null +++ b/doc/examples/plot_lena_tv_denoise.py @@ -0,0 +1,52 @@ +""" +==================================================== +Denoising the picture of Lena using total variation +==================================================== + +In this example, we denoise a noisy version of the picture of Lena +using the total variation denoising filter. The result of this filter +is an image that has a minimal total variation norm, while being as +close to the initial image as possible. The total variation is the L1 +norm of the gradient of the image, and minimizing the total variation +typically produces "posterized" images with flat domains separated by +sharp edges. + +It is possible to change the degree of posterization by controlling +the tradeoff between denoising and faithfulness to the original image. + +""" + +import numpy as np +import matplotlib.pyplot as plt + +from scikits.image import data +from scikits.image.filter import tv_denoise + +l = data.lena() +l = l[230:290, 220:320] + +noisy = l + 0.4*l.std()*np.random.random(l.shape) + +tv_denoised = tv_denoise(noisy, weight=10) + + +plt.figure(figsize=(12,2.8)) + +plt.subplot(131) +plt.imshow(noisy, cmap=plt.cm.gray, vmin=40, vmax=220) +plt.axis('off') +plt.title('noisy', fontsize=20) +plt.subplot(132) +plt.imshow(tv_denoised, cmap=plt.cm.gray, vmin=40, vmax=220) +plt.axis('off') +plt.title('TV denoising', fontsize=20) + +tv_denoised = tv_denoise(noisy, weight=50) +plt.subplot(133) +plt.imshow(tv_denoised, cmap=plt.cm.gray, vmin=40, vmax=220) +plt.axis('off') +plt.title('(more) TV denoising', fontsize=20) + +plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0, left=0, + right=1) +plt.show() diff --git a/doc/ext/LICENSE.txt b/doc/ext/LICENSE.txt index c23b6bfc..3ebab893 100644 --- a/doc/ext/LICENSE.txt +++ b/doc/ext/LICENSE.txt @@ -66,3 +66,46 @@ products or services of Licensee, or any third party. 8. By copying, installing or otherwise using matplotlib 0.98.3, Licensee agrees to be bound by the terms and conditions of this License Agreement. + +The file + +- gen_rst.py + +was taken from the scikit-learn (http://scikit-learn.sourceforge.net), which +has the following license: + +New BSD License + +Copyright (c) 2007 - 2011 The scikit-learn developers. +All rights reserved. + + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + a. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + b. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + c. Neither the name of the Scikit-learn Developers nor the names of + its contributors may be used to endorse or promote products + derived from this software without specific prior written + permission. + + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE. + + + + diff --git a/doc/ext/gen_rst.py b/doc/ext/gen_rst.py new file mode 100644 index 00000000..0cf8cbd2 --- /dev/null +++ b/doc/ext/gen_rst.py @@ -0,0 +1,305 @@ +""" +Example generation for the scikit image. + +Generate the rst files for the examples by iterating over the python +example files. + +Files that generate images should start with 'plot'. + +This code was taken from the scikit-learn. + +""" +import os +import shutil +import traceback +import glob + +import matplotlib +matplotlib.use('Agg') + +import token, tokenize + +rst_template = """ + +.. _example_%(short_fname)s: + +%(docstring)s + +**Python source code:** :download:`%(fname)s <%(fname)s>` + +.. literalinclude:: %(fname)s + :lines: %(end_row)s- + """ + +plot_rst_template = """ + +.. _example_%(short_fname)s: + +%(docstring)s + +%(image_list)s + +**Python source code:** :download:`%(fname)s <%(fname)s>` + +.. literalinclude:: %(fname)s + :lines: %(end_row)s- + """ + +# The following strings are used when we have several pictures: we use +# an html div tag that our CSS uses to turn the lists into horizontal +# lists. +HLIST_HEADER = """ +.. rst-class:: horizontal + +""" + +HLIST_IMAGE_TEMPLATE = """ + * + + .. image:: images/%s + :scale: 50 +""" + +SINGLE_IMAGE = """ +.. image:: images/%s + :align: center +""" + +def extract_docstring(filename): + """ Extract a module-level docstring, if any + """ + lines = file(filename).readlines() + start_row = 0 + if lines[0].startswith('#!'): + lines.pop(0) + start_row = 1 + + docstring = '' + first_par = '' + tokens = tokenize.generate_tokens(lines.__iter__().next) + for tok_type, tok_content, _, (erow, _), _ in tokens: + tok_type = token.tok_name[tok_type] + if tok_type in ('NEWLINE', 'COMMENT', 'NL', 'INDENT', 'DEDENT'): + continue + elif tok_type == 'STRING': + docstring = eval(tok_content) + # If the docstring is formatted with several paragraphs, extract + # the first one: + paragraphs = '\n'.join(line.rstrip() + for line in docstring.split('\n')).split('\n\n') + if len(paragraphs) > 0: + first_par = paragraphs[0] + break + return docstring, first_par, erow+1+start_row + + +def generate_example_rst(app): + """ Generate the list of examples, as well as the contents of + examples. + """ + root_dir = os.path.join(app.builder.srcdir, 'auto_examples') + example_dir = os.path.abspath(app.builder.srcdir + '/../' + 'examples') + try: + plot_gallery = eval(app.builder.config.plot_gallery) + except TypeError: + plot_gallery = bool(app.builder.config.plot_gallery) + if not os.path.exists(example_dir): + os.makedirs(example_dir) + if not os.path.exists(root_dir): + os.makedirs(root_dir) + + # we create an index.rst with all examples + fhindex = file(os.path.join(root_dir, 'index.txt'), 'w') + fhindex.write("""\ + +.. raw:: html + + + +Examples +======== + +.. _examples-index: +""") + # Here we don't use an os.walk, but we recurse only twice: flat is + # better than nested. + generate_dir_rst('.', fhindex, example_dir, root_dir, plot_gallery) + for dir in sorted(os.listdir(example_dir)): + if os.path.isdir(os.path.join(example_dir, dir)): + generate_dir_rst(dir, fhindex, example_dir, root_dir, plot_gallery) + fhindex.flush() + + +def generate_dir_rst(dir, fhindex, example_dir, root_dir, plot_gallery): + """ Generate the rst file for an example directory. + """ + if not dir == '.': + target_dir = os.path.join(root_dir, dir) + src_dir = os.path.join(example_dir, dir) + else: + target_dir = root_dir + src_dir = example_dir + if not os.path.exists(os.path.join(src_dir, 'README.txt')): + print 80*'_' + print ('Example directory %s does not have a README.txt file' + % src_dir) + print 'Skipping this directory' + print 80*'_' + return + fhindex.write(""" + + +%s + + +""" % file(os.path.join(src_dir, 'README.txt')).read()) + if not os.path.exists(target_dir): + os.makedirs(target_dir) + + def sort_key(a): + # put last elements without a plot + if not a.startswith('plot') and a.endswith('.py'): + return 'zz' + a + return a + for fname in sorted(os.listdir(src_dir), key=sort_key): + if fname.endswith('py'): + generate_file_rst(fname, target_dir, src_dir, plot_gallery) + thumb = os.path.join(dir, 'images', 'thumb', fname[:-3] + '.png') + link_name = os.path.join(dir, fname).replace(os.path.sep, '_') + fhindex.write('.. figure:: %s\n' % thumb) + if link_name.startswith('._'): + link_name = link_name[2:] + if dir != '.': + fhindex.write(' :target: ./%s/%s.html\n\n' % (dir, fname[:-3])) + else: + fhindex.write(' :target: ./%s.html\n\n' % link_name[:-3]) + fhindex.write(' :ref:`example_%s`\n\n' % link_name) + fhindex.write(""" +.. raw:: html + +
+ """) # clear at the end of the section + + +def generate_file_rst(fname, target_dir, src_dir, plot_gallery): + """ Generate the rst file for a given example. + """ + base_image_name = os.path.splitext(fname)[0] + image_fname = '%s_%%s.png' % base_image_name + + this_template = rst_template + last_dir = os.path.split(src_dir)[-1] + # to avoid leading . in file names, and wrong names in links + if last_dir == '.' or last_dir == 'examples': + last_dir = '' + else: + last_dir += '_' + short_fname = last_dir + fname + src_file = os.path.join(src_dir, fname) + example_file = os.path.join(target_dir, fname) + shutil.copyfile(src_file, example_file) + + # The following is a list containing all the figure names + figure_list = [] + + image_dir = os.path.join(target_dir, 'images') + thumb_dir = os.path.join(image_dir, 'thumb') + if not os.path.exists(image_dir): + os.makedirs(image_dir) + if not os.path.exists(thumb_dir): + os.makedirs(thumb_dir) + image_path = os.path.join(image_dir, image_fname) + thumb_file = os.path.join(thumb_dir, fname[:-3] + '.png') + if plot_gallery and fname.startswith('plot'): + # generate the plot as png image if file name + # starts with plot and if it is more recent than an + # existing image. + first_image_file = image_path % 1 + + if (not os.path.exists(first_image_file) or + os.stat(first_image_file).st_mtime <= + os.stat(src_file).st_mtime): + # We need to execute the code + print 'plotting %s' % fname + import matplotlib.pyplot as plt + plt.close('all') + cwd = os.getcwd() + try: + # First CD in the original example dir, so that any file created + # by the example get created in this directory + os.chdir(os.path.dirname(src_file)) + execfile(os.path.basename(src_file), {'pl' : plt}) + os.chdir(cwd) + + # In order to save every figure we have two solutions : + # * iterate from 1 to infinity and call plt.fignum_exists(n) + # (this requires the figures to be numbered + # incrementally: 1, 2, 3 and not 1, 2, 5) + # * iterate over [fig_mngr.num for fig_mngr in + # matplotlib._pylab_helpers.Gcf.get_all_fig_managers()] + for fig_num in (fig_mngr.num for fig_mngr in + matplotlib._pylab_helpers.Gcf.get_all_fig_managers()): + # Set the fig_num figure as the current figure as we can't + # save a figure that's not the current figure. + plt.figure(fig_num) + plt.savefig(image_path % fig_num) + figure_list.append(image_fname % fig_num) + except: + print 80*'_' + print '%s is not compiling:' % fname + traceback.print_exc() + print 80*'_' + finally: + os.chdir(cwd) + else: + figure_list = [f[len(image_dir):] + for f in glob.glob(image_path % '[1-9]')] + #for f in glob.glob(image_path % '*')] + + # generate thumb file + this_template = plot_rst_template + from matplotlib import image + if os.path.exists(first_image_file): + image.thumbnail(first_image_file, thumb_file, 0.2) + + if not os.path.exists(thumb_file): + # create something not to replace the thumbnail + shutil.copy('images/blank_image.png', thumb_file) + + docstring, short_desc, end_row = extract_docstring(example_file) + + # Depending on whether we have one or more figures, we're using a + # horizontal list or a single rst call to 'image'. + if len(figure_list) == 1: + figure_name = figure_list[0] + image_list = SINGLE_IMAGE % figure_name.lstrip('/') + else: + image_list = HLIST_HEADER + for figure_name in figure_list: + image_list += HLIST_IMAGE_TEMPLATE % figure_name.lstrip('/') + + f = open(os.path.join(target_dir, fname[:-2] + 'txt'),'w') + f.write(this_template % locals()) + f.flush() + + +def setup(app): + app.connect('builder-inited', generate_example_rst) + app.add_config_value('plot_gallery', True, 'html') diff --git a/doc/source/auto_examples/images/blank_image.png b/doc/source/auto_examples/images/blank_image.png new file mode 100644 index 00000000..4913b99b Binary files /dev/null and b/doc/source/auto_examples/images/blank_image.png differ diff --git a/doc/source/conf.py b/doc/source/conf.py index c9959dd4..d9e20793 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -23,11 +23,17 @@ sys.path.append(os.path.join(curpath, '..', 'ext')) # -- General configuration ----------------------------------------------------- +try: + import gen_rst +except: + pass + + # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. extensions = ['sphinx.ext.autodoc', 'sphinx.ext.pngmath', 'numpydoc', 'sphinx.ext.autosummary', 'sphinx.ext.inheritance_diagram', - 'plot_directive'] + 'plot_directive', 'gen_rst'] # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] diff --git a/doc/source/coverage_generator.py b/doc/source/coverage_generator.py old mode 100644 new mode 100755 index 25ab388a..83fd5d86 --- a/doc/source/coverage_generator.py +++ b/doc/source/coverage_generator.py @@ -40,16 +40,16 @@ def calculate_coverage(reader): total_items += 1 if ":done:" in row[2] or ":done:" in row[3]: done_items += 1 - elif ":partial:" in row[2] or ":partial:" in row[3]: + if ":partial:" in row[2] or ":partial:" in row[3]: partial_items += 1 - elif ":na:" in row[2] or ":na:" in row[3]: + if ":na:" in row[2] or ":na:" in row[3]: na_items += 1 counts = (done_items, partial_items, total_items - (partial_items + done_items + na_items), na_items) - + return list(i / total_items for i in counts) def read_table_titles(reader): @@ -235,7 +235,7 @@ Coverage Bar for item, style in enumerate(('done-bar', 'partial-bar', 'missing-bar', 'na-bar')): - stream.write('XX' % \ + stream.write(' ' % \ (item_counts[item] * 100, style)) stream.write("\n\n") @@ -269,6 +269,6 @@ if __name__ == "__main__": generate_page(csv_path,output) output.close() - print "Generated %s from %s." % (output_path,csv_path) + print("Generated %s from %s." % (output_path,csv_path)) diff --git a/doc/source/index.txt b/doc/source/index.txt index 15b5b57e..2ec428ac 100644 --- a/doc/source/index.txt +++ b/doc/source/index.txt @@ -39,6 +39,11 @@ Sections Conditions on the use and redistribution of this package. + +`Examples `_ + +Introductory examples + Indices and Tables ================== diff --git a/scikits/image/data/__init__.py b/scikits/image/data/__init__.py new file mode 100644 index 00000000..c46a2e37 --- /dev/null +++ b/scikits/image/data/__init__.py @@ -0,0 +1,4 @@ +import scipy.misc as _m + +lena = _m.lena + diff --git a/scikits/image/setup.py b/scikits/image/setup.py index b1831d4e..2d24389d 100644 --- a/scikits/image/setup.py +++ b/scikits/image/setup.py @@ -11,6 +11,7 @@ def configuration(parent_package='', top_path=None): config.add_subpackage('morphology') config.add_subpackage('filter') config.add_subpackage('transform') + config.add_subpackage('data') def add_test_directories(arg, dirname, fnames): if dirname.split(os.path.sep)[-1] == 'tests': diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py index d8adede9..01b8fafb 100644 --- a/scikits/image/transform/__init__.py +++ b/scikits/image/transform/__init__.py @@ -1,5 +1,6 @@ -from hough_transform import * -from finite_radon_transform import * -from radon_transform import * -from project import * - +from .hough_transform import * +from .radon import * +from .finite_radon_transform import * +from .project import * +from ._project import homography as fast_homography +from .integral import * diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx new file mode 100644 index 00000000..3e2b9b7d --- /dev/null +++ b/scikits/image/transform/_project.pyx @@ -0,0 +1,209 @@ +#cython: cdivison=True boundscheck=False + +__all__ = ['homography'] + +cimport cython +cimport numpy as np + +import numpy as np +import cython + +from cython.operator import dereference + +np.import_array() + +cdef extern from "math.h": + double floor(double) + double fmod(double, double) + +cdef double get_pixel(double *image, int rows, int cols, + int r, int c, char mode, double cval=0): + """Get a pixel from the image, taking wrapping mode into consideration. + + Parameters + ---------- + image : *double + Input image. + rows, cols : int + Dimensions of image. + r, c : int + Position at which to get the pixel. + mode : {'C', 'W', 'M'} + Wrapping mode. Constant, Wrap or Mirror. + cval : double + Constant value to use for mode constant. + + """ + if mode == 'C': + if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): + return cval + else: + return image[r * cols + c] + else: + return image[coord_map(rows, r, mode) * cols + + coord_map(cols, c, mode)] + +cdef int coord_map(int dim, int coord, char mode): + """ + Wrap a coordinate, according to a given dimension and mode. + + Parameters + ---------- + dim : int + Maximum coordinate. + coord : int + Coord provided by user. May be < 0 or > dim. + mode : {'W', 'M'} + Whether to wrap or mirror the coordinate if it + falls outside [0, dim). + + """ + dim = dim - 1 + if mode == 'M': # mirror + if (coord < 0): + # How many times times does the coordinate wrap? + if ((-coord / dim) % 2 != 0): + return dim - (-coord % dim) + else: + return (-coord % dim) + elif (coord > dim): + if ((coord / dim) % 2 != 0): + return (dim - (coord % dim)) + else: + return (coord % dim) + elif mode == 'W': # wrap + if (coord < 0): + return (dim - (-coord % dim)) + elif (coord > dim): + return (coord % dim) + + return coord + +cdef tf(double x, double y, double* H, double *x_, double *y_): + """Apply a homography to a coordinate. + + Parameters + ---------- + x, y : double + Input coordinate. + H : (3,3) *double + Transformation matrix. + x_, y_ : *double + Output coordinate. + + """ + cdef double xx, yy, zz + + xx = H[0] * x + H[1] * y + H[2] + yy = H[3] * x + H[4] * y + H[5] + zz = H[6] * x + H[7] * y + H[8] + + xx = xx / zz + yy = yy / zz + + x_[0] = xx + y_[0] = yy + +@cython.boundscheck(False) +def homography(np.ndarray image, np.ndarray H, output_shape=None, + mode='constant', double cval=0): + """ + Projective transformation (homography). + + Perform a projective transformation (homography) of a + floating point image, using bi-linear interpolation. + + For each pixel, given its homogeneous coordinate :math:`\mathbf{x} + = [x, y, 1]^T`, its target position is calculated by multiplying + with the given matrix, :math:`H`, to give :math:`H \mathbf{x}`. + E.g., to rotate by theta degrees clockwise, the matrix should be + + :: + + [[cos(theta) -sin(theta) 0] + [sin(theta) cos(theta) 0] + [0 0 1]] + + or, to translate x by 10 and y by 20, + + :: + + [[1 0 10] + [0 1 20] + [0 0 1 ]]. + + Parameters + ---------- + image : 2-D array + Input image. + H : array of shape ``(3, 3)`` + Transformation matrix H that defines the homography. + output_shape : tuple (rows, cols) + Shape of the output image generated. + order : int + Order of splines used in interpolation. + mode : {'constant', 'mirror', 'wrap'} + How to handle values outside the image borders. + cval : string + Used in conjunction with mode 'C' (constant), the value + outside the image boundaries. + + """ + + cdef np.ndarray[dtype=np.double_t, ndim=2] img = \ + np.asarray(image, dtype=np.double) + cdef np.ndarray[dtype=np.double_t, ndim=2] M = \ + np.ascontiguousarray(np.linalg.inv(H)) + + if mode not in ('constant', 'wrap', 'mirror'): + raise ValueError("Invalid mode specified. Please use " + "`constant`, `wrap` or `mirror`.") + if mode == 'constant': + mode_c = ord('C') + elif mode == 'wrap': + mode_c = ord('W') + elif mode == 'mirror': + mode_c = ord('M') + + cdef int out_r, out_c, columns, rows + if output_shape is None: + out_r = img.shape[0] + out_c = img.shape[1] + else: + out_r = output_shape[0] + out_c = output_shape[1] + + rows = img.shape[0] + columns = img.shape[1] + + cdef np.ndarray[dtype=np.double_t, ndim=2] out = \ + np.zeros((out_r, out_c), dtype=np.double) + + cdef int tfr, tfc, r_int, c_int + cdef double y0, y1, y2, y3 + cdef double r, c, z, t, u + + for tfr in range(out_r): + for tfc in range(out_c): + tf(tfc, tfr, M.data, &c, &r) + r_int = floor(r) + c_int = floor(c) + + t = r - r_int + u = c - c_int + + y0 = get_pixel(img.data, rows, columns, + r_int, c_int, mode_c) + y1 = get_pixel(img.data, rows, columns, + r_int + 1, c_int, mode_c) + y2 = get_pixel(img.data, rows, columns, + r_int + 1, c_int + 1, mode_c) + y3 = get_pixel(img.data, rows, columns, + r_int, c_int + 1, mode_c) + + out[tfr, tfc] = \ + (1 - t) * (1 - u) * y0 + \ + t * (1 - u) * y1 + \ + t * u * y2 + (1 - t) * u * y3; + + return out diff --git a/scikits/image/transform/integral.py b/scikits/image/transform/integral.py new file mode 100644 index 00000000..440d10e6 --- /dev/null +++ b/scikits/image/transform/integral.py @@ -0,0 +1,60 @@ +def integral_image(x): + """Integral image / summed area table. + + The integral image contains the sum of all elements above and to the + left of it, i.e.: + + .. math:: + + S[m, n] = \sum_{i \leq m} \sum_{j \leq n} X[i, j] + + Parameters + ---------- + x : ndarray + Input image. + + Returns + ------- + S : ndarray + Integral image / summed area table. + + References + ---------- + .. [1] F.C. Crow, "Summed-area tables for texture mapping," + ACM SIGGRAPH Computer Graphics, vol. 18, 1984, pp. 207-212. + + """ + return x.cumsum(1).cumsum(0) + +def integrate(ii, r0, c0, r1, c1): + """Use an integral image to integrate over a given window. + + Parameters + ---------- + ii : ndarray + Integral image. + r0, c0 : int + Top-left corner of block to be summed. + r1, c1 : int + Bottom-right corner of block to be summed. + + Returns + ------- + S : int + Integral (sum) over the given window. + + """ + S = 0 + + S += ii[r1, c1] + + if (r0 - 1 >= 0) and (c0 - 1 >= 0): + S += ii[r0 - 1, c0 - 1] + + if (r0 - 1 >= 0): + S -= ii[r0 - 1, c1] + + if (c0 - 1 >= 0): + S -= ii[r1, c0 - 1] + + return S diff --git a/scikits/image/transform/setup.py b/scikits/image/transform/setup.py index f15dd08c..c5629eb0 100644 --- a/scikits/image/transform/setup.py +++ b/scikits/image/transform/setup.py @@ -15,10 +15,14 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['_hough_transform.pyx'], working_path=base_path) + cython(['_project.pyx'], working_path=base_path) config.add_extension('_hough_transform', sources=['_hough_transform.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_project', sources=['_project.c'], + include_dirs=[get_numpy_include_dirs()]) + return config if __name__ == '__main__': diff --git a/scikits/image/transform/tests/test_integral.py b/scikits/image/transform/tests/test_integral.py new file mode 100644 index 00000000..65cdd632 --- /dev/null +++ b/scikits/image/transform/tests/test_integral.py @@ -0,0 +1,28 @@ +import numpy as np +from numpy.testing import * + +from scikits.image.transform import integral_image, integrate + +x = (np.random.random((50, 50)) * 255).astype(np.uint8) +s = integral_image(x) + +def test_validity(): + y = np.arange(12).reshape((4,3)) + + y = (np.random.random((50, 50)) * 255).astype(np.uint8) + assert_equal(integral_image(y)[-1, -1], + y.sum()) + +def test_basic(): + assert_equal(x[12:24, 10:20].sum(), integrate(s, 12, 10, 23, 19)) + assert_equal(x[:20, :20].sum(), integrate(s, 0, 0, 19, 19)) + assert_equal(x[:20, 10:20].sum(), integrate(s, 0, 10, 19, 19)) + assert_equal(x[10:20, :20].sum(), integrate(s, 10, 0, 19, 19)) + +def test_single(): + assert_equal(x[0, 0], integrate(s, 0, 0, 0, 0)) + assert_equal(x[10, 10], integrate(s, 10, 10, 10, 10)) + + +if __name__ == '__main__': + run_module_suite() diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py index 44c0e557..ba8e277e 100644 --- a/scikits/image/transform/tests/test_project.py +++ b/scikits/image/transform/tests/test_project.py @@ -1,7 +1,9 @@ import numpy as np from numpy.testing import assert_array_almost_equal -from scikits.image.transform.project import _stackcopy, homography +from scikits.image.transform.project import _stackcopy +from scikits.image.transform import homography, fast_homography +from scikits.image import data def test_stackcopy(): layers = 4 @@ -19,3 +21,41 @@ def test_homography(): [0, 0, 1]]) x90 = homography(x, M, order=1) assert_array_almost_equal(x90, np.rot90(x)) + +def test_fast_homography(): + img = data.lena() + img = img[:, :100] + + theta = np.deg2rad(30) + scale = 0.5 + tx, ty = 50, 50 + + H = np.eye(3) + S = scale * np.sin(theta) + C = scale * np.cos(theta) + + H[:2, :2] = [[C, -S], [S, C]] + H[:2, 2] = [tx, ty] + + for mode in ('constant', 'mirror', 'wrap'): + print 'Transform mode:', mode + + p0 = homography(img, H, mode=mode, order=1) + p1 = fast_homography(img, H, mode=mode) + p1 = np.round(p1) + + ## import matplotlib.pyplot as plt + ## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4) + ## ax0.imshow(img) + ## ax1.imshow(p0, cmap=plt.cm.gray) + ## ax2.imshow(p1, cmap=plt.cm.gray) + ## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray) + ## plt.show() + + d = np.mean(np.abs(p0 - p1)) + assert d < 0.2 + + +if __name__ == "__main__": + from numpy.testing import run_module_suite + run_module_suite()