Shared Samples

Try out source code samples stored in Bitbucket repositories.
Showing results for 
Search instead for 
Do you mean 

Job imageResample: how to resample an image

by Technical Evangelist on ‎09-08-2015 10:36 AM - edited on ‎03-29-2016 11:30 AM by Anonymous (1,310 Views)

Example Type

Job

 

 

Purpose

 

This example shows how to resample an image

 

Files

 

imageResample.c                    (C source code)

imageResample.eml    (EML script)

imageResample.dsp    (Visual C++ Project)(Windows Only)

 

 

Description

 

The imageResample example performs a map-orienting calibration (transformation), subsetting and decimation simultaneously. This job will use the eimg_LayerStackCreate and the eimg_LayerStackWrite functions to produce the new image. In creating the new resampled transformation, the Exfr_Xform package will be utilized. In summary, the process of creating and new image and modifying the transformation will be demonstrated.

 

 

Command Line Options

 

job geotran -in[put] <inputfilename>  -out[put] <outputfilename> -sub[set] <subset>  -pix[size] <pixelsize>

 

 -in[put]

 The input .img file name.

 

-out[put]

 The output .img file name.

 

-sub[set]

 The upper left x, upper left y,  lower right x,  lower right y (new image size).

 

-pix[size]

 The new pixel size of the image.

 

 

Interface

 demo.jpg

 

 

Output

The output image is obtained by resampling the specified subset of the input image according to the given pixel size.

 

Source Code:

https://bitbucket.org/HGDN/hgdn_erdas_imagine/src/e1ade3b7d148b0f2f1aa6427044dfa4bce6226ee/Job_image...

 

/************************************************************************
*************************************************************************
 Job: imageResample 
 
 Purpose: This example shows how to resample an image.

 Calling syntax: 

 imageresample -in[put] <inputfile> -out[put] <outputfilename> 
 -sub[set] ulx uly lrx lry -pix[size] pixelsize -m <0|1>

 e.g.
 imageresample -in "$IMAGINE_HOME/examples/lanier.img" 
 -out "$IMAGINE_HOME/examples/temp.img" -sub 233085 3807070 248415 3791740
 -pix 10 -m

 IMPORTANT NOTE: This example is going to make a critical assumption about
 the input data. The assumption is that all layers of the input image 
 have the same coordinate extents and pixel size. Multiple layer image
 files are often NOT homogenous and as such beware of the effects
 that this may have on either testing this example or recreating this 
 logic in a custom application.
*************************************************************************
*************************************************************************/

/*
** Toolkit Includes
*/
#include <exfr.h>
#include <emsc.h>
#include <eprj.h>
#include <eimg.h>
#include <eirf.h>
#include <eerr.h>
#include <esmg.h>
#include <eeml.h>
#include <earg.h>
#include <estr.h>
#include <efnp.h>
#include <emet.h>



/* set up macro for error checking */

#define EFS1( msg, arg ) EERR_FATAL_SET_1( __LINE__, msg, arg )
#define EFS0( msg ) EERR_FATAL_SET_0( __LINE__, msg )
#define ES0( e, msg ) EERR_SHOW_0( e, EERR_DEBUG, __LINE__, msg )

/* prototype declaration */

static Exfr_XForm *XformResample(
    Exfr_XForm *,
    long *,
    long *,
    double,
    double,
    double,
    double,
    double,
    double,
    Exfr_XForm **,
    Eerr_ErrorReport **
);
/*
** Global Functions 
*/
static void jobMain(int, char **);
static void set_subset(int, char **);
static void set_pixelsize(int, char **);
static void set_input(int, char **);
static void set_output(int, char **);
static void set_meter(int, char **);


static Emet_MeterInfo *meter = NULL;
static char *input = NULL;
static char *output = NULL;
static int pixSizeX;
static int pixSizeY;
static double ss_ulx;
static double ss_uly;
static double ss_lrx;
static double ss_lry;
static Emsc_Boolean pixFlag = EMSC_FALSE;
static Emsc_Boolean ssFlag = EMSC_FALSE;
static	int	rc = -1;

static Earg_Cmd format[] = 
{
        set_pixelsize,    "-pix[size] %d",     "", EARG_NO_ERR,
        set_subset,   "-sub[set] %f %f %f %f", "", EARG_NO_ERR,
        set_input,    "-in[put] %s",        "", EARG_NO_ERR,
        set_output,   "-out[put] %s",       "", EARG_NO_ERR,
        set_meter,    "-m[eter] [%d]",      "", EARG_NO_ERR,
        jobMain,      "imageresample",      "", EARG_NO_ERR,
        EARG_END
};

/*
** Input file 
*/
static void 
set_input(
    int argc,
    char **argv
)
{ 
    input = estr_Duplicate(argv[1]);
    return;
}


/*
** Output file 
*/
static void 
set_output(
    int argc,
    char **argv
)
{
    output = estr_Duplicate(argv[1]);
    return;
}


/*
** Set the pixel output
*/
static void
set_pixelsize(
    int argc,
    char **argv
)
{
    pixSizeY = pixSizeX = atoi(argv[1]);
    pixFlag = EMSC_TRUE;
    return;
}


/* 
** Set the subset 
*/
static void
set_subset(
    int argc,
    char **argv
)
{
    /* read in new min and max */
    ss_ulx = atof(argv[1]);
    ss_uly = atof(argv[2]);
    ss_lrx = atof(argv[3]);
    ss_lry = atof(argv[4]);
    ssFlag = EMSC_TRUE;   
    return;
}           

/*
**  meterinfo for metering
*/
static void
set_meter(
    int argc, 
    char **argv
)

{
    EERR_INIT("set meter", &lclerr, lclerr);
    
    if (argc == 1)
    {
        meter = emet_MeterInfoCreate("Starting...", "imageResample", NULL, &lclerr);
        EERR_CHECK(lclerr, EFS0("Error creating MeterInfo"));
    }
    
    else
    {
        if (estr_Eqic(argv[1], "0"))
        {
            meter = NULL;
        }
        else if (estr_Eqic(argv[1], "1"))
        {
            meter = emet_MeterInfoCreate("Starting...", "imageResample", NULL, &lclerr);
            EERR_CHECK(lclerr, EFS0("Error creating MeterInfo"));
        }
        else 
        {
            meter = NULL;
        }
    }
    
cleanup:

    ES0(lclerr, "Errors reported by set_meter");
    
    return;
}

int main(int argc, char **argv)
{
    EERR_INIT("imageResample: main", &lclerr, lclerr);
    
    /*
    **  Initialize the ERDAS toolkit
    */
    (void)eint_InitToolkit((Emsc_Opaque**)&lclerr);
    EERR_CHECK(lclerr, EFS0("Error initializing the Toolkit"));
    
    /*
    ** Job Initialization
    */
    esmg_JobInit("imageResample", &lclerr);
    EERR_CHECK(lclerr, EFS0("Error with Job Initalization!"));
    
    /*
    **  Parse the command line, and perform the file copy
    */
    earg_DoArgs(format, argc, argv, &lclerr);
    ES0(lclerr, "Error with getting arguments");
    
    /*
    **  Disconnect from the session manager
    */
    esmg_JobEnd(&lclerr);
    ES0(lclerr, "Error ending job");
    
cleanup:
    ES0(lclerr, "Initialization error");
    
    return (rc);
} 

/*
** Geometric Transformation part of the job
*/      
static void
jobMain(
    int argc,
    char **argv
)
{
    EERR_INIT("Georeference", &lclerr, lclerr);
    
    Eint_InitToolkitData *erdinit = eint_GetInit();
    Eimg_LayerNameStack *inlayernames = NULL;
    Eimg_LayerNameStack *outlayernames = NULL;
    Eimg_LayerStack *inlayerstack = NULL;
    Eimg_LayerStack *outlayerstack = NULL;
    Eimg_PixelRectStack *pixRectStack = NULL;
    Eprj_MapInfo *outmapinfo = NULL;
    Eimg_LayerType layertype;
    Egda_PixelType pixeltype; 
    Eprj_MapProjection *projection = NULL;
    Exfr_XForm *resamplingXForm = NULL;
    Exfr_XForm *xform_old = NULL;
    Exfr_XForm *xform_new = NULL;
    char *outlayername = NULL;
    char *proName = NULL;
    char *units = NULL;
    double ulx;
    double uly;
    double lrx;
    double lry;
    long nlayers;
    long outWidth;
    long outHeight;
    long blockx;
    long blocky;
    int layer;
    int i;
    int j;
    long currentblock = 0;
    double percentdone = 0;
    /*
    ** Initialize metering options
    */
    emet_MeterInfoSet(meter, 0.0, 100.0, &lclerr);
    EERR_SHOW(lclerr, EERR_DEBUG);
    
    emet_MeterInfoPrint(meter, 0.0, &lclerr);
    EERR_SHOW(lclerr, EERR_DEBUG);
    
    /*
    ** Check to see if the input and output filename were received.
    */
    EERR_CHECK(!input, EFS0("You did not use -in[put] option to specify the input file!"));
    EERR_CHECK(!output, EFS0("You did not use -out[put] option to specify the output file!"));
    
    /*
    ** Open correct layer(s)
    */
    inlayernames = eimg_LayerGetNames(input, erdinit, &lclerr);
    EERR_CHECK(lclerr, EFS0("Error getting layer!"));
    
    /*
    ** Get Layer Count
    */
    nlayers = inlayernames->count;
    
    /* 
    ** Open Stack layer
    */
    inlayerstack = eimg_LayerStackOpen(inlayernames, erdinit, &lclerr,
        EIMG_LAYER_OPTION_READONLY, EIMG_LAYER_OPTION_END);
    EERR_CHECK(lclerr, EFS0("Error with stack open!"));
    
    /* 
    ** Get the layer, pixel type and sizes
    */
    layertype = inlayerstack->layers[0]->layerType;
    pixeltype = inlayerstack->layers[0]->pixelType;

    /* 
    ** Read the map information 
    */
    eimg_LayerMapInformationRead(inlayerstack->layers[0], &proName,
        &units, &xform_old, &lclerr);
    EERR_CHECK(lclerr, EFS0("Error Read Map Information!"));
    
    
    /*
    ** Read the input layer projection info for use layer
    */      
    eimg_LayerMapProjectionRead(inlayerstack->layers[0],
        &projection, &lclerr);
    EERR_CHECK(lclerr, EFS0("Error in LayerMapProjectionRead!"));  
    
    /* 
    ** Get the mapinformation 
    */
    outmapinfo = eimg_LayerMapInfoGet(inlayerstack->layers[0], &lclerr);
    EERR_CHECK(lclerr, EFS0("Error Getting MapInfo!"));
    
    if (ssFlag) 
    {
    /* 
    ** Calculate new dimensions for XForm 
    */
        ulx = ss_ulx;
        uly = ss_uly;
        lrx = ss_lrx;
        lry = ss_lry;
    }
    else 
    {
        ulx = outmapinfo->upperLeftCenter->x;
        uly = outmapinfo->upperLeftCenter->y;
        lrx = outmapinfo->lowerRightCenter->x;
        lry = outmapinfo->lowerRightCenter->y;
    }
    
    if (!pixFlag)
    {
        pixSizeX = outmapinfo->pixelSize->width;
        pixSizeY = outmapinfo->pixelSize->height;
    }
    
    emet_MeterInfoChangeTaskName(meter, "Computing output X-Form", &lclerr);
    EERR_SHOW(lclerr, EERR_DEBUG);

    /* 
    ** Call the xform computation function 
    */
    resamplingXForm = XformResample(xform_old, &outWidth, &outHeight, ulx, uly, 
        lrx, lry, (double)pixSizeX, (double)pixSizeY, &xform_new, &lclerr);
    EERR_CHECK(lclerr, EFS0("Unable to compute the transformation matrix"));

    /*
    ** This is required. What putting the transformation on the input 
    ** layerstack allows us to do is to read the input in the output 
    ** data space. 
    */
    eimg_LayerStackChangeOptions( inlayerstack, &lclerr,
		EIMG_LAYER_OPTION_TRANSFORMATION, resamplingXForm,
			EIRF_RC_NEAREST_NEIGHBOR,
		EIMG_LAYER_OPTION_USE_SS_LAYERS,
		EIMG_LAYER_OPTION_END );
    
    /*
    ** Get ready to start working on the layerstack
    */

    blockx = inlayerstack->layers[0]->blockWidth;
    blocky = inlayerstack->layers[0]->blockHeight;

    outlayernames = eimg_LayerNameStackCreate(nlayers, &lclerr);
    EERR_CHECK(lclerr, EFS0( "Error creating output layer!" ));
    
    for(layer = 0; nlayers > layer; layer++) 
	{
        /* Create output layer names */
	emsc_Free(outlayername);
        outlayername = estr_Sprintf(NULL, "Layer_%d", &lclerr, layer+1);
        EERR_CHECK(lclerr, EFS0( "estr_Sprintf failed" ) );
        outlayernames->strings[layer] = efnp_FileNodeCreate(NULL, 
            output, NULL, outlayername, &lclerr);
        EERR_CHECK(lclerr, EFS0( "Error creating node name for layer!" )); 
    }
    
    pixRectStack = eimg_PixelRectStackCreate(nlayers, blockx, blocky, pixeltype, &lclerr);
    EERR_CHECK(lclerr, EFS0("Error creating PixelRectStack matrix"));
    
    outlayerstack = eimg_LayerStackCreate(outlayernames, layertype,
        pixeltype, outWidth, outHeight, erdinit, &lclerr,
        EIMG_LAYER_OPTION_COMPUTE_STATS_ON_CLOSE,
        EIMG_LAYER_OPTION_STATS_IGNORE_ZEROS, 
        EIMG_LAYER_OPTION_PROGRESS_METER_ON_CLOSE,meter,
        EIMG_LAYER_OPTION_END);
    EERR_CHECK(lclerr, EFS0("Error creating stack output !"));
    
    /*
    ** Work on the entire stack at one block by block
    */
    emet_MeterInfoChangeTaskName(meter, "Processing data", &lclerr);
    EERR_SHOW(lclerr, EERR_DEBUG);

    /*
    ** An important point to note here is that we can use this logic
    ** of simply looping through the output blocks, reading from the
    ** input and placing directly in the output because we applied
    ** the transformation to the input stack. Just think about what
    ** would happen there was no transform on the input and you 
    ** attempted to simply halve the pixel size. You should get 
    ** twice as many pixels and therefore twice as many blocks. If   
    ** there was no transformation on the input layer stack, you would   
    ** start reading outside the image bounds when you got half way   
    ** through the blocks in each coordinate direction. The   
    ** transformation allows the input to be resampled to the output.
    */

    for (i = 0; i < outHeight; i += blocky)
    {
        for (j = 0; j < outWidth; j += blockx)
        {
            /*
            **  A simple metering scheme here it to figure out how many blocks there are
            **  and then calculate percentage by the current block being worked on.  
            */            
            percentdone = currentblock * 100.0/((outHeight/blocky) * (outWidth/blockx));
	    
	    	emet_MeterInfoPrint(meter, percentdone, &lclerr);
            EERR_SHOW(lclerr, EERR_DEBUG);
            
            eimg_LayerStackRead(inlayerstack, j, i, blockx, blocky, pixRectStack, &lclerr);
            EERR_CHECK(lclerr, EFS0("Error reading block into pixRectStack!"));
            
            eimg_LayerStackWrite(outlayerstack, j, i,
                blockx, blocky, pixRectStack, &lclerr);
            EERR_CHECK(lclerr, EFS0("Error writing pixRectStack to the output layer stack!"));

            currentblock++;
        }
    }
    
    /*
    ** Write the Map Info to the output layer 
    */      
    eimg_LayerStackMapInformationWrite(outlayerstack, proName, units, xform_new, &lclerr);
    EERR_CHECK(lclerr, EFS0("Error in LayerMapInfoWrite!"));
 
    
    /*
    **  Write the projection info to the output layer if it exists
    */     
    if (projection) 
    {
        eimg_LayerStackMapProjectionWrite(outlayerstack, projection, &lclerr);
        EERR_CHECK(lclerr, EFS0("Error in MapProjection!"));     
    }
    
    emet_MeterInfoPrint(meter, 100, &lclerr);
    EERR_SHOW(lclerr, EERR_DEBUG);

    /*
    ** The outlayerstack used EIMG_LAYER_OPTION_PROGRESS_METER_ON_CLOSE which means
    ** that when the program gets to the  eimg_LayerStackClose(outlayerstack, &lclerr);
    ** function below it will meter the calculation of stats and closing of the layer.
    ** With that in mind we should give the meter something pertinent to say.
    */
    emet_MeterInfoChangeTaskName(meter, "Computing stats and closing image", &lclerr);
    EERR_SHOW(lclerr, EERR_DEBUG);
    
	rc = 0;

cleanup:
    ES0(lclerr, "Errors reported by jobMain");
    
    eimg_LayerStackClose(inlayerstack, &lclerr);
    ES0(lclerr,"eimg_LayerStackClose failed");
    
    eimg_LayerStackClose(outlayerstack, &lclerr);
    ES0(lclerr,"eimg_LayerStackClose failed");

    emet_MeterInfoDelete(meter, &lclerr);
    EERR_DELETE(lclerr);
   
    eimg_LayerNameStackDelete(outlayernames, &lclerr);
    EERR_DELETE(lclerr);
    
    eimg_LayerNameStackDelete(inlayernames, &lclerr);
    EERR_DELETE(lclerr);

    exfr_XFormFree(&xform_old, &lclerr);
    EERR_DELETE(lclerr);
    
    exfr_XFormFree(&xform_new, &lclerr);
    EERR_DELETE(lclerr);
    
    exfr_XFormFree(&resamplingXForm, &lclerr);
    EERR_DELETE(lclerr);

	eimg_PixelRectStackDelete( pixRectStack, &lclerr );
	EERR_DELETE(lclerr);
    
    eprj_MapInfoFree(&outmapinfo);  
    eprj_ProjectionFree(&projection);   
    emsc_Free(input);
    emsc_Free(output);
    emsc_Free(outlayername);
    emsc_Free(proName);
    emsc_Free(units);

    return;   
}  


static Exfr_XForm * 
XformResample(
    Exfr_XForm *xform_old,
    long *outWidth,
    long *outHeight,
    double ulx,
    double uly,
    double lrx,
    double lry,
	double pixW,
	double pixH,
    Exfr_XForm **outLayerMapInformation,
    Eerr_ErrorReport **outerr
)
{
    long  newWidth;
    long  newHeight;
    Exfr_XForm *resampledFileToMap = NULL;
    Exfr_XForm *mapToResampledFile = NULL;
    Exfr_XForm *resampledFileToOriginalFile = NULL;
    Exfr_XForm *result = NULL;
    
    /*Macro for error checking*/
    EERR_INIT("XFormCompute", outerr, lclerr);
    
    
    newWidth = (long)((fabs(ulx - lrx))/ pixW) + 1;
    newHeight = (long)((fabs(uly - lry))/ pixH) + 1;
    
    resampledFileToMap = exfr_XFormCreateUsingCoeffs(
        pixW, 0, ulx, 0, -pixH, uly, &lclerr);
    EERR_CHECK(lclerr, EFS0("Error resampling file to map!")); 
    
    mapToResampledFile = exfr_XFormInvert(resampledFileToMap,
        NULL, &lclerr);
    EERR_CHECK(lclerr, EFS0("Error with map to Resampled file!"));
    
    resampledFileToOriginalFile = exfr_XFormCompose(xform_old,
        resampledFileToMap, NULL, &lclerr);    
    EERR_CHECK(lclerr, EFS0("Error Resampled to Original"));
    
    result = resampledFileToOriginalFile;
    resampledFileToOriginalFile = NULL;
    
    if (outLayerMapInformation)
    {
        *outLayerMapInformation = mapToResampledFile;
        mapToResampledFile = NULL;
    }
    if (outHeight)
        *outHeight = newHeight;
    if (outWidth)
        *outWidth = newWidth;
    
cleanup:
    EERR_SHOW(lclerr, EERR_DEBUG);
    
    exfr_XFormFree(&mapToResampledFile, &lclerr);
    EERR_DELETE(lclerr);
    
    exfr_XFormFree(&resampledFileToMap, &lclerr);
    EERR_DELETE(lclerr);
    
    exfr_XFormFree(&resampledFileToOriginalFile, &lclerr);
    EERR_DELETE(lclerr);

    return result;
}                 

 

Overview
Contributors