Code
-using FileIO
-using MeshIO
-
-using stlProcess
-
-using CSV
-using DataFrames
-
-using LinearAlgebra
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b3133b9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +.jupyter_cache +*.log +*.tex \ No newline at end of file diff --git a/.jupyter_cache/__version__.txt b/.jupyter_cache/__version__.txt deleted file mode 100644 index 79a2734..0000000 --- a/.jupyter_cache/__version__.txt +++ /dev/null @@ -1 +0,0 @@ -0.5.0 \ No newline at end of file diff --git a/.jupyter_cache/executed/2f3f010a51c40fce4d85e1fa29ada8c0/base.ipynb b/.jupyter_cache/executed/2f3f010a51c40fce4d85e1fa29ada8c0/base.ipynb deleted file mode 100644 index 5139368..0000000 --- a/.jupyter_cache/executed/2f3f010a51c40fce4d85e1fa29ada8c0/base.ipynb +++ /dev/null @@ -1,248 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "a45c1344", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n", - "(process:28652): GLib-GIO-WARNING **: 12:50:31.472: Unexpectedly, UWP app `KDEe.V.Okular_21.1203.941.0_x64__7vt06qxq7ptv8' (AUMId `KDEe.V.Okular_7vt06qxq7ptv8!KDEe.V.Okular') supports 5 extensions but has no verbs\n" - ] - } - ], - "source": [ - "import IJulia\n", - "\n", - "# The julia kernel has built in support for Revise.jl, so this is the \n", - "# recommended approach for long-running sessions:\n", - "# https://github.com/JuliaLang/IJulia.jl/blob/9b10fa9b879574bbf720f5285029e07758e50a5e/src/kernel.jl#L46-L51\n", - "\n", - "# Users should enable revise within .julia/config/startup_ijulia.jl:\n", - "# https://timholy.github.io/Revise.jl/stable/config/#Using-Revise-automatically-within-Jupyter/IJulia-1\n", - "\n", - "# clear console history\n", - "IJulia.clear_history()\n", - "\n", - "# Intialize Plots w/ default fig width/height\n", - "try\n", - " fig_width = 7\n", - " fig_height = 5\n", - " fig_format = :retina\n", - " fig_dpi = 96\n", - " # no retina format type, use svg for high quality type/marks\n", - " if fig_format == :retina\n", - " fig_format = :svg\n", - " # IJulia doesn't support PDF output so use png (if the DPI \n", - " # remains the default of 300 then set to 96)\n", - " elseif fig_format == :pdf\n", - " fig_format = :png\n", - " fig_dpi = 96\n", - " end\n", - " # convert inches to pixels\n", - " fig_width = fig_width * fig_dpi\n", - " fig_height = fig_height * fig_dpi\n", - " using Plots\n", - " gr(size=(fig_width, fig_height), fmt = fig_format, dpi = fig_dpi)\n", - "catch e\n", - " # @warn \"Plots init\" exception=(e, catch_backtrace())\n", - "end\n", - "\n", - "# Set run_path if specified\n", - "try\n", - " run_path = \"\"\n", - " if !isempty(run_path)\n", - " cd(run_path)\n", - " end\n", - "catch e\n", - " @warn \"Run path init:\" exception=(e, catch_backtrace())\n", - "end\n", - "\n", - "# don't return kernel dependencies (b/c Revise should take care of dependencies)\n", - "nothing\n" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "189db2ad", - "metadata": {}, - "outputs": [], - "source": [ - "#| code-fold: true\n", - "#| output: false\n", - "using FileIO\n", - "using MeshIO\n", - "\n", - "using stlProcess\n", - "\n", - "using CSV\n", - "using DataFrames\n", - "\n", - "using LinearAlgebra" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "5e47d862", - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
0 rows × 8 columns (omitted printing of 1 columns)
surface_area | characteristic_length | sbx | sby | sbz | Ix | Iy | |
---|---|---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 |
8 rows × 7 columns
variable | mean | min | median | max | nmissing | eltype | |
---|---|---|---|---|---|---|---|
Symbol | Float64 | Float64 | Float64 | Float64 | Int64 | DataType | |
1 | surface_area | 25.2002 | 5.60865 | 13.3338 | 159.406 | 0 | Float64 |
2 | characteristic_length | 79.5481 | 0.158521 | 1.55816 | 1582.23 | 0 | Float64 |
3 | sbx | 1.40222 | 0.0417367 | 0.967078 | 10.0663 | 0 | Float64 |
4 | sby | 3.3367 | 0.0125824 | 2.68461 | 9.68361 | 0 | Float64 |
5 | sbz | 3.91184 | 0.29006 | 1.8185 | 14.7434 | 0 | Float64 |
6 | Ix | 1.58725 | 0.0311782 | 0.23401 | 11.1335 | 0 | Float64 |
7 | Iy | 3.74345 | 0.178598 | 1.01592 | 24.6735 | 0 | Float64 |
8 | Iz | 5.20207 | 0.178686 | 1.742 | 32.0083 | 0 | Float64 |
0 rows × 8 columns (omitted printing of 1 columns)
surface_area | characteristic_length | sbx | sby | sbz | Ix | Iy | |
---|---|---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 |
8 rows × 7 columns
variable | mean | min | median | max | nmissing | eltype | |
---|---|---|---|---|---|---|---|
Symbol | Float64 | Float64 | Float64 | Float64 | Int64 | DataType | |
1 | surface_area | 25.2002 | 5.60865 | 13.3338 | 159.406 | 0 | Float64 |
2 | characteristic_length | 79.5481 | 0.158521 | 1.55816 | 1582.23 | 0 | Float64 |
3 | sbx | 1.40222 | 0.0417367 | 0.967078 | 10.0663 | 0 | Float64 |
4 | sby | 3.3367 | 0.0125824 | 2.68461 | 9.68361 | 0 | Float64 |
5 | sbz | 3.91184 | 0.29006 | 1.8185 | 14.7434 | 0 | Float64 |
6 | Ix | 1.58725 | 0.0311782 | 0.23401 | 11.1335 | 0 | Float64 |
7 | Iy | 3.74345 | 0.178598 | 1.01592 | 24.6735 | 0 | Float64 |
8 | Iz | 5.20207 | 0.178686 | 1.742 | 32.0083 | 0 | Float64 |
using FileIO
-using MeshIO
-
-using stlProcess
-
-using CSV
-using DataFrames
-
-using LinearAlgebra
Orbital debris is a form of pollution that is growing at an exponential pace and puts current and future space infrastructure at risk. Satellites are critical to military, commercial, and civil operations. Unfortunately, the space that debris occupies is increasingly becoming more crowded and dangerous, potentially leading to a cascade event that could turn orbit around the Earth into an unusable wasteland for decades unless proper mitigation is not introduced. Existing models employed by NASA rely on a dataset created from 2D images and are missing many crucial features required for correctly modeling the space debris environment. This approach aims to use high-resolution 3D scanning to fully capture the geometry of a piece of debris and allow a more advanced analysis of each piece. Coupled with machine learning methods, the scans will allow advances to the current cutting edge. Physical and photograph-based measurements are time-consuming, hard to replicate, and lack precision. 3D scanning allows much more advanced and accurate analysis of each debris sample, focusing on properties such as moment of inertia, cross-section, and drag. Once the characteristics of space debris are more thoroughly understood, we can begin mitigating the creation and danger of future space debris by implementing improved satellite construction methods and more advanced debris avoidance measures.
+This project aims to fix very difficult issues, and although great progress has been made there is still plenty of work to be done. Currently, algorithms have been made that are capable of getting many key features from solid 1 models that are in the stl
format. The algorithm for processing the 3D meshes is implemented in the Julia programming language. Syntactically the language is very similar to Python and Matlab. Julia was chosen because it is nearly as performant as compiled languages like C, while still having tooling geared towards engineers and scientists. The code produces a struct with all the calculated properties as follows:
1 A mesh with a surface that is fully closed and has no holes in its geometry.
struct Properties
+# Volume of the mesh
+ ::Float64
+ volume# Center of gravity, meshes are not always center at [0,0,0]
+ ::Vector{Float64}
+ center_of_gravity# Moment of inertia tensor
+ ::Matrix{Float64}
+ inertia# Surface area of mesh
+ ::Float64
+ surface_area# Average orthogonal dimension of the mesh
+ ::Float64
+ characteristic_length# Projected length of farthest two points in [x,y,z] directions
+ ::Vector{Float64}
+ solidbody_valuesend
# local path to https://gitlab.com/orbital-debris-research/fake-satellite-dataset
= raw"C:\Coding\fake-satellite-dataset"
dataset_path
-= ["1_5U", "assembly1", "cubesat"]
- folders
-= DataFrame(;
- df =Float64[],
- surface_area=Float64[],
- characteristic_length=Float64[],
- sbx=Float64[],
- sby=Float64[],
- sbz=Float64[],
- Ix=Float64[],
- Iy=Float64[],
- Iz )
for path in dataset_path * "\\" .* folders
-println("Processing Path: ", path)
- Threads.@threads for file in readdir(path)
- = load(path * "\\" * file)
- stl = find_scale(stl)
- scale = get_mass_properties(stl; scale=scale)
- props
-= eigvals(props.inertia)
- eigs = sortperm(eigs)
- sort_index = eigs[sort_index]
- Ix, Iy, Iz = props.sb_values[sort_index]
- sbx, sby, sbz
-push!(
-
- df,
- [props.surface_area, props.characteristic_length, sbx, sby, sbz, Ix, Iy, Iz],
- )end
- end
Processing Path: C:\Coding\fake-satellite-dataset\1_5U
-Processing Path: C:\Coding\fake-satellite-dataset\assembly1
-Processing Path: C:\Coding\fake-satellite-dataset\cubesat
-┌ Warning: Characteristic Length Algorithm failed to converge, this usually means stl is flat. Setting length in dir to 0.
-└ @ stlProcess C:\Users\albig\.julia\packages\stlProcess\8rsc7\src\stlProcess.jl:153
-describe(df)
8 rows × 7 columns
variable | mean | min | median | max | nmissing | eltype | |
---|---|---|---|---|---|---|---|
Symbol | Float64 | Float64 | Float64 | Float64 | Int64 | DataType | |
1 | surface_area | 25.2002 | 5.60865 | 13.3338 | 159.406 | 0 | Float64 |
2 | characteristic_length | 79.5481 | 0.158521 | 1.55816 | 1582.23 | 0 | Float64 |
3 | sbx | 1.40222 | 0.0417367 | 0.967078 | 10.0663 | 0 | Float64 |
4 | sby | 3.3367 | 0.0125824 | 2.68461 | 9.68361 | 0 | Float64 |
5 | sbz | 3.91184 | 0.29006 | 1.8185 | 14.7434 | 0 | Float64 |
6 | Ix | 1.58725 | 0.0311782 | 0.23401 | 11.1335 | 0 | Float64 |
7 | Iy | 3.74345 | 0.178598 | 1.01592 | 24.6735 | 0 | Float64 |
8 | Iz | 5.20207 | 0.178686 | 1.742 | 32.0083 | 0 | Float64 |
for path in dataset_path * "\\" .* folders
+Threads.@threads for file in readdir(path)
+ = load(path * "\\" * file)
+ stl = find_scale(stl)
+ scale = get_mass_properties(stl; scale=scale)
+ props
+= eigvals(props.inertia)
+ eigs = sortperm(eigs)
+ sort_index = eigs[sort_index]
+ I1, I2, I3 = props.sb_values[sort_index]
+ sbx, sby, sbz
+push!(
+
+ df,
+ [props.surface_area, props.characteristic_length, sbx, sby, sbz, I3, I2, I1],
+ )end
+ end
8 rows × 7 columns
variable | mean | min | median | max | nmissing | eltype | |
---|---|---|---|---|---|---|---|
Symbol | Float64 | Float64 | Float64 | Float64 | Int64 | DataType | |
1 | surface_area | 17.6375 | 6.62244 | 14.4423 | 52.2496 | 0 | Float64 |
2 | characteristic_length | 35.3451 | 0.2662 | 1.50102 | 434.102 | 0 | Float64 |
3 | sbx | 1.69789 | 0.198322 | 1.48948 | 3.60402 | 0 | Float64 |
4 | sby | 2.30204 | 0.135151 | 1.00749 | 8.0472 | 0 | Float64 |
5 | sbz | 2.25255 | 0.29006 | 1.48957 | 7.76479 | 0 | Float64 |
6 | Ix | 1.79204 | 0.178686 | 1.57634 | 6.83016 | 0 | Float64 |
7 | Iy | 1.50921 | 0.178598 | 0.844538 | 6.54812 | 0 | Float64 |
8 | Iz | 0.507267 | 0.0311782 | 0.279339 | 2.19511 | 0 | Float64 |
= cov(Matrix(df))
+ S
+= eigvals(S);
+ eig_vals
+# sorting eigenvalues from largest to smallest
+= sortperm(eig_vals; rev=true)
+ sort_index
+= eig_vals[sort_index]
+ lambda = names(df)[sort_index]
+ names_sorted
+= cumsum(lambda) ./ sum(lambda)
+ lambda_ratio
+plot(lambda_ratio, marker=:x)
+xticks!(sort_index,names(df), xrotation = 15)
To get started on the project before any scans of the actual debris are made available, I opted to find 3D models online and process them as if they were data collected by my team. GrabCAD is an excellent source of high-quality 3D models, and all the models have, at worst, a non-commercial license making them suitable for this study. The current dataset uses three separate satellite assemblies found on GrabCAD, below is an example of one of the satellites that was used.
@@ -2607,16 +2626,16 @@ Murray, James, Heather Cowardin, J-C Liou, Marlon Sorge, Norman Fitz-Coy, and ToThe algorithm’s speed is critical not only for the eventual large number of debris pieces that have to be processed, but many of the data science algorithms we plan on performing on the compiled data need the data to be normalized. For the current dataset and properties, it makes the most sense to normalize the dataset based on volume. Volume was chosen for multiple reasons, namely because it was easy to implement an efficient algorithm to calculate volume, and currently, volume produces the least amount of variation out of the current set of properties calculated. Unfortunately, scaling a model to a specific volume is an iterative process, but can be done very efficiently using derivative-free numerical root-finding algorithms. The current implementation can scale and process all the properties using only 30% more time than getting the properties without first scaling.
-
- Row │ variable mean min median max
- ─────┼───────────────────────────────────────────────────────────────────
- 1 │ surface_area 25.2002 5.60865 13.3338 159.406
- 2 │ characteristic_length 79.5481 0.158521 1.55816 1582.23
- 3 │ sbx 1.40222 0.0417367 0.967078 10.0663
- 4 │ sby 3.3367 0.0125824 2.68461 9.68361
- 5 │ sbz 3.91184 0.29006 1.8185 14.7434
- 6 │ Ix 1.58725 0.0311782 0.23401 11.1335
- 7 │ Iy 3.74345 0.178598 1.01592 24.6735 8 │ Iz 5.20207 0.178686 1.742 32.0083
+ Row │ variable mean min median max
+ ─────┼───────────────────────────────────────────────────────────────────
+ 1 │ surface_area 25.2002 5.60865 13.3338 159.406
+ 2 │ characteristic_length 79.5481 0.158521 1.55816 1582.23
+ 3 │ sbx 1.40222 0.0417367 0.967078 10.0663
+ 4 │ sby 3.3367 0.0125824 2.68461 9.68361
+ 5 │ sbz 3.91184 0.29006 1.8185 14.7434
+ 6 │ Ix 1.58725 0.0311782 0.23401 11.1335
+ 7 │ Iy 3.74345 0.178598 1.01592 24.6735 8 │ Iz 5.20207 0.178686 1.742 32.0083
Above is a summary of the current 108 part with scaling. Since all the volumes are the same it is left out of the dataset, the center of gravity is also left out of the dataset since it currently is just an artifact of the stl
file format. There are many ways to determine the ‘center’ of a 3D mesh, but since only one is being implemented at the moment comparisons to other properties doesn’t make sense. The other notable part of the data is the model is rotated so that the magnitudes of Iz
, Iy
, and Ix
are in descending order. This makes sure that the rotation of a model doesn’t matter for characterization. The dataset is available for download here:
The first step toward characterization is to perform a principal component analysis to determine what properties of the data capture the most variation. PCA
also requires that the data is scaled, so as discussed above the dataset that is scaled by volume
will be used. PCA
is implemented manually instead of the Matlab built-in function as shown below:
% covaraince matrix of data points
-S=cov(scaled_data);
-
-% eigenvalues of S
-eig_vals = eig(S);
-
-% sorting eigenvalues from largest to smallest
-lambda, sort_index] = sort(eig_vals,'descend');
- [
-
-lambda_ratio = cumsum(lambda) ./ sum(lambda)
Then plotting lambda_ratio
, which is the cumsum
/sum
produces the following plot:
% covaraince matrix of data points
+S=cov(scaled_data);
+
+% eigenvalues of S
+eig_vals = eig(S);
+
+% sorting eigenvalues from largest to smallest
+lambda, sort_index] = sort(eig_vals,'descend');
+ [
+
+lambda_ratio = cumsum(lambda) ./ sum(lambda)
Then plotting lambda_ratio
, which is the cumsum ./ sum
produces the following plot:
The current dataset can be described incredibly well just by looking at Iz
, which again the models are rotated so that Iz
is the largest moment of inertia. Then including Iy
and Iz
means that a 3D plot of the principle moments of inertia almost capture all the variation in the data.
The next step for characterization is to get only the inertia’s from the dataset. Since the current dataset is so small, the scaled dataset will be used for rest of the characterization process. Once more parts are added to the database it will make sense to start looking at the raw dataset. Now we can proceed to cluster the data using the k-means method of clustering. To properly use k-means a value of k, which is the number of clusters, needs to be determined. This can be done by creating an elbow plot using the following code:
-for ii=1:20
-idx,~,sumd] = kmeans(inertia,ii);
- [J(ii)=norm(sumd);
- end
for ii=1:20
+idx,~,sumd] = kmeans(inertia,ii);
+ [J(ii)=norm(sumd);
+ end
Which produces the following plot: