From: Kevin M. Rosenberg Date: Tue, 13 Mar 2001 08:24:41 +0000 (+0000) Subject: r636: Optimized Rebinning, Added Reconstruct with Rebinning option X-Git-Tag: debian-4.5.3-3~381 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=4433641931aa27fd6a2b5ecd0102e6c5bbbccc46 r636: Optimized Rebinning, Added Reconstruct with Rebinning option --- diff --git a/ChangeLog b/ChangeLog index d723fa9..a9677d8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -21,6 +21,9 @@ * ctsim: Added imagefile export to DICOM files. * ctsim: Added clipboard cut/copy/paste for image files. + + * ctsim: Added Reconstruction with Rebinning for faster + divergent beam reconstructions. * ctsim: Added plot t-theta sampling to projection file menu. diff --git a/include/projections.h b/include/projections.h index b86f448..c88fca7 100644 --- a/include/projections.h +++ b/include/projections.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: projections.h,v 1.32 2001/03/13 04:44:25 kevin Exp $ +** $Id: projections.h,v 1.33 2001/03/13 08:24:41 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -65,7 +65,7 @@ public: THETA_RANGE_FOLD_TO_PI, }; - ParallelRaysums (Projections* pProjections, int iThetaRange); + ParallelRaysums (const Projections* pProjections, int iThetaRange); ~ParallelRaysums (); typedef std::vector CoordinateContainer; @@ -84,6 +84,7 @@ private: CoordinateContainer m_vecpCoordinates; CoordinateContainer m_vecpSortedByT; CoordinateContainer m_vecpSortedByTheta; + ParallelRaysumCoordinate* m_pCoordinates; int m_iNumCoordinates; int m_iNumView; int m_iNumDet; @@ -128,7 +129,7 @@ class Projections bool detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int view_num); bool detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int view_num); - Projections* interpolateToParallel(); + Projections* interpolateToParallel() const; bool convertPolar (ImageFile& rIF, int iInterpolation); bool convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad); diff --git a/include/reconstruct.h b/include/reconstruct.h index 61df525..473fa66 100644 --- a/include/reconstruct.h +++ b/include/reconstruct.h @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: reconstruct.h,v 1.8 2001/03/11 15:27:30 kevin Exp $ +** $Id: reconstruct.h,v 1.9 2001/03/13 08:24:41 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 @@ -45,13 +45,14 @@ struct ReconstructionROI { class Reconstructor { private: - const Projections& m_rProj; + const Projections& m_rOriginalProj; + const Projections* m_pProj; ImageFile& m_rImagefile; ProcessSignal* m_pProcessSignal; Backprojector* m_pBackprojector; int m_nFilteredProjections; int m_iTrace; - + const bool m_bRebinToParallel; bool m_bFail; std::string m_strFailMessage; @@ -61,7 +62,7 @@ class Reconstructor Reconstructor (const Projections& rProj, ImageFile& rIF, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace, - ReconstructionROI* pROI = NULL, SGP* pSGP = NULL); + ReconstructionROI* pROI = NULL, bool bRebinToParallel = false, SGP* pSGP = NULL); ~Reconstructor (); diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 60f9467..3eab99c 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: projections.cpp,v 1.62 2001/03/13 04:44:25 kevin Exp $ +** $Id: projections.cpp,v 1.63 2001/03/13 08:24:41 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 @@ -895,59 +895,53 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa m_dFocalLength = 510; m_dSourceDetectorLength = 890; m_detInc = convertDegreesToRadians (3.06976 / 60); + m_dFanBeamAngle = (iNDets + 1) * m_detInc; m_detStart = -(m_dFanBeamAngle / 2); m_rotInc = TWOPI / static_cast(iNViews); m_rotStart = HALFPI; - m_dFanBeamAngle = (iNDets + 1) * m_detInc; m_dViewDiameter = sin (m_dFanBeamAngle / 2) * m_dFocalLength * 2; 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 dCenter = (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); + pdCosScale[i] = 1. / cos ((i - dCenter) * m_detInc); long lDataPos = 0; for (int iv = 0; iv < iNViews; iv++) { unsigned char* pArgBase = pData + lDataPos; - unsigned char* p = pArgBase+0; - SwapBytes4IfLittleEndian (p); + unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p); long lProjNumber = *reinterpret_cast(p); - p = pArgBase+20; - SwapBytes4IfLittleEndian (p); + p = pArgBase+20; SwapBytes4IfLittleEndian (p); long lEscale = *reinterpret_cast(p); - p = pArgBase+28; - SwapBytes4IfLittleEndian (p); + p = pArgBase+28; SwapBytes4IfLittleEndian (p); long lTime = *reinterpret_cast(p); - p = pArgBase + 4; - SwapBytes4IfLittleEndian (p); + p = pArgBase + 4; SwapBytes4IfLittleEndian (p); double dAlpha = *reinterpret_cast(p) + HALFPI; - p = pArgBase+12; - SwapBytes4IfLittleEndian (p); + p = pArgBase+12; SwapBytes4IfLittleEndian (p); double dAlign = *reinterpret_cast(p); - p = pArgBase + 16; - SwapBytes4IfLittleEndian (p); + p = pArgBase + 16; SwapBytes4IfLittleEndian (p); double dMaxValue = *reinterpret_cast(p); DetectorArray& detArray = getDetectorArray (iv); detArray.setViewAngle (dAlpha); DetectorValue* detval = detArray.detValues(); - double dViewScale = 2294.4871 * ::pow (2.0, -lEscale); + double dViewScale = 1. / (2294.4871 * ::pow (2.0, -lEscale)); lDataPos += 32; for (int id = 0; id < iNDets; id++) { - int iV = pData[lDataPos+1] + 256 * pData[lDataPos]; + int iV = pData[lDataPos+1] + (pData[lDataPos] << 8); if (iV > 32767) // two's complement signed conversion iV = iV - 65536; - detval[id] = iV / (dViewScale * pdCosScale[id]); + detval[id] = iV * dViewScale * pdCosScale[id]; lDataPos += 2; } } @@ -957,10 +951,10 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa } Projections* -Projections::interpolateToParallel () +Projections::interpolateToParallel () const { if (m_geometry == Scanner::GEOMETRY_PARALLEL) - return this; + return const_cast(this); int nDet = m_nDet; int nView = m_nView; @@ -1045,9 +1039,9 @@ Projections::interpolateToParallel () // /////////////////////////////////////////////////////////////////////////////// -ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange) +ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange) : m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()), - m_iThetaRange (iThetaRange) + m_iThetaRange (iThetaRange), m_pCoordinates(NULL) { int iGeometry = pProjections->geometry(); double dDetInc = pProjections->detInc(); @@ -1055,9 +1049,10 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange) double dFocalLength = pProjections->focalLength(); m_iNumCoordinates = m_iNumView * m_iNumDet; + m_pCoordinates = new ParallelRaysumCoordinate [m_iNumCoordinates]; m_vecpCoordinates.reserve (m_iNumCoordinates); for (int i = 0; i < m_iNumCoordinates; i++) - m_vecpCoordinates[i] = new ParallelRaysumCoordinate; + m_vecpCoordinates[i] = m_pCoordinates + i; int iCoordinate = 0; for (int iV = 0; iV < m_iNumView; iV++) { @@ -1096,8 +1091,7 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange) ParallelRaysums::~ParallelRaysums() { - for (int i = 0; i < m_iNumCoordinates; i++) - delete m_vecpCoordinates[i]; + delete m_pCoordinates; } ParallelRaysums::CoordinateContainer& diff --git a/libctsim/reconstruct.cpp b/libctsim/reconstruct.cpp index c3b6e8b..8caed20 100644 --- a/libctsim/reconstruct.cpp +++ b/libctsim/reconstruct.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: reconstruct.cpp,v 1.16 2001/03/11 15:27:30 kevin Exp $ +** $Id: reconstruct.cpp,v 1.17 2001/03/13 08:24:41 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 @@ -53,11 +53,13 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch double filt_param, const char* const filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int iTrace, - ReconstructionROI* pROI, SGP* pSGP) - : m_rProj(rProj), m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), m_iTrace(iTrace), - m_bFail(false), m_adPlotXAxis(0) + ReconstructionROI* pROI, bool bRebinToParallel, SGP* pSGP) + : m_rOriginalProj(rProj), + m_pProj(bRebinToParallel ? m_rOriginalProj.interpolateToParallel() : &m_rOriginalProj), + m_rImagefile(rIF), m_pProcessSignal(0), m_pBackprojector(0), + m_iTrace(iTrace), m_bRebinToParallel(bRebinToParallel), m_bFail(false), m_adPlotXAxis(0) { - m_nFilteredProjections = m_rProj.nDet() * interpFactor; + m_nFilteredProjections = m_pProj->nDet() * interpFactor; #ifdef HAVE_BSPLINE_INTERP int spline_order = 0, zoom_factor = 0; @@ -69,10 +71,10 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch } #endif - double filterBW = 1. / m_rProj.detInc(); - m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_rProj.detInc(), - m_rProj.nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace, - m_rProj.geometry(), m_rProj.focalLength(), m_rProj.sourceDetectorLength(), pSGP); + double filterBW = 1. / m_pProj->detInc(); + m_pProcessSignal = new ProcessSignal (filterName, filterMethodName, filterBW, m_pProj->detInc(), + m_pProj->nDet(), filt_param, "spatial", filterGenerationName, zeropad, interpFactor, iTrace, + m_pProj->geometry(), m_pProj->focalLength(), m_pProj->sourceDetectorLength(), pSGP); if (m_pProcessSignal->fail()) { m_bFail = true; @@ -82,7 +84,7 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch return; } - m_pBackprojector = new Backprojector (m_rProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI); + m_pBackprojector = new Backprojector (*m_pProj, m_rImagefile, backprojectName, interpName, interpFactor, pROI); if (m_pBackprojector->fail()) { m_bFail = true; m_strFailMessage = "Error creating backprojector: "; @@ -93,17 +95,20 @@ Reconstructor::Reconstructor (const Projections& rProj, ImageFile& rIF, const ch } #ifdef HAVE_SGP - m_adPlotXAxis = new double [m_rProj.nDet()]; - double x = - ((m_rProj.nDet() - 1) / 2) * m_rProj.detInc(); - double xInc = m_rProj.detInc(); + m_adPlotXAxis = new double [m_pProj->nDet()]; + double x = - ((m_pProj->nDet() - 1) / 2) * m_pProj->detInc(); + double xInc = m_pProj->detInc(); - for (int i = 0; i < m_rProj.nDet(); i++, x += xInc) + for (int i = 0; i < m_pProj->nDet(); i++, x += xInc) m_adPlotXAxis[i] = x; #endif } Reconstructor::~Reconstructor () { + if (m_bRebinToParallel) + delete m_pProj; + delete m_pBackprojector; delete m_pProcessSignal; delete m_adPlotXAxis; @@ -139,7 +144,7 @@ Reconstructor::plotFilter (SGP* pSGP) void Reconstructor::reconstructAllViews () { - reconstructView (0, m_rProj.nView()); + reconstructView (0, m_pProj->nView()); postProcessing(); } @@ -156,24 +161,24 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool double* adFilteredProj = new double [m_nFilteredProjections]; // filtered projections if (iViewCount <= 0) - iViewCount = m_rProj.nView() - iStartView; + iViewCount = m_pProj->nView() - iStartView; for (int iView = iStartView; iView < (iStartView + iViewCount); iView++) { if (m_iTrace == Trace::TRACE_CONSOLE) - std::cout <<"Reconstructing view " << iView << " (last = " << m_rProj.nView() - 1 << ")\n"; + std::cout <<"Reconstructing view " << iView << " (last = " << m_pProj->nView() - 1 << ")\n"; - const DetectorArray& rDetArray = m_rProj.getDetectorArray (iView); + const DetectorArray& rDetArray = m_pProj->getDetectorArray (iView); const DetectorValue* detval = rDetArray.detValues(); m_pProcessSignal->filterSignal (detval, adFilteredProj); #ifdef HAVE_BSPLINE_INTERP if (interp_type == I_BSPLINE) - bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj); + bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj); #ifdef HAVE_SGP if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE && pSGP) { - bspline (m_rProj.nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj); + bspline (m_pProj->nDet(), zoom_factor, spline_order, adFilteredProj, adFilteredProj); ezplot_1d (adFilteredProj, m_nFilteredProjections); } #endif @@ -201,14 +206,14 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool ezplotProj.ezset ("box."); ezplotProj.ezset ("grid."); #if 0 // workaround c++ optimizer bug, now disabled by using /O1 in code - double* pdDetval = new double [m_rProj.nDet()]; - for (unsigned int id = 0; id < m_rProj.nDet(); id++) { + double* pdDetval = new double [m_pProj->nDet()]; + for (unsigned int id = 0; id < m_pProj->nDet(); id++) { pdDetval[id] = detval[id]; } - ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_rProj.nDet()); + ezplotProj.addCurve (m_adPlotXAxis, pdDetval, m_pProj->nDet()); delete pdDetval; #else - ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet()); + ezplotProj.addCurve (m_adPlotXAxis, detval, m_pProj->nDet()); #endif pSGP->setTextPointSize (12); ezplotProj.plot (pSGP); diff --git a/msvc/ctsim/ctsim.plg b/msvc/ctsim/ctsim.plg index ba431eb..5c51373 100644 --- a/msvc/ctsim/ctsim.plg +++ b/msvc/ctsim/ctsim.plg @@ -6,19 +6,13 @@ --------------------Configuration: ctsim - Win32 Debug--------------------

Command Lines

-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD2.tmp" with contents +Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP15B.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\docs.cpp" +"C:\ctsim\src\views.cpp" ] -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 /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 "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 +Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP15B.tmp" +Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP15C.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 @@ -45,12 +39,10 @@ winmm.lib rpcrt4.lib ws2_32.lib ../libctsim/Debug/libctsim.lib libcmtd.lib ..\.. \wx2.2.5\lib\zlibd.lib \wx2.2.5\lib\tiffd.lib ] -Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSPD4.tmp" +Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP15C.tmp"

Output Window

Compiling... -docs.cpp -Compiling... -graph3dview.cpp +views.cpp Linking... diff --git a/src/ctsim.h b/src/ctsim.h index ef100ad..b7d20a9 100644 --- a/src/ctsim.h +++ b/src/ctsim.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: ctsim.h,v 1.58 2001/03/11 17:55:29 kevin Exp $ +** $Id: ctsim.h,v 1.59 2001/03/13 08:24:41 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 @@ -306,6 +306,7 @@ enum { PJMENU_FILE_PROPERTIES, PJMENU_RECONSTRUCT_FBP, + PJMENU_RECONSTRUCT_FBP_REBIN, PJMENU_RECONSTRUCT_FOURIER, PJMENU_CONVERT_POLAR, PJMENU_CONVERT_FFT_POLAR, diff --git a/src/threadrecon.cpp b/src/threadrecon.cpp index e8deca4..2fde627 100644 --- a/src/threadrecon.cpp +++ b/src/threadrecon.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2001 Kevin Rosenberg ** -** $Id: threadrecon.cpp,v 1.24 2001/03/11 15:27:30 kevin Exp $ +** $Id: threadrecon.cpp,v 1.25 2001/03/13 08:24:41 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 @@ -54,18 +54,23 @@ ReconstructorSupervisorThread::ReconstructorSupervisorThread (ProjectionFileView* pProjView, int iNX, int iNY, const char* pszFilterName, double dFilterParam, const char* pszFilterMethod, int iZeropad, const char* pszFilterGenerationName, const char* pszInterpName, int iInterpParam, - const char* pszBackprojectName, const char* const pszLabel, ReconstructionROI* pROI) + const char* pszBackprojectName, const char* const pszLabel, ReconstructionROI* pROI, bool bRebinToParallel) : SupervisorThread(), m_pProjView(pProjView), m_iNX(iNX), m_iNY(iNY), m_strFilterName(pszFilterName), m_dFilterParam(dFilterParam), m_strFilterMethod(pszFilterMethod), m_iZeropad(iZeropad), m_strFilterGenerationName(pszFilterGenerationName), m_strInterpName(pszInterpName), m_iInterpParam(iInterpParam), m_strBackprojectName(pszBackprojectName), - m_strLabel(pszLabel), m_reconROI(*pROI) + m_strLabel(pszLabel), m_reconROI(*pROI), m_bRebinToParallel(bRebinToParallel) { } wxThread::ExitCode ReconstructorSupervisorThread::Entry() { - ReconstructorSupervisor reconSupervisor (this, m_pProjView, m_iNX, m_iNY, + Projections* pProj = &m_pProjView->GetDocument()->getProjections(); + + if (m_bRebinToParallel) + pProj = pProj->interpolateToParallel(); + + ReconstructorSupervisor reconSupervisor (this, pProj, m_pProjView, m_iNX, m_iNY, m_strFilterName.c_str(), m_dFilterParam, m_strFilterMethod.c_str(), m_iZeropad, m_strFilterGenerationName.c_str(), m_strInterpName.c_str(), m_iInterpParam, m_strBackprojectName.c_str(), m_strLabel.c_str(), &m_reconROI); @@ -86,6 +91,9 @@ ReconstructorSupervisorThread::Entry() reconSupervisor.onDone(); reconSupervisor.deleteWorkers(); + if (m_bRebinToParallel) + delete pProj; + return static_cast(0); } @@ -101,12 +109,13 @@ ReconstructorSupervisorThread::OnExit() // ///////////////////////////////////////////////////////////////////// -ReconstructorSupervisor::ReconstructorSupervisor (SupervisorThread* pThread, ProjectionFileView* pProjView, - int iImageNX, int iImageNY, const char* pszFilterName, double dFilterParam, const char* pszFilterMethod, - int iZeropad, const char* pszFilterGenerationName, const char* pszInterpName, int iInterpParam, - const char* pszBackprojectName, const char* const pszLabel, ReconstructionROI* pROI) +ReconstructorSupervisor::ReconstructorSupervisor (SupervisorThread* pThread, Projections* pProj, + ProjectionFileView* pProjView, int iImageNX, int iImageNY, const char* pszFilterName, double dFilterParam, + const char* pszFilterMethod, int iZeropad, const char* pszFilterGenerationName, + const char* pszInterpName, int iInterpParam, const char* pszBackprojectName, const char* const pszLabel, + ReconstructionROI* pROI) : BackgroundSupervisor (pThread, pProjView->GetFrame(), pProjView->GetDocument(), "Reconstructing", pProjView->GetDocument()->getProjections().nView()), - m_pProjView(pProjView), m_pProjDoc(pProjView->GetDocument()), + m_pProj(pProj), m_pProjView(pProjView), m_pProjDoc(pProjView->GetDocument()), m_iImageNX(iImageNX), m_iImageNY(iImageNY), m_pszFilterName(pszFilterName), m_dFilterParam(dFilterParam), m_pszFilterMethod(pszFilterMethod), m_iZeropad(iZeropad), m_pszFilterGenerationName(pszFilterGenerationName), m_pszInterpName(pszInterpName), @@ -132,8 +141,8 @@ BackgroundWorkerThread* ReconstructorSupervisor::createWorker (int iThread, int iStartUnit, int iNumUnits) { ReconstructorWorker* pThread = new ReconstructorWorker (this, iThread, iStartUnit, iNumUnits); - pThread->SetParameters (m_pProjView, m_vecpChildImageFile[iThread], m_pszFilterName, m_dFilterParam, - m_pszFilterMethod, m_iZeropad, m_pszFilterGenerationName, m_pszInterpName, + pThread->SetParameters (m_pProj, m_pProjView, m_vecpChildImageFile[iThread], m_pszFilterName, + m_dFilterParam, m_pszFilterMethod, m_iZeropad, m_pszFilterGenerationName, m_pszInterpName, m_iInterpParam, m_pszBackprojectName, m_pReconROI); return pThread; @@ -146,7 +155,7 @@ ReconstructorSupervisor::onDone() wxCriticalSectionLocker critsect (doneSection); ImageFile* pImageFile = getImageFile(); - pImageFile->labelAdd (m_pProjView->GetDocument()->getProjections().getLabel()); + pImageFile->labelAdd (m_pProj->getLabel()); pImageFile->labelAdd (m_pszLabel, getTimerEnd()); wxCommandEvent eventLog (wxEVT_COMMAND_MENU_SELECTED, MAINMENU_LOG_EVENT ); @@ -189,11 +198,12 @@ ReconstructorSupervisor::getImageFile() ///////////////////////////////////////////////////////////////////// void -ReconstructorWorker::SetParameters (ProjectionFileView* pProjView, ImageFile* pImageFile, +ReconstructorWorker::SetParameters (const Projections* pProj, ProjectionFileView* pProjView, ImageFile* pImageFile, const char* pszFilterName, double dFilterParam, const char* pszFilterMethod, int iZeropad, const char* pszFilterGenerationName, const char* pszInterpName, int iInterpParam, const char* pszBackprojectName, ReconstructionROI* pROI) { + m_pProj = pProj; m_pProjView = pProjView; m_pImageFile = pImageFile; m_pszFilterName = pszFilterName; @@ -210,10 +220,9 @@ ReconstructorWorker::SetParameters (ProjectionFileView* pProjView, ImageFile* pI wxThread::ExitCode ReconstructorWorker::Entry () { - Reconstructor* pReconstructor = new Reconstructor (m_pProjView->GetDocument()->getProjections(), - *m_pImageFile, m_pszFilterName, m_dFilterParam, m_pszFilterMethod, m_iZeropad, - m_pszFilterGenerationName, m_pszInterpName, m_iInterpParam, m_pszBackprojectName, Trace::TRACE_NONE, - m_pReconROI); + Reconstructor* pReconstructor = new Reconstructor (*m_pProj, *m_pImageFile, m_pszFilterName, + m_dFilterParam, m_pszFilterMethod, m_iZeropad, m_pszFilterGenerationName, m_pszInterpName, + m_iInterpParam, m_pszBackprojectName, Trace::TRACE_NONE, m_pReconROI, false); bool bFail = pReconstructor->fail(); std::string failMsg; diff --git a/src/threadrecon.h b/src/threadrecon.h index 0b56a8c..d0987ca 100644 --- a/src/threadrecon.h +++ b/src/threadrecon.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2001 Kevin Rosenberg ** -** $Id: threadrecon.h,v 1.13 2001/03/11 15:27:30 kevin Exp $ +** $Id: threadrecon.h,v 1.14 2001/03/13 08:24:41 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 @@ -61,12 +61,13 @@ private: const std::string m_strBackprojectName; const std::string m_strLabel; ReconstructionROI m_reconROI; + const bool m_bRebinToParallel; public: ReconstructorSupervisorThread(ProjectionFileView* pProjView, int iNX, int iNY, const char* pszFilterName, double dFilterParam, const char* pszFilterMethod, int iZeropad, const char* pszFilterGenerationName, const char* pszInterpName, int iInterpParam, const char* pszBackprojectName, const char* const pszLabel, - ReconstructionROI* pROI); + ReconstructionROI* pROI, bool bRebinToParallel); virtual wxThread::ExitCode Entry(); @@ -79,6 +80,7 @@ class ReconstructorSupervisor : public BackgroundSupervisor { private: std::vector m_vecpChildImageFile; + const Projections* m_pProj; ProjectionFileView* m_pProjView; ProjectionFileDocument* m_pProjDoc; @@ -97,10 +99,10 @@ private: ReconstructionROI* m_pReconROI; public: - ReconstructorSupervisor (SupervisorThread* pMyThread, ProjectionFileView* pProjView, int iNX, int iNY, const char* pszFilterName, - double dFilterParam, const char* pszFilterMethod, int iZeropad, const char* pszFilterGenerationName, - const char* pszInterpName, int iInterpParam, const char* pszBackprojectName, const char* const pszLabel, - ReconstructionROI* pReconROI); + ReconstructorSupervisor (SupervisorThread* pMyThread, Projections* pProj, ProjectionFileView* pProjView, + int iNX, int iNY, const char* pszFilterName, double dFilterParam, const char* pszFilterMethod, int iZeropad, + const char* pszFilterGenerationName, const char* pszInterpName, int iInterpParam, + const char* pszBackprojectName, const char* const pszLabel, ReconstructionROI* pReconROI); virtual BackgroundWorkerThread* createWorker (int iThread, int iStartUnit, int iNumUnits); @@ -113,10 +115,9 @@ public: }; - - class ReconstructorWorker : public BackgroundWorkerThread { private: + const Projections* m_pProj; ProjectionFileView* m_pProjView; ImageFile* m_pImageFile; const char* m_pszFilterName; @@ -134,7 +135,7 @@ public: : BackgroundWorkerThread (pSupervisor, iThread, iStartView, iNumViews) {} - void SetParameters (ProjectionFileView* pProjFile, ImageFile* pImageFile, + void SetParameters (const Projections* pProj, ProjectionFileView* pProjFile, ImageFile* pImageFile, const char* const pszFilterName, double dFilterParam, const char* const pszFilterMethod, int iZeropad, const char* const pszFilterGenerationName, const char* const pszInterpName, int iInterpParam, const char* pszBackprojectName, ReconstructionROI* pROI); diff --git a/src/views.cpp b/src/views.cpp index dfce190..8b40f35 100644 --- a/src/views.cpp +++ b/src/views.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: views.cpp,v 1.136 2001/03/13 04:44:25 kevin Exp $ +** $Id: views.cpp,v 1.137 2001/03/13 08:24:41 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 @@ -2378,6 +2378,7 @@ IMPLEMENT_DYNAMIC_CLASS(ProjectionFileView, wxView) BEGIN_EVENT_TABLE(ProjectionFileView, wxView) EVT_MENU(PJMENU_FILE_PROPERTIES, ProjectionFileView::OnProperties) EVT_MENU(PJMENU_RECONSTRUCT_FBP, ProjectionFileView::OnReconstructFBP) +EVT_MENU(PJMENU_RECONSTRUCT_FBP_REBIN, ProjectionFileView::OnReconstructFBPRebin) EVT_MENU(PJMENU_RECONSTRUCT_FOURIER, ProjectionFileView::OnReconstructFourier) EVT_MENU(PJMENU_CONVERT_POLAR, ProjectionFileView::OnConvertPolar) EVT_MENU(PJMENU_CONVERT_FFT_POLAR, ProjectionFileView::OnConvertFFTPolar) @@ -2599,10 +2600,23 @@ ProjectionFileView::OnReconstructFourier (wxCommandEvent& event) wxMessageBox ("Fourier Reconstruction is not yet supported", "Unimplemented function"); } +void +ProjectionFileView::OnReconstructFBPRebin (wxCommandEvent& event) +{ + Projections& rProj = GetDocument()->getProjections(); + doReconstructFBP (rProj, true); +} + void ProjectionFileView::OnReconstructFBP (wxCommandEvent& event) { - const Projections& rProj = GetDocument()->getProjections(); + Projections& rProj = GetDocument()->getProjections(); + doReconstructFBP (rProj, false); +} + +void +ProjectionFileView::doReconstructFBP (const Projections& rProj, bool bRebinToParallel) +{ ReconstructionROI defaultROI; defaultROI.m_dXMin = -rProj.phmLen() / 2; defaultROI.m_dXMax = defaultROI.m_dXMin + rProj.phmLen(); @@ -2641,14 +2655,17 @@ ProjectionFileView::OnReconstructFBP (wxCommandEvent& event) std::ostringstream os; os << "Reconstruct " << rProj.getFilename() << ": xSize=" << m_iDefaultNX << ", ySize=" << m_iDefaultNY << ", Filter=" << optFilterName.c_str() << ", FilterParam=" << m_dDefaultFilterParam << ", FilterMethod=" << optFilterMethodName.c_str() << ", FilterGeneration=" << optFilterGenerationName.c_str() << ", Zeropad=" << m_iDefaultZeropad << ", Interpolation=" << optInterpName.c_str() << ", InterpolationParam=" << m_iDefaultInterpParam << ", Backprojection=" << optBackprojectName.c_str(); - + if (bRebinToParallel) + os << "; Interpolate to Parallel"; + Timer timerRecon; ImageFile* pImageFile = NULL; if (m_iDefaultTrace > Trace::TRACE_CONSOLE) { 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(), - optInterpName.c_str(), m_iDefaultInterpParam, optBackprojectName.c_str(), m_iDefaultTrace, &defaultROI); + optInterpName.c_str(), m_iDefaultInterpParam, optBackprojectName.c_str(), m_iDefaultTrace, + &defaultROI, bRebinToParallel); ReconstructDialog* pDlgReconstruct = new ReconstructDialog (*pReconstructor, rProj, *pImageFile, m_iDefaultTrace, getFrameForChild()); for (int iView = 0; iView < rProj.nView(); iView++) { @@ -2671,10 +2688,10 @@ ProjectionFileView::OnReconstructFBP (wxCommandEvent& event) } else { #if HAVE_WXTHREADS if (theApp->getUseBackgroundTasks()) { - ReconstructorSupervisorThread* pReconstructor = new ReconstructorSupervisorThread (this, - m_iDefaultNX, m_iDefaultNY, optFilterName.c_str(), - m_dDefaultFilterParam, optFilterMethodName.c_str(), m_iDefaultZeropad, optFilterGenerationName.c_str(), - optInterpName.c_str(), m_iDefaultInterpParam, optBackprojectName.c_str(), os.str().c_str(), &defaultROI); + ReconstructorSupervisorThread* pReconstructor = new ReconstructorSupervisorThread (this, m_iDefaultNX, + m_iDefaultNY, optFilterName.c_str(), m_dDefaultFilterParam, optFilterMethodName.c_str(), + m_iDefaultZeropad, optFilterGenerationName.c_str(), optInterpName.c_str(), m_iDefaultInterpParam, + optBackprojectName.c_str(), os.str().c_str(), &defaultROI, bRebinToParallel); if (pReconstructor->Create() != wxTHREAD_NO_ERROR) { sys_error (ERR_SEVERE, "Error creating reconstructor thread"); delete pReconstructor; @@ -2687,11 +2704,12 @@ ProjectionFileView::OnReconstructFBP (wxCommandEvent& event) #endif { pImageFile = new ImageFile (m_iDefaultNX, m_iDefaultNY); + wxProgressDialog dlgProgress (wxString("Reconstruction"), wxString("Reconstruction Progress"), rProj.nView() + 1, getFrameForChild(), wxPD_CAN_ABORT ); Reconstructor* pReconstructor = new Reconstructor (rProj, *pImageFile, optFilterName.c_str(), m_dDefaultFilterParam, optFilterMethodName.c_str(), m_iDefaultZeropad, optFilterGenerationName.c_str(), - optInterpName.c_str(), m_iDefaultInterpParam, optBackprojectName.c_str(), m_iDefaultTrace, &defaultROI); + optInterpName.c_str(), m_iDefaultInterpParam, optBackprojectName.c_str(), m_iDefaultTrace, + &defaultROI, bRebinToParallel); - wxProgressDialog dlgProgress (wxString("Reconstruction"), wxString("Reconstruction Progress"), rProj.nView() + 1, getFrameForChild(), wxPD_CAN_ABORT ); for (int iView = 0; iView < rProj.nView(); iView++) { pReconstructor->reconstructView (iView, 1); if (! dlgProgress.Update (iView + 1)) { @@ -2800,6 +2818,7 @@ ProjectionFileView::CreateChildFrame(wxDocument *doc, wxView *view) 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_FBP_REBIN, "Filtered &Backprojection (Rebin to Parallel)...\tCtrl-B", "Reconstruct image using filtered backprojection"); reconstruct_menu->Append (PJMENU_RECONSTRUCT_FOURIER, "&Fourier...\tCtrl-E", "Reconstruct image using inverse Fourier"); reconstruct_menu->Enable (PJMENU_RECONSTRUCT_FOURIER, false); @@ -2821,14 +2840,15 @@ ProjectionFileView::CreateChildFrame(wxDocument *doc, wxView *view) subframe->SetMenuBar(menu_bar); subframe->Centre(wxBOTH); - wxAcceleratorEntry accelEntries[6]; + wxAcceleratorEntry accelEntries[7]; accelEntries[0].Set (wxACCEL_CTRL, static_cast('L'), PJMENU_CONVERT_POLAR); accelEntries[1].Set (wxACCEL_CTRL, static_cast('M'), PJMENU_CONVERT_FFT_POLAR); accelEntries[2].Set (wxACCEL_CTRL, static_cast('R'), PJMENU_RECONSTRUCT_FBP); - accelEntries[3].Set (wxACCEL_CTRL, static_cast('E'), PJMENU_RECONSTRUCT_FOURIER); - accelEntries[4].Set (wxACCEL_CTRL, static_cast('I'), PJMENU_FILE_PROPERTIES); - accelEntries[5].Set (wxACCEL_CTRL, static_cast('T'), PJMENU_PLOT_TTHETA_SAMPLING); - wxAcceleratorTable accelTable (6, accelEntries); + accelEntries[3].Set (wxACCEL_CTRL, static_cast('B'), PJMENU_RECONSTRUCT_FBP_REBIN); + accelEntries[4].Set (wxACCEL_CTRL, static_cast('E'), PJMENU_RECONSTRUCT_FOURIER); + accelEntries[5].Set (wxACCEL_CTRL, static_cast('I'), PJMENU_FILE_PROPERTIES); + accelEntries[6].Set (wxACCEL_CTRL, static_cast('T'), PJMENU_PLOT_TTHETA_SAMPLING); + wxAcceleratorTable accelTable (7, accelEntries); subframe->SetAcceleratorTable (accelTable); return subframe; diff --git a/src/views.h b/src/views.h index 9aca964..69c2a14 100644 --- a/src/views.h +++ b/src/views.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: views.h,v 1.50 2001/03/11 17:55:29 kevin Exp $ +** $Id: views.h,v 1.51 2001/03/13 08:24:41 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 @@ -241,6 +241,7 @@ public: bool OnClose (bool deleteWindow = true); void OnProperties (wxCommandEvent& event); void OnReconstructFBP (wxCommandEvent& event); + void OnReconstructFBPRebin (wxCommandEvent& event); void OnReconstructFourier (wxCommandEvent& event); void OnConvertPolar (wxCommandEvent& event); void OnConvertFFTPolar (wxCommandEvent& event); @@ -248,6 +249,8 @@ public: void OnConvertParallel (wxCommandEvent& event); void OnArtifactReduction (wxCommandEvent& event); + void doReconstructFBP (const Projections& rProj, bool bRebinToParallel); + #if CTSIM_MDI wxDocMDIChildFrame* getFrame() { return m_pFrame; } #else