diff --git a/.gitignore b/.gitignore index 92e34e201..31fa17aa7 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,6 @@ data/berlin/* !data/berlin/images eval + +# Ignore virtualenv files +env/* diff --git a/README.md b/README.md index 9207ac607..5d320de84 100644 --- a/README.md +++ b/README.md @@ -30,7 +30,6 @@ Install OpenCV using brew install opencv brew install homebrew/science/ceres-solver brew install boost-python - sudo pip install -r requirements.txt And install OpenGV using @@ -63,6 +62,10 @@ When running OpenSfM on top of [OpenCV][] 3.0 the [OpenCV Contrib][] modules are ## Building + sudo pip install virtualenv + virtualenv env + source env/bin/activate + pip install -r requirements.txt python setup.py build @@ -72,7 +75,7 @@ An example dataset is available at `data/berlin`. 1. Put some images in `data/DATASET_NAME/images/` 2. Put config.yaml in `data/DATASET_NAME/config.yaml` - 3. Go to the root of the project and run `bin/run_all data/DATASET_NAME` + 3. Go to the root of the project and run `bin/opensfm_run_all data/DATASET_NAME` 4. Start an http server from the root with `python -m SimpleHTTPServer` 5. Browse `http://localhost:8000/viewer/reconstruction.html#file=/data/DATASET_NAME/reconstruction.meshed.json`. diff --git a/bin/run_all b/bin/opensfm_run_all similarity index 100% rename from bin/run_all rename to bin/opensfm_run_all diff --git a/bin/plot_depthmaps b/bin/plot_depthmaps new file mode 100644 index 000000000..965e5cbef --- /dev/null +++ b/bin/plot_depthmaps @@ -0,0 +1,76 @@ +#!/usr/bin/env python + +import argparse +import matplotlib.pyplot as pl +import numpy as np +from textwrap import wrap + +from opensfm import dataset +from opensfm import io + + +def plot_depthmap(im, title, depth, plane, score, nghbr): + ax = pl.subplot2grid((2, 3), (0, 0), rowspan=2) + ax_title = ax.set_title(title) + ax_title.set_y(1.05) + pl.imshow(im) + ax = pl.subplot(2, 3, 2) + ax_title = ax.set_title("Depth") + ax_title.set_y(1.08) + pl.imshow(depth) + pl.colorbar() + ax = pl.subplot(2, 3, 3) + ax_title = ax.set_title("Score") + ax_title.set_y(1.08) + pl.imshow(score) + pl.colorbar() + ax = pl.subplot(2, 3, 5) + ax_title = ax.set_title("Neighbor") + ax_title.set_y(1.08) + pl.imshow(nghbr) + pl.colorbar() + ax = pl.subplot(2, 3, 6) + ax_title = ax.set_title("Plane normal") + ax_title.set_y(1.02) + pl.imshow(color_plane_normals(plane)) + + +def color_plane_normals(plane): + l = np.linalg.norm(plane, axis=2) + normal = plane / l[..., np.newaxis] + normal[..., 1] *= -1 # Reverse Y because it points down + normal[..., 2] *= -1 # Reverse Z because standard colormap does so + return ((normal + 1) * 128).astype(np.uint8) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description='Compute reconstruction') + parser.add_argument('dataset', + help='path to the dataset to be processed') + parser.add_argument('--image', + help='name of the image to show') + parser.add_argument('--save-figs', + help='save figures istead of showing them', + action='store_true') + args = parser.parse_args() + + data = dataset.DataSet(args.dataset) + + images = [args.image] if args.image else data.images() + for image in images: + depth, plane, score, nghbr, nghbrs = data.load_raw_depthmap(image) + + print "Plotting depthmap for {0}".format(image) + pl.close("all") + pl.figure(figsize=(30, 16), dpi=90, facecolor='w', edgecolor='k') + title = "Shot: " + image + "\n" + "\n".join(wrap("Neighbors: " + ', '.join(nghbrs), 80)) + plot_depthmap(data.image_as_array(image), title, depth, plane, score, nghbr) + pl.tight_layout() + + if args.save_figs: + p = args.dataset + '/plot_depthmaps' + io.mkdir_p(p) + pl.savefig(p + '/' + image + '.png') + pl.close() + else: + pl.show() diff --git a/bin/plot_matches b/bin/plot_matches index eebeb161e..776bf39f6 100755 --- a/bin/plot_matches +++ b/bin/plot_matches @@ -1,6 +1,7 @@ #!/usr/bin/env python import argparse +from itertools import combinations import matplotlib.pyplot as pl import networkx as nx @@ -33,7 +34,9 @@ if __name__ == "__main__": parser.add_argument('dataset', help='path to the dataset to be processed') parser.add_argument('--image', - help='name of the image to show') + help='show tracks for a specific') + parser.add_argument('--images', + help='show tracks between a subset of images (separated by commas)') parser.add_argument('--graph', help='display image graph', action='store_true') @@ -65,39 +68,40 @@ if __name__ == "__main__": else: # Plot matches between images if args.image: - toplot = [args.image] + pairs = [(args.image, o) for o in images if o != args.image] + elif args.images: + subset = args.images.split(',') + pairs = combinations(subset, 2) else: - toplot = images + pairs = combinations(images, 2) i = 0 - for im1 in toplot: - for im2 in images: - if im1 != im2: - matches = data.find_matches(im1, im2) - if len(matches) == 0: - continue - print 'plotting matches between', im1, im2 - - p1, f1, c1 = data.load_features(im1) - p2, f2, c2 = data.load_features(im2) - p1 = p1[matches[:, 0]] - p2 = p2[matches[:, 1]] - - pl.figure(figsize=(20, 10)) - pl.title('Images: ' + im1 + ' - ' + im2 + ', matches: ' + - str(matches.shape[0])) - plot_matches(data.image_as_array(im1), - data.image_as_array(im2), p1, p2) - i += 1 - if args.save_figs: - p = args.dataset + '/plot_tracks' - io.mkdir_p(p) - pl.savefig(p + '/' + im1 + '_' + im2 + '.jpg', dpi=100) - pl.close() - else: - if i >= 10: - i = 0 - pl.show() + for im1, im2 in pairs: + matches = data.find_matches(im1, im2) + if len(matches) == 0: + continue + print 'plotting matches between', im1, im2 + + p1, f1, c1 = data.load_features(im1) + p2, f2, c2 = data.load_features(im2) + p1 = p1[matches[:, 0]] + p2 = p2[matches[:, 1]] + + pl.figure(figsize=(20, 10)) + pl.title('Images: ' + im1 + ' - ' + im2 + ', matches: ' + + str(matches.shape[0])) + plot_matches(data.image_as_array(im1), + data.image_as_array(im2), p1, p2) + i += 1 + if args.save_figs: + p = args.dataset + '/plot_tracks' + io.mkdir_p(p) + pl.savefig(p + '/' + im1 + '_' + im2 + '.jpg', dpi=100) + pl.close() + else: + if i >= 10: + i = 0 + pl.show() if not args.save_figs and i > 0: pl.show() diff --git a/bin/plot_tracks b/bin/plot_tracks index 932b6e00c..cc76017ec 100755 --- a/bin/plot_tracks +++ b/bin/plot_tracks @@ -1,6 +1,8 @@ #!/usr/bin/env python import argparse +from itertools import combinations + import matplotlib.pyplot as pl import networkx as nx import numpy as np @@ -30,11 +32,13 @@ def plot_matches(im1, im2, p1, p2): if __name__ == "__main__": - parser = argparse.ArgumentParser(description='Compute reconstruction') + parser = argparse.ArgumentParser(description='Plot tracks') parser.add_argument('dataset', help='path to the dataset to be processed') parser.add_argument('--image', - help='name of the image to show') + help='show tracks for a specific') + parser.add_argument('--images', + help='show tracks between a subset of images (separated by commas)') parser.add_argument('--graph', help='display image graph', action='store_true') @@ -60,32 +64,33 @@ if __name__ == "__main__": else: # Plot matches between images if args.image: - toplot = [args.image] + pairs = [(args.image, o) for o in images if o != args.image] + elif args.images: + subset = args.images.split(',') + pairs = combinations(subset, 2) else: - toplot = images + pairs = combinations(images, 2) i = 0 - for im1 in toplot: - for im2 in images: - if im1 != im2: - t, p1, p2 = matching.common_tracks(graph, im1, im2) - if len(t) >= 10: - pl.figure(figsize=(20, 10)) - pl.title('Images: ' + im1 + ' - ' + im2 + - ', matches: ' + str(len(t))) - plot_matches(data.image_as_array(im1), - data.image_as_array(im2), p1, p2) - i += 1 - if args.save_figs: - p = args.dataset + '/plot_tracks' - io.mkdir_p(p) - pl.savefig(p + '/' + im1 + '_' + im2 + '.jpg', - dpi=100) - pl.close() - else: - if i >= 10: - i = 0 - pl.show() + for im1, im2 in pairs: + t, p1, p2 = matching.common_tracks(graph, im1, im2) + if len(t) >= 10: + pl.figure(figsize=(20, 10)) + pl.title('Images: ' + im1 + ' - ' + im2 + + ', matches: ' + str(len(t))) + plot_matches(data.image_as_array(im1), + data.image_as_array(im2), p1, p2) + i += 1 + if args.save_figs: + p = args.dataset + '/plot_tracks' + io.mkdir_p(p) + pl.savefig(p + '/' + im1 + '_' + im2 + '.jpg', + dpi=100) + pl.close() + else: + if i >= 10: + i = 0 + pl.show() if not args.save_figs and i > 0: pl.show() diff --git a/bin/run_eval b/bin/run_eval index e7e69b812..2ae07724c 100755 --- a/bin/run_eval +++ b/bin/run_eval @@ -17,7 +17,7 @@ from opensfm import types def run_dataset(run_root, name): folder = run_root + '/' + name - call(['bin/run_all', folder]) + call(['bin/opensfm_run_all', folder]) with open(run_root + '/index.html', 'a') as fout: s = ''' diff --git a/config.yaml b/config.yaml index d732cfe20..62c135604 100644 --- a/config.yaml +++ b/config.yaml @@ -29,10 +29,6 @@ akaze_descriptor_channels: 3 # Number of feature channels (1,2,3) hahog_peak_threshold: 0.00001 hahog_edge_threshold: 10 -# Masks for regions that will be ignored for feature extraction -# List of bounding boxes specified as the ratio to image width and height -# masks: [{top: 0.96, bottom: 1.0, left: 0.0, right: 0.15}, {top: 0.95, bottom: 1.0, left: 0, right: 0.05}] - # Params for general matching lowes_ratio: 0.8 # Ratio test for matches preemptive_lowes_ratio: 0.6 # Ratio test for preemptive matches diff --git a/data/berlin/config.yaml b/data/berlin/config.yaml index 8b092a971..863426711 100644 --- a/data/berlin/config.yaml +++ b/data/berlin/config.yaml @@ -5,3 +5,4 @@ processes: 8 # Number of threads to use depthmap_min_consistent_views: 2 # Min number of views that should reconstruct a point for it to be valid depthmap_save_debug_files: yes # Save debug files with partial reconstruction results +feature_process_size: 1024 diff --git a/data/berlin/masks/01.jpg.png b/data/berlin/masks/01.jpg.png new file mode 100644 index 000000000..59cee8b11 Binary files /dev/null and b/data/berlin/masks/01.jpg.png differ diff --git a/data/berlin/masks/02.jpg.png b/data/berlin/masks/02.jpg.png new file mode 100644 index 000000000..587354dc5 Binary files /dev/null and b/data/berlin/masks/02.jpg.png differ diff --git a/data/berlin/masks/03.jpg.png b/data/berlin/masks/03.jpg.png new file mode 100644 index 000000000..4c27e47ea Binary files /dev/null and b/data/berlin/masks/03.jpg.png differ diff --git a/doc/Makefile b/doc/Makefile index 39db9881e..137f03184 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -58,7 +58,7 @@ html: @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." livehtml: html - sphinx-autobuild -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html + sphinx-autobuild -b html -p 8001 $(ALLSPHINXOPTS) $(BUILDDIR)/html dirhtml: $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml diff --git a/doc/source/_static/mathjax_conf.js b/doc/source/_static/mathjax_conf.js index 95e4dc947..3b20bff7d 100644 --- a/doc/source/_static/mathjax_conf.js +++ b/doc/source/_static/mathjax_conf.js @@ -1,6 +1,6 @@ MathJax.Hub.Config({ "HTML-CSS": { availableFonts: ["TeX"], - scale: 90 + scale: 100 } }); diff --git a/doc/source/building.rst b/doc/source/building.rst index 52a65b9e1..e6a82620f 100644 --- a/doc/source/building.rst +++ b/doc/source/building.rst @@ -70,3 +70,12 @@ When running OpenSfM on top of OpenCV version 3.0 the `OpenCV Contrib`_ modules .. _Ceres solver: http://ceres-solver.org/ .. _Boost Python: http://www.boost.org/ .. _Networkx: https://github.com/networkx/networkx + + +Building the documentation +-------------------------- +To build the documentation and browse it locally use +:: + cd doc + make livehtml +and browse `http://localhost:8001/ `_ diff --git a/doc/source/conf.py b/doc/source/conf.py index dadacf9e1..7bf13f8dd 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -54,7 +54,7 @@ # General information about the project. project = u'OpenSfM' -copyright = u'2016, Mapillary' +copyright = u'2017, Mapillary' author = u'Mapillary' # The version info for the project you're documenting, acts as replacement for @@ -62,9 +62,9 @@ # built documents. # # The short X.Y version. -version = u'0.0' +version = u'0.1' # The full version, including alpha/beta/rc tags. -release = u'0.0.0' +release = u'0.1.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/source/dataset.rst b/doc/source/dataset.rst index 809dc5005..8ef6e74fe 100644 --- a/doc/source/dataset.rst +++ b/doc/source/dataset.rst @@ -4,4 +4,21 @@ Dataset Structure ================= -TODO +:: + + project/ + ├── config.yaml + ├── images/ + ├── masks/ + ├── gcp_list.txt + ├── metadata/ + ├── features/ + ├── matches/ + ├── tracks.tsv + ├── reconstruction.json + ├── reconstruction.meshed.json + ├── undistorted/ + ├── undistorted_tracks.json + ├── undistorted_reconstruction.json + └── depthmaps/ + └── merged.ply diff --git a/doc/source/dense.rst b/doc/source/dense.rst index afc2a4f30..c97a1d213 100644 --- a/doc/source/dense.rst +++ b/doc/source/dense.rst @@ -120,7 +120,7 @@ And the linear approximation around :math:`(x_0, y_0)` is Undistortion ------------ -The dense module assumes that images are taken with perspective projection and no radial distortion. For perspective images, undistorted versions can be generated by taking into account the computed distortion parameters, :math:`k1` and :math: `k2`. +The dense module assumes that images are taken with perspective projection and no radial distortion. For perspective images, undistorted versions can be generated by taking into account the computed distortion parameters, :math:`k1` and :math:`k2`. Equirectangular images (360 panoramas) however can not be unwarped into a single persepective view. We need to generate multiple perspective views to cover the field of view of a panorama. diff --git a/doc/source/gcp.rst b/doc/source/gcp.rst new file mode 100644 index 000000000..f3658890d --- /dev/null +++ b/doc/source/gcp.rst @@ -0,0 +1,49 @@ + +Ground Control Points +--------------------- + +When EXIF data contains GPS location, it is used by OpenSfM to georeference the reconstruction. Additionally, it is possible to use ground control points. + +Ground control points (GCP) are landmarks visible on the images for which the geospatial position (latitude, longitude and altitude) is known. A single GCP can be observed in one or more images. + +OpenSfM uses GCP in two steps of the reconstruction process: alignment and bundle adjustment. In the alignment step, points are used to globaly move the reconstruction so that the observed GCP align with their GPS position. Two or more observations for each GCP are required for it to be used during the aligment step. + +In the bundle adjustment step, GCP observations are used as a constraint to refine the reconstruction. In this step, all ground control points are used. No minimum number of observation is required. + +File format +``````````` +GCPs can be specified by adding a text file named ``gcp_list.txt`` at the root folder of the dataset. The format of the file should be as follows. + +- The first line should contain the name of the projection used for the geo coordinates. + +- The following lines should con should contain the data for each ground control point observation. One per line and in the format:: + + + + Where `` `` are the geospatial coordinates of the GCP and `` `` are the pixel coordinates where the GCP is observed. + + +Supported projections +````````````````````` +The geospatial coordinates can be specified in one the following formats. + +- `WGS84`_: This is the standard latitude, longitude coordinates used by most GPS devices. In this case, `` = longitude``, `` = latitude`` and `` = altitude`` + +- `UTM`_: UTM projections can be specified using a string projection string such as ``WGS84 UTM 32N``, where 32 is the region and N is . In this case, `` = E``, `` = N`` and `` = altitude`` + +- `proj4`_: Any valid proj4 format string can be used. For example, for UTM 32N we can use ``+proj=utm +zone=32 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs`` + +.. _WGS84: https://en.wikipedia.org/wiki/World_Geodetic_System +.. _UTM: https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system +.. _proj4: http://proj4.org/ + +Example +``````` +This file defines 2 GCP whose coordinates are specified in the WGS84 standard. The first one is observed in both ``01.jpg`` and ``02.jpg``, while the second one is only observed in ``01.jpg`` :: + + WGS84 + 13.400740745 52.519134104 12.0792090446 2335.0 1416.7 01.jpg + 13.400740745 52.519134104 12.0792090446 2639.1 938.0 02.jpg + 13.400502446 52.519251158 16.7021233002 766.0 1133.1 01.jpg + + diff --git a/doc/source/geometry.rst b/doc/source/geometry.rst index d0ae5f2f5..22209603e 100644 --- a/doc/source/geometry.rst +++ b/doc/source/geometry.rst @@ -10,7 +10,51 @@ TODO Coordinate Systems ------------------ -TODO +Normalized Image Coordinates +```````````````````````````` + +The 2d position of a point in images is stored in what we will call *normalized image coordinates*. The origin is in the middle of the image. The x coordinate grows to the right and y grows downwards. The larger dimension of the image is 1. + +This means, for example, that all the pixels in an image with aspect ratio 4:3 will be contained in the intervals ``[-0.5, 0.5]`` and ``[3/4 * (-0.5), 3/4 * 0.5]`` for the X and Y axis respectively. + +:: + + +-----------------------------+ + | | + | | + | | + | + -------------> + | | (0, 0) | (0.5, 0) + | | | + | | | + +-----------------------------+ + | + v + (0, 0.5) + +Normalized coordinates are independent of the resolution of the image and give better numerical stability for some multi-view geometry algorithms than pixel coordinates. + + +Pixel Coordinates +````````````````` + +Many OpenCV functions that work with images use *pixel coordinates*. In that reference frame, the origin is at the center of the top-left pixel, x grow by one for every pixel to the right and y grows by one for every pixel downwards. The bottom-right pixel is therefore at ``(width -1, height - 1)``. + + +World Coordinates +````````````````` +The position of the reconstructed 3D points is stored in *world coordinates*. In general, this is an arbitrary euclidean reference frame. + +When GPS data is available, a topocentric reference frame is used for the world coordinates reference. This is a reference frame that with the origin somewhere near the ground, the X axis pointing to the east, the Y axis pointing to the north and the Z axis pointing to the zenith. The latitude, longitude, and altitude of the origin are stored in the ``reference_lla.json`` file. + +When GPS data is not available, the reconstruction process makes its best to rotate the world reference frame so that the vertical direction is Z and the ground is near the `z = 0` plane. It does so by assuming that the images are taken from similar altitudes and that the up vector of the images corresponds to the up vector of the world. + + +Camera Coordinates +`````````````````` +The *camera coordinate* reference frame has the origin at the camera's optical center, the X axis is pointing to the right of the camera the Y axis is pointing down and the Z axis is pointing to the front. A point in front of the camera has positive Z camera coordinate. + +The pose of a camera is determined by the rotation and translation that converts world coordinates to camera coordinates. Camera Models diff --git a/doc/source/index.rst b/doc/source/index.rst index 013e2afcc..527ec71ad 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -3,10 +3,6 @@ OpenSfM ======= -Only tests here for now. - -test contents: - .. toctree:: :maxdepth: 2 diff --git a/doc/source/using.rst b/doc/source/using.rst index d40756f89..7b6d11221 100644 --- a/doc/source/using.rst +++ b/doc/source/using.rst @@ -10,7 +10,7 @@ Quickstart An example dataset is available at ``data/berlin``. You can reconstruct it using by running :: - bin/run_all data/berlin + bin/opensfm_run_all data/berlin This will run the entire SfM pipeline and produce the file ``data/berlin/reconstruction.meshed.json`` as output. To visualize the result you can start a HTTP server running :: @@ -79,38 +79,95 @@ Here is the usage page of ``bin/opensfm``, which lists the available commands extract_metadata ```````````````` -TODO + +This commands extracts EXIF metadata from the images an stores them in the ``exif`` folder and the ``camera_models.json`` file. + +The following data is extracted for each image: + +- ``width`` and ``height``: image size in pixels + +- ``gps`` ``latitude``, ``longitude``, ``altitude`` and ``dop``: The GPS coordinates of the camera at capture time and the corresponding Degree Of Precission). This is used to geolocate the reconstruction. + +- ``capture_time``: The capture time. Used to choose candidate matching images when the option ``matching_time_neighbors`` is set. + +- ``camera orientation``: The EXIF orientation tag (see this `exif orientation documentation`_). Used to orient the reconstruction straigh up. + +- ``projection_type``: The camera projection type. It is extracted from the GPano_ metadata and used to determine which projection to use for each camera. Supported types are `perspective`, `equirectangular` and `fisheye`. + +- ``focal_ratio``: The focal length provided by the EXIF metadata divided by the sensor width. This is used as initialization and prior for the camera focal length parameter. + +- ``make`` and ``model``: The camera make and model. Used to build the camera ID. + +- ``camera``: The camera ID string. Used to identify a camera. When multiple images have the same camera ID string, they will be assumed to be taken with the same camera and will share its parameters. + + +Once the metadata for all images has been extracted, a list of camera models is created and stored in ``camera_models.json``. A camera is created for each diferent camera ID string found on the images. + +For each camera the following data is stored: + +- ``width`` and ``height``: image size in pixels +- ``projection_type``: the camera projection type +- ``focal``: The initial estimation of the focal length (as a multiple of the sensor width). +- ``k1`` and ``k2``: The initial estimation of the radial distortion parameters. Only used for `perspective` and `fisheye` projection models. +- ``focal_prior``: The focal length prior. The final estimated focal length will be forced to be similar to it. +- ``k1_prior`` and ``k2_prior``: The radial distortion parameters prior. + + +Providing your own camera parameters +'''''''''''''''''''''''''''''''''''' + +By default, the camera parameters are taken from the EXIF metadata but it is also possible to override the default parameters. To do so, place a file named ``camera_models_override.json`` in the project folder. This file should have the same structure as ``camera_models.json``. When running the ``extract_metadata`` command, the parameters of any camera present in the ``camera_models_overrides.json`` file will be copied to ``camera_models.json`` overriding the default ones. + +Simplest way to create the ``camera_models_overrides.json`` file is to rename ``camera_models.json`` and modify the parameters. You will need to rerun the ``extract_metadata`` command after that. + + +.. _`exif orientation documentation`: http://sylvana.net/jpegcrop/exif_orientation.html +.. _GPano: TODO(pau): link to Google Pano metadata documentation + detect_features ``````````````` -TODO +This command detect feature points in the images and stores them in the `feature` folder. + match_features `````````````` -TODO +This command matches feature points between images and stores them in the `matches` folder. It first determines the list of image pairs to run, and then run the matching process for each pair to find corresponding feature points. + +Since there are a lot of possible image pairs, the process can be very slow. It can be speeded up by restricting the list of pairs to match. The pairs can be restricted by GPS distance, capture time or file name order. + create_tracks ````````````` -TODO +This command links the matches between pairs of images to build feature point tracks. The tracks are stored in the `tracks.csv` file. A track is a set of feature points from different images that have been recognized to correspond to the same pysical point. + reconstruct ``````````` -TODO +This command runs the incremental reconstruction process. The goal of the reconstruction process is to find the 3D position of tracks (the `structure`) together with the position of the cameras (the `motion`). The computed reconstruction is stored in the ``reconstruction.json`` file. + mesh ```` -TODO +This process computes a rough triangular mesh of the scene seen by each images. Such mesh is used for simulating smooth motions between images in the web viewer. The reconstruction with the mesh added is stored in ``reconstruction.meshed.json`` file. + +Note that the only difference between ``reconstruction.json`` and ``reconstruction.meshed.json`` is that the later contains the triangular meshes. If you don't need that, you only need the former file and there's no need to run this command. + undistort ````````` -TODO +This command creates undistorted version of the reconstruction, tracks and images. The undistorted version can later be used for computing depth maps. + compute_depthmaps ````````````````` -TODO +This commands computes a dense point cloud of the scene by computing and merging depthmaps. It requires an undistorted reconstructions. The resulting depthmaps are stored in the ``depthmaps`` folder and the merged point cloud is stored in ``depthmaps/merged.ply`` Configuration ------------- TODO explain config.yaml and the available parameters + + +.. include:: gcp.rst diff --git a/opensfm/commands/detect_features.py b/opensfm/commands/detect_features.py index cbb6a10e6..cc2879c3e 100644 --- a/opensfm/commands/detect_features.py +++ b/opensfm/commands/detect_features.py @@ -40,10 +40,14 @@ def detect(args): image, data = args logger.info('Extracting {} features for image {}'.format( data.feature_type().upper(), image)) + if not data.feature_index_exists(image): + mask = data.mask_as_array(image) + if mask is not None: + logger.info('Found mask to apply for image {}'.format(image)) preemptive_max = data.config.get('preemptive_max', 200) p_unsorted, f_unsorted, c_unsorted = features.extract_features( - data.image_as_array(image), data.config) + data.image_as_array(image), data.config, mask) if len(p_unsorted) == 0: return @@ -56,5 +60,7 @@ def detect(args): f_pre = f_sorted[-preemptive_max:] data.save_features(image, p_sorted, f_sorted, c_sorted) data.save_preemptive_features(image, p_pre, f_pre) - index = features.build_flann_index(f_sorted, data.config) - data.save_feature_index(image, index) + + if data.config.get('matcher_type', 'FLANN') == 'FLANN': + index = features.build_flann_index(f_sorted, data.config) + data.save_feature_index(image, index) diff --git a/opensfm/commands/match_features.py b/opensfm/commands/match_features.py index 76497d5a2..17ab9bfbb 100644 --- a/opensfm/commands/match_features.py +++ b/opensfm/commands/match_features.py @@ -108,6 +108,7 @@ def match_candidates_from_metadata(images, exifs, data): max_distance = data.config['matching_gps_distance'] max_neighbors = data.config['matching_gps_neighbors'] max_time_neighbors = data.config['matching_time_neighbors'] + max_order_neighbors = data.config['matching_order_neighbors'] if not all(map(has_gps_info, exifs.values())) and max_neighbors != 0: logger.warn("Not all images have GPS info. " @@ -115,24 +116,30 @@ def match_candidates_from_metadata(images, exifs, data): max_neighbors = 0 pairs = set() - for im1 in images: + images.sort() + for index1, im1 in enumerate(images): distances = [] timediffs = [] - for im2 in images: + indexdiffs = [] + for index2, im2 in enumerate(images): if im1 != im2: dx = distance_from_exif(exifs[im1], exifs[im2]) dt = timediff_from_exif(exifs[im1], exifs[im2]) + di = abs(index1 - index2) if dx <= max_distance: distances.append((dx, im2)) timediffs.append((dt, im2)) + indexdiffs.append((di, im2)) distances.sort() timediffs.sort() + indexdiffs.sort() - if max_neighbors or max_time_neighbors: + if max_neighbors or max_time_neighbors or max_order_neighbors: distances = distances[:max_neighbors] timediffs = timediffs[:max_time_neighbors] + indexdiffs = indexdiffs[:max_order_neighbors] - for d, im2 in distances + timediffs: + for d, im2 in distances + timediffs + indexdiffs: if im1 < im2: pairs.add((im1, im2)) else: diff --git a/opensfm/commands/mesh.py b/opensfm/commands/mesh.py index bbee53f5b..7cae57ae9 100644 --- a/opensfm/commands/mesh.py +++ b/opensfm/commands/mesh.py @@ -31,7 +31,8 @@ def run(self, args): shot.mesh.faces = faces data.save_reconstruction(reconstructions, - filename='reconstruction.meshed.json') + filename='reconstruction.meshed.json', + minify=True) end = time.time() with open(data.profile_log(), 'a') as fout: diff --git a/opensfm/commands/undistort.py b/opensfm/commands/undistort.py index fb350a9fc..e6666119d 100644 --- a/opensfm/commands/undistort.py +++ b/opensfm/commands/undistort.py @@ -34,16 +34,26 @@ def undistort_images(self, graph, reconstruction, data): for shot in reconstruction.shots.values(): if shot.camera.projection_type == 'perspective': + image = data.image_as_array(shot.id) + undistorted = undistort_image(image, shot.camera) + data.save_undistorted_image(shot.id, undistorted) + urec.add_camera(shot.camera) urec.add_shot(shot) - + elif shot.camera.projection_type == 'fisheye': image = data.image_as_array(shot.id) - undistorted = undistort_image(image, shot) + undistorted = undistort_fisheye_image(image, shot.camera) data.save_undistorted_image(shot.id, undistorted) + + shot.camera = perspective_camera_from_fisheye(shot.camera) + urec.add_camera(shot.camera) + urec.add_shot(shot) elif shot.camera.projection_type in ['equirectangular', 'spherical']: original = data.image_as_array(shot.id) - image = cv2.resize(original, (2048, 1024), interpolation=cv2.INTER_AREA) - shots = perspective_views_of_a_panorama(shot) + width = 4 * int(data.config['depthmap_resolution']) + height = width / 2 + image = cv2.resize(original, (width, height), interpolation=cv2.INTER_AREA) + shots = perspective_views_of_a_panorama(shot, width) for subshot in shots: urec.add_camera(subshot.camera) urec.add_shot(subshot) @@ -56,21 +66,40 @@ def undistort_images(self, graph, reconstruction, data): data.save_undistorted_reconstruction([urec]) -def undistort_image(image, shot): +def undistort_image(image, camera): """Remove radial distortion from a perspective image.""" - camera = shot.camera height, width = image.shape[:2] K = camera.get_K_in_pixel_coordinates(width, height) distortion = np.array([camera.k1, camera.k2, 0, 0]) return cv2.undistort(image, K, distortion) -def perspective_views_of_a_panorama(spherical_shot): +def undistort_fisheye_image(image, camera): + """Remove radial distortion from a perspective image.""" + height, width = image.shape[:2] + K = camera.get_K_in_pixel_coordinates(width, height) + distortion = np.array([camera.k1, camera.k2, 0, 0]) + return cv2.fisheye.undistortImage(image, K, distortion, K) + + +def perspective_camera_from_fisheye(fisheye): + """Create a perspective camera from a fisheye.""" + camera = types.PerspectiveCamera() + camera.id = fisheye.id + camera.width = fisheye.width + camera.height = fisheye.height + camera.focal = fisheye.focal + camera.focal_prior = fisheye.focal_prior + camera.k1 = camera.k1_prior = camera.k2 = camera.k2_prior = 0.0 + return camera + + +def perspective_views_of_a_panorama(spherical_shot, width): """Create 6 perspective views of a panorama.""" camera = types.PerspectiveCamera() camera.id = 'perspective_panorama_camera' - camera.width = 640 - camera.height = 640 + camera.width = width + camera.height = width camera.focal = 0.5 camera.focal_prior = camera.focal camera.k1 = camera.k1_prior = camera.k2 = camera.k2_prior = 0.0 @@ -125,19 +154,19 @@ def render_perspective_view_of_a_panorama(image, panoshot, perspectiveshot): src_pixels = np.column_stack([src_x.ravel(), src_y.ravel()]) src_pixels_denormalized = features.denormalized_image_coordinates( - src_pixels, - image.shape[1], - image.shape[0]) + src_pixels, image.shape[1], image.shape[0]) # Sample color - colors = image[src_pixels_denormalized[:, 1].astype(int), - src_pixels_denormalized[:, 0].astype(int)] + colors = cv2.remap(image, + src_pixels_denormalized[:, 0].astype(np.float32), + src_pixels_denormalized[:, 1].astype(np.float32), + cv2.INTER_LINEAR) colors.shape = dst_shape + (-1,) return colors def add_subshot_tracks(graph, panoshot, perspectiveshot): - """Add edges betwene subshots and visible tracks.""" + """Add edges between subshots and visible tracks.""" graph.add_node(perspectiveshot.id, bipartite=0) for track in graph[panoshot.id]: edge = graph[panoshot.id][track] diff --git a/opensfm/config.py b/opensfm/config.py index 2101451a9..a3b1c45d0 100644 --- a/opensfm/config.py +++ b/opensfm/config.py @@ -33,10 +33,6 @@ hahog_peak_threshold: 0.00001 hahog_edge_threshold: 10 -# Masks for regions that will be ignored for feature extraction -# List of bounding boxes specified as the ratio to image width and height -# masks: [{top: 0.96, bottom: 1.0, left: 0.0, right: 0.15}, {top: 0.95, bottom: 1.0, left: 0, right: 0.05}] - # Params for general matching lowes_ratio: 0.8 # Ratio test for matches preemptive_lowes_ratio: 0.6 # Ratio test for preemptive matches @@ -51,6 +47,7 @@ matching_gps_distance: 150 # Maximum gps distance between two images for matching matching_gps_neighbors: 0 # Number of images to match selected by GPS distance. Set to 0 to use no limit matching_time_neighbors: 0 # Number of images to match selected by time taken. Set to 0 to use no limit +matching_order_neighbors: 0 # Number of images to match selected by image name. Set to 0 to use no limit preemptive_max: 200 # Number of features to use for preemptive matching preemptive_threshold: 0 # If number of matches passes the threshold -> full feature matching @@ -102,6 +99,7 @@ nav_rotation_threshold: 30 # Maximum general rotation in degrees between cameras for steps # Params for depth estimation +depthmap_method: PATCH_MATCH # Raw depthmap computationg algorithm (PATCH_MATCH, BRUTE_FORCE, PATCH_MATCH_SAMPLE) depthmap_resolution: 640 # Resolution of the depth maps depthmap_num_neighbors: 10 # Number of neighboring views depthmap_num_matching_views: 2 # Number of neighboring views used for each depthmaps diff --git a/opensfm/data/sensor_data.json b/opensfm/data/sensor_data.json index 68a5ab14f..7d1ee637f 100644 --- a/opensfm/data/sensor_data.json +++ b/opensfm/data/sensor_data.json @@ -701,8 +701,10 @@ "Contax U4R": 5.33, "Contax i4R": 5.33, "DJI FC300S": 6.16, + "DJI FC300C": 6.31, "DJI FC300X": 6.2, "DJI FC350": 6.17, + "DJI FC330": 6.25, "Epson L-500V": 5.75, "Epson PhotoPC 3000 Zoom": 7.11, "Epson PhotoPC 3100 Zoom": 7.11, diff --git a/opensfm/dataset.py b/opensfm/dataset.py index 61077252d..bbf40125a 100644 --- a/opensfm/dataset.py +++ b/opensfm/dataset.py @@ -41,6 +41,15 @@ def __init__(self, data_path): else: self.set_image_path(os.path.join(self.data_path, 'images')) + # Load list of masks if they exist. + mask_list_file = os.path.join(self.data_path, 'mask_list.txt') + if os.path.isfile(mask_list_file): + with open(mask_list_file) as fin: + lines = fin.read().splitlines() + self.set_mask_list(lines) + else: + self.set_mask_path(os.path.join(self.data_path, 'masks')) + def _load_config(self): config_file = os.path.join(self.data_path, 'config.yaml') self.config = config.load_config(config_file) @@ -57,12 +66,11 @@ def __image_file(self, image): return self.image_files[image] def load_image(self, image): - return open(self.__image_file(image)) + return open(self.__image_file(image), 'rb') def image_as_array(self, image): """Return image pixels as 3-dimensional numpy array (R G B order)""" - IMREAD_COLOR = cv2.IMREAD_COLOR if context.OPENCV3 else cv2.CV_LOAD_IMAGE_COLOR - return cv2.imread(self.__image_file(image), IMREAD_COLOR)[:,:,::-1] # Turn BGR to RGB + return io.imread(self.__image_file(image)) def _undistorted_image_path(self): return os.path.join(self.data_path, 'undistorted') @@ -73,13 +81,28 @@ def _undistorted_image_file(self, image): def undistorted_image_as_array(self, image): """Undistorted image pixels as 3-dimensional numpy array (R G B order)""" - IMREAD_COLOR = cv2.IMREAD_COLOR if context.OPENCV3 else cv2.CV_LOAD_IMAGE_COLOR - return cv2.imread(self._undistorted_image_file(image), IMREAD_COLOR)[:,:,::-1] # Turn BGR to RGB + return io.imread(self._undistorted_image_file(image)) def save_undistorted_image(self, image, array): io.mkdir_p(self._undistorted_image_path()) cv2.imwrite(self._undistorted_image_file(image), array[:, :, ::-1]) + def masks(self): + """Return list of file names of all masks in this dataset""" + return self.mask_list + + def mask_as_array(self, image): + """Given an image, returns the associated mask as an array if it exists, otherwise returns None""" + mask_name = image + '.png' + if mask_name in self.masks(): + mask_path = self.mask_files[mask_name] + mask = cv2.imread(mask_path) + if len(mask.shape) == 3: + mask = mask.max(axis=2) + else: + mask = None + return mask + def _depthmap_path(self): return os.path.join(self.data_path, 'depthmaps') @@ -90,14 +113,14 @@ def _depthmap_file(self, image, suffix): def raw_depthmap_exists(self, image): return os.path.isfile(self._depthmap_file(image, 'raw.npz')) - def save_raw_depthmap(self, image, depth, plane, score): + def save_raw_depthmap(self, image, depth, plane, score, nghbr, nghbrs): io.mkdir_p(self._depthmap_path()) filepath = self._depthmap_file(image, 'raw.npz') - np.savez_compressed(filepath, depth=depth, plane=plane, score=score) + np.savez_compressed(filepath, depth=depth, plane=plane, score=score, nghbr=nghbr, nghbrs=nghbrs) def load_raw_depthmap(self, image): o = np.load(self._depthmap_file(image, 'raw.npz')) - return o['depth'], o['plane'], o['score'] + return o['depth'], o['plane'], o['score'], o['nghbr'], o['nghbrs'] def clean_depthmap_exists(self, image): return os.path.isfile(self._depthmap_file(image, 'clean.npz')) @@ -116,7 +139,7 @@ def __is_image_file(filename): return filename.split('.')[-1].lower() in {'jpg', 'jpeg', 'png', 'tif', 'tiff', 'pgm', 'pnm', 'gif'} def set_image_path(self, path): - """Set image path and find the all images in there""" + """Set image path and find all images in there""" self.image_list = [] self.image_files = {} if os.path.exists(path): @@ -134,6 +157,29 @@ def set_image_list(self, image_list): self.image_list.append(name) self.image_files[name] = path + @staticmethod + def __is_mask_file(filename): + return DataSet.__is_image_file(filename) + + def set_mask_path(self, path): + """Set mask path and find all masks in there""" + self.mask_list = [] + self.mask_files = {} + if os.path.exists(path): + for name in os.listdir(path): + if self.__is_mask_file(name): + self.mask_list.append(name) + self.mask_files[name] = os.path.join(path, name) + + def set_mask_list(self, mask_list): + self.mask_list = [] + self.mask_files = {} + for line in mask_list: + path = os.path.join(self.data_path, line) + name = os.path.basename(path) + self.mask_list.append(name) + self.mask_files[name] = path + def __exif_path(self): """Return path of extracted exif directory""" return os.path.join(self.data_path, 'exif') @@ -159,12 +205,12 @@ def load_exif(self, image): :param image: Image name, with extension (i.e. 123.jpg) """ - with open(self.__exif_file(image), 'r') as fin: + with open(self.__exif_file(image), 'rb') as fin: return json.load(fin) def save_exif(self, image, data): io.mkdir_p(self.__exif_path()) - with open(self.__exif_file(image), 'w') as fout: + with open(self.__exif_file(image), 'wb') as fout: io.json_dump(data, fout) def feature_type(self): @@ -174,27 +220,16 @@ def feature_type(self): if self.config.get('feature_root', False): feature_name = 'root_' + feature_name return feature_name - def descriptor_type(self): - """Return the type of the descriptor (if exists) - """ - if self.feature_type() == 'akaze': - return self.config.get('akaze_descriptor', '') - else: - return '' - def __feature_path(self): """Return path of feature descriptors and FLANN indices directory""" - __feature_path = self.feature_type() - if len(self.descriptor_type()) > 0: - __feature_path += '_' + self.descriptor_type() - return os.path.join(self.data_path, __feature_path) + return os.path.join(self.data_path, "features") def __feature_file(self, image): """ Return path of feature file for specified image :param image: Image name, with extension (i.e. 123.jpg) """ - return os.path.join(self.__feature_path(), image + '.' + self.feature_type() + '.npz') + return os.path.join(self.__feature_path(), image + '.npz') def __save_features(self, filepath, image, points, descriptors, colors=None): io.mkdir_p(self.__feature_path()) @@ -232,7 +267,7 @@ def __feature_index_file(self, image): Return path of FLANN index file for specified image :param image: Image name, with extension (i.e. 123.jpg) """ - return os.path.join(self.__feature_path(), image + '.' + self.feature_type() + '.flann') + return os.path.join(self.__feature_path(), image + '.flann') def load_feature_index(self, image, features): index = cv2.flann.Index() if context.OPENCV3 else cv2.flann_Index() @@ -248,7 +283,7 @@ def __preemptive_features_file(self, image): for specified image :param image: Image name, with extension (i.e. 123.jpg) """ - return os.path.join(self.__feature_path(), image + '_preemptive.' + self.feature_type() + '.npz') + return os.path.join(self.__feature_path(), image + '_preemptive' + '.npz') def load_preemtive_features(self, image): s = np.load(self.__preemptive_features_file(image)) @@ -257,16 +292,6 @@ def load_preemtive_features(self, image): def save_preemptive_features(self, image, points, descriptors): self.__save_features(self.__preemptive_features_file(image), image, points, descriptors) - def matcher_type(self): - """Return the type of matcher - """ - matcher_type = self.config.get('matcher_type', 'BruteForce') - if 'BruteForce' in matcher_type: - if self.feature_type() == 'akaze' and (self.config.get('akaze_descriptor', 5) >= 4): - matcher_type = 'BruteForce-Hamming' - self.config['matcher_type'] = matcher_type - return matcher_type # BruteForce, BruteForce-L1, BruteForce-Hamming - def __matches_path(self): """Return path of matches directory""" return os.path.join(self.data_path, 'matches') @@ -328,9 +353,9 @@ def load_reconstruction(self, filename=None): reconstructions = io.reconstructions_from_json(json.load(fin)) return reconstructions - def save_reconstruction(self, reconstruction, filename=None): + def save_reconstruction(self, reconstruction, filename=None, minify=False): with open(self.__reconstruction_file(filename), 'w') as fout: - io.json_dump(io.reconstructions_to_json(reconstruction), fout) + io.json_dump(io.reconstructions_to_json(reconstruction), fout, minify) def load_undistorted_reconstruction(self): return self.load_reconstruction( diff --git a/opensfm/dense.py b/opensfm/dense.py index 1668328e5..6dac650ad 100644 --- a/opensfm/dense.py +++ b/opensfm/dense.py @@ -56,23 +56,34 @@ def parallel_run(function, arguments, num_processes): def compute_depthmap(arguments): """Compute depthmap for a single shot.""" data, reconstruction, neighbors, min_depth, max_depth, shot = arguments + method = data.config['depthmap_method'] if data.raw_depthmap_exists(shot.id): logger.info("Using precomputed raw depthmap {}".format(shot.id)) return - logger.info("Computing depthmap for image {}".format(shot.id)) + logger.info("Computing depthmap for image {0} with {1}".format(shot.id, method)) de = csfm.DepthmapEstimator() de.set_depth_range(min_depth, max_depth, 100) de.set_patchmatch_iterations(data.config['depthmap_patchmatch_iterations']) de.set_min_patch_sd(data.config['depthmap_min_patch_sd']) add_views_to_depth_estimator(data, reconstruction, neighbors[shot.id], de) - depth, plane, score = de.compute_patch_match() + + if (method == 'BRUTE_FORCE'): + depth, plane, score, nghbr = de.compute_brute_force() + elif (method == 'PATCH_MATCH'): + depth, plane, score, nghbr = de.compute_patch_match() + elif (method == 'PATCH_MATCH_SAMPLE'): + depth, plane, score, nghbr = de.compute_patch_match_sample() + else: + raise ValueError('Unknown depthmap method type ' \ + '(must be BRUTE_FORCE, PATCH_MATCH or PATCH_MATCH_SAMPLE)') + good_score = score > data.config['depthmap_min_correlation_score'] depth = depth * (depth < max_depth) * good_score # Save and display results - data.save_raw_depthmap(shot.id, depth, plane, score) + data.save_raw_depthmap(shot.id, depth, plane, score, nghbr, neighbors[shot.id][1:]) if data.config['depthmap_save_debug_files']: image = data.undistorted_image_as_array(shot.id) @@ -83,16 +94,21 @@ def compute_depthmap(arguments): if data.config.get('interactive'): import matplotlib.pyplot as plt - plt.subplot(2, 2, 1) + plt.figure() + plt.suptitle("Shot: " + shot.id + ", neighbors: " + ', '.join(neighbors[shot.id][1:])) + plt.subplot(2, 3, 1) plt.imshow(image) - plt.subplot(2, 2, 2) + plt.subplot(2, 3, 2) plt.imshow(color_plane_normals(plane)) - plt.subplot(2, 2, 3) + plt.subplot(2, 3, 3) plt.imshow(depth) plt.colorbar() - plt.subplot(2, 2, 4) + plt.subplot(2, 3, 4) plt.imshow(score) plt.colorbar() + plt.subplot(2, 3, 5) + plt.imshow(nghbr) + plt.colorbar() plt.show() @@ -112,7 +128,7 @@ def clean_depthmap(arguments): depth = dc.clean() # Save and display results - raw_depth, raw_plane, raw_score = data.load_raw_depthmap(shot.id) + raw_depth, raw_plane, raw_score, raw_nghbr, nghbrs = data.load_raw_depthmap(shot.id) data.save_clean_depthmap(shot.id, depth, raw_plane, raw_score) if data.config['depthmap_save_debug_files']: @@ -124,6 +140,8 @@ def clean_depthmap(arguments): if data.config.get('interactive'): import matplotlib.pyplot as plt + plt.figure() + plt.suptitle("Shot: " + shot.id) plt.subplot(2, 2, 1) plt.imshow(raw_depth) plt.colorbar() @@ -190,7 +208,7 @@ def add_views_to_depth_cleaner(data, reconstruction, neighbors, dc): shot = reconstruction.shots[neighbor] if not data.raw_depthmap_exists(shot.id): continue - depth, plane, score = data.load_raw_depthmap(shot.id) + depth, plane, score, nghbr, nghbrs = data.load_raw_depthmap(shot.id) height, width = depth.shape K = shot.camera.get_K_in_pixel_coordinates(width, height) R = shot.pose.get_rotation_matrix() @@ -212,7 +230,7 @@ def compute_depth_range(graph, reconstruction, shot): def find_neighboring_images(shot, common_tracks, reconstruction, num_neighbors=5): - """Find neighbouring images based on common tracks.""" + """Find neighboring images based on common tracks.""" theta_min = np.pi / 60 theta_max = np.pi / 6 ns = [] diff --git a/opensfm/exif.py b/opensfm/exif.py index f20fa41c8..e20d29bc7 100644 --- a/opensfm/exif.py +++ b/opensfm/exif.py @@ -31,7 +31,10 @@ def get_float_tag(tags, key): def get_frac_tag(tags, key): if key in tags: - return eval_frac(tags[key].values[0]) + try: + return eval_frac(tags[key].values[0]) + except ZeroDivisionError: + return None else: return None diff --git a/opensfm/features.py b/opensfm/features.py index 54cf1dde5..949e04000 100644 --- a/opensfm/features.py +++ b/opensfm/features.py @@ -1,13 +1,8 @@ -# -*- coding: utf-8 -*- +"""Tools to extract features.""" -import os, sys -import tempfile import time import logging -from subprocess import call import numpy as np -import json -import uuid import cv2 import csfm @@ -18,14 +13,17 @@ def resized_image(image, config): - feature_process_size = config.get('feature_process_size', -1) - size = np.array(image.shape[0:2]) - if 0 < feature_process_size < size.max(): - new_size = size * feature_process_size / size.max() - return cv2.resize(image, dsize=(new_size[1], new_size[0])) + """Resize image to feature_process_size.""" + max_size = config.get('feature_process_size', -1) + h, w, _ = image.shape + size = max(w, h) + if 0 < max_size < size: + dsize = w * max_size / size, h * max_size / size + return cv2.resize(image, dsize=dsize, interpolation=cv2.INTER_AREA) else: return image + def root_feature(desc, l2_normalization=False): if l2_normalization: s2 = np.linalg.norm(desc, axis=1) @@ -34,6 +32,7 @@ def root_feature(desc, l2_normalization=False): desc = np.sqrt(desc.T/s).T return desc + def root_feature_surf(desc, l2_normalization=False, partial=False): """ Experimental square root mapping of surf-like feature, only work for 64-dim surf now @@ -54,6 +53,7 @@ def root_feature_surf(desc, l2_normalization=False, partial=False): desc[:, ii] = desc_sub*desc_sub_sign return desc + def normalized_image_coordinates(pixel_coords, width, height): size = max(width, height) p = np.empty((len(pixel_coords), 2)) @@ -61,6 +61,7 @@ def normalized_image_coordinates(pixel_coords, width, height): p[:, 1] = (pixel_coords[:, 1] + 0.5 - height / 2.0) / size return p + def denormalized_image_coordinates(norm_coords, width, height): size = max(width, height) p = np.empty((len(norm_coords), 2)) @@ -68,23 +69,27 @@ def denormalized_image_coordinates(norm_coords, width, height): p[:, 1] = norm_coords[:, 1] * size - 0.5 + height / 2.0 return p -def mask_and_normalize_features(points, desc, colors, width, height, config): - masks = np.array(config.get('masks',[])) - for mask in masks: - top = mask['top'] * height - left = mask['left'] * width - bottom = mask['bottom'] * height - right = mask['right'] * width - ids = np.invert ( (points[:,1] > top) * - (points[:,1] < bottom) * - (points[:,0] > left) * - (points[:,0] < right) ) + +def mask_and_normalize_features(points, desc, colors, width, height, mask=None): + """Remove features outside the mask and normalize image coordinates.""" + + if mask is not None: + ids = np.array([_in_mask(point, width, height, mask) for point in points]) points = points[ids] desc = desc[ids] colors = colors[ids] + points[:, :2] = normalized_image_coordinates(points[:, :2], width, height) return points, desc, colors + +def _in_mask(point, width, height, mask): + """Check if a point is inside a binary mask.""" + u = mask.shape[1] * (point[0] + 0.5) / width + v = mask.shape[0] * (point[1] + 0.5) / height + return mask[int(v), int(u)] != 0 + + def extract_features_sift(image, config): sift_edge_threshold = config.get('sift_edge_threshold', 10) sift_peak_threshold = float(config.get('sift_peak_threshold', 0.1)) @@ -124,6 +129,7 @@ def extract_features_sift(image, config): points = np.array([(i.pt[0], i.pt[1], i.size, i.angle) for i in points]) return points, desc + def extract_features_surf(image, config): surf_hessian_threshold = config.get('surf_hessian_threshold', 3000) if context.OPENCV3: @@ -167,6 +173,7 @@ def extract_features_surf(image, config): points = np.array([(i.pt[0], i.pt[1], i.size, i.angle) for i in points]) return points, desc + def akaze_descriptor_type(name): d = csfm.AkazeDescriptorType.__dict__ if name in d: @@ -175,6 +182,7 @@ def akaze_descriptor_type(name): logger.debug('Wrong akaze descriptor type') return d['MSURF'] + def extract_features_akaze(image, config): options = csfm.AKAZEOptions() options.omax = config.get('akaze_omax', 4) @@ -202,6 +210,7 @@ def extract_features_akaze(image, config): points = points.astype(float) return points, desc + def extract_features_hahog(image, config): t = time.time() points, desc = csfm.hahog(image.astype(np.float32) / 255, # VlFeat expects pixel values between 0, 1 @@ -222,7 +231,8 @@ def extract_features_hahog(image, config): logger.debug('Found {0} points in {1}s'.format( len(points), time.time()-t )) return points, desc -def extract_features(color_image, config): + +def extract_features(color_image, config, mask=None): assert len(color_image.shape) == 3 color_image = resized_image(color_image, config) image = cv2.cvtColor(color_image, cv2.COLOR_RGB2GRAY) @@ -243,7 +253,7 @@ def extract_features(color_image, config): ys = points[:,1].round().astype(int) colors = color_image[ys, xs] - return mask_and_normalize_features(points, desc, colors, image.shape[1], image.shape[0], config) + return mask_and_normalize_features(points, desc, colors, image.shape[1], image.shape[0], mask) def build_flann_index(features, config): diff --git a/opensfm/io.py b/opensfm/io.py index 821dbfc24..b2eb17da8 100644 --- a/opensfm/io.py +++ b/opensfm/io.py @@ -1,5 +1,6 @@ import errno import json +import logging import os import cv2 @@ -9,6 +10,10 @@ from opensfm import features from opensfm import geo from opensfm import types +from opensfm import context + + +logger = logging.getLogger(__name__) def camera_from_json(key, obj): @@ -28,6 +33,18 @@ def camera_from_json(key, obj): camera.k1_prior = obj.get('k1_prior', camera.k1) camera.k2_prior = obj.get('k2_prior', camera.k2) return camera + elif pt == 'fisheye': + camera = types.FisheyeCamera() + camera.id = key + camera.width = obj.get('width', 0) + camera.height = obj.get('height', 0) + camera.focal = obj['focal'] + camera.k1 = obj.get('k1', 0.0) + camera.k2 = obj.get('k2', 0.0) + camera.focal_prior = obj.get('focal_prior', camera.focal) + camera.k1_prior = obj.get('k1_prior', camera.k1) + camera.k2_prior = obj.get('k2_prior', camera.k2) + return camera elif pt in ['equirectangular', 'spherical']: camera = types.SphericalCamera() camera.id = key @@ -153,6 +170,18 @@ def camera_to_json(camera): 'k1_prior': camera.k1_prior, 'k2_prior': camera.k2_prior } + elif camera.projection_type == 'fisheye': + return { + 'projection_type': camera.projection_type, + 'width': camera.width, + 'height': camera.height, + 'focal': camera.focal, + 'k1': camera.k1, + 'k2': camera.k2, + 'focal_prior': camera.focal_prior, + 'k1_prior': camera.k1_prior, + 'k2_prior': camera.k2_prior + } elif camera.projection_type in ['equirectangular', 'spherical']: return { 'projection_type': camera.projection_type, @@ -345,14 +374,35 @@ def mkdir_p(path): raise -def json_dump(data, fout, indent=4, codec='utf-8'): - return json.dump(data, fout, indent=indent, ensure_ascii=False, encoding=codec) +def json_dump(data, fout, minify=False, codec='utf-8'): + if minify: + indent, separators = None, (',',':') + else: + indent, separators = 4, None + return json.dump(data, fout, indent=indent, ensure_ascii=False, encoding=codec, separators=separators) def json_loads(text, codec='utf-8'): return json.loads(text.decode(codec)) +def imread(filename): + """Load image as an RGB array ignoring EXIF orientation.""" + if context.OPENCV3: + flags = cv2.IMREAD_COLOR + try: + flags |= cv2.IMREAD_IGNORE_ORIENTATION + except AttributeError: + logger.warning( + "OpenCV version {} does not support loading images without " + "rotating them according to EXIF. Please upgrade OpenCV to " + "version 3.2 or newer.".format(cv2.__version__)) + else: + flags = cv2.CV_LOAD_IMAGE_COLOR + bgr = cv2.imread(filename, flags) + return bgr[:, :, ::-1] # Turn BGR to RGB + + # Bundler def export_bundler(image_list, reconstructions, track_graph, bundle_file_path, @@ -473,7 +523,7 @@ def import_bundler(data_path, bundle_file, list_file, track_file, focal, k1, k2 = map(float, lines[offset].rstrip('\n').split(' ')) if focal > 0: - im = cv2.imread(os.path.join(data_path, image_list[i])) + im = imread(os.path.join(data_path, image_list[i])) height, width = im.shape[0:2] camera = types.PerspectiveCamera() camera.id = 'camera_' + str(i) diff --git a/opensfm/mesh.py b/opensfm/mesh.py index 260bcfd0d..3afc5db0d 100644 --- a/opensfm/mesh.py +++ b/opensfm/mesh.py @@ -18,6 +18,8 @@ def triangle_mesh(shot_id, r, graph, data): if shot.camera.projection_type == 'perspective': return triangle_mesh_perspective(shot_id, r, graph) + elif shot.camera.projection_type == 'fisheye': + return triangle_mesh_fisheye(shot_id, r, graph) elif shot.camera.projection_type in ['equirectangular', 'spherical']: return triangle_mesh_equirectangular(shot_id, r, graph) else: @@ -76,6 +78,53 @@ def back_project_no_distortion(shot, pixel, depth): return shot.pose.transform_inverse(p) +def triangle_mesh_fisheye(shot_id, r, graph): + shot = r.shots[shot_id] + + bearings = [] + vertices = [] + + # Add boundary vertices + num_circle_points = 20 + for i in range(num_circle_points): + a = 2 * np.pi * float(i) / num_circle_points + point = 30 * np.array([np.cos(a), np.sin(a), 0]) + bearing = point / np.linalg.norm(point) + point = shot.pose.transform_inverse(point) + vertices.append(point.tolist()) + bearings.append(bearing) + + # Add a single vertex in front of the camera + point = 30 * np.array([0, 0, 1]) + bearing = 0.3 * point / np.linalg.norm(point) + point = shot.pose.transform_inverse(point) + vertices.append(point.tolist()) + bearings.append(bearing) + + # Add reconstructed points + for track_id, edge in graph[shot_id].items(): + if track_id in r.points: + point = r.points[track_id].coordinates + vertices.append(point) + direction = shot.pose.transform(point) + pixel = direction / np.linalg.norm(direction) + bearings.append(pixel.tolist()) + + # Triangulate + tri = scipy.spatial.ConvexHull(bearings) + faces = tri.simplices.tolist() + + # Remove faces having only boundary vertices + def good_face(face): + return (face[0] >= num_circle_points or + face[1] >= num_circle_points or + face[2] >= num_circle_points) + + faces = filter(good_face, faces) + + return vertices, faces + + def triangle_mesh_equirectangular(shot_id, r, graph): shot = r.shots[shot_id] @@ -100,4 +149,5 @@ def triangle_mesh_equirectangular(shot_id, r, graph): tri = scipy.spatial.ConvexHull(bearings) faces = tri.simplices.tolist() + return vertices, faces diff --git a/opensfm/reconstruction.py b/opensfm/reconstruction.py index addacb4ee..de4a8cce6 100644 --- a/opensfm/reconstruction.py +++ b/opensfm/reconstruction.py @@ -32,7 +32,11 @@ def bundle(graph, reconstruction, gcp, config, fix_cameras=False): str(camera.id), camera.focal, camera.k1, camera.k2, camera.focal_prior, camera.k1_prior, camera.k2_prior, fix_cameras) - + elif camera.projection_type == 'fisheye': + ba.add_fisheye_camera( + str(camera.id), camera.focal, camera.k1, camera.k2, + camera.focal_prior, camera.k1_prior, camera.k2_prior, + fix_cameras) elif camera.projection_type in ['equirectangular', 'spherical']: ba.add_equirectangular_camera(str(camera.id)) @@ -60,14 +64,14 @@ def bundle(graph, reconstruction, gcp, config, fix_cameras=False): if config['bundle_use_gps']: for shot in reconstruction.shots.values(): g = shot.metadata.gps_position - ba.add_position_prior(shot.id, g[0], g[1], g[2], + ba.add_position_prior(str(shot.id), g[0], g[1], g[2], shot.metadata.gps_dop) if config['bundle_use_gcp'] and gcp: for observation in gcp: if observation.shot_id in reconstruction.shots: ba.add_ground_control_point_observation( - observation.shot_id, + str(observation.shot_id), observation.coordinates[0], observation.coordinates[1], observation.coordinates[2], @@ -96,6 +100,11 @@ def bundle(graph, reconstruction, gcp, config, fix_cameras=False): camera.focal = c.focal camera.k1 = c.k1 camera.k2 = c.k2 + elif camera.projection_type == 'fisheye': + c = ba.get_fisheye_camera(str(camera.id)) + camera.focal = c.focal + camera.k1 = c.k1 + camera.k2 = c.k2 for shot in reconstruction.shots.values(): s = ba.get_shot(str(shot.id)) @@ -123,6 +132,10 @@ def bundle_single_view(graph, reconstruction, shot_id, config): ba.add_perspective_camera( str(camera.id), camera.focal, camera.k1, camera.k2, camera.focal_prior, camera.k1_prior, camera.k2_prior, True) + elif camera.projection_type == 'fisheye': + ba.add_fisheye_camera( + str(camera.id), camera.focal, camera.k1, camera.k2, + camera.focal_prior, camera.k1_prior, camera.k2_prior, True) elif camera.projection_type in ['equirectangular', 'spherical']: ba.add_equirectangular_camera(str(camera.id)) @@ -145,7 +158,7 @@ def bundle_single_view(graph, reconstruction, shot_id, config): if config['bundle_use_gps']: g = shot.metadata.gps_position - ba.add_position_prior(shot.id, g[0], g[1], g[2], + ba.add_position_prior(str(shot.id), g[0], g[1], g[2], shot.metadata.gps_dop) ba.set_loss_function(config.get('loss_function', 'SoftLOneLoss'), @@ -216,7 +229,7 @@ def bundle_local(graph, reconstruction, gcp, central_shot_id, config): for shot_id in interior: shot = reconstruction.shots[shot_id] g = shot.metadata.gps_position - ba.add_position_prior(shot.id, g[0], g[1], g[2], + ba.add_position_prior(str(shot.id), g[0], g[1], g[2], shot.metadata.gps_dop) if config['bundle_use_gcp'] and gcp: @@ -831,7 +844,7 @@ def grow_reconstruction(data, graph, reconstruction, images, gcp): if data.config.get('save_partial_reconstructions', False): paint_reconstruction(data, graph, reconstruction) data.save_reconstruction( - reconstruction, 'reconstruction.{}.json'.format( + [reconstruction], 'reconstruction.{}.json'.format( datetime.datetime.now().isoformat().replace(':', '_'))) common_tracks = reconstructed_points_for_images(graph, reconstruction, diff --git a/opensfm/src/CMakeLists.txt b/opensfm/src/CMakeLists.txt index e0f9a1625..a275b0d4e 100644 --- a/opensfm/src/CMakeLists.txt +++ b/opensfm/src/CMakeLists.txt @@ -77,6 +77,10 @@ target_link_libraries(akaze ${OpenCV_LIBS}) # VLFeat include_directories(third_party/vlfeat) file(GLOB VLFEAT_SRCS third_party/vlfeat/vl/*.c third_party/vlfeat/vl/*.h) +if (NOT CMAKE_SYSTEM_PROCESSOR MATCHES + "(x86)|(X86)|(x86_64)|(X86_64)|(amd64)|(AMD64)") + add_definitions(-DVL_DISABLE_SSE2) +endif () add_definitions(-DVL_DISABLE_AVX) add_library(vl ${VLFEAT_SRCS}) diff --git a/opensfm/src/akaze.cc b/opensfm/src/akaze.cc index 749a2100d..22af173b8 100644 --- a/opensfm/src/akaze.cc +++ b/opensfm/src/akaze.cc @@ -42,8 +42,14 @@ bp::object akaze(PyObject *image, bp::list retn; npy_intp keys_shape[2] = {keys.rows, keys.cols}; retn.append(bpn_array_from_data(2, keys_shape, keys.ptr(0))); - npy_intp desc_shape[2] = {desc.rows, desc.cols}; - retn.append(bpn_array_from_data(2, desc_shape, desc.ptr(0))); + + if (options.descriptor == MLDB_UPRIGHT || options.descriptor == MLDB) { + npy_intp desc_shape[2] = {desc.rows, desc.cols}; + retn.append(bpn_array_from_data(2, desc_shape, desc.ptr(0))); + } else { + npy_intp desc_shape[2] = {desc.rows, desc.cols}; + retn.append(bpn_array_from_data(2, desc_shape, desc.ptr(0))); + } return retn; } diff --git a/opensfm/src/bundle.h b/opensfm/src/bundle.h index 83d83b5b1..fc78ebd4e 100644 --- a/opensfm/src/bundle.h +++ b/opensfm/src/bundle.h @@ -18,6 +18,7 @@ extern "C" { enum BACameraType { BA_PERSPECTIVE_CAMERA, + BA_FISHEYE_CAMERA, BA_EQUIRECTANGULAR_CAMERA }; @@ -50,6 +51,21 @@ struct BAPerspectiveCamera : public BACamera{ void SetK2(double v) { parameters[BA_CAMERA_K2] = v; } }; +struct BAFisheyeCamera : public BACamera{ + double parameters[BA_CAMERA_NUM_PARAMS]; + double focal_prior; + double k1_prior; + double k2_prior; + + BACameraType type() { return BA_FISHEYE_CAMERA; } + double GetFocal() { return parameters[BA_CAMERA_FOCAL]; } + double GetK1() { return parameters[BA_CAMERA_K1]; } + double GetK2() { return parameters[BA_CAMERA_K2]; } + void SetFocal(double v) { parameters[BA_CAMERA_FOCAL] = v; } + void SetK1(double v) { parameters[BA_CAMERA_K1] = v; } + void SetK2(double v) { parameters[BA_CAMERA_K2] = v; } +}; + struct BAEquirectangularCamera : public BACamera { BACameraType type() { return BA_EQUIRECTANGULAR_CAMERA; } }; @@ -231,6 +247,63 @@ struct PerspectiveReprojectionError { double scale_; }; +template +void FisheyeProject(const T* const camera, + const T point[3], + T projection[2]) { + const T& focal = camera[BA_CAMERA_FOCAL]; + const T& k1 = camera[BA_CAMERA_K1]; + const T& k2 = camera[BA_CAMERA_K2]; + const T &x = point[0]; + const T &y = point[1]; + const T &z = point[2]; + + T l = sqrt(x * x + y * y); + T theta = atan2(l, z); + T theta2 = theta * theta; + T theta_d = theta * (T(1.0) + theta2 * (k1 + theta2 * k2)); + T s = focal * theta_d / l; + + projection[0] = s * x; + projection[1] = s * y; +} + + +struct FisheyeReprojectionError { + FisheyeReprojectionError(double observed_x, double observed_y, double std_deviation) + : observed_x_(observed_x) + , observed_y_(observed_y) + , scale_(1.0 / std_deviation) + {} + + template + bool operator()(const T* const camera, + const T* const shot, + const T* const point, + T* residuals) const { + T camera_point[3]; + WorldToCameraCoordinates(shot, point, camera_point); + + if (camera_point[2] <= T(0.0)) { + residuals[0] = residuals[1] = T(99.0); + return true; + } + + T predicted[2]; + FisheyeProject(camera, camera_point, predicted); + + // The error is the difference between the predicted and observed position. + residuals[0] = T(scale_) * (predicted[0] - T(observed_x_)); + residuals[1] = T(scale_) * (predicted[1] - T(observed_y_)); + + return true; + } + + double observed_x_; + double observed_y_; + double scale_; +}; + struct EquirectangularReprojectionError { EquirectangularReprojectionError(double observed_x, double observed_y, double std_deviation) : scale_(1.0 / std_deviation) @@ -515,6 +588,10 @@ class BundleAdjuster { return *(BAPerspectiveCamera *)cameras_[id].get(); } + BAFisheyeCamera GetFisheyeCamera(const std::string &id) { + return *(BAFisheyeCamera *)cameras_[id].get(); + } + BAEquirectangularCamera GetEquirectangularCamera(const std::string &id) { return *(BAEquirectangularCamera *)cameras_[id].get(); } @@ -548,6 +625,27 @@ class BundleAdjuster { c.k2_prior = k2_prior; } + void AddFisheyeCamera( + const std::string &id, + double focal, + double k1, + double k2, + double focal_prior, + double k1_prior, + double k2_prior, + bool constant) { + cameras_[id] = std::unique_ptr(new BAFisheyeCamera()); + BAFisheyeCamera &c = static_cast(*cameras_[id]); + c.id = id; + c.parameters[BA_CAMERA_FOCAL] = focal; + c.parameters[BA_CAMERA_K1] = k1; + c.parameters[BA_CAMERA_K2] = k2; + c.constant = constant; + c.focal_prior = focal_prior; + c.k1_prior = k1_prior; + c.k2_prior = k2_prior; + } + void AddEquirectangularCamera( const std::string &id) { cameras_[id] = std::unique_ptr(new BAEquirectangularCamera()); @@ -761,6 +859,13 @@ class BundleAdjuster { problem.SetParameterBlockConstant(c.parameters); break; } + case BA_FISHEYE_CAMERA: + { + BAFisheyeCamera &c = static_cast(*i.second); + problem.AddParameterBlock(c.parameters, BA_CAMERA_NUM_PARAMS); + problem.SetParameterBlockConstant(c.parameters); + break; + } case BA_EQUIRECTANGULAR_CAMERA: // No parameters for now break; @@ -799,6 +904,22 @@ class BundleAdjuster { observation.point->coordinates); break; } + case BA_FISHEYE_CAMERA: + { + BAFisheyeCamera &c = static_cast(*observation.camera); + ceres::CostFunction* cost_function = + new ceres::AutoDiffCostFunction( + new FisheyeReprojectionError(observation.coordinates[0], + observation.coordinates[1], + reprojection_error_sd_)); + + problem.AddResidualBlock(cost_function, + loss, + c.parameters, + observation.shot->parameters, + observation.point->coordinates); + break; + } case BA_EQUIRECTANGULAR_CAMERA: { BAEquirectangularCamera &c = static_cast(*observation.camera); @@ -878,6 +999,11 @@ class BundleAdjuster { observation.shot->parameters); break; } + case BA_FISHEYE_CAMERA: + { + std::cerr << "NotImplemented: GCP for fisheye cameras\n"; + break; + } case BA_EQUIRECTANGULAR_CAMERA: { ceres::CostFunction* cost_function = @@ -912,6 +1038,21 @@ class BundleAdjuster { c.parameters); break; } + case BA_FISHEYE_CAMERA: + { + BAFisheyeCamera &c = static_cast(*i.second); + + ceres::CostFunction* cost_function = + new ceres::AutoDiffCostFunction( + new InternalParametersPriorError(c.focal_prior, focal_prior_sd_, + c.k1_prior, k1_sd_, + c.k2_prior, k2_sd_)); + + problem.AddResidualBlock(cost_function, + NULL, + c.parameters); + break; + } case BA_EQUIRECTANGULAR_CAMERA: break; } @@ -1005,6 +1146,23 @@ class BundleAdjuster { std::max(observations_[i].point->reprojection_error, error); break; } + case BA_FISHEYE_CAMERA: + { + BAFisheyeCamera &c = static_cast(*observations_[i].camera); + + FisheyeReprojectionError pre(observations_[i].coordinates[0], + observations_[i].coordinates[1], + 1.0); + double residuals[2]; + pre(c.parameters, + observations_[i].shot->parameters, + observations_[i].point->coordinates, + residuals); + double error = sqrt(residuals[0] * residuals[0] + residuals[1] * residuals[1]); + observations_[i].point->reprojection_error = + std::max(observations_[i].point->reprojection_error, error); + break; + } case BA_EQUIRECTANGULAR_CAMERA: { BAEquirectangularCamera &c = static_cast(*observations_[i].camera); diff --git a/opensfm/src/csfm.cc b/opensfm/src/csfm.cc index d3110852d..c7223a7db 100644 --- a/opensfm/src/csfm.cc +++ b/opensfm/src/csfm.cc @@ -83,10 +83,12 @@ BOOST_PYTHON_MODULE(csfm) { class_("BundleAdjuster") .def("run", &BundleAdjuster::Run) .def("get_perspective_camera", &BundleAdjuster::GetPerspectiveCamera) + .def("get_fisheye_camera", &BundleAdjuster::GetFisheyeCamera) .def("get_equirectangular_camera", &BundleAdjuster::GetEquirectangularCamera) .def("get_shot", &BundleAdjuster::GetShot) .def("get_point", &BundleAdjuster::GetPoint) .def("add_perspective_camera", &BundleAdjuster::AddPerspectiveCamera) + .def("add_fisheye_camera", &BundleAdjuster::AddFisheyeCamera) .def("add_equirectangular_camera", &BundleAdjuster::AddEquirectangularCamera) .def("add_shot", &BundleAdjuster::AddShot) .def("add_point", &BundleAdjuster::AddPoint) @@ -119,6 +121,15 @@ BOOST_PYTHON_MODULE(csfm) { .def_readwrite("id", &BAPerspectiveCamera::id) ; + class_("BAFisheyeCamera") + .add_property("focal", &BAFisheyeCamera::GetFocal, &BAFisheyeCamera::SetFocal) + .add_property("k1", &BAFisheyeCamera::GetK1, &BAFisheyeCamera::SetK1) + .add_property("k2", &BAFisheyeCamera::GetK2, &BAFisheyeCamera::SetK2) + .def_readwrite("constant", &BAFisheyeCamera::constant) + .def_readwrite("focal_prior", &BAFisheyeCamera::focal_prior) + .def_readwrite("id", &BAFisheyeCamera::id) + ; + class_("BAShot") .add_property("rx", &BAShot::GetRX, &BAShot::SetRX) .add_property("ry", &BAShot::GetRY, &BAShot::SetRY) @@ -154,6 +165,7 @@ BOOST_PYTHON_MODULE(csfm) { .def("set_min_patch_sd", &csfm::DepthmapEstimatorWrapper::SetMinPatchSD) .def("add_view", &csfm::DepthmapEstimatorWrapper::AddView) .def("compute_patch_match", &csfm::DepthmapEstimatorWrapper::ComputePatchMatch) + .def("compute_patch_match_sample", &csfm::DepthmapEstimatorWrapper::ComputePatchMatchSample) .def("compute_brute_force", &csfm::DepthmapEstimatorWrapper::ComputeBruteForce) ; diff --git a/opensfm/src/depthmap.cc b/opensfm/src/depthmap.cc index c39cc0dc2..debe28722 100644 --- a/opensfm/src/depthmap.cc +++ b/opensfm/src/depthmap.cc @@ -1,7 +1,7 @@ +#include #include - namespace csfm { bool IsInsideImage(const cv::Mat &image, int i, int j) { @@ -135,7 +135,6 @@ float UniformRand(float a, float b) { return a + (b - a) * float(rand()) / RAND_MAX; } - class DepthmapEstimator { public: DepthmapEstimator() @@ -145,6 +144,8 @@ class DepthmapEstimator { , num_depth_planes_(50) , patchmatch_iterations_(3) , min_patch_variance_(5 * 5) + , rng_{std::random_device{}()} + , uni_(0, 0) {} void AddView(const double *pK, @@ -160,6 +161,10 @@ class DepthmapEstimator { Qs_.emplace_back(Rs_.back() * Rs_.front().t()); as_.emplace_back(Qs_.back() * ts_.front() - ts_.back()); images_.emplace_back(cv::Mat(height, width, CV_8U, (void *)pimage).clone()); + std::size_t size = images_.size(); + int a = (size > 1) ? 1 : 0; + int b = (size > 1) ? size - 1 : 0; + uni_.param(std::uniform_int_distribution::param_type(a, b)); } void SetDepthRange(double min_depth, double max_depth, int num_depth_planes) { @@ -176,10 +181,8 @@ class DepthmapEstimator { min_patch_variance_ = sd * sd; } - void ComputeBruteForce(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score) { - *best_depth = cv::Mat(images_[0].rows, images_[0].cols, CV_32F, 0.0f); - *best_plane = cv::Mat(images_[0].rows, images_[0].cols, CV_32FC3, 0.0f); - *best_score = cv::Mat(images_[0].rows, images_[0].cols, CV_32F, 0.0f); + void ComputeBruteForce(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr) { + AssignMatrices(best_depth, best_plane, best_score, best_nghbr); int hpz = (patch_size_ - 1) / 2; for (int i = hpz; i < best_depth->rows - hpz; ++i) { @@ -188,50 +191,68 @@ class DepthmapEstimator { float depth = 1 / (1 / min_depth_ + d * (1 / max_depth_ - 1 / min_depth_) / (num_depth_planes_ - 1)); cv::Vec3f normal(0, 0, -1); cv::Vec3f plane = PlaneFromDepthAndNormal(j, i, Ks_[0], depth, normal); - CheckPlaneCandidate(best_depth, best_plane, best_score, i, j, plane); + CheckPlaneCandidate(best_depth, best_plane, best_score, best_nghbr, i, j, plane); } } } } - void ComputePatchMatch(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score) { - *best_depth = cv::Mat(images_[0].rows, images_[0].cols, CV_32F, 0.0f); - *best_plane = cv::Mat(images_[0].rows, images_[0].cols, CV_32FC3, 0.0f); - *best_score = cv::Mat(images_[0].rows, images_[0].cols, CV_32F, 0.0f); + void ComputePatchMatch(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr) { + AssignMatrices(best_depth, best_plane, best_score, best_nghbr); + RandomInitialization(best_depth, best_plane, best_score, best_nghbr, false); + ComputeIgnoreMask(best_depth, best_plane, best_score, best_nghbr); - RandomInitialization(best_depth, best_plane, best_score); - ComputeIgnoreMask(best_depth, best_plane, best_score); + for (int i = 0; i < patchmatch_iterations_; ++i) { + PatchMatchForwardPass(best_depth, best_plane, best_score, best_nghbr, false); + PatchMatchBackwardPass(best_depth, best_plane, best_score, best_nghbr, false); + } + } + + void ComputePatchMatchSample(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr) { + AssignMatrices(best_depth, best_plane, best_score, best_nghbr); + RandomInitialization(best_depth, best_plane, best_score, best_nghbr, true); + ComputeIgnoreMask(best_depth, best_plane, best_score, best_nghbr); for (int i = 0; i < patchmatch_iterations_; ++i) { - PatchMatchForwardPass(best_depth, best_plane, best_score); - PatchMatchBackwardPass(best_depth, best_plane, best_score); + PatchMatchForwardPass(best_depth, best_plane, best_score, best_nghbr, true); + PatchMatchBackwardPass(best_depth, best_plane, best_score, best_nghbr, true); } } - void RandomInitialization(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score) { + void AssignMatrices(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr) { + *best_depth = cv::Mat(images_[0].rows, images_[0].cols, CV_32F, 0.0f); + *best_plane = cv::Mat(images_[0].rows, images_[0].cols, CV_32FC3, 0.0f); + *best_score = cv::Mat(images_[0].rows, images_[0].cols, CV_32F, 0.0f); + *best_nghbr = cv::Mat(images_[0].rows, images_[0].cols, CV_32S, cv::Scalar(0)); + } + + void RandomInitialization(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr, bool sample) { int hpz = (patch_size_ - 1) / 2; for (int i = hpz; i < best_depth->rows - hpz; ++i) { for (int j = hpz; j < best_depth->cols - hpz; ++j) { float depth = UniformRand(min_depth_, max_depth_); cv::Vec3f normal(UniformRand(-1, 1), UniformRand(-1, 1), -1); cv::Vec3f plane = PlaneFromDepthAndNormal(j, i, Ks_[0], depth, normal); - float score = ComputePlaneScore(i, j, plane); - best_depth->at(i, j) = depth; - best_plane->at(i, j) = plane; - best_score->at(i, j) = score; + int nghbr; + float score; + if (sample) { + nghbr = uni_(rng_); + score = ComputePlaneImageScore(i, j, plane, nghbr); + } else { + ComputePlaneScore(i, j, plane, &score, &nghbr); + } + AssignPixel(best_depth, best_plane, best_score, best_nghbr, i, j, depth, plane, score, nghbr); } } } - void ComputeIgnoreMask(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score) { + void ComputeIgnoreMask(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr) { int hpz = (patch_size_ - 1) / 2; for (int i = hpz; i < best_depth->rows - hpz; ++i) { for (int j = hpz; j < best_depth->cols - hpz; ++j) { float variance = PatchVariance(i, j); if (variance < min_patch_variance_) { - best_depth->at(i, j) = 0; - best_plane->at(i, j) = cv::Vec3f(0, 0, 0); - best_score->at(i, j) = 0; + AssignPixel(best_depth, best_plane, best_score, best_nghbr, i, j, 0.0f, cv::Vec3f(0, 0, 0), 0.0f, 0); } } } @@ -250,43 +271,61 @@ class DepthmapEstimator { } - void PatchMatchForwardPass(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score) { - int neighbors[2][2] = {{-1, 0}, {0, -1}}; + void PatchMatchForwardPass(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr, + bool sample) { + int adjacent[2][2] = {{-1, 0}, {0, -1}}; int hpz = (patch_size_ - 1) / 2; for (int i = hpz; i < best_depth->rows - hpz; ++i) { for (int j = hpz; j < best_depth->cols - hpz; ++j) { - PatchMatchUpdatePixel(best_depth, best_plane, best_score, i, j, neighbors); + PatchMatchUpdatePixel(best_depth, best_plane, best_score, best_nghbr, i, j, adjacent, sample); } } } - void PatchMatchBackwardPass(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score) { - int neighbors[2][2] = {{0, 1}, {1, 0}}; + void PatchMatchBackwardPass(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr, + bool sample) { + int adjacent[2][2] = {{0, 1}, {1, 0}}; int hpz = (patch_size_ - 1) / 2; for (int i = best_depth->rows - hpz - 1; i >= hpz; --i) { for (int j = best_depth->cols - hpz - 1; j >= hpz; --j) { - PatchMatchUpdatePixel(best_depth, best_plane, best_score, i, j, neighbors); + PatchMatchUpdatePixel(best_depth, best_plane, best_score, best_nghbr, i, j, adjacent, sample); } } } - void PatchMatchUpdatePixel(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, + void PatchMatchUpdatePixel(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr, int i, int j, - int neighbors[2][2]) { + int adjacent[2][2], + bool sample) { // Ignore pixels with depth == 0. if (best_depth->at(i, j) == 0.0f) { return; } - // Check neighbor's planes. + // Check neighbors and their planes for adjacent pixels. for (int k = 0; k < 2; ++k) { - cv::Vec3f plane = best_plane->at(i + neighbors[k][0], j + neighbors[k][1]); - CheckPlaneCandidate(best_depth, best_plane, best_score, i, j, plane); + int i_adjacent = i + adjacent[k][0]; + int j_adjacent = j + adjacent[k][1]; + + // Do not propagate ignored adjacent pixels. + if (best_depth->at(i_adjacent, j_adjacent) == 0.0f) { + continue; + } + + cv::Vec3f plane = best_plane->at(i_adjacent, j_adjacent); + + if (sample) { + int nghbr = best_nghbr->at(i_adjacent, j_adjacent); + CheckPlaneImageCandidate(best_depth, best_plane, best_score, best_nghbr, i, j, plane, nghbr); + } else { + CheckPlaneCandidate(best_depth, best_plane, best_score, best_nghbr, i, j, plane); + } } - // Check random planes. + // Check random planes for current neighbor. float depth_range = (1 / max_depth_ - 1 / min_depth_) / 20; float normal_range = 0.5; + int current_nghbr = best_nghbr->at(i, j); for (int k = 0; k < 6; ++k) { float current_depth = best_depth->at(i, j); float depth = 1 / (1 / current_depth + UniformRand(-depth_range, depth_range)); @@ -297,32 +336,68 @@ class DepthmapEstimator { -1.0f); cv::Vec3f plane = PlaneFromDepthAndNormal(j, i, Ks_[0], depth, normal); - CheckPlaneCandidate(best_depth, best_plane, best_score, i, j, plane); + if (sample) { + CheckPlaneImageCandidate(best_depth, best_plane, best_score, best_nghbr, i, j, plane, current_nghbr); + } else { + CheckPlaneCandidate(best_depth, best_plane, best_score, best_nghbr, i, j, plane); + } depth_range *= 0.5; normal_range *= 0.5; } + + if (!sample || images_.size() <= 2) { + return; + } + + // Check random other neighbor for current plane. + int other_nghbr = uni_(rng_); + while (other_nghbr == current_nghbr) { + other_nghbr = uni_(rng_); + } + + cv::Vec3f plane = best_plane->at(i, j); + CheckPlaneImageCandidate(best_depth, best_plane, best_score, best_nghbr, i, j, plane, other_nghbr); } - void CheckPlaneCandidate(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, + void CheckPlaneCandidate(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr, int i, int j, const cv::Vec3f &plane) { - float score = ComputePlaneScore(i, j, plane); + float score; + int nghbr; + ComputePlaneScore(i, j, plane, &score, &nghbr); if (score > best_score->at(i, j)) { - best_score->at(i, j) = score; - best_plane->at(i, j) = plane; - best_depth->at(i, j) = DepthOfPlaneBackprojection(j, i, Ks_[0], plane); + float depth = DepthOfPlaneBackprojection(j, i, Ks_[0], plane); + AssignPixel(best_depth, best_plane, best_score, best_nghbr, i, j, depth, plane, score, nghbr); + } + } + + void CheckPlaneImageCandidate(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr, + int i, int j, const cv::Vec3f &plane, int nghbr) { + float score = ComputePlaneImageScore(i, j, plane, nghbr); + if (score > best_score->at(i, j)) { + float depth = DepthOfPlaneBackprojection(j, i, Ks_[0], plane); + AssignPixel(best_depth, best_plane, best_score, best_nghbr, i, j, depth, plane, score, nghbr); } } - float ComputePlaneScore(int i, int j, const cv::Vec3f &plane) { - float best_score = -1.0f; + void AssignPixel(cv::Mat *best_depth, cv::Mat *best_plane, cv::Mat *best_score, cv::Mat *best_nghbr, + int i, int j, const float depth, const cv::Vec3f &plane, const float score, const int nghbr) { + best_depth->at(i, j) = depth; + best_plane->at(i, j) = plane; + best_score->at(i, j) = score; + best_nghbr->at(i, j) = nghbr; + } + + void ComputePlaneScore(int i, int j, const cv::Vec3f &plane, float *score, int *nghbr) { + *score = -1.0f; + *nghbr = 0; for (int other = 1; other < images_.size(); ++other) { - float score = ComputePlaneImageScore(i, j, plane, other); - if (score > best_score) { - best_score = score; + float image_score = ComputePlaneImageScore(i, j, plane, other); + if (image_score > *score) { + *score = image_score; + *nghbr = other; } } - return best_score; } float ComputePlaneImageScoreUnoptimized(int i, int j, @@ -390,6 +465,8 @@ class DepthmapEstimator { int num_depth_planes_; int patchmatch_iterations_; float min_patch_variance_; + std::mt19937 rng_; + std::uniform_int_distribution uni_; }; diff --git a/opensfm/src/depthmap_wrapper.cc b/opensfm/src/depthmap_wrapper.cc index 8fc0694cd..c027ba359 100644 --- a/opensfm/src/depthmap_wrapper.cc +++ b/opensfm/src/depthmap_wrapper.cc @@ -32,26 +32,34 @@ class DepthmapEstimatorWrapper { } bp::object ComputePatchMatch() { - cv::Mat depth, plane, score; - de_.ComputePatchMatch(&depth, &plane, &score); - return ComputeReturnValues(depth, plane, score); + cv::Mat depth, plane, score, nghbr; + de_.ComputePatchMatch(&depth, &plane, &score, &nghbr); + return ComputeReturnValues(depth, plane, score, nghbr); + } + + bp::object ComputePatchMatchSample() { + cv::Mat depth, plane, score, nghbr; + de_.ComputePatchMatchSample(&depth, &plane, &score, &nghbr); + return ComputeReturnValues(depth, plane, score, nghbr); } bp::object ComputeBruteForce() { - cv::Mat depth, plane, score; - de_.ComputeBruteForce(&depth, &plane, &score); - return ComputeReturnValues(depth, plane, score); + cv::Mat depth, plane, score, nghbr; + de_.ComputeBruteForce(&depth, &plane, &score, &nghbr); + return ComputeReturnValues(depth, plane, score, nghbr); } bp::object ComputeReturnValues(const cv::Mat &depth, const cv::Mat &plane, - const cv::Mat &score) { + const cv::Mat &score, + const cv::Mat &nghbr) { bp::list retn; npy_intp shape[2] = {depth.rows, depth.cols}; npy_intp plane_shape[3] = {depth.rows, depth.cols, 3}; retn.append(bpn_array_from_data(2, shape, depth.ptr(0))); retn.append(bpn_array_from_data(3, plane_shape, plane.ptr(0))); retn.append(bpn_array_from_data(2, shape, score.ptr(0))); + retn.append(bpn_array_from_data(2, shape, nghbr.ptr(0))); return retn; } diff --git a/opensfm/src/openmvs_exporter.h b/opensfm/src/openmvs_exporter.h index 8dd7108d7..ce277ff61 100644 --- a/opensfm/src/openmvs_exporter.h +++ b/opensfm/src/openmvs_exporter.h @@ -73,7 +73,7 @@ class OpenMVSExporter { } void Export(std::string filename) { - ARCHIVE::SerializeSave(scene_, filename); + MVS::ARCHIVE::SerializeSave(scene_, filename); } private: diff --git a/opensfm/src/third_party/openmvs/Interface.h b/opensfm/src/third_party/openmvs/Interface.h index 5fa6dc4b8..eefe98aec 100644 --- a/opensfm/src/third_party/openmvs/Interface.h +++ b/opensfm/src/third_party/openmvs/Interface.h @@ -9,6 +9,14 @@ // D E F I N E S /////////////////////////////////////////////////// +#define MVSI_PROJECT_ID "MVSI" // identifies the project stream +#define MVSI_PROJECT_VER ((uint32_t)2) // identifies the version of a project stream + +// set a default namespace name is none given +#ifndef _INTERFACE_NAMESPACE +#define _INTERFACE_NAMESPACE MVS +#endif + // uncomment to enable custom OpenCV data types // (should be uncommented if OpenCV is not available) #if !defined(_USE_OPENCV) && !defined(_USE_CUSTOM_CV) @@ -27,46 +35,58 @@ namespace cv { // simple cv::Matx -template +template class Matx { public: - typedef _Tp Type; + typedef Type value_type; inline Matx() {} #ifdef _USE_EIGEN EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Type,m*n) typedef Eigen::Matrix1?Eigen::RowMajor:Eigen::Default)> EMat; typedef Eigen::Map CEMatMap; typedef Eigen::Map EMatMap; - inline Matx(const EMat& rhs) { operator EMatMap () = rhs; } - inline Matx& operator = (const EMat& rhs) { operator EMatMap () = rhs; return *this; } + template + inline Matx(const Eigen::EigenBase& rhs) { operator EMatMap () = rhs; } + template + inline Matx& operator = (const Eigen::EigenBase& rhs) { operator EMatMap () = rhs; return *this; } inline operator CEMatMap() const { return CEMatMap((const Type*)val); } inline operator EMatMap () { return EMatMap((Type*)val); } #endif + static Matx eye() { + Matx M; + memset(M.val, 0, sizeof(Type)*m*n); + const int shortdim(m < n ? m : n); + for (int i = 0; i < shortdim; i++) + M(i,i) = 1; + return M; + } Type operator()(int r, int c) const { return val[r*n+c]; } Type& operator()(int r, int c) { return val[r*n+c]; } public: - _Tp val[m*n]; + Type val[m*n]; }; // simple cv::Matx -template +template class Point3_ { public: - typedef _Tp Type; + typedef Type value_type; inline Point3_() {} #ifdef _USE_EIGEN EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Type,3) typedef Eigen::Matrix EVec; typedef Eigen::Map EVecMap; - inline Point3_(const EVec& rhs) { operator EVecMap () = rhs; } - inline Point3_& operator = (const EVec& rhs) { operator EVecMap () = rhs; return *this; } + template + inline Point3_(const Eigen::EigenBase& rhs) { operator EVecMap () = rhs; } + template + inline Point3_& operator = (const Eigen::EigenBase& rhs) { operator EVecMap () = rhs; return *this; } inline operator const EVecMap () const { return EVecMap((Type*)this); } inline operator EVecMap () { return EVecMap((Type*)this); } #endif public: - _Tp x, y, z; + Type x, y, z; }; } // namespace cv @@ -74,62 +94,110 @@ class Point3_ /*----------------------------------------------------------------*/ +namespace _INTERFACE_NAMESPACE { + // custom serialization namespace ARCHIVE { -struct ArchiveSave; -struct ArchiveLoad; - -template -bool Save(ArchiveSave& a, const _Tp& obj) { - const_cast<_Tp&>(obj).serialize(a, 0); - return true; -} -template -bool Load(ArchiveLoad& a, _Tp& obj) { - obj.serialize(a, 0); - return true; -} - - // Basic serialization types struct ArchiveSave { std::ostream& stream; - ArchiveSave(std::ostream& _stream) : stream(_stream) {} + uint32_t version; + ArchiveSave(std::ostream& _stream, uint32_t _version) + : stream(_stream), version(_version) {} template - ArchiveSave& operator & (const _Tp& obj) { - Save(*this, obj); - return *this; - } + ArchiveSave& operator & (const _Tp& obj); }; struct ArchiveLoad { std::istream& stream; - ArchiveLoad(std::istream& _stream) : stream(_stream) {} + uint32_t version; + ArchiveLoad(std::istream& _stream, uint32_t _version) + : stream(_stream), version(_version) {} template - ArchiveLoad& operator & (_Tp& obj) { - Load(*this, obj); - return *this; - } + ArchiveLoad& operator & (_Tp& obj); }; +template +bool Save(ArchiveSave& a, const _Tp& obj) { + const_cast<_Tp&>(obj).serialize(a, a.version); + return true; +} +template +bool Load(ArchiveLoad& a, _Tp& obj) { + obj.serialize(a, a.version); + return true; +} + +template +ArchiveSave& ArchiveSave::operator & (const _Tp& obj) { + Save(*this, obj); + return *this; +} +template +ArchiveLoad& ArchiveLoad::operator & (_Tp& obj) { + Load(*this, obj); + return *this; +} // Main exporter & importer template -bool SerializeSave(const _Tp& obj, const std::string& fileName) { +bool SerializeSave(const _Tp& obj, const std::string& fileName, uint32_t version=MVSI_PROJECT_VER) { + // open the output stream std::ofstream stream(fileName, std::ofstream::binary); if (!stream.is_open()) return false; - ARCHIVE::ArchiveSave serializer(stream); + // write header + if (version > 0) { + // save project ID + stream.write(MVSI_PROJECT_ID, 4); + // save project version + stream.write((const char*)&version, sizeof(uint32_t)); + // reserve some bytes + const uint32_t reserved(0); + stream.write((const char*)&reserved, sizeof(uint32_t)); + } + // serialize out the current state + ARCHIVE::ArchiveSave serializer(stream, version); serializer & obj; return true; } template -bool SerializeLoad(_Tp& obj, const std::string& fileName) { +bool SerializeLoad(_Tp& obj, const std::string& fileName, uint32_t* pVersion=NULL) { + // open the input stream std::ifstream stream(fileName, std::ifstream::binary); if (!stream.is_open()) return false; - ARCHIVE::ArchiveLoad serializer(stream); + // read header + uint32_t version(0); + // load project header ID + char szHeader[4]; + stream.read(szHeader, 4); + if (!stream) + return false; + if (strncmp(szHeader, MVSI_PROJECT_ID, 4) != 0) { + // try to load as the first version that didn't have a header + const size_t size(fileName.size()); + if (size <= 4) + return false; + std::string ext(fileName.substr(size-4)); + std::transform(ext.begin(), ext.end(), ext.begin(), ::towlower); + if (ext != ".mvs") + return false; + stream.seekg(0, std::ifstream::beg); + } else { + // load project version + stream.read((char*)&version, sizeof(uint32_t)); + if (!stream || version > MVSI_PROJECT_VER) + return false; + // skip reserved bytes + uint32_t reserved; + stream.read((char*)&reserved, sizeof(uint32_t)); + } + // serialize in the current state + ARCHIVE::ArchiveLoad serializer(stream, version); serializer & obj; + if (pVersion) + *pVersion = version; return true; } @@ -221,17 +289,16 @@ bool Load(ArchiveLoad& a, std::vector<_Tp>& v) { /*----------------------------------------------------------------*/ -namespace MVS { - // interface used to export/import MVS input data; -// MAX(width,height) is used for normalization +// - MAX(width,height) is used for normalization +// - row-major order is used for storing the matrices struct Interface { typedef cv::Point3_ Pos3f; typedef cv::Point3_ Pos3d; typedef cv::Matx Mat33d; - typedef uint8_t Color; - typedef cv::Point3_ Col3; // x=B, y=G, z=R + typedef cv::Matx Mat44d; + typedef cv::Point3_ Col3; // x=B, y=G, z=R /*----------------------------------------------------------------*/ // structure describing a mobile platform with cameras attached to it @@ -239,13 +306,22 @@ struct Interface // structure describing a camera mounted on a platform struct Camera { std::string name; // camera's name - Mat33d K; // camera's normalized intrinsics matrix + uint32_t width, height; // image resolution in pixels for all images sharing this camera (optional) + Mat33d K; // camera's intrinsics matrix (normalized if image resolution not specified) Mat33d R; // camera's rotation matrix relative to the platform Pos3d C; // camera's translation vector relative to the platform + Camera() : width(0), height(0) {} + bool HasResolution() const { return width > 0 && height > 0; } + bool IsNormalized() const { return !HasResolution(); } + template - void serialize(Archive& ar, const unsigned int /*version*/) { + void serialize(Archive& ar, const unsigned int version) { ar & name; + if (version > 0) { + ar & width; + ar & height; + } ar & K; ar & R; ar & C; @@ -325,8 +401,37 @@ struct Interface typedef std::vector VertexArr; /*----------------------------------------------------------------*/ + // structure describing a 3D line + struct Line { + // structure describing one view for a given 3D feature + struct View { + uint32_t imageID; // image ID corresponding to this view + float confidence; // view's confidence (0 - not available) + + template + void serialize(Archive& ar, const unsigned int /*version*/) { + ar & imageID; + ar & confidence; + } + }; + typedef std::vector ViewArr; + + Pos3f pt1; // 3D line segment end-point + Pos3f pt2; // 3D line segment end-point + ViewArr views; // list of all available views for this 3D feature + + template + void serialize(Archive& ar, const unsigned int /*version*/) { + ar & pt1; + ar & pt2; + ar & views; + } + }; + typedef std::vector LineArr; + /*----------------------------------------------------------------*/ + // structure describing a 3D point's normal (optional) - struct VertexNormal { + struct Normal { Pos3f n; // 3D feature normal template @@ -334,11 +439,11 @@ struct Interface ar & n; } }; - typedef std::vector VertexNormalArr; + typedef std::vector NormalArr; /*----------------------------------------------------------------*/ // structure describing a 3D point's color (optional) - struct VertexColor { + struct Color { Col3 c; // 3D feature color template @@ -346,26 +451,40 @@ struct Interface ar & c; } }; - typedef std::vector VertexColorArr; + typedef std::vector ColorArr; /*----------------------------------------------------------------*/ PlatformArr platforms; // array of platforms ImageArr images; // array of images VertexArr vertices; // array of reconstructed 3D points - VertexNormalArr verticesNormal; // array of reconstructed 3D points' normal (optional) - VertexColorArr verticesColor; // array of reconstructed 3D points' color (optional) + NormalArr verticesNormal; // array of reconstructed 3D points' normal (optional) + ColorArr verticesColor; // array of reconstructed 3D points' color (optional) + LineArr lines; // array of reconstructed 3D lines + NormalArr linesNormal; // array of reconstructed 3D lines' normal (optional) + ColorArr linesColor; // array of reconstructed 3D lines' color (optional) + Mat44d transform; // transformation used to convert from absolute to relative coordinate system (optional) + + Interface() : transform(Mat44d::eye()) {} template - void serialize(Archive& ar, const unsigned int /*version*/) { + void serialize(Archive& ar, const unsigned int version) { ar & platforms; ar & images; ar & vertices; ar & verticesNormal; ar & verticesColor; + if (version > 0) { + ar & lines; + ar & linesNormal; + ar & linesColor; + if (version > 1) { + ar & transform; + } + } } }; /*----------------------------------------------------------------*/ -} // namespace MVS +} // namespace _INTERFACE_NAMESPACE #endif // _INTERFACE_MVS_H_ diff --git a/opensfm/src/third_party/openmvs/README.opensfm b/opensfm/src/third_party/openmvs/README.opensfm index 82c73570d..036a4e399 100644 --- a/opensfm/src/third_party/openmvs/README.opensfm +++ b/opensfm/src/third_party/openmvs/README.opensfm @@ -1,7 +1,7 @@ Project: openMVS URL: https://github.com/cdcseacave/openMVS License: BSD -Version: 0.7, 2016-04-06 +Version: a3b360016660a1397f6eb6c070c2c19bbb4c7590, 2017-05-01 Local modifications: * Imported only the Interface.h file diff --git a/opensfm/src/third_party/vlfeat/vl/host.c b/opensfm/src/third_party/vlfeat/vl/host.c index ef97cdcae..499db6484 100644 --- a/opensfm/src/third_party/vlfeat/vl/host.c +++ b/opensfm/src/third_party/vlfeat/vl/host.c @@ -441,6 +441,7 @@ _vl_cpuid (vl_int32* info, int function) #endif +#if defined(HAS_CPUID) void _vl_x86cpu_info_init (VlX86CpuInfo *self) { @@ -463,6 +464,7 @@ _vl_x86cpu_info_init (VlX86CpuInfo *self) self->hasAVX = info[2] & (1 << 28) ; } } +#endif char * _vl_x86cpu_info_to_string_copy (VlX86CpuInfo const *self) diff --git a/opensfm/test/test_types.py b/opensfm/test/test_types.py index a074b6e89..c14eef0e7 100644 --- a/opensfm/test/test_types.py +++ b/opensfm/test/test_types.py @@ -1,5 +1,6 @@ import numpy as np +from opensfm import context from opensfm import types """ @@ -101,6 +102,22 @@ def test_perspective_camera_projection(): assert np.allclose(pixel, projected) +def test_fisheye_camera_projection(): + """Test fisheye projection--backprojection loop.""" + if not context.OPENCV3: + return + camera = types.FisheyeCamera() + camera.width = 800 + camera.height = 600 + camera.focal = 0.6 + camera.k1 = -0.1 + camera.k2 = 0.01 + pixel = [0.1, 0.2] + bearing = camera.pixel_bearing(pixel) + projected = camera.project(bearing) + assert np.allclose(pixel, projected) + + def test_spherical_camera_projection(): """Test spherical projection--backprojection loop.""" camera = types.SphericalCamera() diff --git a/opensfm/types.py b/opensfm/types.py index c25f73c6a..48aa247f2 100644 --- a/opensfm/types.py +++ b/opensfm/types.py @@ -238,6 +238,90 @@ def get_K_in_pixel_coordinates(self, width=None, height=None): [0, 0, 1.0]]) +class FisheyeCamera(Camera): + """Define a fisheye camera. + + Attributes: + widht (int): image width. + height (int): image height. + focal (real): estimated focal lenght. + k1 (real): estimated first distortion parameter. + k2 (real): estimated second distortion parameter. + focal_prior (real): prior focal lenght. + k1_prior (real): prior first distortion parameter. + k2_prior (real): prior second distortion parameter. + """ + + def __init__(self): + """Defaut constructor.""" + self.id = None + self.projection_type = 'fisheye' + self.width = None + self.height = None + self.focal = None + self.k1 = None + self.k2 = None + self.focal_prior = None + self.k1_prior = None + self.k2_prior = None + + def project(self, point): + """Project a 3D point in camera coordinates to the image plane.""" + x, y, z = point + l = np.sqrt(x**2 + y**2) + theta = np.arctan2(l, z) + theta_d = theta * (1.0 + theta**2 * (self.k1 + theta**2 * self.k2)) + s = self.focal * theta_d / l + return np.array([s * x, s * y]) + + def pixel_bearing(self, pixel): + """Unit vector pointing to the pixel viewing direction.""" + point = np.asarray(pixel).reshape((1, 1, 2)) + distortion = np.array([self.k1, self.k2, 0., 0.]) + x, y = cv2.fisheye.undistortPoints(point, self.get_K(), distortion).flat + l = np.sqrt(x * x + y * y + 1.0) + return np.array([x / l, y / l, 1.0 / l]) + + def pixel_bearings(self, pixels): + """Unit vector pointing to the pixel viewing directions.""" + points = pixels.reshape((-1, 1, 2)).astype(np.float64) + distortion = np.array([self.k1, self.k2, 0., 0.]) + up = cv2.fisheye.undistortPoints(points, self.get_K(), distortion) + up = up.reshape((-1, 2)) + x = up[:, 0] + y = up[:, 1] + l = np.sqrt(x * x + y * y + 1.0) + return np.column_stack((x / l, y / l, 1.0 / l)) + + def back_project(self, pixel, depth): + """Project a pixel to a fronto-parallel plane at a given depth.""" + bearing = self.pixel_bearing(pixel) + scale = depth / bearing[2] + return scale * bearing + + def get_K(self): + """The calibration matrix.""" + return np.array([[self.focal, 0., 0.], + [0., self.focal, 0.], + [0., 0., 1.]]) + + def get_K_in_pixel_coordinates(self, width=None, height=None): + """The calibration matrix that maps to pixel coordinates. + + Coordinates (0,0) correspond to the center of the top-left pixel, + and (width - 1, height - 1) to the center of bottom-right pixel. + + You can optionally pass the width and height of the image, in case + you are using a resized versior of the original image. + """ + w = width or self.width + h = height or self.height + f = self.focal * max(w, h) + return np.array([[f, 0, 0.5 * (w - 1)], + [0, f, 0.5 * (h - 1)], + [0, 0, 1.0]]) + + class SphericalCamera(Camera): """A spherical camera generating equirectangular projections. diff --git a/requirements.txt b/requirements.txt index 2289cf936..2ef094b72 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,9 +1,10 @@ -exifread -gpxpy -networkx +exifread==2.1.2 +gpxpy==1.1.2 +networkx==1.11 numpy -pytest -python-dateutil -PyYAML -xmltodict -pyproj +pyproj==1.9.5.1 +pytest==3.0.7 +python-dateutil==2.6.0 +PyYAML==3.12 +scipy +xmltodict==0.10.2 diff --git a/setup.cfg b/setup.cfg index ef7356404..ee6817871 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,3 @@ -[pytest] +[tool:pytest] testpaths = opensfm addopts = --doctest-modules diff --git a/setup.py b/setup.py index 35216ea13..b319c8736 100644 --- a/setup.py +++ b/setup.py @@ -31,7 +31,8 @@ def mkdir_p(path): url='https://github.com/mapillary/OpenSfM', author='Mapillary', license='BSD', - packages=['opensfm'], + packages=['opensfm', 'opensfm.commands'], + scripts=['bin/opensfm_run_all', 'bin/opensfm'], package_data={ 'opensfm': ['csfm.so', 'data/sensor_data.json'] }, diff --git a/viewer/reconstruction.html b/viewer/reconstruction.html index 07e255088..f59d7628f 100644 --- a/viewer/reconstruction.html +++ b/viewer/reconstruction.html @@ -169,20 +169,39 @@ - @@ -714,15 +733,14 @@ } function imageVertexShader(cam) { - if (cam.projection_type == 'equirectangular' || cam.projection_type == 'spherical') - return $('#vertexshader_equirectangular').text(); - else - return $('#vertexshader').text(); + return $('#vertexshader').text(); } function imageFragmentShader(cam) { if (cam.projection_type == 'equirectangular' || cam.projection_type == 'spherical') return $('#fragmentshader_equirectangular').text(); + else if (cam.projection_type == 'fisheye') + return $('#fragmentshader_fisheye').text(); else return $('#fragmentshader').text(); }