Inpaint gmca (#54)

* new option descreasing lambda * new option descreasing lambda * inpaint 1D Lambda last iter
CosmoStat · Nov 13, 2024 · 9d03f0f · 9d03f0f
1 parent 811bf17
commit 9d03f0f
Show file tree

Hide file tree

Showing 3 changed files with 33 additions and 6 deletions.
diff --git a/src/mc/mcmain2d/mr_gmca.cc b/src/mc/mcmain2d/mr_gmca.cc
@@ -165,7 +165,8 @@ Bool UseKnownColomn = False;
 Bool UseMask = False;
 int NbrKnownColumn = 0;   // Number of column which are known in the mixing matrix (a priori)
 fltarray MatKnownColumn; // Matrix containing the known column
-type_gmca_threshold_decrease TypeThresholdDecrease = GMCA_THRES_MAD; // Select the decreasing threshold strategy 
+// enum type_gmca_threshold_decrease {GMCA_THRES_LINEAR, GMCA_THRES_EXP, GMCA_THRES_MAD, GMCA_THRES_MOM};
+type_gmca_threshold_decrease TypeThresholdDecrease = GMCA_THRES_MAD; // Select the decreasing threshold strategy
 Bool Inpainting = False;            // If true, GMCA consider missing data, and apply an inpainting technique
 fltarray Mask;              // if Inpainting==True, Mask contains the available data Mask(i) = 1 or 0 (0 for no data).
 Bool PositiveMatrix = False;        // if True, the mixing matrice is assume to be positive
@@ -238,6 +239,13 @@ static void usage(char *argv[])
     fprintf(OUTMAN, "         [-K Last K-Mad]\n");
     fprintf(OUTMAN, "             Last value of K for K-Mad Thresholding.  \n");
     fprintf(OUTMAN, "         [-G Global Thresholding]\n");
+    fprintf(OUTMAN, "         [-E ThresholdDecreasingMethod]\n");
+    fprintf(OUTMAN, "             1: Linear Thresholding.  \n");
+    fprintf(OUTMAN, "             2: Exponential Thresholding.  \n");
+    fprintf(OUTMAN, "             3: Residual MAD Thresholding.  \n");
+    fprintf(OUTMAN, "             4: MOM Thresholding.  \n");
+    fprintf(OUTMAN, "             Default is 3. \n");
+
 //    fprintf(OUTMAN, "         [-O]\n");
 //    fprintf(OUTMAN, "                 Orthogonalization of the spectra\n");
     verbose_usage();
@@ -260,7 +268,7 @@ static void transinit(int argc, char *argv[])
 #endif   
 
     /* get options */
-    while ((c = GetOpt(argc,argv,"i:N:t:n:MS:PpI:A:ldm:L:DK:GUOvzZ:")) != -1) 
+    while ((c = GetOpt(argc,argv,"E:i:N:t:n:MS:PpI:A:ldm:L:DK:GUOvzZ:")) != -1)
     {
 	switch (c) 
         {
@@ -338,6 +346,22 @@ static void transinit(int argc, char *argv[])
 	            exit(-1);
  		}                
  		break;
+            case 'E':
+              /* -d <type> type of transform */
+              if (sscanf(OptArg,"%d",&c ) != 1)
+                      {
+                  fprintf(OUTMAN, "bad type of threshold decreasing method: %s\n", OptArg);
+                      exit(-1);
+
+              }
+                      if ((c > 0) && (c <= NRB_GMCA_THRES_DEC+1))
+                          TypeThresholdDecrease = (type_gmca_threshold_decrease) (c-1);
+                      else
+                      {
+                  fprintf(OUTMAN, "bad type of threshold decreasing method: %s\n", OptArg);
+                      exit(-1);
+               }
+               break;
        case 'M': Normalize = (Normalize == True) ? False: True; break;
  	   case 'v': Verbose = True; break;
 	   case 'n':
@@ -374,6 +398,7 @@ static void transinit(int argc, char *argv[])
 		    exit(-1);
 		}
  		break;
+
 #ifdef LARGE_BUFF
 	    case 'z':
 	        if (OptZ == True)

diff --git a/src/mga/libmga1d/MCA1D.cc b/src/mga/libmga1d/MCA1D.cc
@@ -418,13 +418,14 @@ void MCA1D::decomposition (fltarray& Signal)
 
     // if (LambdaTV > 0)   LambdaTV *= sqrt((double) 2.) / Nbr_Iter;
 
-   for (int Iter=0; Iter < Nbr_Iter; Iter++) 
+   for (int Iter=0; Iter < Nbr_Iter; Iter++)
    {
        if (Iter>0) 
        {
-          if (Linear == True) Lambda -= StepL; 
-          else Lambda =  LastSoftThreshold 
-                          + DeltaThreshold  *(1.-erf(2.8*Iter/ Nbr_Iter));
+          if (Linear == True)
+               Lambda =  LastSoftThreshold + (FirstSoftThreshold - LastSoftThreshold ) *  (Nbr_Iter-1-Iter) / (float) (Nbr_Iter-1);
+          else Lambda =  LastSoftThreshold
+                          + DeltaThreshold  *(1.-erf(2.8*(Iter-1)/ (double)(Nbr_Iter-1))) * cos((double)Iter /(double)(Nbr_Iter-1)*PI/2.);
           if (Lambda < LastSoftThreshold) Lambda = LastSoftThreshold;
        }
 

diff --git a/src/sparse/libsparse2d/GMCA.h b/src/sparse/libsparse2d/GMCA.h
@@ -41,6 +41,7 @@
 
 #define DEF_GMCA_MAX_ITER 100
 
+#define NRB_GMCA_THRES_DEC 4
 enum type_gmca_threshold_decrease {GMCA_THRES_LINEAR, GMCA_THRES_EXP, GMCA_THRES_MAD, GMCA_THRES_MOM};
 
 #define DEF_GMCA_DECREASE GMCA_THRES_LINEAR