Changeset 2339 for trunk/psLib/src/dataManip/psStats.c
- Timestamp:
- Nov 11, 2004, 10:13:57 AM (22 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/dataManip/psStats.c (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/dataManip/psStats.c
r2329 r2339 9 9 * @author GLG, MHPCC 10 10 * 11 * @version $Revision: 1. 89$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-11-1 0 23:22:32$11 * @version $Revision: 1.90 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-11-11 20:13:57 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 53 53 /*****************************************************************************/ 54 54 psVector* p_psConvertToF32(psVector* in); 55 56 55 /*****************************************************************************/ 57 56 /* GLOBAL VARIABLES */ … … 124 123 125 124 /****************************************************************************** 126 125 MISC PRIVATE STATISTICAL FUNCTIONS 127 126 128 NOTE: it is assumed that any call to these statistical functions will 129 have been preceded by a call to the psVectorStats() function. Various 130 sanity tests will only be performed in psVectorStats(). 131 Is the mask vector the same length as the data vector? 132 Is the mask vector of type PS_TYPE_U8? 133 Is the stats data structure NULL? 134 Is the in data structure NULL? 135 Is the in data structure of type PS_TYPE_F32? 127 NOTE: it is assumed that any call to these statistical functions will have 128 been preceded by a call to the psVectorStats() function. Various sanity tests 129 will only be performed in psVectorStats(). 136 130 131 XXX: Should we perform the sanity checks in each routine anyway? 132 133 XXX: For many of these private stats routines, what should be done if there 134 are now acceptable elements in the input vector (if no elements lie within 135 range, or there are no unmasked elements, or the input vector is NULL)? 136 Currently we set the value to NAN. 137 138 XXX: Optimization: many routines have an "empty" boolean variable which keeps 139 track of whether or not the vector has any valid elements. This code can 140 possibly be optimized away. 137 141 *****************************************************************************/ 138 142 /****************************************************************************** … … 146 150 Returns 147 151 NULL 148 *****************************************************************************/149 152 153 *****************************************************************************/ 150 154 void p_psVectorSampleMean(const psVector* restrict myVector, 151 155 const psVector* restrict maskVector, … … 156 160 float mean = 0.0; // The mean 157 161 psS32 count = 0; // # of points in this mean 158 float rangeMin = 0.0; // Exclude data below this159 float rangeMax = 0.0; // Exclude date above this160 162 161 163 // If PS_STAT_USE_RANGE is requested, then we enter a slightly different 162 164 // loop. 163 165 if (stats->options & PS_STAT_USE_RANGE) { 164 rangeMin = stats->min;165 rangeMax = stats->max;166 166 if (maskVector != NULL) { 167 167 for (i = 0; i < myVector->n; i++) { 168 168 // Check if the data is with the specified range 169 169 if (!(maskVal & maskVector->data.U8[i]) && 170 ( rangeMin <= myVector->data.F32[i]) &&171 (myVector->data.F32[i] <= rangeMax)) {170 (stats->min <= myVector->data.F32[i]) && 171 (myVector->data.F32[i] <= stats->max)) { 172 172 mean += myVector->data.F32[i]; 173 173 count++; … … 182 182 } else { 183 183 for (i = 0; i < myVector->n; i++) { 184 if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 184 if ((stats->min <= myVector->data.F32[i]) && 185 (myVector->data.F32[i] <= stats->max)) { 185 186 mean += myVector->data.F32[i]; 186 187 count++; … … 234 235 { 235 236 psS32 i = 0; // Loop index variable 236 float max = -PS_MAX_F32; // The calculated maximum 237 float rangeMin = 0.0; // Exclude data below this 238 float rangeMax = 0.0; // Exclude date above this 237 float max = -PS_MAX_F32; // The calculated maximum 238 psS32 empty = true; // Does this vector have valid elements? 239 239 240 240 // If PS_STAT_USE_RANGE is requested, then we enter a different loop. 241 241 if (stats->options & PS_STAT_USE_RANGE) { 242 rangeMin = stats->min;243 rangeMax = stats->max;244 242 if (maskVector != NULL) { 245 243 for (i = 0; i < myVector->n; i++) { 246 244 if (!(maskVal & maskVector->data.U8[i])) { 247 245 if ((myVector->data.F32[i] > max) && 248 (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 246 (stats->min <= myVector->data.F32[i]) && 247 (myVector->data.F32[i] <= stats->max)) { 249 248 max = myVector->data.F32[i]; 249 empty = false; 250 250 } 251 251 } … … 254 254 for (i = 0; i < myVector->n; i++) { 255 255 if ((myVector->data.F32[i] > max) && 256 (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 256 (stats->min <= myVector->data.F32[i]) && 257 (myVector->data.F32[i] <= stats->max)) { 257 258 max = myVector->data.F32[i]; 259 empty = false; 258 260 } 259 261 } … … 265 267 if (myVector->data.F32[i] > max) { 266 268 max = myVector->data.F32[i]; 269 empty = false; 267 270 } 268 271 } … … 272 275 if (myVector->data.F32[i] > max) { 273 276 max = myVector->data.F32[i]; 274 } 275 } 276 } 277 } 278 279 stats->max = max; 277 empty = false; 278 } 279 } 280 } 281 } 282 if (empty == false) { 283 stats->max = max; 284 } else { 285 stats->max = NAN; 286 } 280 287 } 281 288 … … 298 305 psS32 i = 0; // Loop index variable 299 306 float min = PS_MAX_F32; // The calculated maximum 300 float rangeMin = 0.0; // Exclude data below this 301 float rangeMax = 0.0; // Exclude date above this 307 psS32 empty = true; // Does this vector have valid elements? 302 308 303 309 if (stats->options & PS_STAT_USE_RANGE) { 304 rangeMin = stats->min;305 rangeMax = stats->max;306 310 if (maskVector != NULL) { 307 311 for (i = 0; i < myVector->n; i++) { 308 312 if (!(maskVal & maskVector->data.U8[i])) { 309 313 if ((myVector->data.F32[i] < min) && 310 (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 314 (stats->min <= myVector->data.F32[i]) && 315 (myVector->data.F32[i] <= stats->max)) { 311 316 min = myVector->data.F32[i]; 317 empty = false; 312 318 } 313 319 } … … 316 322 for (i = 0; i < myVector->n; i++) { 317 323 if ((myVector->data.F32[i] < min) && 318 (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 324 (stats->min <= myVector->data.F32[i]) && 325 (myVector->data.F32[i] <= stats->max)) { 319 326 min = myVector->data.F32[i]; 327 empty = false; 320 328 } 321 329 } … … 327 335 if (myVector->data.F32[i] < min) { 328 336 min = myVector->data.F32[i]; 337 empty = false; 329 338 } 330 339 } … … 334 343 if (myVector->data.F32[i] < min) { 335 344 min = myVector->data.F32[i]; 336 } 337 } 338 } 339 } 340 341 stats->min = min; 345 empty = false; 346 } 347 } 348 } 349 } 350 351 if (empty == false) { 352 stats->min = min; 353 } else { 354 stats->min = NAN; 355 } 342 356 } 343 357 … … 345 359 p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the 346 360 number of non-masked pixels in the vector that fall within the min/max 347 range, if given.361 range, if specified. 348 362 Inputs 349 363 myVector … … 363 377 psS32 i = 0; // Loop index variable 364 378 psS32 numData = 0; // The number of data points 365 float rangeMin = 0.0; // Exclude data below this366 float rangeMax = 0.0; // Exclude date above this367 379 368 380 if (stats->options & PS_STAT_USE_RANGE) { 369 rangeMin = stats->min;370 rangeMax = stats->max;371 381 if (maskVector != NULL) { 372 382 for (i = 0; i < myVector->n; i++) { 373 383 if (!(maskVal & maskVector->data.U8[i]) && 374 (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 384 (stats->min <= myVector->data.F32[i]) && 385 (myVector->data.F32[i] <= stats->max)) { 375 386 numData++; 376 387 } … … 378 389 } else { 379 390 for (i = 0; i < myVector->n; i++) { 380 if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 391 if ((stats->min <= myVector->data.F32[i]) && 392 (myVector->data.F32[i] <= stats->max)) { 381 393 numData++; 382 394 } … … 394 406 } 395 407 } 408 396 409 return (numData); 397 410 } … … 408 421 NULL 409 422 *****************************************************************************/ 410 voidp_psVectorSampleMedian(const psVector* restrict myVector,423 bool p_psVectorSampleMedian(const psVector* restrict myVector, 411 424 const psVector* restrict maskVector, 412 425 psU32 maskVal, … … 418 431 psS32 count = 0; // # of points in this mean? 419 432 psS32 nValues = 0; // # of points in vector 420 float rangeMin = 0.0; // Exclude data below this421 float rangeMax = 0.0; // Exclude date above this422 433 423 434 // Determine how many data points fit inside this min/max range … … 425 436 nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); 426 437 438 // XXX: Is a warning message appropriate? What value should we set the 439 // sampleMedian to? Should we generate an error? 440 if (nValues <= 0) { 441 printf("WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n"); 442 return(true); 443 stats->sampleMedian = 0; 444 } 445 427 446 // Allocate temporary vectors for the data. 428 447 unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); 429 sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);430 448 431 449 // Determine if we must only use data points within a min/max range. 432 450 if (stats->options & PS_STAT_USE_RANGE) { 433 rangeMin = stats->min;434 rangeMax = stats->max;435 436 451 // Store all non-masked data points within the min/max range 437 452 // into the temporary vectors. … … 440 455 for (i = 0; i < myVector->n; i++) { 441 456 if (!(maskVal & maskVector->data.U8[i]) && 442 (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 443 unsortedVector->data.F32[count++] = maskVector->data.F32[i]; 457 (stats->min <= myVector->data.F32[i]) && 458 (myVector->data.F32[i] <= stats->max)) { 459 unsortedVector->data.F32[count] = myVector->data.F32[i]; 460 count++; 444 461 } 445 462 } 446 463 } else { 447 464 for (i = 0; i < myVector->n; i++) { 448 if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) { 449 unsortedVector->data.F32[count++] = myVector->data.F32[i]; 465 if ((stats->min <= myVector->data.F32[i]) && 466 (myVector->data.F32[i] <= stats->max)) { 467 unsortedVector->data.F32[count] = myVector->data.F32[i]; 468 count++; 450 469 } 451 470 } … … 457 476 for (i = 0; i < myVector->n; i++) { 458 477 if (!(maskVal & maskVector->data.U8[i])) { 459 unsortedVector->data.F32[count++] = myVector->data.F32[i]; 478 unsortedVector->data.F32[count] = myVector->data.F32[i]; 479 count++; 460 480 } 461 481 } … … 467 487 } 468 488 // Sort the temporary vectors. 469 psVectorSort(sortedVector, unsortedVector); 489 sortedVector = psVectorSort(sortedVector, unsortedVector); 490 if (sortedVector == NULL) { 491 psError(PS_ERR_UNEXPECTED_NULL, 492 false, 493 PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM); 494 return(false); 495 } 470 496 471 497 // Calculate the median exactly. … … 480 506 psFree(unsortedVector); 481 507 psFree(sortedVector); 508 509 // Return "true" on success. 510 return(true); 482 511 } 483 512 … … 485 514 p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input 486 515 robustHistogram with a Gaussian of width sigma. 516 517 XXX: Only PS_TYPE_F32 is supported. 518 519 XXX: Write a general routine which smoothes a psVector. This routine should 520 call that. Is that possible? 487 521 *****************************************************************************/ 488 522 psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, 489 523 float sigma) 490 524 { 525 PS_PTR_CHECK_NULL(robustHistogram, NULL); 526 PS_PTR_CHECK_NULL(robustHistogram->bounds, NULL); 527 491 528 psS32 i = 0; // Loop index variable 492 529 psS32 j = 0; // Loop index variable … … 504 541 x.type.type = PS_TYPE_F32; 505 542 for (i = 0; i < numBins; i++) { 543 // Determine the midpoint of bin i. 506 544 iMid = (robustHistogram->bounds->data.F32[i] + 507 545 robustHistogram->bounds->data.F32[i+1]) / 2.0; 546 547 548 // We determine the bin numbers corresponding to a range of data 549 // values surrounding iMid. The ranges is of size 550 // s*PS_GAUSS_WIDTH*sigma 508 551 509 552 // YYY: The p_psVectorBinDisect() routine does much of the work of … … 516 559 if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { 517 560 jMin = p_psVectorBinDisect(robustHistogram->bounds, &x); 561 if (jMin < 0) { 562 psError(PS_ERR_UNEXPECTED_NULL, 563 false, 564 PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); 565 return(NULL); 566 } 518 567 } else if (x.data.F32 <= firstBound) { 519 568 jMin = 0;
Note:
See TracChangeset
for help on using the changeset viewer.
