Merge pull request #85 from PauloCarvalhoRJ/FFT2D

Fft2d

Author: PauloCarvalhoRJ

GammaRay.pro

@@ -142,7 +142,8 @@ SOURCES += main.cpp\
     viewer3d/view3dlistrecord.cpp \
     viewer3d/view3dviewdata.cpp \
     viewer3d/view3dconfigwidgets/v3dcfgwidforattributeinmapcartesiangrid.cpp \
-    domain/auxiliary/dataloader.cpp
+    domain/auxiliary/dataloader.cpp \
+    array3d.cpp
 
 HEADERS  += mainwindow.h \
     domain/project.h \
@@ -266,7 +267,8 @@ HEADERS  += mainwindow.h \
     viewer3d/view3dlistrecord.h \
     viewer3d/view3dviewdata.h \
     viewer3d/view3dconfigwidgets/v3dcfgwidforattributeinmapcartesiangrid.h \
-    domain/auxiliary/dataloader.h
+    domain/auxiliary/dataloader.h \
+    array3d.h
 
 FORMS    += mainwindow.ui \
     gslib/gslibparams/widgets/widgetgslibpardouble.ui \
@@ -362,12 +364,13 @@ LIBS        += -lvtkGUISupportQt$$VTK_VERSION_SUFFIX \
                -lvtkImagingSources$$VTK_VERSION_SUFFIX \
                -lvtkImagingCore$$VTK_VERSION_SUFFIX \
                -lvtkFiltersCore$$VTK_VERSION_SUFFIX \
-               -lvtkFiltersExtraction$$VTK_VERSION_SUFFIX
+               -lvtkFiltersExtraction$$VTK_VERSION_SUFFIX \
+               -lvtkImagingFourier$$VTK_VERSION_SUFFIX
 
 # The application version
 # Don't forget to update the Util::importSettingsFromPreviousVersion() method to
 # enable the import of registry/user settings of previous versions.
-VERSION = 2.1
+VERSION = 2.2
 
 # Define a preprocessor macro so we can get the application version in application code.
 DEFINES += APP_VERSION=\\\"$$VERSION\\\"

array3d.cpp

@@ -0,0 +1,3 @@
+#include "array3d.h"
+
+

array3d.h

@@ -0,0 +1,43 @@
+#ifndef ARRAY3D_H
+#define ARRAY3D_H
+
+#include <vector>
+#include <cstdlib>
+#include <stdexcept>
+
+/**
+ * The Array3D class is a generic array of any object type with rank 3 (3 dimensions, IJK).
+ * It wraps an unidimensional array, hiding index arithmetic, allowing a more readable code (e.g. a(i,j,k))
+ * while preserving memory locality, leading to better performance (e.g. more cache hits), since multidimensional
+ * arrays tend to scatter in memory.
+ * The Array3D can be used as a two-dimensional array in code (e.g. a(i,j) assumes k = 0)
+ */
+template <class TYPE>
+class Array3D
+{
+public:
+    Array3D(size_t nI, size_t nJ, size_t nK = 1):
+        _nI(nI), _nJ(nJ), _nK(nK), _vec( nI * nJ * nK )
+    {}
+
+    TYPE & operator()( size_t i, size_t j, size_t k = 0 ) {
+        if ( i < _nI && j < _nJ && k < _nK ) {
+            return _vec[ i + j * _nI + k * _nJ * _nI ];
+        }
+        throw std::out_of_range("Indexes out of range");
+    }
+
+    TYPE operator()( size_t i, size_t j, size_t k = 0 ) const {
+        if ( i < _nI && j < _nJ && k < _nK ) {
+            return _vec[ i + j * _nI + k * _nJ * _nI ];
+        }
+        throw std::out_of_range("Indexes out of range");
+    }
+
+private:
+    size_t _nI;
+    size_t _nJ;
+    size_t _nK;
+    std::vector<TYPE> _vec;
+};
+#endif // ARRAY3D_H

docs/GammaRayManual.docx

@@ -180,6 +180,31 @@ void CartesianGrid::setInfoFromGridParameter(GSLibParGrid *pg)
     this->setInfo( x0, y0, z0, dx, dy, dz, nx, ny, nz, rot, nreal, ndv, empty, empty2);
 }
 
+double CartesianGrid::dataIJK(uint column, uint i, uint j, uint k)
+{
+    uint dataRow = i + j*_nx + k*_ny*_nx;
+    return data( dataRow, column );
+}
+
+std::vector<std::complex<double> > CartesianGrid::getArray(int indexColumRealPart, int indexColumImaginaryPart)
+{
+    std::vector< std::complex<double> > result( _nx * _ny * _nz ); //[_nx][_ny][_nz]
+
+    for( uint k = 0; k < _nz; ++k)
+        for( uint j = 0; j < _ny; ++j)
+            for( uint i = 0; i < _nx; ++i){
+                double real = 0.0d;
+                double im = 0.0d;
+                if( indexColumRealPart >= 0 )
+                    real = dataIJK( indexColumRealPart, i, j, k);
+                if( indexColumImaginaryPart >= 0 )
+                    im = dataIJK( indexColumImaginaryPart, i, j, k);
+                result[i + j*_nx + k*_ny*_nx] = std::complex<double>(real, im);
+            }
+
+    return result;
+}
+
 bool CartesianGrid::canHaveMetaData()
 {
     return true;

domain/cartesiangrid.cpp

@@ -2,6 +2,7 @@
 #define CARTESIANGRID_H
 
 #include "datafile.h"
+#include <complex>
 
 class GSLibParGrid;
 
@@ -55,6 +56,23 @@ class CartesianGrid : public DataFile
     uint getNReal(){ return _nreal; }
     //@}
 
+    /**
+     * Returns a value from the data column (0 = 1st column) given a grid topological coordinate (IJK).
+     * @param i must be between 0 and NX-1.
+     * @param j must be between 0 and NY-1.
+     * @param k must be between 0 and NZ-1.
+     */
+    double dataIJK(uint column, uint i, uint j, uint k);
+
+    /** Creates a vector of complex numbers with the values taken from data columns.
+     *  Specify -1 to omit a column, which causes the repective part to be filled with zeros.
+     *  getArray(-1,-1) returns an array filled with zeroes.  The dimension of the array is that
+     *  of the Cartesian grid (getNX(), getNY(), getNZ()).
+     *  @param indexColumRealPart Column index (starting with 0) with the values for the real part.
+     *  @param indexColumRealPart Column index (starting with 0) with the values for the imaginary part.
+     */
+    std::vector< std::complex<double> > getArray( int indexColumRealPart, int indexColumImaginaryPart = -1 );
+
     //DataFile interface
 public:
     /** Cartesian grids never have declustering weights.  At least they are not supposed to be. */

domain/cartesiangrid.h

@@ -137,7 +137,7 @@ MainWindow::MainWindow(QWidget *parent) :
     setAcceptDrops( true );
 
     //show the 3D view widget
-    ui->frmContent->layout()->addWidget( new View3DWidget() );
+    ui->frmContent->layout()->addWidget( new View3DWidget( this ) );
 }
 
 MainWindow::~MainWindow()
@@ -449,6 +449,9 @@ void MainWindow::onProjectContextMenu(const QPoint &mouse_location)
                 _projectContextMenu->addMenu( m_subMenuClassifyWith );
                 _projectContextMenu->addAction("Soft indicator calibration...", this, SLOT(onSoftIndicatorCalib()) );
             }
+            if( parent_file->getFileType() == "CARTESIANGRID"  ){
+                _projectContextMenu->addAction("FFT", this, SLOT(onFFT()));
+            }
         }
     //two items were selected.  The context menu depends on the combination of items.
     } else if ( selected_indexes.size() == 2 ) {
@@ -1293,6 +1296,53 @@ void MainWindow::onFreeLoadedData()
     Application::instance()->getProject()->freeLoadedData();
 }
 
+void MainWindow::onFFT()
+{
+    //propose a name for the new grid to contain the FFT image
+    QString proposed_name = _right_clicked_attribute->getName() + "_FFT.dat";
+
+    //user enters the name for the new grid with FFT image
+    QString new_cg_name = QInputDialog::getText(this, "Name the new grid",
+                                             "Name for the grid with FFT image:", QLineEdit::Normal,
+                                             proposed_name );
+
+    //if the user canceled the input box
+    if ( new_cg_name.isEmpty() ){
+        //abort
+        return;
+    }
+
+    //the parent file is surely a CartesianGrid.
+    CartesianGrid *cg = (CartesianGrid*)_right_clicked_attribute->getContainingFile();
+
+    //get the array containing the data
+    std::vector< std::complex<double> > array = cg->getArray( _right_clicked_attribute->getAttributeGEOEASgivenIndex()-1 );
+
+    //run FFT
+    Util::fft3D( cg->getNX(),
+                 cg->getNY(),
+                 cg->getNZ(),
+                 array,
+                 FFTComputationMode::DIRECT);
+
+    //make a tmp file path
+    QString tmp_file_path = Application::instance()->getProject()->generateUniqueTmpFilePath("dat");
+
+    //crate a new cartesian grid pointing to the tmp path
+    CartesianGrid * new_cg = new CartesianGrid( tmp_file_path );
+
+    //set the geometry info based on the original grid
+    new_cg->setInfoFromOtherCG( cg, false );
+
+    //save the results in the project's tmp directory
+    Util::createGEOEASGrid( "Real part", "Imaginary part", array, tmp_file_path);
+
+    //import the saved file to the project
+    Application::instance()->getProject()->importCartesianGrid( new_cg, new_cg_name );
+
+    Application::instance()->logInfo("FFT 2D completed.");
+}
+
 void MainWindow::onCreateCategoryDefinition()
 {
     CategoryDefinition *cd = new CategoryDefinition("");

mainwindow.cpp

@@ -143,6 +143,7 @@ private slots:
     void onMapAs();
     void onSoftIndicatorCalib();
     void onFreeLoadedData();
+    void onFFT();
 
 private:
     /**