/src/dcmtk/dcmimgle/libsrc/didispfn.cc

Source
/*
 *
 *  Copyright (C) 1999-2024, OFFIS e.V.
 *  All rights reserved.  See COPYRIGHT file for details.
 *
 *  This software and supporting documentation were developed by
 *
 *    OFFIS e.V.
 *    R&D Division Health
 *    Escherweg 2
 *    D-26121 Oldenburg, Germany
 *
 *
 *  Module:  dcmimgle
 *
 *  Author:  Joerg Riesmeier
 *
 *  Purpose: DicomDisplayFunction (Source)
 *
 */


#include "dcmtk/config/osconfig.h"

#include "dcmtk/ofstd/ofstd.h"
#include "dcmtk/ofstd/ofbmanip.h"
#include "dcmtk/dcmimgle/didispfn.h"
#include "dcmtk/dcmimgle/displint.h"
#include "dcmtk/dcmimgle/dicrvfit.h"
#include "dcmtk/dcmimgle/didislut.h"
#include "dcmtk/ofstd/ofstream.h"

/*----------------------------*
 *  constant initializations  *
 *----------------------------*/

const int DiDisplayFunction::MinBits = 2;
const int DiDisplayFunction::MaxBits = 16;


/*----------------*
 *  constructors  *
 *----------------*/

DiDisplayFunction::DiDisplayFunction(const char *filename,
                                     const E_DeviceType deviceType,
                                     const signed int ord)
  : Valid(0),
    DeviceType(deviceType),
    ValueCount(0),
    MaxDDLValue(0),
    Order(0),
    AmbientLight(0),
    Illumination(0),
    MinDensity(-1),
    MaxDensity(-1),
    DDLValue(NULL),
    LODValue(NULL),
    MinValue(0),
    MaxValue(0)
{
    OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
    if (readConfigFile(filename))
    {
        /* overwrite file setting for polynomial order */
        if (ord >= 0)
            Order = ord;
        Valid = createSortedTable(DDLValue, LODValue) && calculateMinMax() && interpolateValues();
    }
}


DiDisplayFunction::DiDisplayFunction(const double *val_tab,             // UNTESTED !!
                                     const unsigned long count,
                                     const Uint16 max,
                                     const E_DeviceType deviceType,
                                     const signed int ord)
  : Valid(0),
    DeviceType(deviceType),
    ValueCount(count),
    MaxDDLValue(max),
    Order(ord),
    AmbientLight(0),
    Illumination(0),
    MinDensity(-1),
    MaxDensity(-1),
    DDLValue(NULL),
    LODValue(NULL),
    MinValue(0),
    MaxValue(0)
{
    OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
    /* check number of entries */
    if ((ValueCount > 0) && (ValueCount == OFstatic_cast(unsigned long, MaxDDLValue) + 1))
    {
        /* copy value table */
        DDLValue = new Uint16[ValueCount];
        LODValue = new double[ValueCount];
        if ((DDLValue != NULL) && (LODValue != NULL))
        {
            unsigned int i;
            for (i = 0; i <= MaxDDLValue; ++i)
            {
                DDLValue[i] = OFstatic_cast(Uint16, i);     // set DDL values
                LODValue[i] = val_tab[i];                   // copy table
            }
            Valid = calculateMinMax();
        }
    }
}


DiDisplayFunction::DiDisplayFunction(const Uint16 *ddl_tab,             // UNTESTED !!
                                     const double *val_tab,
                                     const unsigned long count,
                                     const Uint16 max,
                                     const E_DeviceType deviceType,
                                     const signed int ord)
  : Valid(0),
    DeviceType(deviceType),
    ValueCount(count),
    MaxDDLValue(max),
    Order(ord),
    AmbientLight(0),
    Illumination(0),
    MinDensity(-1),
    MaxDensity(-1),
    DDLValue(NULL),
    LODValue(NULL),
    MinValue(0),
    MaxValue(0)
{
    OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
    /* check for maximum number of entries */
    if (ValueCount <= MAX_TABLE_ENTRY_COUNT)
        Valid = createSortedTable(ddl_tab, val_tab) && calculateMinMax() && interpolateValues();
}


DiDisplayFunction::DiDisplayFunction(const double val_min,
                                     const double val_max,
                                     const unsigned long count,
                                     const E_DeviceType deviceType,
                                     const signed int ord)
  : Valid(0),
    DeviceType(deviceType),
    ValueCount(count),
    MaxDDLValue(0),
    Order(ord),
    AmbientLight(0),
    Illumination(0),
    MinDensity(-1),
    MaxDensity(-1),
    DDLValue(NULL),
    LODValue(NULL),
    MinValue(val_min),
    MaxValue(val_max)
{
    OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
    /* check parameters */
    if ((ValueCount > 1) && (ValueCount <= MAX_TABLE_ENTRY_COUNT) && (MinValue < MaxValue))
    {
        /* create value tables */
        MaxDDLValue = OFstatic_cast(Uint16, count - 1);
        DDLValue = new Uint16[ValueCount];
        LODValue = new double[ValueCount];
        if ((DDLValue != NULL) && (LODValue != NULL))
        {
            Uint16 i;
            const double min = ((DeviceType == EDT_Printer) || (DeviceType == EDT_Scanner)) ? val_max : val_min;
            const double max = ((DeviceType == EDT_Printer) || (DeviceType == EDT_Scanner)) ? val_min : val_max;
            const double val = (max - min) / OFstatic_cast(double, MaxDDLValue);
            DDLValue[0] = 0;
            LODValue[0] = min;
            for (i = 1; i < MaxDDLValue; ++i)
            {
                DDLValue[i] = i;                            // set DDL values
                LODValue[i] = LODValue[i - 1] + val;        // compute luminance/OD value
            }
            DDLValue[MaxDDLValue] = MaxDDLValue;
            LODValue[MaxDDLValue] = max;
            Valid = 1;
        }
    }
}


/*--------------*
 *  destructor  *
 *--------------*/

DiDisplayFunction::~DiDisplayFunction()
{
    delete[] DDLValue;
    delete[] LODValue;
    int i;
    for (i = 0; i < MAX_NUMBER_OF_TABLES; ++i)
        delete LookupTable[i];
}


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


double DiDisplayFunction::getValueforDDL(const Uint16 ddl) const
{
    if ((LODValue != NULL) && (ddl < ValueCount))
        return LODValue[ddl];
    return -1;
}


Uint16 DiDisplayFunction::getDDLforValue(const double value) const
{
    if ((LODValue != NULL) && (ValueCount > 0))
    {
        unsigned long j = 0;
        /* search for closest index, assuming monotony */
        if ((DeviceType == EDT_Printer) || (DeviceType == EDT_Scanner))
        {
            /* hardcopy device: descending values */
            while ((j + 1 < ValueCount) && (LODValue[j] > value))
                ++j;
        } else {
            /* softcopy device: ascending values */
            while ((j + 1 < ValueCount) && (LODValue[j] < value))
                ++j;
        }
        /* check which value is closer, the upper or the lower */
        if ((j > 0) && (fabs(LODValue[j - 1] - value) < fabs(LODValue[j] - value)))
            --j;
        return OFstatic_cast(Uint16, j);
    }
    return 0;
}


const DiDisplayLUT *DiDisplayFunction::getLookupTable(const int bits,
                                                      unsigned long count)
{
    if (Valid && (bits >= MinBits) && (bits <= MaxBits))
    {
        const int idx = bits - MinBits;
        /* automatically compute number of entries */
        if (count == 0)
            count = DicomImageClass::maxval(bits, 0);
        /* check whether existing LUT is still valid */
        if ((LookupTable[idx] != NULL) && ((count != LookupTable[idx]->getCount()) ||
            (AmbientLight != LookupTable[idx]->getAmbientLightValue()) ||
            (Illumination != LookupTable[idx]->getIlluminationValue())))
        {
            delete LookupTable[idx];
            LookupTable[idx] = NULL;
        }
        if (LookupTable[idx] == NULL)                             // first calculation of this LUT
            LookupTable[idx] = getDisplayLUT(count);
        return LookupTable[idx];
    }
    return NULL;
}


int DiDisplayFunction::deleteLookupTable(const int bits)
{
    if (bits == 0)
    {
        /* delete all LUTs */
        int i;
        for (i = 0; i < MAX_NUMBER_OF_TABLES; ++i)
        {
            delete LookupTable[i];
            LookupTable[i] = NULL;
        }
        return 1;
    }
    else if ((bits >= MinBits) && (bits <= MaxBits))
    {
        /* delete the specified LUT */
        const int idx = bits - MinBits;
        if (LookupTable[idx] != NULL)
        {
            delete LookupTable[idx];
            LookupTable[idx] = NULL;
            return 1;
        }
        return 2;
    }
    return 0;
}


int DiDisplayFunction::setAmbientLightValue(const double value)
{
    if (value >= 0)
    {
        AmbientLight = value;
        return 1;
    }
    return 0;
}


int DiDisplayFunction::setIlluminationValue(const double value)
{
    if (value >= 0)
    {
        Illumination = value;
        return 1;
    }
    return 0;
}


int DiDisplayFunction::setMinDensityValue(const double value)
{
    MinDensity = value;
    return (value < 0) ? 2 : 1;
}


int DiDisplayFunction::setMaxDensityValue(const double value)
{
    MaxDensity = value;
    return (value < 0) ? 2 : 1;
}


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


int DiDisplayFunction::readConfigFile(const char *filename)
{
    if ((filename != NULL) && (strlen(filename) > 0))
    {
        STD_NAMESPACE ifstream file(filename, OFopenmode_in_nocreate);
        if (file)
        {
            char c;
            while (file.get(c))
            {
                if (c == '#')                                               // comment character
                {
                    while (file.get(c) && (c != '\n') && (c != '\r'));      // skip comments
                }
                else if (!OFStandard::isspace(c))                           // skip whitespaces
                {
                    file.putback(c);
                    if (MaxDDLValue == 0)                                   // read maxvalue
                    {
                        char str[4];
                        file.get(str, sizeof(str));
                        if (strcmp(str, "max") == 0)                        // check for key word: max
                        {
                            file >> MaxDDLValue;
                            if (MaxDDLValue > 0)
                            {
                                DDLValue = new Uint16[OFstatic_cast(unsigned long, MaxDDLValue) + 1];
                                LODValue = new double[OFstatic_cast(unsigned long, MaxDDLValue) + 1];
                                if ((DDLValue == NULL) || (LODValue == NULL))
                                    return 0;
                            } else {
                                DCMIMGLE_ERROR("invalid or missing value for maximum DDL value in DISPLAY file");
                                return 0;                                   // abort
                            }
                        } else {
                            DCMIMGLE_ERROR("missing keyword 'max' for maximum DDL value in DISPLAY file");
                            return 0;                                       // abort
                        }
                    }
                    else if ((AmbientLight == 0.0) && (c == 'a'))           // read ambient light value (optional)
                    {
                        char str[4];
                        file.get(str, sizeof(str));
                        if (strcmp(str, "amb") == 0)                        // check for key word: amb
                        {
                            file >> AmbientLight;
                            if (AmbientLight < 0)
                            {
                                DCMIMGLE_WARN("invalid value for ambient light in DISPLAY file ... ignoring");
                                AmbientLight = 0;
                            }
                        } else {
                            DCMIMGLE_ERROR("invalid DISPLAY file ... ignoring");
                            return 0;                                       // abort
                        }
                    }
                    else if ((Illumination == 0.0) && (c == 'l'))           // read ambient light value (optional)
                    {
                        char str[4];
                        file.get(str, sizeof(str));
                        if (strcmp(str, "lum") == 0)                        // check for key word: lum
                        {
                            file >> Illumination;
                            if (Illumination < 0)
                            {
                                DCMIMGLE_WARN("invalid value for illumination in DISPLAY file ... ignoring");
                                Illumination = 0;
                            }
                        } else {
                            DCMIMGLE_ERROR("invalid DISPLAY file ... ignoring");
                            return 0;                                       // abort
                        }
                    }
                    else if ((Order == 0) && (c == 'o'))                    // read polynomial order (optional)
                    {
                        char str[4];
                        file.get(str, sizeof(str));
                        if (strcmp(str, "ord") == 0)                        // check for key word: ord
                        {
                            file >> Order;
                            if (Order < 0)
                            {
                                DCMIMGLE_WARN("invalid value for polynomial order in DISPLAY file ... ignoring");
                                Order = 0;
                            }
                        } else {
                            DCMIMGLE_ERROR("invalid DISPLAY file ... ignoring");
                            return 0;                                       // abort
                        }
                    } else {
                        if (ValueCount <= OFstatic_cast(unsigned long, MaxDDLValue))
                        {
                            file >> DDLValue[ValueCount];                   // read DDL value
                            file >> LODValue[ValueCount];                   // read luminance/OD value
                            if (file.fail())
                            {
                                DCMIMGLE_WARN("missing luminance/OD value in DISPLAY file ... ignoring last entry");
                            }
                            else if (DDLValue[ValueCount] > MaxDDLValue)
                            {
                                DCMIMGLE_WARN("DDL value (" << DDLValue[ValueCount] << ") exceeds maximum value ("
                                    << MaxDDLValue << ") in DISPLAY file ... ignoring value");
                            } else
                                ++ValueCount;
                        } else {
                            DCMIMGLE_WARN("too many values in DISPLAY file ... ignoring last line(s)");
                            return 2;
                        }
                    }
                }
            }
            if ((MaxDDLValue > 0) && (ValueCount > 0))
                return ((DDLValue != NULL) && (LODValue != NULL));
            else {
                DCMIMGLE_WARN("invalid DISPLAY file ... ignoring");
            }
        } else {
            DCMIMGLE_WARN("can't open DISPLAY file ... ignoring");
        }
    }
    return 0;
}


int DiDisplayFunction::createSortedTable(const Uint16 *ddl_tab,
                                         const double *val_tab)
{
    int status = 0;
    Uint16 *old_ddl = DDLValue;
    double *old_val = LODValue;
    if ((ValueCount > 0) && (ddl_tab != NULL) && (val_tab != NULL))
    {
        const unsigned long count = OFstatic_cast(unsigned long, MaxDDLValue) + 1;
        DDLValue = new Uint16[ValueCount];
        LODValue = new double[ValueCount];
        Sint32 *sort_tab = new Sint32[count];                                       // auxiliary array (temporary)
        if ((DDLValue != NULL) && (LODValue != NULL) && (sort_tab != NULL))
        {
            OFBitmanipTemplate<Sint32>::setMem(sort_tab, -1, count);                // initialize array
            unsigned long i;
            for (i = 0; i < ValueCount; ++i)
            {
                if (ddl_tab[i] <= MaxDDLValue)                                      // calculate sort table
                    sort_tab[ddl_tab[i]] = i;
            }
            ValueCount = 0;
            for (i = 0; i <= MaxDDLValue; ++i)                                      // sort ascending
            {
                if (sort_tab[i] >= 0)
                {
                    DDLValue[ValueCount] = ddl_tab[sort_tab[i]];
                    LODValue[ValueCount] = (val_tab[sort_tab[i]] > 0) ? val_tab[sort_tab[i]] : 0;
                    ++ValueCount;                                                   // re-count to ignore values exceeding max
                }
            }
            i = 1;
            if ((DeviceType == EDT_Printer) || (DeviceType == EDT_Scanner))
            {
                /* hardcopy device: check for monotonous descending OD values */
                while ((i < ValueCount) && (LODValue[i - 1] >= LODValue[i]))
                    ++i;
                if (i < ValueCount)
                {
                    DCMIMGLE_WARN("OD values (ordered by DDLs) don't descend monotonously");
                }
            } else {
                /* softcopy device: check for monotonous ascending luminance values */
                while ((i < ValueCount) && (LODValue[i - 1] <= LODValue[i]))
                    ++i;
                if (i < ValueCount)
                {
                    DCMIMGLE_WARN("luminance values (ordered by DDLs) don't ascend monotonously");
                }
            }
            status = (ValueCount > 0);
        }
        delete[] sort_tab;
    }
    delete[] old_ddl;
    delete[] old_val;
    return status;
}


int DiDisplayFunction::interpolateValues()
{
    if (ValueCount <= OFstatic_cast(unsigned long, MaxDDLValue))        // interpolation necessary ?
    {
        int status = 0;
        if (Order > 0)
        {
            /* use polynomial curve fitting */
            double *coeff = new double[Order + 1];
            /* compute coefficients */
            if ((coeff != NULL) && DiCurveFitting<Uint16, double>::calculateCoefficients(DDLValue, LODValue,
                    OFstatic_cast(unsigned int, ValueCount), Order, coeff))
            {
                /* delete old data arrays */
                delete[] DDLValue;
                delete[] LODValue;
                /* create new data arrays */
                ValueCount = OFstatic_cast(unsigned long, MaxDDLValue) + 1;
                DDLValue = new Uint16[ValueCount];
                LODValue = new double[ValueCount];
                if ((DDLValue != NULL) && (LODValue != NULL))
                {
                    /* set x values linearly */
                    unsigned int i;
                    for (i = 0; i <= MaxDDLValue; ++i)
                        DDLValue[i] = OFstatic_cast(Uint16, i);
                    /* compute new y values */
                    status = DiCurveFitting<Uint16, double>::calculateValues(0, MaxDDLValue, LODValue,
                        OFstatic_cast(unsigned int, ValueCount), Order, coeff);
                }
            }
            delete[] coeff;
        } else {
            /* use cubic spline interpolation */
            double *spline = new double[ValueCount];
            if ((spline != NULL) &&
                DiCubicSpline<Uint16, double>::Function(DDLValue, LODValue, OFstatic_cast(unsigned int, ValueCount), spline))
            {
                /* save old values */
                const unsigned long count = ValueCount;
                Uint16 *old_ddl = DDLValue;
                double *old_val = LODValue;
                /* create new data arrays */
                ValueCount = OFstatic_cast(unsigned long, MaxDDLValue) + 1;
                DDLValue = new Uint16[ValueCount];
                LODValue = new double[ValueCount];
                if ((DDLValue != NULL) && (LODValue != NULL))
                {
                    /* set x values linearly */
                    unsigned int i;
                    for (i = 0; i <= MaxDDLValue; ++i)
                        DDLValue[i] = OFstatic_cast(Uint16, i);
                    /* compute new y values */
                    status = DiCubicSpline<Uint16, double>::Interpolation(old_ddl, old_val, spline,
                        OFstatic_cast(unsigned int, count), DDLValue, LODValue, OFstatic_cast(unsigned int, ValueCount));
                }
                /* delete old data arrays */
                delete[] old_ddl;
                delete[] old_val;
            }
            delete[] spline;
        }
        return status;
    }
    return 2;
}


int DiDisplayFunction::calculateMinMax()
{
    if ((LODValue != NULL) && (ValueCount > 0))
    {
        MinValue = LODValue[0];
        MaxValue = LODValue[0];
        unsigned long i;
        for (i = 1; i < ValueCount; ++i)
        {
            if (LODValue[i] < MinValue)
                MinValue = LODValue[i];
            if (LODValue[i] > MaxValue)
                MaxValue = LODValue[i];
        }
        return 1;
    }
    return 0;
}


int DiDisplayFunction::checkMinMaxDensity() const
{
    if ((MinDensity >= 0) && (MaxDensity >= 0) && (MinDensity >= MaxDensity))
    {
        DCMIMGLE_WARN("invalid optical density range (Dmin = " << MinDensity << ", Dmax = " << MaxDensity << ")");
        return 0;
    }
    return 1;
}


double DiDisplayFunction::getMinLuminanceValue() const
{
    /* Dmax = -1 means unspecified */
    return (MaxDensity < 0) ? -1 : convertODtoLum(MaxDensity);
}


double DiDisplayFunction::getMaxLuminanceValue() const
{
    /* Dmin = -1 means unspecified */
    return (MinDensity < 0) ? -1 : convertODtoLum(MinDensity);
}


double *DiDisplayFunction::convertODtoLumTable(const double *od_tab,
                                               const unsigned long count,
                                               const OFBool useAmb)
{
    double *lum_tab = NULL;
    if ((od_tab != NULL) && (count > 0))
    {
        /* create a new table for the luminance values */
        lum_tab = new double[count];
        if (lum_tab != NULL)
        {
            /* compute luminance values from optical density */
            unsigned int i;
            if (useAmb)
            {
                for (i = 0; i < count; ++i)
                    lum_tab[i] = AmbientLight + Illumination * pow(OFstatic_cast(double, 10), -od_tab[i]);
            } else {
                /* ambient light is added later */
                for (i = 0; i < count; ++i)
                    lum_tab[i] = Illumination * pow(OFstatic_cast(double, 10), -od_tab[i]);
            }
        }
    }
    return lum_tab;
}


double DiDisplayFunction::convertODtoLum(const double value,
                                         const OFBool useAmb) const
{
    return (useAmb) ? convertODtoLum(value, AmbientLight, Illumination) :
                      convertODtoLum(value, 0, Illumination);
}


double DiDisplayFunction::convertODtoLum(const double value,
                                         const double ambient,
                                         const double illum)
{
    /* formula from DICOM PS3.14: L = La + L0 * 10^-D */
    return (value >= 0) && (ambient >= 0) && (illum >= 0) ?
        ambient + illum * pow(OFstatic_cast(double, 10), -value) : -1 /*invalid*/;
}

Coverage Report

Created: 2026-09-28 07:04

Line	Count	Source
1		/*
2		*
3		* Copyright (C) 1999-2024, OFFIS e.V.
4		* All rights reserved. See COPYRIGHT file for details.
5		*
6		* This software and supporting documentation were developed by
7		*
8		* OFFIS e.V.
9		* R&D Division Health
10		* Escherweg 2
11		* D-26121 Oldenburg, Germany
12		*
13		*
14		* Module: dcmimgle
15		*
16		* Author: Joerg Riesmeier
17		*
18		* Purpose: DicomDisplayFunction (Source)
19		*
20		*/
21
22
23		#include "dcmtk/config/osconfig.h"
24
25		#include "dcmtk/ofstd/ofstd.h"
26		#include "dcmtk/ofstd/ofbmanip.h"
27		#include "dcmtk/dcmimgle/didispfn.h"
28		#include "dcmtk/dcmimgle/displint.h"
29		#include "dcmtk/dcmimgle/dicrvfit.h"
30		#include "dcmtk/dcmimgle/didislut.h"
31		#include "dcmtk/ofstd/ofstream.h"
32
33		/----------------------------
34		* constant initializations *
35		----------------------------/
36
37		const int DiDisplayFunction::MinBits = 2;
38		const int DiDisplayFunction::MaxBits = 16;
39
40
41		/----------------
42		* constructors *
43		----------------/
44
45		DiDisplayFunction::DiDisplayFunction(const char *filename,
46		const E_DeviceType deviceType,
47		const signed int ord)
48	0	: Valid(0),
49	0	DeviceType(deviceType),
50	0	ValueCount(0),
51	0	MaxDDLValue(0),
52	0	Order(0),
53	0	AmbientLight(0),
54	0	Illumination(0),
55	0	MinDensity(-1),
56	0	MaxDensity(-1),
57	0	DDLValue(NULL),
58	0	LODValue(NULL),
59	0	MinValue(0),
60	0	MaxValue(0)
61	0	{
62	0	OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
63	0	if (readConfigFile(filename))
64	0	{
65		/* overwrite file setting for polynomial order */
66	0	if (ord >= 0)
67	0	Order = ord;
68	0	Valid = createSortedTable(DDLValue, LODValue) && calculateMinMax() && interpolateValues();
69	0	}
70	0	}
71
72
73		DiDisplayFunction::DiDisplayFunction(const double *val_tab, // UNTESTED !!
74		const unsigned long count,
75		const Uint16 max,
76		const E_DeviceType deviceType,
77		const signed int ord)
78	0	: Valid(0),
79	0	DeviceType(deviceType),
80	0	ValueCount(count),
81	0	MaxDDLValue(max),
82	0	Order(ord),
83	0	AmbientLight(0),
84	0	Illumination(0),
85	0	MinDensity(-1),
86	0	MaxDensity(-1),
87	0	DDLValue(NULL),
88	0	LODValue(NULL),
89	0	MinValue(0),
90	0	MaxValue(0)
91	0	{
92	0	OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
93		/* check number of entries */
94	0	if ((ValueCount > 0) && (ValueCount == OFstatic_cast(unsigned long, MaxDDLValue) + 1))
95	0	{
96		/* copy value table */
97	0	DDLValue = new Uint16[ValueCount];
98	0	LODValue = new double[ValueCount];
99	0	if ((DDLValue != NULL) && (LODValue != NULL))
100	0	{
101	0	unsigned int i;
102	0	for (i = 0; i <= MaxDDLValue; ++i)
103	0	{
104	0	DDLValue[i] = OFstatic_cast(Uint16, i); // set DDL values
105	0	LODValue[i] = val_tab[i]; // copy table
106	0	}
107	0	Valid = calculateMinMax();
108	0	}
109	0	}
110	0	}
111
112
113		DiDisplayFunction::DiDisplayFunction(const Uint16 *ddl_tab, // UNTESTED !!
114		const double *val_tab,
115		const unsigned long count,
116		const Uint16 max,
117		const E_DeviceType deviceType,
118		const signed int ord)
119	0	: Valid(0),
120	0	DeviceType(deviceType),
121	0	ValueCount(count),
122	0	MaxDDLValue(max),
123	0	Order(ord),
124	0	AmbientLight(0),
125	0	Illumination(0),
126	0	MinDensity(-1),
127	0	MaxDensity(-1),
128	0	DDLValue(NULL),
129	0	LODValue(NULL),
130	0	MinValue(0),
131	0	MaxValue(0)
132	0	{
133	0	OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
134		/* check for maximum number of entries */
135	0	if (ValueCount <= MAX_TABLE_ENTRY_COUNT)
136	0	Valid = createSortedTable(ddl_tab, val_tab) && calculateMinMax() && interpolateValues();
137	0	}
138
139
140		DiDisplayFunction::DiDisplayFunction(const double val_min,
141		const double val_max,
142		const unsigned long count,
143		const E_DeviceType deviceType,
144		const signed int ord)
145	0	: Valid(0),
146	0	DeviceType(deviceType),
147	0	ValueCount(count),
148	0	MaxDDLValue(0),
149	0	Order(ord),
150	0	AmbientLight(0),
151	0	Illumination(0),
152	0	MinDensity(-1),
153	0	MaxDensity(-1),
154	0	DDLValue(NULL),
155	0	LODValue(NULL),
156	0	MinValue(val_min),
157	0	MaxValue(val_max)
158	0	{
159	0	OFBitmanipTemplate<DiDisplayLUT *>::zeroMem(LookupTable, MAX_NUMBER_OF_TABLES);
160		/* check parameters */
161	0	if ((ValueCount > 1) && (ValueCount <= MAX_TABLE_ENTRY_COUNT) && (MinValue < MaxValue))
162	0	{
163		/* create value tables */
164	0	MaxDDLValue = OFstatic_cast(Uint16, count - 1);
165	0	DDLValue = new Uint16[ValueCount];
166	0	LODValue = new double[ValueCount];
167	0	if ((DDLValue != NULL) && (LODValue != NULL))
168	0	{
169	0	Uint16 i;
170	0	const double min = ((DeviceType == EDT_Printer) \|\| (DeviceType == EDT_Scanner)) ? val_max : val_min;
171	0	const double max = ((DeviceType == EDT_Printer) \|\| (DeviceType == EDT_Scanner)) ? val_min : val_max;
172	0	const double val = (max - min) / OFstatic_cast(double, MaxDDLValue);
173	0	DDLValue[0] = 0;
174	0	LODValue[0] = min;
175	0	for (i = 1; i < MaxDDLValue; ++i)
176	0	{
177	0	DDLValue[i] = i; // set DDL values
178	0	LODValue[i] = LODValue[i - 1] + val; // compute luminance/OD value
179	0	}
180	0	DDLValue[MaxDDLValue] = MaxDDLValue;
181	0	LODValue[MaxDDLValue] = max;
182	0	Valid = 1;
183	0	}
184	0	}
185	0	}
186
187
188		/--------------
189		* destructor *
190		--------------/
191
192		DiDisplayFunction::~DiDisplayFunction()
193	0	{
194	0	delete[] DDLValue;
195	0	delete[] LODValue;
196	0	int i;
197	0	for (i = 0; i < MAX_NUMBER_OF_TABLES; ++i)
198	0	delete LookupTable[i];
199	0	}
200
201
202		/********************************************************************/
203
204
205		double DiDisplayFunction::getValueforDDL(const Uint16 ddl) const
206	0	{
207	0	if ((LODValue != NULL) && (ddl < ValueCount))
208	0	return LODValue[ddl];
209	0	return -1;
210	0	}
211
212
213		Uint16 DiDisplayFunction::getDDLforValue(const double value) const
214	0	{
215	0	if ((LODValue != NULL) && (ValueCount > 0))
216	0	{
217	0	unsigned long j = 0;
218		/* search for closest index, assuming monotony */
219	0	if ((DeviceType == EDT_Printer) \|\| (DeviceType == EDT_Scanner))
220	0	{
221		/* hardcopy device: descending values */
222	0	while ((j + 1 < ValueCount) && (LODValue[j] > value))
223	0	++j;
224	0	} else {
225		/* softcopy device: ascending values */
226	0	while ((j + 1 < ValueCount) && (LODValue[j] < value))
227	0	++j;
228	0	}
229		/* check which value is closer, the upper or the lower */
230	0	if ((j > 0) && (fabs(LODValue[j - 1] - value) < fabs(LODValue[j] - value)))
231	0	--j;
232	0	return OFstatic_cast(Uint16, j);
233	0	}
234	0	return 0;
235	0	}
236
237
238		const DiDisplayLUT *DiDisplayFunction::getLookupTable(const int bits,
239		unsigned long count)
240	0	{
241	0	if (Valid && (bits >= MinBits) && (bits <= MaxBits))
242	0	{
243	0	const int idx = bits - MinBits;
244		/* automatically compute number of entries */
245	0	if (count == 0)
246	0	count = DicomImageClass::maxval(bits, 0);
247		/* check whether existing LUT is still valid */
248	0	if ((LookupTable[idx] != NULL) && ((count != LookupTable[idx]->getCount()) \|\|
249	0	(AmbientLight != LookupTable[idx]->getAmbientLightValue()) \|\|
250	0	(Illumination != LookupTable[idx]->getIlluminationValue())))
251	0	{
252	0	delete LookupTable[idx];
253	0	LookupTable[idx] = NULL;
254	0	}
255	0	if (LookupTable[idx] == NULL) // first calculation of this LUT
256	0	LookupTable[idx] = getDisplayLUT(count);
257	0	return LookupTable[idx];
258	0	}
259	0	return NULL;
260	0	}
261
262
263		int DiDisplayFunction::deleteLookupTable(const int bits)
264	0	{
265	0	if (bits == 0)
266	0	{
267		/* delete all LUTs */
268	0	int i;
269	0	for (i = 0; i < MAX_NUMBER_OF_TABLES; ++i)
270	0	{
271	0	delete LookupTable[i];
272	0	LookupTable[i] = NULL;
273	0	}
274	0	return 1;
275	0	}
276	0	else if ((bits >= MinBits) && (bits <= MaxBits))
277	0	{
278		/* delete the specified LUT */
279	0	const int idx = bits - MinBits;
280	0	if (LookupTable[idx] != NULL)
281	0	{
282	0	delete LookupTable[idx];
283	0	LookupTable[idx] = NULL;
284	0	return 1;
285	0	}
286	0	return 2;
287	0	}
288	0	return 0;
289	0	}
290
291
292		int DiDisplayFunction::setAmbientLightValue(const double value)
293	0	{
294	0	if (value >= 0)
295	0	{
296	0	AmbientLight = value;
297	0	return 1;
298	0	}
299	0	return 0;
300	0	}
301
302
303		int DiDisplayFunction::setIlluminationValue(const double value)
304	0	{
305	0	if (value >= 0)
306	0	{
307	0	Illumination = value;
308	0	return 1;
309	0	}
310	0	return 0;
311	0	}
312
313
314		int DiDisplayFunction::setMinDensityValue(const double value)
315	0	{
316	0	MinDensity = value;
317	0	return (value < 0) ? 2 : 1;
318	0	}
319
320
321		int DiDisplayFunction::setMaxDensityValue(const double value)
322	0	{
323	0	MaxDensity = value;
324	0	return (value < 0) ? 2 : 1;
325	0	}
326
327
328		/********************************************************************/
329
330
331		int DiDisplayFunction::readConfigFile(const char *filename)
332	0	{
333	0	if ((filename != NULL) && (strlen(filename) > 0))
334	0	{
335	0	STD_NAMESPACE ifstream file(filename, OFopenmode_in_nocreate);
336	0	if (file)
337	0	{
338	0	char c;
339	0	while (file.get(c))
340	0	{
341	0	if (c == '#') // comment character
342	0	{
343	0	while (file.get(c) && (c != '\n') && (c != '\r')); // skip comments
344	0	}
345	0	else if (!OFStandard::isspace(c)) // skip whitespaces
346	0	{
347	0	file.putback(c);
348	0	if (MaxDDLValue == 0) // read maxvalue
349	0	{
350	0	char str[4];
351	0	file.get(str, sizeof(str));
352	0	if (strcmp(str, "max") == 0) // check for key word: max
353	0	{
354	0	file >> MaxDDLValue;
355	0	if (MaxDDLValue > 0)
356	0	{
357	0	DDLValue = new Uint16[OFstatic_cast(unsigned long, MaxDDLValue) + 1];
358	0	LODValue = new double[OFstatic_cast(unsigned long, MaxDDLValue) + 1];
359	0	if ((DDLValue == NULL) \|\| (LODValue == NULL))
360	0	return 0;
361	0	} else {
362	0	DCMIMGLE_ERROR("invalid or missing value for maximum DDL value in DISPLAY file");
363	0	return 0; // abort
364	0	}
365	0	} else {
366	0	DCMIMGLE_ERROR("missing keyword 'max' for maximum DDL value in DISPLAY file");
367	0	return 0; // abort
368	0	}
369	0	}
370	0	else if ((AmbientLight == 0.0) && (c == 'a')) // read ambient light value (optional)
371	0	{
372	0	char str[4];
373	0	file.get(str, sizeof(str));
374	0	if (strcmp(str, "amb") == 0) // check for key word: amb
375	0	{
376	0	file >> AmbientLight;
377	0	if (AmbientLight < 0)
378	0	{
379	0	DCMIMGLE_WARN("invalid value for ambient light in DISPLAY file ... ignoring");
380	0	AmbientLight = 0;
381	0	}
382	0	} else {
383	0	DCMIMGLE_ERROR("invalid DISPLAY file ... ignoring");
384	0	return 0; // abort
385	0	}
386	0	}
387	0	else if ((Illumination == 0.0) && (c == 'l')) // read ambient light value (optional)
388	0	{
389	0	char str[4];
390	0	file.get(str, sizeof(str));
391	0	if (strcmp(str, "lum") == 0) // check for key word: lum
392	0	{
393	0	file >> Illumination;
394	0	if (Illumination < 0)
395	0	{
396	0	DCMIMGLE_WARN("invalid value for illumination in DISPLAY file ... ignoring");
397	0	Illumination = 0;
398	0	}
399	0	} else {
400	0	DCMIMGLE_ERROR("invalid DISPLAY file ... ignoring");
401	0	return 0; // abort
402	0	}
403	0	}
404	0	else if ((Order == 0) && (c == 'o')) // read polynomial order (optional)
405	0	{
406	0	char str[4];
407	0	file.get(str, sizeof(str));
408	0	if (strcmp(str, "ord") == 0) // check for key word: ord
409	0	{
410	0	file >> Order;
411	0	if (Order < 0)
412	0	{
413	0	DCMIMGLE_WARN("invalid value for polynomial order in DISPLAY file ... ignoring");
414	0	Order = 0;
415	0	}
416	0	} else {
417	0	DCMIMGLE_ERROR("invalid DISPLAY file ... ignoring");
418	0	return 0; // abort
419	0	}
420	0	} else {
421	0	if (ValueCount <= OFstatic_cast(unsigned long, MaxDDLValue))
422	0	{
423	0	file >> DDLValue[ValueCount]; // read DDL value
424	0	file >> LODValue[ValueCount]; // read luminance/OD value
425	0	if (file.fail())
426	0	{
427	0	DCMIMGLE_WARN("missing luminance/OD value in DISPLAY file ... ignoring last entry");
428	0	}
429	0	else if (DDLValue[ValueCount] > MaxDDLValue)
430	0	{
431	0	DCMIMGLE_WARN("DDL value (" << DDLValue[ValueCount] << ") exceeds maximum value ("
432	0	<< MaxDDLValue << ") in DISPLAY file ... ignoring value");
433	0	} else
434	0	++ValueCount;
435	0	} else {
436	0	DCMIMGLE_WARN("too many values in DISPLAY file ... ignoring last line(s)");
437	0	return 2;
438	0	}
439	0	}
440	0	}
441	0	}
442	0	if ((MaxDDLValue > 0) && (ValueCount > 0))
443	0	return ((DDLValue != NULL) && (LODValue != NULL));
444	0	else {
445	0	DCMIMGLE_WARN("invalid DISPLAY file ... ignoring");
446	0	}
447	0	} else {
448	0	DCMIMGLE_WARN("can't open DISPLAY file ... ignoring");
449	0	}
450	0	}
451	0	return 0;
452	0	}
453
454
455		int DiDisplayFunction::createSortedTable(const Uint16 *ddl_tab,
456		const double *val_tab)
457	0	{
458	0	int status = 0;
459	0	Uint16 *old_ddl = DDLValue;
460	0	double *old_val = LODValue;
461	0	if ((ValueCount > 0) && (ddl_tab != NULL) && (val_tab != NULL))
462	0	{
463	0	const unsigned long count = OFstatic_cast(unsigned long, MaxDDLValue) + 1;
464	0	DDLValue = new Uint16[ValueCount];
465	0	LODValue = new double[ValueCount];
466	0	Sint32 *sort_tab = new Sint32[count]; // auxiliary array (temporary)
467	0	if ((DDLValue != NULL) && (LODValue != NULL) && (sort_tab != NULL))
468	0	{
469	0	OFBitmanipTemplate<Sint32>::setMem(sort_tab, -1, count); // initialize array
470	0	unsigned long i;
471	0	for (i = 0; i < ValueCount; ++i)
472	0	{
473	0	if (ddl_tab[i] <= MaxDDLValue) // calculate sort table
474	0	sort_tab[ddl_tab[i]] = i;
475	0	}
476	0	ValueCount = 0;
477	0	for (i = 0; i <= MaxDDLValue; ++i) // sort ascending
478	0	{
479	0	if (sort_tab[i] >= 0)
480	0	{
481	0	DDLValue[ValueCount] = ddl_tab[sort_tab[i]];
482	0	LODValue[ValueCount] = (val_tab[sort_tab[i]] > 0) ? val_tab[sort_tab[i]] : 0;
483	0	++ValueCount; // re-count to ignore values exceeding max
484	0	}
485	0	}
486	0	i = 1;
487	0	if ((DeviceType == EDT_Printer) \|\| (DeviceType == EDT_Scanner))
488	0	{
489		/* hardcopy device: check for monotonous descending OD values */
490	0	while ((i < ValueCount) && (LODValue[i - 1] >= LODValue[i]))
491	0	++i;
492	0	if (i < ValueCount)
493	0	{
494	0	DCMIMGLE_WARN("OD values (ordered by DDLs) don't descend monotonously");
495	0	}
496	0	} else {
497		/* softcopy device: check for monotonous ascending luminance values */
498	0	while ((i < ValueCount) && (LODValue[i - 1] <= LODValue[i]))
499	0	++i;
500	0	if (i < ValueCount)
501	0	{
502	0	DCMIMGLE_WARN("luminance values (ordered by DDLs) don't ascend monotonously");
503	0	}
504	0	}
505	0	status = (ValueCount > 0);
506	0	}
507	0	delete[] sort_tab;
508	0	}
509	0	delete[] old_ddl;
510	0	delete[] old_val;
511	0	return status;
512	0	}
513
514
515		int DiDisplayFunction::interpolateValues()
516	0	{
517	0	if (ValueCount <= OFstatic_cast(unsigned long, MaxDDLValue)) // interpolation necessary ?
518	0	{
519	0	int status = 0;
520	0	if (Order > 0)
521	0	{
522		/* use polynomial curve fitting */
523	0	double *coeff = new double[Order + 1];
524		/* compute coefficients */
525	0	if ((coeff != NULL) && DiCurveFitting<Uint16, double>::calculateCoefficients(DDLValue, LODValue,
526	0	OFstatic_cast(unsigned int, ValueCount), Order, coeff))
527	0	{
528		/* delete old data arrays */
529	0	delete[] DDLValue;
530	0	delete[] LODValue;
531		/* create new data arrays */
532	0	ValueCount = OFstatic_cast(unsigned long, MaxDDLValue) + 1;
533	0	DDLValue = new Uint16[ValueCount];
534	0	LODValue = new double[ValueCount];
535	0	if ((DDLValue != NULL) && (LODValue != NULL))
536	0	{
537		/* set x values linearly */
538	0	unsigned int i;
539	0	for (i = 0; i <= MaxDDLValue; ++i)
540	0	DDLValue[i] = OFstatic_cast(Uint16, i);
541		/* compute new y values */
542	0	status = DiCurveFitting<Uint16, double>::calculateValues(0, MaxDDLValue, LODValue,
543	0	OFstatic_cast(unsigned int, ValueCount), Order, coeff);
544	0	}
545	0	}
546	0	delete[] coeff;
547	0	} else {
548		/* use cubic spline interpolation */
549	0	double *spline = new double[ValueCount];
550	0	if ((spline != NULL) &&
551	0	DiCubicSpline<Uint16, double>::Function(DDLValue, LODValue, OFstatic_cast(unsigned int, ValueCount), spline))
552	0	{
553		/* save old values */
554	0	const unsigned long count = ValueCount;
555	0	Uint16 *old_ddl = DDLValue;
556	0	double *old_val = LODValue;
557		/* create new data arrays */
558	0	ValueCount = OFstatic_cast(unsigned long, MaxDDLValue) + 1;
559	0	DDLValue = new Uint16[ValueCount];
560	0	LODValue = new double[ValueCount];
561	0	if ((DDLValue != NULL) && (LODValue != NULL))
562	0	{
563		/* set x values linearly */
564	0	unsigned int i;
565	0	for (i = 0; i <= MaxDDLValue; ++i)
566	0	DDLValue[i] = OFstatic_cast(Uint16, i);
567		/* compute new y values */
568	0	status = DiCubicSpline<Uint16, double>::Interpolation(old_ddl, old_val, spline,
569	0	OFstatic_cast(unsigned int, count), DDLValue, LODValue, OFstatic_cast(unsigned int, ValueCount));
570	0	}
571		/* delete old data arrays */
572	0	delete[] old_ddl;
573	0	delete[] old_val;
574	0	}
575	0	delete[] spline;
576	0	}
577	0	return status;
578	0	}
579	0	return 2;
580	0	}
581
582
583		int DiDisplayFunction::calculateMinMax()
584	0	{
585	0	if ((LODValue != NULL) && (ValueCount > 0))
586	0	{
587	0	MinValue = LODValue[0];
588	0	MaxValue = LODValue[0];
589	0	unsigned long i;
590	0	for (i = 1; i < ValueCount; ++i)
591	0	{
592	0	if (LODValue[i] < MinValue)
593	0	MinValue = LODValue[i];
594	0	if (LODValue[i] > MaxValue)
595	0	MaxValue = LODValue[i];
596	0	}
597	0	return 1;
598	0	}
599	0	return 0;
600	0	}
601
602
603		int DiDisplayFunction::checkMinMaxDensity() const
604	0	{
605	0	if ((MinDensity >= 0) && (MaxDensity >= 0) && (MinDensity >= MaxDensity))
606	0	{
607	0	DCMIMGLE_WARN("invalid optical density range (Dmin = " << MinDensity << ", Dmax = " << MaxDensity << ")");
608	0	return 0;
609	0	}
610	0	return 1;
611	0	}
612
613
614		double DiDisplayFunction::getMinLuminanceValue() const
615	0	{
616		/* Dmax = -1 means unspecified */
617	0	return (MaxDensity < 0) ? -1 : convertODtoLum(MaxDensity);
618	0	}
619
620
621		double DiDisplayFunction::getMaxLuminanceValue() const
622	0	{
623		/* Dmin = -1 means unspecified */
624	0	return (MinDensity < 0) ? -1 : convertODtoLum(MinDensity);
625	0	}
626
627
628		double DiDisplayFunction::convertODtoLumTable(const double od_tab,
629		const unsigned long count,
630		const OFBool useAmb)
631	0	{
632	0	double *lum_tab = NULL;
633	0	if ((od_tab != NULL) && (count > 0))
634	0	{
635		/* create a new table for the luminance values */
636	0	lum_tab = new double[count];
637	0	if (lum_tab != NULL)
638	0	{
639		/* compute luminance values from optical density */
640	0	unsigned int i;
641	0	if (useAmb)
642	0	{
643	0	for (i = 0; i < count; ++i)
644	0	lum_tab[i] = AmbientLight + Illumination * pow(OFstatic_cast(double, 10), -od_tab[i]);
645	0	} else {
646		/* ambient light is added later */
647	0	for (i = 0; i < count; ++i)
648	0	lum_tab[i] = Illumination * pow(OFstatic_cast(double, 10), -od_tab[i]);
649	0	}
650	0	}
651	0	}
652	0	return lum_tab;
653	0	}
654
655
656		double DiDisplayFunction::convertODtoLum(const double value,
657		const OFBool useAmb) const
658	0	{
659	0	return (useAmb) ? convertODtoLum(value, AmbientLight, Illumination) :
660	0	convertODtoLum(value, 0, Illumination);
661	0	}
662
663
664		double DiDisplayFunction::convertODtoLum(const double value,
665		const double ambient,
666		const double illum)
667	0	{
668		/* formula from DICOM PS3.14: L = La + L0 * 10^-D */
669	0	return (value >= 0) && (ambient >= 0) && (illum >= 0) ?
670	0	ambient + illum * pow(OFstatic_cast(double, 10), -value) : -1 /invalid/;
671	0	}