|
|
Erstellen einer Höhenmatrix aus den GIS-Daten
Da uns nur Satellitendaten zur Verfügung standen, die entweder Shape Type 1 oder 3 hatten, konnten wir nicht einfach die Vektoren
aus der Shapefile erstellen, sondern mussten aus der Shapefile den 2D Teil der Vektoren ziehen und schließlich das Attribut "Höhe" hinzufügen.
Wie im Teil "Shapelib" beschrieben haben wir also zuerst einmal die dBASE, die SHP und die SHX File aus dm HIGHT Ordner
in C++ eingelesen und hatten darauf dann eine Punktwolke von 21000 Punkten in einem Feld mit insgesamt 504000000 Punkten.
Dass damit eine Grafikengine, die nur Höhen setzt und fehlende Höhen einfach auf NULL setzt nichts anfangen kann, ist einleuchtend, aber was tun?
Wir haben also nach einer Möglichkeit gesucht, alle Höhen zu bekommen und sind auf das das Praktikum "Terraindigitalisierung" gestoßen.
In
dem ursprünglichen Programm wurden nun die Teile 2 und drei sinnlos. Unser Programm startet also so:
#define MAXLINESIZE 300
#define MAX_CORRECTION_RADIUS 20
#define POINT_DENSITY_AT_BORDER 10
#include <fstream>
#include <iostream>
#include "tiffio.h"
#include <string.h>
#include <math.h>
#include <time.h>
#include "terrain.h"
#include "shapefil.h"
using namespace std;
#ifdef WIN32 /* Windows */
#define TRIANGLE_EXECUTE "triangle.exe vertex -I -N"
#endif
#ifdef LINUX /* Linux */
#define TRIANGLE_EXECUTE "./triangle vertex -I -N"
#endif
#ifdef MACOS /* MacOS */
#define TRIANGLE_EXECUTE "./triangle vertex -I -N"
#endif
int cmpxminfirst(const void*, const void*);
inline unsigned short seekIndex(unsigned short*,unsigned
short*,unsigned long,uint32);
void computeBoundingBox(T_triangle*);
int checkParameter(const char*, int, char**);
int main(int argc, char *argv[])
{
TIFF *output;
uint16 compression, *buffer16;
char lineBuffer[MAXLINESIZE];
unsigned long vertexA, vertexB, vertexC, zeile, spalte;
unsigned short *index_maxx, *index_count, countmin, cache_count;
unsigned char *freq;
double parameter1, parameter2, temp;
long count, pixelNumber, countTriangles, countVertices,i;
bool triangleFound;
FILE *vertexnode, *vertexele, *heightmapasc;
T_point* mypoints;
T_triangle* mytriangles;
clock_t prgStart, prgEnd, startTime, endTime;
prgStart=clock();
Bis dahin ist alles, bis auf einige Variablen genau gleich zum ursprünglichen Code, also nehme ich an, dass bisher nichts zu erklären ist.
Aber nun steigen wir erst in Part 4 des ursprünglichen Programms ein:
const char* name="Contour";
double min[4],max[4];
int entities;
int cp;
cout << "opening shapefiles...\n";
SHPHandle shape=SHPOpen(name,"rb");
printf("reading shapefiles...\n");
SHPGetInfo(shape,&entities,NULL,min,max);
printf("entities = %d.\n",entities);
cp=0;
DBFHandle dbf=DBFOpen(name,"rb");
int x,y;
double r=4;
double* z;
int width = (int)((max[0]-min[0])/r)+1, length = (int)((max[1]-min[1])/r)+1;
pixelNumber = width * length;
Wir greifen auf die Shapelib zurück, um die Shapefile Contour einzulesen und
definieren die zwei Arrays min[] und max[] der Länge 4, in die später die
Grenzen der Bounding Boxes geschrieben werden, d.h. min[]={Xmin, Ymin, Zmin, Mmin} und max[] entsprechend.
Die Daten kann man ganz einfach aus dem Header der SHP-Datei auslesen.
Da wir eine Shapefile des Typs 3 (Polyline) haben, werden in den letzten beiden Einträgen Nuller stehen.
entities enthält die Anzahl der Einträge, die die Datei hat.
Zuletzt wird noch
die Dimension des Terrains berechnet, wobei r einfach ein Faktor ist, um die Rechenzeit zu verkürzen.
(Bei r=2 hatten wir schon eine Rechenzeit von 7,4 h!!)
if((z = (double *) malloc(pixelNumber*sizeof(double)))
== NULL){
fprintf(stderr, "Could not allocate enough memory for the
uncompressed image\n");
system ("PAUSE");
exit(EXIT_FAILURE);
}
if((freq = (unsigned char *) malloc(pixelNumber)) == NULL){
fprintf(stderr, "Could not allocate enough memory for freq\n");
system ("PAUSE");
exit(EXIT_FAILURE);
}
for(i=0;i<pixelNumber;i++)
{
z[i]=0;
freq[i]=0;
}
Hier wird der Speicher für z reserviert, was später die Höhenmatrix
ist. Desweiteren
wird für freq Speicher reserviert, was später die Häufigkeit beschreibt, wie oft verschiedene Punkte auf einen durch das
Runterskalieren wegen den Integerquotienten 1/r vorkommen. Die letzte For-Schleife initialisiert die beiden Speicherbereiche.
for(i = 0; i < entities; i++)
{
SHPObject* object=SHPReadObject(shape,i);
for(int j=0; j<object->nVertices;j++)
{
x=(int)((object->padfX[j]-min[0])/r);
y=(int)((object->padfY[j]-min[1])/r);
// für die Heightmap muss man hier (max[1]-object->padfY[j])/r
nehmen, denn sonst erscheint die Heightmap vertikal gespiegelt (Weil
bei tiff der Ursprung oben links ist und nicht unten rechts)
if(x>=0 && x<width && y>=0 && y<width)
{
freq[y*width+x]++;
z[y*width+x]+=DBFReadDoubleAttribute(dbf,i,0);
}
}
SHPDestroyObject(object);
}
SHPClose(shape);
DBFClose(dbf);
Nun werden einfach die Koordinaten x,y und die
zugehörige Höhe ausgelesen, wobei bei überschneidung von Punkten
freq[] erhöht wird, so dass freq[] genau die Anzahl der
ursprünglichen Punkte beschreibt, die durch den Integerquotienten
1/r aufeinander fallen.
z[] enthält zunächst die Summe der Höhen aller
aufeinander liegender Punkte. for(i = 0; i < pixelNumber; i++)
if(z[i]!=0) cp++;
if((vertexnode = fopen("vertex.node","wt")) == NULL) {
fprintf(stderr, "Coult not open vertex.node for writing data");
system ("PAUSE");
exit(EXIT_FAILURE);
}
startTime = clock();
printf("writing vertex.node ... ");
fprintf(vertexnode,"# <Number of vertices> <dimensions> <number
attributes> <number of boundary markers>\n");
fprintf(vertexnode,"%lu 2 1 0\n",cp);
Zunächst noch eine Initialisierung und ein Abfangen von möglichen Fehlern und dann wird
zuerst vertex.node für das Programm triangle.exe
beschrieben. Die Datei muss aus vier Spalten bestehen. Die erste
Spalte muss die Anzahl der zu triangularisierenden Eckpunkte
beschreiben, wir schreiben hier cp, was einfach für "count points" steht.
Die nächsten Spalten enthalten allgemeine Informationen (Dimension, # Attribute und # Grenzmarkierungen). Da wir ein zwei dimensionales Feld haben, das ein Attribut enthält (die Höhe) und keine Grenzen, schreiben wir also die Werte "Zwei, Eins, Null".
Diese Informationen bilden quasi den Head von vertex.node.
fprintf(vertexnode,"#\n# <PointNumber> <x-coordinate [pixel]>
<y-coordinate [pixel]>\n# pixel count starts with index 1\n#\n");
if((buffer16 = (uint16 *) malloc(pixelNumber*2)) == NULL){
fprintf(stderr, "Could not allocate enough memory for the
uncompressed image\n");
system ("PAUSE");
exit(EXIT_FAILURE);
}
cp=0;
for(i = 0; i < pixelNumber; i++)
{
if(z[i]!=0)
{
cp++; buffer16[i]=(uint16)((z[i]/freq[i])*64);
fprintf(vertexnode,"%lu\t%lu\t%lu\t%hu\n", cp, (i%width)+1, i/width+1,
buffer16[i]);
}
else buffer16[i]=0;
}
free(z);
free(freq);
fclose(vertexnode);
endTime = clock();
Die eigentlichen Daten werden darunter geschrieben, in die Spalten "Pointnumber, x-Koord., y-Koord. und Attribut"
Wir schreiben diese Daten in der folgenden For-Schleife. Zuvor muss
aber noch einmal Speicher reserviert werden und zwar für buffer16[],
dies wird später die Höhenmatrix mit den richtigen Höhen.
Das Array buffer16[] ist wie folgt aufgebaut: Da wir ein zweidimensionales Feld
beschreiben müssen (jedem Punkt soll eine Höhe zugeordnet werden), gehen wir zeilenweise vor.
Wir erhalten folgende Form:
buffer16[]={h((0,0)),h((1,0)),h((2,0)),...,h((with,0)),h((0,1)),h((1,1)),...,h((width,length))}
Daher kommen auch die beiden Formeln für width und heigth.
Der Rest des Programms ist gegenüber dem Original nicht weiter abgeändert und
somit verweise ich auf die Praktikumsseite für weitere Informationen.
|
