change interace function from t_cg to p_cg to 1) solve the efficiency problem for windows; 2) more straightforward to use probability than the chi-square critical value to users

Author: SuYe99

src/cxx/cold_flex.c

@@ -44,7 +44,8 @@ int cold_flex(
     int tmask_b2,               /* I: the band id used for tmask */
     int valid_num_scenes,       /* I: number of valid scenes  */
     int pos,                    /* I: the position id of pixel */
-    double pcg,                 /* I: probability threshold of change threshold  */
+    double t_cg,                /* I: threshold for identfying breaks */
+    double max_t_cg,            /* I: threshold for identfying outliers */
     int conse,                  /* I: consecutive observation number   */
     bool b_outputCM,            /* I: indicate if outputting change magnitudes for object-based cold, for cold only, it is the false */
     int starting_date,          /* I: (optional) the starting date of the whole dataset to enable reconstruct CM_date, all pixels for a tile should have the same date, only for b_outputCM is True */
@@ -67,8 +68,6 @@ int cold_flex(
     int i;
     char FUNC_NAME[] = "cold_flex";
     int result;
-    double t_cg = X2(nbands, pcg);
-
     if (valid_num_scenes == 0)
     {
         return (SUCCESS);
@@ -136,7 +135,7 @@ int cold_flex(
         /* standard_procedure for CCD                                 */
         /*                                                            */
         /**************************************************************/
-        result = stand_procedure_flex(valid_num_scenes, valid_date_array, ts_data, fmask_buf, nbands, id_range, t_cg, conse, b_outputCM, starting_date, rec_cg, num_fc, CM_OUTPUT_INTERVAL, CM_outputs,
+        result = stand_procedure_flex(valid_num_scenes, valid_date_array, ts_data, fmask_buf, nbands, id_range, t_cg, max_t_cg, conse, b_outputCM, starting_date, rec_cg, num_fc, CM_OUTPUT_INTERVAL, CM_outputs,
                                       CM_outputs_date, gap_days, tmask_b1, tmask_b2);
         //        result = stand_procedure_fixeddays(valid_num_scenes, valid_date_array, buf_b, buf_g, buf_r, buf_n, buf_s1, buf_s2, buf_t, fmask_buf, id_range,
         //                                 tcg, conse, b_outputCM, starting_date, rec_cg, num_fc, CM_OUTPUT_INTERVAL, CM_outputs,
@@ -183,7 +182,8 @@ int stand_procedure_flex(
     int64_t *fmask_buf, /* I:  mask-based time series  */
     int nbands,         /* I: input band number */
     int *id_range,
-    double tcg,            /* I: threshold of change threshold  */
+    double t_cg,           /* I: threshold for identfying breaks */
+    double max_t_cg,       /* I: threshold for identfying outliers */
     int conse,             /* I: consecutive observation number   */
     bool b_outputCM,       /* I: indicate if cold is running as the first step of object-based cold*/
     int starting_date,     /* I: the starting date of the whole dataset to enable reconstruct CM_date, all pixels for a tile should have the same date, only for b_outputCM is True */
@@ -271,7 +271,6 @@ int stand_procedure_flex(
     int current_CM_n;
     float prob_angle; // change probability for angle
     int i_span_skip = 0;
-    double t_cg_max = X2(nbands, 0.99999);
 
     fit_cft = (float **)allocate_2d_array(nbands, LASSO_COEFFS,
                                           sizeof(float));
@@ -445,7 +444,7 @@ int stand_procedure_flex(
     //        adj_TCG = tcg;
     //    }
     adj_conse = conse;
-    adj_TCG = tcg;
+    adj_TCG = t_cg;
 
     v_dif_mag = (float **)allocate_2d_array(nbands, adj_conse,
                                             sizeof(float));
@@ -1106,7 +1105,7 @@ int stand_procedure_flex(
                             }
                             break;
                         }
-                        else if (vec_magg[0] > t_cg_max) /* false change */
+                        else if (vec_magg[0] > max_t_cg) /* false change */
                         {
                             for (k = i_ini; k < end - 1; k++)
                             {
@@ -1851,7 +1850,7 @@ int stand_procedure_flex(
                     i_span_skip = 0;
                 }
 
-                else if (vec_mag[0] > t_cg_max) /*false change*/
+                else if (vec_mag[0] > max_t_cg) /*false change*/
                 {
                     /**********************************************/
                     /*                                            */

src/cxx/cold_flex.h

@@ -14,7 +14,8 @@ int cold_flex(
     int tmask_b2,               /* I: the band id used for tmask */
     int valid_num_scenes,       /* I: number of valid scenes  */
     int pos,                    /* I: the position id of pixel */
-    double pcg,                 /* I: probability threshold of change threshold  */
+    double t_cg,                /* I: threshold for identfying breaks */
+    double max_t_cg,            /* I: threshold for identfying outliers */
     int conse,                  /* I: consecutive observation number   */
     bool b_outputCM,            /* I: indicate if outputting change magnitudes for object-based cold, for cold only, it is the false */
     int starting_date,          /* I: (optional) the starting date of the whole dataset to enable reconstruct CM_date, all pixels for a tile should have the same date, only for b_outputCM is True */
@@ -33,7 +34,8 @@ int stand_procedure_flex(
     int64_t *fmask_buf, /* I:  mask-based time series  */
     int nbands,         /* I: input band number */
     int *id_range,
-    double tcg,            /* I: threshold of change threshold  */
+    double t_cg,           /* I: threshold for identfying breaks */
+    double max_t_cg,       /* I: threshold for identfying outliers */
     int conse,             /* I: consecutive observation number   */
     bool b_outputCM,       /* I: indicate if cold is running as the first step of object-based cold*/
     int starting_date,     /* I: the starting date of the whole dataset to enable reconstruct CM_date, all pixels for a tile should have the same date, only for b_outputCM is True */

src/cxx/sccd-desktop.c

@@ -773,8 +773,10 @@ int main(int argc, char *argv[])
             }
             else if (method == COLD_FLEX)
             {
-                int nbands = 5;
+                int nbands = 4;
                 rec_cg_flex = malloc(NUM_FC * sizeof(Output_t_flex));
+                double t_cg = X2(nbands, pcg);
+                double max_t_cg = X2(nbands, 0.99999);
                 if (rec_cg_flex == NULL)
                 {
                     RETURN_ERROR("ERROR allocating rec_cg_flex",
@@ -792,7 +794,7 @@ int main(int argc, char *argv[])
                     }
                 }
 
-                result = cold_flex(buf_stack, fmask_buf, sdate, nbands, 1, 4, valid_scene_count, i_col + 1, pcg, conse, b_outputCM, starting_date, rec_cg_flex, &num_fc, CM_OUTPUT_INTERVAL, CM_outputs,
+                result = cold_flex(buf_stack, fmask_buf, sdate, nbands, 1, 4, valid_scene_count, i_col + 1, t_cg, max_t_cg, conse, b_outputCM, starting_date, rec_cg_flex, &num_fc, CM_OUTPUT_INTERVAL, CM_outputs,
                                    CM_outputs_date, gap_days);
 
                 // snprintf (msg_str, sizeof(msg_str), "pixel %d COLD calculation finished\n", i_col+1);

src/python/pycold/_colds_cython.pyx

@@ -15,6 +15,7 @@ from libc.stdint cimport int32_t, int64_t, uint32_t, uint64_t
 from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
 from .common import reccg_dt, reccg_dt_flex, sccd_dt, nrtqueue_dt, nrtmodel_dt, pinpoint_dt
 np.import_array()
+from scipy.stats import chi2
 
 try:
     import typing
@@ -103,11 +104,12 @@ cdef extern from "../../cxx/cold.h":
                   short int *cm_outputs_date, double gap_days);
 
 cdef extern from "../../cxx/cold_flex.h":
-    cdef int32_t cold_flex(int64_t *ts_stack, int64_t *fmask_buf, int64_t *valid_date_array, int nbands, int tmask_b1, int tmask_b1, 
-                            int32_t valid_num_scenes, int32_t pos, double tcg, int32_t conse, 
-                            bool b_output_cm, int32_t starting_date, Output_t_flex *rec_cg,
-                            int32_t *num_fc, int32_t cm_output_interval, short int *cm_outputs,
-                            short int *cm_outputs_date, double gap_days);
+    cdef int32_t cold_flex(int64_t *ts_stack, int64_t *fmask_buf, int64_t *valid_date_array, int nbands, 
+                           int tmask_b1, int tmask_b1, int32_t valid_num_scenes, int32_t pos, 
+                           double tcg, double max_tcg, int32_t conse, bool b_output_cm, 
+                           int32_t starting_date, Output_t_flex *rec_cg, int32_t *num_fc, 
+                           int32_t cm_output_interval, short int *cm_outputs, 
+                           short int *cm_outputs_date, double gap_days);
 
 
 cdef extern from "../../cxx/cold.h":
@@ -601,7 +603,7 @@ cpdef _sccd_update(sccd_pack,
 
 cpdef _cold_detect_flex(np.ndarray[np.int64_t, ndim=1, mode='c'] dates, np.ndarray[np.int64_t, ndim=1, mode='c'] ts_stack,
                         np.ndarray[np.int64_t, ndim=1, mode='c'] qas, int32_t valid_num_scenes, int32_t nbands,
-                        double p_cg = 0.99, int32_t conse=6, int32_t pos=1, bint b_output_cm=False, 
+                        double t_cg, double max_t_cg, int32_t conse=6, int32_t pos=1, bint b_output_cm=False, 
                         int32_t starting_date=0, int32_t n_cm=0, int32_t cm_output_interval=0, 
                         double gap_days=365.25, int32_t tmask_b1=1, int32_t tmask_b2=1):
     """
@@ -614,7 +616,8 @@ cpdef _cold_detect_flex(np.ndarray[np.int64_t, ndim=1, mode='c'] dates, np.ndarr
         qas: 1d array, the QA cfmask bands. '0' - clear; '1' - water; '2' - shadow; '3' - snow; '4' - cloud
         valid_num_scenes: the number of valid observations
         nbands: the number of inputted bands
-        p_cg: threshold of change magnitude, default is chi2.ppf(0.99,5)
+        t_cg: threshold for change magnitude
+        max_t_cg: threshold for identifying outlier
         pos: position id of the pixel
         conse: consecutive observation number
         b_output_cm: bool, 'True' means outputting change magnitude and change magnitude dates, only for object-based COLD
@@ -657,9 +660,9 @@ cpdef _cold_detect_flex(np.ndarray[np.int64_t, ndim=1, mode='c'] dates, np.ndarr
     cdef short [:] cm_outputs_view = cm_outputs  # memory view
     cdef short [:] cm_outputs_date_view = cm_outputs_date  # memory view
     result = cold_flex(&ts_stack_view[0], &qas_view[0], &dates_view[0], nbands, tmask_b1, tmask_b2, 
-                        valid_num_scenes, pos, p_cg, conse, b_output_cm, starting_date, &rec_cg_view[0], 
-                        &num_fc, cm_output_interval, &cm_outputs_view[0], &cm_outputs_date_view[0], 
-                        gap_days)
+                        valid_num_scenes, pos, t_cg, max_t_cg, conse, b_output_cm, starting_date, 
+                        &rec_cg_view[0], &num_fc, cm_output_interval, &cm_outputs_view[0], 
+                        &cm_outputs_date_view[0], gap_days)
     if result != 0:
         raise RuntimeError("cold function fails for pos = {} ".format(pos))
     else: