diff --git a/components/isceobj/Sensor/ALOS.py b/components/isceobj/Sensor/ALOS.py index fc494e8..4a34295 100755 --- a/components/isceobj/Sensor/ALOS.py +++ b/components/isceobj/Sensor/ALOS.py @@ -25,7 +25,12 @@ # Author: Walter Szeliga #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - +#******************************************************************************** +#This program has been upgraded to handle the ALOS-1 PRF change issue. +# +#Cunren Liang, 12-December-2019 +#California Institute of Technology, Pasadena, CA +#******************************************************************************** import os @@ -114,8 +119,8 @@ class ALOS(Sensor): complex(-6.297074e-3,8.026685e-3), complex(7.217117e-1,-2.367683e-2)) - constants = Constants(iBias=15.5, - qBias=15.5, + constants = Constants(iBias=63.5, + qBias=63.5, pointingDirection=-1, antennaLength=8.9) @@ -499,7 +504,7 @@ class ALOS(Sensor): outputNow = self.output + appendStr if not (self._resampleFlag == ''): filein = self.output + '__tmp__' - self.imageFile.extractImage(filein) + self.imageFile.extractImage(filein, i) #image number start with 0 self.populateMetadata() objResample = None if(self._resampleFlag == 'single2dual'): @@ -513,7 +518,7 @@ class ALOS(Sensor): objResample.updateFrame(self.frame) os.remove(filein) else: - self.imageFile.extractImage(outputNow) + self.imageFile.extractImage(outputNow, i) #image number start with 0 self.populateMetadata() width = self.frame.getImage().getWidth() # self.readOrbitPulse(self._leaderFile,outputNow,width) @@ -721,7 +726,7 @@ class ImageFile(object): return None - def extractImage(self,output=None): + def extractImage(self,output=None, image_i=0): """For now, call a wrapped version of ALOS_pre_process""" productLevel = float(self.parent.leaderFile.sceneHeaderRecord.metadata[ 'Product level code']) @@ -731,13 +736,13 @@ class ImageFile(object): elif productLevel == 1.1: self.extractSLC(output) elif productLevel == 1.0: - self.extractRaw(output) + self.extractRaw(output, image_i) #image number start with 0 else: raise ValueError(productLevel) return None @use_api - def extractRaw(self,output=None): + def extractRaw(self,output=None, image_i=0): #if (self.numberOfSarChannels == 1): # print "Single Pol Data Found" # self.extractSinglePolImage(output=output) @@ -748,15 +753,16 @@ class ImageFile(object): if self.parent.leaderFile.sceneHeaderRecord.metadata[ 'Processing facility identifier'] == 'ERSDAC': prmDict = alos.alose_Py(self.parent._leaderFile, - self.parent._imageFile, output) + self.parent._imageFile, output, image_i) #image number start with 0 else: prmDict = alos.alos_Py(self.parent._leaderFile, - self.parent._imageFile, output) + self.parent._imageFile, output, image_i) #image number start with 0 pass # updated 07/24/2012 self.width = prmDict['NUMBER_BYTES_PER_LINE'] - 2 * prmDict['FIRST_SAMPLE'] - self.length = self.imageFDR.metadata['Number of lines per data set'] + #self.length = self.imageFDR.metadata['Number of lines per data set'] + self.length = prmDict['NUMBER_LINES'] self.prefix = self.imageFDR.metadata[ 'Number of bytes of prefix data per record'] self.suffix = self.imageFDR.metadata[ diff --git a/components/isceobj/Sensor/bindings/alosmodule.cpp b/components/isceobj/Sensor/bindings/alosmodule.cpp index 9378d0b..0bcaf0a 100644 --- a/components/isceobj/Sensor/bindings/alosmodule.cpp +++ b/components/isceobj/Sensor/bindings/alosmodule.cpp @@ -25,7 +25,8 @@ // Author: Giangi Sacco //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 #include @@ -68,11 +69,12 @@ PyInit_alos() PyObject *alos_C(PyObject* self,PyObject *args) { char *imageFile,*leaderFile,*outFile; + int image_i; struct PRM inputPRM; struct PRM outputPRM; struct GLOBALS globals; - if(!PyArg_ParseTuple(args,"sss",&leaderFile,&imageFile,&outFile)) + if(!PyArg_ParseTuple(args,"sssi",&leaderFile,&imageFile,&outFile, &image_i)) { return NULL; } @@ -96,7 +98,7 @@ PyObject *alos_C(PyObject* self,PyObject *args) globals.dopp = 0; // Are we calculating a doppler? globals.tbias = 0.0; // Is there a time bias to fix poor orbits? - ALOS_pre_process(inputPRM,&outputPRM,globals); + ALOS_pre_process(inputPRM,&outputPRM,globals,image_i); PyObject * dict = PyDict_New(); createDictionaryOutput(&outputPRM,dict); @@ -106,11 +108,12 @@ PyObject *alos_C(PyObject* self,PyObject *args) PyObject *alose_C(PyObject* self,PyObject *args) { char *imageFile,*leaderFile,*outFile; + int image_i; struct PRM inputPRM; struct PRM outputPRM; struct GLOBALS globals; - if(!PyArg_ParseTuple(args,"sss",&leaderFile,&imageFile,&outFile)) + if(!PyArg_ParseTuple(args,"sssi",&leaderFile,&imageFile,&outFile, &image_i)) { return NULL; } @@ -134,7 +137,7 @@ PyObject *alose_C(PyObject* self,PyObject *args) globals.dopp = 0; // Are we calculating a doppler? globals.tbias = 0.0; // Is there a time bias to fix poor orbits? - ALOS_pre_process(inputPRM,&outputPRM,globals); + ALOS_pre_process(inputPRM,&outputPRM,globals,image_i); PyObject * dict = PyDict_New(); createDictionaryOutput(&outputPRM,dict); @@ -184,6 +187,14 @@ PyObject * createDictionaryOutput(struct PRM * prm, PyObject * dict) intVal = PyLong_FromLong((long) prm->good_bytes); PyDict_SetItemString(dict,"NUMBER_GOOD_BYTES",intVal); Py_XDECREF(intVal); + + + + intVal = PyLong_FromLong((long) prm->num_lines); + PyDict_SetItemString(dict,"NUMBER_LINES",intVal); + Py_XDECREF(intVal); + + intVal = PyLong_FromLong((long) prm->num_rng_bins); PyDict_SetItemString(dict,"NUMBER_RANGE_BIN",intVal); Py_XDECREF(intVal); diff --git a/components/isceobj/Sensor/include/alosmodule.h b/components/isceobj/Sensor/include/alosmodule.h index 157e758..0e6840a 100644 --- a/components/isceobj/Sensor/include/alosmodule.h +++ b/components/isceobj/Sensor/include/alosmodule.h @@ -25,7 +25,8 @@ // Author: Giangi Sacco //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 #ifndef alosmodule_h @@ -42,7 +43,7 @@ extern "C" PyObject *alose_C(PyObject *self,PyObject *args); PyObject *createDictionaryOutput(struct PRM *prm,PyObject *dict); int ALOS_pre_process(struct PRM inputPRM, struct PRM *outputPRM, - struct GLOBALS globals); + struct GLOBALS globals, int image_i); } static PyMethodDef alos_methods[] = diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/ALOS_pre_process.c b/components/isceobj/Sensor/src/ALOS_pre_process/ALOS_pre_process.c index 332fe0b..790a95b 100644 --- a/components/isceobj/Sensor/src/ALOS_pre_process/ALOS_pre_process.c +++ b/components/isceobj/Sensor/src/ALOS_pre_process/ALOS_pre_process.c @@ -17,9 +17,25 @@ * added write_roi * *****************************************************************************/ +/******************************************************************************** +This program has been upgraded to handle the ALOS-1 PRF change issue. + +Cunren Liang, 12-December-2019 +California Institute of Technology, Pasadena, CA +*********************************************************************************/ + + #include "alosglobals.h" #include"image_sio.h" #include"lib_functions.h" +#include "resamp.h" + +//ALOS I or Q mean = 15.5, so get 15 or 16 randomly here +//#define ZERO_VALUE (char)(15 + rand() % 2) +//I changed the dynamic range when reading data +//ALOS I or Q mean = 63.5, so get 63 or 64 randomly here +#define ZERO_VALUE (char)(63 + rand() % 2) + char *USAGE = "\n\nUsage: ALOS_pre_process imagefile LEDfile [-near near_range] [-radius RE] [-swap] [-V] [-debug] [-quiet] \n" "\ncreates data.raw and writes out parameters (PRM format) to stdout\n" @@ -44,9 +60,8 @@ char *USAGE = "\n\nUsage: ALOS_pre_process imagefile LEDfile [-near near_rang "-tbias tbias correct the clock bias (positive value means plus)\n" "Example:\n" "ALOS_pre_process IMG-HH-ALPSRP050420840-H1.0__A LED-ALPSRP050420840-H1.0__A \n"; - -long read_ALOS_data (FILE *, FILE *, struct PRM *, long *); -long read_ALOSE_data (FILE *, FILE *, struct PRM *, long *); +long read_ALOS_data (FILE *, FILE *, struct PRM *, long *, struct resamp_info *, int); +long read_ALOSE_data (FILE *, FILE *, struct PRM *, long *, struct resamp_info *, int); void parse_ALOS_commands(int, char **, char *, struct PRM *); void set_ALOS_defaults(struct PRM *); void print_ALOS_defaults(struct PRM *); @@ -58,19 +73,54 @@ int write_roi(char *, FILE *, struct PRM, struct ALOS_ORB, char *); // ISCE stuff void init_from_PRM(struct PRM inPRM, struct PRM *prm); +int resamp_azimuth(char *slc2, char *rslc2, int nrg, int naz1, int naz2, double prf, double *dopcoeff, double *azcoef, int n, double beta); + int -ALOS_pre_process(struct PRM inputPRM, struct PRM *outputPRM,struct GLOBALS globals) +ALOS_pre_process(struct PRM inputPRM, struct PRM *outputPRM,struct GLOBALS globals, int image_i) //image number starts with 0!!! { FILE *imagefile, *ldrfile; -FILE *rawfile[11];//*prmfile[11]; -//char prmfilename[128]; +FILE *rawfile[11], *prmfile[11]; +char prmfilename[128]; int nPRF; long byte_offset; struct PRM prm; struct ALOS_ORB orb; char date[8]; + +////////////////////////////////////////////// + FILE *resampinfofile; + struct resamp_info rspi; + struct resamp_info rspi_new; + struct resamp_info rspi_pre[100];//maximum number of frames: 100 + int i, j, k; + double SC_clock_start; + double SC_clock_start_resamp; + double d2s = 24.0 * 3600.0; + double line_number_first; + int num_lines_out; + int gap_flag; + + double prf_all[200];//maximum number of prfs: 200 + int frame_counter_start_all[200];//maximum number of prfs: 200 + int nPRF_all;//maximum number of prfs: 200 + + double dopcoeff[4]; + double azcoef[2]; + int num_lines_max, j_max; + char outputfile[256]; + char *data; + FILE *first_prf_fp; + FILE *next_prf_fp; + int num_lines_append; + //int num_lines_gap; + int ret; +////////////////////////////////////////////// + + + //if (argc < 3) die (USAGE,""); + printf("reading image: %d\n", image_i); /* set flags */ dopp = globals.dopp; @@ -91,19 +141,21 @@ char date[8]; init_from_PRM(inputPRM,&prm); //parse_ALOS_commands(argc, argv, USAGE, &prm); - //if (verbose) print_ALOS_defaults(&prm); + /* apply an additional timing bias based on corner reflector analysis */ + //tbias = tbias - 0.0020835; + + if (verbose) print_ALOS_defaults(&prm); if (is_big_endian_() == -1) {swap = 1;fprintf(stderr,".... swapping bytes\n");} else {swap = 0;} /* IMG and LED files should exist already */ - if ((rawfile[0] = fopen(prm.input_file,"w")) == NULL) die("can't open ",prm.input_file); if ((imagefile = fopen(globals.imagefilename, "r")) == NULL) die ("couldn't open Level 1.0 IMG file \n",globals.imagefilename); if ((ldrfile = fopen(inputPRM.led_file, "r")) == NULL) die ("couldn't open LED file \n",inputPRM.led_file); /* if it exists, copy to prm structure */ - //strcpy(prm.led_file,leaderFilename); + strcpy(prm.led_file,inputPRM.led_file); /* name and open output files and header files for raw data (but input for later processing) */ - //get_files(&prm, &rawfile[nPRF], &prmfile[nPRF], prmfilename, argv[1], nPRF); + get_files(&prm, &rawfile[nPRF], &prmfile[nPRF], prmfilename, prm.input_file, nPRF); /* read sarleader; put info into prm; write log file if specified */ read_ALOS_sarleader(ldrfile, &prm, &orb); @@ -125,7 +177,7 @@ char date[8]; /* if prf changes, create new prm and data files */ if (nPRF > 0 ) { if (verbose) fprintf(stderr,"creating multiple files due to PRF change (*.%d) \n",nPRF+1); - //get_files(&prm, &rawfile[nPRF], &prmfile[nPRF], prmfilename, argv[1], nPRF); + get_files(&prm, &rawfile[nPRF], &prmfile[nPRF], prmfilename, prm.input_file, nPRF); } /* set the chirp extension to 500 if FBD fs = 16000000 */ @@ -141,21 +193,26 @@ char date[8]; returns byte offset if the PRF changes */ /* calculate parameters from orbit */ if (ALOS_format == 0) { - byte_offset = read_ALOS_data(imagefile, rawfile[nPRF], &prm, &byte_offset); + byte_offset = read_ALOS_data(imagefile, rawfile[nPRF], &prm, &byte_offset, &rspi, nPRF); } /* ERSDAC - use read_ALOSE_data */ if (ALOS_format == 1) { - byte_offset = read_ALOSE_data(imagefile, rawfile[nPRF], &prm, &byte_offset); + byte_offset = read_ALOSE_data(imagefile, rawfile[nPRF], &prm, &byte_offset, &rspi, nPRF); } - //fclose(rawfile[nPRF]); - // should work for AUIG and ERSDAC ALOS_ldr_orbit(&orb, &prm); /* calculate doppler from raw file */ + dopp=1;//always compute doppler for doing prf resampling if (dopp == 1) calc_dop(&prm); + //prf as a function of range in Hz + rspi.fd1[nPRF] = prm.fd1; + rspi.fdd1[nPRF] = prm.fdd1; + rspi.fddd1[nPRF] = prm.fddd1; + //rspi.input_file[nPRF] = prm.input_file; + strcpy(rspi.input_file[nPRF], prm.input_file); /* divide prf in half for quad_pol */ /* fix chirp slope */ @@ -172,7 +229,7 @@ char date[8]; if (force_slope == 1) prm.chirp_slope = forced_slope; /* write ascii output, SIO format */ - //put_sio_struct(prm, prmfile[nPRF]); + put_sio_struct(prm, prmfile[nPRF]); /* write roi_pac output */ if (roi) { @@ -184,6 +241,322 @@ char date[8]; nPRF++; } + rspi.nPRF=nPRF; + + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////// + printf("\nPRF details of frame: %d\n", image_i); + printf("+++++++++++++++++++++++++++++++++++++++++++++++\n"); + printf("number of PRF: %d\n", rspi.nPRF); + for (i = 0; i < rspi.nPRF; i++){ + printf("PRF %d prf (Hz): %f\n", i+1, rspi.prf[i]); + printf("PRF %d start time (days): %20.12f\n", i+1, rspi.SC_clock_start[i]); + printf("PRF %d frame_counter_start: %d\n", i+1, rspi.frame_counter_start[i]); + printf("PRF %d frame_counter_end: %d\n", i+1, rspi.frame_counter_end[i]); + printf("PRF %d number of lines: %d\n\n", i+1, rspi.frame_counter_end[i]-rspi.frame_counter_start[i]+1); + } + + //open parameter file for doing time adjustment and interpolation + if (image_i == 0){ + if((resampinfofile = fopen("resampinfo.bin", "wb")) == NULL) + die("couldn't open resampinfo file","resampinfo.bin"); + } + else{ + //open the file for reading and appending + if((resampinfofile = fopen("resampinfo.bin", "ab+")) == NULL) + die("couldn't open resampinfo file","resampinfo.bin"); + rewind(resampinfofile); + for(i=0; i < image_i; i++){ + if((fread((void *) &rspi_pre[i],sizeof(struct resamp_info), 1, resampinfofile)) != 1) + die("couldn't read from file","resampinfo.bin"); + } + } + + //get parameter from this image + memcpy(&rspi_pre[image_i], &rspi, sizeof(struct resamp_info)); + + //initialize rspi_new with resamp_info from reading the image, put the adjusted time in it + memcpy(&rspi_new, &rspi, sizeof(struct resamp_info)); + + //adjust start time + //unified PRF of the full track: first prf of first image + //start time of the full track: first line of first image + //only adjust time when the format is not ERSDAC format, becasue ERSDAC format does not have sdr.frame_counter. + printf("adjust start times\n"); + if(ALOS_format == 0){ + + if(image_i==0){ + //adjust start time of prf file i, no need to adjust for first prf + for(i = 1; i < rspi_pre[0].nPRF; i++){ + //time of the line just before the first line of first prf file + SC_clock_start = rspi_pre[0].SC_clock_start[0] - (1.0/rspi_pre[0].prf[0]) / d2s; + //time of the last line of each prf file + for(j = 0; j < i; j++){ + if(rspi_pre[0].num_lines[j] != rspi_pre[0].frame_counter_end[j] - rspi_pre[0].frame_counter_start[j] + 1) + fprintf(stderr, "\n\nWARNING: in image %d prf file %d, \ + number of lines in file: %d is not equal to that computed from frame_counter: %d\n\n", \ + 0, j, rspi_pre[0].num_lines[j], rspi_pre[0].frame_counter_end[j] - rspi_pre[0].frame_counter_start[j] + 1); + SC_clock_start += (rspi_pre[0].frame_counter_end[j]-rspi_pre[0].frame_counter_start[j]+1) * (1.0/rspi_pre[0].prf[j]) / d2s; + } + //time of the first line of current prf file + SC_clock_start += (1.0/rspi_pre[0].prf[i]) / d2s; + + printf("time adjustment result for image %d, prf %d:\n", image_i, i); + printf("+++++++++++++++++++++++++++++++++++++++++++++++\n"); + printf("original start time: %20.12f\n", rspi_pre[0].SC_clock_start[i]); + printf("adjusted start time: %20.12f\n", SC_clock_start); + printf("original - adjusted: %f (number of PRI)\n\n", (rspi_pre[0].SC_clock_start[i]-SC_clock_start)*d2s/(1.0/rspi_pre[0].prf[i])); + //update + rspi_new.SC_clock_start[i] = SC_clock_start; + } + } + else{ + //1. check to see if there is gap between images + gap_flag = 0; + for(i = 0; i < image_i; i++){ + if (rspi_pre[i].frame_counter_end[rspi_pre[i].nPRF-1] - rspi_pre[i+1].frame_counter_start[0] <= -2){ + fprintf(stderr, "\n\nWARNING: there are gaps between image %d and image: %d\n", i, i+1); + fprintf(stderr, "since we don't know the prf of these gap lines, we are not able to adjust starting time\n\n"); + gap_flag = 1; + } + } + //2. adjust start time + if(gap_flag == 0){ + //2.1 count the number of prf chunks in the full track including this image + nPRF_all = 0; + for(i = 0; i < image_i+1; i++){ + for(j = 0; j < rspi_pre[i].nPRF; j++){ + if((i==0) && (j==0)){ + prf_all[nPRF_all] = rspi_pre[i].prf[j]; + frame_counter_start_all[nPRF_all] = rspi_pre[i].frame_counter_start[j]; + nPRF_all += 1; + } + else{ + if((rspi_pre[i].frame_counter_start[j]>frame_counter_start_all[nPRF_all-1]) && (rspi_pre[i].prf[j]!=prf_all[nPRF_all-1])){ + prf_all[nPRF_all] = rspi_pre[i].prf[j]; + frame_counter_start_all[nPRF_all] = rspi_pre[i].frame_counter_start[j]; + nPRF_all += 1; + } + } + } + } + printf("number of prfs including this image: %d\n", nPRF_all); + printf("list of prfs:\n"); + for(i = 0; i < nPRF_all; i++){ + printf("frame_counter: %d, prf: %f\n", frame_counter_start_all[i], prf_all[i]); + } + + //2.2 adjust start time + for(i = 0; i < rspi_pre[image_i].nPRF; i++){ + //time of the line just before the first line of first prf file + //because the unite is day, the errors caused can be 0.042529743164777756 lines, should remove the integer or year part of SC_clock_start, or + //use second as unit in the future + SC_clock_start = rspi_pre[0].SC_clock_start[0] - (1.0/rspi_pre[0].prf[0]) / d2s; + //if there is only one PRF (no prf changes across all images) + if(nPRF_all == 1){ + SC_clock_start += (rspi_pre[image_i].frame_counter_start[0] - rspi_pre[0].frame_counter_start[0] + 1) * (1.0/rspi_pre[0].prf[0]) / d2s; + } + else{ + //find its position among the prfs, start from the second prf + for(j = 1; j < nPRF_all; j++){ + if(rspi_pre[image_i].frame_counter_start[i] < frame_counter_start_all[j]){ + //time of the last line of each prf chuck + for(k = 1; k < j; k++) + SC_clock_start += (frame_counter_start_all[k]-frame_counter_start_all[k-1]) * (1.0/prf_all[k-1]) / d2s; + SC_clock_start += (rspi_pre[image_i].frame_counter_start[i] - frame_counter_start_all[j-1] + 1) * (1.0/prf_all[j-1]) / d2s; + break; + } + else if(rspi_pre[image_i].frame_counter_start[i] == frame_counter_start_all[j]){ + //time of the last line of each prf chuck + for(k = 1; k < j; k++) + SC_clock_start += (frame_counter_start_all[k]-frame_counter_start_all[k-1]) * (1.0/prf_all[k-1]) / d2s; + SC_clock_start += (rspi_pre[image_i].frame_counter_start[i] - frame_counter_start_all[j-1] + 1) * (1.0/prf_all[j-1]) / d2s; + //extra pri of j-1 above, so remove it and add the pri of j + SC_clock_start += (1.0/prf_all[j]) / d2s - (1.0/prf_all[j-1]) / d2s; + break; + } + else{ + if(j == nPRF_all - 1){ + for(k = 1; k < j+1; k++) + SC_clock_start += (frame_counter_start_all[k]-frame_counter_start_all[k-1]) * (1.0/prf_all[k-1]) / d2s; + SC_clock_start += (rspi_pre[image_i].frame_counter_start[i] - frame_counter_start_all[j] + 1) * (1.0/prf_all[j]) / d2s; + break; + } + else{ + continue; + } + } + } + } + + //time of the first line of current prf file + printf("time adjustment result for image %d, prf %d:\n", image_i, i); + printf("+++++++++++++++++++++++++++++++++++++++++++++++\n"); + printf("original start time: %20.12f\n", rspi_pre[image_i].SC_clock_start[i]); + printf("adjusted start time: %20.12f\n", SC_clock_start); + printf("original - adjusted: %f (number of PRI)\n\n", (rspi_pre[image_i].SC_clock_start[i]-SC_clock_start)*d2s/(1.0/rspi_pre[image_i].prf[i])); + + //update + rspi_new.SC_clock_start[i] = SC_clock_start; + } + } + } + + } + + + // use parameters from rspi_pre[image_i], instead of rspi_new (to be updated) + //except rspi_new.SC_clock_start[i], since it was updated (more accurate) above. + printf("azimuth resampling\n"); + for(i = 0; i < rspi_pre[image_i].nPRF; i++){ + if((image_i==0)&&(i==0)) + continue; + //convention: line numbers start with zero + //line number of first line of first prf of first image: 0 + //line number of first line of this prf file + line_number_first = (rspi_new.SC_clock_start[i] - rspi_pre[0].SC_clock_start[0]) * d2s / (1.0 / rspi_pre[0].prf[0]); + //unit: pri of first prf of first image + num_lines_out = (rspi_pre[image_i].frame_counter_end[i] - rspi_pre[image_i].frame_counter_start[i] + 1) * (1.0/rspi_pre[image_i].prf[i]) / (1.0/rspi_pre[0].prf[0]); + + if((fabs(roundfi(line_number_first)-line_number_first)<0.1) && (rspi_pre[image_i].prf[i]==rspi_pre[0].prf[0])) + continue; + + //time of first line of the resampled image + SC_clock_start_resamp = rspi_pre[0].SC_clock_start[0] + roundfi(line_number_first) * (1.0 / rspi_pre[0].prf[0]) / d2s; + //compute offset parameters + //azcoef[0] + azpos * azcoef[1] + azcoef[0] = (SC_clock_start_resamp - rspi_new.SC_clock_start[i]) * d2s / (1.0/rspi_pre[image_i].prf[i]); + azcoef[1] = (1.0/rspi_pre[0].prf[0]) / (1.0/rspi_pre[image_i].prf[i]) - 1.0; + + //use doppler centroid frequency estimated from prf with maximum number of lines in this image + num_lines_max = -1; + j_max = -1; + for(j = 0; j < rspi_pre[image_i].nPRF; j++){ + if(rspi_pre[image_i].num_lines[j] >= num_lines_max){ + num_lines_max = rspi_pre[image_i].num_lines[j]; + j_max = j; + } + } + dopcoeff[0] = rspi_pre[image_i].fd1[j_max]; //average prf for alos-1 is good enough (calc_dop.c). + dopcoeff[1] = 0.0; + dopcoeff[2] = 0.0; + dopcoeff[3] = 0.0; + + //The filenames of all three files created for each prf, are from prm.input_file + //PRM: prm.input_file.PRM + (.prfno_start_from_1, if not first prf) + //data: prm.input_file + (.prfno_start_from_1, if not first prf) + //data after resampling: prm.input_file + (.prfno_start_from_1, if not first prf) + .interp + + sprintf(outputfile,"%s.interp", rspi_pre[image_i].input_file[i]); + //start interpolation + resamp_azimuth(rspi_pre[image_i].input_file[i], outputfile, rspi_pre[image_i].num_bins[i], num_lines_out, rspi_pre[image_i].num_lines[i], rspi_pre[image_i].prf[i], dopcoeff, azcoef, 9, 5.0); + + //update parameters + rspi_new.SC_clock_start[i] = SC_clock_start_resamp; + rspi_new.num_lines[i] = num_lines_out; + rspi_new.prf[i] = rspi_pre[0].prf[0]; + rspi_new.fd1[i] = dopcoeff[0]; + rspi_new.fdd1[i]= dopcoeff[1]; + rspi_new.fddd1[i]=dopcoeff[2]; + strcpy(rspi_new.input_file[i], outputfile); + } + + + //concatenate prfs: put all prfs to the first prf + // use parameters from rspi_new (updated), instead of rspi_pre[image_i] + if(rspi_new.nPRF > 1){ + + //prepare for appending subsequent prfs to first prf: open files and allocate memory + if((first_prf_fp = fopen(rspi_new.input_file[0], "ab")) == NULL) + die("can't open", rspi_new.input_file[0]); + //number of range samples in each prf is asummed to be same + if((data = (char *)malloc(2*sizeof(char)*rspi_new.num_bins[0])) == NULL) + die("can't allocate memory for data", ""); + + //append prf i + for(i = 1; i < rspi_new.nPRF; i++){ + //number of lines to be appended between frames if there are gaps + num_lines_append = (rspi_new.SC_clock_start[i] - rspi_new.SC_clock_start[0]) * d2s / (1.0/rspi_pre[0].prf[0]) - rspi_new.num_lines[0]; + if(num_lines_append >= 1){ + for(j = 0; j < num_lines_append; j++){ + for(k = 0; k < 2*rspi_new.num_bins[i]; k++) + data[k] = ZERO_VALUE; + if(fwrite((char *)data, 2*sizeof(char)*rspi_new.num_bins[i], 1, first_prf_fp) != 1) + die("can't write data to", rspi_new.input_file[0]); + } + rspi_new.num_lines[0] += num_lines_append; + } + + //append data from rspi_new.input_file[i] + if((next_prf_fp = fopen(rspi_new.input_file[i], "rb")) == NULL) + die("can't open", rspi_new.input_file[i]); + num_lines_append = 0; + for(j = 0; j < rspi_new.num_lines[i]; j++){ + if((rspi_new.SC_clock_start[i] + j * (1.0/rspi_pre[0].prf[0]) / d2s - rspi_new.SC_clock_start[0]) * d2s / (1.0/rspi_pre[0].prf[0]) >= rspi_new.num_lines[0]){ + if(fread((char *)data, 2*sizeof(char)*rspi_new.num_bins[i], 1, next_prf_fp) != 1) + die("can't read data from", rspi_new.input_file[i]); + if(fwrite((char *)data, 2*sizeof(char)*rspi_new.num_bins[i], 1, first_prf_fp) != 1) + die("can't write data to", rspi_new.input_file[0]); + num_lines_append += 1; + } + else{ + fseek(next_prf_fp, 2*sizeof(char)*rspi_new.num_bins[i], SEEK_CUR); + } + } + rspi_new.num_lines[0] += num_lines_append; + fclose(next_prf_fp); + } + free(data); + fclose(first_prf_fp); + } + + + //tidy up intermediate files + for(i = 0; i < rspi_pre[image_i].nPRF; i++){ + //if Return value = 0 then it indicates str1 is equal to str2. + ret = strcmp(rspi_new.input_file[i], rspi_pre[image_i].input_file[i]); + if(i == 0){ + if(ret != 0){ + //remove original + if(remove(rspi_pre[image_i].input_file[i]) != 0) + die("can't delete file", rspi_pre[image_i].input_file[i]); + //keep resampled and appended + if(rename(rspi_new.input_file[i], rspi_pre[image_i].input_file[i]) != 0) + die("can't rename file", rspi_new.input_file[i]); + } + } + else{ + //remove original + if(remove(rspi_pre[image_i].input_file[i]) != 0) + die("can't delete file", rspi_pre[image_i].input_file[i]); + //remove resampled + if(ret != 0){ + if(remove(rspi_new.input_file[i]) != 0) + die("can't delete file", rspi_new.input_file[i]); + } + } + } + + + //update prm + prm.prf = rspi_new.prf[0]; + prm.num_lines = rspi_new.num_lines[0]; + prm.SC_clock_start = rspi_new.SC_clock_start[0]; + prm.SC_clock_stop = prm.SC_clock_start + (prm.num_lines - 1) * (1.0/prm.prf) / d2s; + prm.fd1 = rspi_pre[image_i].fd1[j_max]; //average prf for alos-1 is good enough (calc_dop.c). + prm.fdd1 = 0.0; + prm.fddd1 =0.0; + + prm.xmi = 63.5; + prm.xmq = 63.5; + + //write to resampinfo.bin + if((fwrite((void *)&rspi_pre[image_i], sizeof(struct resamp_info), 1, resampinfofile)) != 1 ) + die("couldn't write to file", "resampinfo.bin"); + fclose(resampinfofile); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////// + if (orb.points != NULL) { @@ -198,12 +571,19 @@ void get_files(struct PRM *prm, FILE **rawfile, FILE **prmfile, char *prmfilenam /* name and open output file for raw data (but input for later processing) */ /* if more than 1 set of output files, append an integer (beginning with 2) */ - if (n == 0) { - sprintf(prm->input_file,"%s.raw", name); + //if (n == 0) { + // sprintf(prm->input_file,"%s.raw", name); + // sprintf(prmfilename,"%s.PRM", name); + //} else { + // sprintf(prm->input_file,"%s.raw.%d",name,n+1); + // sprintf(prmfilename,"%s.PRM.%d", name, n+1); + //} + if (n==0) { sprintf(prmfilename,"%s.PRM", name); + sprintf(prm->input_file,"%s",name); } else { - sprintf(prm->input_file,"%s.raw.%d",name,n+1); sprintf(prmfilename,"%s.PRM.%d", name, n+1); + sprintf(prm->input_file,"%s.%d",name,n+1); } /* now open the files */ @@ -212,4 +592,4 @@ void get_files(struct PRM *prm, FILE **rawfile, FILE **prmfile, char *prmfilenam if ((*prmfile = fopen(prmfilename, "w")) == NULL) die ("couldn't open output PRM file \n",prmfilename); } -/*------------------------------------------------------*/ + diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/SConscript b/components/isceobj/Sensor/src/ALOS_pre_process/SConscript index 5ed0b70..12ea6f5 100644 --- a/components/isceobj/Sensor/src/ALOS_pre_process/SConscript +++ b/components/isceobj/Sensor/src/ALOS_pre_process/SConscript @@ -6,11 +6,11 @@ Import('envSensorSrc1') package = envSensorSrc1['PACKAGE'] project = envSensorSrc1['PROJECT'] install = envSensorSrc1['PRJ_LIB_DIR'] -headerFiles = ['data_ALOS.h','data_ALOSE.h','image_sio.h','orbit_ALOS.h','sarleader_ALOS.h','sarleader_fdr.h','siocomplex.h'] -sourceFiles = ['ALOSE_orbits_utils.c','ALOS_ldr_orbit.c','ALOS_pre_process.c','calc_dop.c','hermite_c.c','init_from_PRM.c', - 'interpolate_ALOS_orbit.c','null_sio_struct.c','parse_ALOS_commands.c','polyfit.c','read_ALOSE_data.c', +headerFiles = ['data_ALOS.h','data_ALOSE.h','image_sio.h','orbit_ALOS.h','sarleader_ALOS.h','sarleader_fdr.h','siocomplex.h', 'resamp.h'] +sourceFiles = ['ALOSE_orbits_utils.c','ALOS_ldr_orbit.c','ALOS_pre_process.c','calc_dop.c','get_sio_struct.c','hermite_c.c','init_from_PRM.c', + 'interpolate_ALOS_orbit.c','null_sio_struct.c','parse_ALOS_commands.c','polyfit.c','put_sio_struct.c','read_ALOSE_data.c', 'read_ALOS_data.c','read_ALOS_sarleader.c','roi_utils.c','set_ALOS_defaults.c','siocomplex.c', - 'swap_ALOS_data_info.c','utils.c','write_ALOS_prm.c','readOrbitPulse.f','readOrbitPulseState.f','readOrbitPulseSetState.f'] + 'swap_ALOS_data_info.c','utils.c','write_ALOS_prm.c','readOrbitPulse.f','readOrbitPulseState.f','readOrbitPulseSetState.f', 'lib_array.c', 'lib_cpx.c', 'lib_file.c', 'lib_func.c', 'resamp_azimuth.c'] lib = envSensorSrc1.Library(target = 'alos', source = sourceFiles) envSensorSrc1.Install(install,lib) envSensorSrc1.Alias('install',install) diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/calc_dop.c b/components/isceobj/Sensor/src/ALOS_pre_process/calc_dop.c index 4a2d3db..039d097 100644 --- a/components/isceobj/Sensor/src/ALOS_pre_process/calc_dop.c +++ b/components/isceobj/Sensor/src/ALOS_pre_process/calc_dop.c @@ -12,6 +12,9 @@ * Date: * * *****************************************************************************/ +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 + #include "image_sio.h" #include "lib_functions.h" #include "siocomplex.h" @@ -23,8 +26,8 @@ void calc_dop(struct PRM *prm) long n; float *xr, *ac, *sg; double sumd; - fcomplex *ai, *bi, *ab; - fcomplex ctmp; + fcomplex_sio *ai, *bi, *ab; + fcomplex_sio ctmp; FILE *fin; fprintf(stderr,".... calculating doppler for %s\n",prm->input_file); @@ -40,9 +43,15 @@ void calc_dop(struct PRM *prm) ac = (float *) malloc(n*sizeof(float)); sg = (float *) malloc(n*sizeof(float)); - ai = (fcomplex *) malloc(n*sizeof(fcomplex)); - bi = (fcomplex *) malloc(n*sizeof(fcomplex)); - ab = (fcomplex *) malloc(2*n*sizeof(fcomplex)); + ai = (fcomplex_sio *) malloc(n*sizeof(fcomplex_sio)); + bi = (fcomplex_sio *) malloc(n*sizeof(fcomplex_sio)); + ab = (fcomplex_sio *) malloc(2*n*sizeof(fcomplex_sio)); + + + for(i = 0; i< n;i++){ + ab[i].r = 0; + ab[i].i = 0; + } /* read a line of data from fin (input file, chars) to ai (complex floats) */ fread(indata, sizeof(unsigned char), prm->bytes_per_line, fin); @@ -52,7 +61,7 @@ void calc_dop(struct PRM *prm) /* inefficient; could put loops inside each other */ for (i=prm->first_line; inum_lines-1; i++){ - if (i/2000 == i/2000.0) fprintf(stderr," Working on line %d \n",i); + //if (i/2000 == i/2000.0) fprintf(stderr," Working on line %d \n",i); fread(indata, sizeof(unsigned char), prm->bytes_per_line, fin); @@ -87,9 +96,10 @@ void calc_dop(struct PRM *prm) free(xr); free(ac); free(sg); free(ai); free(bi); free(ab); free(indata); + fprintf(stderr,"done\n"); } /*---------------------------------------------------*/ -void read_data(fcomplex *data, unsigned char *indata, int i, struct PRM *prm) +void read_data(fcomplex_sio *data, unsigned char *indata, int i, struct PRM *prm) { int ii ; diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/get_sio_struct.c b/components/isceobj/Sensor/src/ALOS_pre_process/get_sio_struct.c new file mode 100644 index 0000000..88afe26 --- /dev/null +++ b/components/isceobj/Sensor/src/ALOS_pre_process/get_sio_struct.c @@ -0,0 +1,201 @@ +/*--------------------------------------------------------------------*/ +/* + Read parameters into PRM structure from PRM file + Based on get_params by Evelyn J. Price + Modified by RJM +*/ +/*--------------------------------------------------------------------*/ + +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 + +#include "image_sio.h" +#include "lib_functions.h" + +/* +void get_sio_struct(FILE *, struct PRM *); +void get_string(char *, char *, char *, char *); +void get_int(char *, char *, char *, int *); +void get_double(char *, char *, char *, double *); +*/ + +void get_sio_struct(FILE *fh, struct PRM *s) { + char name[256], value[256]; + + debug = 0; + if (debug) { + fprintf(stderr, "get_sio_struct:\n"); + fprintf(stderr, "PRMname (PRM value) interpreted value\n"); + } + + while (fscanf(fh, "%s = %s \n", name, value) != EOF) { + + /* strings */ + if (strcmp(name, "input_file") == 0) + get_string(name, "input_file", value, s->input_file); + if (strcmp(name, "led_file") == 0) + get_string(name, "led_file", value, s->led_file); + if (strcmp(name, "out_amp_file") == 0) + get_string(name, "out_amp_file", value, s->out_amp_file); + if (strcmp(name, "out_data_file") == 0) + get_string(name, "out_data_file", value, s->out_data_file); + if (strcmp(name, "scnd_rng_mig") == 0) + get_string(name, "scnd_rng_mig", value, s->srm); + if (strcmp(name, "deskew") == 0) + get_string(name, "deskew", value, s->deskew); + if (strcmp(name, "Flip_iq") == 0) + get_string(name, "Flip_iq", value, s->iqflip); + if (strcmp(name, "offset_video") == 0) + get_string(name, "offset_video", value, s->offset_video); + if (strcmp(name, "ref_file") == 0) + get_string(name, "ref_file", value, s->ref_file); + if (strcmp(name, "SLC_file") == 0) + get_string(name, "SLC_file", value, s->SLC_file); + if (strcmp(name, "orbdir") == 0) + get_string(name, "orbdir", value, s->orbdir); + //if (strcmp(name, "lookdir") == 0) + // get_string(name, "lookdir", value, s->lookdir); + if (strcmp(name, "date") == 0) + get_string(name, "date", value, s->date); + + /* integers */ + if (strcmp(name, "nrows") == 0) + get_int(name, "nrows", value, &s->nrows); + if (strcmp(name, "num_lines") == 0) + get_int(name, "num_lines", value, &s->num_lines); + if (strcmp(name, "bytes_per_line") == 0) + get_int(name, "bytes_per_line", value, &s->bytes_per_line); + if (strcmp(name, "good_bytes_per_line") == 0) + get_int(name, "good_bytes_per_line", value, &s->good_bytes); + if (strcmp(name, "first_line") == 0) + get_int(name, "first_line", value, &s->first_line); + if (strcmp(name, "num_patches") == 0) + get_int(name, "num_patches", value, &s->num_patches); + if (strcmp(name, "first_sample") == 0) + get_int(name, "first_sample", value, &s->first_sample); + if (strcmp(name, "num_valid_az") == 0) + get_int(name, "num_valid_az", value, &s->num_valid_az); + if (strcmp(name, "SC_identity") == 0) + get_int(name, "SC_identity", value, &s->SC_identity); + if (strcmp(name, "chirp_ext") == 0) + get_int(name, "chirp_ext", value, &s->chirp_ext); + if (strcmp(name, "st_rng_bin") == 0) + get_int(name, "st_rng_bin", value, &s->st_rng_bin); + if (strcmp(name, "num_rng_bins") == 0) + get_int(name, "num_rng_bins", value, &s->num_rng_bins); + if (strcmp(name, "ref_identity") == 0) + get_int(name, "ref_identity", value, &s->ref_identity); + if (strcmp(name, "nlooks") == 0) + get_int(name, "nlooks", value, &s->nlooks); + if (strcmp(name, "rshift") == 0) + get_int(name, "rshift", value, &s->rshift); + if (strcmp(name, "ashift") == 0) + get_int(name, "ashift", value, &s->ashift); + /* backwards compatibility for xshift/rshift yshift/ashift */ + if (strcmp(name, "xshift") == 0) + get_int(name, "rshift", value, &s->rshift); + if (strcmp(name, "yshift") == 0) + get_int(name, "ashift", value, &s->ashift); + if (strcmp(name, "SLC_format") == 0) + get_int(name, "SLC_format", value, &s->SLC_format); + + /* doubles */ + if (strcmp(name, "SC_clock_start") == 0) + get_double(name, "SC_clock_start", value, &s->SC_clock_start); + if (strcmp(name, "SC_clock_stop") == 0) + get_double(name, "SC_clock_stop", value, &s->SC_clock_stop); + if (strcmp(name, "icu_start") == 0) + get_double(name, "icu_start", value, &s->icu_start); + //if (strcmp(name, "clock_start") == 0) + // get_double(name, "clock_start", value, &s->clock_start); + //if (strcmp(name, "clock_stop") == 0) + // get_double(name, "clock_stop", value, &s->clock_stop); + if (strcmp(name, "caltone") == 0) + get_double(name, "caltone", value, &s->caltone); + if (strcmp(name, "earth_radius") == 0) + get_double(name, "earth_radius", value, &s->RE); + if (strcmp(name, "equatorial_radius") == 0) + get_double(name, "equatorial_radius", value, &s->ra); + if (strcmp(name, "polar_radius") == 0) + get_double(name, "polar_radius", value, &s->rc); + if (strcmp(name, "SC_vel") == 0) + get_double(name, "SC_vel", value, &s->vel); + if (strcmp(name, "SC_height") == 0) + get_double(name, "SC_height", value, &s->ht); + if (strcmp(name, "SC_height_start") == 0) + get_double(name, "SC_height_start", value, &s->ht_start); + if (strcmp(name, "SC_height_end") == 0) + get_double(name, "SC_height_end", value, &s->ht_end); + if (strcmp(name, "near_range") == 0) + get_double(name, "near_range", value, &s->near_range); + if (strcmp(name, "PRF") == 0) + get_double(name, "PRF", value, &s->prf); + if (strcmp(name, "I_mean") == 0) + get_double(name, "I_mean", value, &s->xmi); + if (strcmp(name, "Q_mean") == 0) + get_double(name, "Q_mean", value, &s->xmq); + if (strcmp(name, "az_res") == 0) + get_double(name, "az_res", value, &s->az_res); + if (strcmp(name, "rng_samp_rate") == 0) + get_double(name, "rng_samp_rate", value, &s->fs); + if (strcmp(name, "chirp_slope") == 0) + get_double(name, "chirp_slope", value, &s->chirp_slope); + if (strcmp(name, "pulse_dur") == 0) + get_double(name, "pulse_dur", value, &s->pulsedur); + if (strcmp(name, "radar_wavelength") == 0) + get_double(name, "radar_wavelength", value, &s->lambda); + if (strcmp(name, "rng_spec_wgt") == 0) + get_double(name, "rng_spec_wgt", value, &s->rhww); + if (strcmp(name, "rm_rng_band") == 0) + get_double(name, "rm_rng_band", value, &s->pctbw); + if (strcmp(name, "rm_az_band") == 0) + get_double(name, "rm_az_band", value, &s->pctbwaz); + if (strcmp(name, "fd1") == 0) + get_double(name, "fd1", value, &s->fd1); + if (strcmp(name, "fdd1") == 0) + get_double(name, "fdd1", value, &s->fdd1); + if (strcmp(name, "fddd1") == 0) + get_double(name, "fddd1", value, &s->fddd1); + if (strcmp(name, "sub_int_r") == 0) + get_double(name, "sub_int_r", value, &s->sub_int_r); + if (strcmp(name, "sub_int_a") == 0) + get_double(name, "sub_int_a", value, &s->sub_int_a); + if (strcmp(name, "stretch_r") == 0) + get_double(name, "stretch_r", value, &s->stretch_r); + if (strcmp(name, "stretch_a") == 0) + get_double(name, "stretch_a", value, &s->stretch_a); + if (strcmp(name, "a_stretch_r") == 0) + get_double(name, "a_stretch_r", value, &s->a_stretch_r); + if (strcmp(name, "a_stretch_a") == 0) + get_double(name, "a_stretch_a", value, &s->a_stretch_a); + if (strcmp(name, "baseline_start") == 0) + get_double(name, "baseline_start", value, &s->baseline_start); + if (strcmp(name, "alpha_start") == 0) + get_double(name, "alpha_start", value, &s->alpha_start); + if (strcmp(name, "baseline_end") == 0) + get_double(name, "baseline_end", value, &s->baseline_end); + if (strcmp(name, "alpha_end") == 0) + get_double(name, "alpha_end", value, &s->alpha_end); + //if (strcmp(name, "SLC_scale") == 0) + // get_double(name, "SLC_scale", value, &s->SLC_scale); + } +} +/*--------------------------------------------------------------------------------*/ +void get_string(char *s1, char *name, char *value, char *s2) { + strcpy(s2, value); + if (debug == 1) + fprintf(stderr, " %s (%s) = %s\n", s1, name, value); +} +/*--------------------------------------------------------------------------------*/ +void get_int(char *s1, char *name, char *value, int *iparam) { + *iparam = atoi(value); + if (debug == 1) + fprintf(stderr, " %s (%s) = %s (%d)\n", s1, name, value, *iparam); +} +/*--------------------------------------------------------------------------------*/ +void get_double(char *s1, char *name, char *value, double *param) { + *param = atof(value); + if (debug == 1) + fprintf(stderr, " %s (%s) = %s (%lf)\n", s1, name, value, *param); +} +/*--------------------------------------------------------------------------------*/ diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/image_sio.h b/components/isceobj/Sensor/src/ALOS_pre_process/image_sio.h index 8437a45..c562d69 100644 --- a/components/isceobj/Sensor/src/ALOS_pre_process/image_sio.h +++ b/components/isceobj/Sensor/src/ALOS_pre_process/image_sio.h @@ -1,3 +1,6 @@ +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 + /* taken from soi.h */ #include #include @@ -33,9 +36,9 @@ #define NULL_DOUBLE -99999.9999 #define NULL_CHAR "XXXXXXXX" -typedef struct SCOMPLEX {short r,i;} scomplex; -typedef struct FCOMPLEX {float r,i;} fcomplex; -typedef struct DCOMPLEX {double r,i;} dcomplex; +typedef struct SCOMPLEX_SIO {short r,i;} scomplex_sio; +typedef struct FCOMPLEX_SIO {float r,i;} fcomplex_sio; +typedef struct DCOMPLEX_SIO {double r,i;} dcomplex_sio; struct PRM { char input_file[256]; @@ -121,6 +124,20 @@ struct PRM { double bpara; /* parallel baseline - added by RJM */ double bperp; /* perpendicular baseline - added by RJM */ }; +struct resamp_info { + //we assume there are no more than 20 prfs per image + int nPRF; //number of prfs, start with 1 + int frame_counter_start[20]; + int frame_counter_end[20]; + int num_lines[20]; + int num_bins[20]; + double prf[20]; + double SC_clock_start[20]; /* YYDDD.DDDD */ + double fd1[20]; + double fdd1[20]; + double fddd1[20]; + char input_file[20][256]; //we assume there are no more than 256 characters in the file name +}; /* offset_video off_vid chirp_ext nextend diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/lib_array.c b/components/isceobj/Sensor/src/ALOS_pre_process/lib_array.c new file mode 100644 index 0000000..b420da6 --- /dev/null +++ b/components/isceobj/Sensor/src/ALOS_pre_process/lib_array.c @@ -0,0 +1,575 @@ +////////////////////////////////////// +// Cunren Liang, NASA JPL/Caltech +// Copyright 2017 +////////////////////////////////////// + + +#include "resamp.h" + +/****************************************************************/ +/* allocating arrays */ +/****************************************************************/ + +signed char *vector_char(long nl, long nh) +/* allocate a signed char vector with subscript range v[nl..nh] */ +{ + signed char *v; + + v=(signed char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(signed char))); + if (!v){ + fprintf(stderr,"Error: cannot allocate 1-D vector\n"); + exit(1); + } + + return v-nl+NR_END; +} + +void free_vector_char(signed char *v, long nl, long nh) +/* free a signed char vector allocated with vector() */ +{ + free((FREE_ARG) (v+nl-NR_END)); +} + +unsigned char *vector_unchar(long nl, long nh) +/* allocate a unsigned char vector with subscript range v[nl..nh] */ +{ + unsigned char *v; + + v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); + if (!v){ + fprintf(stderr,"Error: cannot allocate 1-D vector\n"); + exit(1); + } + + return v-nl+NR_END; +} + +void free_vector_unchar(unsigned char *v, long nl, long nh) +/* free a unsigned char vector allocated with vector() */ +{ + free((FREE_ARG) (v+nl-NR_END)); +} + +int *vector_int(long nl, long nh) +/* allocate an int vector with subscript range v[nl..nh] */ +{ + int *v; + + v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); + if (!v) nrerror("Error: cannot allocate vector_int()"); + return v-nl+NR_END; +} + +void free_vector_int(int *v, long nl, long nh) +/* free an int vector allocated with ivector() */ +{ + free((FREE_ARG) (v+nl-NR_END)); +} + +float *vector_float(long nl, long nh) +/* allocate a float vector with subscript range v[nl..nh] */ +{ + float *v; + + v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); + if (!v){ + fprintf(stderr,"Error: cannot allocate 1-D vector\n"); + exit(1); + } + + return v-nl+NR_END; +} + +void free_vector_float(float *v, long nl, long nh) +/* free a float vector allocated with vector() */ +{ + free((FREE_ARG) (v+nl-NR_END)); +} + +double *vector_double(long nl, long nh) +/* allocate a double vector with subscript range v[nl..nh] */ +{ + double *v; + + v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); + if (!v){ + fprintf(stderr,"Error: cannot allocate 1-D vector\n"); + exit(1); + } + + return v-nl+NR_END; +} + +void free_vector_double(double *v, long nl, long nh) +/* free a double vector allocated with vector() */ +{ + free((FREE_ARG) (v+nl-NR_END)); +} + +fcomplex *vector_fcomplex(long nl, long nh) +/* allocate a fcomplex vector with subscript range v[nl..nh] */ +{ + fcomplex *v; + + v=(fcomplex *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(fcomplex))); + if (!v) nrerror("cannot allocate fcvector()"); + return v-nl+NR_END; +} + +void free_vector_fcomplex(fcomplex *v, long nl, long nh) +/* free a fcomplex vector allocated with fcvector() */ +{ + free((FREE_ARG) (v+nl-NR_END)); +} + +signed char **matrix_char(long nrl, long nrh, long ncl, long nch) +/* allocate a signed char matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + signed char **m; + + /* allocate pointers to rows */ + m=(signed char **) malloc((size_t)((nrow+NR_END)*sizeof(signed char*))); + if (!m) nrerror("Error: cannot allocate vector2d_float()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(signed char *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(signed char))); + if (!m[nrl]) nrerror("Error: cannot allocate vector2d_float()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +void free_matrix_char(signed char **m, long nrl, long nrh, long ncl, long nch) +/* free a signed char matrix allocated by matrix() */ +{ + free((FREE_ARG) (m[nrl]+ncl-NR_END)); + free((FREE_ARG) (m+nrl-NR_END)); +} + +unsigned char **matrix_unchar(long nrl, long nrh, long ncl, long nch) +/* allocate a unsigned char matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + unsigned char **m; + + /* allocate pointers to rows */ + m=(unsigned char **) malloc((size_t)((nrow+NR_END)*sizeof(unsigned char*))); + if (!m) nrerror("Error: cannot allocate vector2d_float()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(unsigned char *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(unsigned char))); + if (!m[nrl]) nrerror("Error: cannot allocate vector2d_float()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +void free_matrix_unchar(unsigned char **m, long nrl, long nrh, long ncl, long nch) +/* free a unsigned char matrix allocated by matrix() */ +{ + free((FREE_ARG) (m[nrl]+ncl-NR_END)); + free((FREE_ARG) (m+nrl-NR_END)); +} + +float **matrix_float(long nrl, long nrh, long ncl, long nch) +/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + float **m; + + /* allocate pointers to rows */ + m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); + if (!m) nrerror("Error: cannot allocate vector2d_float()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))); + if (!m[nrl]) nrerror("Error: cannot allocate vector2d_float()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +void free_matrix_float(float **m, long nrl, long nrh, long ncl, long nch) +/* free a float matrix allocated by matrix() */ +{ + free((FREE_ARG) (m[nrl]+ncl-NR_END)); + free((FREE_ARG) (m+nrl-NR_END)); +} + +double **matrix_double(long nrl, long nrh, long ncl, long nch) +/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + double **m; + + /* allocate pointers to rows */ + m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); + if (!m) nrerror("Error: cannot allocate vector2d_double()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); + if (!m[nrl]) nrerror("Error: cannot allocate vector2d_double()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +void free_matrix_double(double **m, long nrl, long nrh, long ncl, long nch) +/* free a double matrix allocated by matrix() */ +{ + free((FREE_ARG) (m[nrl]+ncl-NR_END)); + free((FREE_ARG) (m+nrl-NR_END)); +} + + + +/****************************************************************/ +/* allocating C-style arrays */ +/****************************************************************/ + +FILE **array1d_FILE(long nc){ + + FILE **fv; + + fv = (FILE **)malloc(nc * sizeof(FILE *)); + if(!fv){ + fprintf(stderr,"Error: cannot allocate 1-D FILE array\n"); + exit(1); + } + + return fv; +} + +void free_array1d_FILE(FILE **fv){ + free(fv); +} + +signed char *array1d_char(long nc){ + + signed char *fv; + + fv = (signed char*) malloc(nc * sizeof(signed char)); + if(!fv){ + fprintf(stderr,"Error: cannot allocate 1-D signed char vector\n"); + exit(1); + } + + return fv; +} + +void free_array1d_char(signed char *fv){ + free(fv); +} + +unsigned char *array1d_unchar(long nc){ + + unsigned char *fv; + + fv = (unsigned char*) malloc(nc * sizeof(unsigned char)); + if(!fv){ + fprintf(stderr,"Error: cannot allocate 1-D unsigned char vector\n"); + exit(1); + } + + return fv; +} + +void free_array1d_unchar(unsigned char *fv){ + free(fv); +} + +int *array1d_int(long nc){ + + int *fv; + + fv = (int*) malloc(nc * sizeof(int)); + if(!fv){ + fprintf(stderr,"Error: cannot allocate 1-D int array\n"); + exit(1); + } + + return fv; +} + +void free_array1d_int(int *fv){ + free(fv); +} + +float *array1d_float(long nc){ + + float *fv; + + fv = (float*) malloc(nc * sizeof(float)); + if(!fv){ + fprintf(stderr,"Error: cannot allocate 1-D float vector\n"); + exit(1); + } + + return fv; +} + +void free_array1d_float(float *fv){ + free(fv); +} + +double *array1d_double(long nc){ + + double *fv; + + fv = (double*) malloc(nc * sizeof(double)); + if(!fv){ + fprintf(stderr,"Error: cannot allocate 1-D double vector\n"); + exit(1); + } + + return fv; +} + +void free_array1d_double(double *fv){ + free(fv); +} + +fcomplex *array1d_fcomplex(long nc){ + + fcomplex *fcv; + + fcv = (fcomplex*) malloc(nc * sizeof(fcomplex)); + if(!fcv){ + fprintf(stderr,"Error: cannot allocate 1-D float complex vector\n"); + exit(1); + } + + return fcv; + +} + +void free_array1d_fcomplex(fcomplex *fcv){ + free(fcv); +} + +dcomplex *array1d_dcomplex(long nc){ + + dcomplex *fcv; + + fcv = (dcomplex*) malloc(nc * sizeof(dcomplex)); + if(!fcv){ + fprintf(stderr,"Error: cannot allocate 1-D double complex vector\n"); + exit(1); + } + + return fcv; + +} + +void free_array1d_dcomplex(dcomplex *fcv){ + free(fcv); +} + +signed char **array2d_char(long nl, long nc){ +/* allocate a signed char 2-D matrix */ + + signed char **m; + int i; + + /* allocate pointers to rows */ + m = (signed char **) malloc(nl * sizeof(signed char *)); + if(!m){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* allocate rows */ + m[0] = (signed char*) malloc(nl * nc * sizeof(signed char)); + if(!m[0]){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* set pointers */ + for(i = 1; i < nl; i++){ + m[i] = m[i-1] + nc; + } + + return m; +} + +void free_array2d_char(signed char **m){ +/* free a signed char matrix allocated by farray2d() */ + free(m[0]); + free(m); +} + +unsigned char **array2d_unchar(long nl, long nc){ +/* allocate a unsigned char 2-D matrix */ + + unsigned char **m; + int i; + + /* allocate pointers to rows */ + m = (unsigned char **) malloc(nl * sizeof(unsigned char *)); + if(!m){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* allocate rows */ + m[0] = (unsigned char*) malloc(nl * nc * sizeof(unsigned char)); + if(!m[0]){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* set pointers */ + for(i = 1; i < nl; i++){ + m[i] = m[i-1] + nc; + } + + return m; +} + +void free_array2d_unchar(unsigned char **m){ +/* free a signed unchar matrix allocated by farray2d() */ + free(m[0]); + free(m); +} + +float **array2d_float(long nl, long nc){ +/* allocate a float 2-D matrix */ + + float **m; + int i; + + /* allocate pointers to rows */ + m = (float **) malloc(nl * sizeof(float *)); + if(!m){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* allocate rows */ + m[0] = (float*) malloc(nl * nc * sizeof(float)); + if(!m[0]){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* set pointers */ + for(i = 1; i < nl; i++){ + m[i] = m[i-1] + nc; + } + + return m; +} + +void free_array2d_float(float **m){ +/* free a float matrix allocated by farray2d() */ + free(m[0]); + free(m); +} + +double **array2d_double(long nl, long nc){ +/* allocate a double 2-D matrix */ + + double **m; + int i; + + /* allocate pointers to rows */ + m = (double **) malloc(nl * sizeof(double *)); + if(!m){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* allocate rows */ + m[0] = (double*) malloc(nl * nc * sizeof(double)); + if(!m[0]){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* set pointers */ + for(i = 1; i < nl; i++){ + m[i] = m[i-1] + nc; + } + + return m; +} + +void free_array2d_double(double **m){ +/* free a double matrix allocated by farray2d() */ + free(m[0]); + free(m); +} + +fcomplex **array2d_fcomplex(long nl, long nc){ +/* allocate a fcomplex 2-D matrix */ + + fcomplex **m; + int i; + + /* allocate pointers to rows */ + m = (fcomplex **) malloc(nl * sizeof(fcomplex *)); + if(!m){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* allocate rows */ + m[0] = (fcomplex*) malloc(nl * nc * sizeof(fcomplex)); + if(!m[0]){ + fprintf(stderr,"Error: cannot allocate 2-D matrix\n"); + exit(1); + } + + /* set pointers */ + for(i = 1; i < nl; i++){ + m[i] = m[i-1] + nc; + } + + return m; +} + +void free_array2d_fcomplex(fcomplex **m){ +/* free a fcomplex matrix allocated by fcarray2d() */ + free(m[0]); + free(m); +} + + +/****************************************************************/ +/* handling error */ +/****************************************************************/ + +void nrerror(char error_text[]) +/* Numerical Recipes standard error handler */ +{ + fprintf(stderr,"Numerical Recipes run-time error...\n"); + fprintf(stderr,"%s\n",error_text); + fprintf(stderr,"...now exiting to system...\n"); + exit(1); +} diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/lib_cpx.c b/components/isceobj/Sensor/src/ALOS_pre_process/lib_cpx.c new file mode 100644 index 0000000..956abb6 --- /dev/null +++ b/components/isceobj/Sensor/src/ALOS_pre_process/lib_cpx.c @@ -0,0 +1,72 @@ +////////////////////////////////////// +// Cunren Liang, NASA JPL/Caltech +// Copyright 2017 +////////////////////////////////////// + + +#include "resamp.h" + +// complex operations +fcomplex cmul(fcomplex a, fcomplex b) +{ + fcomplex c; + c.re=a.re*b.re-a.im*b.im; + c.im=a.im*b.re+a.re*b.im; + return c; +} + +fcomplex cconj(fcomplex z) +{ + fcomplex c; + c.re=z.re; + c.im = -z.im; + return c; +} + +fcomplex cadd(fcomplex a, fcomplex b) +{ + fcomplex c; + c.re=a.re+b.re; + c.im=a.im+b.im; + return c; +} + +float xcabs(fcomplex z) +{ + float x,y,ans,temp; + x=fabs(z.re); + y=fabs(z.im); + if (x == 0.0) + ans=y; + else if (y == 0.0) + ans=x; + else if (x > y) { + temp=y/x; + ans=x*sqrt(1.0+temp*temp); + } else { + temp=x/y; + ans=y*sqrt(1.0+temp*temp); + } + return ans; +} + +float cphs(fcomplex z){ + float ans; + + if(z.re == 0.0 && z.im == 0.0) + ans = 0.0; + else + ans = atan2(z.im, z.re); + + return ans; +//it seems that there is no need to add the if clause +//do a test: +// printf("%12.4f, %12.4f, %12.4f, %12.4f, %12.4f\n", \ +// atan2(0.0, 1.0), atan2(1.0, 0.0), atan2(0.0, -1.0), atan2(-1.0, 0.0), atan2(0.0, 0.0)); +//output: +// 0.0000, 1.5708, 3.1416, -1.5708, 0.0000 +} + + + + diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/lib_file.c b/components/isceobj/Sensor/src/ALOS_pre_process/lib_file.c new file mode 100644 index 0000000..c62a953 --- /dev/null +++ b/components/isceobj/Sensor/src/ALOS_pre_process/lib_file.c @@ -0,0 +1,43 @@ +////////////////////////////////////// +// Cunren Liang, NASA JPL/Caltech +// Copyright 2017 +////////////////////////////////////// + + +#include "resamp.h" + +FILE *openfile(char *filename, char *pattern){ + FILE *fp; + + fp=fopen(filename, pattern); + if (fp==NULL){ + fprintf(stderr,"Error: cannot open file: %s\n", filename); + exit(1); + } + + return fp; +} + +void readdata(void *data, size_t blocksize, FILE *fp){ + if(fread(data, blocksize, 1, fp) != 1){ + fprintf(stderr,"Error: cannot read data\n"); + exit(1); + } +} + +void writedata(void *data, size_t blocksize, FILE *fp){ + if(fwrite(data, blocksize, 1, fp) != 1){ + fprintf(stderr,"Error: cannot write data\n"); + exit(1); + } +} + +long file_length(FILE* fp, long width, long element_size){ + long length; + + fseeko(fp,0L,SEEK_END); + length = ftello(fp) / element_size / width; + rewind(fp); + + return length; +} diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/lib_func.c b/components/isceobj/Sensor/src/ALOS_pre_process/lib_func.c new file mode 100644 index 0000000..762b66b --- /dev/null +++ b/components/isceobj/Sensor/src/ALOS_pre_process/lib_func.c @@ -0,0 +1,275 @@ +////////////////////////////////////// +// Cunren Liang, NASA JPL/Caltech +// Copyright 2017 +////////////////////////////////////// + + +#include "resamp.h" + +long next_pow2(long a){ + long i; + long x; + + x = 2; + while(x < a){ + x *= 2; + } + + return x; +} + +void circ_shift(fcomplex *in, int na, int nc){ + + int i; + int ncm; + + ncm = nc%na; + + if(ncm < 0){ + for(i = 0; i < abs(ncm); i++) + left_shift(in, na); + } + else if(ncm > 0){ + for(i = 0; i < ncm; i++) + right_shift(in, na); + } + else{ //ncm == 0, no need to shift + i = 0; + } +} + +void left_shift(fcomplex *in, int na){ + + int i; + fcomplex x; + + if(na < 1){ + fprintf(stderr, "Error: array size < 1\n\n"); + exit(1); + } + else if(na > 1){ + x.re = in[0].re; + x.im = in[0].im; + for(i = 0; i <= na - 2; i++){ + in[i].re = in[i+1].re; + in[i].im = in[i+1].im; + } + in[na-1].re = x.re; + in[na-1].im = x.im; + } + else{ //na==1, no need to shift + i = 0; + } +} + +void right_shift(fcomplex *in, int na){ + + int i; + fcomplex x; + + if(na < 1){ + fprintf(stderr, "Error: array size < 1\n\n"); + exit(1); + } + else if(na > 1){ + x.re = in[na-1].re; + x.im = in[na-1].im; + for(i = na - 1; i >= 1; i--){ + in[i].re = in[i-1].re; + in[i].im = in[i-1].im; + } + in[0].re = x.re; + in[0].im = x.im; + } + else{ //na==1, no need to shift + i = 0; + } +} + +int roundfi(float a){ + int b; + + if(a > 0) + b = (int)(a + 0.5); + else if (a < 0) + b = (int)(a - 0.5); + else + b = a; + + return b; +} + +void sinc(int n, int m, float *coef){ + + int i; + int hmn; + + hmn = n * m / 2; + + for(i=-hmn; i<=hmn; i++){ + if(i != 0){ + coef[i] = sin(PI * i / m) / (PI * i / m); + //coef[i] = sin(pi * i / m) / (pi * i / m); + } + else{ + coef[i] = 1.0; + } + } + +} + +//kaiser() is equivalent to kaiser2() +//it is created to just keep the same style of sinc(). +void kaiser(int n, int m, float *coef, float beta){ + + int i; + int hmn; + float a; + + hmn = n * m / 2; + + for(i = -hmn; i <= hmn; i++){ + a = 1.0 - 4.0 * i * i / (n * m) / (n * m); + coef[i] = bessi0(beta * sqrt(a)) / bessi0(beta); + } +} + +void kaiser2(float beta, int n, float *coef){ + + int i; + int hn; + float a; + + hn = (n - 1) / 2; + + for(i = -hn; i<=hn; i++){ + a = 1.0 - 4.0 * i * i / (n - 1.0) / (n - 1.0); + coef[i] = bessi0(beta * sqrt(a)) / bessi0(beta); + } +} + +void bandpass_filter(float bw, float bc, int n, int nfft, int ncshift, float beta, fcomplex *filter){ + + int i; + float *kw; + int hn; + fcomplex bwx, bcx; + + hn = (n-1)/2; + + if(n > nfft){ + fprintf(stderr, "Error: fft length too small!\n\n"); + exit(1); + } + if(abs(ncshift) > nfft){ + fprintf(stderr, "Error: fft length too small or shift too big!\n\n"); + exit(1); + } + + //set all the elements to zero + for(i = 0; i < nfft; i++){ + filter[i].re = 0.0; + filter[i].im = 0.0; + } + + //calculate kaiser window + kw = vector_float(-hn, hn); + kaiser2(beta, n, kw); + + //calculate filter + for(i = -hn; i <= hn; i++){ + bcx.re = cos(bc * 2.0 * PI * i); + bcx.im = sin(bc * 2.0 * PI * i); + + if(i == 0){ + bwx.re = 1.0; + bwx.im = 0.0; + } + else{ + bwx.re = sin(bw * PI * i) / (bw * PI * i); + bwx.im = 0.0; + } + + filter[i+hn] = cmul(bcx, bwx); + + filter[i+hn].re = bw * kw[i] * filter[i+hn].re; + filter[i+hn].im = bw * kw[i] * filter[i+hn].im; + } + + //circularly shift filter, we shift the filter to left. + ncshift = -abs(ncshift); + circ_shift(filter, nfft, ncshift); + + free_vector_float(kw, -hn, hn); +} + + +float bessi0(float x) +{ + float ax,ans; + double y; + + if ((ax=fabs(x)) < 3.75) { + y=x/3.75; + y*=y; + ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492 + +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2))))); + } else { + y=3.75/ax; + ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1 + +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2 + +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1 + +y*0.392377e-2)))))))); + } + return ans; +} + +#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr +void four1(float data[], unsigned long nn, int isign) +{ + unsigned long n,mmax,m,j,istep,i; + double wtemp,wr,wpr,wpi,wi,theta; + float tempr,tempi; + + n=nn << 1; + j=1; + for (i=1;i i) { + SWAP(data[j],data[i]); + SWAP(data[j+1],data[i+1]); + } + m=nn; + while (m >= 2 && j > m) { + j -= m; + m >>= 1; + } + j += m; + } + mmax=2; + while (n > mmax) { + istep=mmax << 1; + theta=isign*(6.28318530717959/mmax); + wtemp=sin(0.5*theta); + wpr = -2.0*wtemp*wtemp; + wpi=sin(theta); + wr=1.0; + wi=0.0; + for (m=1;m 127) ? 127 : (((A) < 0) ? 0 : A)) /* #define znew (int) (z=36969*(z&65535)+(z>>16)) typedef unsigned long UL; @@ -54,11 +63,12 @@ int assign_sardata_params_ALOSE(struct PRM *, int, int *, int *); void swap_ALOS_data_info(struct sardata_info_ALOSE *sdr); void settable(unsigned long); void print_params(struct PRM *prm); -int check_shift(struct PRM *, int *, int *, int *, int); +int check_shift(struct PRM *, int *, int *, int *, int, int); int set_file_position(FILE *, long *, int); int reset_params(struct PRM *prm, long *, int *, int *); int fill_shift_data(int, int, int, int, int, char *, char *, FILE *); int handle_prf_change_ALOSE(struct PRM *, FILE *, long *, int); +void change_dynamic_range(char *data, long length); struct sardata_record r1; struct sardata_descriptor_ALOSE dfd; @@ -74,11 +84,13 @@ struct sardata_info_ALOSE SARDATA__WCS_ALOSE SARDATA_RVL_ALOSE(SP) */ -long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte_offset) { +long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte_offset, struct resamp_info *rspi, int nPRF) { char *data_fbd, *data, *shift_data; int record_length0; /* length of record read at start of file */ int record_length1; /* length of record read in file */ + int start_sdr_rec_len = 0; /* sdr record length for fisrt record */ + int slant_range_old = 0; /* slant range of previous record */ int line_suffix_size; /* number of bytes after data */ int data_length; /* bytes of data */ int k, n, m, ishift, shift, shift0; @@ -91,6 +103,13 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt if (verbose) fprintf(stderr,".... reading header \n"); + //here we still get sdr from the first data line no matter whether prf changes. + //this sdr is used to initialize record_length0 in assign_sardata_params, which + //is used at line 152 to check if record_length changed. + //I think we should get sdr from first prf-change data line for the output of prf-change file. + //Cunren Liang. 02-DEC-2019 + + /* read header information */ read_sardata_info_ALOSE(imagefile, prm, &header_size, &line_prefix_size); if (verbose) fprintf(stderr,".... reading header %d %d\n", header_size, line_prefix_size); @@ -115,7 +134,7 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt shift0 = 0; n = 1; - m = 0; + m = 2;//first line sequence_number /* read the rest of the file */ while ( (fread((void *) &sdr,sizeof(struct sardata_info_ALOSE), 1, imagefile)) == 1 ) { @@ -124,9 +143,26 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt /* checks for little endian/ big endian */ if (swap) swap_ALOS_data_info(&sdr); + + if (n == 2) + //rspi->frame_counter_start[nPRF] = sdr.frame_counter; + //unfortunately restec format does not have this info, so we are not able to adjust time + rspi->frame_counter_start[nPRF] = 0; + + + /* if this is partway through the file due to prf change, reset sequence, PRF, and near_range */ + if (n == 2) + start_sdr_rec_len = sdr.record_length; + if ((*byte_offset > 0) && (n == 2)) reset_params(prm, byte_offset, &n, &m); + if (sdr.record_length != start_sdr_rec_len) { + printf(" ***** warning sdr.record_length error %d \n", sdr.record_length); + sdr.record_length = start_sdr_rec_len; + sdr.PRF = prm->prf; + sdr.slant_range = slant_range_old; + } if (sdr.sequence_number != n) printf(" missing line: n, seq# %d %d \n", n, sdr.sequence_number); /* check for changes in record_length and PRF */ @@ -136,11 +172,15 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt /* if prf changes, close file and set byte_offset */ if ((sdr.PRF) != prm->prf) { handle_prf_change_ALOSE(prm, imagefile, byte_offset, n); + n-=1; break; } + //rspi->frame_counter_end[nPRF] = sdr.frame_counter; + //unfortunately restec format does not have this info, so we are not able to adjust time + rspi->frame_counter_end[nPRF] = 0; /* check shift to see if it varies from beginning or from command line value */ - check_shift(prm, &shift, &ishift, &shift0, record_length1); + check_shift(prm, &shift, &ishift, &shift0, record_length1, 1); if ((verbose) && (n/2000.0 == n/2000)) { fprintf(stderr," Working on line %d prf %f record length %d slant_range %d \n" @@ -151,6 +191,7 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt if ( fread ((char *) data, record_length1, (size_t) 1, imagefile) != 1 ) break; data_length = record_length1; + slant_range_old = sdr.slant_range; /* write line header to output data */ /* PSA - turning off headers @@ -165,6 +206,7 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt } /* write fbd data */ if (shift == 0) { + change_dynamic_range(data_fbd, data_length/2); fwrite((char *) data_fbd, data_length/2, 1, outfile); } else if (shift != 0) { fill_shift_data(shift, ishift, data_length/2, line_suffix_size, record_length1, data_fbd, shift_data, outfile); @@ -173,18 +215,27 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt else { /* write fbs data */ if (shift == 0) { + change_dynamic_range(data, data_length); fwrite((char *) data, data_length, 1, outfile); } else if (shift != 0) { fill_shift_data(shift, ishift, data_length, line_suffix_size, record_length1, data, shift_data, outfile); } } } - + + //we are not writing out line prefix data, need to correct these parameters + //as they are used in doppler computation. + prm->first_sample = 0; + prm->bytes_per_line -= line_prefix_size; + prm->good_bytes -= line_prefix_size; + + //this is the sdr of the first prf-change data line, should seek back to get last sdr to be used here. /* calculate end time */ prm->SC_clock_stop = get_clock_ALOSE(sdr, tbias); /* m is non-zero only in the event of a prf change */ - prm->num_lines = n - m - 1; + //not correct if PRF changes, so I updated it here. + prm->num_lines = n - m + 1; prm->num_patches = (int)((1.0*n)/(1.0*prm->num_valid_az)); if (prm->num_lines == 0) prm->num_lines = 1; @@ -194,6 +245,15 @@ long read_ALOSE_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byt prm->prf = 1.e3/pri; + prm->xmi = 63.5; + prm->xmq = 63.5; + + rspi->prf[nPRF] = prm->prf; + rspi->SC_clock_start[nPRF] = prm->SC_clock_start; + rspi->num_lines[nPRF] = prm->num_lines; + rspi->num_bins[nPRF] = prm->bytes_per_line/(2*sizeof(char)); + + if (verbose) print_params(prm); free(data); @@ -321,14 +381,14 @@ double get_clock(); /***************************************************************************/ int handle_prf_change_ALOSE(struct PRM *prm, FILE *imagefile, long *byte_offset, int n) { - prm->num_lines = n; + //prm->num_lines = n; fseek(imagefile, -1*sizeof(struct sardata_info_ALOSE), SEEK_CUR); *byte_offset = ftell(imagefile); - printf(" *** PRF changed from %lf to %lf at line %d (byte %ld)\n", (0.001*prm->prf),(0.001*sdr.PRF), n, *byte_offset); - printf(" end: PRF changed from %lf to %lf at line %d \n", (0.001*prm->prf),(0.001*sdr.PRF), n); + printf(" *** PRF changed from %lf to %lf at line %d (byte %ld)\n", (0.001*prm->prf),(0.001*sdr.PRF), n-1, *byte_offset); + // printf(" end: PRF changed from %lf to %lf at line %d \n", (0.001*prm->prf),(0.001*sdr.PRF), n); return(EXIT_SUCCESS); } diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/read_ALOS_data.c b/components/isceobj/Sensor/src/ALOS_pre_process/read_ALOS_data.c index 87d7465..f3f049a 100644 --- a/components/isceobj/Sensor/src/ALOS_pre_process/read_ALOS_data.c +++ b/components/isceobj/Sensor/src/ALOS_pre_process/read_ALOS_data.c @@ -33,6 +33,13 @@ * 07/17/08 reformatted; added functions RJM * ***************************************************************************/ +/******************************************************************************** +This program has been upgraded to handle the ALOS-1 PRF change issue. + +Cunren Liang, 12-December-2019 +California Institute of Technology, Pasadena, CA +*********************************************************************************/ + /* the data header information is read into the structure dfd the line prefix information is read into sdr @@ -47,7 +54,8 @@ SC_clock_start SC_clock_stop #include "image_sio.h" #include "lib_functions.h" - +#define ZERO_VALUE (char)(63 + rand() % 2) +#define clip127(A) (((A) > 127) ? 127 : (((A) < 0) ? 0 : A)) #define znew (int) (z=36969*(z&65535)+(z>>16)) typedef unsigned long UL; static UL z=362436069, t[256]; @@ -61,21 +69,24 @@ void swap_ALOS_data_info(struct sardata_info *sdr); long read_sardata_info(FILE *, struct PRM *, int *, int *); void print_params(struct PRM *prm); int assign_sardata_params(struct PRM *, int, int *, int *); -int check_shift(struct PRM *, int *, int *, int *, int); +int check_shift(struct PRM *, int *, int *, int *, int, int); int set_file_position(FILE *, long *, int); int reset_params(struct PRM *prm, long *, int *, int *); int fill_shift_data(int, int, int, int, int, char *, char *, FILE *); int handle_prf_change(struct PRM *, FILE *, long *, int); +void change_dynamic_range(char *data, long length); struct sardata_record r1; struct sardata_descriptor dfd; struct sardata_info sdr; -long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte_offset) { +long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte_offset, struct resamp_info *rspi, int nPRF) { char *data, *shift_data; int record_length0; /* length of record read at start of file */ int record_length1; /* length of record read in file */ + int start_sdr_rec_len = 0; /* sdr record length for fisrt record */ + int slant_range_old = 0; /* slant range of previous record */ int line_suffix_size; /* number of bytes after data */ int data_length; /* bytes of data */ int n, m, ishift, shift, shift0, npatch_max; @@ -87,6 +98,13 @@ long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte if (debug) fprintf(stderr,".... reading header \n"); + //here we still get sdr from the first data line no matter whether prf changes. + //this sdr is used to initialize record_length0 in assign_sardata_params, which + //is used at line 152 to check if record_length changed. + //I think we should get sdr from first prf-change data line for the output of prf-change file. + //Cunren Liang. 02-DEC-2019 + + /* read header information */ read_sardata_info(imagefile, prm, &header_size, &line_prefix_size); @@ -111,7 +129,7 @@ long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte shift0 = 0; n = 1; - m = 0; + m = 2;//first line sequence_number /* read the rest of the file */ while ( (fread((void *) &sdr,sizeof(struct sardata_info), 1, imagefile)) == 1 ) { @@ -120,10 +138,29 @@ long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte /* checks for little endian/ big endian */ if (swap) swap_ALOS_data_info(&sdr); - /* if this is partway through the file due to prf change, reset sequence, PRF, and near_range */ - if ((*byte_offset > 0) && (n == 2)) reset_params(prm, byte_offset, &n, &m); - if (sdr.sequence_number != n) printf(" missing line: n, seq# %d %d \n", n, sdr.sequence_number); + if (n == 2) + rspi->frame_counter_start[nPRF] = sdr.frame_counter; + + + + /* if this is partway through the file due to prf change, reset sequence, + * PRF, and near_range */ + if (n == 2) + start_sdr_rec_len = sdr.record_length; + + if ((*byte_offset > 0) && (n == 2)) + reset_params(prm, byte_offset, &n, &m); + + if (sdr.record_length != start_sdr_rec_len) { + printf(" ***** warning sdr.record_length error %d \n", sdr.record_length); + sdr.record_length = start_sdr_rec_len; + sdr.PRF = prm->prf; + sdr.slant_range = slant_range_old; + } + if (sdr.sequence_number != n) + printf(" missing line: n, seq# %d %d \n", n, sdr.sequence_number); + /* check for changes in record_length and PRF */ record_length1 = sdr.record_length - line_prefix_size; @@ -132,11 +169,13 @@ long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte /* if prf changes, close file and set byte_offset */ if ((sdr.PRF) != prm->prf) { handle_prf_change(prm, imagefile, byte_offset, n); + n-=1; break; } + rspi->frame_counter_end[nPRF] = sdr.frame_counter; /* check shift to see if it varies from beginning or from command line value */ - check_shift(prm, &shift, &ishift, &shift0, record_length1); + check_shift(prm, &shift, &ishift, &shift0, record_length1, 0); if ((verbose) && (n/2000.0 == n/2000)) { fprintf(stderr," Working on line %d prf %f record length %d slant_range %d \n" @@ -147,27 +186,37 @@ long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte if ( fread ((char *) data, record_length1, (size_t) 1, imagefile) != 1 ) break; data_length = record_length1; + slant_range_old = sdr.slant_range; /* write line header to output data */ - /* PSA turning off headers - fwrite((void *) &sdr, line_prefix_size, 1, outfile); */ + //header is not written to output + //fwrite((void *) &sdr, line_prefix_size, 1, outfile); /* write data */ if (shift == 0) { + change_dynamic_range(data, data_length); fwrite((char *) data, data_length, 1, outfile); /* if data is shifted, fill in with data values of NULL_DATA at start or end*/ } else if (shift != 0) { fill_shift_data(shift, ishift, data_length, line_suffix_size, record_length1, data, shift_data, outfile); } } + + //we are not writing out line prefix data, need to correct these parameters + //as they are used in doppler computation. + prm->first_sample = 0; + prm->bytes_per_line -= line_prefix_size; + prm->good_bytes -= line_prefix_size; /* calculate end time and fix prf */ prm->prf = 0.001*prm->prf; + //this is the sdr of the first prf-change data line, should seek back to get last sdr to be used here. prm->SC_clock_stop = get_clock(sdr, tbias); /* m is non-zero only in the event of a prf change */ - prm->num_lines = n - m - 1; + //not correct if PRF changes, so I updated it here. + prm->num_lines = n - m + 1; /* calculate the maximum number of patches and use that if the default is set to 1000 */ npatch_max = (int)((1.0*n)/(1.0*prm->num_valid_az)); @@ -175,6 +224,15 @@ long read_ALOS_data (FILE *imagefile, FILE *outfile, struct PRM *prm, long *byte if (prm->num_lines == 0) prm->num_lines = 1; + prm->xmi = 63.5; + prm->xmq = 63.5; + + rspi->prf[nPRF] = prm->prf; + rspi->SC_clock_start[nPRF] = prm->SC_clock_start; + rspi->num_lines[nPRF] = prm->num_lines; + rspi->num_bins[nPRF] = prm->bytes_per_line/(2*sizeof(char)); + + if (verbose) print_params(prm); free(data); @@ -293,7 +351,7 @@ double get_clock(); return(EXIT_SUCCESS); } /***************************************************************************/ -int check_shift(struct PRM *prm, int *shift, int *ishift, int *shift0, int record_length1) +int check_shift(struct PRM *prm, int *shift, int *ishift, int *shift0, int record_length1, int ALOS_format) { *shift = 2*floor(0.5 + (sdr.slant_range - prm->near_range)/(0.5*SOL/prm->fs)); *ishift = abs(*shift); @@ -304,7 +362,13 @@ int check_shift(struct PRM *prm, int *shift, int *ishift, int *shift0, int recor } if(*shift != *shift0) { - printf(" near_range, shift = %d %d \n", sdr.slant_range, *shift); + + if(ALOS_format==0) + printf(" near_range, shift = %d %d , at frame_counter: %d, line number: %d\n", sdr.slant_range, *shift, sdr.frame_counter, sdr.sequence_number-1); + if(ALOS_format==1) + printf(" near_range, shift = %d %d\n", sdr.slant_range, *shift); + + *shift0 = *shift; } @@ -324,21 +388,21 @@ int set_file_position(FILE *imagefile, long *byte_offset, int header_size) return(EXIT_SUCCESS); } /***************************************************************************/ -int reset_params(struct PRM *prm, long *byte_offset, int *n, int *m) -{ -double get_clock(); +int reset_params(struct PRM *prm, long *byte_offset, int *n, int *m) { + double get_clock(); prm->SC_clock_start = get_clock(sdr, tbias); prm->prf = sdr.PRF; - prm->near_range = sdr.slant_range; + //comment out so that all data files with different prfs can be aligned at the same starting range + //prm->near_range = sdr.slant_range; *n = sdr.sequence_number; *m = *n; *byte_offset = 0; if (verbose) { - fprintf(stderr," new parameters: \n sequence number %d \n PRF %f\n near_range %lf\n", - *n, 0.001*prm->prf,prm->near_range); - } - return(EXIT_SUCCESS); + fprintf(stderr, " new parameters: \n sequence number %d \n PRF %f\n near_range %lf\n", *n, 0.001 * prm->prf, + prm->near_range); + } + return (EXIT_SUCCESS); } /***************************************************************************/ int fill_shift_data(int shift, int ishift, int data_length, @@ -359,6 +423,7 @@ int k; } /* write the shifted data out */ + change_dynamic_range(shift_data, data_length); fwrite((char *) shift_data, data_length, 1, outfile); return(EXIT_SUCCESS); @@ -366,7 +431,7 @@ int k; /***************************************************************************/ int handle_prf_change(struct PRM *prm, FILE *imagefile, long *byte_offset, int n) { - prm->num_lines = n; + //prm->num_lines = n; /* skip back to beginning of the line */ fseek(imagefile, -1*sizeof(struct sardata_info), SEEK_CUR); @@ -375,9 +440,34 @@ int handle_prf_change(struct PRM *prm, FILE *imagefile, long *byte_offset, int n *byte_offset = ftell(imagefile); /* tell the world */ - printf(" *** PRF changed from %lf to %lf at line %d (byte %ld)\n", (0.001*prm->prf),(0.001*sdr.PRF), n, *byte_offset); - printf(" end: PRF changed from %lf to %lf at line %d \n", (0.001*prm->prf),(0.001*sdr.PRF), n); + printf(" *** PRF changed from %lf to %lf at line %d (byte %ld)\n", (0.001*prm->prf),(0.001*sdr.PRF), n-1, *byte_offset); + // printf(" end: PRF changed from %lf to %lf at line %d \n", (0.001*prm->prf),(0.001*sdr.PRF), n); return(EXIT_SUCCESS); } /***************************************************************************/ + + +void change_dynamic_range(char *data, long length){ + + long i; + + for(i = 0; i < length; i++) + //THIS SHOULD NOT AFFECT DOPPLER COMPUTATION (SUCH AS IN calc_dop.c), BECAUSE + // 1. IQ BIAS IS REMOVED BEFORE COMPUTATION OF DOPPLER. + // 2. 2.0 WILL BE CANCELLED OUT IN atan2f(). + // 3. actual computation results also verified this (even if there is a difference, it is about 0.* Hz) + //data[i] = (unsigned char)clip127(rintf(2. * (data[i] - 15.5) + 63.5)); + data[i] = (unsigned char)clip127(rintf(2.0 * (data[i] - 15.5) + ZERO_VALUE)); + +} + + + + + + + + + + diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/resamp.h b/components/isceobj/Sensor/src/ALOS_pre_process/resamp.h new file mode 100644 index 0000000..7b6858a --- /dev/null +++ b/components/isceobj/Sensor/src/ALOS_pre_process/resamp.h @@ -0,0 +1,106 @@ +////////////////////////////////////// +// Cunren Liang, NASA JPL/Caltech +// Copyright 2017 +////////////////////////////////////// + + +#include +#include +#include +#include + + +#define NR_END 1 +#define FREE_ARG char* +#define PI 3.1415926535897932384626433832795028841971693993751058 + +typedef struct { + float re; + float im; +} fcomplex; + +typedef struct { + double re; + double im; +} dcomplex; + +//allocate arrays +signed char *vector_char(long nl, long nh); +void free_vector_char(signed char *v, long nl, long nh); +unsigned char *vector_unchar(long nl, long nh); +void free_vector_unchar(unsigned char *v, long nl, long nh); +int *vector_int(long nl, long nh); +void free_vector_int(int *v, long nl, long nh); +float *vector_float(long nl, long nh); +void free_vector_float(float *v, long nl, long nh); +double *vector_double(long nl, long nh); +void free_vector_double(double *v, long nl, long nh); +fcomplex *vector_fcomplex(long nl, long nh); +void free_vector_fcomplex(fcomplex *v, long nl, long nh); +signed char **matrix_char(long nrl, long nrh, long ncl, long nch); +void free_matrix_char(signed char **m, long nrl, long nrh, long ncl, long nch); +unsigned char **matrix_unchar(long nrl, long nrh, long ncl, long nch); +void free_matrix_unchar(unsigned char **m, long nrl, long nrh, long ncl, long nch); +float **matrix_float(long nrl, long nrh, long ncl, long nch); +void free_matrix_float(float **m, long nrl, long nrh, long ncl, long nch); +double **matrix_double(long nrl, long nrh, long ncl, long nch); +void free_matrix_double(double **m, long nrl, long nrh, long ncl, long nch); + + +//allocate C-style arrays +FILE **array1d_FILE(long nc); +void free_array1d_FILE(FILE **fv); +signed char *array1d_char(long nc); +void free_array1d_char(signed char *fv); +unsigned char *array1d_unchar(long nc); +void free_array1d_unchar(unsigned char *fv); +int *array1d_int(long nc); +void free_array1d_int(int *fv); +float *array1d_float(long nc); +void free_array1d_float(float *fv); +double *array1d_double(long nc); +void free_array1d_double(double *fv); +fcomplex *array1d_fcomplex(long nc); +void free_array1d_fcomplex(fcomplex *fcv); +dcomplex *array1d_dcomplex(long nc); +void free_array1d_dcomplex(dcomplex *fcv); +signed char **array2d_char(long nl, long nc); +void free_array2d_char(signed char **m); +unsigned char **array2d_unchar(long nl, long nc); +void free_array2d_unchar(unsigned char **m); +float **array2d_float(long nl, long nc); +void free_array2d_float(float **m); +double **array2d_double(long nl, long nc); +void free_array2d_double(double **m); +fcomplex **array2d_fcomplex(long nl, long nc); +void free_array2d_fcomplex(fcomplex **m); + +//handling error +void nrerror(char error_text[]); + +//complex operations +fcomplex cmul(fcomplex a, fcomplex b); +fcomplex cconj(fcomplex z); +fcomplex cadd(fcomplex a, fcomplex b); +float xcabs(fcomplex z); +float cphs(fcomplex z); + +//functions +long next_pow2(long a); +void circ_shift(fcomplex *in, int na, int nc); +void left_shift(fcomplex *in, int na); +void right_shift(fcomplex *in, int na); +int roundfi(float a); +void sinc(int n, int m, float *coef); +void kaiser(int n, int m, float *coef, float beta); +void kaiser2(float beta, int n, float *coef); +void bandpass_filter(float bw, float bc, int n, int nfft, int ncshift, float beta, fcomplex *filter); +float bessi0(float x); +void four1(float data[], unsigned long nn, int isign); + +//file operations +FILE *openfile(char *filename, char *pattern); +void readdata(void *data, size_t blocksize, FILE *fp); +void writedata(void *data, size_t blocksize, FILE *fp); +long file_length(FILE* fp, long width, long element_size); + diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/resamp_azimuth.c b/components/isceobj/Sensor/src/ALOS_pre_process/resamp_azimuth.c new file mode 100644 index 0000000..3e1fd61 --- /dev/null +++ b/components/isceobj/Sensor/src/ALOS_pre_process/resamp_azimuth.c @@ -0,0 +1,246 @@ +////////////////////////////////////// +// Cunren Liang +// California Institute of Technology +// Copyright 2019 +////////////////////////////////////// + +//this program is tested against resamp.c, the outputs of the two are exactly the same. + +#include "resamp.h" + +//ALOS I or Q mean = 15.5, so get 15 or 16 randomly here +//#define ZERO_VALUE (char)(15 + rand() % 2) +//I changed the dynamic range when reading data +//ALOS I or Q mean = 63.5, so get 63 or 64 randomly here +#define ZERO_VALUE (char)(63 + rand() % 2) +typedef struct { + char re; + char im; +} char_complex; +char_complex *array1d_char_complex(long nc); +void free_array1d_char_complex(char_complex *fcv); +void normalize_kernel(float *kernel, long start_index, long end_index); +int resamp_azimuth(char *slc2, char *rslc2, int nrg, int naz1, int naz2, double prf, double *dopcoeff, double *azcoef, int n, double beta){ + int i; + int verbose = 0; + if(verbose){ + printf("\n\ninput parameters:\n"); + printf("slc2: %s\n", slc2); + printf("rslc2: %s\n", rslc2); + printf("nrg: %d\n", nrg); + printf("naz1: %d\n", naz1); + printf("naz2: %d\n\n", naz2); + printf("prf: %f\n\n", prf); + for(i = 0; i < 4; i++){ + printf("dopcoeff[%d]: %e\n", i, dopcoeff[i]); + } + printf("\n"); + for(i = 0; i < 2; i++){ + printf("azcoef[%d]: %e\n", i, azcoef[i]); + } + printf("\n"); + } + FILE *slc2fp; + FILE *rslc2fp; + int m; + int interp_method; + int edge_method; + float azpos; + float azoff; + float az2; + int azi2; + float azf; + int azfn; + int hnm; + int hn; + float *sincc; + float *kaiserc; + float *kernel; + float *azkernel; + fcomplex *azkernel_fc; + fcomplex *rgrs; + fcomplex *azca; + fcomplex *rgrsb; + fcomplex *azrs; + char_complex *inl; + char_complex *outl; + float *dop; + float dopx; + fcomplex **inb; + int j, k, k1, k2; + int tmp1, tmp2; + int zero_flag; + float ftmp1, ftmp2; + fcomplex fctmp1, fctmp2; + m = 10000; + interp_method = 0; + edge_method = 2; + if((n % 2 == 0) || (n < 3)){ + fprintf(stderr, "number of samples to be used in the resampling must be odd, and larger or equal to than 3\n"); + exit(1); + } + slc2fp = openfile(slc2, "rb"); + rslc2fp = openfile(rslc2, "wb"); + hn = n / 2; + hnm = n * m / 2; + sincc = vector_float(-hnm, hnm); + kaiserc = vector_float(-hnm, hnm); + kernel = vector_float(-hnm, hnm); + azkernel = vector_float(-hn, hn); + azkernel_fc = vector_fcomplex(-hn, hn); + rgrs = vector_fcomplex(-hn, hn); + azca = vector_fcomplex(-hn, hn); + rgrsb = vector_fcomplex(-hn, hn); + azrs = array1d_fcomplex(nrg); + inl = array1d_char_complex(nrg); + outl = array1d_char_complex(nrg); + dop = array1d_float(nrg); + inb = array2d_fcomplex(naz2, nrg); + sinc(n, m, sincc); + kaiser(n, m, kaiserc, beta); + for(i = -hnm; i <= hnm; i++) + kernel[i] = kaiserc[i] * sincc[i]; + for(i = 0; i < nrg; i++){ + dop[i] = dopcoeff[0] + dopcoeff[1] * i + dopcoeff[2] * i * i + dopcoeff[3] * i * i * i; + if(verbose){ + if(i % 500 == 0) + printf("range sample: %5d, doppler centroid frequency: %8.2f Hz\n", i, dop[i]); + } + } + for(i = 0; i < naz2; i++){ + readdata((char_complex *)inl, (size_t)nrg * sizeof(char_complex), slc2fp); + for(j =0; j < nrg; j++){ + inb[i][j].re = inl[j].re; + inb[i][j].im = inl[j].im; + } + } + for(i = 0; i < naz1; i++){ + if((i + 1) % 100 == 0) + fprintf(stderr,"processing line: %6d of %6d\r", i+1, naz1); + for(j = 0; j < nrg; j++){ + azrs[j].re = 0.0; + azrs[j].im = 0.0; + } + azpos = i; + azoff = azcoef[0] + azpos * azcoef[1]; + az2 = i + azoff; + azi2 = roundfi(az2); + azf = az2 - azi2; + azfn = roundfi(azf * m); + if(edge_method == 0){ + if(azi2 < hn || azi2 > naz2 - 1 - hn){ + for(j = 0; j < nrg; j++){ + outl[j].re = ZERO_VALUE; + outl[j].im = ZERO_VALUE; + } + writedata((char_complex *)outl, (size_t)nrg * sizeof(char_complex), rslc2fp); + continue; + } + } + else if(edge_method == 1){ + if(azi2 < 0 || azi2 > naz2 - 1){ + for(j = 0; j < nrg; j++){ + outl[j].re = ZERO_VALUE; + outl[j].im = ZERO_VALUE; + } + writedata((char_complex *)outl, (size_t)nrg * sizeof(char_complex), rslc2fp); + continue; + } + } + else{ + if(azi2 < -hn || azi2 > naz2 - 1 + hn){ + for(j = 0; j < nrg; j++){ + outl[j].re = ZERO_VALUE; + outl[j].im = ZERO_VALUE; + } + writedata((char_complex *)outl, (size_t)nrg * sizeof(char_complex), rslc2fp); + continue; + } + } + for(k = -hn; k <= hn; k++){ + tmp2 = k * m - azfn; + if(tmp2 > hnm) tmp2 = hnm; + if(tmp2 < -hnm) tmp2 = -hnm; + azkernel[k] = kernel[tmp2]; + } + normalize_kernel(azkernel, -hn, hn); + for(j = 0; j < nrg; j++){ + for(k1 = -hn; k1 <= hn; k1++){ + if((azi2 + k1 >= 0)&&(azi2 + k1 <= naz2-1)){ + rgrs[k1].re = inb[azi2 + k1][j].re; + rgrs[k1].im = inb[azi2 + k1][j].im; + } + else{ + rgrs[k1].re = ZERO_VALUE; + rgrs[k1].im = ZERO_VALUE; + } + } + dopx = dop[j]; + for(k = -hn; k <= hn; k++){ + ftmp1 = 2.0 * PI * dopx * k / prf; + azca[k].re = cos(ftmp1); + azca[k].im = sin(ftmp1); + if(interp_method == 0){ + rgrsb[k] = cmul(rgrs[k], cconj(azca[k])); + azrs[j].re += rgrsb[k].re * azkernel[k]; + azrs[j].im += rgrsb[k].im * azkernel[k]; + } + else{ + azkernel_fc[k].re = azca[k].re * azkernel[k]; + azkernel_fc[k].im = azca[k].im * azkernel[k]; + azrs[j] = cadd(azrs[j], cmul(rgrs[k], azkernel_fc[k])); + } + } + if(interp_method == 0){ + ftmp1 = 2.0 * PI * dopx * azf / prf; + fctmp1.re = cos(ftmp1); + fctmp1.im = sin(ftmp1); + azrs[j] = cmul(azrs[j], fctmp1); + } + } + for(j = 0; j < nrg; j++){ + outl[j].re = roundfi(azrs[j].re); + outl[j].im = roundfi(azrs[j].im); + } + writedata((char_complex *)outl, (size_t)nrg * sizeof(char_complex), rslc2fp); + } + fprintf(stderr,"processing line: %6d of %6d\n", naz1, naz1); + free_vector_float(sincc, -hnm, hnm); + free_vector_float(kaiserc, -hnm, hnm); + free_vector_float(kernel, -hnm, hnm); + free_vector_float(azkernel, -hn, hn); + free_vector_fcomplex(azkernel_fc, -hn, hn); + free_vector_fcomplex(rgrs, -hn, hn); + free_vector_fcomplex(azca, -hn, hn); + free_vector_fcomplex(rgrsb, -hn, hn); + free_array1d_fcomplex(azrs); + free_array1d_char_complex(inl); + free_array1d_char_complex(outl); + free_array1d_float(dop); + free_array2d_fcomplex(inb); + fclose(slc2fp); + fclose(rslc2fp); + return 0; +} +char_complex *array1d_char_complex(long nc){ + char_complex *fcv; + fcv = (char_complex*) malloc(nc * sizeof(char_complex)); + if(!fcv){ + fprintf(stderr,"Error: cannot allocate 1-D char complex array\n"); + exit(1); + } + return fcv; +} +void free_array1d_char_complex(char_complex *fcv){ + free(fcv); +} +void normalize_kernel(float *kernel, long start_index, long end_index){ + double sum; + long i; + sum = 0.0; + for(i = start_index; i <= end_index; i++) + sum += kernel[i]; + if(sum!=0) + for(i = start_index; i <= end_index; i++) + kernel[i] /= sum; +} diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.c b/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.c index 47b6e7c..8bffd80 100644 --- a/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.c +++ b/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.c @@ -1,48 +1,51 @@ +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 + #include "image_sio.h" #include "siocomplex.h" #include -fcomplex Cmul(fcomplex x, fcomplex y) +fcomplex_sio Cmul(fcomplex_sio x, fcomplex_sio y) { - fcomplex z; + fcomplex_sio z; z.r = x.r*y.r - x.i*y.i; z.i = x.i*y.r + x.r*y.i; return z; } -fcomplex Cexp(float theta) +fcomplex_sio Cexp(float theta) { - fcomplex z; + fcomplex_sio z; z.r = cos(theta); z.i = sin(theta); return z; } -fcomplex Conjg(fcomplex z) +fcomplex_sio Conjg(fcomplex_sio z) { - fcomplex x; + fcomplex_sio x; x.r = z.r; x.i = -z.i; return x; } -fcomplex RCmul(float a, fcomplex z) +fcomplex_sio RCmul(float a, fcomplex_sio z) { - fcomplex x; + fcomplex_sio x; x.r = a*z.r; x.i = a*z.i; return x; } -fcomplex Cadd(fcomplex x, fcomplex y) +fcomplex_sio Cadd(fcomplex_sio x, fcomplex_sio y) { - fcomplex z; + fcomplex_sio z; z.r = x.r + y.r; z.i = x.i + y.i; return z; } -float Cabs(fcomplex z) +float Cabs(fcomplex_sio z) { return hypot(z.r, z.i); } diff --git a/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.h b/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.h index 63b5e0c..8600c3b 100644 --- a/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.h +++ b/components/isceobj/Sensor/src/ALOS_pre_process/siocomplex.h @@ -1,11 +1,14 @@ +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 + #ifndef _COMPLEX_H #define _COMPLEX_H -fcomplex Cmul(fcomplex x, fcomplex y); -fcomplex Cexp(float theta); -fcomplex Conjg(fcomplex z); -fcomplex RCmul(float a, fcomplex z); -fcomplex Cadd(fcomplex x, fcomplex y); -float Cabs(fcomplex z); +fcomplex_sio Cmul(fcomplex_sio x, fcomplex_sio y); +fcomplex_sio Cexp(float theta); +fcomplex_sio Conjg(fcomplex_sio z); +fcomplex_sio RCmul(float a, fcomplex_sio z); +fcomplex_sio Cadd(fcomplex_sio x, fcomplex_sio y); +float Cabs(fcomplex_sio z); #endif /* _COMPLEX_H */ diff --git a/components/stdproc/alosreformat/ALOS_fbd2fbs/bindings/ALOS_fbd2fbsmodule.c b/components/stdproc/alosreformat/ALOS_fbd2fbs/bindings/ALOS_fbd2fbsmodule.c index 6dbd139..58fc7e2 100644 --- a/components/stdproc/alosreformat/ALOS_fbd2fbs/bindings/ALOS_fbd2fbsmodule.c +++ b/components/stdproc/alosreformat/ALOS_fbd2fbs/bindings/ALOS_fbd2fbsmodule.c @@ -25,7 +25,8 @@ // Author: Giangi Sacco //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 #include @@ -219,8 +220,8 @@ PyObject * ALOS_fbd2fbs_C(PyObject* self, PyObject* args) i = j + r.first_sample; /* increase dynamic range by 2 and set the mean value to 63.5 */ - rtest = rintf(2.*cout[j].r+63.5); - itest = rintf(2.*cout[j].i+63.5); + rtest = rintf(cout[j].r+r.xmi); + itest = rintf(cout[j].i+r.xmq); /* sometimes the range can exceed 0-127 so clip the numbers to be in the correct range */ diff --git a/components/stdproc/alosreformat/ALOS_fbs2fbd/bindings/ALOS_fbs2fbdmodule.c b/components/stdproc/alosreformat/ALOS_fbs2fbd/bindings/ALOS_fbs2fbdmodule.c index d035c97..a6fd4fa 100644 --- a/components/stdproc/alosreformat/ALOS_fbs2fbd/bindings/ALOS_fbs2fbdmodule.c +++ b/components/stdproc/alosreformat/ALOS_fbs2fbd/bindings/ALOS_fbs2fbdmodule.c @@ -25,7 +25,8 @@ // Author: Giangi Sacco //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - +//program updated to handle PRF change issue of ALOS-1 +//Cunren Liang, December 2019 #include @@ -197,10 +198,10 @@ PyObject * ALOS_fbs2fbd_C(PyObject* self, PyObject* args) n4 = nffti/4; for(i=0; i