src/tools/bwccmdl.c

@@ -115,6 +115,14 @@
 #define TILES   0x0004
 #define USQLY   0x0002
 
+/*================================================================================================*/
+/**
+ * @details Macros to determine the minimum or maximum of two values.
+ */
+/*===========================================================================|====================*/
+#define MAX(x, y)             (((x) < (y))?(y):(x))     // Returns maximum between two values
+#define MIN(x, y)             (((x) > (y))?(y):(x))     // Returns minimum between two values
+
 /*================================================================================================*/
 /**
  * @details This macro defines a simple operation to remove a supplied deliminator from a string.
@@ -1059,6 +1067,105 @@ parse_opt(int                key,
   return EXIT_SUCCESS;
 }
 
+static uchar
+output_analysis(eas3_data *const ref_data, eas3_data *const org_data)
+{
+  /*-----------------------*\
+  ! DEFINE INT VARIABLES:   !
+  \*-----------------------*/
+  uint64_t            size;
+  uint64_t            i;
+
+  uint8_t             p;
+
+  /*-----------------------*\
+  ! DEFINE FLOAT VARIABLES: !
+  \*-----------------------*/
+  double              MSE, PSNR;
+  double              peakVal;
+  double              sum;
+
+  double              *dOrig, *dRef;
+  float               *fOrig, *fRef;
+
+  bwc_float           minVal, maxVal;
+
+  /*-----------------------*\
+  ! DEFINE STRUCTS:         !
+  \*-----------------------*/
+  eas3_param_names    *param_name;
+
+  dOrig = org_data->field.d;
+  dRef  = ref_data->field.d;
+
+  fOrig = org_data->field.f;
+  fRef  = ref_data->field.f;
+
+  if(org_data->params.ndim1 == ref_data->params.ndim1 &&
+     org_data->params.ndim2 == ref_data->params.ndim2 &&
+     org_data->params.ndim3 == ref_data->params.ndim3 &&
+     org_data->params.nts   == ref_data->params.nts   &&
+     org_data->params.npar  == ref_data->params.npar)
+    {
+      size = (uint64_t)ref_data->params.ndim1 * ref_data->params.ndim2 *
+                       ref_data->params.ndim3 * ref_data->params.nts;
+
+      peakVal = -1.7976931348623157e+308;
+      PSNR    =
+      MSE     = 0;
+
+      if(ref_data->param_names)
+        {
+          param_name = ref_data->param_names->root;
+          p = 0;
+          while(param_name != NULL)
+            {
+              minVal = 1.7976931348623157e+308;
+              maxVal = -1.7976931348623157e+308;
+
+              if(ref_data->params.accuracy == 2)
+                {
+                  for(i = 0; i < size; ++i)
+                    {
+                      minVal  = MIN(minVal, (double)dRef[i + p * size]);
+                      maxVal  = MAX(maxVal, (double)dRef[i + p * size]);
+
+                      sum     = ((double)dRef[i + p * size] - (double)dOrig[i + p * size]);
+
+                      MSE    += sum * sum;
+                    }
+                }
+              else if(ref_data->params.accuracy == 1)
+                {
+                  for(i = 0; i < size; ++i)
+                    {
+                      minVal  = MIN(minVal, (double)fRef[i + p * size]);
+                      maxVal  = MAX(maxVal, (double)fRef[i + p * size]);
+
+                      sum     = ((double)fRef[i + p * size] - (double)fOrig[i + p * size]);
+
+                      MSE    += sum * sum;
+                    }
+                }
+              peakVal = MAX(peakVal, maxVal - minVal);
+              param_name = param_name->next;
+              p++;
+            }
+        }
+
+      MSE /= (double)size * ref_data->params.npar;
+      PSNR = 20 * log10(peakVal/(2 * sqrt(MSE)));
+
+      printf("==============================================================\n");
+      printf("  Mean Square Error:             %*.2e\n", 22, MSE);
+      printf("  Peak Signal-to-Noise Ratio:    %*.2f\n", 22, PSNR);
+      printf("==============================================================\n");
+
+    }
+
+  return EXIT_SUCCESS;
+}
+
 /*================================================================================================*/
 /**
  * @details Initialize the argp struct. used to parse the command line arguments
@@ -1128,6 +1235,7 @@ int main(int argc, char *argv[])
   ! DEFINE STRUCTS:         !
   \*-----------------------*/
   eas3_data              *data = NULL;
+  eas3_data              *ref_data = NULL;
   cli_arguments           arguments = {0};
 
   /* Parse the command line arguments and invoke the appro- *
@@ -1543,7 +1651,67 @@ int main(int argc, char *argv[])
     }
   else if (arguments.mode == cli_anl)
     {
-      printf("Analysis\n");
+      /* Ingest the reference data input. */
+      if ((ref_data = read_eas3(arguments.ref)) == NULL)
+        {
+          error_handle = EXIT_FAILURE;
+          goto OUT;
+        }
+
+      /* Ingest the compressed data input. */
+      if ((fp = fopen(arguments.in, "r")) == NULL)
+        {
+          error_handle = EXIT_FAILURE;
+          printf(FINERROR);
+          goto OUT;
+        }
+
+      root   = ftell(fp);
+      fseek(fp, 0L,   SEEK_END);
+      Lfield = ftell(fp) - root;
+      fseek(fp, root, SEEK_SET);
+
+      /* Read the compressed data from the input file. */
+      input = calloc(Lfield, sizeof(uchar));
+      if (fread(input, sizeof(uchar), Lfield, fp) != Lfield)
+        {
+          error_handle = EXIT_FAILURE;
+          printf(RDERROR);
+          goto OUT;
+        }
+
+      /* Retrieve header information and allocate output buffer. */
+      header = bwc_open_header(input);
+      size = header->info.nX * header->info.nY * header->info.nZ *
+             header->info.nTS * header->info.nPar;
+      if(header->info.data_prec == bwc_precision_double)
+        {
+          output = calloc(size, sizeof(double));
+        }
+      else if(header->info.data_prec == bwc_precision_single)
+        {
+          output = calloc(size, sizeof(float));
+        }
+      bwc_close_header(header);
+
+      /* Initialize and run the decompression. */
+      stream = bwc_init_stream(input, output, comp);
+      coder  = bwc_alloc_decoder();
+      if (bwc_create_decompression(coder, stream, 0) == EXIT_FAILURE)
+        {
+          error_handle = EXIT_FAILURE;
+          goto OUT;
+        }
+      bwc_decompress(coder, stream);
+
+      /* Parse decompressed data into eas3 data structure *
+       * and write to the output file.                    */
+      data = calloc(1, sizeof(eas3_data));
+      bwc_to_eas3(stream, data);
+
+      output_analysis(ref_data, data);
+
+      goto OUT;
     }
   else if (arguments.mode == cli_inf)
     {
@@ -1674,6 +1842,9 @@ OUT:
   if (data != NULL)
     eas3_free_data(data);
 
+  if (ref_data != NULL)
+    eas3_free_data(ref_data);
+
   if (stream !=NULL)
     free(stream);