summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/CompositeGeometryManager.cpp321
-rw-r--r--src/FilteredBackProjectionAlgorithm.cpp2
-rw-r--r--src/Fourier.cpp206
-rw-r--r--src/Utilities.cpp34
4 files changed, 281 insertions, 282 deletions
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index c5b4d3b..7c4f8e6 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -197,6 +197,69 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div
return true;
}
+
+static std::pair<double, double> reduceProjectionVertical(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom)
+{
+ double vmin_g, vmax_g;
+
+ // reduce self to only cover intersection with projection of VolumePart
+ // (Project corners of volume, take bounding box)
+
+ for (int i = 0; i < pProjGeom->getProjectionCount(); ++i) {
+
+ double vol_u[8];
+ double vol_v[8];
+
+ double pixx = pVolGeom->getPixelLengthX();
+ double pixy = pVolGeom->getPixelLengthY();
+ double pixz = pVolGeom->getPixelLengthZ();
+
+ // TODO: Is 0.5 sufficient?
+ double xmin = pVolGeom->getWindowMinX() - 0.5 * pixx;
+ double xmax = pVolGeom->getWindowMaxX() + 0.5 * pixx;
+ double ymin = pVolGeom->getWindowMinY() - 0.5 * pixy;
+ double ymax = pVolGeom->getWindowMaxY() + 0.5 * pixy;
+ double zmin = pVolGeom->getWindowMinZ() - 0.5 * pixz;
+ double zmax = pVolGeom->getWindowMaxZ() + 0.5 * pixz;
+
+ pProjGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]);
+ pProjGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]);
+ pProjGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]);
+ pProjGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]);
+ pProjGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]);
+ pProjGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]);
+ pProjGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]);
+ pProjGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]);
+
+ double vmin = vol_v[0];
+ double vmax = vol_v[0];
+
+ for (int j = 1; j < 8; ++j) {
+ if (vol_v[j] < vmin)
+ vmin = vol_v[j];
+ if (vol_v[j] > vmax)
+ vmax = vol_v[j];
+ }
+
+ if (i == 0 || vmin < vmin_g)
+ vmin_g = vmin;
+ if (i == 0 || vmax > vmax_g)
+ vmax_g = vmax;
+ }
+
+ if (vmin_g < -1.0)
+ vmin_g = -1.0;
+ if (vmax_g > pProjGeom->getDetectorRowCount())
+ vmax_g = pProjGeom->getDetectorRowCount();
+
+ if (vmin_g >= vmax_g)
+ vmin_g = vmax_g = 0.0;
+
+ return std::pair<double, double>(vmin_g, vmax_g);
+
+}
+
+
CCompositeGeometryManager::CPart::CPart(const CPart& other)
{
eType = other.eType;
@@ -237,119 +300,135 @@ size_t CCompositeGeometryManager::CPart::getSize()
}
+static bool testVolumeRange(const std::pair<double, double>& fullRange,
+ const CVolumeGeometry3D *pVolGeom,
+ const CProjectionGeometry3D *pProjGeom,
+ int zmin, int zmax)
+{
+ double pixz = pVolGeom->getPixelLengthZ();
+
+ CVolumeGeometry3D test(pVolGeom->getGridColCount(),
+ pVolGeom->getGridRowCount(),
+ zmax - zmin,
+ pVolGeom->getWindowMinX(),
+ pVolGeom->getWindowMinY(),
+ pVolGeom->getWindowMinZ() + zmin * pixz,
+ pVolGeom->getWindowMaxX(),
+ pVolGeom->getWindowMaxY(),
+ pVolGeom->getWindowMinZ() + zmax * pixz);
+
+
+ std::pair<double, double> subRange = reduceProjectionVertical(&test, pProjGeom);
+
+ // empty
+ if (subRange.first == subRange.second)
+ return true;
+
+ // fully outside of fullRange
+ if (subRange.first >= fullRange.second || subRange.second <= fullRange.first)
+ return true;
+
+ return false;
+}
+
CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other)
{
const CProjectionPart *other = dynamic_cast<const CProjectionPart *>(_other);
assert(other);
- // TODO: Is 0.5 sufficient?
- double umin = -0.5;
- double umax = other->pGeom->getDetectorColCount() + 0.5;
- double vmin = -0.5;
- double vmax = other->pGeom->getDetectorRowCount() + 0.5;
-
- double uu[4];
- double vv[4];
- uu[0] = umin; vv[0] = vmin;
- uu[1] = umin; vv[1] = vmax;
- uu[2] = umax; vv[2] = vmin;
- uu[3] = umax; vv[3] = vmax;
-
- double pixx = pGeom->getPixelLengthX();
- double pixy = pGeom->getPixelLengthY();
- double pixz = pGeom->getPixelLengthZ();
- double xmin = pGeom->getWindowMinX() - 0.5 * pixx;
- double xmax = pGeom->getWindowMaxX() + 0.5 * pixx;
- double ymin = pGeom->getWindowMinY() - 0.5 * pixy;
- double ymax = pGeom->getWindowMaxY() + 0.5 * pixy;
-
- // NB: Flipped
- double zmax = pGeom->getWindowMinZ() - 2.5 * pixz;
- double zmin = pGeom->getWindowMaxZ() + 2.5 * pixz;
-
- // TODO: This isn't as tight as it could be.
- // In particular it won't detect the detector being
- // missed entirely on the u side.
-
- for (int i = 0; i < other->pGeom->getProjectionCount(); ++i) {
- for (int j = 0; j < 4; ++j) {
- double px, py, pz;
-
- other->pGeom->backprojectPointX(i, uu[j], vv[j], xmin, py, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
- other->pGeom->backprojectPointX(i, uu[j], vv[j], xmax, py, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
-
- other->pGeom->backprojectPointY(i, uu[j], vv[j], ymin, px, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
- other->pGeom->backprojectPointY(i, uu[j], vv[j], ymax, px, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
+ std::pair<double, double> fullRange = reduceProjectionVertical(pGeom, other->pGeom);
+
+ int top_slice = 0, bottom_slice = 0;
+
+ if (fullRange.first < fullRange.second) {
+
+
+ // TOP SLICE
+
+ int zmin = 0;
+ int zmax = pGeom->getGridSliceCount()-1; // (Don't try empty region)
+
+ // Setting top slice to zmin is always valid.
+
+ while (zmin < zmax) {
+ int zmid = (zmin + zmax + 1) / 2;
+
+ bool ok = testVolumeRange(fullRange, pGeom, other->pGeom,
+ 0, zmid);
+
+ ASTRA_DEBUG("binsearch min: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much");
+
+ if (ok)
+ zmin = zmid;
+ else
+ zmax = zmid - 1;
}
- }
- //ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax);
+ top_slice = zmin;
- // Clip both zmin and zmax to get rid of extreme (or infinite) values
- // NB: When individual pz values are +/-Inf, the sign is determined
- // by ray direction and on which side of the face the ray passes.
- if (zmin < pGeom->getWindowMinZ() - 2*pixz)
- zmin = pGeom->getWindowMinZ() - 2*pixz;
- if (zmin > pGeom->getWindowMaxZ() + 2*pixz)
- zmin = pGeom->getWindowMaxZ() + 2*pixz;
- if (zmax < pGeom->getWindowMinZ() - 2*pixz)
- zmax = pGeom->getWindowMinZ() - 2*pixz;
- if (zmax > pGeom->getWindowMaxZ() + 2*pixz)
- zmax = pGeom->getWindowMaxZ() + 2*pixz;
- zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz;
- zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz;
+ // BOTTOM SLICE
- int _zmin = (int)floor(zmin);
- int _zmax = (int)ceil(zmax);
+ zmin = top_slice + 1; // (Don't try empty region)
+ zmax = pGeom->getGridSliceCount();
- //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax);
+ // Setting bottom slice to zmax is always valid
+
+ while (zmin < zmax) {
+ int zmid = (zmin + zmax) / 2;
+
+ bool ok = testVolumeRange(fullRange, pGeom, other->pGeom,
+ zmid, pGeom->getGridSliceCount());
+
+ ASTRA_DEBUG("binsearch max: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much");
+
+ if (ok)
+ zmax = zmid;
+ else
+ zmin = zmid + 1;
+
+ }
- if (_zmin < 0)
- _zmin = 0;
- if (_zmax > pGeom->getGridSliceCount())
- _zmax = pGeom->getGridSliceCount();
+ bottom_slice = zmax;
- if (_zmax <= _zmin) {
- _zmin = _zmax = 0;
}
- //ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax);
+
+ ASTRA_DEBUG("found extent: %d - %d", top_slice, bottom_slice);
+
+ top_slice -= 1;
+ if (top_slice < 0)
+ top_slice = 0;
+ bottom_slice += 1;
+ if (bottom_slice >= pGeom->getGridSliceCount())
+ bottom_slice = pGeom->getGridSliceCount();
+
+ ASTRA_DEBUG("adjusted extent: %d - %d", top_slice, bottom_slice);
+
+ double pixz = pGeom->getPixelLengthZ();
CVolumePart *sub = new CVolumePart();
sub->subX = this->subX;
sub->subY = this->subY;
- sub->subZ = this->subZ + _zmin;
+ sub->subZ = this->subZ + top_slice;
sub->pData = pData;
- if (_zmin == _zmax) {
+ if (top_slice == bottom_slice) {
sub->pGeom = 0;
} else {
sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(),
pGeom->getGridRowCount(),
- _zmax - _zmin,
+ bottom_slice - top_slice,
pGeom->getWindowMinX(),
pGeom->getWindowMinY(),
- pGeom->getWindowMinZ() + _zmin * pixz,
+ pGeom->getWindowMinZ() + top_slice * pixz,
pGeom->getWindowMaxX(),
pGeom->getWindowMaxY(),
- pGeom->getWindowMinZ() + _zmax * pixz);
+ pGeom->getWindowMinZ() + bottom_slice * pixz);
}
- ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax);
+ ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d ( %f - %f )", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + top_slice, this->subZ + bottom_slice, pGeom->getWindowMinZ() + top_slice * pixz, pGeom->getWindowMinZ() + bottom_slice * pixz);
return sub;
}
@@ -584,7 +663,9 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -631,7 +712,9 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -678,7 +761,9 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -743,62 +828,16 @@ void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, s
}
+
CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other)
{
const CVolumePart *other = dynamic_cast<const CVolumePart *>(_other);
assert(other);
- double vmin_g, vmax_g;
-
- // reduce self to only cover intersection with projection of VolumePart
- // (Project corners of volume, take bounding box)
-
- for (int i = 0; i < pGeom->getProjectionCount(); ++i) {
-
- double vol_u[8];
- double vol_v[8];
-
- double pixx = other->pGeom->getPixelLengthX();
- double pixy = other->pGeom->getPixelLengthY();
- double pixz = other->pGeom->getPixelLengthZ();
-
- // TODO: Is 0.5 sufficient?
- double xmin = other->pGeom->getWindowMinX() - 0.5 * pixx;
- double xmax = other->pGeom->getWindowMaxX() + 0.5 * pixx;
- double ymin = other->pGeom->getWindowMinY() - 0.5 * pixy;
- double ymax = other->pGeom->getWindowMaxY() + 0.5 * pixy;
- double zmin = other->pGeom->getWindowMinZ() - 0.5 * pixz;
- double zmax = other->pGeom->getWindowMaxZ() + 0.5 * pixz;
-
- pGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]);
- pGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]);
- pGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]);
- pGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]);
- pGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]);
- pGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]);
- pGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]);
- pGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]);
-
- double vmin = vol_v[0];
- double vmax = vol_v[0];
-
- for (int j = 1; j < 8; ++j) {
- if (vol_v[j] < vmin)
- vmin = vol_v[j];
- if (vol_v[j] > vmax)
- vmax = vol_v[j];
- }
-
- if (i == 0 || vmin < vmin_g)
- vmin_g = vmin;
- if (i == 0 || vmax > vmax_g)
- vmax_g = vmax;
- }
-
- // fprintf(stderr, "v extent: %f %f\n", vmin_g, vmax_g);
-
- int _vmin = (int)floor(vmin_g - 1.0f);
- int _vmax = (int)ceil(vmax_g + 1.0f);
+ std::pair<double, double> r = reduceProjectionVertical(other->pGeom, pGeom);
+ // fprintf(stderr, "v extent: %f %f\n", r.first, r.second);
+ int _vmin = (int)floor(r.first - 1.0);
+ int _vmax = (int)ceil(r.second + 1.0);
if (_vmin < 0)
_vmin = 0;
if (_vmax > pGeom->getDetectorRowCount())
@@ -838,7 +877,11 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
+
+ ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
for (int x = -(rem / 2); x < sliceCount; x += blockSize) {
int newsubX = x;
@@ -880,7 +923,11 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
+
+ ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
for (int z = -(rem / 2); z < sliceCount; z += blockSize) {
int newsubZ = z;
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index 90efd52..70462f7 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -275,7 +275,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
float32* pf = new float32[2 * iAngleCount * zpDetector];
- int *ip = new int[int(2+sqrt(zpDetector)+1)];
+ int *ip = new int[int(2+sqrt((float)zpDetector)+1)];
ip[0]=0;
float32 *w = new float32[zpDetector/2];
diff --git a/src/Fourier.cpp b/src/Fourier.cpp
index 5ca22e6..c33f7bd 100644
--- a/src/Fourier.cpp
+++ b/src/Fourier.cpp
@@ -27,6 +27,7 @@ $Id$
*/
#include "astra/Fourier.h"
+#include <cmath>
namespace astra {
@@ -320,11 +321,45 @@ Appendix :
*/
-void cdft(int n, int isgn, float32 *a, int *ip, float32 *w)
+static int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
+static void bitrv208(float32 *a);
+static void bitrv208neg(float32 *a);
+static void bitrv216(float32 *a);
+static void bitrv216neg(float32 *a);
+static void bitrv2conj(int n, int *ip, float32 *a);
+static void bitrv2(int n, int *ip, float32 *a);
+static void cftb040(float32 *a);
+static void cftb1st(int n, float32 *a, float32 *w);
+static void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
+static void cftf040(float32 *a);
+static void cftf081(float32 *a, float32 *w);
+static void cftf082(float32 *a, float32 *w);
+static void cftf161(float32 *a, float32 *w);
+static void cftf162(float32 *a, float32 *w);
+static void cftf1st(int n, float32 *a, float32 *w);
+static void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
+static void cftfx41(int n, float32 *a, int nw, float32 *w);
+static void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
+static void cftmdl1(int n, float32 *a, float32 *w);
+static void cftmdl2(int n, float32 *a, float32 *w);
+static void *cftrec1_th(void *p);
+static void *cftrec2_th(void *p);
+static void cftrec4(int n, float32 *a, int nw, float32 *w);
+static void cftx020(float32 *a);
+static void dctsub(int n, float32 *a, int nc, float32 *c);
+static void dstsub(int n, float32 *a, int nc, float32 *c);
+static void makect(int nc, int *ip, float32 *c);
+static void makeipt(int nw, int *ip);
+static void makewt(int nw, int *ip, float32 *w);
+static void rftbsub(int n, float32 *a, int nc, float32 *c);
+static void rftfsub(int n, float32 *a, int nc, float32 *c);
+#ifdef USE_CDFT_THREADS
+static void cftrec4_th(int n, float32 *a, int nw, float32 *w);
+#endif /* USE_CDFT_THREADS */
+
+
+_AstraExport void cdft(int n, int isgn, float32 *a, int *ip, float32 *w)
{
- void makewt(int nw, int *ip, float32 *w);
- void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
int nw;
nw = ip[0];
@@ -340,14 +375,8 @@ void cdft(int n, int isgn, float32 *a, int *ip, float32 *w)
}
-void rdft(int n, int isgn, float32 *a, int *ip, float32 *w)
+_AstraExport void rdft(int n, int isgn, float32 *a, int *ip, float32 *w)
{
- void makewt(int nw, int *ip, float32 *w);
- void makect(int nc, int *ip, float32 *c);
- void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void rftfsub(int n, float32 *a, int nc, float32 *c);
- void rftbsub(int n, float32 *a, int nc, float32 *c);
int nw, nc;
float32 xi;
@@ -384,15 +413,8 @@ void rdft(int n, int isgn, float32 *a, int *ip, float32 *w)
}
-void ddct(int n, int isgn, float32 *a, int *ip, float32 *w)
+_AstraExport void ddct(int n, int isgn, float32 *a, int *ip, float32 *w)
{
- void makewt(int nw, int *ip, float32 *w);
- void makect(int nc, int *ip, float32 *c);
- void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void rftfsub(int n, float32 *a, int nc, float32 *c);
- void rftbsub(int n, float32 *a, int nc, float32 *c);
- void dctsub(int n, float32 *a, int nc, float32 *c);
int j, nw, nc;
float32 xr;
@@ -440,15 +462,8 @@ void ddct(int n, int isgn, float32 *a, int *ip, float32 *w)
}
-void ddst(int n, int isgn, float32 *a, int *ip, float32 *w)
+_AstraExport void ddst(int n, int isgn, float32 *a, int *ip, float32 *w)
{
- void makewt(int nw, int *ip, float32 *w);
- void makect(int nc, int *ip, float32 *c);
- void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void rftfsub(int n, float32 *a, int nc, float32 *c);
- void rftbsub(int n, float32 *a, int nc, float32 *c);
- void dstsub(int n, float32 *a, int nc, float32 *c);
int j, nw, nc;
float32 xr;
@@ -496,13 +511,8 @@ void ddst(int n, int isgn, float32 *a, int *ip, float32 *w)
}
-void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w)
+_AstraExport void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w)
{
- void makewt(int nw, int *ip, float32 *w);
- void makect(int nc, int *ip, float32 *c);
- void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void rftfsub(int n, float32 *a, int nc, float32 *c);
- void dctsub(int n, float32 *a, int nc, float32 *c);
int j, k, l, m, mh, nw, nc;
float32 xr, xi, yr, yi;
@@ -589,13 +599,8 @@ void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w)
}
-void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w)
+_AstraExport void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w)
{
- void makewt(int nw, int *ip, float32 *w);
- void makect(int nc, int *ip, float32 *c);
- void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
- void rftfsub(int n, float32 *a, int nc, float32 *c);
- void dstsub(int n, float32 *a, int nc, float32 *c);
int j, k, l, m, mh, nw, nc;
float32 xr, xi, yr, yi;
@@ -675,12 +680,8 @@ void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w)
/* -------- initializing routines -------- */
-
-#include <math.h>
-
-void makewt(int nw, int *ip, float32 *w)
+static void makewt(int nw, int *ip, float32 *w)
{
- void makeipt(int nw, int *ip);
int j, nwh, nw0, nw1;
float32 delta, wn4r, wk1r, wk1i, wk3r, wk3i;
@@ -739,7 +740,7 @@ void makewt(int nw, int *ip, float32 *w)
}
-void makeipt(int nw, int *ip)
+static void makeipt(int nw, int *ip)
{
int j, l, m, m2, p, q;
@@ -759,7 +760,7 @@ void makeipt(int nw, int *ip)
}
-void makect(int nc, int *ip, float32 *c)
+static void makect(int nc, int *ip, float32 *c)
{
int j, nch;
float32 delta;
@@ -835,23 +836,8 @@ void makect(int nc, int *ip, float32 *c)
#endif /* USE_CDFT_WINTHREADS */
-void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w)
+static void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w)
{
- void bitrv2(int n, int *ip, float32 *a);
- void bitrv216(float32 *a);
- void bitrv208(float32 *a);
- void cftf1st(int n, float32 *a, float32 *w);
- void cftrec4(int n, float32 *a, int nw, float32 *w);
- void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
- void cftfx41(int n, float32 *a, int nw, float32 *w);
- void cftf161(float32 *a, float32 *w);
- void cftf081(float32 *a, float32 *w);
- void cftf040(float32 *a);
- void cftx020(float32 *a);
-#ifdef USE_CDFT_THREADS
- void cftrec4_th(int n, float32 *a, int nw, float32 *w);
-#endif /* USE_CDFT_THREADS */
-
if (n > 8) {
if (n > 32) {
cftf1st(n, a, &w[nw - (n >> 2)]);
@@ -883,23 +869,8 @@ void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w)
}
-void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w)
+static void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w)
{
- void bitrv2conj(int n, int *ip, float32 *a);
- void bitrv216neg(float32 *a);
- void bitrv208neg(float32 *a);
- void cftb1st(int n, float32 *a, float32 *w);
- void cftrec4(int n, float32 *a, int nw, float32 *w);
- void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
- void cftfx41(int n, float32 *a, int nw, float32 *w);
- void cftf161(float32 *a, float32 *w);
- void cftf081(float32 *a, float32 *w);
- void cftb040(float32 *a);
- void cftx020(float32 *a);
-#ifdef USE_CDFT_THREADS
- void cftrec4_th(int n, float32 *a, int nw, float32 *w);
-#endif /* USE_CDFT_THREADS */
-
if (n > 8) {
if (n > 32) {
cftb1st(n, a, &w[nw - (n >> 2)]);
@@ -931,7 +902,7 @@ void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w)
}
-void bitrv2(int n, int *ip, float32 *a)
+static void bitrv2(int n, int *ip, float32 *a)
{
int j, j1, k, k1, l, m, nh, nm;
float32 xr, xi, yr, yi;
@@ -1278,7 +1249,7 @@ void bitrv2(int n, int *ip, float32 *a)
}
-void bitrv2conj(int n, int *ip, float32 *a)
+static void bitrv2conj(int n, int *ip, float32 *a)
{
int j, j1, k, k1, l, m, nh, nm;
float32 xr, xi, yr, yi;
@@ -1633,7 +1604,7 @@ void bitrv2conj(int n, int *ip, float32 *a)
}
-void bitrv216(float32 *a)
+static void bitrv216(float32 *a)
{
float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i,
@@ -1690,7 +1661,7 @@ void bitrv216(float32 *a)
}
-void bitrv216neg(float32 *a)
+static void bitrv216neg(float32 *a)
{
float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i,
@@ -1760,7 +1731,7 @@ void bitrv216neg(float32 *a)
}
-void bitrv208(float32 *a)
+static void bitrv208(float32 *a)
{
float32 x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
@@ -1783,7 +1754,7 @@ void bitrv208(float32 *a)
}
-void bitrv208neg(float32 *a)
+static void bitrv208neg(float32 *a)
{
float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
x5r, x5i, x6r, x6i, x7r, x7i;
@@ -1819,7 +1790,7 @@ void bitrv208neg(float32 *a)
}
-void cftf1st(int n, float32 *a, float32 *w)
+static void cftf1st(int n, float32 *a, float32 *w)
{
int j, j0, j1, j2, j3, k, m, mh;
float32 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
@@ -2025,7 +1996,7 @@ void cftf1st(int n, float32 *a, float32 *w)
}
-void cftb1st(int n, float32 *a, float32 *w)
+static void cftb1st(int n, float32 *a, float32 *w)
{
int j, j0, j1, j2, j3, k, m, mh;
float32 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
@@ -2242,10 +2213,8 @@ struct cdft_arg_st {
typedef struct cdft_arg_st cdft_arg_t;
-void cftrec4_th(int n, float32 *a, int nw, float32 *w)
+static void cftrec4_th(int n, float32 *a, int nw, float32 *w)
{
- void *cftrec1_th(void *p);
- void *cftrec2_th(void *p);
int i, idiv4, m, nthread;
cdft_thread_t th[4];
cdft_arg_t ag[4];
@@ -2276,11 +2245,8 @@ void cftrec4_th(int n, float32 *a, int nw, float32 *w)
}
-void *cftrec1_th(void *p)
+static void *cftrec1_th(void *p)
{
- int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
- void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
- void cftmdl1(int n, float32 *a, float32 *w);
int isplt, j, k, m, n, n0, nw;
float32 *a, *w;
@@ -2305,11 +2271,8 @@ void *cftrec1_th(void *p)
}
-void *cftrec2_th(void *p)
+static void *cftrec2_th(void *p)
{
- int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
- void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
- void cftmdl2(int n, float32 *a, float32 *w);
int isplt, j, k, m, n, n0, nw;
float32 *a, *w;
@@ -2337,11 +2300,8 @@ void *cftrec2_th(void *p)
#endif /* USE_CDFT_THREADS */
-void cftrec4(int n, float32 *a, int nw, float32 *w)
+static void cftrec4(int n, float32 *a, int nw, float32 *w)
{
- int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
- void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
- void cftmdl1(int n, float32 *a, float32 *w);
int isplt, j, k, m;
m = n;
@@ -2361,8 +2321,6 @@ void cftrec4(int n, float32 *a, int nw, float32 *w)
int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w)
{
- void cftmdl1(int n, float32 *a, float32 *w);
- void cftmdl2(int n, float32 *a, float32 *w);
int i, isplt, m;
if ((k & 3) != 0) {
@@ -2394,15 +2352,8 @@ int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w)
}
-void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w)
+static void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w)
{
- void cftmdl1(int n, float32 *a, float32 *w);
- void cftmdl2(int n, float32 *a, float32 *w);
- void cftf161(float32 *a, float32 *w);
- void cftf162(float32 *a, float32 *w);
- void cftf081(float32 *a, float32 *w);
- void cftf082(float32 *a, float32 *w);
-
if (n == 512) {
cftmdl1(128, a, &w[nw - 64]);
cftf161(a, &w[nw - 8]);
@@ -2459,7 +2410,7 @@ void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w)
}
-void cftmdl1(int n, float32 *a, float32 *w)
+static void cftmdl1(int n, float32 *a, float32 *w)
{
int j, j0, j1, j2, j3, k, m, mh;
float32 wn4r, wk1r, wk1i, wk3r, wk3i;
@@ -2569,7 +2520,7 @@ void cftmdl1(int n, float32 *a, float32 *w)
}
-void cftmdl2(int n, float32 *a, float32 *w)
+static void cftmdl2(int n, float32 *a, float32 *w)
{
int j, j0, j1, j2, j3, k, kr, m, mh;
float32 wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
@@ -2703,13 +2654,8 @@ void cftmdl2(int n, float32 *a, float32 *w)
}
-void cftfx41(int n, float32 *a, int nw, float32 *w)
+static void cftfx41(int n, float32 *a, int nw, float32 *w)
{
- void cftf161(float32 *a, float32 *w);
- void cftf162(float32 *a, float32 *w);
- void cftf081(float32 *a, float32 *w);
- void cftf082(float32 *a, float32 *w);
-
if (n == 128) {
cftf161(a, &w[nw - 8]);
cftf162(&a[32], &w[nw - 32]);
@@ -2724,7 +2670,7 @@ void cftfx41(int n, float32 *a, int nw, float32 *w)
}
-void cftf161(float32 *a, float32 *w)
+static void cftf161(float32 *a, float32 *w)
{
float32 wn4r, wk1r, wk1i,
x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
@@ -2883,7 +2829,7 @@ void cftf161(float32 *a, float32 *w)
}
-void cftf162(float32 *a, float32 *w)
+static void cftf162(float32 *a, float32 *w)
{
float32 wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
x0r, x0i, x1r, x1i, x2r, x2i,
@@ -3066,7 +3012,7 @@ void cftf162(float32 *a, float32 *w)
}
-void cftf081(float32 *a, float32 *w)
+static void cftf081(float32 *a, float32 *w)
{
float32 wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
@@ -3128,7 +3074,7 @@ void cftf081(float32 *a, float32 *w)
}
-void cftf082(float32 *a, float32 *w)
+static void cftf082(float32 *a, float32 *w)
{
float32 wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i,
y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
@@ -3200,7 +3146,7 @@ void cftf082(float32 *a, float32 *w)
}
-void cftf040(float32 *a)
+static void cftf040(float32 *a)
{
float32 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
@@ -3223,7 +3169,7 @@ void cftf040(float32 *a)
}
-void cftb040(float32 *a)
+static void cftb040(float32 *a)
{
float32 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
@@ -3246,7 +3192,7 @@ void cftb040(float32 *a)
}
-void cftx020(float32 *a)
+static void cftx020(float32 *a)
{
float32 x0r, x0i;
@@ -3259,7 +3205,7 @@ void cftx020(float32 *a)
}
-void rftfsub(int n, float32 *a, int nc, float32 *c)
+static void rftfsub(int n, float32 *a, int nc, float32 *c)
{
int j, k, kk, ks, m;
float32 wkr, wki, xr, xi, yr, yi;
@@ -3284,7 +3230,7 @@ void rftfsub(int n, float32 *a, int nc, float32 *c)
}
-void rftbsub(int n, float32 *a, int nc, float32 *c)
+static void rftbsub(int n, float32 *a, int nc, float32 *c)
{
int j, k, kk, ks, m;
float32 wkr, wki, xr, xi, yr, yi;
@@ -3309,7 +3255,7 @@ void rftbsub(int n, float32 *a, int nc, float32 *c)
}
-void dctsub(int n, float32 *a, int nc, float32 *c)
+static void dctsub(int n, float32 *a, int nc, float32 *c)
{
int j, k, kk, ks, m;
float32 wkr, wki, xr;
@@ -3330,7 +3276,7 @@ void dctsub(int n, float32 *a, int nc, float32 *c)
}
-void dstsub(int n, float32 *a, int nc, float32 *c)
+static void dstsub(int n, float32 *a, int nc, float32 *c)
{
int j, k, kk, ks, m;
float32 wkr, wki, xr;
diff --git a/src/Utilities.cpp b/src/Utilities.cpp
index c9740bf..8b0ca94 100644
--- a/src/Utilities.cpp
+++ b/src/Utilities.cpp
@@ -28,10 +28,6 @@ $Id$
#include "astra/Utilities.h"
-#include <boost/algorithm/string.hpp>
-#include <boost/algorithm/string/split.hpp>
-#include <boost/algorithm/string/classification.hpp>
-
#include <sstream>
#include <locale>
#include <iomanip>
@@ -84,18 +80,16 @@ std::vector<double> stringToDoubleVector(const std::string &s)
template<typename T>
std::vector<T> stringToVector(const std::string& s)
{
- // split
- std::vector<std::string> items;
- boost::split(items, s, boost::is_any_of(",;"));
-
- // init list
std::vector<T> out;
- out.resize(items.size());
+ size_t current = 0;
+ size_t next;
+ do {
+ next = s.find_first_of(",;", current);
+ std::string t = s.substr(current, next - current);
+ out.push_back(stringTo<T>(t));
+ current = next + 1;
+ } while (next != std::string::npos);
- // loop elements
- for (unsigned int i = 0; i < items.size(); i++) {
- out[i] = stringTo<T>(items[i]);
- }
return out;
}
@@ -120,6 +114,18 @@ std::string doubleToString(double f)
template<> std::string toString(float f) { return floatToString(f); }
template<> std::string toString(double f) { return doubleToString(f); }
+void splitString(std::vector<std::string> &items, const std::string& s,
+ const char *delim)
+{
+ items.clear();
+ size_t current = 0;
+ size_t next;
+ do {
+ next = s.find_first_of(",;", current);
+ items.push_back(s.substr(current, next - current));
+ current = next + 1;
+ } while (next != std::string::npos);
+}
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 13:53:54 +0000