// template library for multidimensional image processing and analysis
// Copyright (C) 2013-2023 Jiri Janacek
// 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.
#if !defined(jil_included)
#define jil_included
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
#include
#include
#include
#include
#include
#include
#include
#include
#include
#ifdef min
#undef min
#define jil_redefine_min
#endif
#ifdef max
#undef max
#define jil_redefine_max
#endif
//assuming sizeof(size_t)==NumB
//assuming sizeof(ptrdiff_t)==NumB
//assuming Numb==8*NumB
#ifdef _WIN64
#define Numb 64
#define NumB 8
#else
#define Numb 32
#define NumB 4
#endif
const float max_float= std::numeric_limits::max();
const float max_double= std::numeric_limits::max();
//ConstantsXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
//operations
enum UnOpGray {IL_GU_SET, IL_GU_PLUS, IL_GU_MINUS, IL_GU_SUB,
IL_GU_MULT, IL_GU_DIV, IL_GU_MAX, IL_GU_MIN, IL_GU_NOT, IL_GU_SQR, IL_GU_POW};
enum UnOpFloat {IL_FU_MULT, IL_FU_DIV, IL_FU_POW};
enum BinOpGray {IL_GB_PLUS, IL_GB_MINUS, IL_GB_SUB, IL_GB_MAX,
IL_GB_MULT, IL_GB_DIV, IL_GB_MIN};
enum UnOpBin {IL_BU_SET, IL_BU_NOT};
enum BinOpBin {IL_BB_OR, IL_BB_AND, IL_BB_XOR, IL_BB_DIF, IL_BB_NIF, IL_BB_XORN};
enum UnOpObj {IL_OU_SET, IL_OU_NOT, IL_OU_SQR};
//connectivity
enum con2D {con2= 0, con4 = 1, con8};
static inline con2D complement(con2D con) {
if (con!=con2) return con2D(3-con);
else return con;
}
enum con3D {con3= 0, con6 = 1, con18, con26};
static inline con3D complement(con3D con) {
if ((con!=con3)&&(con!=con18)) return con3D(4-con);
else return con;
}
//resampling method
enum resmet {resnear= 0, reslin = 1, rescub = 2};
//thresholding mode
enum thrmod {thrgt, thrge, thrlt, thrle, thre};
//reconstruction mode
enum recmod {recpart, rectot};
//binary mode
enum binmod {binpos, binneg};
//neigbourhood mode
enum ngbhsiz {ngbh2= 0, ngbh3};
//math
#ifndef SQR
static double dsqrarg;
#define SQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
#endif
#ifndef M_PI
const double M_PI = 3.14159265358979323846;
#endif
#ifndef M_PI_2
const double M_PI_2 = 1.57079632679489661923;
#endif
#ifndef M_SQRT2
const double M_SQRT2 = 1.41421356237309504880;
#endif
#ifndef M_SQRT2_2
const double M_SQRT2_2 = 0.70710678118654752440;
#endif
//errors
class Imgerr { }; //general
class Allocerr : public Imgerr { }; //data not allocated
class Dimserr : public Imgerr { }; //dimensions do not match
class NVarerr : public Imgerr { }; //number of variables do not match
class Sizeerr : public Imgerr { }; //sizes do not match
class Argerr : public Imgerr { }; //bad argument
class Memerr : public Imgerr { }; //memory alloc
class Matherr : public Imgerr { }; //bad math operation
class Valerr : public Imgerr { }; //too many labels
class Numerr : public Imgerr { }; //too few objects
//Utilities****************************************************************************
//binary constants
const size_t LB = (size_t)(1)<<(Numb-1);
const size_t RB = 1;
const size_t AB = (size_t)(-1);
const ptrdiff_t narr[0x100]= {
0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8};
static inline ptrdiff_t bitnum(size_t x) {
return narr[0xff & x] + narr[0xff & (x>>8)]
+ narr[0xff & (x>>16)] + narr[0xff & (x>>24)]
#ifdef _WIN64
+ narr[0xff & (x>>32)] + narr[0xff & (x>>40)]
+ narr[0xff & (x>>48)] + narr[0xff & (x>>56)]
#endif
;
}
inline ptrdiff_t DimsProd (int nDim0, ptrdiff_t* dims0, int nVars= 1) {
ptrdiff_t i, p;
for (i= 0, p= 1; i < nDim0; i++)
p*= dims0[i];
return p * nVars;
}
inline ptrdiff_t DimsMax (int nDim0, ptrdiff_t* dims0) {
ptrdiff_t p;
int i;
for (i= 1, p= *dims0; i < nDim0; i++)
if (dims0[i] > p) p= dims0[i];
return p;
}
inline ptrdiff_t DimsProdEx (int nDim0, ptrdiff_t* dims0, int except) {
ptrdiff_t p;
int i;
for (i= 0, p= 1; i < nDim0; i++)
if (i != except) p*= dims0[i];
return p;
}
inline bool DimsEq (int nDim0, ptrdiff_t* dims0, ptrdiff_t* dims1) {
int i;
bool p= true;
for (i= 0; i < nDim0; i++)
p&= (dims0[i] == dims1[i]);
return p;
}
inline bool DimsEqEx (int nDim0, ptrdiff_t* dims0, ptrdiff_t* dims1, int except) {
int i;
bool p= true;
for (i= 0; i < nDim0; i++)
p&= (i == except) || (dims0[i] == dims1[i]);
return p;
}
inline bool DimsEqProj(int nDim0, ptrdiff_t* dims0, ptrdiff_t* dims1, int dimension) {
int i;
int j= 0;
bool p= true;
for (i= 0; i < nDim0; i++)
if (i != dimension)
p&= (dims0[j++] == dims1[i]);
return p;
}
inline void SetBit(size_t *data, ptrdiff_t index) {
ptrdiff_t poloha= (index / Numb);
size_t k= LB;
k>>= (index % Numb);
data[poloha]|= k;
}
inline bool GetBit(size_t *data, ptrdiff_t index) {
ptrdiff_t poloha= (index / Numb);
size_t k= LB;
k>>= (index % Numb);
return ((data[poloha] & k)!= 0);
}
//3d vectors utilities
inline void substvect3d(double x, double y, double z, double *v) {
v[0]= x; v[1]= y; v[2]= z;
}
inline void conv2vect3d(double *x1, double *x2, double a, double *x) {
int j;
for (j= 0; j < 3; j++)
x[j]= a * x1[j] + (1.0 - a) * x2[j];
}
inline void copyvect3d(double *u, double *v) {
int j;
for (j= 0; j < 3; j++)
v[j]= u[j];
}
inline void normalize3d(double *x) {
int j;
double s;
for (j= 0, s= 0.0; j < 3; j++) s+= x[j] * x[j];
s= sqrt(s);
for (j= 0; j < 3; j++) x[j]/= s;
}
inline void subperp3d(double *x, double *y) {
int j;
double s;
for (j= 0, s= 0.0; j < 3; j++) s+= x[j] * y[j];
for (j= 0; j < 3; j++) x[j]-= s * y[j];
}
inline void comb2vect3d(double *x1, double *x2, double a, double b) {
int j;
for (j= 0; j < 3; j++)
x1[j]= a * x1[j] + b * x2[j];
}
inline void sub3d(double *x, double *y) {
int j;
for (j= 0; j < 3; j++) x[j]-= y[j];
}
inline void mul3d(double *x, double *y) {
int j;
for (j= 0; j < 3; j++) x[j]*= y[j];
}
inline void vec3d(double *x, double *y) {
int j;
double s;
double value[3];
for (j= 0, s= 0.0; j < 3; j++) value[j]= x[j];
x[0]= value[1] * y[2] - value[2] * y[1];
x[1]= value[2] * y[0] - value[0] * y[2];
x[2]= value[0] * y[1] - value[1] * y[0];
}
inline int GetIcosaNum(int divs) {
return 5 * divs * divs + 1;
}
template bool GenIcosaSyst(int divs, double **&y) {
const double eps= 1e-6;
int num_nodes, vrch;
double dp, dq;
int i, j, k, p1, r1, s1, p2, r2, p3, r3,
p21, r21, p20, r20, s20, t1, t2, t3;
double **x;
int connec0[60]=
{ 0, 1, 2, 0, 2, 3, 0, 3, 4, 0, 4, 5,
0, 5, 1, 1, 6, 2, 2, 6, 7, 2, 7, 3,
3, 7, 8, 3, 8, 4, 4, 8, 9, 4, 9, 5,
5, 9,10, 5,10, 1, 1,10, 6, 6,11, 7,
7,11, 8, 8,11, 9, 9,11,10,10,11, 6};
int edges[60]=
{ 0, 1, 0, 2, 0, 3, 0, 4, 0, 5,
1, 2, 2, 3, 3, 4, 4, 5, 5, 1,
1,10, 1, 6, 2, 6, 2, 7, 3, 7,
3, 8, 4, 8, 4, 9, 5, 9, 5,10,
6, 7, 7, 8, 8, 9, 9,10,10, 6,
6,11, 7,11, 8,11, 9,11,10,11};
num_nodes= 10 * divs * divs + 2;
typedef double *pdouble;
x= new pdouble[num_nodes];
if (!x) return false;
x[0]= new double[num_nodes * 3];
if (!x[0]) {
delete[] x;
return false;
}
for(i= 1; i < num_nodes; i++) x[i]= x[i-1]+3;
//vertices
substvect3d(0.0, 0.0, 1.0, x[0]);
if (divs <= 0)
substvect3d(0.0, 0.0, -1.0, x[1]);
else {
substvect3d(0.0, 0.0, -1.0, x[11]);
dp= 1.0 / sqrt(5.0); dq= sqrt(1.0 - dp * dp);
for (i= 0; i < 5; i++) {
substvect3d(dq * cos(2.0 * M_PI * i / 5.0),
dq * sin(2.0 * M_PI * i / 5.0),
dp, x[1 + i]);
substvect3d(dq * cos(2.0 * M_PI * (i + 0.5) / 5.0),
dq * sin(2.0 * M_PI * (i + 0.5) / 5.0),
-dp, x[6 + i]);
}
//edges
if (divs > 1) for (i= 0; i < 30; i++) {
vrch= 11 + i * (divs - 1);
for (j= 1; j < divs; j++)
conv2vect3d(x[edges[2 * i + 1]], x[edges[2 * i]], (float)j/(float)divs, x[vrch + j]);
}
vrch= 12 + 30 * (divs - 1);
//planes
if (divs > 2) for (i= 0; i < 20; i++) {
t1= connec0[3 * i];
t2= connec0[3 * i + 1];
t3= connec0[3 * i + 2];
j= 0;
while (!((edges[2 * j] == t1)&&(edges[2 * j + 1] == t2))
&&!((edges[2 * j] == t2)&&(edges[2 * j + 1] == t1))) j++;
p1= j;
if (edges[2 * p1] == t1) {
p3= 1; p2= 12 + p1 * (divs - 1);
}
else {
p3= -1; p2= 12 + (p1 + 1) * (divs - 1) - 1;
}
j= 0;
while (!((edges[2 * j] == t1)&&(edges[2 * j + 1] == t3))
&&!((edges[2 * j] == t3)&&(edges[2 * j + 1] == t1))) j++;
r1= j;
if (edges[2 * r1] == t1) {
r3= 1; r2= 12 + r1 * (divs - 1);
}
else {
r3= -1; r2= 12 + (r1 + 1) * (divs - 1) - 1;
}
j= 0;
while (!((edges[2 * j] == t2)&&(edges[2 * j + 1] == t3))
&&!((edges[2 * j] == t3)&&(edges[2 * j + 1] == t2))) j++;
s1= j;
p20= p2; r20= r2; s20= vrch;
for (j= 2; j < divs; j++) {
p21= p2;
r21= r2;
p2+= p3;
r2+= r3;
conv2vect3d(x[r2], x[p2], 1.0/(float)j, x[vrch]);
vrch++;
for (k= 2; k < j; k++) {
conv2vect3d(x[r2], x[p2], (float)k/(float)j, x[vrch]);
vrch++;
}
}
}
}
//normalize
for (i= 0; i < num_nodes; i++) normalize3d(x[i]);
int n_nodes= num_nodes / 2;
y= new pdouble[n_nodes];
if (!y) return false;
y[0]= new double[n_nodes * 3];
if (!y[0]) {
delete[] y;
return false;
}
for(i= 1; i < n_nodes; i++) y[i]= y[i-1]+3;
//symmetrize
for (i= 0, j= 0; i < num_nodes; i++)
if ((x[i][2]>eps)||((fabs(x[i][2])<=eps)&&(x[i][1]>eps))
||((fabs(x[i][2])<=eps)&&(fabs(x[i][1])<=eps)&&(x[i][0]>0.0))) {
if (j >= (num_nodes / 2)) {
delete[] x[0];
delete[] x;
return false;
}
copyvect3d(x[i],y[j]);
j++;
}
//clean up
delete[] x[0];
delete[] x;
return true;
}
//ClassesXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
class BoundBox {
public:
double* bounds;
int nDim;
BoundBox () : bounds (0), nDim (0) { }
explicit BoundBox (int nDim0) : bounds (0), nDim (nDim0) { }
BoundBox (const BoundBox& x0) {
nDim= x0.nDim;
bounds= new double[2 * nDim];
SetBounds(x0.bounds);
}
BoundBox (int nDim0, double* bounds0) {
nDim= nDim0;
bounds= new double[2 * nDim];
SetBounds(bounds0);
}
void SetBounds (double* bounds0) {
if (!bounds)
bounds= new double[2 * nDim];
int i;
for (i= 0; i < 2 * nDim; i++) bounds[i]= bounds0[i];
}
void SetBounds (float* bounds0) {
if (!bounds)
bounds= new double[2 * nDim];
int i;
for (i= 0; i < 2 * nDim; i++) bounds[i]= bounds0[i];
}
double Size (int ndim) const {
if (!bounds)
return 0.0;
return bounds[2 * ndim + 1] - bounds[2 * ndim];
}
double Diagonal () const {
if (!bounds)
return 0.0;
int i;
double p;
for (i= 0, p= 0.0; i < nDim; i++) p+= SQR(bounds[2 * i + 1] - bounds[2 * i]);
return sqrt(p);
}
double Volume () const {
if (!bounds)
return 0.0;
int i;
double p;
for (i= 0, p= 1.0; i < nDim; i++) p*= bounds[2 * i + 1] - bounds[2 * i];
return p;
}
void Shift(int ndim, double by) {
bounds[2 * ndim]+= by;
bounds[2 * ndim + 1]+= by;
}
void UCrop(int ndim, double by) {
bounds[2 * ndim + 1]-= by;
}
void LCrop(int ndim, double by) {
bounds[2 * ndim]+= by;
}
void Crop(int ndim, double by) {
bounds[2 * ndim]+= by;
bounds[2 * ndim + 1]-= by;
}
~BoundBox () {
if (bounds) delete[] bounds;
}
};
class Frame : public BoundBox {
public:
ptrdiff_t* dims;
Frame () : BoundBox(), dims (0) { }
Frame (const Frame& x0) : BoundBox(x0) {
dims= new ptrdiff_t[nDim];
SetDims(x0.dims);
}
Frame (int nDim0, ptrdiff_t* dims0 = 0, double* bounds0 = 0) {
nDim= nDim0;
bounds= new double[2 * nDim];
dims= new ptrdiff_t[nDim];
int i;
if (dims0) SetDims(dims0);
else for (i= 0; i < nDim; i++)
dims[i]= 1;
if (bounds0) SetBounds(bounds0);
else for (i= 0; i < nDim; i++) {
bounds[2 * i]= 0;
bounds[2 * i + 1]= dims[i];
}
}
double Dx(int i) const {
return (bounds[2 * i + 1] - bounds[2 * i]) / dims[i];
}
void SetDims (ptrdiff_t* dims0) {
if (!dims)
dims= new ptrdiff_t[nDim];
int i;
for (i= 0; i < nDim; i++) dims[i]= dims0[i];
}
void ResetBounds() {
int i;
for (i= 0; i < nDim; i++) {
bounds[2 * i]= 0;
bounds[2 * i + 1]= dims[i];
}
}
~Frame () {
if (dims) delete[] dims;
}
};
template class Image : public Frame {
public:
F* data;
int nVars;
Image () : Frame (), data (0), nVars(1) { }
Image (Image& x0) : Frame (x0), data (0), nVars(x0.nVars) { } //overrides default!
Image (const Frame& x0, int nVars0 = 1) : Frame (x0), data (0), nVars(nVars0) { }
Image (int nDim0, ptrdiff_t* dims0, double* bounds0 = 0, int nVars0 = 1) : Frame(nDim0, dims0, bounds0), data (0), nVars (nVars0) { }
void NewData () {
if (!data) {
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
data= new F[ndat];
if (!data) throw Memerr();
}
}
void DeleteData () {
if (data) delete[] data;
data= 0;
}
templatevoid SetData (G* data0, int step= 1) {
if (!data) NewData();
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
F *p= data;
G *q= data0;
ptrdiff_t i;
for (i= 0; i < ndat; i++, p++, q+= step) *p= *q;
}
templatevoid SetData (G* data0, int step= 1, H mult= 1) {
if (!data) NewData();
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
F *p= data;
G *q= data0;
ptrdiff_t i;
for (i= 0; i < ndat; i++, p++, q+= step) *p= *q * mult;
}
templatevoid SetDataChannel (int channel, G* data0, int step= 1) {
if ((channel < 0)||(channel >= nVars)) throw Argerr();
if (!data) NewData();
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + channel;
G *q= data0;
ptrdiff_t i;
for (i= 0; i < ndat; i++, p+= nVars, q+= step) *p= *q;
}
templatevoid GetData (G* data0, int step= 1) {
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
F *p= data;
G *q= data0;
ptrdiff_t i;
for (i= 0; i < ndat; i++, q+= step, p++) *q= *p;
}
templatevoid GetDataChannel (int channel, G* data0, int step= 1) {
if ((channel < 0)||(channel >= nVars)) throw Argerr();
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + channel;
G *q= data0;
ptrdiff_t i;
for (i= 0; i < ndat; i++, q+= step, p+= nVars) *q= *p;
}
void SetData (UnOpGray op = IL_GU_SET, F val = 0) {
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
F* p= data;
ptrdiff_t i;
int issigned= std::numeric_limits::is_signed;
switch (op) {
case IL_GU_SET:
for (i= 0; i < ndat; i++, p++)
*p= val;
break;
case IL_GU_PLUS:
for (i= 0; i < ndat; i++, p++)
*p+= val;
break;
case IL_GU_MINUS: {
if (issigned)
for (i= 0; i < ndat; i++, p++)
*p-= val;
else
for (i= 0; i < ndat; i++, p++)
if (*p > val) *p-= val;
else *p= 0;
}
break;
case IL_GU_SUB: {
if (issigned)
for (i= 0; i < ndat; i++, p++)
*p= val - *p;
else
for (i= 0; i < ndat; i++, p++)
if (*p < val) *p= val - *p;
else *p= 0;
}
break;
case IL_GU_MULT:
for (i= 0; i < ndat; i++, p++)
*p*= val;
break;
case IL_GU_DIV:
for (i= 0; i < ndat; i++, p++)
*p/= val;
break;
case IL_GU_MIN:
for (i= 0; i < ndat; i++, p++)
if (*p > val) *p= val;
break;
case IL_GU_MAX:
for (i= 0; i < ndat; i++, p++)
if (*p < val) *p= val;
break;
case IL_GU_NOT: {
F maxF= (std::numeric_limits::is_integer)
?(std::numeric_limits::max()):(F)0;
for (i= 0; i < ndat; i++, p++)
*p= maxF - *p;
}
break;
case IL_GU_SQR:
for (i= 0; i < ndat; i++, p++)
*p= SQR(*p);
break;
case IL_GU_POW:
for (i= 0; i < ndat; i++, p++)
*p= pow(double(*p),double(val));
break;
}
}
void SetDataChannel (int channel, UnOpGray op = IL_GU_SET, F val = 0) {
if ((channel < 0)||(channel >= nVars)) throw Argerr();
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + channel;
ptrdiff_t i;
int issigned= std::numeric_limits::is_signed;
switch (op) {
case IL_GU_SET:
for (i= 0; i < ndat; i++, p+= nVars)
*p= val;
break;
case IL_GU_PLUS:
for (i= 0; i < ndat; i++, p+= nVars)
*p+= val;
break;
case IL_GU_MINUS: {
if (issigned)
for (i= 0; i < ndat; i++, p+= nVars)
*p-= val;
else
for (i= 0; i < ndat; i++, p+= nVars)
if (*p > val) *p-= val;
else *p= 0;
}
break;
case IL_GU_SUB: {
if (issigned)
for (i= 0; i < ndat; i++, p+= nVars)
*p= val - *p;
else
for (i= 0; i < ndat; i++, p+= nVars)
if (*p < val) *p= val - *p;
else *p= 0;
}
break;
case IL_GU_MULT:
for (i= 0; i < ndat; i++, p+= nVars)
*p*= val;
break;
case IL_GU_DIV:
for (i= 0; i < ndat; i++, p+= nVars)
*p/= val;
break;
case IL_GU_MIN:
for (i= 0; i < ndat; i++, p+= nVars)
if (*p > val) *p= val;
break;
case IL_GU_MAX:
for (i= 0; i < ndat; i++, p+= nVars)
if (*p < val) *p= val;
break;
case IL_GU_NOT: {
F maxF= (std::numeric_limits::is_integer)
?(std::numeric_limits::max()):(F)0;
for (i= 0; i < ndat; i++, p+= nVars)
*p= maxF - *p;
}
break;
case IL_GU_SQR:
for (i= 0; i < ndat; i++, p+= nVars)
*p= SQR(*p);
break;
case IL_GU_POW:
for (i= 0; i < ndat; i++, p+= nVars)
*p= pow(double(*p),double(val));
break;
}
}
void SetData (UnOpFloat op, double val = 0) {
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
F* p= data;
ptrdiff_t i;
int testmax= std::numeric_limits::is_integer;
F fmax= std::numeric_limits::max();
switch (op) {
case IL_FU_MULT:
for (i= 0; i < ndat; i++, p++) {
double pf= *p * val;
if (testmax && (pf > fmax)) *p= fmax;
else *p= pf;
}
break;
case IL_FU_DIV:
for (i= 0; i < ndat; i++, p++) {
double pf= *p / val;
if (testmax && (pf > fmax)) *p= fmax;
else *p= pf;
}
break;
case IL_FU_POW:
for (i= 0; i < ndat; i++, p++) {
double pf= pow(double(*p),val);
if (testmax && (pf > fmax)) *p= fmax;
else *p= pf;
}
break;
}
}
void SetDataChannel (int channel, UnOpFloat op, double val = 0) {
if ((channel < 0)||(channel >= nVars)) throw Argerr();
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + channel;
ptrdiff_t i;
int testmax= std::numeric_limits::is_integer;
F fmax= std::numeric_limits::max();
switch (op) {
case IL_FU_MULT:
for (i= 0; i < ndat; i++, p+= nVars) {
double pf= *p * val;
if (testmax && (pf > fmax)) *p= fmax;
else *p= pf;
}
break;
case IL_FU_DIV:
for (i= 0; i < ndat; i++, p+= nVars) {
double pf= *p / val;
if (testmax && (pf > fmax)) *p= fmax;
else *p= pf;
}
break;
case IL_FU_POW:
for (i= 0; i < ndat; i++, p+= nVars) {
double pf= pow(double(*p),val);
if (testmax && (pf > fmax)) *p= fmax;
else *p= pf;
}
break;
}
}
templatevoid SetData (BinOpGray op, G* data0, int step= 1) {
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
F *p= data;
G *q= data0;
ptrdiff_t i;
int issigned=std::numeric_limits::is_signed;
switch (op) {
case IL_GB_PLUS:
for (i= 0; i < ndat; i++, p++, q+= step) *p+= *q;
break;
case IL_GB_MINUS: {
if (issigned)
for (i= 0; i < ndat; i++, p++, q+= step)
*p-= *q;
else
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p > *q) *p-= *q;
else *p= 0;
}
break;
case IL_GB_SUB: {
if (issigned)
for (i= 0; i < ndat; i++, p++, q+= step)
*p= *q - *p;
else
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p < *q) *p= *q - *p;
else *p= 0;
}
break;
case IL_GB_MULT:
for (i= 0; i < ndat; i++, p++, q+= step) *p*= *q;
break;
case IL_GB_DIV:
for (i= 0; i < ndat; i++, p++, q+= step) *p/= *q;
break;
case IL_GB_MAX:
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p < *q) *p= *q;
break;
case IL_GB_MIN:
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p > *q) *p= *q;
break;
}
}
templatevoid SetData (BinOpGray op, G* data0, int step= 1, H mult= 1) {
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, nVars);
F *p= data;
G *q= data0;
ptrdiff_t i;
int issigned=std::numeric_limits::is_signed;
switch (op) {
case IL_GB_PLUS:
for (i= 0; i < ndat; i++, p++, q+= step) *p+= *q * mult;
break;
case IL_GB_MINUS: {
if (issigned)
for (i= 0; i < ndat; i++, p++, q+= step)
*p-= *q * mult;
else
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p > (*q * mult)) *p-= *q * mult;
else *p= 0;
}
break;
case IL_GB_SUB: {
if (issigned)
for (i= 0; i < ndat; i++, p++, q+= step)
*p= *q * mult - *p;
else
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p < (*q * mult)) *p= *q * mult - *p;
else *p= 0;
}
break;
case IL_GB_MULT:
for (i= 0; i < ndat; i++, p++, q+= step) *p*= *q * mult;
break;
case IL_GB_DIV:
for (i= 0; i < ndat; i++, p++, q+= step) *p/= *q * mult;
break;
case IL_GB_MAX:
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p < (*q * mult)) *p= *q * mult;
break;
case IL_GB_MIN:
for (i= 0; i < ndat; i++, p++, q+= step)
if (*p > (*q * mult)) *p= *q * mult;
break;
}
}
templatevoid SetDataChannel (int channel, BinOpGray op, G* data0, int step= 1) {
if ((channel < 0)||(channel >= nVars)) throw Argerr();
if (!data) throw Allocerr();
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + channel;
G *q= data0;
ptrdiff_t i;
int issigned=std::numeric_limits::is_signed;
switch (op) {
case IL_GB_PLUS:
for (i= 0; i < ndat; i++, p+= nVars, q+= step) *p+= *q;
break;
case IL_GB_MINUS: {
if (issigned)
for (i= 0; i < ndat; i++, p+= nVars, q+= step)
*p-= *q;
else
for (i= 0; i < ndat; i++, p+= nVars, q+= step)
if (*p > *q) *p-= *q;
else *p= 0;
}
break;
case IL_GB_SUB: {
if (issigned)
for (i= 0; i < ndat; i++, p+= nVars, q+= step)
*p= *q - *p;
else
for (i= 0; i < ndat; i++, p+= nVars, q+= step)
if (*p < *q) *p= *q - *p;
else *p= 0;
}
break;
case IL_GB_MULT:
for (i= 0; i < ndat; i++, p+= nVars, q+= step) *p*= *q;
break;
case IL_GB_DIV:
for (i= 0; i < ndat; i++, p+= nVars, q+= step) *p/= *q;
break;
case IL_GB_MAX:
for (i= 0; i < ndat; i++, p+= nVars, q+= step)
if (*p < *q) *p= *q;
break;
case IL_GB_MIN:
for (i= 0; i < ndat; i++, p+= nVars, q+= step)
if (*p > *q) *p= *q;
break;
}
}
double GetAvg(int chan= 0) {
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + chan;
double s;
ptrdiff_t i;
for (i= 0, s= 0.0; i < ndat; i++, p+= nVars) s+= *p;
return s / ndat;
}
double GetSum(int chan= 0) {
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + chan;
double s;
ptrdiff_t i;
for (i= 0, s= 0.0; i < ndat; i++, p+= nVars) s+= *p;
return s;
}
double GetSumSq(int chan= 0) {
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + chan;
double s;
ptrdiff_t i;
for (i= 0, s= 0.0; i < ndat; i++, p+= nVars) s+= SQR(*p);
return s;
}
double GetMax(int chan= 0) {
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + chan;
double maxval= std::numeric_limits::lowest() ;
ptrdiff_t i;
for (i= 0; i < ndat; i++, p+= nVars)
if (*p > maxval) maxval= *p;
return maxval;
}
double GetMin(int chan= 0) {
ptrdiff_t ndat= DimsProd(nDim, dims, 1);
F *p= data + chan;
double minval= max_double ;
ptrdiff_t i;
for (i= 0; i < ndat; i++, p+= nVars)
if (*p < minval) minval= *p;
return minval;
}
};
template class Img : public Image {
public:
Img () : Image () { }
Img (const Img& x0) : Image (x0) { NewData (); } //overrides default!
Img (Image& x0) : Image (x0) { NewData (); }
Img (const Frame& x0, int nVars0 = 1) : Image (x0, nVars0) { NewData (); }
Img (int nDim0, ptrdiff_t* dims0, double* bounds0 = 0, int nVars0 = 1) : Image (nDim0, dims0, bounds0, nVars0) {
NewData ();
}
~Img () {
if (data) delete[] data;
}
};
class Binary : public Frame {
public:
size_t* data;
Binary () : Frame (), data (0) { }
Binary (Binary& x0) : Frame (x0), data (0) { } //overrides default!
Binary (const Frame& x0) : Frame (x0), data (0) { }
Binary (int nDim0, ptrdiff_t* dims0, double* bounds0 = 0) : Frame(nDim0, dims0, bounds0), data (0) { }
void NewData () {
if (!data) {
ptrdiff_t ndat= (dims[0] + Numb - 1) / Numb * DimsProd(nDim - 1, dims + 1);
data= new size_t[ndat];
if (!data) throw Memerr();
}
}
void DeleteData () {
if (data) delete[] data;
data= 0;
}
void EraseEnd () {
ptrdiff_t i, ints, height;
size_t c, *p;
if(data && (dims[0] % Numb)) {
c= AB;
c>>= (dims[0] % Numb);
c= ~c;
ints= ((dims[0] + Numb - 1) / Numb);
height= DimsProd(nDim - 1, dims + 1);
p= data + ints - 1;
for(i= 0; i < height; i++, p+= ints) *p&= c;
}
}
int GetEmpty () {
ptrdiff_t i, j, ints, height, full;
size_t c, *p;
full= 0;
if(data) {
c= AB;
c>>= ((dims[0] - 1) % Numb) + 1;
c= ~c;
ints= ((dims[0] + Numb - 1) / Numb);
height= DimsProd(nDim - 1, dims + 1);
for (j= 0, p= data; j < height; j++) {
for(i= 1; i < ints; i++, p++) full= full || *p;
full= full || (*p & c); p++;
}
}
return !full;
}
ptrdiff_t GetCount() {
ptrdiff_t i, j, ints, height, count;
size_t c, *p;
count= 0;
if(data) {
c= AB;
c>>= ((dims[0] - 1) % Numb) + 1;
c= ~c;
ints= ((dims[0] + Numb - 1) / Numb);
height= DimsProd(nDim - 1, dims + 1);
for (j= 0, p= data; j < height; j++) {
for(i= 1; i < ints; i++, p++) count+= bitnum(*p);
count+= bitnum(*p & c); p++;
}
}
return count;
}
void GetData (size_t* data0) const {
if (!data) return;
ptrdiff_t ndat= (dims[0] + Numb - 1) / Numb * DimsProd(nDim - 1, dims + 1);
size_t *p, *q;
p= data; q= data0;
ptrdiff_t i;
for (i= 0; i < ndat; i++, p++, q++) *q= *p;
}
void SetData (size_t* data0) {
if (!data) NewData();
ptrdiff_t ndat= (dims[0] + Numb - 1) / Numb * DimsProd(nDim - 1, dims + 1);
size_t *p, *q;
p= data; q= data0;
ptrdiff_t i;
for (i= 0; i < ndat; i++, p++, q++) *p= *q;
}
void SetData (UnOpBin op = IL_BU_SET, int val = 0) {
if (!data) throw Allocerr();
ptrdiff_t ndat= (dims[0] + Numb - 1) / Numb * DimsProd(nDim - 1, dims + 1);
size_t *p;
p= data;
ptrdiff_t i;
switch (op) {
case IL_BU_SET:
if (val) {
for (i= 0; i < ndat; i++) *(p++)= AB;
EraseEnd();
}
else for (i= 0; i < ndat; i++) *(p++)= 0;
break;
case IL_BU_NOT: {
for (i= 0; i < ndat; i++) *(p++)= ~*p;
EraseEnd();
}
}
}
void SetData (BinOpBin op, size_t* data0) {
if (!data) throw Allocerr();
ptrdiff_t ndat= (dims[0] + Numb - 1) / Numb * DimsProd(nDim - 1, dims + 1);
size_t *p, *q;
p= data;
q= data0;
ptrdiff_t i;
switch (op) {
case IL_BB_OR:
for (i= 0; i < ndat; i++, p++, q++) *p|= *q;
break;
case IL_BB_AND:
for (i= 0; i < ndat; i++, p++, q++) *p&= *q;
break;
case IL_BB_XOR: {
for (i= 0; i < ndat; i++, p++, q++) *p^= *q;
EraseEnd();
}
break;
case IL_BB_DIF:
for (i= 0; i < ndat; i++, p++, q++) *p&= ~*q;
break;
case IL_BB_NIF:
for (i= 0; i < ndat; i++, p++, q++) *p|= ~*q;
break;
case IL_BB_XORN: {
for (i= 0; i < ndat; i++, p++, q++) *p^= ~*q;
EraseEnd();
}
}
}
};
class Bin : public Binary {
public:
Bin () : Binary () { }
Bin (const Bin& x0) : Binary (x0) { NewData (); } //overrides default!
Bin (Binary& x0) : Binary (x0) { NewData (); }
Bin (const Frame& x0) : Binary (x0) { NewData (); }
Bin (int nDim0, ptrdiff_t* dims0, double* bounds0 = 0) : Binary (nDim0, dims0, bounds0) {
NewData ();
}
~Bin () {
if (data) delete[] data;
}
};
class Object {
public:
Object() : nVar (0) , data (0) { }
Object(const Object& x0) : nVar (x0.nVar) {
data= new double[nVar];
int i;
for (i= 0; i < nVar; i++) data[i]= x0.data[i];
}
explicit Object(int nVar0) : nVar (nVar0) {
data= new double[nVar];
}
templateObject(int nVar0, G* data0) : nVar (nVar0) {
data= new double[nVar];
double *p;
G *q;
p= data;
q= data0;
int i;
for (i= 0; i < nVar; i++) *(p++)= *(q++);
}
Object& operator=(const Object& x0) {
if (!data || nVar != x0.nVar) {
nVar= x0.nVar;
if (data) delete data;
data= new double[nVar];
}
int i;
for (i= 0; i < nVar; i++) data[i]= x0.data[i];
return *this;
}
templatevoid SetData (G* data0, int num = nVar) {
if (!data) data= new double[nVar];
if (num > nVar) num= nVar;
double *p;
G *q;
p= data;
q= data0;
int i;
for (i= 0; i < num; i++) *(p++)= *(q++);
}
double Dist(const Object& arg, int nDim) {
if (!data) return 0.0;
double sum= 0.0;
int i;
for (i= 0; i < nDim; i++) {
double pom= data[i] - arg.data[i];
sum+= pom * pom;
}
return sqrt(sum);
}
double operator[](int i) const { return data[i]; }
double& operator[](int i) { return data[i]; }
~Object() {
if (data) delete [] data;
}
int nVar;
double* data;
};
typedef std::vector