Subversion Repositories public iLand

#include "lightroom.h"

#include "../tools/helper.h"

#include <QtCore>

LightRoom::LightRoom()

{

    m_roomObject = 0;

    m_aggregationMode = 0;

}

/** setup routine.

  sets up datastructures (3d space, hemi grids)

  @param dimx size of space in x direction [m]

  @param dimy size of space in y direction [m]

  @param dimz size of space in z direction [m]

  @param cellsize metric length of cells (used for all 3 dimensions).

  @param hemigridsize size of hemigrid-cells (in degrees).

  @param latitude lat. in degrees.

  @param diffus_frac fraction [0..1] of diffus radiation of global radiation. */

void LightRoom::setup(const double dimx, const double dimy, const double dimz,

                      const double cellsize, const double hemigridsize,

                      const double latitude, const double diffus_frac)

{

    DebugTimer t1("setup of lightroom");

    m_countX = int(dimx / cellsize);

    m_countY = int(dimy / cellsize);

    m_countZ = int(dimz / cellsize);

    m_cellsize = cellsize;

    if (m_countX%2==0) m_countX++; // make uneven

    if (m_countY%2==0) m_countY++;

    QRectF rect(-m_countX/2.*cellsize, -m_countY/2.*cellsize, m_countX*cellsize, m_countY*cellsize);

    qDebug() << "Lightroom size: " << m_countX << "/" << m_countY << "/" << m_countZ << " elements. rect: " << rect;

    m_2dvalues.setup(rect, cellsize);

    m_2dvalues.initialize(0.);

    // setup hemigrids

    SolarRadiation solar;

    solar.setLatidude(latitude); // austria

    solar.setVegetationPeriod(0,367); // no. veg. period

    solar.setDiffusRadFraction(diffus_frac); // 50% diffuse radiation

    // calculate global radiation values

    DebugTimer t2("calculate solar radiation matrix");

    solar.calculateRadMatrix(hemigridsize, m_solarGrid);

    t2.showElapsed();

    m_shadowGrid.setup(hemigridsize); // setup size

    m_solarrad_factor = 1. / m_solarGrid.sum(RAD(45)); // sum of rad. > 45°

    m_centervalue = 0.;

}

double LightRoom::calculateGridAtPoint(const double p_x, const double p_y, const double p_z, bool fillShadowGrid)

{

    if (!m_roomObject)

        return 0;

    // check feasibility

    if (m_roomObject->noHitGuaranteed(p_x, p_y, p_z)) {

        //qDebug()<<"skipped";

        return 0;

    }

    if (fillShadowGrid)

        m_shadowGrid.clear(0.);

    // start with 45°

    int ie = m_shadowGrid.indexElevation(RAD(45));

    int ia = 0;

    int max_a = m_shadowGrid.matrixCountAzimuth();

    int max_e = m_shadowGrid.matrixCountElevation();

    double elevation, azimuth;

    double solar_sum=0;

    int hit;

    int c_hit = 0;

    int c_test = 0;

    for (;ie<max_e;ie++){

        for (ia=0;ia<max_a;ia++) {

            azimuth = m_shadowGrid.azimuthNorth(ia);

            elevation = m_shadowGrid.elevation(ie);

            hit = m_roomObject->hittest(p_x, p_y, p_z,azimuth,elevation);

            // if inside the crown: do nothing and return.

            // 20090708: if inside the crown: return "totally dark"

            if (hit==-1)

                return 1;

            //qDebug() << "testing azimuth" << GRAD(azimuth) << "elevation" << GRAD(elevation)<<"hit"<<hit;

            c_test++;

            if (hit==1) {

                // retrieve value from solar grid

                // Sum(cells) of solargrid =1 -> the sum of all "shadowed" pixels therefore is already the "ratio" of shaded/total radiation

                solar_sum += m_solarGrid.rGetByIndex(ia, ie);

                if (fillShadowGrid)

                  m_shadowGrid.rGetByIndex(ia, ie) = m_solarGrid.rGetByIndex(ia, ie);

                c_hit++;

            }

        }

    }

    // solar-rad-factor = 1/(sum rad > 45°)

    return solar_sum * m_solarrad_factor;

    //double ratio = c_hit / double(c_test);

    //qDebug() << "tested"<< c_test<<"hit count:" << c_hit<<"ratio"<<c_hit/double(c_test)<<"total sum"<<m_shadowGrid.getSum();

    //return ratio; // TODO: the global radiation is not calculated!!!!!

}

void LightRoom::calculateFullGrid()

{

    float *v = m_2dvalues.begin();

    float *vend = m_2dvalues.end();

    int z;

    QPoint pindex;

    QPointF coord;

    double hit_ratio;

    DebugTimer t("calculate full grid");

    int c=0;

    float maxh = m_roomObject->maxHeight();

    float *values = new float[m_countZ];

    //double h_realized;

    double sum;

    while (v!=vend) {

        pindex = m_2dvalues.indexOf(v);

        coord = m_2dvalues.cellCenterPoint(pindex);

        double hor_distance = sqrt(coord.x()*coord.x() + coord.y()*coord.y());

        for (z=0;z<m_countZ && z*m_cellsize <= maxh;z++) {

            // only calculate values up to the 45° line

            // tan(45°)=1 -> so this is the case when the distance p->tree > height of the tree

            if (hor_distance > maxh-z*m_cellsize)

                break;

            hit_ratio = calculateGridAtPoint(coord.x(), coord.y(), // coords x,y

                                             z*m_cellsize,false); // heigth (z), false: do not clear and fill shadow grid structure

            // stop calculating when return = -1 (e.g. when inside the crown)

            if (hit_ratio==-1)

                break;

            values[z]=hit_ratio; // this value * height of cells

        }

        // calculate average

        sum = 0;

        // 20090708: do not average, but keep the sum

        // aggregate mean for all cells with angles>45° to the tree-top!

        // 20090711: again average it, but now include the tree top.

        // 20090713: go back to sum...

        // 20090812: remove the weighting for the 10m cells again.

        // 20100520: again, use the average value (per meter)

        for(int i=0;i<z;i++)

            sum+=values[i];

        if (m_aggregationMode==0) {

            // average shadow / meter height (within the 45° cone)

            if (z)

                sum/=float(z);

        } else {

            // aggregation: sum

            sum *= m_cellsize;  // multiply with height (shadow * meter)

        }

        *v = sum; // store in matrix

        // sum

        v++;

        c++;

        if (c%1000==0) {

            qDebug() << c << "processed...time: ms: " << t.elapsed();

            QCoreApplication::processEvents();

        }

        // save value of the middle column...

        m_centervalue = m_2dvalues(0.f, 0.f);

    } // while(v!=vend)

}

//////////////////////////////////////////////////////////

// Lightroom Object

//////////////////////////////////////////////////////////

LightRoomObject::~LightRoomObject()

{

    if (m_radiusFormula)

        delete m_radiusFormula;

}

void LightRoomObject::setuptree(const double height, const double crownheight, const QString &formula)

{

    if (m_radiusFormula)

        delete m_radiusFormula;

    m_radiusFormula = new Expression(formula);

    m_baseradius = m_radiusFormula->calculate(crownheight/height);

    m_height = height;

    m_crownheight = crownheight;

    double h=0., r;

    m_treeHeights.clear();

    // preprocess crown radii for each meter step

    while (h<=m_height) {

        if (h<m_crownheight)

            r=0.;

        else

            r = m_radiusFormula->calculate(h/m_height);

        m_treeHeights.push_back(r);

        h++;

    }

}

// Angle-function. return the  difference between

// two angles as a radian value between -pi..pi [-180..180°]

// >0: directionB is "left" (ccw) of directionA, <0: "right", clockwise

// e.g.: result=-10: 10° cw, result=10°: 10° ccw

// result:-180/+180: antiparallel.

double DiffOfAngles(double DirectionA, double DirectionB)

{

    DirectionA = fmod(DirectionA, PI2); // -> -2pi .. 2pi

    if (DirectionA < 0) DirectionA+=PI2; // -> 0..2pi (e.g.: AngleA = -30 -> 330)

    DirectionB = fmod(DirectionB, PI2);

    if (DirectionB < 0) DirectionB+=PI2;

    double Delta = DirectionB - DirectionA;

    if (Delta<0) {

      if (Delta<-M_PI)

         return Delta  + PI2; // ccw

      else

         return Delta; // cw

    } else {

      if (Delta>M_PI)

         return Delta - PI2; // cw

      else

         return Delta;  // ccw

    }

}

// Angle-function. return the absolute difference between

// two angles as a radian value between 0..pi [0..180°]

// e.g. 10° = +10° or -10°; maximum value is 180°

double AbsDiffOfAngles(double AngleA, double AngleB)

{

     return fabs(DiffOfAngles(AngleA, AngleB));

}

/** The tree is located in x/y=0/0.

*/

int LightRoomObject::hittest(const double p_x, const double p_y, const double p_z,

                              const double azimuth_rad, const double elevation_rad)

{

    bool inside_crown=false;

    if (p_z > m_height)

        return 0;

    // Test 1: does the ray (azimuth) direction hit the crown?

    double phi = atan2(-p_y, -p_x); // angle between P and the tree center

    double dist2d = sqrt(p_x*p_x + p_y*p_y); // distance at ground

    //if (dist2d==0)

    //    return true;

    double alpha = DiffOfAngles(phi, azimuth_rad); //  phi - azimuth_rad; // angle between the ray and the center of the tree

    if (dist2d>m_baseradius) { // test only, if p not the crown

        double half_max_angle = asin(m_baseradius / dist2d); // maximum deviation angle from direct connection p - tree where ray hits maxradius

        if (fabs(alpha) > half_max_angle)

            return 0;

    } else {

        inside_crown = true;

        // test if p is inside the crown

        if (p_z<=m_height && p_z>=m_crownheight) {

            double radius_hit = m_radiusFormula->calculate(p_z / m_height);

            if (dist2d <= radius_hit)

                return -1;

        }

    }

    // Test 2: test if the crown-"plate" at the bottom of the crown is hit.

    if (elevation_rad>0. && p_z<m_crownheight) {

        // calc. base distance between p and the point where the height of the ray reaches the bottom of the crown:

        double r_hitbottom = dist2d; // for 90°

        if (elevation_rad < M_PI_2) {

            double d_hitbottom = (m_crownheight - p_z) / tan(elevation_rad);

            // calc. position (projected) of that hit point

            double rx = p_x + cos(azimuth_rad)*d_hitbottom;

            double ry = p_y + sin(azimuth_rad)*d_hitbottom;

            r_hitbottom = sqrt(rx*rx + ry*ry);

        }

        if (r_hitbottom <= m_baseradius)

            return 1;

    }

    // Test 3: test for height steps...

    // distance from p to the plane normal to p-vector through the center of the tree

    // do only when p is

    double rx,ry,rhit;

    double d_center = cos(alpha)*dist2d;

    if (d_center>=0) {

        double h_center = p_z + d_center*tan(elevation_rad);

        if (h_center<=m_height && h_center>=m_crownheight) {

            rx = p_x + cos(azimuth_rad)*d_center;

            ry = p_y + sin(azimuth_rad)*d_center;

            rhit = sqrt(rx*rx + ry*ry);

            double r_h = m_radiusFormula->calculate(h_center / m_height);

            if (rhit < r_h)

                return 1;

        }

    }

    // Test 4: "walk" through crown using 1m steps.

    double h=floor(p_z);

    double d_1m = 1 / tan(elevation_rad); //projected ground distance equivalent of 1m height difference

    double d_cur = 0;

    if (h!=p_z) {

       d_cur += ((h+1)-p_z)*d_1m;

    }

    while (h<=m_height) {

        double r_tree = m_treeHeights[int(h)];

        rx = p_x + cos(azimuth_rad)*d_cur;

        ry = p_y + sin(azimuth_rad)*d_cur;

        rhit = rx*rx + ry*ry;

        // hit if inside of radius

        if (rhit < r_tree*r_tree)

            return 1;

        // no hit: if formerly was inside crown and now left

        if (inside_crown && rhit > m_baseradius*m_baseradius)

            return 0;

        // enter crown

        if (!inside_crown && rhit <=  m_baseradius*m_baseradius)

            inside_crown = true;

        d_cur+=d_1m;

        h++;

    }

    return 0;

}

bool LightRoomObject::noHitGuaranteed(const double p_x, const double p_y, const double p_z)

{

    // 1. simple: compare height...

    if (p_z > m_height)

        return true;

    // 2. 45° test:

    if (p_z > m_height - sqrt(p_x*p_x + p_y*p_y)) // 45°: height = distance from tree center

        return true;

    return false;

}