From 0604739e8d3a27a4642cc9b92aa4696d14195d80 Mon Sep 17 00:00:00 2001 From: Richard Chapman Date: Sat, 13 Dec 2014 16:56:27 +0000 Subject: [PATCH 1/4] Add vertical option Add compensation for longitude Fix error if width and height not equal Signed-off-by: Richard Chapman --- elevstl.cpp | 123 ++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 96 insertions(+), 27 deletions(-) diff --git a/elevstl.cpp b/elevstl.cpp index c09c10e..2c17fff 100644 --- a/elevstl.cpp +++ b/elevstl.cpp @@ -27,6 +27,8 @@ struct triangle{ char endTag[2] = {0,0}; ofstream out; +float longscale = 1.0; +bool vertical; //Determines the normal vector of a triangle from three vertices vertex normalOf(vertex p1, vertex p2, vertex p3){ @@ -46,14 +48,27 @@ vertex normalOf(vertex p1, vertex p2, vertex p3){ //Creates a vertex from x, y, and z coordinates vertex createVertex(float x, float y, float z){ vertex v; - v.x =x; + v.x = x * longscale; v.y = y; v.z = z; return v; } +vertex rotateVertex(vertex v){ + if (!vertical) + return v; + vertex r; + r.x = -v.z; + r.y = v.y; + r.z = v.x; + return r; +} + //Creates a triangle and its normal vector from three vertices -triangle createTriangle(vertex j, vertex k, vertex l){ +triangle createTriangle(vertex _j, vertex _k, vertex _l){ + vertex j = rotateVertex(_j); + vertex k = rotateVertex(_k); + vertex l = rotateVertex(_l); triangle t; t.a = j; t.b = k; @@ -86,10 +101,28 @@ void addTriangle(triangle t){ int width = 40; //default width and length of model int height = 40; +float thickness = 2; vector hList; string savefile = "stls/"; +float baseheight(float below) +{ + if (vertical) + return below - thickness; + else + return 0; +} + +vertex createVertexBelow(vertex v) +{ + vertex r; + r.x = v.x; + r.y = v.y; + r.z = baseheight(v.z); + return r; +} + //Takes a height array height array of variable length and turns it into an STL file void writeSTLfromArray(){ out.open(savefile.c_str(),ios_base::binary); @@ -97,7 +130,10 @@ void writeSTLfromArray(){ uint32_t triangleCount = (width-1)*(height-1)*2; //number of facets in a void-free surface triangleCount += 4*(width-1); //triangle counts for the walls of the model triangleCount += 4*(height-1); - triangleCount += 2; //base triangles + if (vertical) + triangleCount += (width-1)*(height-1)*2; //For bottom surface + else + triangleCount += 2; //base triangles float planarScale = 40/width; if(out.good()){ @@ -112,8 +148,8 @@ void writeSTLfromArray(){ vertex b = createVertex(c-1, 0,hList.at(c-1)); vertex d = createVertex(c-1, 1,hList.at(c+width-1)); - vertex w = createVertex(c,0,0); //used in model walls - vertex z = createVertex(c-1,0,0); + vertex w = createVertexBelow(a); //used in model walls + vertex z = createVertexBelow(b); addTriangle(createTriangle(a,d,b)); addTriangle(createTriangle(b,z,a)); //model walls @@ -129,8 +165,15 @@ void writeSTLfromArray(){ vertex b = createVertex(x,y-1,hList.at((y-1)*width+x)); vertex c = createVertex(x-1,y,hList.at(y*width+x-1)); addTriangle(createTriangle(a,c,b)); + + if(vertical){ + vertex a1 = createVertexBelow(a); + vertex b1 = createVertexBelow(b); + vertex c1 = createVertexBelow(c); + addTriangle(createTriangle(a1,b1,c1)); + } }else{ - triangleCount--; + triangleCount-=vertical?2:1; } } for(int x = 1; x < width; x++){ @@ -139,19 +182,26 @@ void writeSTLfromArray(){ vertex b = createVertex(x-1,y,hList.at(y*width+x-1)); vertex c = createVertex(x-1,y+1,hList.at((y+1)*width+x-1)); addTriangle(createTriangle(a,c,b)); + + if(vertical){ + vertex a1 = createVertexBelow(a); + vertex b1 = createVertexBelow(b); + vertex c1 = createVertexBelow(c); + addTriangle(createTriangle(a1,b1,c1)); + } }else{ - triangleCount--; + triangleCount-=vertical?2:1; } } } for(int x = 1; x < width; x++){ - if((int)hList.at((height-1)*width+x)!=-9 & (int)hList.at((height-2)*width+x)!=-99 & (int)hList.at((height-1)*width+x-1)!=99 ){ + if((int)hList.at((height-1)*width+x)!=-99 & (int)hList.at((height-2)*width+x)!=-99 & (int)hList.at((height-1)*width+x-1)!=99 ){ vertex a = createVertex(x,height-1,hList.at((height-1)*width+x)); //same vertex b = createVertex(x,height-2,hList.at((height-2)*width+x)); vertex c = createVertex(x-1,height-1,hList.at((height-1)*width+x-1)); - vertex w = createVertex(x,height-1,0); //used in model walls - vertex z = createVertex(x-1,height-1,0); + vertex w = createVertexBelow(a); //used in model walls + vertex z = createVertexBelow(c); addTriangle(createTriangle(a,c,b)); addTriangle(createTriangle(c,a,z)); //model walls @@ -161,32 +211,32 @@ void writeSTLfromArray(){ } } - for(int y = 1; y < width; y++){ //adds walls in the y direction for + for(int y = 1; y < height; y++){ //adds walls in the y direction for vertex st = createVertex(0,y,hList.at(y*width)); //for x=0 first vertex sb = createVertex(0,y-1,hList.at((y-1)*width)); - vertex bt = createVertex(0,y,0); - vertex bb = createVertex(0,y-1,0); + vertex bt = createVertexBelow(st); + vertex bb = createVertexBelow(sb); addTriangle(createTriangle(bb,sb,st)); addTriangle(createTriangle(st,bt,bb)); st = createVertex(width-1,y,hList.at(y*width+width-1)); //for x=width next sb = createVertex(width-1,y-1,hList.at(y*width-1)); - bt = createVertex(width-1,y,0); - bb = createVertex(width-1,y-1,0); + bt = createVertexBelow(st); + bb = createVertexBelow(sb); addTriangle(createTriangle(sb,bb,st)); addTriangle(createTriangle(bt,st,bb)); } - - vertex origin = createVertex(0,0,0); - vertex bottomright = createVertex(width-1,0,0); - vertex topleft = createVertex(0,height-1,0); - vertex topright = createVertex(width-1,height-1,0); - addTriangle(createTriangle(origin,topright,bottomright)); - addTriangle(createTriangle(origin,topleft,topright)); - - + + if(!vertical){ + vertex origin = createVertex(0,0,0); + vertex bottomright = createVertex(width-1,0,0); + vertex topleft = createVertex(0,height-1,0); + vertex topright = createVertex(width-1,height-1,0); + addTriangle(createTriangle(origin,topright,bottomright)); + addTriangle(createTriangle(origin,topleft,topright)); + } out.seekp(80); out.write((char *)&triangleCount,4); } @@ -214,12 +264,17 @@ int main(int argc, char **argv) //lat, long, res printf("Bad resolution\n"); return 0; } + vertical = true; res = res/3; //SRTM data already has a resolution of 3 arcseconds printf("Using resolution %i\n",res); //Arc - height = 40*res; - width = 40*res; + + // Mathematical expression: Length of a degree of longitude = cos(latitude) * 111.325 kilometers + longscale = cos(lat*0.0174532925199433); // ratio of long to lat + printf ("Longitude scaling factor (at NW corner): %f\n", longscale); + + height = 40*res; + width = 40*res / longscale; hList.resize(width*height,0); - savefile.append(std::string(argv[4])); printf("Saving to '%s'\n",savefile.c_str()); @@ -286,6 +341,20 @@ int main(int argc, char **argv) //lat, long, res hList.at((height-1-y)*width+x) = h/(verticalscale)+1; //+1 so that the bottom of the model does not bleed through to the top } } + for(int y = 0; y < height; y++){ + for(int x = 0; x < width; x++){ + int index = (height-1-y)*width+x; + h = hList.at(index); + if (h == -99 && x && y) + { + // Ideally you would take average of surrounding 4 (or 8) non-void points + // This actually only looks at previous x/y (which are already corrected). + // Will go a bit wrong if you are unlucky enough to have a void when x or y is zero. + h = (hList.at(index-1) + hList.at(index+width))/2; + hList.at(index) = h; + } + } + } } tile.close(); writeSTLfromArray(); From f788cbd87c037f382299f3f5604ac42a293f186f Mon Sep 17 00:00:00 2001 From: Richard Chapman Date: Sat, 13 Dec 2014 17:55:20 +0000 Subject: [PATCH 2/4] Add argument parsing Separate crust mode and verticall mode Signed-off-by: Richard Chapman --- elevstl.cpp | 93 +++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 68 insertions(+), 25 deletions(-) diff --git a/elevstl.cpp b/elevstl.cpp index 2c17fff..0c191a4 100644 --- a/elevstl.cpp +++ b/elevstl.cpp @@ -29,7 +29,9 @@ char endTag[2] = {0,0}; ofstream out; float longscale = 1.0; bool vertical; - +bool crust; +float minheight = 1e30; +; //Determines the normal vector of a triangle from three vertices vertex normalOf(vertex p1, vertex p2, vertex p3){ vertex u,v,r; @@ -108,10 +110,10 @@ string savefile = "stls/"; float baseheight(float below) { - if (vertical) + if (crust) return below - thickness; else - return 0; + return minheight - thickness; } vertex createVertexBelow(vertex v) @@ -130,7 +132,7 @@ void writeSTLfromArray(){ uint32_t triangleCount = (width-1)*(height-1)*2; //number of facets in a void-free surface triangleCount += 4*(width-1); //triangle counts for the walls of the model triangleCount += 4*(height-1); - if (vertical) + if (crust) triangleCount += (width-1)*(height-1)*2; //For bottom surface else triangleCount += 2; //base triangles @@ -166,14 +168,14 @@ void writeSTLfromArray(){ vertex c = createVertex(x-1,y,hList.at(y*width+x-1)); addTriangle(createTriangle(a,c,b)); - if(vertical){ + if(crust){ vertex a1 = createVertexBelow(a); vertex b1 = createVertexBelow(b); vertex c1 = createVertexBelow(c); addTriangle(createTriangle(a1,b1,c1)); } }else{ - triangleCount-=vertical?2:1; + triangleCount-=crust?2:1; } } for(int x = 1; x < width; x++){ @@ -183,14 +185,14 @@ void writeSTLfromArray(){ vertex c = createVertex(x-1,y+1,hList.at((y+1)*width+x-1)); addTriangle(createTriangle(a,c,b)); - if(vertical){ + if(crust){ vertex a1 = createVertexBelow(a); vertex b1 = createVertexBelow(b); vertex c1 = createVertexBelow(c); addTriangle(createTriangle(a1,b1,c1)); } }else{ - triangleCount-=vertical?2:1; + triangleCount-=crust?2:1; } } } @@ -229,7 +231,7 @@ void writeSTLfromArray(){ addTriangle(createTriangle(bt,st,bb)); } - if(!vertical){ + if(!crust){ vertex origin = createVertex(0,0,0); vertex bottomright = createVertex(width-1,0,0); vertex topleft = createVertex(0,height-1,0); @@ -243,7 +245,26 @@ void writeSTLfromArray(){ cout << triangleCount << "\n"; } -int main(int argc, char **argv) //lat, long, res +float verticalscale = 23.2; // true verticalscale 92.7 meters/arcsecond (at equator) gives models that are too flat to be interesting + +const char ** checkOption(const char *arg, const char **argv) +{ + if (strcmp(arg, "--vertical")==0) + vertical = true; + else if (strcmp(arg, "--crust")==0) + crust = true; + else if (strcmp(arg, "--thickness")==0) + thickness = atoi(*argv++); + else if (strcmp(arg, "--scale")==0) + verticalscale = atof(*argv++); + else{ + printf("Unrecognised argument %s\n", arg); + exit(0); + } + return argv; +} + +int main(int argc, const char **argv) //lat, long, res, outfile { string file; int point; @@ -252,19 +273,41 @@ int main(int argc, char **argv) //lat, long, res int res; int tile_n; int tile_w; - //float true_verticalscale = 92.7; //meters/arcsecond at equator - float verticalscale = 23.2; //true_verticalscale gives models that are too flat to be interesting - - lat = atof(argv[1]); //Latitude of NW corner - printf("Using latitude: %f\n",lat); - lng = atof(argv[2]); //Longitude of NW corner - printf("Using longitude: %f\n",lng); - res = atoi(argv[3]); //arcseconds/tick in model - if(res%3!=0){ //must be a multiple of 3 - printf("Bad resolution\n"); - return 0; + if(argc < 5){ + printf("Not enough arguents\n"); + exit(0); + } + int argno = 1; + argv++; // not interested in argv[0] + while(*argv){ + const char *arg = *argv++; + if(arg[0]=='-' && arg[1]=='-') + argv = checkOption(arg, argv); + else switch(argno++){ + case 1: + lat = atof(arg); //Latitude of NW corner + printf("Using latitude: %f\n",lat); + break; + case 2: + lng = atof(arg); //Longitude of NW corner + printf("Using longitude: %f\n",lng); + break; + case 3: + res = atoi(arg); //arcseconds/tick in model + if(res%3!=0){ //must be a multiple of 3 + printf("Bad resolution\n"); + return 0; + } + break; + case 4: + savefile.append(std::string(arg)); + printf("Saving to '%s'\n",savefile.c_str()); + break; + default: + printf("Too many arguments\n"); + exit(0); + } } - vertical = true; res = res/3; //SRTM data already has a resolution of 3 arcseconds printf("Using resolution %i\n",res); //Arc @@ -275,8 +318,6 @@ int main(int argc, char **argv) //lat, long, res height = 40*res; width = 40*res / longscale; hList.resize(width*height,0); - savefile.append(std::string(argv[4])); - printf("Saving to '%s'\n",savefile.c_str()); //-------get correct tile------------------- if(lat>=0){ //Positive is north @@ -338,7 +379,7 @@ int main(int argc, char **argv) //lat, long, res } //rotate model to correct orientation //hList.at((height-1-y)*width+x) = h/(verticalscale*res); //cast verticalscale to int for COOl effect! - hList.at((height-1-y)*width+x) = h/(verticalscale)+1; //+1 so that the bottom of the model does not bleed through to the top + hList.at((height-1-y)*width+x) = h/(verticalscale)+1; //+1 so that the bottom of the model does not bleed through to the top } } for(int y = 0; y < height; y++){ @@ -353,6 +394,8 @@ int main(int argc, char **argv) //lat, long, res h = (hList.at(index-1) + hList.at(index+width))/2; hList.at(index) = h; } + if (h != -99 && h < minheight) + minheight = h; } } } From 05b726ba52a85ca861994a4ba0014bd9af688880 Mon Sep 17 00:00:00 2001 From: Richard Chapman Date: Sat, 13 Dec 2014 18:11:23 +0000 Subject: [PATCH 3/4] Better hole fixing More options Signed-off-by: Richard Chapman --- elevstl.cpp | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/elevstl.cpp b/elevstl.cpp index 0c191a4..091160a 100644 --- a/elevstl.cpp +++ b/elevstl.cpp @@ -103,7 +103,7 @@ void addTriangle(triangle t){ int width = 40; //default width and length of model int height = 40; -float thickness = 2; +float thickness = 0; vector hList; string savefile = "stls/"; @@ -112,8 +112,10 @@ float baseheight(float below) { if (crust) return below - thickness; - else + else if (thickness) return minheight - thickness; + else + return 0; } vertex createVertexBelow(vertex v) @@ -232,10 +234,10 @@ void writeSTLfromArray(){ } if(!crust){ - vertex origin = createVertex(0,0,0); - vertex bottomright = createVertex(width-1,0,0); - vertex topleft = createVertex(0,height-1,0); - vertex topright = createVertex(width-1,height-1,0); + vertex origin = createVertex(0,0,minheight); + vertex bottomright = createVertex(width-1,0,minheight); + vertex topleft = createVertex(0,height-1,minheight); + vertex topright = createVertex(width-1,height-1,minheight); addTriangle(createTriangle(origin,topright,bottomright)); addTriangle(createTriangle(origin,topleft,topright)); } @@ -308,6 +310,8 @@ int main(int argc, const char **argv) //lat, long, res, outfile exit(0); } } + if (crust && !thickness) + thickness = 2; res = res/3; //SRTM data already has a resolution of 3 arcseconds printf("Using resolution %i\n",res); //Arc @@ -386,15 +390,24 @@ int main(int argc, const char **argv) //lat, long, res, outfile for(int x = 0; x < width; x++){ int index = (height-1-y)*width+x; h = hList.at(index); - if (h == -99 && x && y) + if (h == -99) { // Ideally you would take average of surrounding 4 (or 8) non-void points // This actually only looks at previous x/y (which are already corrected). - // Will go a bit wrong if you are unlucky enough to have a void when x or y is zero. - h = (hList.at(index-1) + hList.at(index+width))/2; + // Will go a bit wrong if you are unlucky enough to have a void when x and y is zero. + if (x && y) + h = (hList.at(index-1) + hList.at(index+width))/2; + else if (x) + h = hList.at(index-1); + else if (y) + h = hList.at(index+width); + else{ + printf("Missing data at NW corner - aborting\n"); + exit(0); + } hList.at(index) = h; } - if (h != -99 && h < minheight) + if (h < minheight) minheight = h; } } From 1fc1273047b716c47427d9a2e127d8158d2a5c6c Mon Sep 17 00:00:00 2001 From: Richard Chapman Date: Sat, 13 Dec 2014 19:50:12 +0000 Subject: [PATCH 4/4] Added support for spanning lat/long tile boundaries Signed-off-by: Richard Chapman --- elevstl.cpp | 170 +++++++++++++++++++++++++++------------------------- 1 file changed, 89 insertions(+), 81 deletions(-) diff --git a/elevstl.cpp b/elevstl.cpp index 091160a..69e036b 100644 --- a/elevstl.cpp +++ b/elevstl.cpp @@ -105,6 +105,11 @@ int width = 40; //default width and length of model int height = 40; float thickness = 0; vector hList; +ifstream tiles[4]; +int ytiles; +int xtiles; +float lat; +float lng; string savefile = "stls/"; @@ -266,15 +271,44 @@ const char ** checkOption(const char *arg, const char **argv) return argv; } -int main(int argc, const char **argv) //lat, long, res, outfile +void checkOpen(int tilex, int tiley) { + int tidx = tilex + tiley*ytiles; + if (tiles[tidx].is_open()) + return; string file; - int point; - float lat; - float lng; + //-------get correct tile------------------- + float tlat = lat - tiley; + float tlng = lng + tilex; + if(tlat>=0){ //Positive is north + file = "N"; + }else{ + file = "S"; + } + file.append( to_string( abs( (int)floor(tlat) ) ) ); + if(file.length()==2) file.insert(1,"0"); + if(tlng>=0){ //Positive is east + file.append("E"); + }else{ + file.append("W"); + } + file.append( to_string( abs( (int)floor(tlng) ) ) ); + if(file.length()==5)file.insert(4,"00"); + if(file.length()==6)file.insert(4,"0"); + + file.append(".hgt"); + file.insert(0,"./hgt_files/"); + printf("Opening %s\n",file.c_str()); + tiles[tidx].open(file.c_str(),ios::in|ios::binary); + if (!tiles[tidx].is_open()){ + printf ("Error opening file %s", file.c_str()); + exit(0); + } +} + +int main(int argc, const char **argv) //lat, long, res, outfile +{ int res; - int tile_n; - int tile_w; if(argc < 5){ printf("Not enough arguents\n"); exit(0); @@ -323,27 +357,6 @@ int main(int argc, const char **argv) //lat, long, res, outfile width = 40*res / longscale; hList.resize(width*height,0); - //-------get correct tile------------------- - if(lat>=0){ //Positive is north - file = "N"; - }else{ - file = "S"; - } - file.append( to_string( abs( (int)floor(lat) ) ) ); - if(file.length()==2) file.insert(1,"0"); - if(lng>=0){ //Positive is east - file.append("E"); - }else{ - file.append("W"); - } - file.append( to_string( abs( (int)floor(lng) ) ) ); - if(file.length()==5)file.insert(4,"00"); - if(file.length()==6)file.insert(4,"0"); - - file.append(".hgt"); - file.insert(0,"./hgt_files/"); - puts(file.c_str()); - //-------Find starting file index--------------- float n = (lat-floor(lat))*3600; float e = (lng-floor(lng))*3600; @@ -351,68 +364,63 @@ int main(int argc, const char **argv) //lat, long, res, outfile int i = 1201-(int)(n/(3)); int j = (int)(e/(3)); - if(i+height>1201){ - perror("Error in y"); - } - if(j+width>1201){ - perror("Error in x"); - } - - point = j+i*1201; //the file index of the NW corner - //------------Open file and read data into array---------------------------- - ifstream tile; - tile.open(file.c_str(),ios::in|ios::binary); int h; - if (!tile.is_open()){ - perror ("Error opening file"); - }else{ - char number [2]; - for(int y = 0; y < height; y++){ - for(int x = 0; x < width; x++){ - tile.seekg((point+x+y*1201)*2,ios::beg); - tile.read(number,2); - h = number[1]; - if(h<0){ - h = h+255; - } - h+= number[0]*256; - //If a void exists, marks it as -100 - if(h<-100){ - h=-verticalscale*100; - } - //rotate model to correct orientation - //hList.at((height-1-y)*width+x) = h/(verticalscale*res); //cast verticalscale to int for COOl effect! - hList.at((height-1-y)*width+x) = h/(verticalscale)+1; //+1 so that the bottom of the model does not bleed through to the top + ytiles = (i+height) / 1200; + xtiles = (j+width) / 1200; + //tiles.resize(ytiles*xtiles); + char number [2]; + for(int y = i; y < i+height; y++){ + int tiley = y / 1200; + int coly = y % 1200; + int hy = y - i; + for(int x = j; x < j+width; x++) { + int tilex = x / 1200; + int colx = x % 1200; + int hx = x - j; + checkOpen(tilex, tiley); + tiles[tiley*ytiles+tilex].seekg((colx+coly*1201)*2,ios::beg); + tiles[tiley*ytiles+tilex].read(number,2); + h = number[1]; + if(h<0){ + h = h+255; + } + h+= number[0]*256; + //If a void exists, marks it as -100 + if(h<-100){ + h=-verticalscale*100; } + //rotate model to correct orientation + //hList.at((height-1-hy)*width+hx) = h/(verticalscale*res); //cast verticalscale to int for COOl effect! + hList.at((height-1-hy)*width+hx) = h/(verticalscale)+1; //+1 so that the bottom of the model does not bleed through to the top } - for(int y = 0; y < height; y++){ - for(int x = 0; x < width; x++){ - int index = (height-1-y)*width+x; - h = hList.at(index); - if (h == -99) - { - // Ideally you would take average of surrounding 4 (or 8) non-void points - // This actually only looks at previous x/y (which are already corrected). - // Will go a bit wrong if you are unlucky enough to have a void when x and y is zero. - if (x && y) - h = (hList.at(index-1) + hList.at(index+width))/2; - else if (x) - h = hList.at(index-1); - else if (y) - h = hList.at(index+width); - else{ - printf("Missing data at NW corner - aborting\n"); - exit(0); - } - hList.at(index) = h; + } + for(int y = 0; y < height; y++){ + for(int x = 0; x < width; x++){ + int index = (height-1-y)*width+x; + h = hList.at(index); + if (h == -99) + { + // Ideally you would take average of surrounding 4 (or 8) non-void points + // This actually only looks at previous x/y (which are already corrected). + // Will go a bit wrong if you are unlucky enough to have a void when x and y is zero. + if (x && y) + h = (hList.at(index-1) + hList.at(index+width))/2; + else if (x) + h = hList.at(index-1); + else if (y) + h = hList.at(index+width); + else{ + printf("Missing data at NW corner - aborting\n"); + exit(0); } - if (h < minheight) - minheight = h; + hList.at(index) = h; } + if (h < minheight) + minheight = h; } } - tile.close(); + // tile.close(); writeSTLfromArray(); return 0; }