Photogrammetry: Reconstructing 3D Scenes using Bit-sliced Indexing
Advisor: Swetha Bhattacharya
Internship Location: San Jose State University
Onsite Mentor: Gheorghi Guzun
To conduct operations during a disaster, emergency responders need to understand the infrastructure of the affected area. To map out the area with drones, I chose to employ photogrammetry using OpenSFM, which calculates the center of reference in a series of images to create a point cloud. In my research, I will introduce a data structure called bit-sliced indexing (BSI) to structure from motion, a type of photogrammetry. OpenSFM generally takes a few hours to run on 30 to 60 images, which generates a very photorealistic 3D model. Speeding up the binary operations will make creating 3D models time-efficient and increase the chances of the 3D model being useful for urgent scenarios. Building on my research with Professor Guzun in implementing BSI, I will implement the mathematical operations necessary to run structure from motion. Since all of BSI’s operations should run in linear time, the purpose of this project is to compete with the constant coefficients in NumPy’s also linear runtime. The added benefit of BSI’s compression scheme could also contribute to the resulting runtime.
Project Blog
Week 1: 2D Greetings
Welcome to my 2D blogs about 3D models! Over the next weeks, I’ll be diving into a journey of clouds, not the type in the sky, but point clouds that will be the foundation of my senior project. As a commemoration of the Disney trip my fellow seniors and I have just come back from and an introduction to what point clouds look like, I will include a video of a dense point cloud generated from our Mickey ears in my next blog. I was going to include it in this blog, but I’m currently getting errors with building the OpenSFM library.
Structure from Motion Photogrammetry
Photogrammetry is the art of using 2D images to get information about its 3D qualities, such as depth, shapes, and sizes of objects in the image.
Structure from motion photogrammetry is a technique for generating 3D models by estimating the X, Y, and Z coordinates of each point in a photo by comparing the variation in different photos of the same scene. The resulting 3D model is called a point cloud because it is comprised of a group of points, which don’t form the continuous surfaces seen in animations. Since a lot of images with different perspectives of the same location are required to accurately reconstruct a 3D model, structure from motion can take several hours. When I ran it on 60 images of a 5 square feet plot of land in my backyard, it took around 2-3 hours.
The variation in photos is caused by the phenomenon known as parallax, which is what causes mountains far away to appear stationary even when you’re driving at 90mph on the highway. The Apple website extensively uses parallax to convince you to make your annual purchase of their software update in an iPhone case with one more camera.
Here’s another example of a website that uses parallax to create a visual novel experience: https://www.sbs.com.au/theboat/.
Purpose
I tackled the optimization of photogrammetry with the main motivation of increasing the relevance of using this technology for disaster response. Since emergency responders are not always trained in interpreting graphs generated with techniques like LIDAR and SLAM, the photorealism (maintained colors and textures from pictures) of photogrammetry makes it very useful. Additionally, you don’t need fancy (and expensive) equipment, unlike with LIDAR. Instead, a simple phone camera and cheap drone are enough.
Goal
I’ll be investigating if replacing NumPy vector operations with Bit-sliced Indexing (BSI) operations as the primary data structure can speed up the algorithm used in OpenSFM, an open source implementation of structure from motion photogrammetry. Since both data structures are capable of running many of their operations in linear time, I’ll be investigating how to diminish the constant factor to turn out better results. Speeding up the binary operations will make creating 3D models time-efficient and increase the chances of the 3D model being useful for urgent scenarios.
Currently, it takes a few hours to run OpenSFM on up to 100 photos. My goal is to reduce the time to maybe an hour or so, enough to demo on the day of senior project presentations. (Early advertisement to come to my presentation if you want a chance to get your own 3D model of anything you can bring to school, including yourself!)
Preview of Upcoming Blog
I cloned and partly set up the starter code with the necessary dependencies for running OpenSFM to make the Mickey ears demo in this blog. Next week, I will finish setting it up, and then it will be all about the math because I’m going to have to figure out how to break down the equations into functions in BSI.
Sources:
Week 2: Reconstructing Understanding
This week I’ve been focusing on two major tasks: setting up OpenSFM on my Mac and reading through the OpenSFM library and documentation.
For context, however, we’re gonna start with some previous research done on this topic.
Related Works
To capitalize on photogrammetry’s photorealism, there is an important barrier to overcome, bad atmospheric conditions. Since the detail in photogrammetry relies on contrast, bad weather can make landscapes appear too symmetric to correctly identify matching features across images. Dark lighting causes errors in camera calibration, which can distort the resulting geometry. As discovered in [1], the increase of uniformly spread ground control points (GCPs), points with exactly known coordinates, increased accuracy of the model near GCPs.
The equipment used to capture images for photogrammetry is much cheaper than that of LIDAR. To evaluate the accuracy of 3D models produced by images from a low-cost UAV and digital camera, [2] conducts an experiment over the surface mine Jastrabá in Slovakia. Out of 237 test points, only 3 points failed the minimum requirements for orthophoto maps in national legislation. The results of this experiment not only demonstrates the utility of photogrammetry for mining, but also for other land operations.
Optimizations can take place during the process of taking photos or during the pre-processing of data. An example is the pipeline for color enhancement, image denoising, color-to-gray conversion and image content enrichment proposed in [3]. CBM3D-new and BID, the proposed methods for image denoising and grayscale reduction, respectively, proved to result in denser and more complete dense point clouds. [4] covers trajectory planning to reduce redundant camera angles. [5] shows that a high camera inclination of 20-40 degrees enhances photogrammetry and can improve the spatial structure of errors more significantly than the loss of accuracy from higher flying height.
Within the calculations of the structure from motion process, there are proposed mathematical models that can increase accuracy of feature matching. The scale invariant feature transform (SIFT) from [6] has been fundamental to many point-based methods like the PCA-SIFT [7].
Algorithm Explained
OpenSFM employs incremental structure from motion: starting from a pair of images and iteratively adding the rest of the images to reconstruct the 3D scene.
Step 1: Choosing a pair of images
Good pairs need enough parallax, i.e. the camera was displaced between 2 shots relatively far compared to distance to scene.
To compute how much parallax is “good enough,” OpenSFM tries to fit a rotation-only camera model to the two images. Using the RANSAC algorithm, an iterative method for estimating parameters of a math model by random sampling, the number of outliers that don’t fit the model is determined. If the number of outliers exceeds 30%, the pair is accepted.
Step 2: Initializing a reconstruction
To find the relative camera motion between the two images, at least 5 corresponding points is needed. For each configuration, a set of inliers is found, and the algorithm triangulates, or in other words, calculates the position of the points. The Rodrigues transformation is an efficient algorithm for applying transformation matrices used to check the camera poses with respect to the images. Using the best inliers with the least error, the reconstruction is initialized with the corresponding poses, triangulated with matches from GCP and GPS data, and bundle adjusted to fit the positions.
Step 3: Adding more images
For each new image added, the process is repeated. The algorithm finds the camera position that makes the reconstructed 3D points project to the correct position in the new image via resectioning. If the resectioning works, the image is added, the features of the new image seen in the reconstructed images are triangulated, the reconstruction is bundle adjusted, and all the features are re-triangulated. Finally, the reconstruction is rigidly moved, as in any affine transformations are applied to the whole reconstruction rather than parts, to adjust to GPS or GCP data.
OpenSFM actually has built-in optimization for large datasets. If the GPS data of the images are known, the dataset can be split into subsets to reconstruct submodels. The algorithm can be run in parallel on all of the submodels. After all the submodels have been reconstructed, the GPS data can be used to find the relative positions of the submodels to each other in order to reconstruct the final model.
Matrix Operations
After going through the files, I’ve compiled almost all the matrix operations I will have to override. I say almost all because I only focused on the NumPy operations, while the library has defined several of its own Matrix classes composed of NumPy arrays, so I’ll have to see how that translates. I won’t be missing any major matrix operations, just some implementation differences. Thus, here’s the noncomprehensive list kind of in the order they first appear when running the algorithm (yes, I traced almost the whole library just to get this):
Dot product
Transpose
Addition and subtraction
Normalize
Rodrigues transform (in cv2 library)
Determinant
Inverse
Multiplication and division by scalars
Matrix multiplication
Nonzero
Drop column
Divide by column
Find homography (in cv2 library)
Argmax
Mean
Eigenvalues and eigenvectors
Sort
Less than, greater than
Sum
Add column
Vstack
Is empty
Identity matrix
Extract diagonal
Standard deviation
Concatenate
Reshape
Flatten
Outer product
Pseudo inverse
Angle between vectors
Rotate matrix
Roll elements
Game Plan
As frustrating as it is, there’s no magic trick to get OpenSFM working (still having trouble building, high hopes for the future because setting up is usually the hardest part), so I’ll grind that for another indefinite period of time. In the mean time, I can get started with creating an auxiliary class for BSI to fit into the functions I identified above. Most of the algorithms have implementations, but I need to check how much quantization they do, if at all, since there will be a lot of floats when dealing with pi and conversions between degrees and radians. Not sure if anyone will believe me, but reading through the OpenSFM library was actually pretty fun. While it is tedious, there are some veeeery brief introductions to the algorithms and optimizations used, and with a little extra Google search, I found out about some pretty cool stuff, like the Rodrigues transformations. If you’re into Linear Algebra, I highly suggest taking a look at it.
Sources:
[1] P. Burdziakowski and K. Bobkowska. “UAV Photogrammetry under Poor Lighting Conditions-Accuracy Considerations.” MDPI, Multidisciplinary Digital Publishing Institute, 19 May 2021, www.mdpi.com/1424-8220/21/10/3531.
[2] B. Kršák, P. Blišťan, A. Pauliková, P. Puškárová, L. Kovanič, J. Palková, and V. Zelizňaková. “Use of Low-Cost UAV Photogrammetry to Analyze the Accuracy of a Digital Elevation Model in a Case Study.” Measurement, Elsevier, 12 May 2016, www.sciencedirect.com/science/article/abs/pii/S0263224116301749.
[3] M. Gaiani, F. Remondino, F. I. Apollonio, and A. Ballabeni. “An Advanced Pre-Processing Pipeline to Improve Automated Photogrammetric Reconstructions of Architectural Scenes.” MDPI, Multidisciplinary Digital Publishing Institute, 25 Feb. 2016, www.mdpi.com/2072-4292/8/3/178.
[4] Q. Li, W. Yu, and S. Jiang. “Optimized Views Photogrammetry: Precision Analysis and a Large-Scale Case Study in Qingdao.” arXiv.Org, 24 June 2022, arxiv.org/abs/2206.12216.
[5] W. Dai, R. Qiu, B. Wang, W. Lu, G. Zheng, S.O.Y. Amankwah, and G. Wang. “Enhancing UAV-SFM Photogrammetry for Terrain Modeling from the Perspective of Spatial Structure of Errors.” MDPI, Multidisciplinary Digital Publishing Institute, 31 Aug. 2023, www.mdpi.com/2072-4292/15/17/4305.
[6] D. G. Lowe. Object recognition from local scale-invariant features. In Proc. of the International Conference on Computer Vision ICCV, Corfu, 1999.
[7] Y. Ke and R. Sukthankar. PCA-SIFT: a more distinctive representation for local image descriptors. In Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004., volume 2, pages II–506–II–513 Vol.2, June 2004.
Week 3: Mickey Docks at Shore
Mickey Ears
Good news, Mickey has finally “docked”, and you can see the long-awaited point cloud!
This week, I have found a way to run OpenSFM commands by using a Docker image kindly provided by Dr. Jeff Liu, technical staff member in the Humanitarian Assistance and Disaster Relief group at MIT Lincoln Lab and my instructor from the Remote Sensing for Disaster Response course at MIT Beaverworks. Docker containers are isolated containers of images dictating the instructions for all the code and libraries for running an application. Since Docker containers are independent of your local system, they can easily be shared with others to solve the problem of “it doesn’t work on my computer.” Here’s how I ran a Docker container, which contains an image with the instructions for compiling an environment with OpenSFM already built. I had to bind my local server to the exposed Docker port (-p 8888:8888) in order to be able to open the Jupyter Notebook link generated, which Professor Guzun gave me the command for.
Once I opened the link, the Jupyter Notebook initially has an empty folder “work.” I added a Python notebook inside and another folder called “mickey_with_ears” with a folder called “images” where I uploaded all of my images.
I put the Mickey ears on top of another Mickey because I thought the images had a problem at first. When I uploaded the images, I didn’t make a folder called “images” inside the “mickey_ears” folder and directly uploaded the images to the “mickey_ears” folder. The opensfm command was returning 0 reconstructions, so I thought maybe the images were too symmetric. Then, I looked back at my work from the summer and realized my mistake, but since I already have these images, might as well use them.
Here’s the command I ran to generate the .ply file. This took about an hour for 60 images.
And here’s Mickey!
Sparse cloud (converted to HTML file and opened in browser):
The red, green, and blue dots show the estimated camera location and orientation for each image.
Dense cloud (imported PLY file on MeshLabs):
(You can faintly make out “BISV” “Class of ‘24”)
Surprise slightly squashed two Mobius strip hearts on the back
As you can see, the point cloud looks a bit like a super low-resolution picture where everything is pixelated, but it’s actually a collection of points rendered with colors in 3D.
However, in order to run OpenSFM with BSI (bit-sliced indexing, the data structure I will be using for matrix operations), I still need to be able to build OpenSFM on my computer and later incorporate the build into the Docker image. Thus, allow me to explain the problems I have been encountering.
Building OpenSFM
These past two weeks, I have been having trouble building OpenSFM on my computer, which means that the project is failing to configure the source code to run on my local system (MacOS). The build process also includes the compile time, which is probably the most commonly known part of it.
Earlier in the process, I fixed some issues with mismatched dependencies/version issues. In case you want to build OpenSFM for yourself, I will go through the errors in the order they appear.
#1:
pip install -r requirements
a.
Error:
Solution:
Add .txt to requirements and run
pip install -r requirements.txt
b.
Error:
Solution:
Cython 3.0.0 has an issue where it’s not compatible with pyyaml 5.4.1, so remove the version number from pyyaml in requirements.txt. If you’re working in the terminal, run
open requirements.txt
to open the text editor for the file.
c.
Error:
Solution:
You can also ignore this because you will install opencv with Homebrew later.
#2: export PYTHONPATH=/usr/local/lib/python3.7/site-packages:$PYTHONPATH
In case you don’t know how to find your Python path to the installed packages, look at the output of pip install requirements.txt.
Where it says something like “/Library/Frameworks…site-packages” is where your path is. Your path might not start with “/Library/Frameworks.”
export PYTHONPATH=/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages:$PYTHONPATH
You can include any of these export lines into your .zshrc or .bashrc, or, better yet, include it in your CMakeLists.txt with slightly different syntax (search up SET CMake documentation).
#3: export C_INCLUDE_PATH=/usr/local/include
export CPLUS_INCLUDE_PATH=/usr/local/include
export DYLD_LIBRARY_PATH=$HOME/local/lib64
To find where your libraries are installed by Homebrew, run brew info opencv (or any library you know you installed).
The “/usr/local/Cellar” is where mine was installed, so I would run
export C_INCLUDE_PATH=/usr/local/Cellar
export CPLUS_INCLUDE_PATH=/usr/local/Cellar
export DYLD_LIBRARY_PATH=$HOME/local/lib64
#4: python3 setup.py build
a.
Error:
Solution:
What you see above is an error due to adding OpenMP as a library, but clang only works with libomp, not libgomp. You need to comment out these lines in opensfm/src/CMakeLists.txt:
And add
if(APPLE)
if(CMAKE_C_COMPILER_ID MATCHES "Clang")
set(OpenMP_C "${CMAKE_C_COMPILER}" CACHE STRING "" FORCE)
set(OpenMP_C_FLAGS "-fopenmp=libomp -Wno-unused-command-line-argument" CACHE STRING "" FORCE)
set(OpenMP_C_LIB_NAMES "libomp" "libgomp" "libiomp5" CACHE STRING "" FORCE)
set(OpenMP_libomp_LIBRARY ${OpenMP_C_LIB_NAMES} CACHE STRING "" FORCE)
set(OpenMP_libgomp_LIBRARY ${OpenMP_C_LIB_NAMES} CACHE STRING "" FORCE)
set(OpenMP_libiomp5_LIBRARY ${OpenMP_C_LIB_NAMES} CACHE STRING "" FORCE)
endif()
if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(OpenMP_CXX "${CMAKE_CXX_COMPILER}" CACHE STRING "" FORCE)
set(OpenMP_CXX_FLAGS "-fopenmp=libomp -Wno-unused-command-line-argument" CACHE STRING "" FORCE)
set(OpenMP_CXX_LIB_NAMES "libomp" "libgomp" "libiomp5" CACHE STRING "" FORCE)
set(OpenMP_libomp_LIBRARY ${OpenMP_CXX_LIB_NAMES} CACHE STRING "" FORCE)
set(OpenMP_libgomp_LIBRARY ${OpenMP_CXX_LIB_NAMES} CACHE STRING "" FORCE)
set(OpenMP_libiomp5_LIBRARY ${OpenMP_CXX_LIB_NAMES} CACHE STRING "" FORCE)
endif()
endif()
Sources:
https://opensfm.org/docs/building.html
https://stackoverflow.com/questions/48825416/missing-openmp-c-flags-openmp-c-lib-names
Week 4: Breaking Down Computers and Teeth
This week, I have not only started to break down the math operations in BSI but also my computer and my teeth. I got my last two wisdom teeth removed (hoping to never see another numbing needle again) and suffered most of the week, so this week’s blog will mostly be to catch you up on my prior research. As for my computer, after restarting the build process several times, I decided to put the build on hold and pivot my attention to tweaking the BSI library, only to discover that somewhere in my build process, I accidentally changed an important environment variable which now makes my BSI build also fail. In conclusion, until I can figure out which variable I changed, I’m stuck with either using a different computer or just focusing on the math part. Sorry, computer.
Building on Windows
After witnessing countless failures on my MacOS related to OpenMP, I decided to try my hand at building OpenSFM on my Windows computer instead to see if it really is more straightforward.
For those of you wondering why OpenMP is such a big issue, Apple disabled Clang support for OpenMP (probably out of pure malice to ruin my efforts), so you have to go the long way and install a version of LLVM based on your Clang version. The shortest explanation I found has been on this project. I haven’t tried the installation on the website yet because I’m not completely sure how reliable it is, but I will probably investigate it further sometime later.
Now, to Windows.
I first tried installing the Windows Subsystem for Linux (WSL) so that I could have Linux because my external advisor was actually able to build OpenSFM on his Linux computer. However, I ran into this error opening it and haven’t yet had the time to figure out what went wrong because it isn’t one of the named errors in their installation FAQ.
On the regular Windows command prompt, one new error which I did see on MacOS but forgot to mention is this:
On MacOS, I installed all the dependencies it was returning an error for, so I didn’t see this error again. On Windows, I decided to ignore the error to see if it would cause any future problems. So far, I am getting some ModuleNotFoundError’s so the next thing I’m gonna try is fixing those dependency errors.
Bit-Slice Indexing (BSI)
Now that we’re here at the math part, let me run through the basic algorithm behind bit-slice indexing. Indexing is a technique to speed up data retrieval and query processing. The two types I’m going to talk about are bit-slice indexing and bitmap indexing.
Bitmap indexing creates integer indices where the bit values represent the presence of each element. For example:
Bitmap Blue Fruit Vegetable
101 1 0 1
001 0 0 1
010 0 1 0
This makes for very fast evaluation of boolean expressions.
Bit-slice indexing uses bitmaps to encode its values, where the dataset is partitioned into slices of its binary representation. For example:
Number Binary Representation Slice 3 Slice 2 Slice 1
1 001 0 0 1
5 101 1 0 1
7 111 1 1 1
Slice 3 is encoded as bit-string 110 = decimal number 6
Slice 2 is encoded as bit-string 100 = decimal number 4
Slice 1 is encoded as bit-string 111 = decimal number 7
Bit-sliced indexing is preferred for operations like addition and multiplication.
Hybrid BSI
Data compression is the technique of reducing file size and saving storage by encoding large amounts of data using fewer bits.
One such technique is the Word Aligned Hybrid (WAH) compression for bitmap indices. WAH represents bitmaps using two types of words, filler and literal, indicated by a single bit at the beginning of each word. Literal words are chunks of length w - 1, where w is the length of a word on the local computer, in the bitmap that are represented verbatim. Filler words compress consecutive groups of zeros and ones, where one bit is used to represent the fill value (zero or one), and the rest of the w - 2 bits represent the run length (how many times the fill value repeats). For example:
If w = 8, and verbatim words were represented by 0, filler words by 1,
10000101
could be encoded as
11000001 10000100 0101
(Word 1: 1 repeated 1 time,
Word 2: 0 repeated 4 times,
Word 3: 101)
Since this framework divides bitmaps into groups of w - 1 bits rather than w bits, this causes some misalignment, which needs to be resolved at query time.
Enhanced WAH (EWAH) solves this issue by introducing a new type of word in place of filler words, a marker word. The upper half of a marker word is used to encode the fill value and the run length. The rest of the marker word is used to represent the number of literal words following the run length. Since this implicitly determines the order of the words, no bits are required to explicitly encode the type of the word. Using the same example from WAH:
10000101
could be encoded as
10010000 01000001 101
However, if the data is not sparse enough and there are no patterns, it is usually better to store it verbatim. In order to perform operations between compressed and uncompressed without built-in support, one would have to decompress the compressed data, which causes a lot of overhead, the runtime cost of decompression.
The problem with current implementations of BSI is that there isn’t support for operations between compressed and uncompressed. To confront this issue, Professor Guzun has proposed hybrid BSI. Using the EWAH compression scheme, hybrid BSI would implement operations to return results without decompression. To do so, when reading in two different BSIs, there are 3 scenarios that need to be handled at each word: compressed and compressed, compressed and uncompressed, uncompressed and uncompressed. Through testing and analysis, we have proved the value in saving overhead through this technique by comparing the runtime of our operations with that of PyTorch tensors, which are optimized for machine learning models to run on GPUs and other specialized hardware. Now it’s time to show another application by comparing its runtime to NumPy arrays, which are optimized for matrix operations on large amounts of data.
Look into the future
Here’s what the schedule for the future looks like:
Week 5: Investigate best way to quantize decimals and multi-dimensional BSI
Week 6: Implement quantization
Week 7: Test quantization
Week 8: Try using same numpy algorithms for BSI for bigger operations (determinant, inverse, etc.)
Week 9: If succeeded last week, try for rest of algorithms, else setup conversion from BSI to numpy
Week 10: Write drafts of smaller operations in BSI (drop col, etc.)
Week 11: Runthrough with OpenSFM and collect some data
Week 12: Make presentation
Sources:
https://mac.r-project.org/openmp/
Guzun, G., Canahuate, G. Hybrid query optimization for hard-to-compress bit-vectors. The VLDB Journal 25, 339–354 (2016). https://doi.org/10.1007/s00778-015-0419-9
Week 5: Going to Win-“dows” Victories
Day N of building OpenSFM, I finally got it working on my WSL. Not Windows PowerShell, Windows Subsystem for Linux. In other words, Linux was the easiest to build OpenSFM on.
Building OpenSFM on Linux
To fix my installation error from last week,
Windows button > search Turn Windows features on or off > enable Virtual Machine Platform > restart computer
Now, I clone the repository, then run the commands from the Dockerfile for Ubuntu by adding “sudo” in front.
sudo apt-get update
sudo apt-get install -y build-essential cmake git libeigen3-dev libopencv-dev libceres-dev python3-dev python3-numpy python3-opencv python3-pip python3-pyproj python3-scipy python3-yaml curl
sudo apt-get clean
sudo rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
sudo pip3 install -r requirements.txt
sudo python3 setup.py build
To find the files stored in WSL, go to File Explorer, scroll down to Linux tab > Ubuntu > home > yourusername > OpenSfM. You can open this in your choice of code editor (VSCode, in my case).
Note that when you edit then save a file, you might get an error that you don’t have permission to do so.
To get access, run this command for your file:
sudo chown czhang /home/czhang/OpenSfM/filename.txt
This is the output from when OpenSFM finishes building. Never seen such beauty in my life.
To test that it worked, I ran it on the sample Berlin dataset.
I linked this project to my original Github fork by running
sudo git remote set-url origin https://username:personalaccesstoken@github.com/Lypho2012/OpenSfM.git
Setting the remote URL to contain your personal access token lets you commit changes without having to enter your personal access token every single time.
Before I get into my work process, there’s a few more tips to note about binary numbers.
Binary - Special Cases
Negative Numbers:
The two conventions for writing negative numbers in binary are sign magnitude and twos complement.
Sign magnitude uses the most significant bit to represent the sign, where 1 represents negative and 0 represents positive. For example:
Decimal Binary
7 0111
-7 1111
Twos complement uses the most significant bit to represent subtracting 2^n, where n is the length of the bit string. You can also convert from the positive to negative form by flipping all the bits and adding 1. For example:
Decimal Binary
7 0111
-7 1001
Floats:
The IEEE-754 is a convention for converting floating point numbers to binary. You can use this WikiHow to learn how it works.
However, there’s an even easier way to work with floats with quantization: multiplying the float by a power of 10 to preserve some digits after the decimal point. This is the method that I’m going to use when performing operations later on.
VSCode Python Debugger
I tried running a test file on its own, but I got an error on the “import opensfm” line that opensfm couldn’t be found. Since the build works and I could even reconstruct the Berlin sample, I suspect the problem might be with permissions since I have to run commands with “sudo” in the Terminal. For now, I decided to try setting up a debugger and running the commands from the Terminal so that I can attach to the process.
Run and Debug on the sidebar (ctrl+shift+d) > dropdown menu at the top > Add Configurations > Python Debugger > Attach Using Process ID
This generates a launch.json file.
Now, put a random input() statement (so that the code process starts but waits for you to attach your debugger) in your code and run it. If you select Python Debugger: Attach Using Process ID from the dropdown menu and click the run button, it should prompt you to pick a process.
For me, this didn’t show up, which I think is a problem with my VSCode, so I will have to either downgrade VSCode or input the process ID manually. There’s a way to see the list of processes by ctrl+shift+p and search Process Explorer, but even when I input the process ID, none of them would connect, so I have to test further.
Another way to use the debugger is to add a Python Debugger: Current File with Arguments configuration. I was trying to test the build statement python3 setup.py build, so I had to add the additional “sudo”: true, “redirectOutput”: true and set “console” to “integratedTerminal.” While this did run, I got an error that the “sphinx” module couldn’t be found, which I think is because this process likely runs the command “python setup.py build” and not “python3 setup.py build.”
What to Expect
For the upcoming week, I’m going to be trying to get the debugger to work so that I can start testing substituting in BSI in the code. I will likely need to create an auxiliary class using BSI as the base because BSI is not suited to work with multi-dimensional matrices. Instead, I will be using the algorithms from NumPy and calling BSI operations instead of the regular vector operations. Of course, since I have figured out the quantization part, this merits a change to the predicted upcoming schedule:
Week 6: Create auxiliary class for BSI arrays and start substituting it in small sections
Week 7: Continue testing and substituting further
Week 8: Try using same numpy algorithms for BSI for bigger operations (determinant, inverse, etc.)
Week 9: If succeeded last week, try for rest of algorithms, else setup conversion from BSI to numpy
Week 10: Write drafts of smaller operations in BSI (drop col, etc.)
Week 11: Runthrough with OpenSFM and collect some data
Week 12: Make presentation
Sources:
https://www.wikihow.com/Convert-a-Number-from-Decimal-to-IEEE-754-Floating-Point-Representation
Week 6: The “Four’s” of Python/C++
Very technical blog today because this week contained a special occasion. On the 4th day of the week, it was the 4th day of the 4th month of 2024. Whatever you were doing at 4:44 is now your lucky activity (or extremely unlucky, possibly death-inducing, if you believe in the unlucky number 4). To celebrate this miraculous string of 4’s, I will explain 4 methods today: 2 ways to use C++ in Python and vice versa and 2 ways to set up a debugger in the Terminal.
Debugging in the Terminal
While VSCode provides extensions to set up debuggers, ultimately all the work happens in the Terminal; VSCode extensions simply provide ease of access. However, in my case, those extensions are giving me issues because of bugs in the version of VSCode I have that hasn't been resolved on Linux yet, so I'm going to try relying only on the Terminal instead.
Method 1: GDB (GNU Debugger)
GDB is a command-line tool for debugging C, C++, Fortran, etc. To use it, compile a C++ executable by running this command to include debugging information
g++ -g -o filename filename.cpp
Then, activate GDB with the command
gdb filename
And run the program
run
If you want to attach it to a running process, you can also activate GDB with this command instead
gdb -p your_process_id
To use GDB with Python files, you need to set up the python-gdb extension. Then, build your Python file with debug mode by adding the option
--with-pydebug
After launching GDB, attach it to your running Python process that uses your extension module. The process is explained in full detail here.
Method 2: Valgrind
Valgrind is a toolset for catching memory leaks. You can set it up on your computer, then run tools like Memcheck by using the command
valgrind --leak-check=yes filename
To use it in tandem with GDB, you can get vgdb, an embedded GDB server. Launch Valgrind with the option
--vgdb=yes --vgdb-error=0
This turns on the embedded GDB server and tells Valgrind to wait for 0 errors to occur before waiting for a connection from GDB.
Then, in another window, launch GDB with your program and run
target remote | vgdb
Combining Python and C++
Method 1: Extension modules
The Python/C API is a well-defined process for either extending Python with C modules (extension modules) or using Python as an embedded component in a larger application. To define an extension module, you need to include this at the top of your C file
#define PY_SSIZE_T_CLEAN
#include <Python.h>
Then, you need an init function for Python to recognize your C file as an extension module.
PyMODINIT_FUNC
PyInit_your_module_name (void) {
// some code
}
In C, Python objects are represented by a single datatype, PyObject. When calling a function from the extension module in Python, the parameters need to be pointers to a PyObject. In order to read the PyObject, there are built-in functions like
PyObject_Length(PyObject *your_arr_name)
PyList_GetItem(PyObject *your_arr_name)
PyLong_AsLong(PyObject *your_int_or_long_name)
You can read more about the API on the official documentation. NumPy even has its own guide to using extension modules.
Method 2: pybind11
pybind11 is a simplified version of Boost.Python based on C++11 features, a library used to expose C++ code to Python and vice versa. The goal of Boost.Python is to minimize boilerplate (near repetitive) code in extension modules. Since pybind11 is a simplified version of Boost.Python, it’s an even more simplified version of extension modules. To use pybind11, you need to include this header file (you might need more header files based on the complexity of your program).
#include <pybind11/pybind11.h>
To make things easy, you can shorten the namespace name by writing
namespace py = pybind11;
Like with extension modules, you need an initialization function for Python to recognize your code.
PYBIND11_MODULE(your_module_name, m) {
// some code
}
To add a function, write it like you normally would above this function, then inside the PYBIND11_MODULE function, add this line
m.def(“your_function_name_in_python”, &your_function_name_in_c++);
To add a class, write the class, then in PYBIND11_MODULE, add this line
py::class_<your_class_name>(m, “your_class_name_in_python”)
.def(py::init())
If you want to add functions to the class, just add another .def().
You can read more about pybind11 in its official documentation.
There is also a way to directly cast C++ data into NumPy arrays using a buffer object.
Sources:
https://docs.python.org/3/howto/gdb_helpers.html#gdb
https://valgrind.org/docs/manual/quick-start.html
https://docs.python.org/3/c-api/intro.html
https://numpy.org/doc/stable/user/c-info.how-to-extend.html
https://pybind11.readthedocs.io/en/stable/basics.html
https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html
Week 7: Mesh of Languages
Last week, I introduced several C++ debuggers and extension module techniques. This week, I didn’t end up using the C++ debuggers because while I was hoping to trace the calls to the C++ code from Python, my gdb had problems connecting to the process, so I traced the Python code using the Python debugger (pdb) instead.
PDB
To use PDB, you can either import pdb at the beginning of the file and put pdb.set_trace() wherever you want a break in your code, or use breakpoint() to skip the import statement. When you run your Python program, it will automatically open the PDB command line.
Some commonly used commands are:
c to continue to next breakpoint
s to step into the current line
whatis some_var to get the data type of some_var
p some_var to print the contents of some_var
until linenumber to continue until that line number; continues to end of current function if the line number is greater than the last line number of the function
To get the process ID of your running Python script, use the os.getpid() function:
Connect GDB to the running process using gdb -p your_pid:
I had a problem with how something was linked as you can see from the colorful message on the bottom, so this didn’t connect GDB, but normally that process would be correct.
The other commands in opensfm_run_all
In the blogs of previous weeks, I’ve only been talking about the reconstruction part of the process, but there’s in fact code to process before and after as well. When you run opensfm_run_all, it executes a bash file with these commands
While reconstruct is one of the lengthiest processes, I needed to pick a command with less action so that I could start out with a test of a complete process with few changes. All of the commands invoke a call to opensfm_main.py, which initializes a DataSet to store the input images, then calls the respective command with the DataSet object. Since the DataSet object stores its data as matrices defined by the C++ Eigen class, which has built-in implementation for conversion to NumPy, I may have to create a list of types in a header file to direct the initialization of these matrices to vectors of BsiAttribute’s.
The current types.h file:
extract_metadata formats all the metadata in each image into attributes of the DataSet object.
detect_features detects features in images using OpenCV algorithms like SIFT (Scale-Invariant Feature Transform), SURF (Sped Up Robust Features), AKAZE (Accelerated KAZE), ORB (Oriented FAST and Rotated BRIEF), and HAHOG (Histogram of Oriented Gradients).
match_features generates similar pairs for matching using strategies based on GPS distance, time, triangulation, Bag of Words, and VLAD.
create_tracks creates a track manager that links matches into tracks using union find.
reconstruct calls the incremental reconstruction process.
mesh adds the Delaunay meshes to the reconstruction using the corresponding perspective transformations.
undistort exports the reconstruction to NVM_V3 format from VisualSfM.
compute_depthmaps finds the angle between neighboring images to merge depthmaps.
Out of these, I decided to test mesh because it has one of the higher level linear algebra operations, inverse, and a lower level dot product, ignoring the operations in the cameras that aren’t perspective, which is what my phone takes pictures in. Since dot product was already implemented in BsiAttribute, I started out with a small testcase using pybind11 to make sure I knew how to use it.
pybind11 example
After writing the function and adding it to pybind11 in your C++ code,
You need to add these lines to created a shared library using your pybind11 module
Here, you can see that pybsi was generated along with the other .so files when running setup.py.
Now you can call it from your Python code.
Now that I know this pybind11 structure works, I need to implement inverse with BsiAttribute. To avoid having to run setup.py every time I want to test whether it works in Python, I’ll start off by testing it separately in C++. I turned to using CLion because I was already familiar with the debugger.
Inverse
NumPy uses the sgesv LAPACK routine for inverse. LAPACK (Linear Algebra Package) is a collection of high level linear algebra operations written in Fortran. It looks a lot like pseudo code because it’s mostly used for mathematical computation.
sgesv is a routine that solves for X, the solution to A * X = B, where A is a n-by-n matrix and X and B are both n-by-m matrices, using LU decomposition with partial pivoting (only row interchanges). This is the equivalent of using Gaussian elimination to find the inverse.
LU decomposition breaks down a matrix A into
A = P * L * U.
In the case of inverse,
P is the permutation matrix that stores the row interchanges,
L is the lower (all 0’s above the diagonal) unit (has 1’s along the diagonal) triangular matrix that stores the elementary elimination to get to row-echelon form, and
U is the upper triangular matrix that stores the elementary elimination to get to reduced row-echelon form.
By solving for the right hand side using an identity matrix, you get the inverse.
So far, I’ve set up a skeleton for the main functions and written and tested a compare function that compares the elements at the same index of a BsiAttribute.
The testing code:
Next week, I’ll probably be continuing to write the inverse function since there are several small functions in between that haven’t been implemented yet, and I need to run testcases on all of them to make sure they’re accurate for development.
Sources:
https://docs.python.org/3/library/pdb.html
https://www.netlib.org/lapack/
Week 8: Breakpoint and Retracing Steps
Before I got into the testing phase of inverse, I was made aware of a huge miscalculation: the application of the matrix inverse I had found in the mesh command only took inputs of 3x3 matrices. When I searched the repository for the other applications, I came up with the same result. In fact, even though there were templates of Eigen matrices of n by 1 and n by n defined in the types.h file, the variable-dimensioned matrices were only used to verify whether an object was a matrix of that type.
This was a pretty big setback for me because that means there isn’t really any room for optimization with inputs that small. While it is disappointing to not be able to apply BSI to such an orthodox linear algebra operation as of yet, it will definitely be something I want to continue in the future because of how perfect BSI is for linear algebra operations consisting only of row operations (partial pivoting). For now, it’s back to the drawing board for me.
compute_depthmaps
To find what operations do take the longest, I put logging statements and timer functions throughout the code in addition to the ones already in place for reports. Here’s an example of some of the outputs from compute_depthmaps, which generates a dense point cloud from depth maps which store information about the location of points relative to each camera position, in microseconds. A lot of the function names in the output are repeated because of parallel computing.
From here, I found my first candidates for implementation, PatchMatchForwardPass and PatchMatchBackwardPass, which apply a function to update all the pixels in a depthmap of dimensions 852 by 640 with patch size 8.
I will have to read the function more to see whether it’s possible to run the algorithm on the whole array of points at once rather than one by one.
reconstruct
The operation that took the next longest was the global bundle, which bundle adjusts all of the cameras/images in the current reconstruction to fit the information from the new image.
Since the operation that takes the longest is the Ceres solve function, I also have to keep in mind that I might have to either add the Ceres library as a subdirectory to the project instead of just including it from my local system so that I can edit it, or write a different object to replace the Ceres Solver.
Other options
My remaining candidates are the NumPy operations in reconstruct that do operate on large vectors/matrices. While they don’t take a large percentage of the total time to run OpenSFM, it could be valuable as a stand-alone testcase.
Many NumPy operations can be done using row operations alone because many of them rely on solving systems of equations. While the LAPACK routine algorithm isn’t entirely suited to BSI because it comprises of breaking the matrix into submatrices and operation on smaller chunks, which BSI is not optimized to do, BSI can do elementary row eliminations in linear time and reduced space complexity. This makes it trivial to convert a matrix to reduced row-echelon form, and performing the same row operations on an identity matrix will result in the inverse. This was the method I opted for when implementing inverse using BSI.
For the upcoming week, I will likely be reading the library to understand my options and hopefully devise a plan to implement them because they are not the conventional linear algebra or vector operations.
Week 9: Automatic differentiation but not automatic comprehension
When I saw that the function taking the longest was from Ceres Solver, I was very surprised to see the familiar name and quite happy that it turned out to be so important after all the hours I was forced to labor for it. The operation that I looked into this week was bundle adjust. To start off, let’s revisit how reconstruction works. First, a pair of images is added to the reconstruction, kicking off the grow_reconstruction function. Using the matching features, the RANSAC random sampling algorithm is used to create a triangulation, a representation of the geometric model using triangles. To increase the accuracy of the triangulation, it is bundle adjusted to the prior reconstruction by minimizing the discrepancy in their projections. If there’s GPS or GCP (ground control points) data as well, those are considered first to be absolutely accurate.
Ceres Solver
Ceres Solver, the library I was struggling with several weeks ago to setup, is a framework for solving optimization problems. To use Ceres Solver, you need to setup the problem with a cost function and a parameter block. A cost function is the customizable measure for difference between expected and actual values. A parameter block represents either a single parameter or a group of scalars that define a parameter, like the 3 components of a translation vector and the 4 components of a quaternion for camera poses. Typically the least squares method is used on the cost function to calculate the residual block. Then, the parameter block is iteratively tweaked until the residual block is minimized. Since we are working with experimental data, which might have outliers, we can introduce loss functions to reduce the impact of the outliers on the residuals.
To calculate the change in the residual block after tweaking the parameter block, we need to know the Jacobian. The Jacobian is a matrix representation of the partial derivatives of a function with respect to each of its parameters, for each of the functions defined with respect to those parameters. Rather than analytically calculating and storing the Jacobian, there is a method called automatic differentiation.
Automatic differentiation is based on the concept of dual numbers. Similar to imaginary numbers, dual numbers introduce another unit e which has the property that e^2 = 0 but e is an infinitesimally small number that approaches but does not equal 0. We can make use of this unit to calculate the derivative of a function by substituting the variable which we are taking the derivative with respect to, with the sum of that variable and e. By solving for the function, the derivative is the coefficient of e. For example, if we had the function
f(x) = x(x + 1)
We know that the derivative of f(x) must be
f’(x) = 2x + 1
If we substitute x + e for x,
f(x + e) = (x + e)(x + e + 1)
= x^2 + (2x + 1)e + x + e^2
In Ceres Solver, these dual numbers are implemented with the class Jet.
Bundle Adjust
The primary function of bundle adjust is to setup the minimization problem. The linear solver used for bundle adjust is usually SPARSE_SCHUR because the solution tends to be sparse due to the lack of correlation between the parameters for different points and cameras. Here is a sample of how much data was in each of the parameter blocks:
There are cost functions for prior errors, relative motion/rotation errors, common position errors, up vector errors, pan/tilt/roll angle errors, and linear motion errors. After adding all the cost functions, parameter blocks, and loss functions, the ceres::Solve() function is called to initialize the process. Each complete run took about 10 seconds on my 60 images.
I didn’t get to do that much this week because I was bedridden with the flu for most of it (never been sick so often in the period of a few months.) For the upcoming week, I’ll be looking at the function in compute_depthmaps, PatchMatchUpdatePixel().
Week 10: Levenberg Marquardt
Before I get into compute_depthmaps, I looked a little deeper into the math behind bundle adjustment and dug up some concepts that heavily rely on Linear Algebra and Multivariable Calculus.
Minimization Problem
Last week, I said that bundle adjust uses a linear solver but never specified what the linear problem was. Bundle adjustment is actually a nonlinear least squares problem that solves the optimization problem of minimizing
where F(x) is vector of m functions in n-dimensional space.
The solution to the vector x can be approximated by taking step 𝚫x each iteration. In other words, the minimization problem can be written as
Gradient Descent Method
The gradient descent method takes steps in the direction of steepest descent using step size determined by a positive scalar ⍺. ⍺ multiplied by the Jacobian matrix is the update to F(x). With the proper selection of ⍺, this method always converges, but may take a long time doing so.
Gauss-Newton Method
For moderately-sized problems, the Gauss-Newton method converges much faster than the gradient descent method, but there’s a chance that it diverges.
Using the Taylor series expansion,
If is small, the higher order terms can be ignored, yielding this linear function
with an error on the order of (x - 𝚫x)2.
is the Jacobian matrix J(x). This final expression is called the Jacobian linearization
Solve for 𝚫x for a number of iterations.
Levenberg-Marquardt Method
To make sure F(x + 𝚫x) converges, the Levenberg-Marquardt method introduces a regularization term to dampen the Hessian matrix (matrix of 2nd degree partial derivatives).
𝜆 is a regularization parameter that is larger at the beginning to make 𝚫x smaller. As the iterations progress and the solution improves, 𝜆 is decreased, and the solution accelerates toward convergence as the Levenberg-Marquardt method becomes more like the Gauss-Newton method.
compute_depthmaps
Now that I’ve finished introducing the minimization problem in bundle adjustment, it’s time to look at the algorithm behind compute_depthmaps. The entire command looks at the data for depth of points in the images and compiles a ply file of the points. To check which image has the most accurate location for each point, PatchMatchUpdatePixel() uses Normalized Cross Correlation (NCC). The equation for NCC is the cosine of the angle between the images.
For each of the pixels, the function checks the neighboring pixels and 6 random planes for the same pixel and takes the position with the highest cross correlation score. Since the function checks several random planes, the position might be slightly different every time compute_depthmaps is run.
Source:
Week 11: Depths of Pixels and Math
After expanding the equations in PatchMatchUpdatePixel(), I have arrived at a formula that I can use to express the position of the pixel being solved for in terms of BSI.
In bundle adjust, the camera poses of all of the images were solved for and stored as rotation and translation matrices, which are extrinsic parameters of the camera. Using those matrices, compute_depthmaps will now calculate the 3D positions of the pixels in the images by randomizing the depth and calculating the other coordinates using some formulas based on the camera extrinsic parameters.
Equations
Since all of the camera extrinsic parameters were calculated on a coordinate system centering around the first image, the position of the pixels are calculated with respect to the first image as well. As such, the homography matrix, a 3x3 matrix representing the transformation matrix from one image to another, needs to be calculated between each image and the first image.
H = 3x3 homography matrix
K = 3x3 projection matrix, how cameras projects onto an image plane based on focal length
R = 3x1 rotation matrices of the first image and the other image combined
t = 3x1 translation matrix
Using this homography matrix, the 3D points of each pixel in the image can be calculated
Similarly, based on the homography matrix, the a pixel in the first image can be projected to a position in the other image.
Since these equations rely purely on H, i, and j, we can store the values of H, i, and j for each iteration and rewrite the equations in terms of a vector of each element. For example,
where multiplication between vectors is defined as
Now that the equations have been rewritten in terms of vectors, it can be easily rewritten in terms of BSI by converting those vectors into BSI.
Implementation
For now, I will only implement the equations for u, v, and w in terms of H because I need to start off somewhere simple. To solve for H in terms of the extrinsic parameters would require 36 additional vectors for each element in each of the matrices. I also can’t implement the Hx0 or dfdx variables yet because I haven’t implemented the framework for division with BSI.
Currently, I have created a vector for all of the H matrices, i, and j values for each iteration of PatchMatchForwardPass and PatchMatchBackwardPass. There are about 4 million H matrices per iteration, and each of the elements of H are floats ranging from 0 to 2 with up to 5 decimal places. Since BSI uses quantization to store floats, I will test with different precisions of H.
For now, my program keeps crashing after some number of iterations, so I need to figure out a way to break up how many values of H I store into even more chunks.
BSI Division
As I mentioned, I haven’t yet implemented the framework for division, but I can explain the algorithm that I plan to use. Division by a number can be represented as multiplying by the inverse of that number.
Fractions can be rewritten such that it is in the form
Since X is stored in binary format, it is easy to divide X by powers of 2 by shifting X. As such, it would be easy to divide X by D if we could instead divide X by a power of 2, then multiply by some constant Y.
In other words, the division can be rewritten as
Depending on how much precision the final answer needs to be in, there may be some multiplication by a positive power of 2 with the BSI performed earlier on, which may be later offset by storing the power of 2 the BSI was multiplied by.
For the upcoming week, I will be fixing the cause of my program crashing, likely due to memory limit exceeded. Once I have fixed that, I will be able to show some numbers for runtime and accuracy on the final presentation!
Week 12: Concluding Notes on Final Presentation
I was able to debug the sum and multiply algorithms in BSI in time for my presentation, run some testcases, and even generate a point cloud of my rabbit. Here’s a link to the video of my presentation where you’ll be able to see my results:
It was a long journey, but it’s nowhere near over. Like I mentioned in my conclusion, I would like to finish my implementation so that the entire compute depthmaps operation relies on BSI so that I don’t have to do a conversion in the middle. My results have shown that BSI runs 50-100 times faster than the linear solution, which spells a huge improvement to the OpenSfM algorithm overall.
I want to end off on the note that photogrammetry has a huge variety of applications, and there are still many areas within it that can be explored, so I would encourage anyone who has been reading my blogs to consider diving into the field. I couldn’t have picked a better project to end off my senior year with.
This ends the last blog of my senior project. Thank you for everyone who has been following my blogs! I wish you the best of luck if you are also pursuing your own senior project.