diff --git a/bento.info b/bento.info index cba403d9..7d21bfa8 100644 --- a/bento.info +++ b/bento.info @@ -161,6 +161,9 @@ Library: Extension: skimage.external.tifffile._tifffile Sources: skimage/external/tifffile/_tifffile.c + Extension: skimage.transform._seam_carving + Sources: + skimage/transform/_seam_carving.pyx Executable: skivi Module: skimage.scripts.skivi diff --git a/doc/examples/plot_seam_carving.py b/doc/examples/plot_seam_carving.py new file mode 100644 index 00000000..f1dcb948 --- /dev/null +++ b/doc/examples/plot_seam_carving.py @@ -0,0 +1,95 @@ +""" +============ +Seam Carving +============ + +This example demonstrates how images can be resized using seam carving [1]_. +Resizing to a new aspect ratio distorts image contents. Seam carving attempts +to resize *without* distortion, by removing regions of an image which are less +important. In this example we are using the Sobel filter to signify the +importance of each pixel. + +.. [1] Shai Avidan and Ariel Shamir + "Seam Carving for Content-Aware Image Resizing" + http://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Avidan07.pdf + +""" +from skimage import data, draw +from skimage import transform, util +import numpy as np +from skimage import filters, color +from matplotlib import pyplot as plt + + +hl_color = np.array([0, 1, 0]) + +img = data.rocket() +img = util.img_as_float(img) +eimg = filters.sobel(color.rgb2gray(img)) + +plt.title('Original Image') +plt.imshow(img) + +""" +.. image:: PLOT2RST.current_figure +""" + +resized = transform.resize(img, (img.shape[0], img.shape[1] - 200)) +plt.figure() +plt.title('Resized Image') +plt.imshow(resized) + + +""" +.. image:: PLOT2RST.current_figure +""" + +out = transform.seam_carve(img, eimg, 'vertical', 200) +plt.figure() +plt.title('Resized using Seam-Carving') +plt.imshow(out) + +""" +.. image:: PLOT2RST.current_figure + +As you can see, resizing as distorted the rocket and the objects around, +whereas seam carving has reszied by removing the empty spaces in between. + +Object Removal +-------------- + +Seam Carving can also be used to remove atrifacts from images. To do that, we +have to ensure that pixels to be removes get less importance. In the following +code I approximately mark the rocket with a mask, and then decrease the +importance of those pixels + +""" + +masked_img = img.copy() + +poly = [(404, 281), (404, 360), (359, 364), (338, 337), (145, 337), (120, 322), + (145, 304), (340, 306), (362, 284)] +pr = np.array([p[0] for p in poly]) +pc = np.array([p[1] for p in poly]) +rr, cc = draw.polygon(pr, pc) + +masked_img[rr, cc, :] = masked_img[rr, cc, :]*0.5 + hl_color*.5 +plt.figure() +plt.title('Object Marked') + +plt.imshow(masked_img) +""" +.. image:: PLOT2RST.current_figure +""" + +eimg[rr, cc] -= 1000 + +plt.figure() +plt.title('Object Removed') +out = transform.seam_carve(img, eimg, 'vertical', 90) +resized = transform.resize(img, out.shape) +plt.imshow(out) +plt.show() +""" +.. image:: PLOT2RST.current_figure +""" diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py index 1711185c..a5f5068d 100644 --- a/skimage/data/__init__.py +++ b/skimage/data/__init__.py @@ -26,6 +26,7 @@ __all__ = ['load', 'chelsea', 'coffee', 'hubble_deep_field', + 'rocket', 'astronaut'] @@ -241,3 +242,22 @@ def hubble_deep_field(): """ return load("hubble_deep_field.jpg") + + +def rocket(): + """Launch photo of DSCOVR on Falcon 9 by SpaceX. + + This is the launch photo of Falcon 9 carrying DSCOVR lifted off from + SpaceX's Launch Complex 40 at Cape Canaveral Air Force Station, FL. + + Notes + ----- + This image was downloaded from + `SpaceX Photos + `__. + + The image was captured by SpaceX and `released in the public domain + `_. + + """ + return load("rocket.jpg") diff --git a/skimage/data/rocket.jpg b/skimage/data/rocket.jpg new file mode 100644 index 00000000..32bdd339 Binary files /dev/null and b/skimage/data/rocket.jpg differ diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 3049ebb0..9b40b3a9 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -12,6 +12,7 @@ from ._geometric import (warp, warp_coords, estimate_transform, from ._warps import swirl, resize, rotate, rescale, downscale_local_mean from .pyramids import (pyramid_reduce, pyramid_expand, pyramid_gaussian, pyramid_laplacian) +from .seam_carving import seam_carve __all__ = ['hough_circle', @@ -43,4 +44,5 @@ __all__ = ['hough_circle', 'pyramid_reduce', 'pyramid_expand', 'pyramid_gaussian', - 'pyramid_laplacian'] + 'pyramid_laplacian', + 'seam_carve'] diff --git a/skimage/transform/_seam_carving.pyx b/skimage/transform/_seam_carving.pyx new file mode 100644 index 00000000..3499514c --- /dev/null +++ b/skimage/transform/_seam_carving.pyx @@ -0,0 +1,226 @@ +# cython: cdivision=True +# cython: boundscheck=False +# cython: nonecheck=False +# cython: wraparound=False +import numpy as np +cimport numpy as cnp + + +cdef cnp.double_t DBL_MAX = np.finfo(np.double).max + +cdef void _preprocess_image(cnp.double_t[:, :, ::1] energy_img, + cnp.double_t[:, ::1] cumulative_img, + cnp.int8_t[:, ::1] track_img, + Py_ssize_t cols) nogil: + """ For each row, compute the lowest seam value for all its columns. + + This function updates `cumulative_img` such that `cumulative_img[r, c]` + is the total energy of the lowest energy seam ending at `(r, c)`. + + Parameters + ---------- + energy_img : (M, N, 1) ndarray + Cost array representing the expense to remove each pixel. Seam carving + tries to avoid pixels with high costs. + cumulative_img : (M, N) ndarray + The array to be updated inplace with the total cost of lowest energy + seams. + track_img : (M, N) ndarray + For each pixel, `track_img` stores the relative column offset in + the previous row which has the lowest value in `cumulative_img`. This + helps in in re-tracing the minimum cost seam. + cols : int + Number of columns to process. + """ + + cdef Py_ssize_t r, c, offset, c_idx + cdef Py_ssize_t rows = energy_img.shape[0] + cdef cnp.double_t min_cost = DBL_MAX + cdef Py_ssize_t colsm1 = cols - 1 + cdef Py_ssize_t rm1 + + for c in range(cols): + cumulative_img[0, c] = energy_img[0, c, 0] + + for r in range(1, rows): + rm1 = r - 1 + for c in range(cols): + min_cost = DBL_MAX + for offset in range(-1, 2): + + c_idx = c + offset + if (c_idx > colsm1) or (c_idx < 0): + continue + + if cumulative_img[rm1, c_idx] < min_cost: + min_cost = cumulative_img[rm1, c_idx] + track_img[r, c] = offset + + cumulative_img[r, c] = min_cost + energy_img[r, c, 0] + +cdef bint _mark_seam(cnp.int8_t[:, ::1] track_img, + Py_ssize_t start_index, + cnp.uint8_t[:, ::1] seam_map, + Py_ssize_t[::1] seam_buffer) nogil: + + """ Re-trace the optimal seam from a given column in the last row. + + This function tries to re-track an optimal seam from `start_index` and + tries to mark it in `seam_map`. If this seam intersects with any existing + seam in `seam_map` the function returns `0` without marking anything. Else + it marks the seam in `seam_map` and returns `1`. + + track_img : (M, N) ndarray + The array of relative column indices as updated by `_preprocess_image`. + start_index : int + The column number of the bottom most row from where to start re-tracing + the seam. + seam_map : (M, N) ndarray + The array used to mark seams. If a pixel is marked as as seam it is set + to `1`, else `0`. + seam_buffer : (M,) ndarray + Buffer used to store the column indices of the seam currently being + checked. This is preallocated to save time. + + Returns + ------- + success : int + `1` if seam was marked, `0` is seam intersects and was not marked. + """ + cdef Py_ssize_t rows = track_img.shape[0] + cdef Py_ssize_t[::1] current_seam_indices = seam_buffer + cdef Py_ssize_t row, col + cdef cnp.int8_t offset + cdef Py_ssize_t seams + + current_seam_indices[rows - 1] = start_index + for row in range(rows - 2, -1, -1): + col = current_seam_indices[row + 1] + offset = track_img[row, col] + col = col + offset + current_seam_indices[row] = col + + if seam_map[row, col]: + return 0 + + for row in range(rows): + col = current_seam_indices[row] + seam_map[row, col] = 1 + + return 1 + +cdef void _remove_seam(cnp.double_t[:, :, ::1] img, + cnp.uint8_t[:, ::1] seam_map, Py_ssize_t cols) nogil: + """ Remove marked seams from an image. + + Parameters + ---------- + img : (M, N, P) ndarray + Input image whose vertical seams are to be removed. + seam_map : (M, N) ndarray + Array with seams to be removed marked by non-zero entries. + cols : int + The number of columns to process. + """ + cdef Py_ssize_t rows = img.shape[0] + cdef Py_ssize_t channels = img.shape[2] + cdef Py_ssize_t r, c, ch, shift + cdef Py_ssize_t c_shift + + for r in range(rows): + shift = 0 + for c in range(cols): + shift += seam_map[r, c] + c_shift = c + shift + for ch in range(channels): + img[r, c, ch] = img[r, c_shift, ch] + + +def _seam_carve_v(img, energy_map, iters, border): + """ Carve vertical seams off an image. + + Carves out vertical seams from an image while using the given energy map to + decide the importance of each pixel.[1]_ + + Parameters + ---------- + img : (M, N) or (M, N, 3) ndarray + Input image whose vertical seams are to be removed. + energy_map : (M, N) ndarray + Cost array denoting importance of each pixel. The algorithm will try to + retain high valued pixels. + iters : int + Number of vertical seams to be removed. + border : int, optional + The number of pixels in the right, left and bottom end of the image + to be excluded from being considered for a seam. This is important as + certain filters just ignore image boundaries and set them to `0`. + By default border is set to `1`. + + Returns + ------- + image : (M, N - iters, 3) ndarray of float + The cropped image with the vertical seams removed. + + References + ---------- + .. [1] Shai Avidan and Ariel Shamir + "Seam Carving for Content-Aware Image Resizing" + http://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Avidan07.pdf + """ + # This reference has been kept to be used for the `np.argsort` call + last_row_obj = np.zeros(img.shape[1], dtype=np.float) + + cdef cnp.double_t[::1] last_row = last_row_obj + cdef Py_ssize_t[::1] sorted_indices + cdef cnp.uint8_t[:, ::1] seam_map = np.zeros(img.shape[0:2], + dtype=np.uint8) + cdef Py_ssize_t cols = img.shape[1] + cdef Py_ssize_t rows = img.shape[0] + cdef Py_ssize_t seams_left = iters + cdef Py_ssize_t seams_removed + cdef Py_ssize_t seam_idx + cdef Py_ssize_t[::1] seam_buffer = np.zeros(rows, dtype=np.int) + + cdef cnp.double_t[:, :, ::1] image = img + cdef cnp.int8_t[:, ::1] track_img = np.zeros(img.shape[0:2], dtype=np.int8) + cdef cnp.double_t[:, ::1] cumulative_img = np.zeros(img.shape[0:2], + dtype=np.float) + cdef cnp.double_t[:, :, ::1] energy_img + + energy_map[:, 0:border] = DBL_MAX + energy_map[:, cols-border:cols] = DBL_MAX + + # Filters often let the boundary be `0`. If all the entries in the last + # row of `energy_img` are equal, the minimum value in the penultimate row + # of `cumulative_img` will result in 3 minimum values in its last row. + # Hence, two successive removals will always intersect as the 3 least seams + # will share the same pixels except they will differ in the last row. + energy_map[rows-border:rows, :] = energy_map[rows-2*border:rows-border, :] + + energy_map = np.ascontiguousarray(energy_map[:, :, np.newaxis]) + energy_img = energy_map + + _preprocess_image(energy_img, cumulative_img, track_img, cols) + last_row[...] = cumulative_img[rows - 1, :] + sorted_indices = np.argsort(last_row_obj) + seam_idx = 0 + + while seams_left > 0: + if _mark_seam(track_img, sorted_indices[seam_idx], seam_map, + seam_buffer): + seams_left -= 1 + cols -= 1 + seam_idx += 1 + else: + seam_idx = 0 + _remove_seam(image, seam_map, cols) + _remove_seam(energy_img, seam_map, cols) + seam_map[...] = 0 + _preprocess_image(energy_img, cumulative_img, track_img, cols) + last_row[:cols] = cumulative_img[rows - 1, :cols] + sorted_indices = np.argsort(last_row_obj) + + _remove_seam(image, seam_map, cols) + + return img[:, 0:cols] diff --git a/skimage/transform/seam_carving.py b/skimage/transform/seam_carving.py new file mode 100644 index 00000000..95e975ce --- /dev/null +++ b/skimage/transform/seam_carving.py @@ -0,0 +1,66 @@ +from ._seam_carving import _seam_carve_v +from .. import util +from .._shared import utils +import numpy as np + + +def seam_carve(img, energy_map, mode, num, border=1, force_copy=True): + """ Carve vertical or horizontal seams off an image. + + Carves out vertical/horizontal seams from an image while using the given + energy map to decide the importance of each pixel. + + Parameters + ---------- + image : (M, N) or (M, N, 3) ndarray + Input image whose seams are to be removed. + energy_map : (M, N) ndarray + The array to decide the importance of each pixel. The higher + the value corresponding to a pixel, the more the algorithm will try + to keep it in the image. + mode : str {'horizontal', 'vertical'} + Indicates whether seams are to be removed vertically or horizontally. + Removing seams horizontally will decrease the height whereas removing + vertically will decrease the width. + num : int + Number of seams are to be removed. + border : int, optional + The number of pixels in the right, left and bottom end of the image + to be excluded from being considered for a seam. This is important as + certain filters just ignore image boundaries and set them to `0`. + By default border is set to `1`. + force_copy : bool, optional + If set, the `image` and `energy_map` are copied before being used by + the method which modifies it in place. Set this to `False` if the + original image and the energy map are no longer needed after + this opetration. + + Returns + ------- + out : ndarray + The cropped image with the seams removed. + + References + ---------- + .. [1] Shai Avidan and Ariel Shamir + "Seam Carving for Content-Aware Image Resizing" + http://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Avidan07.pdf + """ + + utils.assert_nD(img, (2, 3)) + image = util.img_as_float(img, force_copy) + energy_map = util.img_as_float(energy_map, force_copy) + + if image.ndim == 2: + image = image[..., np.newaxis] + + if mode == 'horizontal': + image = np.transpose(image, (1, 0, 2)) + + image = np.ascontiguousarray(image) + out = _seam_carve_v(image, energy_map, num, border) + + if mode == 'horizontal': + out = np.transpose(out, (1, 0, 2)) + + return np.squeeze(out) diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py index 22f31696..ff2e9bc6 100644 --- a/skimage/transform/setup.py +++ b/skimage/transform/setup.py @@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None): cython(['_hough_transform.pyx'], working_path=base_path) cython(['_warps_cy.pyx'], working_path=base_path) cython(['_radon_transform.pyx'], working_path=base_path) + cython(['_seam_carving.pyx'], working_path=base_path) config.add_extension('_hough_transform', sources=['_hough_transform.c'], include_dirs=[get_numpy_include_dirs()]) @@ -27,6 +28,8 @@ def configuration(parent_package='', top_path=None): sources=['_radon_transform.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_seam_carving', sources=['_seam_carving.c'], + include_dirs=[get_numpy_include_dirs()]) return config if __name__ == '__main__': diff --git a/skimage/transform/tests/test_seam_carving.py b/skimage/transform/tests/test_seam_carving.py new file mode 100644 index 00000000..b4eb8977 --- /dev/null +++ b/skimage/transform/tests/test_seam_carving.py @@ -0,0 +1,23 @@ +from skimage import transform +import numpy as np +from numpy import testing + + +def test_seam_carving(): + img = np.array([[0, 0, 1, 0, 0], + [0, 0, 1, 0, 0], + [0, 0, 1, 0, 0], + [0, 1, 0, 0, 0], + [1, 0, 0, 0, 0]], dtype=np.float) + energy = 1 - img + + out = transform.seam_carve(img, energy, 'vertical', 1, border=0) + testing.assert_allclose(out, 0) + + img = img.T + out = transform.seam_carve(img, energy, 'horizontal', 1, border=0) + testing.assert_allclose(out, 0) + + +if __name__ == '__main__': + np.testing.run_module_suite()