pktools 2.6.7
Processing Kernel for geospatial data
pkfilterascii.cc
1/**********************************************************************
2pkfilterascii.cc: program to filter data in an ASCII file
3Copyright (C) 2008-2014 Pieter Kempeneers
4
5This file is part of pktools
6
7pktools is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published by
9the Free Software Foundation, either version 3 of the License, or
10(at your option) any later version.
11
12pktools is distributed in the hope that it will be useful,
13but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15GNU General Public License for more details.
16
17You should have received a copy of the GNU General Public License
18along with pktools. If not, see <http://www.gnu.org/licenses/>.
19***********************************************************************/
20#include <assert.h>
21#include <iostream>
22#include <string>
23#include <fstream>
24#include <math.h>
25#include <sys/types.h>
26#include <stdio.h>
27#include "base/Optionpk.h"
28#include "base/Vector2d.h"
29#include "algorithms/Filter.h"
30#include "fileclasses/FileReaderAscii.h"
31
32extern "C" {
33#include <gsl/gsl_sort.h>
34}
35
36/******************************************************************************/
81using namespace std;
82
83/*------------------
84 Main procedure
85 ----------------*/
86int main(int argc,char **argv) {
87 Optionpk<std::string> input_opt("i","input","input ASCII file");
88 Optionpk<std::string> output_opt("o", "output", "Output ASCII file");
89 Optionpk<int> inputCols_opt("ic", "inputCols", "input columns (e.g., for three dimensional input data in first three columns use: -ic 0 -ic 1 -ic 2");
90 Optionpk<std::string> method_opt("f", "filter", "filter function (to be implemented: dwt, dwti,dwt_cut)");
91 Optionpk<std::string> wavelet_type_opt("wt", "wavelet", "wavelet type: daubechies,daubechies_centered, haar, haar_centered, bspline, bspline_centered", "daubechies");
92 Optionpk<int> family_opt("wf", "family", "wavelet family (vanishing moment, see also http://www.gnu.org/software/gsl/manual/html_node/DWT-Initialization.html)", 4);
93 Optionpk<double> threshold_opt("cut", "cut", "threshold to cut dwt coefficients. Use 0 to keep all.", 0);
94 Optionpk<int> dimZ_opt("dz", "dz", "filter kernel size in z (band or spectral dimension), must be odd (example: 3).. Set dz>0 if 1-D filter must be used in band domain");
95 Optionpk<double> tapZ_opt("tapz", "tapz", "taps used for spectral filtering");
96 Optionpk<double> fwhm_opt("fwhm", "fwhm", "list of full width half to apply spectral filtering (-fwhm band1 -fwhm band2 ...)");
97 Optionpk<std::string> srf_opt("srf", "srf", "list of ASCII files containing spectral response functions (two columns: wavelength response)");
98 Optionpk<int> wavelengthIn_opt("win", "wavelengthIn", "column number of input ASCII file containing wavelengths");
99 Optionpk<double> wavelengthOut_opt("wout", "wavelengthOut", "list of wavelengths in output spectrum (-wout band1 -wout band2 ...)");
100 Optionpk<std::string> interpolationType_opt("interp", "interp", "type of interpolation for spectral filtering (see http://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html)","akima");
101 Optionpk<bool> transpose_opt("t", "transpose", "transpose output with samples in rows and wavelengths in cols", false);
102 Optionpk<short> verbose_opt("v", "verbose", "verbose mode if > 0", 0,2);
103
104 tapZ_opt.setHide(1);
105 fwhm_opt.setHide(1);
106 srf_opt.setHide(1);
107 wavelengthIn_opt.setHide(1);
108 wavelengthOut_opt.setHide(1);
109 interpolationType_opt.setHide(1);
110 transpose_opt.setHide(1);
111 wavelet_type_opt.setHide(1);
112 family_opt.setHide(1);
113 threshold_opt.setHide(1);
114
115 bool doProcess;//stop process when program was invoked with help option (-h --help)
116 try{
117 doProcess=input_opt.retrieveOption(argc,argv);
118 output_opt.retrieveOption(argc,argv);
119 inputCols_opt.retrieveOption(argc,argv);
120 method_opt.retrieveOption(argc,argv);
121 dimZ_opt.retrieveOption(argc,argv);
122 tapZ_opt.retrieveOption(argc,argv);
123 fwhm_opt.retrieveOption(argc,argv);
124 srf_opt.retrieveOption(argc,argv);
125 wavelengthIn_opt.retrieveOption(argc,argv);
126 wavelengthOut_opt.retrieveOption(argc,argv);
127 interpolationType_opt.retrieveOption(argc,argv);
128 transpose_opt.retrieveOption(argc,argv);
129 wavelet_type_opt.retrieveOption(argc,argv);
130 family_opt.retrieveOption(argc,argv);
131 threshold_opt.retrieveOption(argc,argv);
132 verbose_opt.retrieveOption(argc,argv);
133 }
134 catch(string predefinedString){
135 std::cout << predefinedString << std::endl;
136 exit(0);
137 }
138 if(!doProcess){
139 cout << endl;
140 cout << "Usage: pkfilterascii -i input.txt [-ic column]*" << endl;
141 cout << endl;
142 std::cout << "short option -h shows basic options only, use long option --help to show all options" << std::endl;
143 exit(0);//help was invoked, stop processing
144 }
145
146 Vector2d<double> inputData(inputCols_opt.size());
147 Vector2d<double> filteredData(inputCols_opt.size());
148 std::vector<double> wavelengthIn;
149 std::vector<double> wavelengthOut;
150 assert(input_opt.size());
151 FileReaderAscii asciiReader(input_opt[0]);
152 if(wavelengthIn_opt.size())
153 asciiReader.readData(wavelengthIn,wavelengthIn_opt[0]);
154 assert(inputCols_opt.size());
155 asciiReader.readData(inputData,inputCols_opt);
156 if(verbose_opt[0]){
157 std::cout << "wavelengthIn.size(): " << wavelengthIn.size() << std::endl;
158 std::cout << "inputData[0].size(): " << inputData[0].size() << std::endl;
159 }
160 if(wavelengthIn.size())
161 assert(wavelengthIn.size()==inputData[0].size());
162 asciiReader.close();
163 filter::Filter filter1d;
164 if(fwhm_opt.size()){
165 filteredData.resize(inputCols_opt.size(),wavelengthOut_opt.size());
166 assert(wavelengthIn_opt.size());
167 if(verbose_opt[0])
168 std::cout << "spectral filtering to " << fwhm_opt.size() << " bands with provided fwhm " << std::endl;
169 assert(wavelengthOut_opt.size()==fwhm_opt.size());
170 std::vector<double> fwhmData(wavelengthOut_opt.size());
171 for(int icol=0;icol<inputCols_opt.size();++icol)
172 filter1d.applyFwhm<double>(wavelengthIn,inputData[icol], wavelengthOut_opt,fwhm_opt, interpolationType_opt[0], filteredData[icol],verbose_opt[0]);
173 if(verbose_opt[0])
174 std::cout << "spectra filtered to " << wavelengthOut_opt.size() << " bands" << std::endl;
175 wavelengthOut=wavelengthOut_opt;
176 }
177 else if(srf_opt.size()){
178 wavelengthOut.resize(srf_opt.size());
179 filteredData.resize(inputCols_opt.size(),srf_opt.size());
180 Vector2d<double> srfData(srf_opt.size(),inputCols_opt.size());//transposed output
181 if(verbose_opt[0])
182 std::cout << "spectral filtering to " << srf_opt.size() << " bands with provided SRF " << std::endl;
183 std::vector< Vector2d<double> > srf(srf_opt.size());//[0] srf_nr, [1]: wavelength, [2]: response
184 ifstream srfFile;
185 for(int isrf=0;isrf<srf_opt.size();++isrf){
186 srf[isrf].resize(2);
187 srfFile.open(srf_opt[isrf].c_str());
188 double v;
189 //add 0 to make sure srf is 0 at boundaries after interpolation step
190 srf[isrf][0].push_back(0);
191 srf[isrf][1].push_back(0);
192 srf[isrf][0].push_back(1);
193 srf[isrf][1].push_back(0);
194 while(srfFile >> v){
195 srf[isrf][0].push_back(v);
196 srfFile >> v;
197 srf[isrf][1].push_back(v);
198 }
199 srfFile.close();
200 //add 0 to make sure srf[isrf] is 0 at boundaries after interpolation step
201 srf[isrf][0].push_back(srf[isrf][0].back()+1);
202 srf[isrf][1].push_back(0);
203 srf[isrf][0].push_back(srf[isrf][0].back()+1);
204 srf[isrf][1].push_back(0);
205 if(verbose_opt[0])
206 cout << "srf file details: " << srf[isrf][0].size() << " wavelengths defined" << endl;
207 if(verbose_opt[0]>1){
208 for(int iw=0;iw<srf[isrf][0].size();++iw)
209 std::cout << srf[isrf][0][iw] << " " << srf[isrf][1][iw] << std::endl;
210 }
211 }
212 double centreWavelength=0;
213 for(int icol=0;icol<inputCols_opt.size();++icol)
214 filteredData[icol].resize(srf.size());
215
216 for(int isrf=0;isrf<srf.size();++isrf){
217 double delta=1.0;
218 bool normalize=true;
219 centreWavelength=filter1d.applySrf<double>(wavelengthIn,inputData,srf[isrf], interpolationType_opt[0], srfData[isrf], delta, normalize,1,true, verbose_opt[0]);
220 if(verbose_opt[0])
221 std::cout << "centre wavelength srf " << isrf << ": " << centreWavelength << std::endl;
222 wavelengthOut[isrf]=static_cast<int>(centreWavelength+0.5);
223 }
224 srfData.transpose(filteredData);
225 if(verbose_opt[0])
226 std::cout << "spectra filtered to " << srf.size() << " bands" << std::endl;
227 }
228 else{//no filtering
229 if(verbose_opt[0])
230 std::cout << "no filtering selected" << std::endl;
231 for(int icol=0;icol<inputCols_opt.size();++icol)
232 filteredData[icol]=inputData[icol];
233 }
234
235 if(method_opt.size()){
236 wavelengthOut=wavelengthIn;
237 for(int icol=0;icol<inputCols_opt.size();++icol){
238 switch(filter::Filter::getFilterType(method_opt[0])){
239 case(filter::dwt):
240 filter1d.dwtForward(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
241 break;
242 case(filter::dwti):
243 filter1d.dwtInverse(filteredData[icol],wavelet_type_opt[0],family_opt[0]);
244 break;
245 case(filter::dwt_cut):
246 filter1d.dwtCut(filteredData[icol],wavelet_type_opt[0],family_opt[0],threshold_opt[0]);
247 break;
248 case(filter::smooth):
249 if(tapZ_opt.size()){
250 filter1d.setTaps(tapZ_opt);
251 filter1d.filter(inputData[icol],filteredData[icol]);
252 }
253 else{
254 assert(dimZ_opt.size());
255 filter1d.smooth(inputData[icol],filteredData[icol],dimZ_opt[0]);
256 }
257 break;
258 default:
259 assert(tapZ_opt.size());
260 filter1d.filter(inputData[icol],filteredData[icol]);
261 // if(verbose_opt[0])
262 // std::cout << "method to be implemented" << std::endl;
263 // exit(1);
264 break;
265 }
266 }
267 }
268 ofstream outputStream;
269 if(!output_opt.empty())
270 outputStream.open(output_opt[0].c_str(),ios::out);
271
272 if(verbose_opt[0])
273 std::cout << "stream to output" << std::endl;
274 if(transpose_opt[0]){
275 for(int icol=0;icol<inputCols_opt.size();++icol){
276 for(int iband=0;iband<filteredData[icol].size();++iband){
277 if(!output_opt.empty()){
278 outputStream << filteredData[icol][iband];
279 if(iband<filteredData[icol].size()-1)
280 outputStream << " ";
281 else
282 outputStream << std::endl;
283 }
284 else{
285 std::cout << filteredData[icol][iband];
286 if(iband<filteredData[icol].size()-1)
287 std::cout << " ";
288 else
289 std::cout << std::endl;
290 }
291 }
292 }
293 }
294 else{
295 // int nband=wavelengthOut.size()? wavelengthOut.size() : filteredData[0].size();
296 int nband=0;
297 if(method_opt.size()){
298 switch(filter::Filter::getFilterType(method_opt[0])){
299 case(filter::dwt):
300 nband=filteredData[0].size();
301 break;
302 case(filter::dwti):
303 nband=filteredData[0].size();
304 break;
305 case(filter::dwt_cut):
306 nband=filteredData[0].size();
307 break;
308 default:
309 nband=wavelengthOut.size();
310 break;
311 }
312 }
313 else
314 nband=wavelengthOut.size();
315 if(verbose_opt[0]){
316 std::cout << "number of bands: " << nband << std::endl;
317 std::cout << "wavelengthOut.size(): " << wavelengthOut.size() << std::endl;
318 }
319 for(int iband=0;iband<nband;++iband){
320 if(!output_opt.empty()){
321 if(wavelengthOut.size())
322 outputStream << wavelengthOut[iband] << " ";
323 else if(wavelengthIn_opt.size())
324 outputStream << iband << " ";
325 }
326 else{
327 if(wavelengthOut.size())
328 std::cout << wavelengthOut[iband] << " ";
329 else if(wavelengthIn_opt.size())
330 std::cout << iband << " ";
331 }
332 for(int icol=0;icol<inputCols_opt.size();++icol){
333 if(!output_opt.empty()){
334 outputStream << filteredData[icol][iband];
335 if(icol<inputCols_opt.size()-1)
336 outputStream << " ";
337 else
338 outputStream << std::endl;
339 }
340 else{
341 std::cout << filteredData[icol][iband];
342 if(icol<inputCols_opt.size()-1)
343 std::cout << " ";
344 else
345 std::cout << std::endl;
346 }
347 }
348 }
349 }
350 if(!output_opt.empty())
351 outputStream.close();
352 return 0;
353}