TerraLib - Development Source: TePDILevelRemap.cpp Source File

00001 #include "TePDILevelRemap.hpp"
00002 
00003 #include "../kernel/TeAgnostic.h"
00004 #include "TePDIUtils.hpp"
00005 #include "TePDIStatistic.hpp"
00006 #include "TePDIMatrix.hpp"
00007 
00008 #include "../kernel/TeUtils.h"
00009 #include "../kernel/TeDefines.h"
00010 
00011 #include <math.h>
00012 #include <float.h>
00013 
00014 TePDILevelRemap::TePDILevelRemap()
00015 {
00016 }
00017 
00018 
00019 TePDILevelRemap::~TePDILevelRemap()
00020 {
00021 }
00022 
00023 
00024 
00025 void TePDILevelRemap::ResetState( const TePDIParameters& params )
00026 {
00027   TePDIParameters dummy_params = params;
00028 
00029   histo_cache_.clear();
00030 }
00031 
00032 
00033 bool TePDILevelRemap::RemapLevels(
00034   TePDITypes::TePDIRasterPtrType& inRaster,
00035   remap_func_3_ptr_type remapp_func,
00036   int in_channel,
00037   int out_channel,
00038   double gain,
00039   double offset,
00040   bool normalize_output,
00041   TePDITypes::TePDIRasterPtrType& outRaster )
00042 {
00043   TEAGN_TRUE_OR_THROW( inRaster.isActive(),
00044     "inRaster inactive" );
00045   TEAGN_TRUE_OR_THROW( outRaster.isActive(),
00046     "outRaster inactive" );
00047   TEAGN_TRUE_OR_THROW( 
00048     ( inRaster->params().status_ != TeRasterParams::TeNotReady ),
00049     "inRaster not ready" );
00050   TEAGN_TRUE_OR_THROW( 
00051     ( outRaster->params().status_ != TeRasterParams::TeNotReady ),
00052     "outRaster not ready" );
00053   TEAGN_CHECK_EQUAL( inRaster->params().nlines_,
00054     outRaster->params().nlines_,
00055     "Lines number mismatch between input and output image" );
00056   TEAGN_CHECK_EQUAL( inRaster->params().ncols_,
00057     outRaster->params().ncols_,
00058     "Columns number mismatch between input and output image" );
00059   TEAGN_TRUE_OR_THROW( in_channel < inRaster->nBands(), "Invalid input band" );
00060   TEAGN_TRUE_OR_THROW( out_channel < outRaster->nBands(), "Invalid output band" );
00061   
00062   if( normalize_output ) {
00063     TEAGN_TRUE_OR_THROW( 
00064       ( outRaster->params().dataType_[ out_channel ] != TeDOUBLE ),
00065       "Cannot apply normalization for TeDOUBLE pixel type output bands" );
00066   }
00067   
00068   TeRaster* inRasterNPtr = inRaster.nakedPointer();
00069   TeRaster* outRasterNPtr = outRaster.nakedPointer();
00070   const int in_raster_lines = inRaster->params().nlines_;
00071   const int in_raster_columns = inRaster->params().ncols_;  
00072   
00073   /* Guessing channel bounds */
00074   
00075   double in_channel_min_level = 0;
00076   double in_channel_max_level = 0;
00077 
00078   TEAGN_TRUE_OR_THROW( TePDIUtils::TeGetRasterMinMaxBounds(
00079     inRaster, in_channel, in_channel_min_level,
00080     in_channel_max_level ), 
00081     "Unable to get raster input channel level bounds" );  
00082 
00083   double out_channel_min_level = 0;
00084   double out_channel_max_level = 0;
00085 
00086   TEAGN_TRUE_OR_THROW( TePDIUtils::TeGetRasterMinMaxBounds(
00087     outRaster, out_channel, out_channel_min_level,
00088     out_channel_max_level ), 
00089     "Unable to get raster output channel level bounds" );
00090     
00091   /* Guessing dummy use */
00092   
00093   bool inRaster_uses_dummy = inRaster->params().useDummy_;
00094   double inRaster_dummy = 0;
00095   if( inRaster_uses_dummy ) {
00096     inRaster_dummy = inRaster->params().dummy_[ in_channel ];
00097   }
00098   
00099   bool outRaster_uses_dummy = outRaster->params().useDummy_;
00100   double outRaster_dummy = 0;
00101   if( outRaster_uses_dummy ) {
00102     outRaster_dummy = outRaster->params().dummy_[ out_channel ];
00103   }
00104   
00105   /* Level remapping */
00106   
00107   if( normalize_output ) {
00108     /* Guessing normalization factors */
00109   
00110     double norm_min_level = DBL_MAX;
00111     double norm_max_level = -1.0 * norm_min_level;
00112     double norm_factor = 0;
00113     double norm_off = 0;
00114     
00115     double current_level = 0;
00116     
00117     TePDIMatrix< double > mapped_matrix;
00118     TEAGN_TRUE_OR_THROW( mapped_matrix.Reset( in_raster_lines, 
00119       in_raster_columns, TePDIMatrix< double >::AutoMemPol ),
00120       "Mapped levels matrix reset error" );
00121     
00122     StartProgInt( "Guessing normalization factors...", in_raster_lines );
00123     
00124     int line = 0;
00125     int column = 0;
00126       
00127     for( line = 0 ; line < in_raster_lines ; ++line ) {
00128       TEAGN_FALSE_OR_RETURN( UpdateProgInt( line ), "Canceled by the user" );
00129     
00130       for( column = 0 ; column < in_raster_columns ; ++column ) {
00131         if( ! inRasterNPtr->getElement( column, line, current_level, 
00132           in_channel ) ) {
00133           
00134           TEAGN_TRUE_OR_RETURN( inRaster_uses_dummy, "Raster read error" );
00135           
00136           current_level = inRaster_dummy;
00137         }
00138         
00139         current_level = remapp_func( current_level, gain, offset ); 
00140           
00141         mapped_matrix( line, column ) = current_level;
00142           
00143         if( current_level < norm_min_level ) {
00144           norm_min_level = current_level;
00145         } 
00146         if( current_level > norm_max_level ) {
00147           norm_max_level = current_level;
00148         }                  
00149       }
00150     }
00151        
00152     StopProgInt();
00153     
00154     norm_factor = ( out_channel_max_level - out_channel_min_level ) /
00155       ( norm_max_level - norm_min_level );
00156     norm_off = out_channel_min_level - norm_min_level;
00157     
00158     /* Output Level remapping */
00159     
00160     double output_level = 0;
00161     double matrix_level = 0;
00162     
00163     StartProgInt( "Remapping Levels...", in_raster_lines );
00164   
00165     for( line = 0 ; line < in_raster_lines ; ++line ) {
00166       TEAGN_FALSE_OR_RETURN( UpdateProgInt( line ), "Canceled by the user" );
00167     
00168       for( column = 0 ; column < in_raster_columns ; ++column ) {
00169         matrix_level = mapped_matrix( line, column );
00170         output_level = ( matrix_level + norm_off ) * norm_factor;
00171         
00172         TEAGN_TRUE_OR_RETURN( outRasterNPtr->setElement( column, line,
00173           output_level, out_channel ),
00174           "Level remmaping error at " + Te2String( line ) +
00175           "," + Te2String( column ) );            
00176       }
00177     }
00178     
00179     StopProgInt();
00180   } else {
00181     /* Level remapping */
00182   
00183     double current_level;
00184     double output_level;
00185     
00186     StartProgInt( "Remapping Levels...", in_raster_lines );
00187   
00188     int line = 0;
00189     int column = 0;
00190     for( line = 0 ; line < in_raster_lines ; ++line ) {
00191       TEAGN_FALSE_OR_RETURN( UpdateProgInt( line ), "Canceled by the user" );
00192     
00193       for( column = 0 ; column < in_raster_columns ; ++column ) {
00194         if( inRasterNPtr->getElement( column, line, current_level,
00195             in_channel ) ) {
00196   
00197           /* Finding mapped level by using the level remapping function */
00198   
00199           output_level = remapp_func( current_level, gain, offset );
00200   
00201           /* Level range filtering */
00202           
00203           if( output_level < out_channel_min_level ) {
00204             output_level = out_channel_min_level;
00205           } else if( output_level > out_channel_max_level ) {
00206             output_level = out_channel_max_level;
00207           }
00208   
00209           TEAGN_TRUE_OR_RETURN( outRasterNPtr->setElement( column, line,
00210             output_level, out_channel ),
00211             "Level remmaping error at " + Te2String( line ) +
00212             "," + Te2String( column ) );
00213         } else {
00214           TEAGN_TRUE_OR_RETURN( inRaster_uses_dummy, "Raster read error" );
00215             
00216           TEAGN_TRUE_OR_RETURN( outRasterNPtr->setElement( column, line,
00217             outRaster_dummy, out_channel ),
00218             "Level remmaping error at " + Te2String( line ) +
00219             "," + Te2String( column ) );          
00220         }
00221       }
00222     }
00223     
00224     StopProgInt();
00225   }
00226   
00227   return true;
00228 }
00229 
00230 
00231 bool TePDILevelRemap::RemapLevels(
00232   TePDITypes::TePDIRasterPtrType& inRaster1,
00233   TePDITypes::TePDIRasterPtrType& inRaster2,
00234   remap_func_4_ptr_type remapp_func,
00235   int in_channel1,
00236   int in_channel2,
00237   int out_channel,
00238   double gain,
00239   double offset,
00240   bool normalize_output,
00241   TePDITypes::TePDIRasterPtrType& outRaster )
00242 {
00243   TEAGN_TRUE_OR_THROW( inRaster1.isActive(),
00244     "inRaster1 inactive" );
00245   TEAGN_TRUE_OR_THROW( inRaster2.isActive(),
00246     "inRaster2 inactive" );
00247   TEAGN_TRUE_OR_THROW( outRaster.isActive(),
00248     "outRaster inactive" );
00249   TEAGN_TRUE_OR_THROW( 
00250     ( inRaster1->params().status_ != TeRasterParams::TeNotReady ),
00251     "inRaster1 not ready" );
00252   TEAGN_TRUE_OR_THROW( 
00253     ( inRaster2->params().status_ != TeRasterParams::TeNotReady ),
00254     "inRaster2 not ready" );
00255   TEAGN_TRUE_OR_THROW( 
00256     ( outRaster->params().status_ != TeRasterParams::TeNotReady ),
00257     "outRaster not ready" );
00258   TEAGN_CHECK_EQUAL( inRaster1->params().nlines_,
00259     outRaster->params().nlines_,
00260     "Lines number mismatch between input1 and output image" );
00261   TEAGN_CHECK_EQUAL( inRaster1->params().ncols_,
00262     outRaster->params().ncols_,
00263     "Columns number mismatch between input1 and output image" );
00264   TEAGN_CHECK_EQUAL( inRaster2->params().nlines_,
00265     outRaster->params().nlines_,
00266     "Lines number mismatch between input2 and output image" );
00267   TEAGN_CHECK_EQUAL( inRaster2->params().ncols_,
00268     outRaster->params().ncols_,
00269     "Columns number mismatch between input2 and output image" );
00270   TEAGN_TRUE_OR_THROW( in_channel1 < inRaster1->nBands(), "Invalid input1 band" );
00271   TEAGN_TRUE_OR_THROW( in_channel2 < inRaster2->nBands(), "Invalid input2 band" );
00272   TEAGN_TRUE_OR_THROW( out_channel < outRaster->nBands(), "Invalid output band" );
00273   
00274   if( normalize_output ) {
00275     TEAGN_TRUE_OR_THROW( 
00276       ( outRaster->params().dataType_[ out_channel ] != TeDOUBLE ),
00277       "Cannot apply normalization for TeDOUBLE pixel type output bands" );
00278   }
00279   
00280   TeRaster* inRaster1NPtr = inRaster1.nakedPointer();
00281   TeRaster* inRaster2NPtr = inRaster2.nakedPointer();
00282   TeRaster* outRasterNPtr = outRaster.nakedPointer();
00283   const int in_raster_lines = inRaster1->params().nlines_;
00284   const int in_raster_columns = inRaster1->params().ncols_;    
00285   
00286   /* Guessing dummy use */
00287   
00288   bool inRaster1_uses_dummy = inRaster1->params().useDummy_;
00289   bool inRaster2_uses_dummy = inRaster2->params().useDummy_;
00290   bool outRaster_uses_dummy = outRaster->params().useDummy_;
00291   double outRaster_dummy = 0;
00292   if( outRaster_uses_dummy ) {
00293     outRaster_dummy = outRaster->params().dummy_[ out_channel ];
00294   }   
00295   
00296   /* Guessing channel bounds */
00297   
00298   double in_channel1_min_level = 0;
00299   double in_channel1_max_level = 0;
00300 
00301   TEAGN_TRUE_OR_THROW( TePDIUtils::TeGetRasterMinMaxBounds(
00302     inRaster1, in_channel1, in_channel1_min_level,
00303     in_channel1_max_level ), 
00304     "Unable to get raster1 input channel level bounds" );  
00305     
00306   double in_channel2_min_level = 0;
00307   double in_channel2_max_level = 0;
00308 
00309   TEAGN_TRUE_OR_THROW( TePDIUtils::TeGetRasterMinMaxBounds(
00310     inRaster2, in_channel2, in_channel2_min_level,
00311     in_channel2_max_level ), 
00312     "Unable to get raster2 input channel level bounds" );  
00313 
00314   double out_channel_min_level = 0;
00315   double out_channel_max_level = 0;
00316 
00317   TEAGN_TRUE_OR_THROW( TePDIUtils::TeGetRasterMinMaxBounds(
00318     outRaster, out_channel, out_channel_min_level,
00319     out_channel_max_level ), 
00320     "Unable to get raster output channel level bounds" );
00321     
00322   /* Level remapping */
00323     
00324   if( normalize_output ) {
00325     /* Guessing normalization factors */
00326     
00327     double norm_factor = 0;
00328     double norm_off = 0;
00329     
00330     TePDIMatrix< double > mapped_matrix;
00331     TEAGN_TRUE_OR_THROW( mapped_matrix.Reset( in_raster_lines, 
00332       in_raster_columns, TePDIMatrix< double >::AutoMemPol ),
00333       "Mapped levels matrix reset error" );    
00334   
00335     StartProgInt( "Guessing normalization factors...", in_raster_lines );
00336   
00337     double norm_min_level = DBL_MAX;
00338     double norm_max_level = -1.0 * norm_min_level;
00339     int line = 0;
00340     int column = 0;
00341     double value1 = 0;
00342     double value2= 0;
00343     bool got_element1 = false;
00344     bool got_element2 = false;
00345     double mapped_level = 0;
00346     
00347     for( line = 0 ; line < in_raster_lines ; ++line ) {
00348       TEAGN_FALSE_OR_RETURN( UpdateProgInt( line ), "Canceled by the user" );
00349   
00350       for( column = 0 ; column < in_raster_columns ; ++column ) {
00351         got_element1 = inRaster1NPtr->getElement( column, line, value1, 
00352           in_channel1 );
00353         got_element2 = inRaster2NPtr->getElement( column, line, value2, 
00354           in_channel2 );
00355           
00356         if( got_element1 && got_element2 ) {
00357           mapped_level = remapp_func( value1, value2 , gain, offset );    
00358         
00359           mapped_matrix( line, column ) = mapped_level;
00360             
00361           if( mapped_level < norm_min_level ) {
00362             norm_min_level = mapped_level;
00363           } 
00364           if( mapped_level > norm_max_level ) {
00365             norm_max_level = mapped_level;
00366           }
00367         } else {
00368           if( got_element1 ) {
00369             TEAGN_TRUE_OR_RETURN( inRaster2_uses_dummy, 
00370               "Raster2 read error" );
00371           } else {
00372             TEAGN_TRUE_OR_RETURN( inRaster1_uses_dummy, 
00373               "Raster1 read error" );
00374           }        
00375         
00376           mapped_matrix( line, column ) = outRaster_dummy;
00377         }
00378       }
00379     }
00380      
00381     StopProgInt();
00382     
00383     norm_factor = ( out_channel_max_level - out_channel_min_level ) /
00384       ( norm_max_level - norm_min_level );
00385     norm_off = out_channel_min_level - norm_min_level;
00386     
00387     /* Remapping levels */
00388     
00389     double output_level = 0;
00390     
00391     StartProgInt( "Remapping Levels...", in_raster_lines );
00392     
00393     for( line = 0 ; line < in_raster_lines ; ++line ) {
00394       TEAGN_FALSE_OR_RETURN( UpdateProgInt( line ), "Canceled by the user" );
00395     
00396       for( column = 0 ; column < in_raster_columns ; ++column ) {
00397         output_level = mapped_matrix( line, column );
00398               
00399         output_level = ( output_level + norm_off ) * norm_factor;
00400         
00401         TEAGN_TRUE_OR_RETURN( outRasterNPtr->setElement( column, line,
00402           output_level, out_channel ),
00403           "Level remmaping error at " + Te2String( line ) +
00404           "," + Te2String( column ) );        
00405       }  
00406     }
00407   } else {
00408     int line = 0;
00409     int column = 0;  
00410     double current_level1 = 0;
00411     double current_level2 = 0;
00412     double output_level = 0;
00413     bool got_element1 = false;
00414     bool got_element2 = false;
00415     
00416     StartProgInt( "Remapping Levels...", in_raster_lines );
00417   
00418     for( line = 0 ; line < in_raster_lines ; ++line ) {
00419       TEAGN_FALSE_OR_RETURN( UpdateProgInt( line ), "Canceled by the user" );
00420     
00421       for( column = 0 ; column < in_raster_columns ; ++column ) {
00422         got_element1 = inRaster1NPtr->getElement( column, line, 
00423           current_level1, in_channel1 );
00424         got_element2 = inRaster2NPtr->getElement( column, line, 
00425           current_level2, in_channel2 );
00426               
00427         if( got_element1 && got_element2 ) {
00428   
00429           /* Finding mapped level by using the level remapping function */
00430   
00431           output_level = remapp_func( current_level1, current_level2 , gain, 
00432             offset );
00433   
00434           /* Level range filtering */
00435           
00436           if( output_level < out_channel_min_level ) {
00437             output_level = out_channel_min_level;
00438           } else if( output_level > out_channel_max_level ) {
00439             output_level = out_channel_max_level;
00440           }
00441 
00442           TEAGN_TRUE_OR_RETURN( outRasterNPtr->setElement( column, line,
00443             output_level, out_channel ),
00444             "Level remmaping error at " + Te2String( line ) +
00445             "," + Te2String( column ) );
00446         } else {
00447           if( got_element1 ) {
00448             TEAGN_TRUE_OR_RETURN( inRaster2_uses_dummy, 
00449               "Raster2 read error" );
00450           } else {
00451             TEAGN_TRUE_OR_RETURN( inRaster1_uses_dummy, 
00452               "Raster1 read error" );
00453           }
00454             
00455           TEAGN_TRUE_OR_RETURN( outRaster->setElement( column, line,
00456             outRaster_dummy, out_channel ),
00457             "Level remmaping error at " + Te2String( line ) +
00458             "," + Te2String( column ) );
00459         }
00460       }  
00461     }
00462   }
00463   
00464   return true;
00465 }
00466 
00467 
00468 void TePDILevelRemap::BuildHistograms(
00469   TePDITypes::TePDIRasterPtrType& inRaster,
00470   unsigned int histo_levels,
00471   std::vector< int >& channels,
00472   bool force )
00473 {
00474   TEAGN_TRUE_OR_THROW( inRaster.isActive(),
00475     "inRaster inactive" );
00476   TEAGN_TRUE_OR_THROW( inRaster->params().status_ != TeRasterParams::TeNotReady,
00477     "inRaster not ready" );
00478 
00479   std::pair< TeRaster*, unsigned int  > mapkey;
00480   mapkey.first = inRaster.nakedPointer();
00481 
00482   for( unsigned int channels_index = 0 ; channels_index < channels.size() ;
00483        ++channels_index ) {
00484 
00485     TEAGN_TRUE_OR_THROW( channels[ channels_index ] < inRaster->nBands(),
00486       "Trying to creat histogram from an invalid band" );
00487 
00488     mapkey.second = channels_index;
00489 
00490     if( ( histo_cache_.find( mapkey ) ==  histo_cache_.end() ) || force ) {
00491 
00492       TePDIHistogram::pointer temp_hist( new TePDIHistogram );
00493 
00494       TEAGN_TRUE_OR_THROW( temp_hist->reset( inRaster,
00495         channels[ channels_index ], histo_levels, 
00496         TeBoxPixelIn ),
00497         "Histogram Generation Error" );
00498 
00499       histo_cache_[ mapkey ] = temp_hist;
00500     }
00501   }
00502 }
00503 
00504 
00505 TePDIHistogram::pointer TePDILevelRemap::GetHistRef(
00506   TePDITypes::TePDIRasterPtrType& inRaster,
00507   unsigned int band )
00508 {
00509   std::pair< TeRaster*, unsigned int  > mapkey;
00510   mapkey.first = inRaster.nakedPointer();
00511   mapkey.second = band;
00512 
00513   TEAGN_TRUE_OR_THROW(
00514     ( histo_cache_.find( mapkey ) != histo_cache_.end() ),
00515     "Histogram not stored inside cache" );
00516 
00517   return histo_cache_[ mapkey ];
00518 }
00519 
00520