//add all lines from poLineLayer to collection      OGRGeometryCollection *gc = new OGRGeometryCollection();      poLineLayer->ResetReading();      while (OGRFeature *feature = poLineLayer->GetNextFeature())          gc->addGeometry(feature->GetGeometryRef());      //polygonize lines      CPLDebug( "OGR_ILI", "Polygonizing layer %s with %d multilines", poAreaLayer->GetLayerDefn()->GetName(), gc->getNumGeometries());      OGRMultiPolygon* polys = Polygonize( gc );      //associate polygon feature with data row according to centroid      int i;      OGRPolygon emptyPoly;#if defined(POLYGONIZE_AREAS)      GEOSGeom *ahInGeoms = NULL;      ahInGeoms = (GEOSGeom *) CPLCalloc(sizeof(void*),polys->getNumGeometries());      for( i = 0; i < polys->getNumGeometries(); i++ )      {          ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS();          if (!GEOSisValid(ahInGeoms[i])) ahInGeoms[i] = NULL;      }      poAreaLayer->ResetReading();      while (OGRFeature *feature = poAreaLayer->GetNextFeature())      {        GEOSGeom point = (GEOSGeom)feature->GetGeometryRef()->exportToGEOS();        for (i = 0; i < polys->getNumGeometries(); i++ )        {          if (ahInGeoms[i] && GEOSWithin(point, ahInGeoms[i]))          {            feature->SetGeometry( polys->getGeometryRef(i) );            break;          }        }        if (i == polys->getNumGeometries())        {          CPLDebug( "OGR_ILI", "Association between area and point failed.");          feature->SetGeometry( &emptyPoly );        }        GEOSGeom_destroy( point );      }      poAreaLayer->GetLayerDefn()->SetGeomType(wkbPolygon);      for( i = 0; i < polys->getNumGeometries(); i++ )          GEOSGeom_destroy( ahInGeoms[i] );      CPLFree( ahInGeoms );#endif    }}void ILI1Reader::JoinSurfaceLayers(){    for(int iLayer = 0; iLayer < nSurfaceLayers; iLayer++ )    {      OGRLayer *poSurfaceLayer = papoSurfaceLayers[iLayer];      OGRLayer *poPolyLayer = papoSurfacePolyLayers[iLayer];      poSurfaceLayer->GetLayerDefn()->SetGeomType(poPolyLayer->GetLayerDefn()->GetGeomType());      poSurfaceLayer->ResetReading();      poPolyLayer->ResetReading();      while (OGRFeature *feature = poSurfaceLayer->GetNextFeature()) {        //Assume same sequence -> Dangerous?        OGRFeature *polyfeature = poPolyLayer->GetNextFeature();        if (polyfeature) { //&& EQUAL(feature->GetFieldAsString(0), polyfeature->GetFieldAsString(0))          feature->SetGeometry(polyfeature->GetGeometryRef());        }      }    }}int ILI1Reader::ReadTable() {        char **tokens = NULL;    const char *firsttok = NULL;    int ret = TRUE;    int warned = FALSE;    int fIndex;    OGRFeatureDefn *featureDef = curLayer->GetLayerDefn();    OGRFieldDefn *fieldDef = NULL;    OGRFeature *feature = NULL;    long fpos = VSIFTell(fpItf);    while (ret && (tokens = ReadParseLine()))    {      firsttok = CSLGetField(tokens, 0);      if (EQUAL(firsttok, "OBJE"))      {        if (featureDef->GetFieldCount() == 0)        {          CPLDebug( "OGR_ILI", "No field definition found for table: %s", featureDef->GetName() );          //Model not read - use heuristics          for (fIndex=1; fIndex<CSLCount(tokens); fIndex++)          {            fieldDef = new OGRFieldDefn(CPLStrdup("Field00"), OFTString);            *(char *)(fieldDef->GetNameRef()+strlen(fieldDef->GetNameRef())-2) = '0'+fIndex/10;            *(char *)(fieldDef->GetNameRef()+strlen(fieldDef->GetNameRef())-1) = '0'+fIndex%10;            featureDef->AddFieldDefn(fieldDef);          }        }        feature = new OGRFeature(featureDef);        int fieldno = 0;        for (fIndex=1; fIndex<CSLCount(tokens) && fieldno < featureDef->GetFieldCount(); fIndex++, fieldno++)        {          if (!EQUAL(tokens[fIndex], "@")) {            //CPLDebug( "OGR_ILI", "Adding Field %d: %s", fieldno, tokens[fIndex]);            feature->SetField(fieldno, CPLStrdup(tokens[fIndex]));            if (featureDef->GetFieldDefn(fieldno)->GetType() == OFTReal                && fieldno > 0                && featureDef->GetFieldDefn(fieldno-1)->GetType() == OFTReal                && featureDef->GetGeomType() == wkbPoint) {              //add Point geometry              OGRPoint *ogrPoint = new OGRPoint(atof(tokens[fIndex-1]), atof(tokens[fIndex]));              feature->SetGeometry(ogrPoint);            }          }        }        if (!warned && featureDef->GetFieldCount() != CSLCount(tokens)-1) {          CPLDebug( "OGR_ILI", "Field count doesn't match. %d declared, %d found", featureDef->GetFieldCount(), CSLCount(tokens)-1);          warned = TRUE;        }        curLayer->AddFeature(feature);      }      else if (EQUAL(firsttok, "STPT"))      {        OGRGeometry *geom = ReadGeom(tokens, featureDef->GetGeomType());        if (EQUAL(featureDef->GetFieldDefn(featureDef->GetFieldCount()-1)->GetNameRef(), "ILI_Geometry"))        {          AddIliGeom(feature, featureDef->GetFieldCount()-1, fpos);        }        feature->SetGeometry(geom);      }      else if (EQUAL(firsttok, "ELIN"))      {        //empty geom      }      else if (EQUAL(firsttok, "EDGE"))      {        tokens = ReadParseLine(); //STPT        OGRGeometry *geom = ReadGeom(tokens, wkbMultiLineString);        feature->SetGeometry(geom);        if (EQUAL(featureDef->GetFieldDefn(featureDef->GetFieldCount()-1)->GetNameRef(), "ILI_Geometry"))        {          AddIliGeom(feature, featureDef->GetFieldCount()-1, fpos);        }      }      else if (EQUAL(firsttok, "PERI"))      {      }      else if (EQUAL(firsttok, "ETAB"))      {        return TRUE;      }      else      {        CPLDebug( "OGR_ILI", "Unexpected token: %s", firsttok );      }            CSLDestroy(tokens);      fpos = VSIFTell(fpItf);    }        return ret;}OGRGeometry *ILI1Reader::ReadGeom(char **stgeom, OGRwkbGeometryType eType) {        char **tokens = NULL;    const char *firsttok = NULL;    int end = FALSE;    OGRGeometry *ogrGeom = NULL;    OGRLineString *ogrLine = NULL; //current line    int isArc = FALSE;    OGRPoint ogrPoint, arcPoint, endPoint; //points for arc interpolation    OGRMultiLineString *ogrMultiLine = NULL; //current multi line        //tokens = ["STPT", "1111", "22222"]    ogrPoint.setX(atof(stgeom[1])); ogrPoint.setY(atof(stgeom[2]));    ogrLine = new OGRLineString();    ogrLine->addPoint(&ogrPoint);    if (eType == wkbMultiLineString || eType == wkbGeometryCollection)    {      ogrMultiLine = new OGRMultiLineString();    }    while (!end && (tokens = ReadParseLine()))    {      firsttok = CSLGetField(tokens, 0);      if (EQUAL(firsttok, "LIPT"))      {        if (isArc) {          endPoint.setX(atof(tokens[1])); endPoint.setY(atof(tokens[2]));          interpolateArc(ogrLine, &ogrPoint, &arcPoint, &endPoint, arcIncr);        }        ogrPoint.setX(atof(tokens[1])); ogrPoint.setY(atof(tokens[2])); isArc = FALSE;        ogrLine->addPoint(&ogrPoint);      }      else if (EQUAL(firsttok, "ARCP"))      {        isArc = TRUE;        arcPoint.setX(atof(tokens[1])); arcPoint.setY(atof(tokens[2]));      }      else if (EQUAL(firsttok, "ELIN"))      {        if (ogrMultiLine)        {          ogrMultiLine->addGeometryDirectly(ogrLine);        }        if (eType != wkbGeometryCollection) end = TRUE;      }      else if (EQUAL(firsttok, "STPT"))      {        //AREA lines spread over mutltiple objects        ogrPoint.setX(atof(tokens[1])); ogrPoint.setY(atof(tokens[2])); isArc = FALSE;        ogrLine = new OGRLineString();        ogrLine->addPoint(&ogrPoint);      }      else if (EQUAL(firsttok, "EEDG"))      {        end = TRUE;      }      else if (EQUAL(firsttok, "LATT"))      {        //Line Attributes (ignored)      }      else if (EQUAL(firsttok, "EFLA"))      {        end = TRUE;      }      else if (EQUAL(firsttok, "ETAB"))      {        end = TRUE;      }      else if (EQUAL(firsttok, "OBJE"))      {        //AREA lines spread over mutltiple objects      }      else      {        CPLDebug( "OGR_ILI", "Unexpected token: %s", firsttok );      }          CSLDestroy(tokens);    }        if (eType == wkbPolygon)    {      OGRPolygon *gc = new OGRPolygon();      gc->addRing((OGRLinearRing *)ogrLine);      ogrGeom = gc;    }    else if (ogrMultiLine)    {      ogrGeom = ogrMultiLine;    }    else    {      ogrGeom = ogrLine;    }    return ogrGeom;}/************************************************************************//*                              AddLayer()                              *//************************************************************************/void ILI1Reader::AddLayer( OGRLayer * poNewLayer ){    papoLayers = (OGRLayer **)        CPLRealloc( papoLayers, sizeof(void*) * ++nLayers );        papoLayers[nLayers-1] = poNewLayer;}/************************************************************************//*                              AddAreaLayer()                              *//************************************************************************/void ILI1Reader::AddAreaLayer( OGRLayer * poAreaLayer,  OGRLayer * poLineLayer ){    ++nAreaLayers;    papoAreaLayers = (OGRLayer **)        CPLRealloc( papoAreaLayers, sizeof(void*) * nAreaLayers );    papoAreaLineLayers = (OGRLayer **)        CPLRealloc( papoAreaLineLayers, sizeof(void*) * nAreaLayers );        papoAreaLayers[nAreaLayers-1] = poAreaLayer;    papoAreaLineLayers[nAreaLayers-1] = poLineLayer;}/************************************************************************//*                              AddSurfaceLayer()                              *//************************************************************************/void ILI1Reader::AddSurfaceLayer( OGRLayer * poDataLayer,  OGRLayer * poPolyLayer ){    ++nSurfaceLayers;    papoSurfaceLayers = (OGRLayer **)        CPLRealloc( papoSurfaceLayers, sizeof(void*) * nSurfaceLayers );    papoSurfacePolyLayers = (OGRLayer **)        CPLRealloc( papoSurfacePolyLayers, sizeof(void*) * nSurfaceLayers );        papoSurfaceLayers[nSurfaceLayers-1] = poDataLayer;    papoSurfacePolyLayers[nSurfaceLayers-1] = poPolyLayer;}/************************************************************************//*                              GetLayer()                              *//************************************************************************/OGRLayer *ILI1Reader::GetLayer( int iLayer ){    if( iLayer < 0 || iLayer >= nLayers )        return NULL;    else        return papoLayers[iLayer];}OGRLayer *ILI1Reader::GetLayerByName( const char* pszLayerName ){    for(int iLayer = 0; iLayer < nLayers; iLayer++ )    {        if( EQUAL(pszLayerName,                  papoLayers[iLayer]->GetLayerDefn()->GetName()) )            return papoLayers[iLayer];    }    return NULL;}/************************************************************************//*                           GetLayerCount()                            *//************************************************************************/int ILI1Reader::GetLayerCount(){    return nLayers;}/************************************************************************//*     Read one logical line, and return split into fields.  The return *//*     result is a stringlist, in the sense of the CSL functions.       *//************************************************************************/char ** ILI1Reader::ReadParseLine(){    const char  *pszLine;    char **tokens;    char **conttok;    char *token;    CPLAssert( fpItf != NULL );    if( fpItf == NULL )        return( NULL );        pszLine = CPLReadLine( fpItf );    if( pszLine == NULL )        return( NULL );        if (strlen(pszLine) == 0) return NULL;          tokens = CSLTokenizeString2( pszLine, " ", CSLT_PRESERVEESCAPES );    token = tokens[CSLCount(tokens)-1];        //Append CONT lines    while (strlen(pszLine) && EQUALN(token, "\\", 2))    {       //remove last token      CPLFree(tokens[CSLCount(tokens)-1]);      tokens[CSLCount(tokens)-1] = NULL;            pszLine = CPLReadLine( fpItf );      conttok = CSLTokenizeString2( pszLine, " ", CSLT_PRESERVEESCAPES );      if (!conttok || !EQUAL(conttok[0], "CONT")) break;            //append      tokens = CSLInsertStrings(tokens, -1, &conttok[1]);      token = tokens[CSLCount(tokens)-1];            CSLDestroy(conttok);    }    return tokens;}IILI1Reader *CreateILI1Reader() {    return new ILI1Reader();}