GLIMS: GLIMSView Glacier Digitization Software

00001 #include "imageformat.h"
00002 #include <iostream>
00003 #include "TraceLog.h"
00004 
00005 ImageFormat::ImageFormat( ) {
00006   // mParms = new (char*)[15];
00007 }
00008 
00009 ImageFormat::~ImageFormat( ) {
00010   // delete [] mParms;
00011 }
00012 
00013 /*
00014 void ImageFormat::getImageInf( GlimsDataset::ImageInf &imginf )
00015 {
00016 
00017 }
00018 */
00019 
00020 std::vector<ImageFormat::BandInf>& ImageFormat::getBandInfSet( ) {
00021   return mBandSet;
00022 }
00023 
00024 const ImageFormat::BandInf& ImageFormat::getBandInf( int iband ) {
00025   return mBandSet[iband];
00026 }
00027 
00028 std::string ImageFormat::getName( ) const {
00029   return mFName;
00030 }
00031 
00032 std::string ImageFormat::dtypeToString( ImageFormat::DataType dt ) {
00033   std::string dtype;
00034 
00035   switch ( dt ) {
00036     case ImageFormat::NOT_DEFINED:
00037       dtype = "UNSUPPORTED";
00038       break;
00039 
00040     case ImageFormat::INT8:
00041       dtype = "INT8";
00042       break;
00043 
00044     case ImageFormat::UINT8:
00045       dtype = "UINT8";
00046       break;
00047 
00048     case ImageFormat::INT16:
00049       dtype = "INT16";
00050       break;
00051 
00052     case ImageFormat::UINT16:
00053       dtype = "UINT16";
00054       break;
00055 
00056     case ImageFormat::INT32:
00057       dtype = "INT32";
00058       break;
00059 
00060     case ImageFormat::UINT32:
00061       dtype = "UINT32";
00062       break;
00063 
00064     case ImageFormat::INT64:
00065       dtype = "INT64";
00066       break;
00067 
00068     case ImageFormat::UINT64:
00069       dtype = "UINT64";
00070       break;
00071 
00072     case ImageFormat::FLOAT16:
00073       dtype = "FLOAT16";
00074       break;
00075 
00076     case ImageFormat::FLOAT32:
00077       dtype = "FLOAT32";
00078       break;
00079 
00080     case ImageFormat::FLOAT64:
00081       dtype = "FLOAT64";
00082       break;
00083 
00084     default:
00085       dtype = "UNDEF";
00086       break;
00087   }
00088 
00089   return dtype;
00090 }
00091 
00092 int ImageFormat::dtypeSize( ImageFormat::DataType dt ) {
00093   int sz_dtype;
00094 
00095   switch ( dt ) {
00096     case ImageFormat::NOT_DEFINED:
00097       sz_dtype = 0;
00098       break;
00099 
00100     case ImageFormat::INT8:
00101       sz_dtype = 1;
00102       break;
00103 
00104     case ImageFormat::UINT8:
00105       sz_dtype = 1;
00106       break;
00107 
00108     case ImageFormat::INT16:
00109       sz_dtype = 2;
00110       break;
00111 
00112     case ImageFormat::UINT16:
00113       sz_dtype = 2;
00114       break;
00115 
00116     case ImageFormat::INT32:
00117       sz_dtype = 4;
00118       break;
00119 
00120     case ImageFormat::UINT32:
00121       sz_dtype = 4;
00122       break;
00123 
00124     case ImageFormat::INT64:
00125       sz_dtype = 8;
00126       break;
00127 
00128     case ImageFormat::UINT64:
00129       sz_dtype = 8;
00130       break;
00131 
00132     case ImageFormat::FLOAT16:
00133       sz_dtype = 2;
00134       break;
00135 
00136     case ImageFormat::FLOAT32:
00137       sz_dtype = 4;
00138       break;
00139 
00140     case ImageFormat::FLOAT64:
00141       sz_dtype = 8;
00142       break;
00143 
00144     default:
00145       sz_dtype = 0;
00146       break;
00147   }
00148 
00149   return sz_dtype;
00150 }
00151 
00152 bool ImageFormat::isFloat( ImageFormat::DataType dt ) {
00153   switch ( dt ) {
00154     case ImageFormat::NOT_DEFINED:
00155       return false;
00156       break;
00157 
00158     case ImageFormat::INT8:
00159       return false;
00160       break;
00161 
00162     case ImageFormat::UINT8:
00163       return false;
00164       break;
00165 
00166     case ImageFormat::INT16:
00167       return false;
00168       break;
00169 
00170     case ImageFormat::UINT16:
00171       return false;
00172       break;
00173 
00174     case ImageFormat::INT32:
00175       return false;
00176       break;
00177 
00178     case ImageFormat::UINT32:
00179       return false;
00180       break;
00181 
00182     case ImageFormat::INT64:
00183       return false;
00184       break;
00185 
00186     case ImageFormat::UINT64:
00187       return false;
00188       break;
00189 
00190     default:
00191       return true;
00192       break;
00193   }
00194 
00195   return true;
00196 }
00197 
00198 char* ImageFormat::scaleData( char **orig,
00199                               int ow,
00200                               int oh,
00201                               int nw,
00202                               int nh,
00203                               DataType dt,
00204                               ScaleType st,
00205                               bool del ) const {
00206   return scaleData( orig, 0, 0, ow, oh, ow, oh, nw, nh, dt, st, del );
00207 }
00208 
00209 char* ImageFormat::scaleData( char **orig,
00210                               int ox,
00211                               int oy,
00212                               int ow,
00213                               int oh,
00214                               int fw,
00215                               int fh,
00216                               int nw,
00217                               int nh,
00218                               DataType dt,
00219                               ScaleType st,
00220                               bool del ) const {
00221   char *newData = NULL ;
00222   bool gotNew = true;
00223 
00224   //\todo handle SCALE_BC in switch statements to get rid of compile warning
00225   switch ( dt ) {
00226     case INT8:
00227       switch ( st ) {
00228         case SCALE_NN:
00229           newData = (char*)scl_nearestneighbor( (char*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00230           break;
00231 
00232         case SCALE_BI:
00233           newData = (char*)scl_bilinear( (char*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00234           break;
00235         default: ;
00236       }
00237       // newData = (char*)scaleDT( *orig, ox, oy, ow, oh, fw, fh, nw, nh, st );
00238       break;
00239 
00240     case INT16:
00241       switch ( st ) {
00242         case SCALE_NN:
00243           newData = (char*)scl_nearestneighbor( (short*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00244           break;
00245 
00246         case SCALE_BI:
00247           newData = (char*)scl_bilinear( (short*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00248           break;
00249         default: ;
00250       }
00251       break;
00252 
00253     case UINT8:
00254       switch ( st ) {
00255         case SCALE_NN:
00256           newData = (char*)scl_nearestneighbor( (unsigned char*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00257           break;
00258 
00259         case SCALE_BI:
00260           newData = (char*)scl_bilinear( (unsigned char*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00261           break;
00262         default: ;
00263       }
00264       break;
00265 
00266     case UINT16:
00267       switch ( st ) {
00268         case SCALE_NN:
00269           newData = (char*)scl_nearestneighbor( (unsigned short*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00270           break;
00271 
00272         case SCALE_BI:
00273           newData = (char*)scl_bilinear( (unsigned short*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00274           break;
00275         default: ;
00276       }
00277       // newData = (char*)scaleDT( (unsigned short*)*orig, ox, oy, ow, oh, fw, fh, nw, nh, st );
00278       break;
00279 
00280     case FLOAT32:
00281       switch ( st ) {
00282         case SCALE_NN:
00283           newData = (char*)scl_nearestneighbor( (float*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00284           break;
00285 
00286         case SCALE_BI:
00287           newData = (char*)scl_bilinear( (float*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00288           break;
00289         default: ;
00290       }
00291       // newData = (char*)scaleDT( (float*)*orig, ox, oy, ow, oh, fw, fh, nw, nh, st );
00292       break;
00293 
00294     case FLOAT64:
00295       switch ( st ) {
00296         case SCALE_NN:
00297           newData = (char*)scl_nearestneighbor( (double*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00298           break;
00299 
00300         case SCALE_BI:
00301           newData = (char*)scl_bilinear( (double*)*orig, ox, oy, ow, oh, fw, fh, nw, nh );
00302           break;
00303         default: ;
00304       }
00305       // newData = (char*)scaleDT( (double*)*orig, ox, oy, ow, oh, fw, fh, nw, nh, st );
00306       break;
00307 
00308     default:
00309       newData = *orig;
00310       gotNew = false;
00311       break;
00312   }
00313 
00314   if ( gotNew && del )
00315     delete [] *orig;
00316 
00317   *orig = newData;
00318   return newData;
00319 }
00320 
00321 template<class DType> DType*
00322 ImageFormat::scaleDT( DType* orig,
00323                       int ox,
00324                       int oy,
00325                       int ow,
00326                       int oh,
00327                       int fw,
00328                       int fh,
00329                       int nw,
00330                       int nh,
00331                       ScaleType st ) const {
00332   char *newData = NULL;
00333   switch ( st ) {
00334     case SCALE_NN:
00335       //newData = (char*)scl_nearestneighbor( orig, ox, oy, ow, oh, fw, fh, nw, nh );
00336       break;
00337 
00338     case SCALE_BI:
00339       //newData = (char*)scl_bilinear( (char*)orig, ox, oy, ow, oh, fw, fh, nw, nh );
00340       break;
00341 
00342     case SCALE_BC:
00343       // newData = scl_bicubic( orig, ow, oh, nw, nh, 0 );
00344       break;
00345 
00346     default:
00347       return orig;
00348   }
00349 
00350   return(DType*)newData;
00351 }
00352 
00353 template<typename DType> DType*
00354 ImageFormat::scl_nearestneighbor( DType* orig,
00355                                   int ox,
00356                                   int oy,
00357                                   int ow,
00358                                   int , // oh,
00359                                   int fw,
00360                                   int fh,   // fh,
00361                                   int nw,
00362                                   int nh ) const {
00363   if ( nw == 0 ) return NULL;
00364 
00365   DType* newData = new DType[ nw*nh ];
00366   double x1 = 0.0, y1= 0.0, x2=0.0, y2=0.0;
00367   int    ind1, ind2, yind1, yind2;
00368   double scl;
00369 
00370   scl = (double)ow / (double)nw;
00371 
00372   for ( y2=0.0; (int)y2 < nh; y2++ ) {
00373     y1 = y2 * scl + oy;
00374     y1 = ((y1 > 0.0 ) ? (int)(y1+0.5) : (int)(y1-0.5));
00375     yind1 = (int)(y1 * fw);
00376     yind2 = (int)(y2 * nw);
00377 
00378     for ( x2=0.0; (int)x2 < nw; x2++ ) {
00379       x1 = x2 * scl + ox;
00380       x1 = ((x1 > 0.0 ) ? (int)(x1+0.5) : (int)(x1-0.5));
00381 
00382       // COULD POSSIBLY BE TAKEN OUT
00383       /*
00384       if( x1 < 0.0 ) x1 = 0.0;
00385       if( x1 >= fw ) x1 = (double)(fw-1);
00386       if( y1 < 0.0 ) y1 = 0.0; */
00387       if ( y1 >= fh ) y1 = (double)(fh-1);
00388 
00389       ind1 = (int)(yind1 + x1);
00390       ind2 = (int)(yind2 + x2);
00391       newData[ind2] = orig[ind1];
00392     }
00393   }
00394 
00395   return newData;
00396 }
00397 
00398 template<class DType> DType*
00399 ImageFormat::scl_bilinear( DType* orig,
00400                            int ox,
00401                            int oy,
00402                            int ow,
00403                            int oh,
00404                            int fw,
00405                            int , // fh,
00406                            int nw,
00407                            int nh ) const {
00408   if ( nw == 0 || nh == 0 ) return NULL;
00409 
00410   DType   dn1, dn2, dn3, dn4;
00411   double  ix, iy, tx, ty;
00412   double  scl, mx, my;
00413   int     fox, foy, cox, coy;
00414   int     oicx, oicy, oifx, oify;
00415   DType*  newData = new DType[ nw*nh ];
00416 
00417   scl = (double)ow / (double)nw;
00418 
00419   for ( int y=0; y < nh; y++ ) {
00420     iy = y * scl;
00421     for ( int x=0; x < nw; x++ ) {
00422       ix = x * scl;
00423       fox = (int)ix;
00424       foy = (int)iy;
00425       cox = fox + 1;
00426       coy = foy + 1;
00427 
00428 
00429       if ( coy >= oh )
00430         coy = (int)(oh-1);
00431 
00432       oify = (oy+foy) * fw;
00433       oicy = (oy+coy) * fw;
00434       oifx = ox+fox;
00435       oicx = ox+cox;
00436 
00437       dn1 = orig[ oify + oifx ];
00438       dn2 = orig[ oify + oicx ];
00439       dn3 = orig[ oicy + oifx ];
00440       dn4 = orig[ oicy + oicx ];
00441 
00442       tx = ix - (int)ix;
00443       ty = iy - (int)iy;
00444 
00445       mx = 1.0 - tx;
00446       my = 1.0 - ty;
00447       newData[ y*nw+x ] = (DType)( mx*my*dn1 + tx*my*dn2 + mx*ty*dn3 + tx*ty*dn4 );
00448     }
00449   }
00450 
00451   return newData;
00452 }
00453 
00454 unsigned char* ImageFormat::scl_bicubic( unsigned char*,
00455                                          int,
00456                                          int,
00457                                          int,
00458                                          int,
00459                                          int ) const {
00460   return NULL;
00461 }
00462 
00463 void ImageFormat::setCornerPnts( int iband, ImageInf &inf ) {
00464   Node p1, p2, p3, p4;
00465 
00466   p1.x = 0.0;
00467   p1.y = 0.0;
00468   p1.z = 0.0;
00469 
00470   p2.x = mBandSet[iband].width;
00471   p2.y = 0.0;
00472   p2.z = 0.0;
00473 
00474   p3.x = mBandSet[iband].width;
00475   p3.y = mBandSet[iband].height;
00476   p3.z = 0.0;
00477 
00478   p4.x = 0.0;
00479   p4.y = mBandSet[iband].height;
00480   p4.z = 0.0;
00481 
00482   getLL( iband, p1.x, p1.y );
00483   getLL( iband, p2.x, p2.y );
00484   getLL( iband, p3.x, p3.y );
00485   getLL( iband, p4.x, p4.y );
00486 
00487   inf.mP1 = p1;
00488   inf.mP2 = p2;
00489   inf.mP3 = p3;
00490   inf.mP4 = p4;
00491 }
00492 
00493 void ImageFormat::xyToMeter( int iband, double &x, double &y ) {
00494   double angle;
00495   double r;
00496   BandInf &binf = mBandSet[iband];
00497 
00498   if ( iband < 0 || iband >= (int)mBandSet.size() ) return;
00499 
00500   if ( binf.proj == "" )
00501     return;
00502 
00503   if ( binf.rot != 0.0 && (x || y)) {
00504     r = sqrt( x*x + y*y );
00505 
00506     if ( x < 0 )
00507       angle = acos( (x * -1.0) / r );
00508     else {
00509       angle = acos( x / r );
00510     }
00511 
00512     // ADJUST FOR NON IMAGE QUADRANT
00513     if ( x < 0 && y > 0 )
00514       angle = M_PI - angle;
00515     else if ( x < 0 && y < 0 ) {
00516       angle += M_PI ;
00517     }
00518     else if ( x > 0 && y < 0 ) {
00519       angle = M_PI + (M_PI-angle) ;
00520     }
00521 
00522     // APPLY BAND ROTATION
00523     angle += (binf.rot * DEG_TO_RAD) ;
00524 
00525     x = r * ( cos( angle ) ) ;
00526     y = r * ( sin( angle ) ) ;
00527   }
00528 
00529   x = binf.ulx + ( x * binf.res );
00530   y = binf.uly - ( y * binf.res );
00531 }
00532 
00533 void ImageFormat::meterToLL( int iband, double &x, double &y ) {
00534   projXY xy;
00535   projPJ pj;
00536 
00537   if ( iband < 0 || iband >= (int)mBandSet.size() ) return;
00538 
00539   if ( !(pj = createProj( iband )) )
00540     return;
00541 
00542   xy.u = x;
00543   xy.v = y;
00544   xy = pj_inv( xy, pj );
00545   x = xy.u * RAD_TO_DEG;
00546   y = xy.v * RAD_TO_DEG;
00547   pj_free( pj );
00548 }
00549 
00550 void ImageFormat::llToMeter( int iband, double &lon, double &lat ) {
00551 
00552   projXY xy;
00553   projPJ pj;
00554 
00555   if ( iband < 0 || iband >= (int)mBandSet.size() ) return;
00556 
00557   if ( !(pj = createProj( iband )) )
00558     return;
00559 
00560   xy.u = lon * DEG_TO_RAD;
00561   xy.v = lat * DEG_TO_RAD;
00562   xy = pj_fwd( xy, pj );
00563   lon = xy.u;
00564   lat = xy.v;
00565   pj_free( pj );
00566 }
00567 
00568 void ImageFormat::llToXY( int iband, double &lon, double &lat ) {
00569   double x, y;
00570   double r;
00571   double angle;
00572 
00573   if ( iband < 0 || iband >= (int)mBandSet.size() ) return;
00574   llToMeter( iband, lon, lat );
00575   double ulx = mBandSet[iband].ulx;
00576   double uly = mBandSet[iband].uly;
00577 
00578   x = (lon - ulx) / mBandSet[iband].res;
00579   y = uly - lat;
00580 
00581   /*
00582   if( uly > 0 )
00583       y = uly - lat;
00584   else
00585       y = lat - uly;
00586   */
00587 
00588   y /= mBandSet[iband].res;
00589 
00590   // APPLY A ROTATION AROUND [0,0]
00591   r = sqrt( x*x + y*y );
00592   if ( r ) { // divide by zero is bad
00593     if ( x < 0 )
00594       angle = acos( (x * -1.0) / r );
00595     else
00596       angle = acos( x / r );
00597   }
00598   else {
00599     return ;
00600   }
00601 
00602   if ( x < 0 && y > 0 )
00603     angle = M_PI - angle;
00604   else if ( x < 0 && y < 0 )
00605     angle += M_PI;
00606   else if ( x > 0 && y < 0 )
00607     angle = M_PI + (M_PI - angle);
00608 
00609   angle -= (mBandSet[iband].rot * DEG_TO_RAD);
00610 
00611   x = r * ( cos( angle ) );
00612   y = r * ( sin( angle ) );
00613 
00614   lon = x;
00615   lat = y;
00616 }
00617 
00618 void ImageFormat::xyToLL( int iband, double &x, double &y ) {
00619   if ( iband < 0 || iband >= (int)mBandSet.size() ) return;
00620   xyToMeter( iband, x, y );
00621 
00622   meterToLL( iband, x, y );
00623 
00624   if ( x < 0 ) {
00625     x += 360 ;
00626   }
00627 }
00628 
00629 projPJ ImageFormat::createProj( int iband ) {
00630   projPJ pj = NULL;
00631   char   projstr[50];
00632   int i = 0;
00633 
00634   if ( iband < 0 || iband >= (int)mBandSet.size() ) return NULL;
00635 
00636   BandInf &binf = mBandSet[iband];
00637 
00638   if ( binf.proj == "utm" ) {
00639     char zone[50];
00640     sprintf( zone, "zone=%d", binf.UTMZ );
00641 
00642     mParms[0] = "proj=utm";
00643     mParms[1] = &zone[0];
00644     for ( i=0; i < (int)binf.parms.size(); i++ )
00645       mParms[i+2] = (char*)binf.parms[i].c_str();
00646 
00647 
00648     pj = pj_init( i+2, mParms );
00649   } else {
00650     sprintf( projstr, "proj=%s", binf.proj.c_str() );
00651 
00652     mParms[0] = projstr;
00653     for ( i=0; i < (int)binf.parms.size(); i++ )
00654       mParms[i+1] = (char*)binf.parms[i].c_str();
00655 
00656     pj = pj_init( i+1, mParms );
00657 
00658   }
00659 
00660   return pj;
00661 }
00662 
00663 void ImageFormat::getULCorners( ) {
00664 
00665   //double x, y;
00666   for ( int iband = 0 ; iband < (int)mBandSet.size() ; iband++ ) {
00667     BandInf &binf = mBandSet[iband] ;
00668     if ( binf.proj == "" ) continue;
00669 
00670     //x = y = 0 ;
00671     // MUST USE NATIVE COORDINATE SYSTEM
00672     //getLL( iband, x, y, true );
00673 
00674     llToMeter( iband, binf.ulx, binf.uly );
00675   }
00676 }
00677
imageformat.cpp
