Sunday, June 4, 2017

Mars, Middle School, and 3D Models

"Dad, can you print out a 3D model of NE Syrtis Major?"
"Where's that?"
"On Mars. NASA might land the next rover there. Can't you just find a model and print it? I need it for science tomorrow. We're doing group projects."

A quick check on Thingiverse revealed no such model - not surprising. This presents quite the dilemma: Let the boy learn a lesson not to procrastinate or work on a really cool 3D printing project? The project was too cool to pass up but I have a tiny bit of guilt enabling the boy's poor planning.

Two problems to tackle first. Where is Syrtis Major and where do I get elevation data? The latter turns out to be really easy. I remember when NASA mapped Mars with a laser altimeter. A little googling and I find a digital elevation model downloadable as a TIFF file. Next problem is locating Syrtis Major. Some more googling and I download a nice paper with an image and latitude/longitude values. The candidate landing site in NE Srytis is topographically boring (i.e., flat) for obvious reasons. But the whole Srytis Major province is perched 3 miles above Isidis Planitia and a 2000-mile diagonal map is dramatic. That's what we'll target. 

I wrote a Matlab script (see below) to read in the TIFF file, subset out the area, smooth and decimate the surface down to something manageable. Someone wrote a nice mesh to STL exporter. My first crack at it would take 32 hours to print and since this is a fire fighting exercise, I smoothed it down even more and made the model only 6" on a side to get a 14 hour predicted print time. Right about when the school bus comes. The horizontal scale is 16,000,000:1 and the vertical is 220,000:1, which is about a 70x exaggeration in the elevation.

In my experience, many 3D models have errors so they are not manifold. Windows 10 has a built in app called 3D Builder that is great at correcting errors. I usually run my STL files through it before trying to slice them for printing. Here's a rendering of the model in 3D Builder:
The STL file is up on Thingiverse. I sliced it with the 0.2-mm normal setting in Prusa 3D Sli3er, loaded to my Octoprint server and started the job. Here's the outcome. Not bad!


Matlab script to generate the source STL file:

%%
a=imread('Mars_MGS_MOLA_DEM_mosaic_global_463m.tiff');;
%% maps
% https://astrogeology.usgs.gov/search/details/Mars/GlobalSurveyor/MOLA/Mars_MGS_MOLA_DEM_mosaic_global_463m/cub
% Minimum Latitude -90
% Maximum Latitude 90
% Minimum Longitude -180
% Maximum Longitude 180
% Direct Spatial Reference Method Raster
% Object Type Pixel
% Lines (pixels) 23040
% Samples (pixels) 46080
% Bit Type 16
% Radius A 3396000
% Radius C 3396000
% Bands 1
% Pixel Resolution (meters/pixel) 463.0836
% Scale (pixels/degree) 128
% Horizontal Coordinate System Units Meters
% Map Projection Name Simple Cylindrical
% Latitude Type Planetocentric
% Longitude Direction Positive East
% Longitude Domain -180 to 180
imshow(a(1:100:end,1:100:end))
%% subset
% http://onlinelibrary.wiley.com/doi/10.1029/2003JE002143/abstract
% http://www.planetary.brown.edu/pdfs/2763.pdf
% Figure 1: The map covers an area from -10S to 30N and from 270W to 315W.
% 270 W is -90 E and 315 W is -45 E
lat=linspace(-90,90,size(a,1));
lon=linspace(-180,180,size(a,2));
x=interp1(lat,1:size(a,1),-[30 -10],'nearest');
y=interp1(lon,1:size(a,2),-[-45 -90],'nearest');
b=a(x(1):x(2),y(1):y(2));
%%
figure
N=5;
pcolor(flipud(b(1:N:end,1:N:end)))
shading flat
%% smooth
S=25;
c = double(imgaussfilt(b,S));
%% plot
figure
N=2*S;
mesh((c(1:N:end,1:N:end)))
shading flat
%% make surface to export
N=2*S;
d=c(1:N:end,1:N:end);
d=d-min(d(:));
d(d<0)=0;
d(1:end,1)=0;
d(1:end,end)=0;
d(1,1:end)=0;
d(end,1:end)=0;
d=d*35/max(d(:));
[u,v]=meshgrid(1:size(d,2),1:size(d,1));
u=u*150/max(u(:));
v=v*150/max(v(:));
%%
mesh(u,v,d); axis equal
%%  export
% https://www.mathworks.com/matlabcentral/fileexchange/20922-stlwrite-filename--varargin-
stlwrite('mars4.stl',u,v,d)