enable gaussian fitting of multiple snapshots iteratively

gaussian_fit/input.cpp

@@ -2,6 +2,8 @@
 
 Input INPUT;
 ofstream ofs_running;
+ofstream ofs_local_max;
+ofstream ofs_predict;
 
 Input::Input()
 {
@@ -10,6 +12,7 @@ Input::Input()
     intensity_thre = 0.0;
     write_local_max = 0;
     write_predicted_intensity = 0;
+    top_pctg = 0.0;
 }
 
 Input::~Input(){}
@@ -26,6 +29,10 @@ void Input::read()
         else if (strcmp(keyword, "intensity_thre") == 0) read_value(ifs, intensity_thre);
         else if (strcmp(keyword, "write_local_max") == 0) read_value(ifs, write_local_max);
         else if (strcmp(keyword, "write_predicted_intensity") == 0) read_value(ifs, write_predicted_intensity);
+        else if (strcmp(keyword, "intensity_file") == 0) read_value(ifs, intensity_file);
+        else if (strcmp(keyword, "nsnapshot") == 0) read_value(ifs, nsnapshot);
+        else if (strcmp(keyword, "top_pctg") == 0) read_value(ifs, top_pctg);
+
     }
     return;
 }

gaussian_fit/input.h

@@ -25,6 +25,9 @@ class Input
     double intensity_thre;
     int write_local_max;
     int write_predicted_intensity;
+    string intensity_file;
+    int nsnapshot;
+    double top_pctg; // The top percentage will be used to calculate the intensity threshold.
 
 
 
@@ -40,5 +43,7 @@ class Input
 
 extern Input INPUT;
 extern ofstream ofs_running;
+extern ofstream ofs_local_max;
+extern ofstream ofs_predict;
 
 #endif

gaussian_fit/main.cpp

@@ -6,17 +6,42 @@ int main(int argc, char **argv)
 {
     INPUT.read();
     ofs_running.open("running.log");
+    ofs_local_max.open("local_max.txt");
+    ofs_predict.open("predicted_intensity.txt");
 
     Spot spot;
     ifstream ifs;
-    ifs.open("outIntensity.dat");
-    cout << INPUT.intensity_thre << endl;
-    spot.readin(ifs);
-    spot.identify_local_max();
-    spot.fit_gaussian();
-    spot.predict();
-    spot.clean();
+    ifs.open(INPUT.intensity_file);
+    //cout << INPUT.intensity_thre << endl;
+    for (int isnapshot=0; isnapshot<INPUT.nsnapshot; isnapshot++)
+    {
+        spot.readin(ifs);
+        if (INPUT.write_local_max > 0)
+        {
+            ofs_running << "isnapshot=" << isnapshot << endl;
+            cout << "isnapshot=" << isnapshot << endl;
+        }
+        if (INPUT.top_pctg > 0)
+        {
+            spot.calc_quantile();
+        }
+        else
+        {
+            spot.quantile = INPUT.intensity_thre;
+        }
+        ofs_running << "intensity_thre=" << spot.quantile << endl;
+        spot.identify_local_max();
+        if (INPUT.write_local_max)
+        {
+            ofs_local_max << "isnapshot=" << isnapshot << " " << spot.nlocal_max << endl;
+        }
+        spot.fit_gaussian();
+        spot.predict();
+        spot.clean();
+    }
     ofs_running.close();
+    ofs_local_max.close();
+    ofs_predict.close();
 
     return 0;
 }

gaussian_fit/spot.cpp

@@ -21,7 +21,7 @@ void Spot::identify_local_max()
     {
         for (int ix=window_width; ix<width-window_width; ix++)
         {
-            if (value[iy][ix] < this->local_max_thre)
+            if (value[iy][ix] < this->quantile)
             {
                 continue;
             }
@@ -47,7 +47,7 @@ void Spot::identify_local_max()
                 nlocal_max++;
                 if (INPUT.write_local_max > 0)
                 {
-                    ofs_running << "local_max " << nlocal_max-1 << ", ix=" << ix << ", iy=" << iy << endl;
+                    ofs_running << nlocal_max-1 << " " << ix << " " << iy << endl;
                     for (int iy1=-l; iy1<l+1; iy1++)
                     {
                         for (int ix1=-l; ix1<l+1; ix1++)
@@ -73,6 +73,7 @@ void Spot::fit_gaussian()
     {
         local_window[iy] = new double[window_width];
     }
+
     for (int imax=0; imax<nlocal_max; imax++)
     {
         for (int iy=-l; iy<l+1; iy++)
@@ -83,6 +84,11 @@ void Spot::fit_gaussian()
             }
         }
         gaussian[imax].fit(local_window, window_width);
+        if (INPUT.write_local_max > 0)
+        {
+            ofs_local_max << imax << " " << local_max_coord[imax][1] << " " << local_max_coord[imax][0] << " " << gaussian[imax].I << " "
+             << gaussian[imax].x0 << " " << gaussian[imax].y0 << " " << gaussian[imax].sigmax2 << " " << gaussian[imax].sigmay2 << " " << gaussian[imax].alpha << endl;
+        }
     }
     fit_ = true;
     for (int iy=0; iy<this->window_width; iy++)
@@ -109,6 +115,7 @@ void Spot::predict()
     for (int imax=0; imax<nlocal_max; imax++)
     {
         if (gaussian[imax].sigmax2 < 0 or gaussian[imax].sigmay2 < 0) continue;
+        if (abs(gaussian[imax].x0) > 2 or abs(gaussian[imax].y0) > 2) continue;
         for (int iy=0; iy<width; iy++)
         {
             for (int ix=0; ix<width; ix++)
@@ -119,17 +126,17 @@ void Spot::predict()
     }
     if (INPUT.write_predicted_intensity > 0)
     {
-        ofstream ofs("fitted_intensity.txt");
-        ofs << setprecision(8);
+        //ofstream ofs("fitted_intensity.txt");
+        ofs_predict << setprecision(8);
         for (int iy=0; iy<width; iy++)
         {
             for (int ix=0; ix<width; ix++)
             {
-                ofs << fitted_value[iy][ix] << " ";
+                ofs_predict << fitted_value[iy][ix] << " ";
             }
-            ofs << endl;
+            ofs_predict << endl;
         }
-        ofs.close();
+        //ofs.close();
     }
 
     return;
@@ -165,17 +172,18 @@ void Spot::clean()
         }
         delete[] fitted_value;
     }
+
     return;
 }
 
 void Spot::readin(ifstream &ifs)
 {
     assert(INPUT.width > 0);
     assert(INPUT.window_width > 0);
-    assert(INPUT.intensity_thre > 0);
+    assert(INPUT.intensity_thre > 0 or INPUT.top_pctg > 0);
     this->width = INPUT.width;
     this->window_width = INPUT.window_width;
-    this->local_max_thre = INPUT.intensity_thre;
+    //this->local_max_thre = INPUT.intensity_thre;
     this->value = new double*[width];
     for (int iy=0; iy<width; iy++) this->value[iy] = new double[width];
     for (int iy=0; iy<width; iy++)
@@ -187,4 +195,44 @@ void Spot::readin(ifstream &ifs)
     }
     read_ = true;
     return;
+}
+
+void Spot::calc_quantile()
+{
+    // This is an incomplete bubble sorting, since we only have to
+    // sort the first INPUT.top_pctg*width*width times to get the 
+    // quantile.
+    assert(this->read_);
+    int length = this->width*this->width;
+    double* intensity1D = new double[length];
+    for (int iy=0; iy<this->width; iy++)
+    {
+        for (int ix=0; ix<this->width; ix++)
+        {
+            intensity1D[iy*this->width + ix] = this->value[iy][ix];
+        }
+    }
+
+    int Ntop_pctg = int(INPUT.top_pctg * length)+1;
+    //cout << "Ntop_pctg=" << Ntop_pctg << endl;
+    for (int ii=0; ii<Ntop_pctg; ii++)
+    {
+        //cout << "ii=" << ii << endl;
+        for (int jj=length-1; jj>ii; jj--)
+        {
+            if (intensity1D[jj] > intensity1D[jj-1])
+            {
+                double tmp=intensity1D[jj-1];
+                intensity1D[jj-1] = intensity1D[jj];
+                intensity1D[jj] = tmp;
+            }
+        }
+    }
+    //cout << "choosed order = " << Ntop_pctg << endl;
+    this->quantile = intensity1D[Ntop_pctg-1];
+    INPUT.intensity_thre = this->quantile;
+    //cout << "intensity_thre=" << this->quantile << endl;
+
+    delete[] intensity1D;
+    return;
 }

gaussian_fit/spot.h

@@ -11,11 +11,12 @@ class Spot
 
     double** value;
     double** fitted_value;
+    double quantile;
     int width;
     int window_width;
     int nlocal_max;
     int nlocal_max_ub;
-    double local_max_thre;
+    //double local_max_thre;
     int** local_max_coord;
     Gaussian* gaussian;
 
@@ -24,11 +25,13 @@ class Spot
     bool fit_;
     bool predict_;
 
+
     void readin(ifstream &ifs);
     void identify_local_max();
     void fit_gaussian();
     void predict();
     void clean();
+    void calc_quantile();
 };
 
 #endif

gaussian_fit/test/Makefile → gaussian_fit/test/gaussian/Makefile

@@ -1,10 +1,9 @@
-VPATH=./\
-../\
+
 
 OBJS = test_gaussian_fit.o\
-../input.o\
-../gaussian.o\
-../spot.o\
+../../input.o\
+../../gaussian.o\
+../../spot.o\
 
 CC=mpiicc
 

gaussian_fit/test/quantile/INPUT

@@ -0,0 +1,4 @@
+width 1024
+window_width 3
+top_pctg 0.01
+intensity_file ../gaussian/outIntensity.dat

gaussian_fit/test/quantile/Makefile

@@ -0,0 +1,18 @@
+VPATH=./\
+../../\
+
+OBJS = test_quantile.o\
+../../input.o\
+../../gaussian.o\
+../../spot.o\
+
+CC=mpiicc
+
+.cpp.o:
+		${CC} ${OPTION} -c -O3 $< -o $@ -I../../ -I/usr/src/gtest/include/gtest -L/usr/src/gtest -lgtest -lpthread 
+
+quantile_test:${OBJS}
+		${CC} ${OPTION} -o quantile_test.exe ${OBJS} -I/usr/src/gtest/include/gtest -L/usr/src/gtest -lgtest -lpthread
+
+clean:
+		rm -f *.o *.exe

gaussian_fit/test/quantile/test.py

@@ -0,0 +1,3 @@
+import numpy as np
+intensity = np.genfromtxt("../gaussian/outIntensity.dat")
+print(np.max(intensity))

gaussian_fit/test/quantile/test_quantile.cpp

@@ -0,0 +1,26 @@
+#include "gaussian.h"
+#include "input.h"
+#include "spot.h"
+#include <gtest/gtest.h>
+
+TEST(Quantile_test, quantile)
+{
+    cout << 0 << endl;
+    INPUT.read();
+    Spot spot;
+    ifstream ifs;
+    cout << 1 << endl;
+    ifs.open(INPUT.intensity_file);
+    spot.readin(ifs);
+    cout << 2 << endl;
+    ifs.close();
+    spot.calc_quantile();
+    EXPECT_FLOAT_EQ(spot.quantile, 9.86382);
+}
+
+
+int main(int argc, char **argv)
+{
+    testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS(); 
+}