* sgp.cpp: Fixed bug in drawCircle
+ * filter.cpp: Fixed Hanning parameter to be 0.5 rather than 0.54
+
3.0.3 - Released 2/20/01
* ctsim: Fixed core dump on Linux with OpenGL
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: fnetorderstream.h,v 1.8 2001/01/28 19:10:18 kevin Exp $
+** $Id: fnetorderstream.h,v 1.9 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
p[4] = temp;
}
+inline void
+SwapBytes2IfLittleEndian (void* buffer)
+{
+#ifndef WORDS_BIGENDIAN
+ SwapBytes2 (buffer);
+#endif
+}
+
+inline void
+SwapBytes4IfLittleEndian (void* buffer)
+{
+#ifndef WORDS_BIGENDIAN
+ SwapBytes4 (buffer);
+#endif
+}
+
+inline void
+SwapBytes8IfLittleEndian (void* buffer)
+{
+#ifndef WORDS_BIGENDIAN
+ SwapBytes8 (buffer);
+#endif
+}
void ConvertNetworkOrder (void* buffer, size_t bytes);
void ConvertReverseNetworkOrder (void* buffer, size_t bytes);
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.h,v 1.31 2001/03/11 12:37:34 kevin Exp $
+** $Id: projections.h,v 1.32 2001/03/13 04:44:25 kevin Exp $
**
**
** This program is free software; you can redistribute it and/or modify
double m_dTheta; // perpendicular angle to origin
double m_dRaysum;
- bool lessThanT (ParallelRaysumCoordinate& rCompare)
- {return m_dT < rCompare.m_dT; }
- bool lessThanTheta (ParallelRaysumCoordinate& rCompare)
- {return m_dTheta < rCompare.m_dTheta; }
-
static bool compareByTheta (ParallelRaysumCoordinate* a, ParallelRaysumCoordinate* b)
{ return a->m_dTheta == b->m_dTheta ? b->m_dT > a->m_dT : b->m_dTheta > a->m_dTheta; }
class ParallelRaysums {
public:
- ParallelRaysums (Projections* pProjections);
+
+ enum {
+ THETA_RANGE_UNCONSTRAINED = 0,
+ THETA_RANGE_NORMALIZE_TO_TWOPI,
+ THETA_RANGE_FOLD_TO_PI,
+ };
+
+ ParallelRaysums (Projections* pProjections, int iThetaRange);
~ParallelRaysums ();
typedef std::vector<ParallelRaysumCoordinate*> CoordinateContainer;
void getLimits (double* dMinT, double* dMaxT, double* dMinTheta, double* dMaxTheta) const;
CoordinateContainer& getSortedByT();
CoordinateContainer& getSortedByTheta();
- double interpolate (double* pdX, double* pdY, int n, double dXValue);
+ double interpolate (double* pdX, double* pdY, int n, double dXValue, int* iLastFloor = NULL);
void getThetaAndRaysumsForT (int iT, double* pdTheta, double* pdRaysum);
void getDetPositions (double* pdDetPos);
int m_iNumCoordinates;
int m_iNumView;
int m_iNumDet;
+ int m_iThetaRange;
};
** This is part of the CTSim program
** Copyright (c) 1983-2000 Kevin Rosenberg
**
-** $Id: filter.cpp,v 1.38 2001/02/22 18:22:40 kevin Exp $
+** $Id: filter.cpp,v 1.39 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
const int SignalFilter::DOMAIN_INVALID = -1;
const int SignalFilter::DOMAIN_FREQUENCY = 0;
const int SignalFilter::DOMAIN_SPATIAL = 1;
-
+
const char* const SignalFilter::s_aszDomainName[] = {
{"frequency"},
{"spatial"},
/* NAME
- * SignalFilter::SignalFilter Construct a signal
- *
- * SYNOPSIS
- * f = SignalFilter (filt_type, bw, filterMin, filterMax, n, param, domain, analytic)
- * double f Generated filter vector
- * int filt_type Type of filter wanted
- * double bw Bandwidth of filter
- * double filterMin, filterMax Filter limits
- * int nFilterPoints Number of points in signal
- * double param General input parameter to filters
- * int domain FREQUENCY or SPATIAL domain wanted
- */
+* SignalFilter::SignalFilter Construct a signal
+*
+* SYNOPSIS
+* f = SignalFilter (filt_type, bw, filterMin, filterMax, n, param, domain, analytic)
+* double f Generated filter vector
+* int filt_type Type of filter wanted
+* double bw Bandwidth of filter
+* double filterMin, filterMax Filter limits
+* int nFilterPoints Number of points in signal
+* double param General input parameter to filters
+* int domain FREQUENCY or SPATIAL domain wanted
+*/
SignalFilter::SignalFilter (const char* szFilterName, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const char* szDomainName)
- : m_adFilter(NULL), m_fail(false)
+: m_adFilter(NULL), m_fail(false)
{
m_idFilter = convertFilterNameToID (szFilterName);
if (m_idFilter == FILTER_INVALID) {
}
SignalFilter::SignalFilter (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain)
- : m_adFilter(NULL), m_fail(false)
+: m_adFilter(NULL), m_fail(false)
{
init (idFilter, dFilterMinimum, dFilterMaximum, nFilterPoints, dBandwidth, dFilterParam, idDomain);
}
SignalFilter::SignalFilter (const char* szFilterName, const char* szDomainName, double dBandwidth, double dFilterParam)
- : m_adFilter(NULL), m_fail(false)
+: m_adFilter(NULL), m_fail(false)
{
m_nFilterPoints = 0;
m_dBandwidth = dBandwidth;
m_failMessage = " less than 2";
return;
}
-
+
m_nameFilter = convertFilterIDToName (m_idFilter);
m_nameDomain = convertDomainIDToName (m_idDomain);
m_nFilterPoints = nFilterPoints;
m_dBandwidth = dBandwidth;
m_dFilterMin = dFilterMinimum;
m_dFilterMax = dFilterMaximum;
-
+
m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
m_adFilter = new double [m_nFilterPoints];
-
+
if (m_idDomain == DOMAIN_FREQUENCY)
- createFrequencyFilter (m_adFilter);
+ createFrequencyFilter (m_adFilter);
else if (m_idDomain == DOMAIN_SPATIAL)
- createSpatialFilter (m_adFilter);
+ createSpatialFilter (m_adFilter);
}
SignalFilter::~SignalFilter (void)
{
- delete [] m_adFilter;
+ delete [] m_adFilter;
}
void
int center = (m_nFilterPoints - 1) / 2;
int sidelen = center;
m_adFilter[center] = 4. / (a * a);
-
+
for (int i = 1; i <= sidelen; i++ )
m_adFilter [center + i] = m_adFilter [center - i] = c / (4 * (i * i) - 1);
} else {
double x = m_dFilterMin;
for (int i = 0; i < m_nFilterPoints; i++, x += m_dFilterInc) {
if (haveAnalyticSpatial(m_idFilter))
- m_adFilter[i] = spatialResponseAnalytic (x);
+ m_adFilter[i] = spatialResponseAnalytic (x);
else
- m_adFilter[i] = spatialResponseCalc (x);
+ m_adFilter[i] = spatialResponseCalc (x);
}
}
}
SignalFilter::convertFilterNameToID (const char *filterName)
{
int filterID = FILTER_INVALID;
-
+
for (int i = 0; i < s_iFilterCount; i++)
if (strcasecmp (filterName, s_aszFilterName[i]) == 0) {
filterID = i;
break;
}
-
- return (filterID);
+
+ return (filterID);
}
const char *
SignalFilter::convertFilterIDToName (const int filterID)
{
static const char *name = "";
-
+
if (filterID >= 0 && filterID < s_iFilterCount)
- return (s_aszFilterName [filterID]);
-
+ return (s_aszFilterName [filterID]);
+
return (name);
}
SignalFilter::convertFilterIDToTitle (const int filterID)
{
static const char *title = "";
-
+
if (filterID >= 0 && filterID < s_iFilterCount)
- return (s_aszFilterTitle [filterID]);
-
+ return (s_aszFilterTitle [filterID]);
+
return (title);
}
-
+
int
SignalFilter::convertDomainNameToID (const char* const domainName)
{
int dID = DOMAIN_INVALID;
-
+
for (int i = 0; i < s_iDomainCount; i++)
- if (strcasecmp (domainName, s_aszDomainName[i]) == 0) {
+ if (strcasecmp (domainName, s_aszDomainName[i]) == 0) {
dID = i;
break;
}
-
- return (dID);
+
+ return (dID);
}
const char *
SignalFilter::convertDomainIDToName (const int domainID)
{
static const char *name = "";
-
+
if (domainID >= 0 && domainID < s_iDomainCount)
- return (s_aszDomainName [domainID]);
-
+ return (s_aszDomainName [domainID]);
+
return (name);
}
SignalFilter::convertDomainIDToTitle (const int domainID)
{
static const char *title = "";
-
+
if (domainID >= 0 && domainID < s_iDomainCount)
- return (s_aszDomainTitle [domainID]);
-
+ return (s_aszDomainTitle [domainID]);
+
return (title);
}
SignalFilter::response (double x)
{
double response = 0;
-
+
if (m_idDomain == DOMAIN_SPATIAL)
response = spatialResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam);
else if (m_idDomain == DOMAIN_FREQUENCY)
response = frequencyResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam);
-
+
return (response);
}
void
SignalFilter::copyFilterData (double* pdFilter, const int iStart, const int nPoints) const
{
- int iFirst = clamp (iStart, 0, m_nFilterPoints - 1);
- int iLast = clamp (iFirst + nPoints - 1, 0, m_nFilterPoints - 1);
-
- for (int i = iFirst; i <= iLast; i++)
- pdFilter[i - iFirst] = m_adFilter[i];
+ int iFirst = clamp (iStart, 0, m_nFilterPoints - 1);
+ int iLast = clamp (iFirst + nPoints - 1, 0, m_nFilterPoints - 1);
+
+ for (int i = iFirst; i <= iLast; i++)
+ pdFilter[i - iFirst] = m_adFilter[i];
}
/* NAME
- * filter_spatial_response_calc Calculate filter by discrete inverse fourier
- * transform of filters's frequency
- * response
- *
- * SYNOPSIS
- * y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
- * double y Filter's response in spatial domain
- * int filt_type Type of filter (definitions in ct.h)
- * double x Spatial position to evaluate filter
- * double m_bw Bandwidth of window
- * double param General parameter for various filters
- * int n Number of points to calculate integrations
- */
+* filter_spatial_response_calc Calculate filter by discrete inverse fourier
+* transform of filters's frequency
+* response
+*
+* SYNOPSIS
+* y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
+* double y Filter's response in spatial domain
+* int filt_type Type of filter (definitions in ct.h)
+* double x Spatial position to evaluate filter
+* double m_bw Bandwidth of window
+* double param General parameter for various filters
+* int n Number of points to calculate integrations
+*/
double
SignalFilter::spatialResponseCalc (double x) const
SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n)
{
double zmin, zmax;
-
+
if (filterID == FILTER_TRIANGLE) {
zmin = 0;
zmax = bw;
zmax = bw / 2;
}
double zinc = (zmax - zmin) / (n - 1);
-
+
double z = zmin;
double* q = new double [n];
for (int i = 0; i < n; i++, z += zinc)
double y = 2 * integrateSimpson (zmin, zmax, q, n);
delete q;
-
+
return (y);
}
/* NAME
- * filter_frequency_response Return filter frequency response
- *
- * SYNOPSIS
- * h = filter_frequency_response (filt_type, u, m_bw, param)
- * double h Filters frequency response at u
- * int filt_type Type of filter
- * double u Frequency to evaluate filter at
- * double m_bw Bandwidth of filter
- * double param General input parameter for various filters
- */
+* filter_frequency_response Return filter frequency response
+*
+* SYNOPSIS
+* h = filter_frequency_response (filt_type, u, m_bw, param)
+* double h Filters frequency response at u
+* int filt_type Type of filter
+* double u Frequency to evaluate filter at
+* double m_bw Bandwidth of filter
+* double param General input parameter for various filters
+*/
double
SignalFilter::frequencyResponse (double u) const
{
double q;
double au = fabs (u);
-
+ double abw = fabs (bw);
+
switch (filterID) {
case FILTER_BANDLIMIT:
- if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
+ if (au >= (abw / 2) + F_EPSILON)
q = 0.;
else
q = 1;
break;
case FILTER_ABS_BANDLIMIT:
- if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
+ if (au >= (abw / 2) + F_EPSILON)
q = 0.;
else
q = au;
break;
case FILTER_TRIANGLE:
- if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
+ if (au >= (abw / 2) + F_EPSILON)
q = 0;
else
- q = 1 - au / bw;
+ q = 1 - au / abw;
break;
case FILTER_COSINE:
- if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
+ if (au >= (abw / 2) + F_EPSILON)
q = 0;
else
- q = cos(PI * u / bw);
+ q = cos(PI * au / abw);
break;
case FILTER_ABS_COSINE:
- if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
+ if (au >= (abw / 2) + F_EPSILON)
q = 0;
else
- q = au * cos(PI * u / bw);
+ q = au * cos(PI * au / abw);
break;
case FILTER_SINC:
- q = bw * sinc (PI * bw * u, 1.);
+ q = abw * sinc (PI * abw * au, 1.);
break;
case FILTER_ABS_SINC:
- q = au * bw * sinc (PI * bw * u, 1.);
+ if (au >= (abw / 2) + F_EPSILON)
+ q = 0;
+ else
+ q = au * abw * sinc (PI * abw * au, 1.);
break;
case FILTER_HANNING:
- param = 0.54;
+ param = 0.5;
// follow through to G_HAMMING
case FILTER_G_HAMMING:
- if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
+ if (au >= (abw / 2) + F_EPSILON)
q = 0;
else
- q = param + (1 - param) * cos (TWOPI * u / bw);
+ q = param + (1 - param) * cos (TWOPI * au / abw);
break;
case FILTER_ABS_HANNING:
- param = 0.54;
+ param = 0.5;
// follow through to ABS_G_HAMMING
case FILTER_ABS_G_HAMMING:
- if (fabs(au) >= fabs(bw / 2) + F_EPSILON)
+ if (au >= (abw / 2) + F_EPSILON)
q = 0;
else
- q = au * (param + (1 - param) * cos(TWOPI * u / bw));
+ q = au * (param + (1 - param) * cos(TWOPI * au / abw));
break;
default:
q = 0;
sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
break;
}
+
return (q);
}
/* NAME
- * filter_spatial_response_analytic Calculate filter by analytic inverse fourier
- * transform of filters's frequency
- * response
- *
- * SYNOPSIS
- * y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
- * double y Filter's response in spatial domain
- * int filt_type Type of filter (definitions in ct.h)
- * double x Spatial position to evaluate filter
- * double m_bw Bandwidth of window
- * double param General parameter for various filters
- */
+* filter_spatial_response_analytic Calculate filter by analytic inverse fourier
+* transform of filters's frequency
+* response
+*
+* SYNOPSIS
+* y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
+* double y Filter's response in spatial domain
+* int filt_type Type of filter (definitions in ct.h)
+* double x Spatial position to evaluate filter
+* double m_bw Bandwidth of window
+* double param General parameter for various filters
+*/
double
SignalFilter::spatialResponseAnalytic (double x) const
SignalFilter::haveAnalyticSpatial (int filterID)
{
bool haveAnalytic = false;
-
+
switch (filterID) {
case FILTER_BANDLIMIT:
case FILTER_TRIANGLE:
default:
break;
}
-
+
return (haveAnalytic);
}
double w = bw / 2;
double b = PI / bw;
double b2 = TWOPI / bw;
-
+
switch (filterID) {
case FILTER_BANDLIMIT:
q = bw * sinc(u * w, 1.0);
q = sinc(b-u,w) + sinc(b+u,w);
break;
case FILTER_HANNING:
- param = 0.54;
+ param = 0.5;
// follow through to G_HAMMING
case FILTER_G_HAMMING:
q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
break;
case FILTER_ABS_HANNING:
- param = 0.54;
+ param = 0.5;
// follow through to ABS_G_HAMMING
case FILTER_ABS_G_HAMMING:
q = 2 * param * integral_abscos(u,w) +
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.61 2001/03/11 18:52:03 kevin Exp $
+** $Id: projections.cpp,v 1.62 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
{
init (iNViews, iNDets);
m_geometry = Scanner::GEOMETRY_EQUIANGULAR;
- m_dFanBeamAngle = iNDets * convertDegreesToRadians (3.06976 / 60);
m_dFocalLength = 510;
m_dSourceDetectorLength = 890;
m_detInc = convertDegreesToRadians (3.06976 / 60);
- m_detStart = -m_dFanBeamAngle / 2;
+ m_detStart = -(m_dFanBeamAngle / 2);
m_rotInc = TWOPI / static_cast<double>(iNViews);
m_rotStart = HALFPI;
+ m_dFanBeamAngle = (iNDets + 1) * m_detInc;
m_dViewDiameter = sin (m_dFanBeamAngle / 2) * m_dFocalLength * 2;
- if (iNDets != 1024)
- return false;
- bool bValid = (iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) || (iNViews == 1500 && lDataLength == 3120000);
- if (! bValid)
+ if (! ((iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L)
+ || (iNViews == 1500 && lDataLength == 3120000)))
return false;
+ int iCenter = (iNDets - 1) / 2; // change from (Nm+1)/2 because of 0 vs. 1 indexing
+ double* pdCosScale = new double [iNDets];
+ for (int i = 0; i < iNDets; i++)
+ pdCosScale[i] = cos ((i - iCenter) * m_detInc);
+
long lDataPos = 0;
for (int iv = 0; iv < iNViews; iv++) {
unsigned char* pArgBase = pData + lDataPos;
unsigned char* p = pArgBase+0;
-#ifndef WORDS_BIGENDIAN
- SwapBytes4 (p);
-#endif
+ SwapBytes4IfLittleEndian (p);
long lProjNumber = *reinterpret_cast<long*>(p);
p = pArgBase+20;
-#ifndef WORDS_BIGENDIAN
- SwapBytes4 (p);
-#endif
+ SwapBytes4IfLittleEndian (p);
long lEscale = *reinterpret_cast<long*>(p);
p = pArgBase+28;
-#ifndef WORDS_BIGENDIAN
- SwapBytes4 (p);
-#endif
+ SwapBytes4IfLittleEndian (p);
long lTime = *reinterpret_cast<long*>(p);
p = pArgBase + 4;
-#ifndef WORDS_BIGENDIAN
- SwapBytes4 (p);
-#endif
+ SwapBytes4IfLittleEndian (p);
double dAlpha = *reinterpret_cast<float*>(p) + HALFPI;
p = pArgBase+12;
-#ifndef WORDS_BIGENDIAN
- SwapBytes4 (p);
-#endif
+ SwapBytes4IfLittleEndian (p);
double dAlign = *reinterpret_cast<float*>(p);
p = pArgBase + 16;
-#ifndef WORDS_BIGENDIAN
- SwapBytes4 (p);
-#endif
+ SwapBytes4IfLittleEndian (p);
double dMaxValue = *reinterpret_cast<float*>(p);
- lDataPos += 32;
- double dEScale = pow (2.0, -lEscale);
- double dBetaInc = convertDegreesToRadians (3.06976 / 60);
- int iCenter = (iNDets + 1) / 2;
-
- DetectorArray& detArray = getDetectorArray( iv );
+ DetectorArray& detArray = getDetectorArray (iv);
detArray.setViewAngle (dAlpha);
DetectorValue* detval = detArray.detValues();
- double dTempScale = 2294.4871 * dEScale;
+ double dViewScale = 2294.4871 * ::pow (2.0, -lEscale);
+ lDataPos += 32;
for (int id = 0; id < iNDets; id++) {
int iV = pData[lDataPos+1] + 256 * pData[lDataPos];
if (iV > 32767) // two's complement signed conversion
iV = iV - 65536;
- double dCosScale = cos ((id + 1 - iCenter) * dBetaInc);
- detval[id] = iV / (dTempScale * dCosScale);
+ detval[id] = iV / (dViewScale * pdCosScale[id]);
lDataPos += 2;
}
}
+ delete pdCosScale;
return true;
}
if (nDet % 2 == 0) // even
pProjNew->m_detInc = m_dViewDiameter / (nDet - 1);
- ParallelRaysums parallel (this);
+ ParallelRaysums parallel (this, ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI);
double* pdThetaValuesForT = new double [pProjNew->nView()];
double* pdRaysumsForT = new double [pProjNew->nView()];
parallel.getThetaAndRaysumsForT (iD, pdThetaValuesForT, pdRaysumsForT);
double dViewAngle = m_rotStart;
+ int iLastFloor = -1;
for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues();
- detValues[iD] = parallel.interpolate (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), dViewAngle);
+ detValues[iD] = parallel.interpolate (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), dViewAngle, &iLastFloor);
}
}
delete pdThetaValuesForT;
pdDetValueCopy[i] = detValues[i];
double dDetPos = pProjNew->m_detStart;
+ int iLastFloor = -1;
for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) {
- detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos);
+ detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos, &iLastFloor);
}
}
delete pdDetValueCopy;
//
///////////////////////////////////////////////////////////////////////////////
-ParallelRaysums::ParallelRaysums (Projections* pProjections)
-: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet())
+ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange)
+: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()),
+ m_iThetaRange (iThetaRange)
{
int iGeometry = pProjections->geometry();
double dDetInc = pProjections->detInc();
ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++];
if (iGeometry == Scanner::GEOMETRY_PARALLEL) {
- pC->m_dTheta = normalizeAngle (dViewAngle);
+ pC->m_dTheta = dViewAngle;
pC->m_dT = dDetPos;
} else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) {
double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength());
- pC->m_dTheta = normalizeAngle (dViewAngle + dFanAngle);
+ pC->m_dTheta = dViewAngle + dFanAngle;
pC->m_dT = dFocalLength * sin(dFanAngle);
} else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
// fan angle is same as dDetPos
- pC->m_dTheta = normalizeAngle (dViewAngle + dDetPos);
+ pC->m_dTheta = dViewAngle + dDetPos;
pC->m_dT = dFocalLength * sin (dDetPos);
}
-#ifdef CONVERT_PARALLEL_PI
- if (pC->m_dTheta >= PI) { // convert T/Theta to 0-PI interval
- pC->m_dTheta -= PI;
- pC->m_dT = -pC->m_dT;
+ if (m_iThetaRange != THETA_RANGE_UNCONSTRAINED) {
+ pC->m_dTheta = normalizeAngle (pC->m_dTheta);
+ if (m_iThetaRange == THETA_RANGE_FOLD_TO_PI && pC->m_dTheta >= PI) {
+ pC->m_dTheta -= PI;
+ pC->m_dT = -pC->m_dT;
+ }
}
-#endif
pC->m_dRaysum = detValues[iD];
dDetPos += dDetInc;
}
}
// locate by bisection, O(log2(n))
+// iLastFloor is used when sequential calls to interpolate have monotonically increasing dX
double
-ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX)
+ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX, int* iLastFloor)
{
int iLower = -1;
int iUpper = n;
+ if (iLastFloor && *iLastFloor >= 0 && pdX[*iLastFloor] < dX)
+ iLower = *iLastFloor;
while (iUpper - iLower > 1) {
int iMiddle = (iUpper + iLower) >> 1;
return 0;
}
+ if (iLastFloor)
+ *iLastFloor = iLower;
return pdY[iLower] + (pdY[iUpper] - pdY[iLower]) * ((dX - pdX[iLower]) / (pdX[iUpper] - pdX[iLower]));
}
<pre>
<h1>Build Log</h1>
<h3>
---------------------Configuration: libctsim - Win32 Debug--------------------
+--------------------Configuration: ctsim - Win32 Debug--------------------
</h3>
<h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A7C.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD2.tmp" with contents
[
-/nologo /G6 /MTd /W3 /Gm /Gi /GR /GX /Zi /Od /Gy /I "..\..\..\wx2.2.5\src\png" /I "..\..\..\wx2.2.5\src\zlib" /I "..\..\INCLUDE" /I "..\..\getopt" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /I "..\..\..\wx2.2.5\include" /I "\dicom\ctn\include" /D "_DEBUG" /D "HAVE_WXWIN" /D "HAVE_STRING_H" /D "HAVE_GETOPT_H" /D "WIN32" /D "_MBCS" /D "_LIB" /D "MSVC" /D "HAVE_FFTW" /D "HAVE_PNG" /D "HAVE_SGP" /D "HAVE_WXWINDOWS" /D "__WXMSW__" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "HAVE_CTN_DICOM" /D VERSION=\"3.1.0\" /FR"Debug/" /Fp"Debug/libctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c
-"D:\ctsim\libctsim\ctndicom.cpp"
+/nologo /G6 /MTd /W3 /Gm /Gi /GR /GX /Zi /Od /Gy /I "\wx2.2.5\include" /I "..\..\..\fftw-2.1.3\fftw" /I "\wx2.2.5\src\png" /I "\wx2.2.5\src\zlib" /I "..\..\include" /I "..\..\getopt" /I "..\..\..\fftw-2.1.3\rfftw" /I "\dicom\ctn\include" /D VERSION=\"3.0.0beta1\" /D "_DEBUG" /D "__WXMSW__" /D "HAVE_SGP" /D "HAVE_PNG" /D "HAVE_WXWINDOWS" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "HAVE_STRING_H" /D "HAVE_FFTW" /D "HAVE_RFFTW" /D "HAVE_GETOPT_H" /D "MSVC" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "STRICT" /D CTSIMVERSION=\"3.1.0\" /D "HAVE_CTN_DICOM" /D CTSIMVERSION=\"3.0.0alpha5\" /FR"Debug/" /Fp"Debug/ctsim.pch" /YX"ctsim.h" /Fo"Debug/" /Fd"Debug/" /FD /GZ /c
+"C:\ctsim\src\docs.cpp"
]
-Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A7C.tmp"
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A7D.tmp" with contents
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD2.tmp"
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD3.tmp" with contents
[
-/nologo /out:"Debug\libctsim.lib"
-.\Debug\array2dfile.obj
-.\Debug\backprojectors.obj
-.\Debug\clip.obj
-.\Debug\consoleio.obj
-.\Debug\ctndicom.obj
-.\Debug\dlgezplot.obj
-.\Debug\ezplot.obj
-.\Debug\ezset.obj
-.\Debug\ezsupport.obj
-.\Debug\filter.obj
-.\Debug\fnetorderstream.obj
-.\Debug\fourier.obj
-.\Debug\getopt.obj
-.\Debug\getopt1.obj
-.\Debug\globalvars.obj
-.\Debug\hashtable.obj
-.\Debug\imagefile.obj
-.\Debug\interpolator.obj
-.\Debug\mathfuncs.obj
-.\Debug\phantom.obj
-.\Debug\plotfile.obj
-.\Debug\pol.obj
-.\Debug\procsignal.obj
-.\Debug\projections.obj
-.\Debug\reconstruct.obj
-.\Debug\scanner.obj
-.\Debug\sgp.obj
-.\Debug\strfuncs.obj
-.\Debug\syserror.obj
-.\Debug\trace.obj
-.\Debug\transformmatrix.obj
-.\Debug\xform.obj
+/nologo /G6 /MTd /W3 /Gm /Gi /GR /GX /Zi /Od /Gy /I "\wx2.2.5\include" /I "..\..\..\fftw-2.1.3\fftw" /I "\wx2.2.5\src\png" /I "\wx2.2.5\src\zlib" /I "..\..\include" /I "..\..\getopt" /I "..\..\..\fftw-2.1.3\rfftw" /I "\dicom\ctn\include" /D VERSION=\"3.0.0beta1\" /D "_DEBUG" /D "__WXMSW__" /D "HAVE_SGP" /D "HAVE_PNG" /D "HAVE_WXWINDOWS" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "HAVE_STRING_H" /D "HAVE_FFTW" /D "HAVE_RFFTW" /D "HAVE_GETOPT_H" /D "MSVC" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "STRICT" /D CTSIMVERSION=\"3.1.0\" /D "HAVE_CTN_DICOM" /FR"Debug/" /Fp"Debug/ctsim.pch" /YX"ctsim.h" /Fo"Debug/" /Fd"Debug/" /FD /GZ /c
+"C:\ctsim\src\graph3dview.cpp"
]
-Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A7D.tmp"
-<h3>Output Window</h3>
-Compiling...
-ctndicom.cpp
-Creating library...
-<h3>
---------------------Configuration: ctsim - Win32 Debug--------------------
-</h3>
-<h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A7E.tmp" with contents
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD3.tmp"
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD4.tmp" with contents
[
winmm.lib rpcrt4.lib ws2_32.lib ../libctsim/Debug/libctsim.lib libcmtd.lib ..\..\..\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib ..\..\..\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib wxd.lib xpmd.lib tiffd.lib zlibd.lib pngd.lib comctl32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib opengl32.lib glu32.lib htmlhelp.lib ctn_lib.lib /nologo /subsystem:windows /incremental:yes /pdb:"Debug/ctsim.pdb" /debug /machine:I386 /out:"Debug/ctsim.exe" /pdbtype:sept /libpath:"\wx2.2.5\lib" /libpath:"\dicom\ctn\winctn\ctn_lib\Debug"
.\Debug\backgroundmgr.obj
\wx2.2.5\lib\zlibd.lib
\wx2.2.5\lib\tiffd.lib
]
-Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A7E.tmp"
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD4.tmp"
<h3>Output Window</h3>
+Compiling...
+docs.cpp
+Compiling...
+graph3dview.cpp
Linking...
** This is part of the CTSim program
** Copyright (C) 1983-2001 Kevin Rosenberg
**
-** $Id: backgroundmgr.cpp,v 1.18 2001/03/11 17:55:29 kevin Exp $
+** $Id: backgroundmgr.cpp,v 1.19 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
resizeWindow();
if (m_iNumTasks == 1) {
m_pCanvas->SetFocus();
+ pGauge->SetFocus();
Show(true);
}
}
#define IDH_DLG_COMPARISON 8611
#define IDH_DLG_PREFERENCES 8612
#define IDH_DLG_POLAR 8613
+
+// Need to add to .tex file
#define IDH_DLG_IMPORT 8614
+#define IDH_DLG_THETA_RANGE 8615
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: dialogs.cpp,v 1.51 2001/03/11 17:55:29 kevin Exp $
+** $Id: dialogs.cpp,v 1.52 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
#include "docs.h"
#include "views.h"
#include "imagefile.h"
+#include "projections.h"
#if defined(MSVC) || HAVE_SSTREAM
#include <sstream>
}
+///////////////////////////////////////////////////////////////////////
+// CLASS IMPLEMENTATION
+// DialogGetThetaRange
+///////////////////////////////////////////////////////////////////////
+
+DialogGetThetaRange::DialogGetThetaRange (wxWindow* pParent, int iDefaultThetaRange)
+: wxDialog (pParent, -1, _T("Select Phantom"), wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION)
+{
+ wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL);
+
+ pTopSizer->Add (new wxStaticText (this, -1, "Select Theta Range"), 0, wxCENTER | wxALL, 5);
+
+ pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);
+
+ wxString asTitle[] = {"Unconstrained", "Normalized to 2pi", "Fold to pi"};
+
+ m_pRadioBoxThetaRange = new wxRadioBox (this, -1, _T("Theta Range"), wxDefaultPosition, wxDefaultSize, 3, asTitle, 1, wxRA_SPECIFY_COLS);
+ if (iDefaultThetaRange == ParallelRaysums::THETA_RANGE_UNCONSTRAINED)
+ m_pRadioBoxThetaRange->SetSelection (0);
+ else if (iDefaultThetaRange == ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI)
+ m_pRadioBoxThetaRange->SetSelection (1);
+ else if (iDefaultThetaRange == ParallelRaysums::THETA_RANGE_FOLD_TO_PI)
+ m_pRadioBoxThetaRange->SetSelection (2);
+
+ pTopSizer->Add (m_pRadioBoxThetaRange, 0, wxALL | wxALIGN_CENTER);
+
+ pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5);
+
+ wxBoxSizer* pButtonSizer = new wxBoxSizer (wxHORIZONTAL);
+ wxButton* pButtonOk = new wxButton (this, wxID_OK, "Okay");
+ pButtonSizer->Add (pButtonOk, 0, wxEXPAND | wxALL, 10);
+ wxButton* pButtonCancel = new wxButton (this, wxID_CANCEL, "Cancel");
+ pButtonSizer->Add (pButtonCancel, 0, wxEXPAND | wxALL, 10);
+ CTSimHelpButton* pButtonHelp = new CTSimHelpButton (this, IDH_DLG_THETA_RANGE);
+ pButtonSizer->Add (pButtonHelp, 0, wxEXPAND | wxALL, 10);
+
+ pTopSizer->Add (pButtonSizer, 0, wxALIGN_CENTER);
+ pButtonOk->SetDefault();
+
+ SetAutoLayout (true);
+ SetSizer (pTopSizer);
+ pTopSizer->Fit (this);
+ pTopSizer->SetSizeHints (this);
+}
+
+int
+DialogGetThetaRange::getThetaRange()
+{
+ int iSelection = m_pRadioBoxThetaRange->GetSelection();
+ if (iSelection == 0)
+ return ParallelRaysums::THETA_RANGE_UNCONSTRAINED;
+ else if (iSelection == 1)
+ return ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI;
+ else
+ return ParallelRaysums::THETA_RANGE_FOLD_TO_PI;
+}
+
+
///////////////////////////////////////////////////////////////////////
// CLASS IMPLEMENTATION
// DialogGetComparisonImage
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: dialogs.h,v 1.34 2001/03/11 15:27:30 kevin Exp $
+** $Id: dialogs.h,v 1.35 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
StringValueAndTitleRadioBox* m_pRadioBoxPhantom;
};
+class DialogGetThetaRange : public wxDialog
+{
+ public:
+ DialogGetThetaRange (wxWindow* pParent, int iDefaultThetaRange = ParallelRaysums::THETA_RANGE_UNCONSTRAINED);
+ virtual ~DialogGetThetaRange () {}
+
+ int getThetaRange ();
+
+ private:
+ wxRadioBox* m_pRadioBoxThetaRange;
+};
+
#include <vector>
class ImageFileDocument;
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: docs.cpp,v 1.35 2001/03/11 18:52:03 kevin Exp $
+** $Id: docs.cpp,v 1.36 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
Graph3dFileDocument::~Graph3dFileDocument()
{
-// delete [] m_pVertices;
-// delete [] m_pNormals;
}
bool
bool
Graph3dFileDocument::createFromImageFile (const ImageFile& rImageFile)
{
-// delete [] m_pVertices;
-// delete [] m_pNormals;
-
-
m_nx = rImageFile.nx();
m_ny = rImageFile.ny();
m_array = rImageFile.getArray();
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: graph3dview.cpp,v 1.18 2001/03/08 17:37:09 kevin Exp $
+** $Id: graph3dview.cpp,v 1.19 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
-#if 1
+#if 0
GLfloat eyep[3], lookp[3], up[3];
eyep[0] = -nx/2; eyep[1] = 0; eyep[2] = -ny/2;
lookp[0] = 0; lookp[1] = 0, lookp[2] = 0;
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: views.cpp,v 1.135 2001/03/11 18:52:03 kevin Exp $
+** $Id: views.cpp,v 1.136 2001/03/13 04:44:25 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
edit_menu->Append(IFMENU_EDIT_COPY, "Copy\tCtrl-C");
edit_menu->Append(IFMENU_EDIT_CUT, "Cut\tCtrl-X");
edit_menu->Append(IFMENU_EDIT_PASTE, "Paste\tCtrl-V");
-
+
wxMenu *view_menu = new wxMenu;
view_menu->Append(IFMENU_VIEW_SCALE_MINMAX, "Display Scale S&et...\tCtrl-E");
view_menu->Append(IFMENU_VIEW_SCALE_AUTO, "Display Scale &Auto...\tCtrl-A");
ImageFileView::OnEditCopy (wxCommandEvent& event)
{
wxBitmapDataObject *pBitmapObject = new wxBitmapDataObject;
-
+
pBitmapObject->SetBitmap (m_bitmap);
-
+
if (wxTheClipboard->Open()) {
wxTheClipboard->SetData (pBitmapObject);
wxTheClipboard->Close();
ImageFileView::OnEditPaste (wxCommandEvent& event)
{
ImageFile& rIF = GetDocument()->getImageFile();
-
+
if (wxTheClipboard->Open()) {
wxBitmap bitmap;
if (wxTheClipboard->IsSupported (wxDF_BITMAP)) {
bitmap = bitmapObject.GetBitmap ();
}
wxTheClipboard->Close();
-
+
int nx = rIF.nx();
int ny = rIF.ny();
if (bitmap.Ok() == true && bitmap.GetWidth() == nx && bitmap.GetHeight() == ny) {
pScaledDoc->Modify (true);
pScaledDoc->UpdateAllViews (this);
pScaledDoc->getView()->getFrame()->Show(true);
- pScaledDoc->Activate();
+ pScaledDoc->Activate();
}
}
GetDocumentManager()->ActivateView (pGraph3d->getView(), true, false);
::wxYield();
pGraph3d->getView()->getCanvas()->SetFocus();
- pGraph3d->Activate();
}
#endif
pPlotDoc->Modify (true);
pPlotDoc->getView()->getFrame()->Show(true);
pPlotDoc->UpdateAllViews ();
- pPlotDoc->Activate();
+ pPlotDoc->Activate();
}
}
pPlotDoc->Modify (true);
pPlotDoc->getView()->getFrame()->Show(true);
pPlotDoc->UpdateAllViews ();
- pPlotDoc->Activate();
+ pPlotDoc->Activate();
}
}
pPlotDoc->Modify (true);
pPlotDoc->getView()->getFrame()->Show(true);
pPlotDoc->UpdateAllViews ();
- pPlotDoc->Activate();
+ pPlotDoc->Activate();
}
}
pPlotDoc->Modify (true);
pPlotDoc->getView()->getFrame()->Show(true);
pPlotDoc->UpdateAllViews ();
- pPlotDoc->Activate();
+ pPlotDoc->Activate();
}
}
#endif
return;
} else
#endif // HAVE_WXTHREADS
- {
+ {
pProj = new Projections;
pProj->initFromScanner (theScanner);
wxProgressDialog dlgProgress (wxString("Projection"), wxString("Projection Progress"), pProj->nView() + 1, getFrameForChild(), wxPD_CAN_ABORT );
os << "Rasterize Phantom " << rPhantom.name() << ": XSize=" << m_iDefaultRasterNX << ", YSize="
<< m_iDefaultRasterNY << ", ViewRatio=" << m_dDefaultRasterViewRatio << ", nSamples="
<< m_iDefaultRasterNSamples;;
-
+
#if HAVE_WXTHREADS
if (theApp->getUseBackgroundTasks()) {
RasterizerSupervisorThread* pThread = new RasterizerSupervisorThread (this, m_iDefaultRasterNX, m_iDefaultRasterNY,
return;
}
}
-
+
ImageFileDocument* pRasterDoc = theApp->newImageDoc();
if (! pRasterDoc) {
sys_error (ERR_SEVERE, "Unable to create image file");
rasterView->getFrame()->SetFocus();
rasterView->OnUpdate (rasterView, NULL);
}
- pRasterDoc->Activate();
+ pRasterDoc->Activate();
}
}
*theApp->getLog() << "Error converting to Polar\n";
return;
}
-
+
pPolarDoc = theApp->newImageDoc ();
if (! pPolarDoc) {
sys_error (ERR_SEVERE, "Unable to create image file");
pPolarDoc->Modify (true);
pPolarDoc->getView()->getFrame()->Show(true);
pPolarDoc->UpdateAllViews ();
- pPolarDoc->Activate();
+ pPolarDoc->Activate();
}
}
pPolarDoc->Modify (true);
pPolarDoc->getView()->getFrame()->Show(true);
pPolarDoc->UpdateAllViews ();
- pPolarDoc->Activate();
+ pPolarDoc->Activate();
}
}
void
ProjectionFileView::OnPlotTThetaSampling (wxCommandEvent& event)
{
- Projections& rProj = GetDocument()->getProjections();
- ParallelRaysums parallel (&rProj);
- PlotFileDocument* pPlotDoc = theApp->newPlotDoc();
- PlotFile& rPlot = pPlotDoc->getPlotFile();
- ParallelRaysums::CoordinateContainer& coordContainer = parallel.getCoordinates();
- double* pdT = new double [parallel.getNumCoordinates()];
- double* pdTheta = new double [parallel.getNumCoordinates()];
-
- for (int i = 0; i < parallel.getNumCoordinates(); i++) {
- pdT[i] = coordContainer[i]->m_dT;
- pdTheta[i] = coordContainer[i]->m_dTheta;
- }
- rPlot.setCurveSize (2, parallel.getNumCoordinates(), true);
- rPlot.addEzsetCommand ("title T-Theta Sampling");
- rPlot.addEzsetCommand ("xlabel T");
- rPlot.addEzsetCommand ("ylabel Theta");
- rPlot.addEzsetCommand ("curve 1");
- if (rProj.nDet() < 50 && rProj.nView() < 50)
- rPlot.addEzsetCommand ("symbol 1"); // x symbol
- else
- rPlot.addEzsetCommand ("symbol 6"); // point symbol
- rPlot.addEzsetCommand ("noline");
- rPlot.addColumn (0, pdT);
- rPlot.addColumn (1, pdTheta);
- delete pdT;
- delete pdTheta;
- if (theApp->getAskDeleteNewDocs())
- pPlotDoc->Modify (true);
- pPlotDoc->getView()->getFrame()->Show(true);
- pPlotDoc->UpdateAllViews ();
- pPlotDoc->Activate();
+ DialogGetThetaRange dlgTheta (this->getFrame(), ParallelRaysums::THETA_RANGE_UNCONSTRAINED);
+ if (dlgTheta.ShowModal() != wxID_OK)
return;
+
+ int iThetaRange = dlgTheta.getThetaRange();
+
+ Projections& rProj = GetDocument()->getProjections();
+ ParallelRaysums parallel (&rProj, iThetaRange);
+ PlotFileDocument* pPlotDoc = theApp->newPlotDoc();
+ PlotFile& rPlot = pPlotDoc->getPlotFile();
+ ParallelRaysums::CoordinateContainer& coordContainer = parallel.getCoordinates();
+ double* pdT = new double [parallel.getNumCoordinates()];
+ double* pdTheta = new double [parallel.getNumCoordinates()];
+
+ for (int i = 0; i < parallel.getNumCoordinates(); i++) {
+ pdT[i] = coordContainer[i]->m_dT;
+ pdTheta[i] = coordContainer[i]->m_dTheta;
+ }
+ rPlot.setCurveSize (2, parallel.getNumCoordinates(), true);
+ rPlot.addEzsetCommand ("title T-Theta Sampling");
+ rPlot.addEzsetCommand ("xlabel T");
+ rPlot.addEzsetCommand ("ylabel Theta");
+ rPlot.addEzsetCommand ("curve 1");
+ if (rProj.nDet() < 50 && rProj.nView() < 50)
+ rPlot.addEzsetCommand ("symbol 1"); // x symbol
+ else
+ rPlot.addEzsetCommand ("symbol 6"); // point symbol
+ rPlot.addEzsetCommand ("noline");
+ rPlot.addColumn (0, pdT);
+ rPlot.addColumn (1, pdTheta);
+ delete pdT;
+ delete pdTheta;
+ if (theApp->getAskDeleteNewDocs())
+ pPlotDoc->Modify (true);
+ pPlotDoc->getView()->getFrame()->Show(true);
+ pPlotDoc->UpdateAllViews ();
+ pPlotDoc->Activate();
}
void
Projections* pProjNew = rProj.interpolateToParallel();
ProjectionFileDocument* pProjDocNew = theApp->newProjectionDoc();
pProjDocNew->setProjections (pProjNew);
-
+
if (ProjectionFileView* projView = pProjDocNew->getView()) {
projView->OnUpdate (projView, NULL);
if (projView->getCanvas())
defaultROI.m_dXMax = defaultROI.m_dXMin + rProj.phmLen();
defaultROI.m_dYMin = -rProj.phmLen() / 2;
defaultROI.m_dYMax = defaultROI.m_dYMin + rProj.phmLen();
-
+
DialogGetReconstructionParameters dialogReconstruction (getFrameForChild(), m_iDefaultNX, m_iDefaultNY,
m_iDefaultFilter, m_dDefaultFilterParam, m_iDefaultFilterMethod, m_iDefaultFilterGeneration,
m_iDefaultZeropad, m_iDefaultInterpolation, m_iDefaultInterpParam, m_iDefaultBackprojector,
m_iDefaultBackprojector = Backprojector::convertBackprojectNameToID (optBackprojectName.c_str());
m_iDefaultTrace = dialogReconstruction.getTrace();
dialogReconstruction.getROI (&defaultROI);
-
+
if (m_iDefaultNX <= 0 && m_iDefaultNY <= 0)
return;
return;
} else
#endif
- {
+ {
pImageFile = new ImageFile (m_iDefaultNX, m_iDefaultNY);
Reconstructor* pReconstructor = new Reconstructor (rProj, *pImageFile, optFilterName.c_str(),
m_dDefaultFilterParam, optFilterMethodName.c_str(), m_iDefaultZeropad, optFilterGenerationName.c_str(),
convert_menu->Append (PJMENU_CONVERT_FFT_POLAR, "&FFT->Polar Image...\tCtrl-M");
convert_menu->AppendSeparator();
convert_menu->Append (PJMENU_CONVERT_PARALLEL, "&Interpolate to Parallel");
-
+
wxMenu* filter_menu = new wxMenu;
filter_menu->Append (PJMENU_ARTIFACT_REDUCTION, "&Artifact Reduction");
-
+
wxMenu* analyze_menu = new wxMenu;
analyze_menu->Append (PJMENU_PLOT_TTHETA_SAMPLING, "&Plot T-Theta Sampling\tCtrl-T");
-
+
wxMenu *reconstruct_menu = new wxMenu;
reconstruct_menu->Append (PJMENU_RECONSTRUCT_FBP, "&Filtered Backprojection...\tCtrl-R", "Reconstruct image using filtered backprojection");
reconstruct_menu->Append (PJMENU_RECONSTRUCT_FOURIER, "&Fourier...\tCtrl-E", "Reconstruct image using inverse Fourier");
const int iNColumns = rPlotFile.getNumColumns();
const int iNRecords = rPlotFile.getNumRecords();
const bool bScatterPlot = rPlotFile.getIsScatterPlot();
-
+
if (iNColumns > 0 && iNRecords > 0) {
if (m_pEZPlot)
delete m_pEZPlot;
m_pEZPlot->ezset("box");
m_pEZPlot->ezset("grid");
-
+
double* pdX = new double [iNRecords];
double* pdY = new double [iNRecords];
if (! bScatterPlot) {
rPlotFile.getColumn (0, pdX);
-
+
for (int iCol = 1; iCol < iNColumns; iCol++) {
rPlotFile.getColumn (iCol, pdY);
m_pEZPlot->addCurve (pdX, pdY, iNRecords);
projects / ctsim.git / commitdiff