From b56c370dfaf9a3203259ca4ae26959bea9ca9674 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sat, 17 Sep 2011 09:53:08 -0700 Subject: [PATCH 01/14] Add summed area table, using Nicolas Pinto's suggestion for a pure Python implementation. --- scikits/image/transform/__init__.py | 8 +-- scikits/image/transform/sum_area_table.py | 61 +++++++++++++++++++++++ scikits/image/transform/tests/test_sat.py | 28 +++++++++++ 3 files changed, 93 insertions(+), 4 deletions(-) create mode 100644 scikits/image/transform/sum_area_table.py create mode 100644 scikits/image/transform/tests/test_sat.py diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py index 91711c6e..407ab0d4 100644 --- a/scikits/image/transform/__init__.py +++ b/scikits/image/transform/__init__.py @@ -1,4 +1,4 @@ -from hough_transform import * -from finite_radon_transform import * -from project import * - +from .hough_transform import * +from .finite_radon_transform import * +from .project import * +from .sum_area_table import * diff --git a/scikits/image/transform/sum_area_table.py b/scikits/image/transform/sum_area_table.py new file mode 100644 index 00000000..6811cdc5 --- /dev/null +++ b/scikits/image/transform/sum_area_table.py @@ -0,0 +1,61 @@ +def sat(x): + """Summed area table / integral image. + + 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 of uint + Input image. + + Returns + ------- + S : ndarray + 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 sat_sum(sat, r0, c0, r1, c1): + """Using a summed area table / integral image, calculate the sum + over a given window. + + Parameters + ---------- + sat : ndarray of uint64 + Summed area table / 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 + Sum over the given window. + + """ + S = 0 + + S += sat[r1, c1] + + if (r0 - 1 >= 0) and (c0 - 1 >= 0): + S += sat[r0 - 1, c0 - 1] + + if (r0 - 1 >= 0): + S -= sat[r0 - 1, c1] + + if (c0 - 1 >= 0): + S -= sat[r1, c0 - 1] + + return S diff --git a/scikits/image/transform/tests/test_sat.py b/scikits/image/transform/tests/test_sat.py new file mode 100644 index 00000000..481d8e28 --- /dev/null +++ b/scikits/image/transform/tests/test_sat.py @@ -0,0 +1,28 @@ +import numpy as np +from numpy.testing import * + +from scikits.image.transform import sat, sat_sum + +x = (np.random.random((50, 50)) * 255).astype(np.uint8) +s = sat(x) + +def test_validity(): + y = np.arange(12).reshape((4,3)) + + y = (np.random.random((50, 50)) * 255).astype(np.uint8) + assert_equal(sat(y)[-1, -1], + y.sum()) + +def test_basic(): + assert_equal(x[12:24, 10:20].sum(), sat_sum(s, 12, 10, 23, 19)) + assert_equal(x[:20, :20].sum(), sat_sum(s, 0, 0, 19, 19)) + assert_equal(x[:20, 10:20].sum(), sat_sum(s, 0, 10, 19, 19)) + assert_equal(x[10:20, :20].sum(), sat_sum(s, 10, 0, 19, 19)) + +def test_single(): + assert_equal(x[0, 0], sat_sum(s, 0, 0, 0, 0)) + assert_equal(x[10, 10], sat_sum(s, 10, 10, 10, 10)) + + +if __name__ == '__main__': + run_module_suite() From b2be624de1ed957358326cc42e7ea098fe9685c0 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 18 Sep 2011 14:20:38 -0700 Subject: [PATCH 02/14] DOC: Update sat docs after feedback from Brian Holt. --- scikits/image/transform/sum_area_table.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/scikits/image/transform/sum_area_table.py b/scikits/image/transform/sum_area_table.py index 6811cdc5..6496d7d3 100644 --- a/scikits/image/transform/sum_area_table.py +++ b/scikits/image/transform/sum_area_table.py @@ -10,7 +10,7 @@ def sat(x): Parameters ---------- - X : ndarray of uint + X : ndarray Input image. Returns @@ -27,12 +27,11 @@ def sat(x): return x.cumsum(1).cumsum(0) def sat_sum(sat, r0, c0, r1, c1): - """Using a summed area table / integral image, calculate the sum - over a given window. + """Use a summed area table / integral image to sum over a given window. Parameters ---------- - sat : ndarray of uint64 + sat : ndarray Summed area table / integral image. r0, c0 : int Top-left corner of block to be summed. From ac0fc3aba8950ef43206562577d0194eb894b067 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Mon, 19 Sep 2011 21:45:55 -0500 Subject: [PATCH 03/14] Fix display of coverage bar --- doc/source/coverage_generator.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) mode change 100644 => 100755 doc/source/coverage_generator.py 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)) From 5ba4ebf30ba3e086214ba3719df0381e76aa74ba Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Tue, 20 Sep 2011 08:43:54 -0700 Subject: [PATCH 04/14] ENH: Rename sat and sat_sum to be in-line with CV community convention. --- scikits/image/transform/__init__.py | 2 +- .../{sum_area_table.py => integral.py} | 26 ++++++++--------- .../image/transform/tests/test_integral.py | 28 +++++++++++++++++++ scikits/image/transform/tests/test_sat.py | 28 ------------------- 4 files changed, 42 insertions(+), 42 deletions(-) rename scikits/image/transform/{sum_area_table.py => integral.py} (65%) create mode 100644 scikits/image/transform/tests/test_integral.py delete mode 100644 scikits/image/transform/tests/test_sat.py diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py index 407ab0d4..341671d1 100644 --- a/scikits/image/transform/__init__.py +++ b/scikits/image/transform/__init__.py @@ -1,4 +1,4 @@ from .hough_transform import * from .finite_radon_transform import * from .project import * -from .sum_area_table import * +from .integral import * diff --git a/scikits/image/transform/sum_area_table.py b/scikits/image/transform/integral.py similarity index 65% rename from scikits/image/transform/sum_area_table.py rename to scikits/image/transform/integral.py index 6496d7d3..440d10e6 100644 --- a/scikits/image/transform/sum_area_table.py +++ b/scikits/image/transform/integral.py @@ -1,5 +1,5 @@ -def sat(x): - """Summed area table / integral image. +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.: @@ -10,13 +10,13 @@ def sat(x): Parameters ---------- - X : ndarray + x : ndarray Input image. Returns ------- S : ndarray - Summed area table. + Integral image / summed area table. References ---------- @@ -26,13 +26,13 @@ def sat(x): """ return x.cumsum(1).cumsum(0) -def sat_sum(sat, r0, c0, r1, c1): - """Use a summed area table / integral image to sum over a given window. +def integrate(ii, r0, c0, r1, c1): + """Use an integral image to integrate over a given window. Parameters ---------- - sat : ndarray - Summed area table / integral image. + ii : ndarray + Integral image. r0, c0 : int Top-left corner of block to be summed. r1, c1 : int @@ -41,20 +41,20 @@ def sat_sum(sat, r0, c0, r1, c1): Returns ------- S : int - Sum over the given window. + Integral (sum) over the given window. """ S = 0 - S += sat[r1, c1] + S += ii[r1, c1] if (r0 - 1 >= 0) and (c0 - 1 >= 0): - S += sat[r0 - 1, c0 - 1] + S += ii[r0 - 1, c0 - 1] if (r0 - 1 >= 0): - S -= sat[r0 - 1, c1] + S -= ii[r0 - 1, c1] if (c0 - 1 >= 0): - S -= sat[r1, c0 - 1] + S -= ii[r1, c0 - 1] return S 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_sat.py b/scikits/image/transform/tests/test_sat.py deleted file mode 100644 index 481d8e28..00000000 --- a/scikits/image/transform/tests/test_sat.py +++ /dev/null @@ -1,28 +0,0 @@ -import numpy as np -from numpy.testing import * - -from scikits.image.transform import sat, sat_sum - -x = (np.random.random((50, 50)) * 255).astype(np.uint8) -s = sat(x) - -def test_validity(): - y = np.arange(12).reshape((4,3)) - - y = (np.random.random((50, 50)) * 255).astype(np.uint8) - assert_equal(sat(y)[-1, -1], - y.sum()) - -def test_basic(): - assert_equal(x[12:24, 10:20].sum(), sat_sum(s, 12, 10, 23, 19)) - assert_equal(x[:20, :20].sum(), sat_sum(s, 0, 0, 19, 19)) - assert_equal(x[:20, 10:20].sum(), sat_sum(s, 0, 10, 19, 19)) - assert_equal(x[10:20, :20].sum(), sat_sum(s, 10, 0, 19, 19)) - -def test_single(): - assert_equal(x[0, 0], sat_sum(s, 0, 0, 0, 0)) - assert_equal(x[10, 10], sat_sum(s, 10, 10, 10, 10)) - - -if __name__ == '__main__': - run_module_suite() From cce3d7bcc5f285b17e2accce5661cdf8c0c84e92 Mon Sep 17 00:00:00 2001 From: Emmanuelle Gouillart Date: Fri, 23 Sep 2011 23:53:47 +0200 Subject: [PATCH 05/14] Example generation from scripts, taken from scikit learn (ext/gen_rst.py) --- doc/examples/README.txt | 6 + doc/examples/plot_lena_tv_denoise.py | 51 +++++ doc/ext/gen_rst.py | 303 +++++++++++++++++++++++++++ doc/source/conf.py | 8 +- doc/source/index.txt | 5 + examples/plot_lena_tv_denoise.py | 51 +++++ 6 files changed, 423 insertions(+), 1 deletion(-) create mode 100644 doc/examples/README.txt create mode 100644 doc/examples/plot_lena_tv_denoise.py create mode 100644 doc/ext/gen_rst.py create mode 100644 examples/plot_lena_tv_denoise.py 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..c4000d38 --- /dev/null +++ b/doc/examples/plot_lena_tv_denoise.py @@ -0,0 +1,51 @@ +""" +==================================================== +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 scipy +from scipy import ndimage +import matplotlib.pyplot as plt +from scikits.image.filter import tv_denoise + +l = scipy.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) + diff --git a/doc/ext/gen_rst.py b/doc/ext/gen_rst.py new file mode 100644 index 00000000..cb71a605 --- /dev/null +++ b/doc/ext/gen_rst.py @@ -0,0 +1,303 @@ +""" +Example generation for the scikit learn + +Generate the rst files for the examples by iterating over the python +example files. + +Files that generate images should start with 'plot' + +""" +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/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/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/examples/plot_lena_tv_denoise.py b/examples/plot_lena_tv_denoise.py new file mode 100644 index 00000000..c4000d38 --- /dev/null +++ b/examples/plot_lena_tv_denoise.py @@ -0,0 +1,51 @@ +""" +==================================================== +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 scipy +from scipy import ndimage +import matplotlib.pyplot as plt +from scikits.image.filter import tv_denoise + +l = scipy.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) + From 4174d5fef00e5f1b0cbf2c99c44954e12dd1b0e4 Mon Sep 17 00:00:00 2001 From: Emmanuelle Gouillart Date: Sat, 24 Sep 2011 15:05:16 +0200 Subject: [PATCH 06/14] A few fixes: give credit to scikit-learn guys, fix examples... --- CONTRIBUTORS.txt | 6 ++- doc/examples/plot_lena_tv_denoise.py | 4 +- doc/ext/LICENSE.txt | 43 +++++++++++++++ doc/ext/gen_rst.py | 6 ++- .../auto_examples/images/blank_image.png | Bin 0 -> 1508 bytes examples/plot_lena_tv_denoise.py | 51 ------------------ 6 files changed, 54 insertions(+), 56 deletions(-) create mode 100644 doc/source/auto_examples/images/blank_image.png delete mode 100644 examples/plot_lena_tv_denoise.py 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/plot_lena_tv_denoise.py b/doc/examples/plot_lena_tv_denoise.py index c4000d38..4797aaed 100644 --- a/doc/examples/plot_lena_tv_denoise.py +++ b/doc/examples/plot_lena_tv_denoise.py @@ -21,7 +21,7 @@ from scipy import ndimage import matplotlib.pyplot as plt from scikits.image.filter import tv_denoise -l = scipy.lena() +l = scipy.misc.lena() l = l[230:290, 220:320] noisy = l + 0.4*l.std()*np.random.random(l.shape) @@ -48,4 +48,4 @@ 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 index cb71a605..0cf8cbd2 100644 --- a/doc/ext/gen_rst.py +++ b/doc/ext/gen_rst.py @@ -1,10 +1,12 @@ """ -Example generation for the scikit learn +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' +Files that generate images should start with 'plot'. + +This code was taken from the scikit-learn. """ import os 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 0000000000000000000000000000000000000000..4913b99b19fe0e6f91c068f27ac80831ea473bf2 GIT binary patch literal 1508 zcmVPx#24YJ`L;(K){{a7>y{D4^000SaNLh0L01ejw01ejxLMWSf00007bV*G`2ipk( z2{tlaE*H!I00m@8L_t(|+U?!HZ(CIu#_{J|w@K^*{UHP$2ryOZKm_7PtVp0r-9Ux^ zqNE6kKLCb`g^>kXrvnl~NbC_iYAHzzLmgPi)Z#Eunl^FeYrAsp%V5WG9AC$doNFKR z{Zy)|*ij?<$&>e-b7M;Y00000000000000000000001y>f!q1_=9Tk2An@PxsR_ON z%1OQUVy5}A@YkHA8tQTel66* z;}a_9vP$b{y{)g#jr`Z2;weO{YRz;DyQWlg zPA&*k$z?S)e_NGYmO#S~`9evH`A53@pt5)2<(w)DSxwC^E4(&Kpy5YLFWz~m>#N1i zcQB{&LROiJw-iML8XmZ}Ue%4&;;xq7$+OAHDa>W{Qfs88veBOV*XLv0Gbgl%}gKD)X|`=sU?$B zEmKnZ!(|#8ek>68ZsWBKAFh6ySp9lV;+(wUX|w zSL5$w8p@G4fyRNX6vCLr;dWKK<`sGb8Y8dU>0y2JFqtuOz_VF;Wi5=WXgNvYwX7<+ zSpp3|T6(2WirK=eU*{1-($xGd0uBEchfc1(jVH-2pCiy1gT-Ov>S=fRCjt#m#*Vcx zcI=Qzl>{$VBG4E>OPAtXx_;FXfrdZwg_3@c-R($r?~XuYAdADh>($t!a5<^6FiW7} z(=1(y-R&?rMG0L8N1!p7KtHH<*}gNSEF#eGsiET!S>$t>=-BaBw`-&cG<+L7I&R-N0u5h79oF^)+LYX)@CZ)<25F;Q(`I@1 z?-M#RIc8{#%fiqaFTqPgV=NYiF2U0JN_-!G%>%#$bs?j$KPS%|0jXQp@_`% z?alj&zB$K(0soypdz=RW}C_8=d0W~QMr z63uO9U0d1AdI~v#Mk0ZY3uC8`ZwGw8j}d4LeYEK1_FmhVKx3#QMW+yVYT|=@QxIq* z9_n6$+i7dB*Hc6Gp^->C_p<-oDd^juOG9I1+qs_^xl@?a=`;<%YYl8Fm5Occ=6!hBk@U!M(EIi>Pt(vyJjih)MBE`uw|SZw z8Y3I(UKWu11)7FNqT&J>1i`*X-o0$y(a=bwp Date: Sat, 24 Sep 2011 13:21:44 -0700 Subject: [PATCH 07/14] ENH: Expose data module. Currently just provides scipy's lena in a more convenient location. --- scikits/image/data/__init__.py | 4 ++++ scikits/image/setup.py | 1 + 2 files changed, 5 insertions(+) create mode 100644 scikits/image/data/__init__.py 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': From 707c4c801a0c35a1503575a6bd8c82fed6c589b6 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sat, 24 Sep 2011 14:45:33 -0700 Subject: [PATCH 08/14] DOC: Update tv example to use data module. Rewrap some text. --- doc/examples/plot_lena_tv_denoise.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/doc/examples/plot_lena_tv_denoise.py b/doc/examples/plot_lena_tv_denoise.py index 4797aaed..55f167e1 100644 --- a/doc/examples/plot_lena_tv_denoise.py +++ b/doc/examples/plot_lena_tv_denoise.py @@ -3,25 +3,26 @@ 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. +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. +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 scipy -from scipy import ndimage import matplotlib.pyplot as plt + +from scikits.image import data from scikits.image.filter import tv_denoise -l = scipy.misc.lena() +l = data.lena() l = l[230:290, 220:320] noisy = l + 0.4*l.std()*np.random.random(l.shape) @@ -46,6 +47,6 @@ 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.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0, left=0, + right=1) plt.show() From 9fb46e1e441e375497cd5c8c9148ef95e52705e2 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 25 Sep 2011 15:03:01 -0700 Subject: [PATCH 09/14] Add Cython homography implementation. Not well optimised yet. --- scikits/image/transform/_project.pyx | 150 +++++++++++++++++++++++++++ scikits/image/transform/setup.py | 4 + 2 files changed, 154 insertions(+) create mode 100644 scikits/image/transform/_project.pyx diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx new file mode 100644 index 00000000..5cbbb608 --- /dev/null +++ b/scikits/image/transform/_project.pyx @@ -0,0 +1,150 @@ +#cython: cdivison=True boundscheck=False + +cimport cython +cimport numpy as np + +import numpy as np +import cython + +np.import_array() + +cdef extern from "math.h": + double floor(double) + double fmod(double, double) + +cdef double get_pixel(np.ndarray image, int r, int c, char mode, + double cval=0): + cdef np.ndarray[dtype=np.double_t, ndim=2] img = image + cdef int rows = img.shape[0] + cdef int cols = img.shape[1] + + if mode == 'C': + if (r < 0) or (r >= cols) or (c < 0) or (c >= cols): + return cval + else: + return img[r, c] + else: + return img[coord_map(rows, r, mode), + coord_map(cols, c, mode)] + +cdef int coord_map(int dim, int coord, char mode): + dim = dim - 1 + if mode == 'M': # mirror + if (coord < 0): + return (-coord % dim) + elif (coord > dim): + return (dim - (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, H): + cdef np.ndarray[np.double_t, ndim=2] M = H + cdef double xx, yy, zz + + xx = M[0, 0] * x + M[0, 1] * y + M[0, 2] + yy = M[1, 0] * x + M[1, 1] * y + M[1, 2] + zz = M[2, 0] * x + M[2, 1] * y + M[2, 2] + + xx /= zz + yy /= zz + + return xx, yy + +@cython.boundscheck(False) +def homography(np.ndarray image, np.ndarray H, output_shape=None, + mode='C', 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 : {'C', 'M', 'W'} + How to handle values outside the image borders. + Constant, Mirror or Wrap. + 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 = image + cdef np.ndarray[dtype=np.double_t, ndim=2] M = np.linalg.inv(H) + + if mode not in ('C', 'W', 'M'): + raise ValueError("Invalid mode specified. Please use " + "C [constant], W [wrap] or M [mirror].") + cdef char mode_c = ord(mode[0]) + + 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.float64) + + 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): + c, r = tf(tfc, tfr, M) + r_int = floor(r) + c_int = floor(c) + + t = r - r_int + u = c - c_int + + y0 = get_pixel(img, r_int, c_int, mode_c) + y1 = get_pixel(img, r_int + 1, c_int, mode_c) + y2 = get_pixel(img, r_int + 1, c_int + 1, mode_c) + y3 = get_pixel(img, 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/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__': From 167b43fc33a6fb1f5cf25cd3ac1f2707c6073247 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 25 Sep 2011 16:45:52 -0700 Subject: [PATCH 10/14] ENH: fast_homography: 20x speed improvement. --- scikits/image/transform/_project.pyx | 71 ++++++++++++++++------------ 1 file changed, 42 insertions(+), 29 deletions(-) diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx index 5cbbb608..06bc53cb 100644 --- a/scikits/image/transform/_project.pyx +++ b/scikits/image/transform/_project.pyx @@ -6,26 +6,24 @@ 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(np.ndarray image, int r, int c, char mode, - double cval=0): - cdef np.ndarray[dtype=np.double_t, ndim=2] img = image - cdef int rows = img.shape[0] - cdef int cols = img.shape[1] - +cdef double get_pixel(double *image, int rows, int cols, + int r, int c, char mode, double cval=0): if mode == 'C': if (r < 0) or (r >= cols) or (c < 0) or (c >= cols): return cval else: - return img[r, c] + return image[r * cols + c] else: - return img[coord_map(rows, r, mode), - coord_map(cols, c, mode)] + return image[coord_map(rows, r, mode) * cols + + coord_map(cols, c, mode)] cdef int coord_map(int dim, int coord, char mode): dim = dim - 1 @@ -42,22 +40,28 @@ cdef int coord_map(int dim, int coord, char mode): return coord -cdef tf(double x, double y, H): - cdef np.ndarray[np.double_t, ndim=2] M = H +cdef tf(double x, double y, double* H, double *x_, double *y_): cdef double xx, yy, zz - xx = M[0, 0] * x + M[0, 1] * y + M[0, 2] - yy = M[1, 0] * x + M[1, 1] * y + M[1, 2] - zz = M[2, 0] * x + M[2, 1] * y + M[2, 2] + ## print + ## print H[0], H[1], H[2] + ## print H[3], H[4], H[5] + ## print H[6], H[7], H[8] + ## print - xx /= zz - 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] - return xx, yy + 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='C', double cval=0): + mode='constant', double cval=0): """ Projective transformation (homography). @@ -93,9 +97,8 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, Shape of the output image generated. order : int Order of splines used in interpolation. - mode : {'C', 'M', 'W'} + mode : {'constant', 'mirror', 'wrap'} How to handle values outside the image borders. - Constant, Mirror or Wrap. cval : string Used in conjunction with mode 'C' (constant), the value outside the image boundaries. @@ -103,12 +106,18 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, """ cdef np.ndarray[dtype=np.double_t, ndim=2] img = image - cdef np.ndarray[dtype=np.double_t, ndim=2] M = np.linalg.inv(H) + cdef np.ndarray[dtype=np.double_t, ndim=2] M = \ + np.ascontiguousarray(np.linalg.inv(H)) - if mode not in ('C', 'W', 'M'): + if mode not in ('constant', 'wrap', 'mirror'): raise ValueError("Invalid mode specified. Please use " - "C [constant], W [wrap] or M [mirror].") - cdef char mode_c = ord(mode[0]) + "`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: @@ -130,17 +139,21 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, for tfr in range(out_r): for tfc in range(out_c): - c, r = tf(tfc, tfr, M) + 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, r_int, c_int, mode_c) - y1 = get_pixel(img, r_int + 1, c_int, mode_c) - y2 = get_pixel(img, r_int + 1, c_int + 1, mode_c) - y3 = get_pixel(img, r_int, c_int + 1, mode_c) + 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 + \ From 69ae2f022b4a1f8ae48301e68689094b10b06954 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 25 Sep 2011 16:53:55 -0700 Subject: [PATCH 11/14] BUG: fast_homography: Fix mirror mode. --- scikits/image/transform/_project.pyx | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx index 06bc53cb..e2d08a02 100644 --- a/scikits/image/transform/_project.pyx +++ b/scikits/image/transform/_project.pyx @@ -17,7 +17,7 @@ cdef extern from "math.h": cdef double get_pixel(double *image, int rows, int cols, int r, int c, char mode, double cval=0): if mode == 'C': - if (r < 0) or (r >= cols) or (c < 0) or (c >= cols): + if (r < 0) or (r > cols - 1) or (c < 0) or (c > cols - 1): return cval else: return image[r * cols + c] @@ -31,7 +31,10 @@ cdef int coord_map(int dim, int coord, char mode): if (coord < 0): return (-coord % dim) elif (coord > dim): - return (dim - (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)) From 2d94a273f6c77aa0c2d0f70092fd8966f141b068 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 25 Sep 2011 17:22:34 -0700 Subject: [PATCH 12/14] Fix mirror mode. Update docstrings. Update tests --- scikits/image/transform/_project.pyx | 61 ++++++++++++++++--- scikits/image/transform/tests/test_project.py | 36 +++++++++++ 2 files changed, 88 insertions(+), 9 deletions(-) diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx index e2d08a02..d1788ab8 100644 --- a/scikits/image/transform/_project.pyx +++ b/scikits/image/transform/_project.pyx @@ -1,5 +1,7 @@ #cython: cdivison=True boundscheck=False +__all__ = ['homography'] + cimport cython cimport numpy as np @@ -16,6 +18,22 @@ cdef extern from "math.h": 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 > cols - 1) or (c < 0) or (c > cols - 1): return cval @@ -26,10 +44,28 @@ cdef double get_pixel(double *image, int rows, int 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): - return (-coord % dim) + # 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)) @@ -44,13 +80,19 @@ cdef int coord_map(int dim, int coord, char mode): return coord cdef tf(double x, double y, double* H, double *x_, double *y_): - cdef double xx, yy, zz + """Apply a homography to a coordinate. - ## print - ## print H[0], H[1], H[2] - ## print H[3], H[4], H[5] - ## print H[6], H[7], H[8] - ## print + 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] @@ -108,7 +150,8 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, """ - cdef np.ndarray[dtype=np.double_t, ndim=2] img = image + 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)) @@ -134,7 +177,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, columns = img.shape[1] cdef np.ndarray[dtype=np.double_t, ndim=2] out = \ - np.zeros((out_r, out_c), dtype=np.float64) + np.zeros((out_r, out_c), dtype=np.double) cdef int tfr, tfc, r_int, c_int cdef double y0, y1, y2, y3 diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py index 44c0e557..7ae24049 100644 --- a/scikits/image/transform/tests/test_project.py +++ b/scikits/image/transform/tests/test_project.py @@ -2,6 +2,8 @@ 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 homography as fast_homography +from scikits.image import data def test_stackcopy(): layers = 4 @@ -19,3 +21,37 @@ 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) + p1 = fast_homography(img, H, mode=mode) + p1 = np.round(p1) + + ## import matplotlib.pyplot as plt + ## plt.imshow(np.abs(p0 - p1), cmap=plt.cm.gray) + ## plt.show() + + d = np.mean(np.abs(p0 - p1)) + assert d < 0.1 + + +if __name__ == "__main__": + from numpy.testing import run_module_suite + run_module_suite() From 2b71fd9a0eebb038c41fe61432bc4019dcbe52e6 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 25 Sep 2011 17:29:38 -0700 Subject: [PATCH 13/14] BUG: fast_homography: Fix constant mode. --- scikits/image/transform/_project.pyx | 2 +- scikits/image/transform/tests/test_project.py | 10 +++++++--- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx index d1788ab8..3e2b9b7d 100644 --- a/scikits/image/transform/_project.pyx +++ b/scikits/image/transform/_project.pyx @@ -35,7 +35,7 @@ cdef double get_pixel(double *image, int rows, int cols, """ if mode == 'C': - if (r < 0) or (r > cols - 1) or (c < 0) or (c > cols - 1): + if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): return cval else: return image[r * cols + c] diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py index 7ae24049..b7c25fdf 100644 --- a/scikits/image/transform/tests/test_project.py +++ b/scikits/image/transform/tests/test_project.py @@ -40,16 +40,20 @@ def test_fast_homography(): for mode in ('constant', 'mirror', 'wrap'): print 'Transform mode:', mode - p0 = homography(img, H, 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 - ## plt.imshow(np.abs(p0 - p1), cmap=plt.cm.gray) + ## 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.1 + assert d < 0.2 if __name__ == "__main__": From 3df7df263ea02d3589a26a395c955c33798bb2b8 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Sun, 25 Sep 2011 17:37:04 -0700 Subject: [PATCH 14/14] Import fast homography as transform.fast_homography. --- scikits/image/transform/__init__.py | 1 + scikits/image/transform/tests/test_project.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py index 341671d1..e13d3ac8 100644 --- a/scikits/image/transform/__init__.py +++ b/scikits/image/transform/__init__.py @@ -1,4 +1,5 @@ from .hough_transform 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/tests/test_project.py b/scikits/image/transform/tests/test_project.py index b7c25fdf..ba8e277e 100644 --- a/scikits/image/transform/tests/test_project.py +++ b/scikits/image/transform/tests/test_project.py @@ -1,8 +1,8 @@ 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 homography as fast_homography +from scikits.image.transform.project import _stackcopy +from scikits.image.transform import homography, fast_homography from scikits.image import data def test_stackcopy():