** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: procsignal.h,v 1.15 2001/03/01 07:30:49 kevin Exp $
+** $Id: procsignal.h,v 1.16 2001/03/13 14:53:43 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
static void finiteFourierTransform (const std::complex<double> input[], std::complex<double> output[], const int n, const int direction);
static void finiteFourierTransform (const std::complex<double> input[], double output[], const int n, const int direction);
+ static int addZeropadFactor (int n, int iZeropad);
private:
std::string m_nameFilterMethod;
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.h,v 1.33 2001/03/13 08:24:41 kevin Exp $
+** $Id: projections.h,v 1.34 2001/03/13 14:53:44 kevin Exp $
**
**
** This program is free software; you can redistribute it and/or modify
bool convertPolar (ImageFile& rIF, int iInterpolation);
bool convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad);
- bool calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet);
+ void calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet,
+ int iNumDetWithZeros, double dZeropadRatio);
void interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue,
- double** ppdDet, double** ppdView, unsigned int nView, unsigned int nDet, int iInterpolate);
+ double** ppdDet, double** ppdView, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros,
+ int iInterpolate);
bool reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* frequencyFilterName, const char* const interpName, int interp_param, const char* const backprojName, const int trace) const;
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: ctndicom.cpp,v 1.11 2001/03/13 10:35:06 kevin Exp $
+** $Id: ctndicom.cpp,v 1.12 2001/03/13 14:53:44 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
char szIDImageType[] = "ORIGINAL";
char szIDManufacturerModel[65] = "";
- std::ostringstream osComments;
- m_pImageFile->printLabelsBrief (osComments);
- size_t sizeIDComments = osComments.str().length();
- char* pszIDComments = new char [sizeIDComments+1];
- strncpy (pszIDComments, osComments.str().c_str(), sizeIDComments);
+ std::ostringstream osPatComments;
+ m_pImageFile->printLabelsBrief (osPatComments);
+ size_t sizePatComments = osPatComments.str().length();
+ char* pszPatComments = new char [sizePatComments+1];
+ strncpy (pszPatComments, osPatComments.str().c_str(), sizePatComments);
snprintf (szIDSOPInstanceUID, sizeof(szIDSOPInstanceUID), "%s.2.1.6.1", szCTSimRoot);
snprintf (szRelStudyInstanceUID, sizeof(szRelStudyInstanceUID), "%s.2.1.6.1.1", szCTSimRoot);
{DCM_IDMODALITY, DCM_CS, "", 1, strlen(szModality), szModality},
{DCM_IDSOPCLASSUID, DCM_UI, "", 1, strlen(szSOPClassUID), szSOPClassUID},
{DCM_IDMANUFACTURERMODEL, DCM_LO, "", 1, strlen(szIDManufacturerModel), szIDManufacturerModel},
- {DCM_PATCOMMENTS, DCM_LT, "", 1, strlen(pszIDComments), pszIDComments},
+ {DCM_PATCOMMENTS, DCM_LT, "", 1, strlen(pszPatComments), pszPatComments},
};
int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
cond = DCM_ModifyElements (&m_pObject, &elemPixelData, 1, NULL, 0, &iUpdateCount);
delete pRawPixels;
- delete pszIDComments;
+ delete pszPatComments;
if (cond != DCM_NORMAL || iUpdateCount != 1) {
m_bFail = true;
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: procsignal.cpp,v 1.27 2001/03/01 07:30:49 kevin Exp $
+** $Id: procsignal.cpp,v 1.28 2001/03/13 14:53:44 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
if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
// calculate number of filter points with zeropadding
- m_nFilterPoints = m_nSignalPoints;
- if (m_iZeropad > 0) {
- double logBase2 = log(m_nFilterPoints) / log(2);
- int nextPowerOf2 = static_cast<int>(floor(logBase2));
- if (logBase2 != floor(logBase2))
- nextPowerOf2++;
- nextPowerOf2 += (m_iZeropad - 1);
- m_nFilterPoints = 1 << nextPowerOf2;
-#if defined(DEBUG) || defined(_DEBUG)
- if (m_traceLevel >= Trace::TRACE_CONSOLE)
- sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints);
-#endif
- }
+ m_nFilterPoints = addZeropadFactor (m_nSignalPoints, m_iZeropad);
m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
if (m_nFilterPoints % 2) { // Odd
dlgEZPlot.ShowModal();
}
#endif
-
+
// This works fairly well. I'm not sure why since scaling for geometries is done on
// frequency filter rather than spatial filter as it should be.
// It gives values slightly off than freq/inverse filtering
dlgEZPlot.ShowModal();
}
#endif
-
+
// FILTERING: FREQUENCY - INVERSE FOURIER
-
+
} else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
// calculate number of filter points with zeropadding
int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
}
}
+
+int
+ProcessSignal::addZeropadFactor (int n, int iZeropad)
+{
+ if (iZeropad > 0) {
+ double dLogBase2 = log(n) / log(2);
+ int iLogBase2 = static_cast<int>(floor (dLogBase2));
+ if (dLogBase2 != floor(dLogBase2))
+ iLogBase2++; // raise up to next power of 2
+ n = 1 << (iLogBase2 + (iZeropad - 1));
+ }
+
+ return n;
+}
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.65 2001/03/13 10:35:06 kevin Exp $
+** $Id: projections.cpp,v 1.66 2001/03/13 14:53:44 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
double** ppdView = adView.getArray();
double** ppdDet = adDet.getArray();
- if (! pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet))
- return false;
+ pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, m_nDet, 1.);
std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
unsigned int iView;
for (iView = 0; iView < m_nView; iView++) {
ppcDetValue[iView] = new std::complex<double> [m_nDet];
+ DetectorValue* detval = pProj->getDetectorArray (iView).detValues();
for (unsigned int iDet = 0; iDet < m_nDet; iDet++)
- ppcDetValue[iView][iDet] = std::complex<double>(pProj->getDetectorArray (iView).detValues()[iDet], 0);
+ ppcDetValue[iView][iDet] = std::complex<double>(detval[iDet], 0);
}
- pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, iInterpolationID);
+ pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, pProj->m_nDet, iInterpolationID);
for (iView = 0; iView < m_nView; iView++)
delete [] ppcDetValue[iView];
bool
Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
{
+#ifndef HAVE_FFT
+ return false;
+#else
unsigned int nx = rIF.nx();
unsigned int ny = rIF.ny();
ImageFileArray v = rIF.getArray();
return false;
}
-#ifndef HAVE_FFT
- return false;
-#else
Array2d<double> adView (nx, ny);
Array2d<double> adDet (nx, ny);
double** ppdView = adView.getArray();
double** ppdDet = adDet.getArray();
std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
- unsigned int iView;
- double* pdDet = new double [m_nDet];
- fftw_complex* pcIn = new fftw_complex [m_nDet];
- fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE);
- for (iView = 0; iView < m_nView; iView++) {
- unsigned int iDet;
- for (iDet = 0; iDet < m_nDet; iDet++) {
- pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet];
+ int iNumDetWithZeros = ProcessSignal::addZeropadFactor (m_nDet, iZeropad);
+ double dZeropadRatio = iNumDetWithZeros / static_cast<double>(m_nDet);
+
+ double* pdDet = new double [iNumDetWithZeros];
+ fftw_complex* pcIn = new fftw_complex [iNumDetWithZeros];
+ fftw_plan plan = fftw_create_plan (iNumDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE);
+
+ for (unsigned int iView = 0; iView < m_nView; iView++) {
+ DetectorValue* detval = getDetectorArray(iView).detValues();
+ for (unsigned int iDet = 0; iDet < m_nDet; iDet++) {
+ pcIn[iDet].re = detval[iDet];
pcIn[iDet].im = 0;
}
+ for (unsigned int iDet2 = m_nDet; iDet2 < iNumDetWithZeros; iDet2++)
+ pcIn[iDet2].re = pcIn[iDet2].im = 0;
+
fftw_one (plan, pcIn, NULL);
- ppcDetValue[iView] = new std::complex<double> [m_nDet];
- for (iDet = 0; iDet < m_nDet; iDet++)
- ppcDetValue[iView][iDet] = std::complex<double> (pcIn[iDet].re, pcIn[iDet].im);
- Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet);
+ ppcDetValue[iView] = new std::complex<double> [iNumDetWithZeros];
+ for (unsigned int iD = 0; iD < iNumDetWithZeros; iD++)
+ ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re, pcIn[iD].im);
+//ppcDetValue[iView][iD] = std::complex<double> (std::abs(std::complex<double>(pcIn[iD].re, pcIn[iD].im)), 0);
+ Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumDetWithZeros);
}
fftw_destroy_plan (plan);
delete [] pcIn;
- bool bError = calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
-
- if (! bError)
- interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID);
+ calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumDetWithZeros, dZeropadRatio);
+ interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumDetWithZeros,
+ iInterpolationID);
- for (iView = 0; iView < m_nView; iView++)
- delete [] ppcDetValue[iView];
+ for (int i = 0; i < m_nView; i++)
+ delete [] ppcDetValue[i];
delete [] ppcDetValue;
- return bError;
+ return true;
#endif
}
-bool
-Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet)
+void
+Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet,
+ int iNumDetWithZeros, double dZeropadRatio)
{
double xMin = -phmLen() / 2;
double xMax = xMin + phmLen();
double yMin = -phmLen() / 2;
double yMax = yMin + phmLen();
-
+ double xCent = (xMin + xMax) / 2;
+ double yCent = (yMin + yMax) / 2;
+
+ xMin = (xMin - xCent) * dZeropadRatio + xCent;
+ xMax = (xMax - xCent) * dZeropadRatio + xCent;
+ yMin = (yMin - yCent) * dZeropadRatio + yCent;
+ yMax = (yMax - yCent) * dZeropadRatio + yCent;
+
double xInc = (xMax - xMin) / nx; // size of cells
double yInc = (yMax - yMin) / ny;
-
- int iDetCenter = (m_nDet - 1) / 2; // index refering to L=0 projection
+
+ // +1 is correct for frequency data, ndet-1 is correct for projections
+ int iDetCenter = (iNumDetWithZeros - 1) / 2; // index refering to L=0 projection
+ if (iNumDetWithZeros % 2 == 0)
+ iDetCenter = (iNumDetWithZeros + 1) / 2;
// Calculates polar coordinates (view#, det#) for each point on phantom grid
double x = xMin + xInc / 2; // Rectang coords of center of pixel
double r = ::sqrt (x * x + y * y);
double phi = atan2 (y, x);
+ if (phi < 0)
+ phi += TWOPI;
+
if (phi >= PI) {
phi -= PI;
- } else if (phi < 0) {
- phi += PI;
- } else
r = -r;
+ }
ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc;
ppdDet[ix][iy] = (r / m_detInc) + iDetCenter;
}
}
-
- return true;
}
void
Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
- unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue,
- double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, int iInterpolationID)
+ unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue, double** ppdView,
+ double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID)
{
+
for (unsigned int ix = 0; ix < ny; ix++) {
for (unsigned int iy = 0; iy < ny; iy++) {
if (iInterpolationID == POLAR_INTERP_NEAREST) {
iView = 0;
// iDet = m_nDet - iDet;
}
- if (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) {
+ if (iDet >= 0 && iDet < nDetWithZeros && iView >= 0 && iView < nView) {
v[ix][iy] = ppcDetValue[iView][iDet].real();
if (vImag)
vImag[ix][iy] = ppcDetValue[iView][iDet].imag();
v2 = ppcDetValue[iFloorView + 1][iFloorDet];
else
v2 = ppcDetValue[0][iFloorDet];
- if (iFloorDet < nDet - 1)
+ if (iFloorDet < nDetWithZeros - 1)
v4 = ppcDetValue[iFloorView][iFloorDet+1];
else
v4 = v1;
- if (iFloorView < nView - 1 && iFloorDet < nDet - 1)
+ if (iFloorView < nView - 1 && iFloorDet < nDetWithZeros - 1)
v3 = ppcDetValue [iFloorView+1][iFloorDet+1];
else if (iFloorView < nView - 1)
v3 = v2;
###############################################################################
-Project: "IsosurfVC"=..\..\wx2.2.5\samples\opengl\isosurf\IsosurfVC.dsp - Package Owner=<4>
-
-Package=<5>
-{{{
-}}}
-
-Package=<4>
-{{{
-}}}
-
-###############################################################################
-
Project: "RFFTW2st"="..\..\fftw-2.1.3\Win32\RFFTW2st\RFFTW2st.dsp" - Package Owner=<4>
Package=<5>
<pre>
<h1>Build Log</h1>
<h3>
---------------------Configuration: ctsim - Win32 Debug--------------------
+--------------------Configuration: libctsim - Win32 Debug--------------------
</h3>
<h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1C2.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP282.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
+"C:\ctsim\libctsim\projections.cpp"
+]
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP282.tmp"
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP283.tmp" with contents
[
-/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\views.cpp"
+/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
]
-Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1C2.tmp"
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1C3.tmp" with contents
+Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP283.tmp"
+<h3>Output Window</h3>
+Compiling...
+projections.cpp
+Creating library...
+<h3>
+--------------------Configuration: ctsim - Win32 Debug--------------------
+</h3>
+<h3>Command Lines</h3>
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP284.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\RSP1C3.tmp"
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP284.tmp"
<h3>Output Window</h3>
-Compiling...
-views.cpp
Linking...