Changeset 14774

Show
Ignore:
Timestamp:
06/27/08 21:48:32 (5 months ago)
Author:
warmerdam
Message:

reimplement algorithm a completely different way. Much faster!

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/gdal/alg/gdalproximity.cpp

    r14772 r14774  
    3434CPL_CVSID("$Id: gdalchecksum.cpp 13893 2008-02-28 21:08:37Z rouault $"); 
    3535 
    36 static CPLErr LoadTargetLine( GUInt32 *panTargetMask, int iLine, 
    37                               GDALRasterBandH hSrcBand,  
    38                               int nTargetValues, int *panTargetValues ); 
    39 static CPLErr  
    40 LineProximitySearch( int nXSize, int nYSize, int iLine, int nMaxDist, 
    41                      float *pafLastProximityLine, 
    42                      float *pafThisProximityLine, 
    43                      int iFirstTargetLine, int iFirstTargetLineLoc, 
    44                      int nTargetMaskLines, int nTargetMaskLineWords, 
    45                      GUInt32 *panTargetMaskBuffer ); 
    46  
     36static CPLErr 
     37ProcessProximityLine( GInt32 *panSrcScanline, int *panNearX, int *panNearY,  
     38                      int bForward, int iLine, int nXSize, int nMaxDist, 
     39                      float *pafProximity, 
     40                      int nTargetValues, int *panTargetValues ); 
    4741 
    4842/************************************************************************/ 
    …… 
    170164/*      (the current and last line).                                    */ 
    171165/* -------------------------------------------------------------------- */ 
    172     float *pafThisProximityLine; 
    173     float *pafLastProximityLine; 
    174     float *pafWriteBuffer; 
    175  
    176     pafThisProximityLine = (float *) VSIMalloc(4 * nXSize); 
    177     pafLastProximityLine = (float *) VSIMalloc(4 * nXSize); 
    178     pafWriteBuffer = (float *) VSIMalloc(4 * nXSize); 
    179  
    180     if( pafLastProximityLine == NULL ) 
     166    int   *panNearX, *panNearY; 
     167    float *pafProximity; 
     168    GInt32 *panSrcScanline; 
     169 
     170    pafProximity = (float *) VSIMalloc(4 * nXSize); 
     171    panNearX = (int *) VSIMalloc(sizeof(int) * nXSize); 
     172    panNearY = (int *) VSIMalloc(sizeof(int) * nXSize); 
     173    panSrcScanline = (GInt32 *) VSIMalloc(4 * nXSize); 
     174 
     175    if( pafProximity== NULL  
     176        || panNearX == NULL  
     177        || panNearY == NULL ) 
    181178    { 
    182179        CPLError( CE_Failure, CPLE_OutOfMemory,  
    …… 
    186183    } 
    187184 
     185/* -------------------------------------------------------------------- */ 
     186/*      Loop forward over all the lines.                                */ 
     187/* -------------------------------------------------------------------- */ 
     188    int iLine; 
     189    CPLErr eErr = CE_None; 
     190 
    188191    for( i = 0; i < nXSize; i++ ) 
    189     { 
    190         pafThisProximityLine[i] = -2.0; 
    191         pafLastProximityLine[i] = -2.0; 
    192     } 
    193  
    194 /* -------------------------------------------------------------------- */ 
    195 /*      Allocate the target identification buffer.  We use one bit      */ 
    196 /*      per pixel and allocate a number of lines determined by the      */ 
    197 /*      maximum distance we need to search.                             */ 
    198 /* -------------------------------------------------------------------- */ 
    199     GUInt32 *panTargetMaskBuffer; 
    200     int      nTargetMaskLines; 
    201     int      nTargetMaskLineWords; 
    202      
    203     nTargetMaskLines = MIN(nMaxDist * 2 + 1,nYSize); 
    204     nTargetMaskLineWords = (nXSize+31)/32; 
    205     panTargetMaskBuffer = (GUInt32 *)  
    206         VSIMalloc2(4*nTargetMaskLineWords,nTargetMaskLines); 
    207     
    208     if( panTargetMaskBuffer == NULL ) 
    209     { 
    210         CPLError( CE_Failure, CPLE_OutOfMemory,  
    211                   "Out of memory allocating %g byte buffer.",  
    212                   nTargetMaskLineWords * (double) nTargetMaskLines ); 
    213         return CE_Failure; 
    214     } 
    215  
    216     memset( panTargetMaskBuffer, 0,  
    217             4 * nTargetMaskLineWords * nTargetMaskLines ); 
    218  
    219 /* -------------------------------------------------------------------- */ 
    220 /*      Prefill the target buffer.                                      */ 
    221 /* -------------------------------------------------------------------- */ 
    222     int iLine; 
    223     int iFirstTargetLine = 0; 
    224     int iFirstTargetLineLoc = 0; 
    225     CPLErr eErr; 
    226  
    227     for( iLine = 0; iLine < nTargetMaskLines; iLine++ ) 
    228     { 
    229         eErr =  
    230             LoadTargetLine( panTargetMaskBuffer + iLine * nTargetMaskLineWords,  
    231                             iLine, hSrcBand, nTargetValues, panTargetValues ); 
    232         if( eErr != CE_None ) 
    233         { 
    234             CPLFree( panTargetMaskBuffer ); 
    235             return eErr; 
    236         } 
    237     } 
    238  
    239 /* -------------------------------------------------------------------- */ 
    240 /*      Loop over the lines to process.                                 */ 
    241 /* -------------------------------------------------------------------- */ 
     192        panNearX[i] = panNearY[i] = -1; 
     193 
    242194    for( iLine = 0; eErr == CE_None && iLine < nYSize; iLine++ ) 
    243195    { 
    244         // Switch this/last proximity buffers. 
    245         { 
    246             float *pafTemp = pafThisProximityLine; 
    247             pafThisProximityLine = pafLastProximityLine; 
    248             pafLastProximityLine = pafTemp; 
    249         } 
    250  
    251         // Do we need to load a new target line? 
    252         if( (iFirstTargetLine + nTargetMaskLines) < nYSize 
    253             && (iFirstTargetLine + nTargetMaskLines) < iLine + nMaxDist ) 
    254         { 
    255             eErr =  
    256                 LoadTargetLine( panTargetMaskBuffer  
    257                                 + iFirstTargetLineLoc * nTargetMaskLineWords,  
    258                                 iFirstTargetLine + nTargetMaskLines, 
    259                                 hSrcBand, nTargetValues, panTargetValues ); 
    260  
    261             iFirstTargetLine++; 
    262             iFirstTargetLineLoc = (iFirstTargetLineLoc+1) % nTargetMaskLines; 
    263         } 
    264  
    265         // Perform proximity search for this scanline. 
    266         eErr =  
    267             LineProximitySearch( nXSize, nYSize, iLine, nMaxDist, 
    268                                  pafLastProximityLine, pafThisProximityLine, 
    269                                  iFirstTargetLine, iFirstTargetLineLoc, 
    270                                  nTargetMaskLines, nTargetMaskLineWords, 
    271                                  panTargetMaskBuffer ); 
     196        eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iLine, nXSize, 1,  
     197                             panSrcScanline, nXSize, 1, GDT_Int32, 0, 0 ); 
    272198        if( eErr != CE_None ) 
    273199            break; 
    274200 
    275         // We need to turn the -1.0's into >maxdist value. 
    276201        for( i = 0; i < nXSize; i++ ) 
    277         { 
    278             pafWriteBuffer[i] = pafThisProximityLine[i]; 
    279             if( pafWriteBuffer[i] < 0.0 ) 
    280                 pafWriteBuffer[i] = nMaxDist; 
    281         } 
     202            pafProximity[i] = -1.0; 
     203 
     204        ProcessProximityLine( panSrcScanline, panNearX, panNearY,  
     205                              TRUE, iLine, nXSize, nMaxDist, 
     206                              pafProximity, nTargetValues, panTargetValues ); 
     207 
     208        ProcessProximityLine( panSrcScanline, panNearX, panNearY,  
     209                              FALSE, iLine, nXSize, nMaxDist, 
     210                              pafProximity, nTargetValues, panTargetValues ); 
    282211 
    283212        // Write out results. 
    284213        eErr =  
    285214            GDALRasterIO( hProximityBand, GF_Write, 0, iLine, nXSize, 1,  
    286                           pafWriteBuffer, nXSize, 1, GDT_Float32, 0, 0 ); 
    287  
    288          
     215                          pafProximity, nXSize, 1, GDT_Float32, 0, 0 ); 
     216 
    289217        if( eErr != CE_None ) 
    290218            break; 
    291219 
    292         if( !pfnProgress( (iLine+1) / (double) nYSize, "", pProgressArg ) ) 
     220        if( !pfnProgress( 0.5 * (iLine+1) / (double) nYSize,  
     221                          "", pProgressArg ) ) 
    293222        { 
    294223            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" ); 
    …… 
    298227 
    299228/* -------------------------------------------------------------------- */ 
     229/*      Loop backward over all the lines.                               */ 
     230/* -------------------------------------------------------------------- */ 
     231    for( i = 0; i < nXSize; i++ ) 
     232        panNearX[i] = panNearY[i] = -1; 
     233 
     234    for( iLine = nYSize-1; eErr == CE_None && iLine >= 0; iLine-- ) 
     235    { 
     236        // Read first pass proximity 
     237        eErr =  
     238            GDALRasterIO( hProximityBand, GF_Read, 0, iLine, nXSize, 1,  
     239                          pafProximity, nXSize, 1, GDT_Float32, 0, 0 ); 
     240 
     241        if( eErr != CE_None ) 
     242            break; 
     243 
     244        // Read pixel values. 
     245 
     246        eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iLine, nXSize, 1,  
     247                             panSrcScanline, nXSize, 1, GDT_Int32, 0, 0 ); 
     248        if( eErr != CE_None ) 
     249            break; 
     250 
     251        // Process backwards. 
     252        ProcessProximityLine( panSrcScanline, panNearX, panNearY,  
     253                              FALSE, iLine, nXSize, nMaxDist, 
     254                              pafProximity, nTargetValues, panTargetValues ); 
     255 
     256        // Process backwards. 
     257        ProcessProximityLine( panSrcScanline, panNearX, panNearY,  
     258                              TRUE, iLine, nXSize, nMaxDist, 
     259                              pafProximity, nTargetValues, panTargetValues ); 
     260 
     261        // Write out results. 
     262        eErr =  
     263            GDALRasterIO( hProximityBand, GF_Write, 0, iLine, nXSize, 1,  
     264                          pafProximity, nXSize, 1, GDT_Float32, 0, 0 ); 
     265 
     266        if( eErr != CE_None ) 
     267            break; 
     268 
     269        if( !pfnProgress( 0.5 + 0.5 * (nYSize-iLine) / (double) nYSize,  
     270                          "", pProgressArg ) ) 
     271        { 
     272            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" ); 
     273            eErr = CE_Failure; 
     274        } 
     275    } 
     276 
     277/* -------------------------------------------------------------------- */ 
    300278/*      Cleanup                                                         */ 
    301279/* -------------------------------------------------------------------- */ 
    302     CPLFree( pafThisProximityLine ); 
    303     CPLFree( pafLastProximityLine ); 
    304     CPLFree( panTargetMaskBuffer ); 
    305     CPLFree( pafWriteBuffer ); 
     280    CPLFree( panNearX ); 
     281    CPLFree( panNearY ); 
     282    CPLFree( panSrcScanline ); 
     283    CPLFree( pafProximity ); 
    306284 
    307285    return eErr; 
    …… 
    309287 
    310288/************************************************************************/ 
    311 /*                           LoadTargetLine()                           */ 
     289/*                        ProcessProximityLine()                        */ 
    312290/************************************************************************/ 
    313 static CPLErr LoadTargetLine( GUInt32 *panTargetMask, int iLine, 
    314                               GDALRasterBandH hSrcBand,  
    315                               int nTargetValues, int *panTargetValues ) 
     291                       
     292static CPLErr 
     293ProcessProximityLine( GInt32 *panSrcScanline, int *panNearX, int *panNearY,  
     294                      int bForward, int iLine, int nXSize, int nMaxDist, 
     295                      float *pafProximity, 
     296                      int nTargetValues, int *panTargetValues ) 
    316297 
    317298{ 
    318     int i; 
    319     int nXSize = GDALGetRasterXSize( hSrcBand ); 
    320     GInt32 *panScanline; 
    321     CPLErr eErr; 
    322  
    323     memset( panTargetMask, 0, 4 * ((nXSize+31)/32) ); 
    324     panScanline = (GInt32 *) CPLMalloc(4 * nXSize); 
    325  
    326     eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iLine, nXSize, 1,  
    327                          panScanline, nXSize, 1, GDT_Int32, 0, 0 ); 
    328     for( i = 0; i < nXSize; i++ ) 
     299    int iStart, iEnd, iStep, iPixel; 
     300 
     301    if( bForward ) 
     302    { 
     303        iStart = 0; 
     304        iEnd = nXSize;  
     305        iStep = 1; 
     306    } 
     307    else 
     308    { 
     309        iStart = nXSize-1; 
     310        iEnd = -1;  
     311        iStep = -1; 
     312    } 
     313 
     314    for( iPixel = iStart; iPixel != iEnd; iPixel += iStep ) 
    329315    { 
    330316        int bIsTarget; 
    331317 
    332         if( panTargetValues == NULL ) 
    333             bIsTarget = (panScanline[i] != 0); 
     318/* -------------------------------------------------------------------- */ 
     319/*      Is the current pixel a target pixel?                            */ 
     320/* -------------------------------------------------------------------- */ 
     321        if( nTargetValues == 0 ) 
     322            bIsTarget = (panSrcScanline[iPixel] != 0); 
    334323        else 
    335324        { 
    336             int j; 
    337             bIsTarget = FALSE; 
    338             for( j = 0; j < nTargetValues; j++ ) 
     325            int i; 
     326            for( i = 0; i < nTargetValues; i++ ) 
    339327            { 
    340                 bIsTarget |= (panScanline[i] == panTargetValues[i]); 
     328                if( panSrcScanline[iPixel] == panTargetValues[i] ) 
     329                    bIsTarget = TRUE; 
    341330            } 
    342331        } 
    343332 
    344333        if( bIsTarget ) 
    345             panTargetMask[i>>5] |= (0x1 << (i & 0x1f)); 
    346     } 
    347  
    348     CPLFree( panScanline ); 
    349  
    350     return eErr; 
    351 } 
    352  
    353 #define MASK_LINE_INDEX(iReqLine) \ 
    354     ((iReqLine - iFirstTargetLine + iFirstTargetLineLoc) % nTargetMaskLines) 
    355  
    356 #define MASK_LINE(iReqLine) \ 
    357     (panTargetMaskBuffer + nTargetMaskLineWords * MASK_LINE_INDEX(iReqLine)) 
    358  
    359 #define MASK_BIT(panMaskLine,iPixel) \ 
    360     (panMaskLine[iPixel>>5] & (0x1 << (iPixel & 0x1f))) 
    361  
    362 #define MASK_TEST(iPixel,iLine) MASK_BIT(MASK_LINE(iLine),iPixel) 
    363    
    364 /************************************************************************/ 
    365 /*                        LineProximitySearch()                         */ 
    366 /************************************************************************/ 
    367  
    368 static CPLErr  
    369 LineProximitySearch( int nXSize, int nYSize, int iLine, int nMaxDist, 
    370                      float *pafLastProximityLine, 
    371                      float *pafThisProximityLine, 
    372                      int iFirstTargetLine, int iFirstTargetLineLoc, 
    373                      int nTargetMaskLines, int nTargetMaskLineWords, 
    374                      GUInt32 *panTargetMaskBuffer ) 
    375      
    376 { 
    377     int iPixel; 
    378     GUInt32 *panThisTargetMaskLine = MASK_LINE(iLine); 
    379  
    380     for( iPixel = 0; iPixel < nXSize; iPixel++ ) 
    381     { 
    382         // Is this pixel a target pixel?  Distance zero. 
    383         if( MASK_BIT(panThisTargetMaskLine,iPixel) ) 
    384         { 
    385             pafThisProximityLine[iPixel] = 0.0; 
     334        { 
     335            pafProximity[iPixel] = 0.0; 
     336            panNearX[iPixel] = iPixel; 
     337            panNearY[iPixel] = iLine; 
    386338            continue; 
    387339        } 
    388340 
    389         // Were the left and above pixels outside our last search region? 
    390         // If so, we only need to check one new pixel in the bottom right  
    391         // corner of our search area.  
    392         if( iPixel > 0  
    393             && pafThisProximityLine[iPixel-1] == -1.0 
    394             && iLine > 1  
    395             && pafLastProximityLine[iPixel] == -1.0 ) 
    396         { 
    397             if( iPixel < nXSize - nMaxDist 
    398                 && iLine < nYSize - nMaxDist  
    399                 && MASK_TEST(iPixel+nMaxDist,iLine+nMaxDist) ) 
     341/* -------------------------------------------------------------------- */ 
     342/*      Are we near(er) to the closest target to the above (below)      */ 
     343/*      pixel?                                                          */ 
     344/* -------------------------------------------------------------------- */ 
     345        float fNearDistSq = nXSize * nXSize * 2; 
     346        float fDistSq; 
     347 
     348        if( panNearX[iPixel] != -1 ) 
     349        { 
     350            fDistSq =  
     351                (panNearX[iPixel] - iPixel) * (panNearX[iPixel] - iPixel) 
     352                + (panNearY[iPixel] - iLine) * (panNearY[iPixel] - iLine); 
     353 
     354            if( fDistSq < fNearDistSq ) 
    400355            { 
    401                 pafThisProximityLine[iPixel] =  
    402                     sqrt( nMaxDist * nMaxDist * 2 ); 
     356                fNearDistSq = fDistSq; 
    403357            } 
    404358            else 
    405                 pafThisProximityLine[iPixel] = -1.0; 
    406         } 
    407  
    408         // Add case where we search the whole right column and bottom row.  
    409 #ifdef optimization 
    410 #endif 
    411  
    412         // We will do a full search, in concentric boxes out from our center  
    413         // pixel.  This is the most expensive case.  
    414         float fLeastDistSq = (nMaxDist*2) * (nMaxDist*2); 
    415         int iLeastDist = nMaxDist*2; 
    416         int iDistOut; 
    417  
    418         for( iDistOut = 1;  
    419              iDistOut < nMaxDist && iDistOut < iLeastDist*1.42 /*sqrt(2)*/; 
    420              iDistOut++ ) 
    421         { 
    422             int nEdgeLen = iDistOut*2 + 1; 
    423             int iBaseX, iBaseY, iOff; 
    424             float fDistSq; 
    425  
    426             iBaseX = iPixel - iDistOut; 
    427             iBaseY = iLine - iDistOut; 
    428              
    429             for( iOff = 0; iOff < nEdgeLen; iOff++ ) 
    430359            { 
    431                 if( iBaseX+iOff >= 0 && iBaseX+iOff < nXSize ) 
    432                 { 
    433                     // top edge. 
    434                     if( iBaseY >= 0 ) 
    435                     { 
    436                         if( MASK_TEST(iBaseX+iOff,iBaseY) ) 
    437                         { 
    438                             iLeastDist = MIN(iLeastDist,iDistOut); 
    439                             fDistSq = (iOff - iDistOut) * (iOff - iDistOut)  
    440                                 + iDistOut * iDistOut; 
    441                             if( fDistSq < fLeastDistSq ) 
    442                                 fLeastDistSq = fDistSq; 
    443                         } 
    444                     } 
    445                     // bottom edge 
    446                     if( (iBaseY + nEdgeLen - 1) < nYSize ) 
    447                     { 
    448                         if( MASK_TEST(iBaseX+iOff,iBaseY+nEdgeLen-1) ) 
    449                         { 
    450                             iLeastDist = MIN(iLeastDist,iDistOut); 
    451                             fDistSq = (iOff - iDistOut) * (iOff - iDistOut)  
    452                                 + iDistOut * iDistOut; 
    453                             if( fDistSq < fLeastDistSq ) 
    454                                 fLeastDistSq = fDistSq; 
    455                         } 
    456                     } 
    457                 } 
    458  
    459                 if( iBaseY+iOff >= 0 && iBaseY+iOff < nYSize ) 
    460                 { 
    461                     // left edge 
    462                     if( iBaseX >= 0 ) 
    463                     { 
    464                         if( MASK_TEST(iBaseX,iBaseY+iOff) ) 
    465                         { 
    466                             iLeastDist = MIN(iLeastDist,iDistOut); 
    467                             fDistSq = (iOff - iDistOut) * (iOff - iDistOut)  
    468                                 + iDistOut * iDistOut; 
    469                             if( fDistSq < fLeastDistSq ) 
    470                                 fLeastDistSq = fDistSq; 
    471                         } 
    472                     } 
    473                     // right edge. 
    474                     if( (iBaseX + nEdgeLen - 1) < nXSize ) 
    475                     { 
    476                         if( MASK_TEST(iBaseX+nEdgeLen-1,iBaseY+iOff) ) 
    477                         { 
    478                             iLeastDist = MIN(iLeastDist,iDistOut); 
    479                             fDistSq = (iOff - iDistOut) * (iOff - iDistOut)  
    480                                 + iDistOut * iDistOut; 
    481                             if( fDistSq < fLeastDistSq ) 
    482                                 fLeastDistSq = fDistSq; 
    483                         } 
    484                     } 
    485                 } 
     360                panNearX[iPixel] = -1; 
     361                panNearY[iPixel] = -1; 
    486362            } 
    487         } /* next ring (box) out in search */ 
    488  
    489         if( fLeastDistSq <= nMaxDist * nMaxDist ) 
    490             pafThisProximityLine[iPixel] = sqrt(fLeastDistSq); 
    491         else if( iLeastDist < nMaxDist*2 ) 
    492             pafThisProximityLine[iPixel] = -2.0; // hit in box, but >maxdist 
    493         else 
    494             pafThisProximityLine[iPixel] = -1.0; // no hit in box. 
    495     } /* next pixel on scanline */ 
     363        } 
     364 
     365/* -------------------------------------------------------------------- */ 
     366/*      Are we near(er) to the closest target to the left (right)       */ 
     367/*      pixel?                                                          */ 
     368/* -------------------------------------------------------------------- */ 
     369        int iLast = iPixel-iStep; 
     370 
     371        if( iPixel != iStart && panNearX[iLast] != -1 ) 
     372        { 
     373            fDistSq =  
     374                (panNearX[iLast] - iPixel) * (panNearX[iLast] - iPixel) 
     375                + (panNearY[iLast] - iLine) * (panNearY[iLast] - iLine); 
     376 
     377            if( fDistSq < fNearDistSq ) 
     378            { 
     379                fNearDistSq = fDistSq; 
     380                panNearX[iPixel] = panNearX[iLast]; 
     381                panNearY[iPixel] = panNearY[iLast]; 
     382            } 
     383        } 
     384 
     385/* -------------------------------------------------------------------- */ 
     386/*      Are we near(er) to the closest target to the topright           */ 
     387/*      (bottom left) pixel?                                            */ 
     388/* -------------------------------------------------------------------- */ 
     389        int iTR = iPixel+iStep; 
     390 
     391        if( iTR != iEnd && panNearX[iTR] != -1 ) 
     392        { 
     393            fDistSq =  
     394                (panNearX[iTR] - iPixel) * (panNearX[iTR] - iPixel) 
     395                + (panNearY[iTR] - iLine) * (panNearY[iTR] - iLine); 
     396 
     397            if( fDistSq < fNearDistSq ) 
     398            { 
     399                fNearDistSq = fDistSq; 
     400                panNearX[iPixel] = panNearX[iTR]; 
     401                panNearY[iPixel] = panNearY[iTR]; 
     402            } 
     403        } 
     404 
     405/* -------------------------------------------------------------------- */ 
     406/*      Update our proximity value.                                     */ 
     407/* -------------------------------------------------------------------- */ 
     408        if( panNearX[iPixel] != -1  
     409            && fNearDistSq <= nMaxDist * nMaxDist 
     410            && (pafProximity[iPixel] < 0  
     411                || fNearDistSq < pafProximity[iPixel] * pafProximity[iPixel]) ) 
     412            pafProximity[iPixel] = sqrt(fNearDistSq); 
     413    } 
    496414 
    497415    return CE_None; 
    498416} 
    499