UltraScan III
us_spectrodata.cpp
Go to the documentation of this file.
1 
3 #include "us_spectrodata.h"
4 #include "us_defines.h"
5 #include "us_settings.h"
6 
7 #define LO_DTERM 0.2500 // low decay-term point (1/4)
8 
9 // Provides raster data for QwtSpectrogram
10 US_SpectrogramData::US_SpectrogramData() : QwtRasterData()
11 {
12  rdata.clear();
13  xreso = 300.0;
14  yreso = 300.0;
15  resol = 90.0;
16  zfloor = 0.20;
17  nxpsc = qRound( xreso );
18  nyscn = qRound( yreso );
19  nxypt = nxpsc * nyscn;
20  dbg_level = US_Settings::us_debug();
21 }
22 
23 QwtRasterData* US_SpectrogramData::copy() const
24 {
25  return new US_SpectrogramData();
26 }
27 
28 QwtDoubleInterval US_SpectrogramData::range() const
29 {
30  return QwtDoubleInterval( zmin, zmax );
31 }
32 
33 // Initialize raster: get x,y ranges and image pixel size
34 void US_SpectrogramData::initRaster( QwtDoubleRect& drect, QSize& rsiz )
35 {
36  drect = QwtDoubleRect( xmin, ymin, xrng, yrng );
37  rsiz = QSize( nxpsc, nyscn );
38 }
39 
40 // Get x range
41 QwtDoubleInterval US_SpectrogramData::xrange()
42 {
43  return QwtDoubleInterval( xmin, xmax );
44 }
45 
46 // Get y range
47 QwtDoubleInterval US_SpectrogramData::yrange()
48 {
49  return QwtDoubleInterval( ymin, ymax );
50 }
51 
52 // Necessary method to set up raster ranges: resolutions, floor fraction
53 void US_SpectrogramData::setRastRanges( double a_xres, double a_yres,
54  double a_reso, double a_zfloor, QwtDoubleRect a_drecti )
55 {
56  xreso = a_xres;
57  yreso = a_yres;
58  resol = a_reso;
59  zfloor = a_zfloor;
60  zfloor = ( zfloor < 1.0 ) ? zfloor : ( zfloor / 100.0 );
61  nxpsc = qRound( xreso );
62  nyscn = qRound( yreso );
63  nxypt = nxpsc * nyscn;
64  drecti = a_drecti;
65 }
66 
67 // Set constant Z range for use with manual or looping displays
68 void US_SpectrogramData::setZRange( double a_zmin, double a_zmax )
69 {
70  zmin = a_zmin;
71  zmax = a_zmax;
72 }
73 
74 // Method called by QwtPlotSpectrogram for each raster point.
75 // This version gets or interpolates a point from the raster buffer.
76 double US_SpectrogramData::value(double x, double y) const
77 {
78  double rx = ( x - xmin ) * xinc; // real x pixel position
79  double ry = ( ymax - y ) * yinc; // real y pixel position
80  int jx = (int)rx; // integral x pixel
81  int jy = (int)ry; // integral y pixel
82  int jr1 = jy * nxpsc + jx; // overall raster position
83  int mxr = nxypt - 1; // max raster position index
84 
85  if ( jr1 < 0 || jr1 > mxr )
86  { // for x,y outside raster, return z-minimum
87  return zmin;
88  }
89 
90  double rv1 = rdata.at( jr1 ); // possible output value
91 
92  double dx = rx - (double)jx; // fractional part of x
93  double dy = ry - (double)jy; // fractional part of y
94 
95  if ( dx > 0.1 || dy > 0.1 )
96  { // if either fraction significant, interpolate output value
97  int jr2 = jr1 + 1; // x+1,y
98  int jr3 = jr1 + nxpsc; // x,y+1
99  int jr4 = jr3 + 1; // x+1,y+1
100  jr2 = ( jr2 > mxr ) ? jr1 : jr2;
101  jr3 = ( jr3 > mxr ) ? jr2 : jr3;
102  jr4 = ( jr4 > mxr ) ? jr3 : jr4;
103  double rv2 = rdata.at( jr2 ); // surrounding input values
104  double rv3 = rdata.at( jr3 );
105  double rv4 = rdata.at( jr4 );
106 
107  // interpolate between the four points
108  rv2 = ( dx * rv2 ) + ( ( 1.0 - dx ) * rv1 );
109  rv4 = ( dx * rv4 ) + ( ( 1.0 - dx ) * rv3 );
110  rv1 = ( dy * rv4 ) + ( ( 1.0 - dy ) * rv2 );
111  }
112 
113  return rv1; // return raster value or interpolated value
114 }
115 
116 // Set up raster values from solution distribution set.
117 void US_SpectrogramData::setRaster( QList< S_Solute >* solu )
118 {
119  double xval;
120  double yval;
121  double zval;
122 
123  int nsol = solu->size(); // number Solute points
124  if ( nsol < 1 )
125  return;
126 
127  xmin = drecti.left ();
128  xmax = drecti.right();
129 
130  if ( xmin == xmax )
131  { // Auto-limits: calculate x,y ranges
132  xmin = solu->at( 0 ).s; // initial minima,maxima
133  ymin = solu->at( 0 ).k;
134  zmin = solu->at( 0 ).c;
135  xmax = xmin;
136  ymax = ymin;
137  zmax = zmin;
138 
139  // scan solute distribution for ranges
140 
141  for ( int ii = 1; ii < nsol; ii++ )
142  {
143  xval = solu->at( ii ).s;
144  yval = solu->at( ii ).k;
145  zval = solu->at( ii ).c;
146  xmin = qMin( xval, xmin );
147  xmax = qMax( xval, xmax );
148  ymin = qMin( yval, ymin );
149  ymax = qMax( yval, ymax );
150  zmin = qMin( zval, zmin );
151  zmax = qMax( zval, zmax );
152  }
153 
154  xrng = xmax - xmin; // initial ranges and pixel/data ratios
155  xinc = ( xreso - 1.0 ) / xrng;
156  yrng = ymax - ymin;
157  if ( yrng == 0.0 )
158  {
159  yrng = 0.02;
160  ymin = (double)( qFloor( ymin * 100.0 ) ) * 0.01;
161  ymax = ymin + yrng;
162  }
163  yinc = ( yreso - 1.0 ) / yrng;
164 
165  xmin -= ( 4.0 / xinc ); // adjust for padding, then recalc ranges
166  xmax += ( 4.0 / xinc );
167  ymin -= ( 4.0 / yinc );
168  ymax += ( 4.0 / yinc );
169  }
170 
171  else
172  { // Given x,y ranges
173  ymin = drecti.top ();
174  ymax = drecti.bottom();
175  }
176 
177  if ( nsol == 1 )
178  {
179  zmin *= zfloor;
180  xmin *= 0.90;
181  xmax *= 1.10;
182  ymin *= 0.90;
183  ymax *= 1.10;
184  }
185 
186  xrng = xmax - xmin;
187  yrng = ymax - ymin;
188  xinc = ( xreso - 1.0 ) / xrng;
189  yinc = ( yreso - 1.0 ) / yrng;
190  zminr = zmin - ( ( zmax - zmin ) * zfloor );
191 
192  // set bounding rectangle for raster plot
193  setBoundingRect( QwtDoubleRect( xmin, ymin, xrng, yrng ) );
194 
195  // initialize raster to zmin
196  rdata.clear();
197 
198  for ( int ii = 0; ii < nxypt; ii++ )
199  {
200  rdata.append( zminr );
201  }
202 
203  // Populate raster with z values derived from a Gaussian distribution
204  // around each distribution point.
205 
206  double ssig = xrng / ( 2.0 * resol ); // sigma values
207  double fsig = yrng / ( 2.0 * resol );
208  double sssc = -1.0 / ( 4.0 * ssig * ssig ); // sig scale factors
209  double fssc = -1.0 / ( 4.0 * fsig * fsig );
210 
211  // calculate radius of decay-to-zero in Gaussian distribution
212  double dmin = sqrt( -log( 1.0e-5 ) ) * 2.0; // terms very low
213  double xdif = dmin * ssig * xinc;
214  double ydif = dmin * fsig * yinc;
215  int nxd = qRound( xdif ) + 2; // radius points with pad
216  int nyd = qRound( ydif ) + 2;
217  int hixd = nxpsc / 4; // maximum radii
218  int hiyd = nyscn / 4;
219  nxd = ( nxd < 10 ) ? 10 : ( ( nxd > hixd ) ? hixd : nxd );
220  nyd = ( nyd < 10 ) ? 10 : ( ( nyd > hiyd ) ? hiyd : nyd );
221 
222  if ( resol != 100.0 )
223  {
224  for ( int kk = 0; kk < nsol; kk++ )
225  { // spread z values for each distribution point
226  xval = solu->at( kk ).s; // x,y,z
227  yval = solu->at( kk ).k;
228  zval = solu->at( kk ).c - zminr; // use z in 0,zrng range
229 
230  int rx = (int)( ( xval - xmin ) * xinc ); // x index of this point
231  int fx = rx - nxd; // first reasonable x
232  int lx = rx + nxd; // last reasonable x
233  fx = ( fx > 0 ) ? fx : 0;
234  lx = ( lx < nxpsc ) ? lx : nxpsc;
235  int ry = (int)( ( ymax - yval ) * yinc ); // y index of this point
236  int fy = ry - nyd; // first reasonable y
237  int ly = ry + nyd; // last reasonable y
238  fy = ( fy > 0 ) ? fy : 0;
239  ly = ( ly < nyscn ) ? ly : nyscn;
240 
241  for ( int ii = fy; ii < ly; ii++ )
242  { // calculate y-term and x range for each y
243  double yras = ymax - ( (double)ii / yinc );
244  double ydif = yras - yval;
245  double yterm = exp( ydif * ydif * fssc ); // y term
246  double zterm = zval * yterm; // combine z-value,y-term
247  int kr = ii * nxpsc + fx; // start rast point index
248 
249  for ( int jj = fx; jj < lx; jj++ )
250  { // calculate x-term then calculate z-value for raster point
251  double xras = (double)jj / xinc + xmin;
252  double xdif = xras - xval;
253  double xterm = exp( xdif * xdif * sssc );
254 
255  double zin = rdata.at( kr ); // current value there
256 
257  // Output value according to Gaussian distribution factor.
258  // Note that the expression below adds zmin back in to a
259  // value that is really:
260  // zval * exp( -pow( xdif, 2.0 ) / pow( 2 * ssigma, 2.0 ) )
261  // * exp( -pow( ydif, 2.0 ) / pow( 2 * fsigma, 2.0 ) )
262  double zout = zmin + zterm * xterm;
263 
264  // only replace input if new value is greater
265  if ( zout > zin )
266  rdata.replace( kr, zout );
267 
268  kr++; // bump rast point index
269  }
270  }
271  }
272  }
273  else
274  { // for resolution=100, make all points in circle have zval value
275  for ( int kk = 0; kk < nsol; kk++ )
276  { // spread z values for each distribution point
277  xval = solu->at( kk ).s; // x,y,z
278  yval = solu->at( kk ).k;
279  zval = solu->at( kk ).c;
280 
281  int rx = (int)( ( xval - xmin ) * xinc ); // x index of this point
282  int fx = rx - nxd; // first reasonable x
283  int lx = rx + nxd; // last reasonable x
284  fx = ( fx > 0 ) ? fx : 0;
285  lx = ( lx < nxpsc ) ? lx : nxpsc;
286  int ry = (int)( ( ymax - yval ) * yinc ); // y index of this point
287  int fy = ry - nyd; // first reasonable y
288  int ly = ry + nyd; // last reasonable y
289  fy = ( fy > 0 ) ? fy : 0;
290  ly = ( ly < nyscn ) ? ly : nyscn;
291 
292  for ( int ii = fy; ii < ly; ii++ )
293  { // calculate y-term and x range for each y
294  double yras = ymax - ( (double)ii / yinc );
295  double ydif = yras - yval;
296  double yterm = exp( ydif * ydif * fssc ); // y term
297  int kr = ii * nxpsc + fx; // start rast point index
298 
299  for ( int jj = fx; jj < lx; jj++ )
300  { // calculate x-term then calculate z-value for raster point
301  double xras = (double)jj / xinc + xmin;
302  double xdif = xras - xval;
303  double xterm = exp( xdif * xdif * sssc );
304 
305  double zin = rdata.at( kr ); // current value there
306  double zdecay = yterm * xterm; // composite decay term
307 
308  // replace input if within distr radius and zout>zin
309  if ( zdecay > LO_DTERM && zval > zin)
310  rdata.replace( kr, zval );
311 
312  kr++; // bump rast point index
313  }
314  }
315  }
316  }
317 }
318 
Generated on Sun Jun 21 2015 11:40:43 for UltraScan III by   doxygen 1.8.3.1-20130324