326 lines
10 KiB
C++
326 lines
10 KiB
C++
/*
|
|
*
|
|
* Copyright (C) 1996-2016, 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: DiCurveFitting (header/implementation)
|
|
*
|
|
*/
|
|
|
|
|
|
#ifndef DICRVFIT_H
|
|
#define DICRVFIT_H
|
|
|
|
#include "dcmtk/config/osconfig.h"
|
|
#include "dcmtk/ofstd/oftypes.h"
|
|
#include "dcmtk/ofstd/ofcast.h"
|
|
|
|
#define INCLUDE_CMATH
|
|
#define INCLUDE_CSTDDEF /* For NULL */
|
|
#include "dcmtk/ofstd/ofstdinc.h"
|
|
|
|
|
|
/*---------------------*
|
|
* macro definitions *
|
|
*---------------------*/
|
|
|
|
// SunCC 4.x does not support default values for template types :-/
|
|
#define T3_ double
|
|
|
|
|
|
/*------------------*
|
|
* template class *
|
|
*------------------*/
|
|
|
|
/** Template class for polynomial curve fitting algorithm
|
|
*/
|
|
template <class T1, class T2 /*, class T3 = double*/>
|
|
class DiCurveFitting
|
|
{
|
|
|
|
public:
|
|
|
|
/** calculate coefficients for resulting polynomial function.
|
|
* T1 = type of x coordinates
|
|
* T2 = type of y coordinates
|
|
* T3_ = type of coefficients (and for internal calculations)
|
|
*
|
|
** @param x array with x coordinates of given points
|
|
* @param y array with y coordinates of given points
|
|
* @param n number of entries in array (= points)
|
|
* @param o order of polynomial function
|
|
* @param c array to store the resulting coefficients (o+1 entries !)
|
|
*
|
|
** @return true if successful, false otherwise
|
|
*/
|
|
static int calculateCoefficients(const T1 *x,
|
|
const T2 *y,
|
|
const unsigned int n,
|
|
const unsigned int o,
|
|
T3_ *c)
|
|
{
|
|
int result = 0;
|
|
if ((x != NULL) && (y != NULL) && (c !=NULL) && (n > 0))
|
|
{
|
|
const unsigned int order = o + 1;
|
|
const unsigned int order2 = order * order;
|
|
T3_ *basis = new T3_[order * n];
|
|
T3_ *alpha = new T3_[order2];
|
|
T3_ *beta = new T3_[order];
|
|
if ((basis != NULL) && (alpha != NULL) && (beta != NULL))
|
|
{
|
|
unsigned int i;
|
|
unsigned int j;
|
|
unsigned int k;
|
|
for (i = 0; i < order; ++i)
|
|
{
|
|
for (j = 0; j < n; ++j)
|
|
{
|
|
k = i + j * order;
|
|
if (i == 0)
|
|
basis[k] = 1;
|
|
else
|
|
basis[k] = OFstatic_cast(T3_, x[j]) * basis[k - 1];
|
|
}
|
|
}
|
|
T3_ sum;
|
|
for (i = 0; i < order; ++i)
|
|
{
|
|
const unsigned int i_order = i * order;
|
|
for (j = 0; j <= i; ++j)
|
|
{
|
|
sum = 0;
|
|
for (k = 0; k < n; ++k)
|
|
sum += basis[i + k * order] * basis[j + k * order];
|
|
alpha[i + j * order] = sum;
|
|
if (i != j)
|
|
alpha[j + i_order] = sum;
|
|
}
|
|
}
|
|
for (i = 0; i < order; ++i)
|
|
{
|
|
sum = 0;
|
|
for (j = 0; j < n; ++j)
|
|
sum += OFstatic_cast(T3_, y[j]) * basis[i + j * order];
|
|
beta[i] = sum;
|
|
}
|
|
if (solve(alpha, beta, order))
|
|
{
|
|
for (i = 0; i < order; ++i)
|
|
c[i] = beta[i];
|
|
result = 1;
|
|
}
|
|
}
|
|
delete[] basis;
|
|
delete[] alpha;
|
|
delete[] beta;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
|
|
/** calculate y coordinates for the given range of x coordinates.
|
|
* The polynomial function is defined by the specified coefficients.
|
|
* T1 = type of x coordinates
|
|
* T2 = type of y coordinates
|
|
* T3_ = type of coefficients (and for internal calculations)
|
|
*
|
|
** @param xs first x coordinate for computation
|
|
* @param xe last x coordinate for computation
|
|
* @param y array to store the resulting y coordinates (n entries !)
|
|
* @param n number of entries in array (= points)
|
|
* @param o order of polynomial function
|
|
* @param c array of coefficients computed by the above method (o+1 entries !)
|
|
*
|
|
** @return true if successful, false otherwise
|
|
*/
|
|
static int calculateValues(const T1 xs,
|
|
const T1 xe,
|
|
T2 *y,
|
|
const unsigned int n,
|
|
const unsigned int o,
|
|
const T3_ *c)
|
|
{
|
|
int result = 0;
|
|
if ((y != NULL) && (c != NULL) && (n > 0) && (xe > xs))
|
|
{
|
|
unsigned int i;
|
|
unsigned int j;
|
|
T3_ x;
|
|
T3_ x2;
|
|
T3_ w;
|
|
const T3_ xo = OFstatic_cast(T3_, xs);
|
|
const T3_ xi = OFstatic_cast(T3_, (OFstatic_cast(T3_, xe) - OFstatic_cast(T3_, xs)) / (n - 1));
|
|
for (i = 0; i < n; ++i)
|
|
{
|
|
x = xo + OFstatic_cast(T3_, i) * xi;
|
|
x2 = 1;
|
|
w = 0;
|
|
for (j = 0; j <= o; ++j)
|
|
{
|
|
w += c[j] * x2;
|
|
x2 *= x;
|
|
}
|
|
convertValue(w, y[i]); // cut value if necessary
|
|
}
|
|
result = 1;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
|
|
private:
|
|
|
|
/** helper routine: convert to unsigned 8 bit value
|
|
*
|
|
** @param input input value to be converted
|
|
* @param output output value (range: 0..255)
|
|
*
|
|
** @return output value
|
|
*/
|
|
static void convertValue(const T3_ input, Uint8 &output)
|
|
{
|
|
output = (input <= 0) ? 0 : ((input >= 255) ? 255 : OFstatic_cast(Uint8, input));
|
|
}
|
|
|
|
/** helper routine: convert to signed 8 bit value
|
|
*
|
|
** @param input input value to be converted
|
|
* @param output output value (range: -128..127)
|
|
*
|
|
** @return output value
|
|
*/
|
|
static void convertValue(const T3_ input, Sint8 &output)
|
|
{
|
|
output = (input <= -128) ? -128 : ((input >= 127) ? 127 : OFstatic_cast(Sint8, input));
|
|
}
|
|
|
|
/** helper routine: convert to unsigned 16 bit value
|
|
*
|
|
** @param input input value to be converted
|
|
* @param output output value (range: 0..65535)
|
|
*
|
|
** @return output value
|
|
*/
|
|
static void convertValue(const T3_ input, Uint16 &output)
|
|
{
|
|
output = (input <= 0) ? 0 : ((input >= 65535) ? 65535 : OFstatic_cast(Uint16, input));
|
|
}
|
|
|
|
/** helper routine: convert to signed 16 bit value
|
|
*
|
|
** @param input input value to be converted
|
|
* @param output output value (range: -32768..32767)
|
|
*
|
|
** @return output value
|
|
*/
|
|
static void convertValue(const T3_ input, Sint16 &output)
|
|
{
|
|
output = (input <= -32768) ? -32768 : ((input >= 32767) ? 32767 : OFstatic_cast(Sint16, input));
|
|
}
|
|
|
|
/** helper routine: convert to floating point value (double precision)
|
|
*
|
|
** @param input input value to be converted
|
|
* @param output output value (double)
|
|
*
|
|
** @return output value
|
|
*/
|
|
static inline void convertValue(const T3_ input, double &output)
|
|
{
|
|
output = OFstatic_cast(double, input);
|
|
}
|
|
|
|
/** solve the equation given by the two matrixes.
|
|
* T3_ = type of coefficients (and for internal calculations)
|
|
*
|
|
** @param a first matrix (array of values)
|
|
* @param b second matrix (array of values)
|
|
* @param n number of entries in array
|
|
*
|
|
** @return true if successful, false otherwise
|
|
*/
|
|
static int solve(T3_ *a,
|
|
T3_ *b,
|
|
const unsigned int n)
|
|
{
|
|
int result = 0;
|
|
if ((a != NULL) && (b != NULL) && (n > 0))
|
|
{
|
|
unsigned int i;
|
|
unsigned int j;
|
|
unsigned int k;
|
|
signed int pivot;
|
|
T3_ mag;
|
|
T3_ mag2;
|
|
T3_ temp;
|
|
for (i = 0; i < n; ++i)
|
|
{
|
|
mag = 0;
|
|
pivot = -1;
|
|
for (j = i; j < n; ++j)
|
|
{
|
|
mag2 = fabs(a[i + j * n]);
|
|
if (mag2 > mag)
|
|
{
|
|
mag = mag2;
|
|
pivot = j;
|
|
}
|
|
}
|
|
if ((pivot == -1) || (mag == 0))
|
|
break;
|
|
else
|
|
{
|
|
const unsigned int piv = OFstatic_cast(unsigned int, pivot);
|
|
const unsigned int i_n = i * n;
|
|
if (piv != i)
|
|
{
|
|
const unsigned int piv_n = piv * n;
|
|
for (j = i; j < n; ++j)
|
|
{
|
|
temp = a[j + i_n];
|
|
a[j + i_n] = a[j + piv_n];
|
|
a[j + piv_n] = temp;
|
|
}
|
|
temp = b[i];
|
|
b[i] = b[piv];
|
|
b[piv] = temp;
|
|
}
|
|
mag = a[i + i_n];
|
|
for (j = i; j < n; ++j)
|
|
a[j + i_n] /= mag;
|
|
b[i] /= mag;
|
|
for (j = 0; j < n; ++j)
|
|
{
|
|
if (i == j)
|
|
continue;
|
|
const unsigned int j_n = j * n;
|
|
mag2 = a[i + j_n];
|
|
for (k = i; k < n; ++k)
|
|
a[k + j_n] -= mag2 * a[k + i_n];
|
|
b[j] -= mag2 * b[i];
|
|
}
|
|
result = 1;
|
|
}
|
|
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
};
|
|
|
|
|
|
#endif
|