+ int nDet = rProj.nDet();
+ int nView = rProj.nView();
+ ImageFile* pIF = new ImageFile (nDet, nView);
+ ImageFileArray v = pIF->getArray();
+ for (int iv = 0; iv < nView; iv++) {
+ DetectorValue* detval = rProj.getDetectorArray(iv).detValues();
+
+ for (int id = 0; id < nDet; id++)
+ v[id][iv] = detval[id];
+ }
+
+ ImageFileDocument* pRectDoc = theApp->newImageDoc ();
+ if (! pRectDoc) {
+ sys_error (ERR_SEVERE, "Unable to create image file");
+ return;
+ }
+ pRectDoc->setImageFile (pIF);
+ pIF->labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime());
+ std::ostringstream os;
+ os << "Convert projection file " << getFrame()->GetTitle().c_str() << " to rectangular image";
+ *theApp->getLog() << wxConvUTF8.cMB2WX(os.str().c_str()) << _T("\n");
+ pIF->labelAdd (os.str().c_str());
+ if (theApp->getAskDeleteNewDocs())
+ pRectDoc->Modify (true);
+ pRectDoc->UpdateAllViews();
+ pRectDoc->getView()->setInitialClientSize();
+ pRectDoc->Activate();
+}
+
+void
+ProjectionFileView::OnConvertPolar (wxCommandEvent& event)
+{
+ Projections& rProj = GetDocument()->getProjections();
+ DialogGetConvertPolarParameters dialogPolar (getFrameForChild(), _T("Convert Polar"), m_iDefaultPolarNX, m_iDefaultPolarNY,
+ m_iDefaultPolarInterpolation, -1, IDH_DLG_POLAR);
+ if (dialogPolar.ShowModal() == wxID_OK) {
+ wxProgressDialog dlgProgress (_T("Convert Polar"), _T("Conversion Progress"), 1, getFrameForChild(), wxPD_APP_MODAL);
+ wxString strInterpolation (dialogPolar.getInterpolationName(), wxConvUTF8);
+ m_iDefaultPolarNX = dialogPolar.getXSize();
+ m_iDefaultPolarNY = dialogPolar.getYSize();
+ ImageFile* pIF = new ImageFile (m_iDefaultPolarNX, m_iDefaultPolarNY);
+ m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.mb_str(wxConvUTF8));
+
+ if (! rProj.convertPolar (*pIF, m_iDefaultPolarInterpolation)) {
+ delete pIF;
+ *theApp->getLog() << _T("Error converting to Polar\n");
+ return;
+ }
+
+ ImageFileDocument* pPolarDoc = theApp->newImageDoc();
+ if (! pPolarDoc) {
+ sys_error (ERR_SEVERE, "Unable to create image file");
+ return;
+ }
+ pPolarDoc->setImageFile (pIF);
+ pIF->labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime());
+ std::ostringstream os;
+ os << "Convert projection file " << getFrame()->GetTitle().c_str() << " to polar image: xSize="
+ << m_iDefaultPolarNX << ", ySize=" << m_iDefaultPolarNY << ", interpolation="
+ << strInterpolation.c_str();
+ *theApp->getLog() << wxConvUTF8.cMB2WX(os.str().c_str()) << _T("\n");
+ pIF->labelAdd (os.str().c_str());
+ if (theApp->getAskDeleteNewDocs())
+ pPolarDoc->Modify (true);
+ pPolarDoc->UpdateAllViews ();
+ pPolarDoc->getView()->setInitialClientSize();
+ pPolarDoc->Activate();
+ }
+}
+
+void
+ProjectionFileView::OnConvertFFTPolar (wxCommandEvent& event)
+{
+ Projections& rProj = GetDocument()->getProjections();
+ DialogGetConvertPolarParameters dialogPolar (getFrameForChild(), _T("Convert to FFT Polar"), m_iDefaultPolarNX, m_iDefaultPolarNY,
+ m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad, IDH_DLG_FFT_POLAR);
+ if (dialogPolar.ShowModal() == wxID_OK) {
+ wxProgressDialog dlgProgress (_T("Convert FFT Polar"), _T("Conversion Progress"), 1, getFrameForChild(), wxPD_APP_MODAL);
+ wxString strInterpolation (dialogPolar.getInterpolationName(), wxConvUTF8);
+ m_iDefaultPolarNX = dialogPolar.getXSize();
+ m_iDefaultPolarNY = dialogPolar.getYSize();
+ m_iDefaultPolarZeropad = dialogPolar.getZeropad();
+ ImageFile* pIF = new ImageFile (m_iDefaultPolarNX, m_iDefaultPolarNY);
+
+ m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.mb_str(wxConvUTF8));
+ if (! rProj.convertFFTPolar (*pIF, m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad)) {
+ delete pIF;
+ *theApp->getLog() << _T("Error converting to polar\n");
+ return;
+ }
+ ImageFileDocument* pPolarDoc = theApp->newImageDoc();
+ if (! pPolarDoc) {
+ sys_error (ERR_SEVERE, "Unable to create image file");
+ return;
+ }
+ pPolarDoc->setImageFile (pIF);
+ pIF->labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime());
+ std::ostringstream os;
+ os << "Convert projection file " << getFrame()->GetTitle().c_str() << " to FFT polar image: xSize="
+ << m_iDefaultPolarNX << ", ySize=" << m_iDefaultPolarNY << ", interpolation="
+ << strInterpolation.c_str() << ", zeropad=" << m_iDefaultPolarZeropad;
+ *theApp->getLog() << wxConvUTF8.cMB2WX(os.str().c_str()) << _T("\n");
+ pIF->labelAdd (os.str().c_str());
+ if (theApp->getAskDeleteNewDocs())
+ pPolarDoc->Modify (true);
+ pPolarDoc->UpdateAllViews (this);
+ pPolarDoc->getView()->setInitialClientSize();
+ pPolarDoc->Activate();
+ }
+}
+
+void
+ProjectionFileView::OnPlotTThetaSampling (wxCommandEvent& event)
+{
+ 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
+ProjectionFileView::OnPlotHistogram (wxCommandEvent& event)
+{
+ Projections& rProj = GetDocument()->getProjections();
+ int nDet = rProj.nDet();
+ int nView = rProj.nView();
+
+ if (nDet < 1 || nView < 1)
+ return;
+
+ PlotFileDocument* pPlotDoc = theApp->newPlotDoc();
+ if (! pPlotDoc) {
+ sys_error (ERR_SEVERE, "Internal error: unable to create Plot file");
+ return;
+ }
+
+ DetectorValue* pdDetval = rProj.getDetectorArray(0).detValues();
+ double dMin = pdDetval[0], dMax = pdDetval[0];
+
+ for (int iv = 0; iv < nView; iv++) {
+ pdDetval = rProj.getDetectorArray(iv).detValues();
+ for (int id = 0; id < nDet; id++) {
+ double dV = pdDetval[id];
+ if (dV < dMin)
+ dMin = dV;
+ else if (dV > dMax)
+ dMax = dV;
+ }
+ }
+
+ double* pX = new double [NUMBER_HISTOGRAM_BINS];
+ double* pY = new double [NUMBER_HISTOGRAM_BINS];
+ double dBinWidth = (dMax - dMin) / NUMBER_HISTOGRAM_BINS;
+
+ for (int i = 0; i < NUMBER_HISTOGRAM_BINS; i++) {
+ pX[i] = dMin + (i + 0.5) * dBinWidth;
+ pY[i] = 0;
+ }
+ for (int j = 0; j < nView; j++) {
+ pdDetval = rProj.getDetectorArray(j).detValues();
+ for (int id = 0; id < nDet; id++) {
+ int iBin = nearest<int> ((pdDetval[id] - dMin) / dBinWidth);
+ if (iBin >= 0 && iBin < NUMBER_HISTOGRAM_BINS)
+ pY[iBin] += 1;
+ }
+ }
+ PlotFile& rPlotFile = pPlotDoc->getPlotFile();
+ std::ostringstream os;
+ os << "Histogram";
+ std::string title("title ");
+ title += os.str();
+ rPlotFile.addEzsetCommand (title.c_str());
+ rPlotFile.addEzsetCommand ("xlabel Detector Value");
+ rPlotFile.addEzsetCommand ("ylabel Count");
+ rPlotFile.addEzsetCommand ("box");
+ rPlotFile.addEzsetCommand ("grid");
+ rPlotFile.setCurveSize (2, NUMBER_HISTOGRAM_BINS);
+ rPlotFile.addColumn (0, pX);
+ rPlotFile.addColumn (1, pY);
+ rPlotFile.addDescription (rProj.remark());
+ os << " plot of " << dynamic_cast<wxFrame*>(GetDocument()->GetFirstView()->GetFrame())->GetTitle().c_str();
+ *theApp->getLog() << wxConvUTF8.cMB2WX(os.str().c_str()) << _T("\n");
+ rPlotFile.addDescription (os.str().c_str());
+ delete pX;
+ delete pY;
+ if (theApp->getAskDeleteNewDocs())
+ pPlotDoc->Modify (true);
+ pPlotDoc->getView()->getFrame()->Show(true);
+ pPlotDoc->UpdateAllViews ();
+ pPlotDoc->Activate();
+}
+
+
+void
+ProjectionFileView::OnConvertParallel (wxCommandEvent& event)
+{
+ Projections& rProj = GetDocument()->getProjections();
+ if (rProj.geometry() == Scanner::GEOMETRY_PARALLEL) {
+ wxMessageBox (_T("Projections are already parallel"), _T("Error"));
+ return;
+ }
+ wxProgressDialog dlgProgress (_T("Convert to Parallel"), _T("Conversion Progress"), 1, getFrameForChild(), wxPD_APP_MODAL);
+ Projections* pProjNew = rProj.interpolateToParallel();
+ ProjectionFileDocument* pProjDocNew = theApp->newProjectionDoc();
+ pProjDocNew->setProjections (pProjNew);
+
+ if (ProjectionFileView* projView = pProjDocNew->getView()) {
+ projView->OnUpdate (projView, NULL);
+ if (projView->getCanvas())
+ projView->getCanvas()->SetClientSize (pProjNew->nDet(), pProjNew->nView());
+ if (wxFrame* pFrame = projView->getFrame()) {
+ pFrame->Show(true);
+ pFrame->SetFocus();
+ pFrame->Raise();
+ }
+ GetDocumentManager()->ActivateView (projView, true);
+ }
+ if (theApp->getAskDeleteNewDocs())
+ pProjDocNew-> Modify(true);
+ pProjDocNew->UpdateAllViews (this);
+ pProjDocNew->getView()->setInitialClientSize();
+ pProjDocNew->Activate();
+}
+
+void
+ProjectionFileView::OnReconstructFourier (wxCommandEvent& event)
+{
+ Projections& rProj = GetDocument()->getProjections();
+ DialogGetConvertPolarParameters dialogPolar (getFrameForChild(), _T("Fourier Reconstruction"), m_iDefaultPolarNX, m_iDefaultPolarNY,
+ m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad, IDH_DLG_RECON_FOURIER);
+ if (dialogPolar.ShowModal() == wxID_OK) {
+ wxProgressDialog dlgProgress (_T("Reconstruction Fourier"), _T("Reconstruction Progress"), 1, getFrameForChild(), wxPD_APP_MODAL);
+ wxString strInterpolation (dialogPolar.getInterpolationName(), wxConvUTF8);
+ m_iDefaultPolarNX = dialogPolar.getXSize();
+ m_iDefaultPolarNY = dialogPolar.getYSize();
+ m_iDefaultPolarZeropad = dialogPolar.getZeropad();
+ ImageFile* pIF = new ImageFile (m_iDefaultPolarNX, m_iDefaultPolarNY);
+
+ m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.mb_str(wxConvUTF8));
+ if (! rProj.convertFFTPolar (*pIF, m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad)) {
+ delete pIF;
+ *theApp->getLog() << _T("Error converting to polar\n");
+ return;
+ }
+#ifdef HAVE_FFT
+ pIF->ifft(*pIF);
+#endif
+ pIF->magnitude(*pIF);
+ Fourier::shuffleFourierToNaturalOrder (*pIF);
+
+ ImageFileDocument* pPolarDoc = theApp->newImageDoc();
+ if (! pPolarDoc) {
+ sys_error (ERR_SEVERE, "Unable to create image file");
+ return;
+ }
+ pPolarDoc->setImageFile (pIF);
+ pIF->labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime());
+ std::ostringstream os;
+ os << "Reconstruct Fourier " << getFrame()->GetTitle().mb_str(wxConvUTF8) << ": xSize="
+ << m_iDefaultPolarNX << ", ySize=" << m_iDefaultPolarNY << ", interpolation="
+ << strInterpolation.mb_str(wxConvUTF8) << ", zeropad=" << m_iDefaultPolarZeropad;
+ *theApp->getLog() << wxConvUTF8.cMB2WX(os.str().c_str()) << _T("\n");
+ pIF->labelAdd (os.str().c_str());
+ if (theApp->getAskDeleteNewDocs())
+ pPolarDoc->Modify (true);
+ pPolarDoc->UpdateAllViews ();
+ pPolarDoc->getView()->setInitialClientSize();
+ pPolarDoc->Activate();
+ }
+}
+
+void
+ProjectionFileView::OnReconstructFBPRebin (wxCommandEvent& event)
+{
+ Projections& rProj = GetDocument()->getProjections();
+ doReconstructFBP (rProj, true);
+}
+
+void
+ProjectionFileView::OnReconstructFBP (wxCommandEvent& event)
+{
+ 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();
+ 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_iDefaultTrace, &defaultROI);
+
+ int retVal = dialogReconstruction.ShowModal();
+ if (retVal != wxID_OK)
+ return;
+
+ m_iDefaultNX = dialogReconstruction.getXSize();
+ m_iDefaultNY = dialogReconstruction.getYSize();
+ wxString optFilterName (dialogReconstruction.getFilterName(), wxConvUTF8);
+ m_iDefaultFilter = SignalFilter::convertFilterNameToID (optFilterName.mb_str(wxConvUTF8));
+ m_dDefaultFilterParam = dialogReconstruction.getFilterParam();
+ wxString optFilterMethodName (dialogReconstruction.getFilterMethodName(), wxConvUTF8);
+ m_iDefaultFilterMethod = ProcessSignal::convertFilterMethodNameToID(optFilterMethodName.mb_str(wxConvUTF8));
+ m_iDefaultZeropad = dialogReconstruction.getZeropad();
+ wxString optFilterGenerationName (dialogReconstruction.getFilterGenerationName(), wxConvUTF8);
+ m_iDefaultFilterGeneration = ProcessSignal::convertFilterGenerationNameToID (optFilterGenerationName.mb_str(wxConvUTF8));
+ wxString optInterpName (dialogReconstruction.getInterpName(), wxConvUTF8);
+ m_iDefaultInterpolation = Backprojector::convertInterpNameToID (optInterpName.mb_str(wxConvUTF8));
+ m_iDefaultInterpParam = dialogReconstruction.getInterpParam();
+ wxString optBackprojectName (dialogReconstruction.getBackprojectName(), wxConvUTF8);
+ m_iDefaultBackprojector = Backprojector::convertBackprojectNameToID (optBackprojectName.mb_str(wxConvUTF8));
+ m_iDefaultTrace = dialogReconstruction.getTrace();
+ dialogReconstruction.getROI (&defaultROI);
+
+ if (m_iDefaultNX <= 0 && m_iDefaultNY <= 0)
+ return;
+
+ std::ostringstream os;
+ os << "Reconstruct " << rProj.getFilename() << ": xSize=" << m_iDefaultNX << ", ySize=" << m_iDefaultNY << ", Filter=" << optFilterName.mb_str(wxConvUTF8) << ", FilterParam=" << m_dDefaultFilterParam << ", FilterMethod=" << optFilterMethodName.mb_str(wxConvUTF8) << ", FilterGeneration=" << optFilterGenerationName.mb_str(wxConvUTF8) << ", Zeropad=" << m_iDefaultZeropad << ", Interpolation=" << optInterpName.mb_str(wxConvUTF8) << ", InterpolationParam=" << m_iDefaultInterpParam << ", Backprojection=" << optBackprojectName.mb_str(wxConvUTF8);
+ 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.mb_str(wxConvUTF8),
+ m_dDefaultFilterParam, optFilterMethodName.mb_str(wxConvUTF8), m_iDefaultZeropad, optFilterGenerationName.mb_str(wxConvUTF8),
+ optInterpName.mb_str(wxConvUTF8), m_iDefaultInterpParam, optBackprojectName.mb_str(wxConvUTF8), m_iDefaultTrace,
+ &defaultROI, bRebinToParallel);
+
+ ReconstructDialog* pDlgReconstruct = new ReconstructDialog (*pReconstructor, rProj, *pImageFile, m_iDefaultTrace, getFrameForChild());
+ for (int iView = 0; iView < rProj.nView(); iView++) {
+ ::wxYield();
+ if (pDlgReconstruct->isCancelled() || ! pDlgReconstruct->reconstructView (iView, true)) {
+ delete pDlgReconstruct;
+ delete pReconstructor;
+ return;