Subversion Repositories public iLand

/********************************************************************************************

**    iLand - an individual based forest landscape and disturbance model

**    http://iland.boku.ac.at

**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl

**

**    This program is free software: you can redistribute it and/or modify

**    it under the terms of the GNU General Public License as published by

**    the Free Software Foundation, either version 3 of the License, or

**    (at your option) any later version.

**

**    This program is distributed in the hope that it will be useful,

**    but WITHOUT ANY WARRANTY; without even the implied warranty of

**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

**    GNU General Public License for more details.

**

**    You should have received a copy of the GNU General Public License

**    along with this program.  If not, see <http://www.gnu.org/licenses/>.

********************************************************************************************/

#include "stampcontainer.h"

#include "globalsettings.h"

#include "exception.h"

//constants

const int StampContainer::cBHDclassWidth=4;

const int StampContainer::cBHDclassLow = 4; ///< bhd classes start with 4cm

const int StampContainer::cBHDclassCount = 70; ///< class count, see getKey(): for lower dbhs classes are smaller

const int StampContainer::cHDclassWidth=10;

const int StampContainer::cHDclassLow = 35; ///< hd classes offset is 35: class 0 = 35-45, class 1 = 45-55

const int StampContainer::cHDclassCount = 16; ///< class count. highest class:  185-195

// static values

Grid<float> StampContainer::m_distance;

StampContainer::StampContainer()

{

    //

    m_lookup.setup(1., // cellsize

                   cBHDclassCount, // count x

                   cHDclassCount); // count y

    m_lookup.initialize(NULL);

    //qDebug() << "grid after init" << gridToString(m_lookup);

    m_maxBhd = -1;

    m_useLookup = true;

}

StampContainer::~StampContainer()

{

    // delete stamps.

    while (!m_stamps.isEmpty())

        delete m_stamps.takeLast().stamp;

}

/// getKey: decodes a floating point piar of dbh and hd-ratio to indices for the

/// lookup table containing pointers to the actual stamps.

inline void StampContainer::getKey(const float dbh, const float hd_value, int &dbh_class, int &hd_class) const

{

    hd_class = int(hd_value - cHDclassLow) / cHDclassWidth;

    // dbh_class = int(dbh - cBHDclassLow) / cBHDclassWidth;

    // fixed scheme: smallest classification scheme for tree-diameters:

    // 1cm width from 4 up to 9cm,

    // 2cm bins from 10 to 18cm

    // 4cm bins starting from 20cm, max DBH=255 (with 70 classes)

    if (dbh < 10.f) {

        dbh_class = qMax(0, int(dbh-4.)); // classes from 0..5

    } else if (dbh<20.f) {

        dbh_class = 6 + int((dbh-10.f) / 2.f); // 10-12cm has index 6

    } else {

        dbh_class = 11 + int((dbh-20.f) / 4.f); // 20-24cm has index 11

    }

}

/** fill up the NULLs in the lookup map */

void StampContainer::finalizeSetup()

{

    if (!m_useLookup)

        return;

    Stamp *s;

    int h;

    int max_size=0;

    for (int b=0;b<cBHDclassCount;b++) {

        // find lowest value...

        for (h=0;h<cHDclassCount;h++) {

            s=m_lookup.valueAtIndex(b,h);

            if (s) {

                // fill up values left from this value

                for (int hfill=0;hfill<h;hfill++)

                    m_lookup.valueAtIndex(b,hfill) = s;

                break;

            }

        }

        // go to last filled cell...

        for (;h<cHDclassCount;h++) {

            if (m_lookup.valueAtIndex(b,h)==0)

                break;

            s=m_lookup.valueAtIndex(b,h);

        }

        // fill up the rest...

        for (;h<cHDclassCount;h++) {

            m_lookup.valueAtIndex(b,h)=s;

        }

        if(s)

            max_size = std::max(max_size, s->dataSize());

        // if no stamps in this dbh-class, copy values (from last row)

        if (!s && b>0) {

           for (h=0;h<cHDclassCount;h++)

               m_lookup.valueAtIndex(b,h) = m_lookup(b-1, h);

        }

    }

    if (!m_lookup.valueAtIndex(0,0)) {

        // first values are missing

        int b=0;

        while (b<cBHDclassCount && m_lookup.valueAtIndex(b,0)==NULL)

            b++;

        for (int fill=0;fill<b;fill++)

            for (h=0;h<cHDclassCount;h++)

                m_lookup.valueAtIndex(fill, h) = m_lookup.valueAtIndex(b,h);

    }

    // distance grid

    if (m_distance.sizeX()<max_size) {

        setupDistanceGrid(max_size);

    }

}

void StampContainer::setupDistanceGrid(const int size)

{

    const float px_size = cPxSize;

    m_distance.setup(px_size, size, size);

    float *p=m_distance.begin();

    QPoint idx;

    for (;p!=m_distance.end();++p) {

        idx = m_distance.indexOf(p);

        *p = sqrt(double(idx.x()*idx.x()) + double(idx.y()*idx.y()))*px_size;

    }

}

 void StampContainer::addStamp(Stamp* stamp, const int cls_dbh, const int cls_hd, const float crown_radius_m, const float dbh, const float hd_value)

 {

    if (m_useLookup) {

        if (cls_dbh<0 || cls_dbh>=cBHDclassCount || cls_hd<0 || cls_hd>=cHDclassCount)

            throw IException(QString("StampContainer::addStamp: Stamp out of range. dbh=%1 hd=%2.").arg(dbh).arg(hd_value));

        m_lookup.valueAtIndex(cls_dbh, cls_hd) = stamp; // save address in look up table

    } // if (useLookup)

    stamp->setCrownRadius(crown_radius_m);

    StampItem si;

    si.dbh = dbh;

    si.hd = hd_value;

    si.crown_radius = crown_radius_m;

    si.stamp = stamp;

    m_stamps.append(si); // store entry in list of stamps

 }

/** add a stamp to the internal storage.

    After loading the function finalizeSetup() must be called to ensure that gaps in the matrix get filled. */

void  StampContainer::addStamp(Stamp* stamp, const float dbh, const float hd_value, const float crown_radius)

{

    int cls_dbh, cls_hd;

    getKey(dbh, hd_value, cls_dbh, cls_hd); // decode dbh/hd-value

    addStamp(stamp, cls_dbh, cls_hd, crown_radius, dbh, hd_value); // dont set crownradius

}

void StampContainer::addReaderStamp(Stamp *stamp, const float crown_radius_m)

{

    double rest = fmod(crown_radius_m, 1.f)+0.0001;

    int cls_hd = int( rest * 10 ); // 0 .. 9.99999999

    if (cls_hd>=cHDclassCount)

        cls_hd=cHDclassCount-1;

    int cls_dbh = int(crown_radius_m);

    //qDebug() << "Readerstamp r="<< crown_radius_m<<" index dbh hd:" << cls_dbh << cls_hd;

    stamp->setCrownRadius(crown_radius_m);

    // prepare special keys for reader stamps

    addStamp(stamp,cls_dbh, cls_hd, crown_radius_m, 0., 0.); // set crownradius, but not dbh/hd

}

/** retrieve a read-out-stamp. Readers depend solely on a crown radius.

Internally, readers are stored in the same lookup-table, but using a encoding/decoding trick.*/

const Stamp* StampContainer::readerStamp(const float crown_radius_m) const

{

    // Readers: from 0..10m in 50 steps???

    int cls_hd = int( (fmod(crown_radius_m, 1.f)+0.0001) * 10 ); // 0 .. 9.99999999

    if (cls_hd>=cHDclassCount)

        cls_hd=cHDclassCount-1;

    int cls_bhd = int(crown_radius_m);

    const Stamp* stamp = m_lookup(cls_bhd, cls_hd);

    if (!stamp)

        qDebug() << "Stamp::readerStamp(): no stamp found for radius" << crown_radius_m;

    return stamp;

}

/** fast access for an individual stamp using a lookup table.

    the dimensions of the lookup table are defined by class-constants.

    If stamp is not found there, the more complete list of stamps is searched. */

const Stamp* StampContainer::stamp(const float bhd_cm, const float height_m) const

{

    float hd_value = 100.f * height_m / bhd_cm;

    int cls_dbh, cls_hd;

    getKey(bhd_cm, hd_value, cls_dbh, cls_hd);

    // check loopup table

    if (cls_dbh<cBHDclassCount && cls_dbh>=0 && cls_hd < cHDclassCount && cls_hd>=0) {

        const Stamp* stamp = m_lookup(cls_dbh, cls_hd);

        if (stamp)

            return stamp;

        if (logLevelDebug())

            qDebug() << "StampContainer::stamp(): not in list: dbh height:" << bhd_cm << height_m << "in" << m_fileName;

    }

    // extra work: search in list...

    // look for a stamp if the HD-ratio is out of range

    if (cls_dbh<cBHDclassCount && cls_dbh>=0) {

        if (logLevelDebug())

            qDebug() << "HD for stamp out of range dbh " << bhd_cm << "and h="<< height_m << "(using smallest/largeset HD)";

        if (cls_hd>=cHDclassCount)

            return m_lookup(cls_dbh, cHDclassCount-1); // highest

        return m_lookup(cls_dbh, 0); // smallest

    }

    // look for a stamp if the DBH is out of range.

    if (cls_hd<cHDclassCount && cls_hd>=0) {

        if (logLevelDebug())

            qDebug() << "DBH for stamp out of range dbh " << bhd_cm << "and h="<< height_m << "-> using largest available DBH.";

        if (cls_dbh>=cBHDclassCount)

            return m_lookup(cBHDclassCount-1, cls_hd); // highest

        return m_lookup(0, cls_hd); // smallest

    }

    // handle the case DBH and HD are out of range

    if (cls_dbh>=cBHDclassCount && cls_hd<0) {

        if (logLevelDebug())

            qDebug() << "DBH AND HD for stamp out of range dbh " << bhd_cm << "and h="<< height_m << "-> using largest available DBH/smallest HD.";

        return m_lookup(cBHDclassCount-1, 0);

    }

    // handle the case that DBH is too high and HD is too high (not very likely)

    if (cls_dbh>=cBHDclassCount && cls_hd>=cHDclassCount) {

        if (logLevelDebug())

            qDebug() << "DBH AND HD for stamp out of range dbh " << bhd_cm << "and h="<< height_m << "-> using largest available DBH.";

        return m_lookup(cBHDclassCount-1, cHDclassCount-1);

    }

    qDebug() << "ERROR: No stamp defined for dbh " << bhd_cm << "and h="<< height_m;

    throw IException("StampContainer:: did not find a valid stamp.");

}

void StampContainer::attachReaderStamps(const StampContainer &source)

{

    int found=0, total=0;

    foreach (const StampItem &si, m_stamps) {

        const Stamp *s = source.readerStamp(si.crown_radius);

        si.stamp->setReader(const_cast<Stamp*>(s));

        if (s) found++;

        total++;

        //si.crown_radius

    }

    if (logLevelInfo()) qDebug() << "attachReaderStamps: found" << found << "stamps of" << total;

}

void StampContainer::invert()

{

    StampItem  si;

    foreach(si, m_stamps) {

        Stamp *s =si.stamp;

        float *p = s->data();

        while (p!=s->end()) {

            *p = 1. - *p;

            p++;

        }

    }

}

/// convenience function that loads stamps directly from a single file.

void StampContainer::load(const QString &fileName)

{

    QFile readerfile(fileName);

    if (!readerfile.exists())

        throw IException(QString("The LIP stampfile %1 cannot be found!").arg(fileName));

    m_fileName = fileName;

    readerfile.open(QIODevice::ReadOnly);

    QDataStream rin(&readerfile);

    qDebug() << "loading stamp file" << fileName;

    load(rin);

    readerfile.close();

}

void StampContainer::load(QDataStream &in)

{

    qint32 type;

    qint32 count;

    float dbh, hdvalue, crownradius;

    quint32 magic;

    in >> magic;

    if (magic!=0xFEED0001)

        throw IException("StampContainer: invalid file type!");

    quint16 version;

    in >> version;

    if (version != 100)

        throw IException(QString("StampContainer: invalid file version: %1").arg(version));

    in.setVersion(QDataStream::Qt_4_5);

    in >> count; // read count of stamps

    if (logLevelInfo()) qDebug() << count << "stamps to read";

    QString desc;

    in >> desc; // read textual description of stamp

    if (logLevelInfo()) qDebug() << "Stamp notes:" << desc;

    m_desc = desc;

    for (int i=0;i<count;i++) {

        in >> type; // read type

        in >> dbh;

        in >> hdvalue;

        in >> crownradius;

        //qDebug() << "stamp bhd hdvalue type readsum dominance type" << bhd << hdvalue << type << readsum << domvalue << type;

        Stamp *stamp = new Stamp(type);

        stamp->load(in);

        if (dbh > 0.f)

            addStamp(stamp, dbh, hdvalue, crownradius);

        else

            addReaderStamp(stamp, crownradius);

    }

    finalizeSetup(); // fill up lookup grid

    if (count==0)

        throw IException("no stamps loaded!");

}

/** Saves all stamps of the container to a binary stream.

  Format: * count of stamps (int32)

          * a string containing a description (free text) (QString)

      for each stamp:

      - type (enum Stamp::StampType, 4, 8, 12, 16, ...)

      - bhd of the stamp (float)

      - hd-value of the tree (float)

      - crownradius of the stamp (float) in [m]

      - the sum of values in the center of the stamp (used for read out)

      - the dominance value of the stamp

      - individual data values (Stamp::save() / Stamp::load())

      -- offset (int) no. of pixels away from center

      -- list of data items (type*type items)

      see also stamp creation (FonStudio application, MainWindow.cpp).

*/

void StampContainer::save(QDataStream &out)

{

    qint32 type;

    qint32 size = m_stamps.count();

    out << (quint32)0xFEED0001; // magic number

    out << (quint16)100; // version

    out.setVersion(QDataStream::Qt_4_5);

    out << size; // count of stamps...

    out << m_desc; // text...

    foreach(StampItem si, m_stamps) {

        type = si.stamp->dataSize();

        out << type;

        out << si.dbh;

        out << si.hd;

        out << si.crown_radius;

        si.stamp->save(out);

    }

}

QString StampContainer::dump()

{

    QString res;

    QString line;

    int x,y;

    int maxidx;

    res = QString("****** Dump of StampContainer %1 **********\r\n").arg(m_fileName);

    foreach (StampItem si, m_stamps) {

        line = QString("%5 -> size: %1 offset: %2 dbh: %3 hd-ratio: %4\r\n")

               .arg(sqrt((double)si.stamp->count())).arg(si.stamp->offset())

               .arg(si.dbh).arg(si.hd).arg((ptrdiff_t)si.stamp, 0, 16);

        // add data....

        maxidx = 2*si.stamp->offset() + 1;

        for (y=0;y<maxidx;++y)  {

            for (x=0;x<maxidx;++x)  {

                line+= QString::number(*si.stamp->data(x,y)) + " ";

            }

            line+="\r\n";

        }

        line+="==============================================\r\n";

        res+=line;

    }

    res+= "Dump of lookup map\r\n=====================\r\n";

    for (Stamp **s = m_lookup.begin(); s!=m_lookup.end(); ++s) {

        if (*s)

         res += QString("P: x/y: %1/%2 addr %3\r\n").arg( m_lookup.indexOf(s).x()).arg(m_lookup.indexOf(s).y()).arg((ptrdiff_t)*s, 0, 16);

    }

    res+="\r\n" + gridToString(m_lookup);

    return res;

}