diff --git a/src/tools/bwccmdl.c b/src/tools/bwccmdl.c index c85ce50..85b139c 100644 --- a/src/tools/bwccmdl.c +++ b/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);