Geografic-Information-System/imagery/i.atcorr/Transform.cpp

#include <math.h>
#include "Transform.h"

void EtmDN(int iwave, float asol, bool before, float &lmin, float &lmax)
{
    if (before)		/* ETM+ digital numbers taken before July 1, 2000 */
    {
	switch(iwave)
	{
	case 61:
	{
	    lmin = -6.2f;
	    lmax = 194.3f;
	    break;
	}

	case 62:
	{
	    lmin = -6.0f;
	    lmax = 202.4f;
	    break;
	}

	case 63:
	{
	    lmin = -4.5f;
	    lmax = 158.6f;
	    break;
	}

	case 64:
	{
	    if (asol < 45.)
	    {
		lmin = -4.5f;
		lmax = 235.0f;
	    }
	    else
	    {
		lmin = -4.5f;
		lmax = 157.5f;
	    }
	    break;
	}

	case 65:
	{
	    lmin = -1.0f;
	    lmax = 31.76f;
	    break;
	}

	case 66:
	{
	    lmin = -0.35f;
	    lmax = 10.932f;
	    break;
	}

	case 67:
	{
	    lmin = -5.0f;
	    lmax = 244.00f;
	    break;
	}
	}
    }
    else		/* ETM+ digital numbers taken after July 1, 2000 */
    {
	switch(iwave)
	{
	case 61:
	{
	    lmin = -6.2f;
	    lmax = 191.6f;
	    break;
	}

	case 62:
	{
	    lmin = -6.4f;
	    lmax = 196.5f;
	    break;
	}

	case 63:
	{
	    lmin = -5.0f;
	    lmax = 152.9f;
	    break;
	}

	case 64:
	{
	    if (asol < 45.)
	    {
		lmin = -5.1f;
		lmax = 241.1f;
	    }
	    else
	    {
		lmin = -5.1f;
		lmax = 157.4f;
	    }
	    break;
	}

	case 65:
	{
	    lmin = -1.0f;
	    lmax = 31.06f;
	    break;
	}

	case 66:
	{
	    lmin = -0.35f;
	    lmax = 10.80f;
	    break;
	}

	case 67:
	{
	    lmin = -4.7f;
	    lmax = 243.1f;
	    break;
	}
	}
    }
}

/* Assuming input value between 0 and 1
   if rad is true, idn should first be converted to a reflectance value
   returns adjusted value also between 0 and 1 */
float transform(const TransformInput ti, InputMask imask, float idn)
{
    /* convert from radiance to reflectance */
    if((imask & ETM_BEFORE) || (imask & ETM_AFTER))
    {
        /* http://ltpwww.gsfc.nas */
        float lmin, lmax;
        EtmDN(ti.iwave, ti.asol, imask & ETM_BEFORE, lmin, lmax);

        /* multiply idn by 255.f to correct precondition that idn lies in [0, 255] */
        idn = (lmax - lmin) / 254.f * (idn * 255.f - 1.f) + lmin;
        if (idn < 0.f) idn = 0.f;
        idn /= 255.f;
    }
    if(imask & RADIANCE) idn += (float)M_PI * idn * 255.f * ti.sb / ti.xmus / ti.seb;
          
    float rapp = idn;
    float ainrpix = ti.ainr[0][0];
    float xa = 0.0f;
    float xb = 0.0f;
    float xc = 0.0f;
    float rog = rapp / ti.tgasm;
    /* The if below was added to avoid ground reflectances lower than
       zero when ainr(1,1) greater than rapp/tgasm
       In such case either the choice of atmospheric model was not
       adequate for that image or the calculated apparent reflectance
       was too low. Run the model again for other conditions.
       The lines below just decrease ainr(1,1)/tgasm to avoid too
       bright pixels in the image. Check the output file to see if that
       has happened. */

    float decrfact = 1.0f;
    if (rog < (ainrpix / ti.tgasm))
    {
	do
	{
	    decrfact = decrfact - 0.1f;
	    ainrpix = decrfact * ainrpix;
	}
	while(rog < (ainrpix / ti.tgasm));
    }

    rog = (rog - ainrpix / ti.tgasm) / ti.sutott / ti.sdtott;
    rog = rog / (1.f + rog * ti.sast);
    xa = (float)M_PI * ti.sb / ti.xmus / ti.seb / ti.tgasm / ti.sutott / ti.sdtott;
    xb = ti.srotot / ti.sutott / ti.sdtott / ti.tgasm;
    xc = ti.sast;

    if (rog > 1) rog = 1;
    if (rog < 0) rog = 0;

    return rog;
}