001/* 002 * This file is part of McIDAS-V 003 * 004 * Copyright 2007-2024 005 * Space Science and Engineering Center (SSEC) 006 * University of Wisconsin - Madison 007 * 1225 W. Dayton Street, Madison, WI 53706, USA 008 * https://www.ssec.wisc.edu/mcidas/ 009 * 010 * All Rights Reserved 011 * 012 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and 013 * some McIDAS-V source code is based on IDV and VisAD source code. 014 * 015 * McIDAS-V is free software; you can redistribute it and/or modify 016 * it under the terms of the GNU Lesser Public License as published by 017 * the Free Software Foundation; either version 3 of the License, or 018 * (at your option) any later version. 019 * 020 * McIDAS-V is distributed in the hope that it will be useful, 021 * but WITHOUT ANY WARRANTY; without even the implied warranty of 022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 023 * GNU Lesser Public License for more details. 024 * 025 * You should have received a copy of the GNU Lesser Public License 026 * along with this program. If not, see https://www.gnu.org/licenses/. 027 */ 028 029package edu.wisc.ssec.mcidasv.adt; 030 031import static java.lang.Math.PI; 032import static java.lang.Math.abs; 033import static java.lang.Math.cos; 034import static java.lang.Math.log; 035import static java.lang.Math.max; 036import static java.lang.Math.sqrt; 037import static java.lang.Math.pow; 038 039import java.io.IOException; 040 041import javax.swing.JOptionPane; 042 043class Remap { 044 /** Block size (bytes) input file */ 045 int in_bfw; 046 047 /** Block size (bytes) output file */ 048 int out_bfw; 049 050 /** Number of splines/line */ 051 int nspl; 052 053 /** Number of splines/elem */ 054 int nspe; 055 056 /** Source blocksize */ 057 int slb; 058 059 /** Dest blocksize */ 060 int dlb; 061 062 /** Number of corners in line */ 063 int ncl; 064 065 /** Number of corners in elem */ 066 int nce; 067} 068 069/** TIFF header */ 070class TiffHeader { 071 /** Byte order */ 072 int order; 073 074 /** Version */ 075 int version; 076 077 /** Pointer */ 078 int point; 079} 080 081/** 082 * TIFF directory record. 083 */ 084class TiffRecord { 085 /** TIFF tag */ 086 int tag; 087 088 /** Data type */ 089 int type; 090 091 /** Length */ 092 int length; 093 094 /** Pointer or value */ 095 int voff; 096} 097 098class DisVars { 099 /** Output number of lines */ 100 double xrectl; 101 102 /** Output number of elems */ 103 double xrecte; 104} 105 106class TiffVars { 107 int nbits; 108 int photo; 109 int unit; 110 int in_lines; 111 int in_elems; 112 int out_lines; 113 int out_elems; 114} 115 116public class Auto { 117 118 private static int[][] MoatMaskFlagField = new int[200][200]; 119 120 private static int[][] BlackWhiteFieldArray = new int[200][200]; 121 122 private static double[][] IRData_Remap_Latitude = new double[200][200]; 123 private static double[][] IRData_Remap_Longitude = new double[200][200]; 124 125 private static double[][] IRData_Remap_Temperature = new double[200][200]; 126 127 private static int IRData_Remap_NumberRows; 128 129 private static int IRData_Remap_NumberColumns; 130 131 private static double[][] NSTempGradientArray = new double[200][200]; 132 133 private static double[][] EWTempGradientArray = new double[200][200]; 134 135 private static double[][] SpiralCenterAnalysisField = new double[50000][3]; 136 137 private static double[][] RingScoreAnalysisField = new double[750000][3]; 138 139 private static int[][] CircleFilterRowArray = new int[100][2000]; 140 141 private static int[][] CircleFilterColumnArray = new int[100][2000]; 142 143 /* I don't know this size for sure */ 144 /** Array containing line coordinates */ 145 private static double[] LineCoordinateArray = new double[5000]; 146 147 /* I don't know this size for sure */ 148 /** Array containing element coordinates */ 149 private static double[] ElementCoordinateArray = new double[5000]; 150 151 private static Remap remap_vars = new Remap(); 152 153 private static TiffVars tiff_vars = new TiffVars(); 154 155 private static double RING_WIDTH = 4.0; 156 157 /** Minimum block size (bytes) for domap */ 158 private static int MINBFW = 500000; 159 160 /** Minimum block size (lines) */ 161 private static int MINBLKSIZ = 10; 162 163 public Auto() { 164 IRData_Remap_NumberRows = 0; 165 IRData_Remap_NumberColumns = 0; 166 } 167 168 /** 169 * Determine storm position at time CurrentTime using NHC/JTWC forecast 170 * discussion products. 171 * 172 * Time and location information from these products are then interpolated 173 * to time in question to derive a estimated storm position. If position 174 * estimation cannot be calculated, a lat/lon position of -99.5/-999.5 will 175 * be returned. Inputs : None Outputs : Latitude_Return - estimated latitude 176 * position Longitude_Return - estimated longitude position 177 * PositioningMethodID_Return - method used to derive storm location Return 178 * : -43 : Error w/ forecast file open and BAD extrapolation -44 : Invalid 179 * forecast file and BAD extrapolation -45 : Error w/ forecast file read and 180 * BAD extrapolation -46 : Error w/ forecast interpolation and BAD 181 * extrapolation 42 : GOOD INTERPOLATION 43 : Error w/ forecast file open 182 * but GOOD EXTRAPOLATION 44 : Invalid forecast file but GOOD extrapolation 183 * 45 : Error w/ forecast file read but GOOD extrapolation 46 : Error w/ 184 * forecast interpolation but GOOD EXTRAPOLATION 0 : Subroutine Error 185 */ 186 public double[] AutoMode1(String ForecastFile, int ForecastFileType) throws IOException { 187 188 int RetID = 0; 189 int PositioningMethodID = 0; 190 int InterpolatedReturnFlag = 0; 191 double InterpolatedLatitude = -99.99; 192 double InterpolatedLongitude = -99.99; 193 double InterpolatedIntensity = -99.99; 194 boolean UseExtrapPositionTF = false; 195 196 /* Read Forecast File */ 197 double ThresholdTime = 24.0; 198 double[] ForecastFileOutput = Forecasts.ReadForecasts(ForecastFile, ForecastFileType, 199 ThresholdTime); 200 201 if (ForecastFileOutput == null) { 202 JOptionPane.showMessageDialog(null, 203 "Unable to process Forecast File, please ensure a valid file is present."); 204 return null; 205 } 206 207 InterpolatedReturnFlag = (int) ForecastFileOutput[0]; 208 InterpolatedLatitude = ForecastFileOutput[1]; 209 InterpolatedLongitude = ForecastFileOutput[2]; 210 InterpolatedIntensity = ForecastFileOutput[3]; 211 /* History.IRCurrentRecord.latitude = ForecastLatitude; */ 212 /* History.IRCurrentRecord.longitude = ForecastLongitude; */ 213 214 /* 215 * System.out.printf("InterpolatedReturnFlag=%d\n",InterpolatedReturnFlag 216 * ); 217 */ 218 if (InterpolatedReturnFlag != 0) { 219 PositioningMethodID = 6; 220 UseExtrapPositionTF = true; 221 /* 222 * -1 will give 45=invalid interpretation -2 will give 44=invalid 223 * file type -4 will give 42=error closing file -5 will give 224 * 41=error opening file 225 */ 226 RetID = 46 + InterpolatedReturnFlag; 227 } else { 228 /* Good interpolated values... use 'em */ 229 PositioningMethodID = 1; 230 UseExtrapPositionTF = false; 231 RetID = 43; 232 } 233 234 /* 235 * try to extrapolate storm location from previous storm locations in 236 * history file 237 */ 238 if (UseExtrapPositionTF) { 239 /* call slopecal to get y-intercept values for lat/lon values */ 240 InterpolatedLatitude = Functions.adt_slopecal(12.0, 3); 241 InterpolatedLongitude = Functions.adt_slopecal(12.0, 4); 242 if ((abs(InterpolatedLatitude) > 90.0) || (abs(InterpolatedLongitude) > 180.0)) { 243 /* invalid interp and extrap... negative error code returns */ 244 PositioningMethodID = 0; 245 RetID = -1 * RetID; 246 } else { 247 /* invalid interp but good extrap... positive error code returns */ 248 PositioningMethodID = 6; 249 RetID = 46; 250 } 251 } 252 253 return new double[] { (double) RetID, InterpolatedLatitude, InterpolatedLongitude, 254 InterpolatedIntensity, (double) PositioningMethodID }; 255 } 256 257 /** 258 * Additional automatic positioning of storm center location using official 259 * forecasts from NHC or JTWC as input. 260 * 261 * Storm location will be estimated using spiral fitting and ring fitting 262 * routines derived by Tony Wimmers in his MatLab routines (which have been 263 * converted). 264 * 265 * The final position will be determined utilizing empirically defined 266 * confidence factors for each method. The final storm position will be 267 * returned along with a position determination flag. 268 * 269 * @param InputLatitudePosition 270 * Storm center latitude. 271 * @param InputLongitudePosition 272 * Storm center longitude. 273 * 274 * Outputs : Latitude_Return - final storm center latitude 275 * position Longitude_Return - final storm center longitude 276 * position PositioningMethodID_Return - method used to derive 277 * storm location 0-error 1-interpolation of operational forecast 278 * 2-Laplacian analysis (not used anymore) 3-Warm Spot location 279 * 4-10^ log spiral analysis 5-Combo method of spiral and ring 280 * analyses 6-linear extrapolation from prior locations Return : 281 * Error flag = 0 282 */ 283 public static double[] AutoMode2(double InputLatitudePosition, double InputLongitudePosition) 284 throws IOException { 285 int XInc; 286 int YInc; 287 double SpiralCenterLatitude = -99.9; 288 double SpiralCenterLongitude = -99.9; 289 double SpiralCenterScore = -1.0; 290 double RingFitLatitude = -99.9; 291 double RingFitLongitude = -99.9; 292 double RingFitScore = -1.0; 293 double FinalAutoLatitude = -99.9; 294 double FinalAutoLongitude = -99.9; 295 double FinalAutoScore = -99.9; 296 int FinalAutoFixMethod = 1; 297 boolean DatelineCrossTF = false; 298 boolean DoRemappingTF = true; 299 300 /* 301 * NEED TO MULTIPLY LONGITUDES BY -1.0 SINCE ADT ROUTINES WERE DEVELOPED 302 * ON McIDAS-X, which uses West/East Longitudes of +1.0/-1.0. McV is 303 * reversed. 304 */ 305 InputLongitudePosition = -1.0 * InputLongitudePosition; // flip for ADT 306 // routines 307 for (YInc = 0; YInc < Data.IRData_NumberRows; YInc++) { 308 for (XInc = 0; XInc < Data.IRData_NumberColumns; XInc++) { 309 Data.IRData_Longitude[YInc][XInc] = (float) -1.0 310 * Data.IRData_Longitude[YInc][XInc]; 311 } 312 } 313 314 if (DoRemappingTF) { 315 System.out.printf("REMAPPING DATA\n"); 316 IRData_Remap_NumberRows = Data.IRData_NumberRows; 317 IRData_Remap_NumberColumns = Data.IRData_NumberColumns; 318 319 double LongitudeIncrement = Data.IRData_Longitude[0][0] - Data.IRData_Longitude[0][1]; 320 double LatitudeIncrement = Data.IRData_Latitude[0][0] - Data.IRData_Latitude[1][1]; 321 322 if (abs(LongitudeIncrement - LatitudeIncrement) < 0.001) { 323 System.out.printf("already remapped\n"); 324 /* data is already remapped */ 325 /* crosses dateline check */ 326 int XSizeMax = Data.IRData_NumberColumns - 1; 327 int YSizeMax = Data.IRData_NumberRows - 1; 328 double NWCornerLongitude = Data.IRData_Longitude[0][0]; 329 double NECornerLongitude = Data.IRData_Longitude[0][XSizeMax]; 330 double SWCornerLongitude = Data.IRData_Longitude[YSizeMax][0]; 331 double SECornerLongitude = Data.IRData_Longitude[YSizeMax][XSizeMax]; 332 /* 333 * if((NWCornerLongitude<NECornerLongitude)||(SWCornerLongitude< 334 * SECornerLongitude)) { DatelineCrossTF = true; } 335 */ 336 if ((NWCornerLongitude > NECornerLongitude) 337 || (SWCornerLongitude > SECornerLongitude)) { 338 DatelineCrossTF = true; 339 NWCornerLongitude = NWCornerLongitude + 360.0; 340 341 } 342 for (YInc = 0; YInc < Data.IRData_NumberRows; YInc++) { 343 for (XInc = 0; XInc < Data.IRData_NumberColumns; XInc++) { 344 IRData_Remap_Longitude[YInc][XInc] = Data.IRData_Longitude[YInc][XInc]; 345 /* check for dateline crossing */ 346 if (DatelineCrossTF && (IRData_Remap_Longitude[YInc][XInc] < 0.0)) { 347 IRData_Remap_Latitude[YInc][XInc] = Data.IRData_Latitude[YInc][XInc] + 360.0; 348 } 349 IRData_Remap_Latitude[YInc][XInc] = Data.IRData_Latitude[YInc][XInc]; 350 IRData_Remap_Temperature[YInc][XInc] = Data.IRData_Temperature[YInc][XInc]; 351 } 352 } 353 IRData_Remap_NumberColumns = Data.IRData_NumberColumns; 354 IRData_Remap_NumberRows = Data.IRData_NumberRows; 355 } else { 356 /* remap data to rectilinear projection */ 357 RemapData(); 358 System.out.printf("COMPLETED REMAPPING DATA\n"); 359 } 360 361 /* perform spiral analysis to determine storm center location */ 362 double SpiralReturn[] = SpiralCenterLowRes(InputLatitudePosition, 363 InputLongitudePosition); 364 SpiralCenterLatitude = SpiralReturn[0]; 365 SpiralCenterLongitude = SpiralReturn[1]; 366 SpiralCenterScore = SpiralReturn[2]; 367 368 /* 369 * System.out.printf("Spiral Score : lat=%f lon=%f score=%f\n", 370 * SpiralCenterLatitude, SpiralCenterLongitude,SpiralCenterScore); 371 */ 372 373 /* 374 * redefine first guess for ring analysis as input storm location 375 * point 376 */ 377 double RingFitFirstGuessLatitude = SpiralCenterLatitude; 378 double RingFitFirstGuessLongitude = SpiralCenterLongitude; 379 380 /* calculate Moat Mask for false eye check */ 381 double MoatMaskTempThreshold = 237.0; 382 double MoatMaskMaxRadiusDegree = 0.50; 383 MoatMaskCalc(MoatMaskTempThreshold, MoatMaskMaxRadiusDegree, 1); 384 385 /* perform ring analysis to determine storm center location */ 386 double RingReturn[] = RingFit(RingFitFirstGuessLatitude, RingFitFirstGuessLongitude); 387 RingFitLatitude = RingReturn[0]; 388 RingFitLongitude = RingReturn[1]; 389 RingFitScore = RingReturn[2]; 390 391 /* 392 * System.out.printf("Ring Score : lat=%f lon=%f score=%f \n", 393 * RingFitLatitude,RingFitLongitude,RingFitScore); 394 */ 395 /* 396 * System.out.printf( 397 * "fglat=%f fglon=%f SpiralCenterLatitude=%f SpiralCenterLongitude=%f rglat=%f rglon=%f\n" 398 * , InputLatitudePosition,InputLongitudePosition, 399 * SpiralCenterLatitude,SpiralCenterLongitude, 400 * RingFitLatitude,RingFitLongitude); 401 */ 402 403 /* caluculate confidence factor for combined spiral/ring analyses */ 404 double ScoresReturn[] = CalcScores(InputLatitudePosition, InputLongitudePosition, 405 SpiralCenterLatitude, SpiralCenterLongitude, SpiralCenterScore, 406 RingFitLatitude, RingFitLongitude, RingFitScore); 407 FinalAutoLatitude = ScoresReturn[0]; 408 FinalAutoLongitude = ScoresReturn[1]; 409 FinalAutoScore = ScoresReturn[2]; 410 FinalAutoFixMethod = (int) ScoresReturn[3]; 411 /* 412 * System.out.printf("SPIRAL CENTER : lat=%f lon=%f SCORE=%f \n", 413 * SpiralCenterLatitude,SpiralCenterLongitude,SpiralCenterScore); 414 * System.out.printf("RING CENTER : lat=%f lon=%f SCORE=%f \n", 415 * RingFitLatitude,RingFitLongitude,RingFitScore); 416 * System.out.printf("FINAL CENTER : lat=%f lon=%f SCORE=%f \n", 417 * FinalAutoLatitude,FinalAutoLongitude,FinalAutoScore); 418 */ 419 } else { 420 SpiralCenterLatitude = -99.9; 421 SpiralCenterLongitude = -99.9; 422 SpiralCenterScore = -1.0; 423 RingFitLatitude = -99.9; 424 RingFitLatitude = -99.9; 425 RingFitScore = -1.0; 426 FinalAutoLatitude = -99.9; 427 FinalAutoLongitude = -99.9; 428 FinalAutoScore = -99.9; 429 FinalAutoFixMethod = 1; 430 } 431 432 int PositioningMethodID = History.IRCurrentRecord.autopos; 433 434 double LocationReturn[] = PickFinalLocation(PositioningMethodID, InputLatitudePosition, 435 InputLongitudePosition, FinalAutoLatitude, FinalAutoLongitude, FinalAutoScore, 436 FinalAutoFixMethod); 437 438 double FinalLatitude = LocationReturn[0]; 439 double FinalLongitude = LocationReturn[1]; 440 PositioningMethodID = (int) LocationReturn[3]; 441 442 if (FinalLongitude > 180.0) { 443 FinalLongitude = FinalLongitude - 360.0; 444 } 445 446 /* 447 * System.out.printf("TAC CENTER : lat=%f lon=%f \n", 448 * InputLatitudePosition,InputLongitudePosition); 449 */ 450 /* 451 * System.out.printf( 452 * "AUTO CENTER : lat=%f lon=%f SCORE=%f METHOD=%d\n", 453 * FinalAutoLatitude,FinalAutoLongitude, 454 * FinalAutoScore,FinalAutoFixMethod); 455 */ 456 /* 457 * System.out.printf( 458 * "FINAL LOCATION: lat=%f lon=%f SCORE=%f METHOD=%d\n", 459 * FinalLatitude,FinalLongitude, FinalScoreValue,PositioningMethodID); 460 */ 461 462 /* 463 * need to flip final Longitude value back to McV format by multiplying 464 * by -1.0 465 */ 466 FinalLongitude = -1.0 * FinalLongitude; 467 468 return new double[] { FinalLatitude, FinalLongitude, (double) PositioningMethodID }; 469 } 470 471 private static double[] SpiralCenterLowRes(double InputLatitude, double InputLongitude) { 472 473 double[][] IRData_NormCoord_Latitude = new double[200][200]; 474 double[][] IRData_NormCoord_Longitude = new double[200][200]; 475 double[][] IRData_NormCoord_Temperature = new double[200][200]; 476 477 int XInc, YInc; 478 int IRData_NormCoord_NumberRows; 479 int IRData_NormCoord_NumberColumns; 480 int NumPoints; 481 int RingWidthLocal = (int) RING_WIDTH; 482 483 double XOff, YOff; 484 double SearchRadius = 0.0; 485 double CrossScoreSum = 0.0; 486 double ProxyValueX = 0.0; 487 double ProxyValueY = 0.0; 488 double DenominatorValue = 0.0; 489 double SpiralValueX = 0.0; 490 double SpiralValueY = 0.0; 491 double RawCrossScore = 0.0; 492 double CrossScoreClean = 0.0; 493 double SpiralMeanCrossScore = 0.0; 494 double SpiralMeanCrossMaxScore = 0.0; 495 double SpiralMeanCrossMaxYLocation = 0.0; 496 double SpiralMeanCrossMaxXLocation = 0.0; 497 498 double SpiralCenterLatitude = -99.99; 499 double SpiralCenterLongitude = -999.99; 500 double SpiralCenterScore = 0.0; 501 502 double Alpha = 5.0 * PI / 180.0; 503 double OutsideFactor = 0.62; 504 505 double OuterSearchRadiusDegree = 1.75; 506 double CourseGridSpacingDegree = 0.2; 507 double FineGridSpacingDegree = 0.1; 508 double FilterDiscSize = pow(OuterSearchRadiusDegree + (2.0 * FineGridSpacingDegree), 2); 509 double AlphaPOWP1 = 1.0 + pow(Alpha, 2); 510 511 int ImageResolution = (int) Data.GetCurrentImageResolution(); 512 int IncAddVal = (ImageResolution > RingWidthLocal) ? 1 513 : (RingWidthLocal - ImageResolution + 1); 514 515 double SignFactor = abs(InputLatitude) / InputLatitude; 516 517 YInc = (IRData_Remap_NumberRows) - 1; 518 if ((InputLongitude < 0.0) 519 && ((IRData_Remap_Longitude[0][0] > 180.0) || (IRData_Remap_Longitude[YInc][0] > 180.0))) { 520 /* 521 * System.out.printf("DATELINE CROSS... changing InputLongitude from %f" 522 * ,InputLongitude); 523 */ 524 InputLongitude = InputLongitude + 360.0; 525 } 526 527 for (YInc = 0; YInc < IRData_Remap_NumberRows; YInc++) { 528 for (XInc = 0; XInc < IRData_Remap_NumberColumns; XInc++) { 529 /* compute normalized coordinate system arrays */ 530 IRData_NormCoord_Longitude[YInc][XInc] = (IRData_Remap_Longitude[YInc][XInc] - InputLongitude) 531 * (cos(PI * InputLatitude / 180.0)); 532 IRData_NormCoord_Latitude[YInc][XInc] = IRData_Remap_Latitude[YInc][XInc] 533 - InputLatitude; 534 IRData_NormCoord_Temperature[YInc][XInc] = IRData_Remap_Temperature[YInc][XInc]; 535 /* 536 * System.out.printf( 537 * "YInc=%d XInc=%d lat=%f lon=%f lat=%f lon=%f temp=%f\n", 538 * YInc,XInc, IRData_Remap_Latitude[YInc][XInc], 539 * IRData_Remap_Longitude[YInc][XInc], 540 * IRData_NormCoord_Latitude[YInc][XInc], 541 * IRData_NormCoord_Longitude[YInc][XInc], 542 * IRData_Remap_Temperature[YInc][XInc]); 543 */ 544 } 545 } 546 IRData_NormCoord_NumberColumns = IRData_Remap_NumberColumns; 547 IRData_NormCoord_NumberRows = IRData_Remap_NumberRows; 548 549 /* determine lat/lon grid increment */ 550 /* W to E gradient */ 551 double LongitudeIncrement = abs(IRData_NormCoord_Longitude[0][0] 552 - IRData_NormCoord_Longitude[0][1]); 553 /* N to S gradient */ 554 double LatitudeIncrement = abs(IRData_NormCoord_Latitude[0][0] 555 - IRData_NormCoord_Latitude[1][0]); 556 /* 557 * This is to determine longitude multiplier factor... original routines 558 * were developed using negative, but McIDAS is positive in WH. So if 559 * LongitudeMultFactor is negative, we know (assume) we are using 560 * non-McIDAS input imagery/values, otherwise make LongitudeMultFactor 561 * positive. This all assumes that image is loaded from NW to SE 562 */ 563 double LongitudeMultFactor = IRData_NormCoord_Longitude[0][0] 564 - IRData_NormCoord_Longitude[0][1]; 565 LongitudeMultFactor = (LongitudeMultFactor < 0.0) ? 1.0 : -1.0; 566 /* 567 * System.out.printf( 568 * "LatitudeIncrement=%f LongitudeIncrement=%f LongitudeMultFactor=%f\n" 569 * , LatitudeIncrement, LongitudeIncrement,LongitudeMultFactor); 570 */ 571 572 /* calculate gradient field */ 573 Gradient(IRData_NormCoord_Temperature, IRData_NormCoord_NumberColumns, 574 IRData_NormCoord_NumberRows, LongitudeIncrement, LatitudeIncrement); 575 576 for (YInc = 0; YInc < IRData_NormCoord_NumberRows - 1; YInc++) { 577 for (XInc = 0; XInc < IRData_NormCoord_NumberColumns - 1; XInc++) { 578 /* 579 * System.out.printf( 580 * "YInc=%d XInc=%d remap:lat=%f lon=%f normcoord:lat=%f lon=%f NSTempGradientArray=%f EWTempGrad=%f\n" 581 * , YInc,XInc, 582 * IRData_Remap_Latitude[YInc][XInc],IRData_Remap_Longitude 583 * [YInc][XInc], 584 * IRData_NormCoord_Latitude[YInc][XInc],IRData_NormCoord_Longitude 585 * [YInc][XInc], 586 * NSTempGradientArray[YInc][XInc],EWTempGradientArray 587 * [YInc][XInc]); 588 */ 589 590 double GradientOriginalMag = sqrt(pow(NSTempGradientArray[YInc][XInc], 2) 591 + pow(EWTempGradientArray[YInc][XInc], 2)); 592 double GradientLOGMag = log(1.0 + GradientOriginalMag); 593 double GradientLOGReduction = 0.0; 594 if (GradientLOGMag != 0.0) { 595 GradientLOGReduction = GradientLOGMag / GradientOriginalMag; 596 } 597 NSTempGradientArray[YInc][XInc] = GradientLOGReduction 598 * NSTempGradientArray[YInc][XInc]; 599 EWTempGradientArray[YInc][XInc] = GradientLOGReduction 600 * EWTempGradientArray[YInc][XInc]; 601 } 602 } 603 604 /* COURSE GRID */ 605 /* calculate cross product score at each grid point */ 606 /* XInc/YInc are "starting point" coordinates */ 607 SpiralMeanCrossMaxScore = -99.0; 608 double SearchRadiusMaximum = pow(OuterSearchRadiusDegree 609 + ((2.0 * CourseGridSpacingDegree) / 3.0), 2); 610 for (XOff = -OuterSearchRadiusDegree; XOff <= OuterSearchRadiusDegree; XOff = XOff 611 + CourseGridSpacingDegree) { 612 for (YOff = -OuterSearchRadiusDegree; YOff <= OuterSearchRadiusDegree; YOff = YOff 613 + CourseGridSpacingDegree) { 614 /* XOff/YOff are offset coordinates from "starting point" */ 615 SearchRadius = pow(XOff, 2) + pow(YOff, 2); 616 if (SearchRadius <= SearchRadiusMaximum) { 617 CrossScoreSum = 0.0; 618 NumPoints = 0; 619 for (YInc = 1; YInc < IRData_NormCoord_NumberRows - 2; YInc = YInc + IncAddVal) { 620 for (XInc = 1; XInc < IRData_NormCoord_NumberColumns - 2; XInc = XInc 621 + IncAddVal) { 622 SearchRadius = pow(IRData_NormCoord_Longitude[YInc][XInc], 2) 623 + pow(IRData_NormCoord_Latitude[YInc][XInc], 2); 624 if (SearchRadius < FilterDiscSize) { 625 ProxyValueX = LongitudeMultFactor 626 * IRData_NormCoord_Longitude[YInc][XInc] - XOff; 627 ProxyValueY = IRData_NormCoord_Latitude[YInc][XInc] - YOff; 628 DenominatorValue = sqrt((AlphaPOWP1) 629 * (pow(ProxyValueX, 2) + pow(ProxyValueY, 2))); 630 SpiralValueX = ((Alpha * ProxyValueX) + (SignFactor * ProxyValueY)) 631 / DenominatorValue; 632 SpiralValueY = ((Alpha * ProxyValueY) - (SignFactor * ProxyValueX)) 633 / DenominatorValue; 634 RawCrossScore = (SpiralValueX * NSTempGradientArray[YInc][XInc]) 635 - (SpiralValueY * EWTempGradientArray[YInc][XInc]); 636 CrossScoreClean = max(0.0, -RawCrossScore) 637 + (OutsideFactor * max(0.0, RawCrossScore)); 638 CrossScoreSum = CrossScoreSum + CrossScoreClean; 639 /* 640 * System.out.printf( 641 * "%d %d : lat=%f lon=%f xoff=%f yoff=%f ", 642 * "ProxyValueX=%f ProxyValueY=%f ", 643 * "SpiralValueX=%f SpiralValueY=%f ", 644 * "rawScore=%f clean=%f\n", 645 * IRData_Remap_Latitude[YInc][XInc], 646 * IRData_Remap_Longitude[YInc][XInc], 647 * XOff,YOff,ProxyValueX,ProxyValueY, 648 * SpiralValueX,SpiralValueY, 649 * RawCrossScore,CrossScoreClean); 650 */ 651 NumPoints++; 652 } 653 } 654 } 655 /* calculate mean of all values in CrossScore array */ 656 SpiralMeanCrossScore = CrossScoreSum / (double) NumPoints; 657 /* store location of maximum score position */ 658 if (SpiralMeanCrossScore > SpiralMeanCrossMaxScore) { 659 SpiralMeanCrossMaxScore = SpiralMeanCrossScore; 660 SpiralMeanCrossMaxYLocation = YOff; 661 SpiralMeanCrossMaxXLocation = XOff; 662 } 663 } 664 } 665 } 666 667 /* 668 * System.out.printf("course grid : y=%f x=%f max=%f\n", 669 * SpiralMeanCrossMaxYLocation,SpiralMeanCrossMaxXLocation, 670 * SpiralMeanCrossMaxScore); 671 */ 672 673 /* FINE GRID */ 674 int SpiralCount = 1; 675 double FineGridXMinimum = SpiralMeanCrossMaxXLocation - CourseGridSpacingDegree; 676 double FineGridXMaximum = SpiralMeanCrossMaxXLocation + CourseGridSpacingDegree; 677 double FineGridYMimimum = SpiralMeanCrossMaxYLocation - CourseGridSpacingDegree; 678 double FineGridYMaximum = SpiralMeanCrossMaxYLocation + CourseGridSpacingDegree; 679 SpiralMeanCrossMaxScore = -99.0; 680 for (XOff = FineGridXMinimum; XOff <= FineGridXMaximum; XOff = XOff + FineGridSpacingDegree) { 681 for (YOff = FineGridYMimimum; YOff <= FineGridYMaximum; YOff = YOff 682 + FineGridSpacingDegree) { 683 /* XOff/YOff are offset coordinates from "starting point" */ 684 CrossScoreSum = 0.0; 685 NumPoints = 0; 686 for (YInc = 1; YInc < IRData_NormCoord_NumberRows - 2; YInc++) { 687 for (XInc = 1; XInc < IRData_NormCoord_NumberColumns - 2; XInc++) { 688 SearchRadius = Math.pow(IRData_NormCoord_Longitude[YInc][XInc], 2) 689 + Math.pow(IRData_NormCoord_Latitude[YInc][XInc], 2); 690 if (SearchRadius < FilterDiscSize) { 691 ProxyValueX = LongitudeMultFactor 692 * IRData_NormCoord_Longitude[YInc][XInc] - XOff; 693 ProxyValueY = IRData_NormCoord_Latitude[YInc][XInc] - YOff; 694 DenominatorValue = Math.sqrt((AlphaPOWP1) 695 * (Math.pow(ProxyValueX, 2) + Math.pow(ProxyValueY, 2))); 696 SpiralValueX = (Alpha * ProxyValueX + (SignFactor * ProxyValueY)) 697 / DenominatorValue; 698 SpiralValueY = (Alpha * ProxyValueY - (SignFactor * ProxyValueX)) 699 / DenominatorValue; 700 RawCrossScore = (SpiralValueX * NSTempGradientArray[YInc][XInc]) 701 - (SpiralValueY * EWTempGradientArray[YInc][XInc]); 702 CrossScoreClean = Math.max(0.0, -RawCrossScore) 703 + (OutsideFactor * Math.max(0.0, RawCrossScore)); 704 CrossScoreSum = CrossScoreSum + CrossScoreClean; 705 NumPoints++; 706 } 707 } 708 } 709 /* calculate mean of all values in CrossScore array */ 710 SpiralMeanCrossScore = CrossScoreSum / (double) NumPoints; 711 /* store location of maximum score position */ 712 if (SpiralMeanCrossScore > SpiralMeanCrossMaxScore) { 713 SpiralMeanCrossMaxScore = SpiralMeanCrossScore; 714 SpiralMeanCrossMaxYLocation = YOff; 715 SpiralMeanCrossMaxXLocation = XOff; 716 } 717 SpiralCenterAnalysisField[SpiralCount][0] = SpiralMeanCrossScore; 718 SpiralCenterAnalysisField[SpiralCount][1] = (double) YOff; 719 SpiralCenterAnalysisField[SpiralCount][2] = (double) XOff; 720 SpiralCenterAnalysisField[0][0] = (double) SpiralCount; 721 SpiralCount++; 722 } 723 } 724 725 /* 726 * System.out.printf("fine grid : y=%f x=%f max=%f\n", 727 * SpiralMeanCrossMaxYLocation 728 * ,SpiralMeanCrossMaxXLocation,SpiralMeanCrossMaxScore); 729 */ 730 731 SpiralCenterAnalysisField[0][1] = 0.0; 732 SpiralCenterAnalysisField[0][2] = 0.0; 733 734 /* determine lat/lon point from x/y coordinates */ 735 SpiralCenterLatitude = SpiralMeanCrossMaxYLocation + InputLatitude; 736 SpiralCenterLongitude = ((LongitudeMultFactor * SpiralMeanCrossMaxXLocation) / (Math.cos(PI 737 * InputLatitude / 180.0))) 738 + InputLongitude; 739 SpiralCenterScore = SpiralMeanCrossMaxScore; 740 for (YInc = 1; YInc <= (int) SpiralCenterAnalysisField[0][0]; YInc++) { 741 SpiralCenterAnalysisField[YInc][1] = SpiralCenterAnalysisField[YInc][1] + InputLatitude; 742 SpiralCenterAnalysisField[YInc][2] = ((LongitudeMultFactor * SpiralCenterAnalysisField[YInc][2]) / (Math 743 .cos(PI * InputLatitude / 180.0))) + InputLongitude; 744 } 745 746 return new double[] { SpiralCenterLatitude, SpiralCenterLongitude, SpiralCenterScore }; 747 } 748 749 private static void Gradient(double TemperatureInputArray[][], int ElementXNumber, 750 int LineYNumber, double LongitudeIncrement, double LatitudeIncrement) { 751 752 int XInc, YInc; 753 754 /* initialize arrays */ 755 for (YInc = 0; YInc < LineYNumber - 1; YInc++) { 756 for (XInc = 0; XInc < ElementXNumber - 1; XInc++) { 757 NSTempGradientArray[YInc][XInc] = 0.0; 758 EWTempGradientArray[YInc][XInc] = 0.0; 759 } 760 } 761 762 for (YInc = 1; YInc < LineYNumber - 1; YInc++) { 763 for (XInc = 1; XInc < ElementXNumber - 1; XInc++) { 764 /* determine N-S gradient at point */ 765 NSTempGradientArray[YInc][XInc] = (TemperatureInputArray[YInc - 1][XInc] - TemperatureInputArray[YInc + 1][XInc]) 766 / (2.0 * LatitudeIncrement); 767 /* determine E-W gradient at point */ 768 EWTempGradientArray[YInc][XInc] = (TemperatureInputArray[YInc][XInc + 1] - TemperatureInputArray[YInc][XInc - 1]) 769 / (2.0 * LongitudeIncrement); 770 } 771 } 772 773 } 774 775 private static double[] RingFit(double RingFitFirstGuessLatitude, 776 double RingFitFirstGuessLongitude) { 777 778 int XInc, YInc, ZInc; 779 int CircleFilterRadii; 780 int CircleFilterPoints; 781 782 int RingWidthLocal = (int) RING_WIDTH; 783 int XLocation = 0; 784 int YLocation = 0; 785 int CircleFilterXLocation = 0; 786 int CircleFilterYLocation = 0; 787 788 double[][] GradientFieldXDirection = new double[200][200]; 789 double[][] GradientFieldYDirection = new double[200][200]; 790 791 double RingFitSearchRadiusDegree = 0.75; 792 double RingFitMinRadiusDegree = 0.06; 793 double RingFitMaxRadiusDegree = 0.40; 794 double CircleFilterDistance = 0.0; 795 796 double DotScoreMaximum = -99999.0; 797 double DotProductXDirection = -999.9; 798 double DotProductYDirection = -999.9; 799 double DotProductXYValue = -999.9; 800 double SignFactor = 0.0; 801 double DotScoreValue = 0.0; 802 double DotScoreValueSum = 0.0; 803 double DotScoreFinal = 0.0; 804 double NANAdjustmentValue = 0.0; 805 806 int ImageResolution = (int) Data.GetCurrentImageResolution(); 807 int IncAddVal = (ImageResolution > RingWidthLocal) ? 1 808 : (RingWidthLocal - ImageResolution + 1); 809 810 /* derive values */ 811 double DegreesPerPixel = Math 812 .abs(IRData_Remap_Latitude[0][0] - IRData_Remap_Latitude[1][0]); 813 int RingFitSearchRadiusPixel = (int) Math 814 .round(RingFitSearchRadiusDegree / DegreesPerPixel); 815 int RingFitMinRadiusPix = (int) Math.max(2.0, 816 Math.round(RingFitMinRadiusDegree / DegreesPerPixel)); 817 int MaximumRadiusSizePixels = (int) Math.round(RingFitMaxRadiusDegree / DegreesPerPixel); 818 819 double LongitudeIncrement = 1.0; 820 double LatitudeIncrement = -1.0; 821 822 /* 823 * System.out.printf("RingFit: %d %f %d %d %d\n",IncAddVal,DegreesPerPixel 824 * , RingFitSearchRadiusPixel,RingFitMinRadiusPix, 825 * MaximumRadiusSizePixels); 826 */ 827 /* 828 * for(YInc=0;YInc<IRData_Remap_NumberRows-1;YInc++) { 829 * for(XInc=0;XInc<IRData_Remap_NumberColumns-1;XInc++) { 830 * System.out.printf 831 * ("Y=%d X=%d lat=%f lon=%f\n",YInc,XInc,IRData_Remap_Latitude 832 * [YInc][XInc],IRData_Remap_Longitude[YInc][XInc]); } } 833 */ 834 835 /* NEED TO PASS BACK BOTH 2D ARRAYS FROM GRADIENT AND CIRCLEFILT */ 836 /* calculate gradient field */ 837 Gradient(IRData_Remap_Temperature, IRData_Remap_NumberColumns, IRData_Remap_NumberRows, 838 LongitudeIncrement, LatitudeIncrement); 839 840 for (YInc = 0; YInc < IRData_Remap_NumberRows - 1; YInc++) { 841 for (XInc = 0; XInc < IRData_Remap_NumberColumns - 1; XInc++) { 842 GradientFieldYDirection[YInc][XInc] = NSTempGradientArray[YInc][XInc]; 843 GradientFieldXDirection[YInc][XInc] = EWTempGradientArray[YInc][XInc]; 844 /* 845 * System.out.printf("GRADIENT Y=%d X=%d Xdir=%f Ydir=%f ",YInc 846 * , 847 * XInc,GradientFieldXDirection[YInc][XInc],GradientFieldYDirection 848 * [YInc][XInc]); 849 */ 850 /* 851 * System.out.printf(" lat=%f lon=%f\n",IRData_Remap_Latitude[YInc 852 * ][XInc],IRData_Remap_Longitude[YInc][XInc]); 853 */ 854 } 855 } 856 /* 857 * System.out.printf( 858 * "RingFit : RingFitFirstGuessLatitude=%f RingFitFirstGuessLongitude=%f\n" 859 * , RingFitFirstGuessLatitude,RingFitFirstGuessLongitude); 860 */ 861 862 /* make matricies of row and column numbers */ 863 int Lalo2IndsReturn[] = Lalo2IndsFloat(RingFitFirstGuessLatitude, 864 RingFitFirstGuessLongitude, IRData_Remap_Latitude, IRData_Remap_Longitude, 865 IRData_Remap_NumberColumns, IRData_Remap_NumberRows); 866 867 int FirstGuessXLocation = Lalo2IndsReturn[0]; 868 int FirstGuessYLocation = Lalo2IndsReturn[1]; 869 870 /* 871 * System.out.printf( 872 * "RingFit : FirstGuessXLocation=%d FirstGuessYLocation=%d\n", 873 * FirstGuessXLocation,FirstGuessYLocation); 874 */ 875 876 /* initialize circle/ring filter arrays */ 877 /* radius in pixels */ 878 879 for (CircleFilterRadii = 0; CircleFilterRadii < 100; CircleFilterRadii++) { 880 /* number of points on circle at radius */ 881 for (CircleFilterPoints = 0; CircleFilterPoints < 2000; CircleFilterPoints++) { 882 CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints] = 0; 883 CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints] = 0; 884 } 885 } 886 887 /* 888 * determine digital pixel coordinates for ring analysis for different 889 * radii sizes 890 */ 891 for (CircleFilterRadii = RingFitMinRadiusPix; CircleFilterRadii <= MaximumRadiusSizePixels; CircleFilterRadii++) { 892 CircleFilt(CircleFilterRadii); 893 } 894 895 /* search image box */ 896 ZInc = 1; 897 /* develop the accumulator */ 898 for (CircleFilterRadii = RingFitMinRadiusPix; CircleFilterRadii < MaximumRadiusSizePixels; CircleFilterRadii++) { 899 int CircleFilterPointCount = CircleFilterRowArray[CircleFilterRadii][0]; 900 /* determine each main point in analysis disc */ 901 for (XInc = 1; XInc < IRData_Remap_NumberColumns - 1; XInc = XInc + IncAddVal) { 902 for (YInc = 1; YInc < IRData_Remap_NumberRows - 1; YInc = YInc + IncAddVal) { 903 /* 904 * System.out.printf("XInc=%d YInc=%d %d %d: %d\n",XInc,YInc 905 * , FirstGuessXLocation,FirstGuessYLocation, 906 * RingFitSearchRadiusPixel); 907 */ 908 909 CircleFilterDistance = (double) ((XInc - FirstGuessXLocation) * (XInc - FirstGuessXLocation)) 910 + ((YInc - FirstGuessYLocation) * (YInc - FirstGuessYLocation)); 911 /* 912 * System.out.printf( 913 * "XInc=%d YInc=%d CircleFilterDistance=%f RingFitSearchRadiusPixel=%d\n" 914 * , 915 * XInc,YInc,CircleFilterDistance,RingFitSearchRadiusPixel); 916 */ 917 if (CircleFilterDistance <= (double) (RingFitSearchRadiusPixel * RingFitSearchRadiusPixel)) { 918 /* 919 * if main point (YInc,XInc) is in disc, calculate 920 * dotproduct for each subpoint on ring around main 921 * point 922 */ 923 DotScoreValueSum = 0.0; 924 int NANCount = 0; 925 boolean FoundMoatRegionTF = false; 926 /* 927 * System.out.printf( 928 * "XInc=%d YInc=%d CircleFilterPointCount=%d %d\n", 929 * XInc,YInc,CircleFilterPointCount,CircleFilterRadii); 930 */ 931 for (CircleFilterPoints = 1; CircleFilterPoints <= CircleFilterPointCount; CircleFilterPoints++) { 932 /* 933 * System.out.printf( 934 * "CircleFilterPoints=%d CircleFilterRadii=%d XInc=%d YInc=%d\n" 935 * , 936 * CircleFilterPoints,CircleFilterRadii,XInc,YInc); 937 */ 938 CircleFilterXLocation = XInc 939 + CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints]; 940 CircleFilterYLocation = YInc 941 + CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints]; 942 /* 943 * System.out.printf("%d : %d %d -- %d : %d %d ", 944 * CircleFilterYLocation,YInc, 945 * CircleFilterRowArray[CircleFilterRadii 946 * ][CircleFilterPoints], 947 * CircleFilterXLocation,XInc, 948 * CircleFilterColumnArray 949 * [CircleFilterRadii][CircleFilterPoints]); 950 */ 951 if ((CircleFilterXLocation < 1) || (CircleFilterYLocation < 1)) { 952 DotProductXDirection = -999.9; 953 DotProductYDirection = -999.9; 954 } else { 955 /* 956 * System.out.printf( 957 * "I=%d x=%d J=%d y=%d : %b %b\n", 958 * CircleFilterXLocation 959 * ,IRData_Remap_NumberColumns, 960 * CircleFilterYLocation 961 * ,IRData_Remap_NumberRows, 962 * (CircleFilterXLocation 963 * >=IRData_Remap_NumberColumns), 964 * (CircleFilterYLocation 965 * >=IRData_Remap_NumberRows)); 966 */ 967 if ((CircleFilterXLocation >= (IRData_Remap_NumberColumns - 1)) 968 || (CircleFilterYLocation >= (IRData_Remap_NumberRows - 1))) { 969 DotProductXDirection = -999.9; 970 DotProductYDirection = -999.9; 971 } else { 972 if (MoatMaskFlagField[CircleFilterYLocation][CircleFilterXLocation] == 1) { 973 FoundMoatRegionTF = true; 974 } 975 /* 976 * System.out.printf( 977 * " mask=%d FoundMoatRegionTF=%b\n" 978 * ,MoatMaskFlagField 979 * [CircleFilterYLocation][ 980 * CircleFilterXLocation 981 * ],FoundMoatRegionTF); 982 */ 983 if (FoundMoatRegionTF) { 984 DotProductXDirection = -999.9; 985 DotProductYDirection = -999.9; 986 } else { 987 DotProductXDirection = ((double) CircleFilterRowArray[CircleFilterRadii][CircleFilterPoints] / (double) CircleFilterRadii) 988 * GradientFieldYDirection[CircleFilterYLocation][CircleFilterXLocation]; 989 /* 990 * System.out.printf( 991 * " XX | %d %d j=%d i=%d j=%d i=%d GradientFieldYDirection=%f GradientFieldXDirection=%f | \n" 992 * , 993 * CircleFilterPoints,CircleFilterRadii, 994 * CircleFilterRowArray 995 * [CircleFilterRadii 996 * ][CircleFilterPoints], 997 * CircleFilterColumnArray 998 * [CircleFilterRadii 999 * ][CircleFilterPoints], 1000 * CircleFilterYLocation, 1001 * CircleFilterXLocation, 1002 * GradientFieldYDirection 1003 * [CircleFilterYLocation 1004 * ][CircleFilterXLocation], 1005 * GradientFieldXDirection 1006 * [CircleFilterYLocation 1007 * ][CircleFilterXLocation]); 1008 */ 1009 DotProductYDirection = ((double) CircleFilterColumnArray[CircleFilterRadii][CircleFilterPoints] / (double) CircleFilterRadii) 1010 * GradientFieldXDirection[CircleFilterYLocation][CircleFilterXLocation]; 1011 } 1012 } 1013 } 1014 if ((DotProductXDirection < -999.0) || (DotProductYDirection < -999.0)) { 1015 NANCount++; 1016 } else { 1017 DotProductXYValue = DotProductXDirection + DotProductYDirection; 1018 if (DotProductXYValue == 0.0) { 1019 SignFactor = 0.0; 1020 } else { 1021 /* return -1/+1 for -/+ value */ 1022 SignFactor = Math.abs(DotProductXYValue) / DotProductXYValue; 1023 } 1024 DotScoreValue = SignFactor 1025 * (Math.log(1.0 + Math.abs(DotProductXYValue))); 1026 DotScoreValueSum = DotScoreValueSum + DotScoreValue; 1027 } 1028 /* 1029 * System.out.printf( 1030 * "dot product X=%f Y=%f XY=%f value=%f sum=%f\n" 1031 * ,DotProductXDirection, 1032 * DotProductYDirection,DotProductXYValue 1033 * ,DotScoreValue,DotScoreValueSum); 1034 */ 1035 } /* if indisk */ 1036 /* 1037 * System.out.printf("dot product Final Sum=%f \n", 1038 * DotScoreValueSum); 1039 */ 1040 /* 1041 * check for missing data and adjust DotScoreFinal 1042 * accordingly 1043 */ 1044 /* 1045 * System.out.printf( 1046 * "Y=%d X=%d foundmoat=%b nancount=%f xyz=%f\n" 1047 * ,YInc,XInc 1048 * ,FoundMoatRegionTF,(double)NANCount,0.575*(double 1049 * )CircleFilterPointCount); 1050 */ 1051 if (FoundMoatRegionTF 1052 || ((double) NANCount) > (0.575 * (double) CircleFilterPointCount)) { 1053 DotScoreFinal = 0.0; 1054 } else { 1055 NANAdjustmentValue = (double) CircleFilterPointCount 1056 / (double) (CircleFilterPointCount - NANCount); 1057 DotScoreFinal = -NANAdjustmentValue * DotScoreValueSum 1058 / Math.sqrt((double) CircleFilterPointCount); 1059 } 1060 /* 1061 * System.out.printf("XX final=%f maximum=%f ",DotScoreFinal 1062 * ,DotScoreMaximum); 1063 */ 1064 /* 1065 * System.out.printf( 1066 * " circlefilterpointcount=%d sqrt=%f " 1067 * ,CircleFilterPointCount 1068 * ,Math.sqrt((double)CircleFilterPointCount)); 1069 */ 1070 if (DotScoreFinal > DotScoreMaximum) { 1071 DotScoreMaximum = DotScoreFinal; 1072 XLocation = XInc; 1073 YLocation = YInc; 1074 /* 1075 * System.out.printf(" xloc=%d yloc=%d MaxRadSize=%d" 1076 * ,XInc,YInc,MaximumRadiusSize); 1077 */ 1078 } 1079 /* System.out.printf("\n"); */ 1080 RingScoreAnalysisField[ZInc][0] = DotScoreFinal; 1081 RingScoreAnalysisField[ZInc][1] = IRData_Remap_Latitude[YInc][XInc]; 1082 RingScoreAnalysisField[ZInc][2] = IRData_Remap_Longitude[YInc][XInc]; 1083 RingScoreAnalysisField[0][0] = (double) ZInc; 1084 /* 1085 * System.out.printf("ZZZZ 0=%f 1=%f 2=%f x=%f\n", 1086 * RingScoreAnalysisField 1087 * [ZInc][0],RingScoreAnalysisField[ZInc][1], 1088 * RingScoreAnalysisField 1089 * [ZInc][2],RingScoreAnalysisField[0][0]); 1090 */ 1091 ZInc++; 1092 } 1093 } /* XInc */ 1094 } /* YInc */ 1095 } /* CircleFilterRadii */ 1096 1097 RingScoreAnalysisField[0][1] = 0.0; 1098 RingScoreAnalysisField[0][2] = 0.0; 1099 1100 /* make matricies of row and column numbers */ 1101 /* System.out.printf("inds2lalo: x=%d y=%d \n",XLocation,YLocation); */ 1102 double Inds2LaloReturn[] = Inds2LaloFloat(XLocation, YLocation, IRData_Remap_Latitude, 1103 IRData_Remap_Longitude, IRData_Remap_NumberColumns, IRData_Remap_NumberRows); 1104 double MaxScoreLatitude = Inds2LaloReturn[0]; 1105 double MaxScoreLongitude = Inds2LaloReturn[1]; 1106 1107 double RingFitLatitude = MaxScoreLatitude; 1108 double RingFitLongitude = MaxScoreLongitude; 1109 double RingFitMaxScore = DotScoreMaximum; 1110 1111 return new double[] { RingFitLatitude, RingFitLongitude, RingFitMaxScore }; 1112 } 1113 1114 private static int[] Lalo2IndsFloat(double LatitudeInput, double LongitudeInput, 1115 double LatitudeArrayInput[][], double LongitudeArrayInput[][], int ElementXNumber, 1116 int LineYNumber) { 1117 1118 double LatitudeMinimum = LatitudeArrayInput[LineYNumber - 1][0]; 1119 double LatitudeMaximum = LatitudeArrayInput[0][0]; 1120 double LongitudeMinimum = LongitudeArrayInput[0][ElementXNumber - 1]; 1121 double LongitudeMaximum = LongitudeArrayInput[0][0]; 1122 1123 int YAxisPositionReturn = (int) ((((double) LineYNumber - 1.0) / (LatitudeMaximum - LatitudeMinimum)) * (LatitudeMaximum - LatitudeInput)); 1124 int XAxisPositionReturn = (int) ((((double) ElementXNumber - 1.0) / (LongitudeMaximum - LongitudeMinimum)) * (LongitudeMaximum - LongitudeInput)); 1125 1126 return new int[] { XAxisPositionReturn, YAxisPositionReturn }; 1127 1128 } 1129 1130 private static double[] Inds2LaloFloat(int XAxisPosition, int YAxisPosition, 1131 double LatitudeArrayInput[][], double LongitudeArrayInput[][], int ElementXNumber, 1132 int LineYNumber) { 1133 1134 double LatitudeMinimum = LatitudeArrayInput[LineYNumber - 1][0]; 1135 double LatitudeMaximum = LatitudeArrayInput[0][0]; 1136 double LongitudeMinimum = LongitudeArrayInput[0][ElementXNumber - 1]; 1137 double LongitudeMaximum = LongitudeArrayInput[0][0]; 1138 1139 double LongitudeReturn = LongitudeMaximum 1140 - ((((double) XAxisPosition) / ((double) ElementXNumber - 1.0)) * (LongitudeMaximum - LongitudeMinimum)); 1141 double LatitudeReturn = LatitudeMaximum 1142 - ((((double) YAxisPosition) / ((double) LineYNumber - 1.0)) * (LatitudeMaximum - LatitudeMinimum)); 1143 1144 return new double[] { LatitudeReturn, LongitudeReturn }; 1145 1146 } 1147 1148 private static void CircleFilt(int RingRadiusInput) { 1149 1150 int XInc, YInc; 1151 int PointCount = 1; 1152 int RadiusPlus1 = RingRadiusInput + 1; 1153 1154 double ThresholdDifference = 0.5 * (double) (((RingRadiusInput + 1) * (RingRadiusInput + 1)) - (RingRadiusInput * RingRadiusInput)); 1155 for (XInc = -RadiusPlus1; XInc <= RadiusPlus1; XInc++) { 1156 for (YInc = -RadiusPlus1; YInc <= RadiusPlus1; YInc++) { 1157 double DifferenceValue = (double) ((YInc * YInc) + (XInc * XInc) - (RingRadiusInput * RingRadiusInput)); 1158 if (Math.abs(DifferenceValue) <= ThresholdDifference) { 1159 CircleFilterRowArray[RingRadiusInput][PointCount] = YInc; 1160 CircleFilterColumnArray[RingRadiusInput][PointCount] = XInc; 1161 PointCount++; 1162 } 1163 } 1164 } 1165 1166 /* number of points on given radius size */ 1167 CircleFilterRowArray[RingRadiusInput][0] = PointCount - 1; 1168 CircleFilterColumnArray[RingRadiusInput][0] = PointCount - 1; 1169 1170 } 1171 1172 private static void MoatMaskCalc(double MoatMaskTempThreshold, double RingFitMaxRadiusDegree, 1173 int MoatSignCheckFlag) { 1174 1175 int ArrayInc, XInc, YInc; 1176 int NumberColumns = IRData_Remap_NumberColumns; 1177 int NumberRows = IRData_Remap_NumberRows; 1178 1179 /* initialize arrays */ 1180 for (YInc = 0; YInc < NumberRows; YInc++) { 1181 for (XInc = 0; XInc < NumberColumns; XInc++) { 1182 MoatMaskFlagField[YInc][XInc] = 0; 1183 BlackWhiteFieldArray[YInc][XInc] = 0; 1184 } 1185 } 1186 1187 MeshGrid(MoatMaskTempThreshold, MoatSignCheckFlag); 1188 int NumberOfObjects = BWImage(8, NumberRows, NumberColumns); 1189 1190 for (ArrayInc = 1; ArrayInc <= NumberOfObjects; ArrayInc++) { 1191 double LatitudeMaximum = -90.0; 1192 double LatitudeMinimum = 90.0; 1193 double LongitudeMaximum = -180.0; 1194 double LongitudeMinimum = 180.0; 1195 for (YInc = 0; YInc < NumberRows; YInc++) { 1196 for (XInc = 0; XInc < NumberColumns; XInc++) { 1197 if (MoatMaskFlagField[YInc][XInc] == ArrayInc) { 1198 double LatitudeValue = IRData_Remap_Latitude[YInc][XInc]; 1199 double LongitudeValue = IRData_Remap_Longitude[YInc][XInc]; 1200 if (LatitudeValue > LatitudeMaximum) { 1201 LatitudeMaximum = LatitudeValue; 1202 } 1203 if (LongitudeValue > LongitudeMaximum) { 1204 LongitudeMaximum = LongitudeValue; 1205 } 1206 if (LatitudeValue < LatitudeMinimum) { 1207 LatitudeMinimum = LatitudeValue; 1208 } 1209 if (LongitudeValue < LongitudeMinimum) { 1210 LongitudeMinimum = LongitudeValue; 1211 } 1212 } 1213 } 1214 } 1215 double AverageLatitude = .5 * (LatitudeMinimum + LatitudeMaximum); 1216 double FeatureLength = Math.sqrt(Math.pow((LatitudeMaximum - LatitudeMinimum), 2) 1217 + Math.pow( 1218 (LongitudeMaximum - LongitudeMinimum) 1219 / Math.cos((PI * AverageLatitude) / 180.0), 2)); 1220 int MoatFlagID = (FeatureLength > (RingFitMaxRadiusDegree * 2.0)) ? 1 : 0; 1221 for (YInc = 0; YInc < NumberRows; YInc++) { 1222 for (XInc = 0; XInc < NumberColumns; XInc++) { 1223 if (MoatMaskFlagField[YInc][XInc] == ArrayInc) { 1224 MoatMaskFlagField[YInc][XInc] = MoatFlagID; 1225 } 1226 } 1227 } 1228 } 1229 1230 } 1231 1232 private static void MeshGrid(double TemperatureThresholdValue, int MoatSignCheckFlag) { 1233 1234 int XInc, YInc; 1235 1236 for (XInc = 0; XInc < IRData_Remap_NumberColumns; XInc++) { 1237 for (YInc = 0; YInc < IRData_Remap_NumberRows; YInc++) { 1238 if (MoatSignCheckFlag == 1) { 1239 BlackWhiteFieldArray[YInc][XInc] = (IRData_Remap_Temperature[YInc][XInc] > TemperatureThresholdValue) ? 1 1240 : 0; 1241 } else { 1242 BlackWhiteFieldArray[YInc][XInc] = (IRData_Remap_Temperature[YInc][XInc] < TemperatureThresholdValue) ? 1 1243 : 0; 1244 } 1245 } 1246 } 1247 1248 } 1249 1250 private static int BWImage(int ConnectednessValue, int NumberRows, int NumberColumns) { 1251 /* 1252 * blackwhitefieldarray = BooleanValueArray; moatmaskflagfield = 1253 * LabelImageArray 1254 */ 1255 1256 int IncVal, XInc, YInc; 1257 int B_Value, C_Value, D_Value, E_Value; 1258 int TableElements = 0; 1259 int LabelValue = 0; 1260 int[] LabelSetArray = new int[NumberRows * NumberColumns]; 1261 1262 LabelSetArray[0] = 0; 1263 1264 for (YInc = 0; YInc < NumberRows; YInc++) { 1265 for (XInc = 0; XInc < NumberColumns; XInc++) { 1266 if (BlackWhiteFieldArray[YInc][XInc] == 1) { 1267 /* if A is an object */ 1268 /* 1269 * get the neighboring pixels B_Value, C_Value, D_Value, and 1270 * E_Value 1271 */ 1272 B_Value = 0; 1273 C_Value = 0; 1274 D_Value = 0; 1275 E_Value = 0; 1276 if (XInc > 0) { 1277 B_Value = Find(LabelSetArray, MoatMaskFlagField[YInc][XInc - 1]); 1278 } 1279 if (YInc > 0) { 1280 C_Value = Find(LabelSetArray, MoatMaskFlagField[YInc - 1][XInc]); 1281 } 1282 if ((YInc > 0) && (XInc > 0)) { 1283 D_Value = Find(LabelSetArray, MoatMaskFlagField[YInc - 1][XInc - 1]); 1284 } 1285 if ((YInc > 0) && (XInc > (NumberColumns - 1))) { 1286 E_Value = Find(LabelSetArray, MoatMaskFlagField[YInc - 1][XInc + 1]); 1287 } 1288 if (ConnectednessValue == 4) { 1289 /* apply 4 connectedness */ 1290 if ((B_Value != 0) && (C_Value != 0)) { 1291 /* B_Value and C_Value are labeled */ 1292 if (B_Value == C_Value) { 1293 MoatMaskFlagField[YInc][XInc] = B_Value; 1294 } else { 1295 LabelSetArray[C_Value] = B_Value; 1296 MoatMaskFlagField[YInc][XInc] = B_Value; 1297 } 1298 } else if (B_Value != 0) { 1299 /* B_Value is object but C_Value is not */ 1300 MoatMaskFlagField[YInc][XInc] = B_Value; 1301 } else if (C_Value != 0) { 1302 /* C_Value is object but B_Value is not */ 1303 MoatMaskFlagField[YInc][XInc] = C_Value; 1304 } else { 1305 /* B_Value, C_Value, D_Value not object - new object */ 1306 /* label and put into table */ 1307 TableElements++; 1308 LabelSetArray[TableElements] = TableElements; 1309 MoatMaskFlagField[YInc][XInc] = LabelSetArray[TableElements]; 1310 } 1311 } else if (ConnectednessValue == 6) { 1312 /* apply 6 connected ness */ 1313 if (D_Value != 0) { 1314 /* D_Value object, copy label and move on */ 1315 MoatMaskFlagField[YInc][XInc] = D_Value; 1316 } else if ((B_Value != 0) && (C_Value != 0)) { 1317 /* B_Value and C_Value are labeled */ 1318 if (B_Value == C_Value) { 1319 MoatMaskFlagField[YInc][XInc] = B_Value; 1320 } else { 1321 LabelValue = Math.min(B_Value, C_Value); 1322 LabelSetArray[B_Value] = LabelValue; 1323 LabelSetArray[C_Value] = LabelValue; 1324 MoatMaskFlagField[YInc][XInc] = LabelValue; 1325 } 1326 } else if (B_Value != 0) { 1327 /* B_Value is object but C_Value is not */ 1328 MoatMaskFlagField[YInc][XInc] = B_Value; 1329 } else if (C_Value != 0) { 1330 /* C_Value is object but B_Value is not */ 1331 MoatMaskFlagField[YInc][XInc] = C_Value; 1332 } else { 1333 /* B_Value, C_Value, D_Value not object - new object */ 1334 /* label and put into table */ 1335 TableElements++; 1336 LabelSetArray[TableElements] = TableElements; 1337 MoatMaskFlagField[YInc][XInc] = LabelSetArray[TableElements]; 1338 } 1339 } else if (ConnectednessValue == 8) { 1340 /* apply 8 connectedness */ 1341 if ((B_Value != 0) || (C_Value != 0) || (D_Value != 0) || (E_Value != 0)) { 1342 LabelValue = B_Value; 1343 if (B_Value != 0) { 1344 LabelValue = B_Value; 1345 } else if (C_Value != 0) { 1346 LabelValue = C_Value; 1347 } else if (D_Value != 0) { 1348 LabelValue = D_Value; 1349 } else { 1350 LabelValue = E_Value; 1351 } 1352 MoatMaskFlagField[YInc][XInc] = LabelValue; 1353 if ((B_Value != 0) && (B_Value != LabelValue)) { 1354 LabelSetArray[B_Value] = LabelValue; 1355 } 1356 if ((C_Value != 0) && (C_Value != LabelValue)) { 1357 LabelSetArray[C_Value] = LabelValue; 1358 } 1359 if ((D_Value != 0) && (D_Value != LabelValue)) { 1360 LabelSetArray[D_Value] = LabelValue; 1361 } 1362 if ((E_Value != 0) && (E_Value != LabelValue)) { 1363 LabelSetArray[E_Value] = LabelValue; 1364 } 1365 } else { 1366 /* label and put into table */ 1367 TableElements++; 1368 LabelSetArray[TableElements] = TableElements; 1369 MoatMaskFlagField[YInc][XInc] = LabelSetArray[TableElements]; 1370 } 1371 } 1372 } else { 1373 /* A is not an object so leave it */ 1374 MoatMaskFlagField[YInc][XInc] = 0; 1375 } 1376 } 1377 } 1378 1379 /* consolidate component table */ 1380 for (IncVal = 0; IncVal <= TableElements; IncVal++) { 1381 LabelSetArray[IncVal] = Find(LabelSetArray, IncVal); 1382 } 1383 1384 /* run image through the look-up table */ 1385 for (YInc = 0; YInc < NumberRows; YInc++) { 1386 for (XInc = 0; XInc < NumberColumns; XInc++) { 1387 MoatMaskFlagField[YInc][XInc] = LabelSetArray[MoatMaskFlagField[YInc][XInc]]; 1388 } 1389 } 1390 1391 /* count up the objects in the image */ 1392 for (IncVal = 0; IncVal <= TableElements; IncVal++) { 1393 LabelSetArray[IncVal] = 0; 1394 } 1395 1396 for (YInc = 0; YInc < NumberRows; YInc++) { 1397 for (XInc = 0; XInc < NumberColumns; XInc++) { 1398 LabelSetArray[MoatMaskFlagField[YInc][XInc]]++; 1399 } 1400 } 1401 1402 /* number the objects from 1 through n objects */ 1403 int NumberObjects = 0; 1404 LabelSetArray[0] = 0; 1405 for (IncVal = 1; IncVal <= TableElements; IncVal++) { 1406 if (LabelSetArray[IncVal] > 0) { 1407 LabelSetArray[IncVal] = ++NumberObjects; 1408 } 1409 } 1410 1411 /* run through the look-up table again */ 1412 for (YInc = 0; YInc < NumberRows; YInc++) { 1413 for (XInc = 0; XInc < NumberColumns; XInc++) { 1414 MoatMaskFlagField[YInc][XInc] = LabelSetArray[MoatMaskFlagField[YInc][XInc]]; 1415 } 1416 } 1417 1418 return NumberObjects; 1419 1420 } 1421 1422 private static int Find(int InputArray[], int InputValue) { 1423 int IncVal = InputValue; 1424 1425 while (InputArray[IncVal] != IncVal) { 1426 IncVal = InputArray[IncVal]; 1427 } 1428 1429 return IncVal; 1430 } 1431 1432 /** 1433 * This routine will determine the confidence scores for the spiral fitting 1434 * and ring fitting routines and calculate the best possible position for 1435 * the storm center based upon a combination of the two methods, if 1436 * available. 1437 * 1438 * If the ring fitting routine does not determine a good candidate position, 1439 * the spiral fitting routine will be used alone. If the spiral fitting 1440 * routine candidate point is not "good", the forecast point will be 1441 * selected. 1442 * 1443 * @param FirstGuessLatitude 1444 * First Guess latitude. 1445 * @param FirstGuessLongitude 1446 * First Guess longitude. 1447 * @param SpiralCenterLatitude 1448 * Spiral Analysis latitude at max location. 1449 * @param SpiralCenterLongitude 1450 * Spiral Analysis longitude at max location. 1451 * @param SpiralCenterScoreValue 1452 * Spiral Analysis Score value. 1453 * @param RingFitLatitude 1454 * Ring Analysis latitude at max score location. 1455 * @param RingFitLongitude 1456 * Ring Analysis longitude at max score location. 1457 * @param RingFitScoreValue 1458 * Ring Analysis Score value. 1459 * 1460 * @return Array of four double values. The values represent: latitude of 1461 * final selected location, longitude of final selected location, 1462 * confidence score of final selected location, and the 1463 * {@literal "method"} used to determine the final selected 1464 * location. Possible values for the method: 1 for first guess, 2 1465 * for enhanced spiral analysis, and 5 for combo ring/spiral 1466 * analysis. 1467 */ 1468 private static double[] CalcScores(double FirstGuessLatitude, double FirstGuessLongitude, 1469 double SpiralCenterLatitude, double SpiralCenterLongitude, 1470 double SpiralCenterScoreValue, double RingFitLatitude, double RingFitLongitude, 1471 double RingFitScoreValue) { 1472 int YInc; 1473 int FinalLocationSelectionMethodReturn = 0; 1474 int NumberOfSpirals = (int) SpiralCenterAnalysisField[0][0] + 1; 1475 1476 double[] DistancePenaltyArray = new double[NumberOfSpirals]; 1477 double[] InitialSpiralScoreArray = new double[NumberOfSpirals]; 1478 double[] DistanceFromSpiralCenterArray = new double[NumberOfSpirals]; 1479 double[] DistanceBonusArray = new double[NumberOfSpirals]; 1480 double[] EnhancedSpiralScoreArray = new double[NumberOfSpirals]; 1481 double[] FinalSpiralScoreArray = new double[NumberOfSpirals]; 1482 1483 double FinalLatitudeReturn = -99.99; 1484 double FinalLongitudeReturn = -999.99; 1485 double FinalScoreReturn = 0.0; 1486 double DistanceValue; 1487 double MaximumForecastErrorDegree = 1.0; /* 1488 * max dist between fcst and 1489 * final 1490 */ 1491 double ExpectedMaximumForecastErrorDegree = MaximumForecastErrorDegree; /* 1492 * expected 1493 * error 1494 */ 1495 1496 /* max displacement */ 1497 double MaximumAllowableDisplacement = ExpectedMaximumForecastErrorDegree * 1.15; 1498 1499 /* Spiral Center analysis weight */ 1500 double SPIRALWEIGHT = 10.0; 1501 1502 /* RF distance penalty */ 1503 double DISTPENALTYWEIGHT = (0.5 / ExpectedMaximumForecastErrorDegree); 1504 1505 /* Ring Fit bonus value */ 1506 double PROXIMITYBONUS = 4.5; 1507 1508 /* RF bonus value threshold dist deg */ 1509 double PROXIMITYTHRESH = 0.25; 1510 1511 /* combination score threshold value */ 1512 double COMBOSCORETHRESH = 15.0; 1513 1514 /* convert distance in km to degrees */ 1515 double KMperDegree = 111.0; 1516 1517 double SpiralCenterIndexMaximum = -999.99; 1518 double SpiralCenterMaximumLatitude = -99.99; 1519 double SpiralCenterMaximumLongitude = -999.99; 1520 double SpiralMaximumDistanceFromGuess = 999.99; 1521 double IntermediateRingScore = 0.0; 1522 double FinalSpiralScoreValue = 0.0; 1523 double MaximumSpiralScoreValue = -99.99; 1524 double MaximumSpiralScoreLatitude = -99.99; 1525 double MaximumSpiralScoreLongitude = -999.99; 1526 1527 /* Spiral Score Calculations */ 1528 double LocalValue[] = Functions.distance_angle(FirstGuessLatitude, FirstGuessLongitude, 1529 SpiralCenterLatitude, SpiralCenterLongitude, 1); 1530 DistanceValue = LocalValue[0]; 1531 double InitialSpiralScoreArrayScore = SpiralCenterScoreValue; 1532 1533 for (YInc = 1; YInc < NumberOfSpirals; YInc++) { 1534 double LocalValue1[] = Functions.distance_angle(FirstGuessLatitude, 1535 FirstGuessLongitude, SpiralCenterAnalysisField[YInc][1], 1536 SpiralCenterAnalysisField[YInc][2], 1); 1537 DistanceValue = LocalValue1[0]; 1538 DistancePenaltyArray[YInc] = -DISTPENALTYWEIGHT * (DistanceValue / KMperDegree); 1539 InitialSpiralScoreArray[YInc] = SPIRALWEIGHT 1540 * (SpiralCenterAnalysisField[YInc][0] - InitialSpiralScoreArrayScore); 1541 double LocalValue2[] = Functions.distance_angle(SpiralCenterLatitude, 1542 SpiralCenterLongitude, SpiralCenterAnalysisField[YInc][1], 1543 SpiralCenterAnalysisField[YInc][2], 1); 1544 DistanceValue = LocalValue2[0]; 1545 DistanceFromSpiralCenterArray[YInc] = DistanceValue / KMperDegree; 1546 if (DistanceFromSpiralCenterArray[YInc] <= PROXIMITYTHRESH) { 1547 DistanceBonusArray[YInc] = PROXIMITYBONUS; 1548 } else { 1549 DistanceBonusArray[YInc] = 0.0; 1550 } 1551 EnhancedSpiralScoreArray[YInc] = InitialSpiralScoreArray[YInc] 1552 + DistancePenaltyArray[YInc]; 1553 /* 1554 * System.out.printf( 1555 * "SpiralCenterAnalysisField : YInc: %d Lat: %f Lon: %f Value: %f " 1556 * ,YInc,SpiralCenterAnalysisField[YInc][1], 1557 * SpiralCenterAnalysisField 1558 * [YInc][2],SpiralCenterAnalysisField[YInc][0]); 1559 */ 1560 /* 1561 * System.out.printf( 1562 * " distPentaly=%f DistanceFromSpiralCenterArray=%f DistanceBonusArray=%f max=%f\n" 1563 * , DistancePenaltyArray[YInc],DistanceFromSpiralCenterArray[YInc], 1564 * DistanceBonusArray[YInc],SpiralCenterIndexMaximum); 1565 */ 1566 if (EnhancedSpiralScoreArray[YInc] > SpiralCenterIndexMaximum) { 1567 SpiralCenterIndexMaximum = EnhancedSpiralScoreArray[YInc]; 1568 SpiralCenterMaximumLatitude = SpiralCenterAnalysisField[YInc][1]; 1569 SpiralCenterMaximumLongitude = SpiralCenterAnalysisField[YInc][2]; 1570 } 1571 FinalSpiralScoreArray[YInc] = InitialSpiralScoreArray[YInc] 1572 + DistancePenaltyArray[YInc] + DistanceBonusArray[YInc]; 1573 } 1574 /* 1575 * System.out.printf( 1576 * "FirstGuessLatitude=%f FirstGuessLongitude=%f SpiralCenterMaximumLatitude=%f SpiralCenterMaximumLongitude=%f\n" 1577 * , FirstGuessLatitude,FirstGuessLongitude,SpiralCenterMaximumLatitude, 1578 * SpiralCenterMaximumLongitude); 1579 */ 1580 double LocalValue3[] = Functions.distance_angle(FirstGuessLatitude, FirstGuessLongitude, 1581 SpiralCenterMaximumLatitude, SpiralCenterMaximumLongitude, 1); 1582 DistanceValue = LocalValue3[0]; 1583 SpiralMaximumDistanceFromGuess = DistanceValue / KMperDegree; 1584 /* 1585 * System.out.printf( 1586 * "SpiralMaximumDistanceFromGuess=%f MaximumAllowableDisplacement=%f\n" 1587 * , SpiralMaximumDistanceFromGuess,MaximumAllowableDisplacement); 1588 */ 1589 if (SpiralMaximumDistanceFromGuess <= MaximumAllowableDisplacement) { 1590 /* Ring Score Calculations */ 1591 double LocalValue4[] = Functions.distance_angle(FirstGuessLatitude, 1592 FirstGuessLongitude, RingFitLatitude, RingFitLongitude, 1); 1593 DistanceValue = LocalValue4[0]; 1594 for (YInc = 1; YInc < NumberOfSpirals; YInc++) { 1595 double[] RingScoreReturn = FindRingScore(SpiralCenterAnalysisField[YInc][1], 1596 SpiralCenterAnalysisField[YInc][2], RingScoreAnalysisField); 1597 int RetErr = (int) RingScoreReturn[0]; 1598 IntermediateRingScore = RingScoreReturn[1]; 1599 if (RetErr == 1) { 1600 FinalSpiralScoreValue = FinalSpiralScoreArray[YInc] + IntermediateRingScore; 1601 /* 1602 * System.out.printf( 1603 * " FOUND NumberOfSpirals=%d lat=%f lon=%f ringScore=%f FinalSpiralScoreValue=%f MaximumSpiralScoreValue=%f\n" 1604 * , YInc,SpiralCenterAnalysisField[YInc][1], 1605 * SpiralCenterAnalysisField[YInc][2], 1606 * IntermediateRingScore, 1607 * FinalSpiralScoreValue,MaximumSpiralScoreValue); 1608 */ 1609 if (FinalSpiralScoreValue > MaximumSpiralScoreValue) { 1610 MaximumSpiralScoreValue = FinalSpiralScoreValue; 1611 MaximumSpiralScoreLatitude = SpiralCenterAnalysisField[YInc][1]; 1612 MaximumSpiralScoreLongitude = SpiralCenterAnalysisField[YInc][2]; 1613 } 1614 } 1615 } 1616 /* 1617 * System.out.printf( 1618 * "MaximumSpiralScoreValue=%f RingPart=%f SpiralCenterIndexMaximum=%f\n" 1619 * , MaximumSpiralScoreValue,IntermediateRingScore, 1620 * SpiralCenterIndexMaximum); 1621 */ 1622 if (MaximumSpiralScoreValue >= COMBOSCORETHRESH) { 1623 /* use Combo Method */ 1624 FinalScoreReturn = MaximumSpiralScoreValue; 1625 FinalLatitudeReturn = MaximumSpiralScoreLatitude; 1626 FinalLongitudeReturn = MaximumSpiralScoreLongitude; 1627 FinalLocationSelectionMethodReturn = 5; 1628 } else { 1629 /* use ESP method */ 1630 FinalScoreReturn = 1.0 + SpiralCenterIndexMaximum; 1631 FinalLatitudeReturn = SpiralCenterMaximumLatitude; 1632 FinalLongitudeReturn = SpiralCenterMaximumLongitude; 1633 FinalLocationSelectionMethodReturn = 4; 1634 } 1635 } else { 1636 /* use Forecast Position */ 1637 FinalScoreReturn = 0.0; 1638 FinalLatitudeReturn = FirstGuessLatitude; 1639 FinalLongitudeReturn = FirstGuessLongitude; 1640 FinalLocationSelectionMethodReturn = 1; 1641 } 1642 1643 return new double[] { FinalLatitudeReturn, FinalLongitudeReturn, FinalScoreReturn, 1644 (double) FinalLocationSelectionMethodReturn }; 1645 } 1646 1647 /** 1648 * Find Ring score at selected location (spiral analysis location) 1649 * 1650 * @param FirstGuessLatitude 1651 * Latitude of search location. 1652 * @param FirstGuessLongitude 1653 * Longitude of search location. 1654 * @param RingScoreAnalysisField 1655 * - Array/Grid of Ring Analysis scores. 1656 * 1657 * @return Array of two doubles. The first value will be either -1 (if 1658 * nothing was found) or 1 (if found). The second value is the value 1659 * at search location in ring analysis grid (if the first value is 1660 * 1). 1661 */ 1662 private static double[] FindRingScore(double FirstGuessLatitude, double FirstGuessLongitude, 1663 double RingScoreAnalysisField[][]) { 1664 int YInc = 1; 1665 int Ret_ID = -1; 1666 double MaxRingDistance = RingScoreAnalysisField[0][0]; 1667 double RingFixScoreReturn = -99.9; 1668 1669 double LatitudeIncrement = Math.abs(IRData_Remap_Latitude[0][0] 1670 - IRData_Remap_Latitude[1][0]) / 2.0; 1671 double LongitudeIncrement = Math.abs(IRData_Remap_Longitude[0][1] 1672 - IRData_Remap_Longitude[0][0]) / 2.0; 1673 while (YInc <= MaxRingDistance) { 1674 if ((Math.abs(RingScoreAnalysisField[YInc][1] - FirstGuessLatitude) <= LatitudeIncrement) 1675 && (Math.abs(RingScoreAnalysisField[YInc][2] - FirstGuessLongitude) <= LongitudeIncrement)) { 1676 if (RingScoreAnalysisField[YInc][0] > RingFixScoreReturn) { 1677 RingFixScoreReturn = RingScoreAnalysisField[YInc][0]; 1678 } 1679 Ret_ID = 1; 1680 } 1681 YInc++; 1682 } 1683 1684 return new double[] { (double) Ret_ID, RingFixScoreReturn }; 1685 } 1686 1687 /** 1688 * Determine method location scheme to use by examining various 1689 * empirically-defined confidence factors. 1690 * 1691 * Confidence factors will be derived, with the "most confident" value used 1692 * as the final automatically determined storm position. 1693 * 1694 * @param ForecastLatitude 1695 * NHC/JTWC interpolated latitude position. 1696 * 1697 * @param ForecastLongitude 1698 * NHC/JTWC interpolated longitude position. 1699 * 1700 * @param RingSpiralLatitude 1701 * Ring/Spiral Analysis latitude position. 1702 * 1703 * @param RingSpiralLongitude 1704 * Ring/Spiral Analysis longitude position. 1705 * 1706 * @param RingSpiralScore 1707 * Ring/Spiral Analysis confidence factor score. 1708 * 1709 * @param RingSpiralSelectionIDValue 1710 * Ring/Spiral Analysis position derivation method. 1711 * 1712 * @return Array of three doubles. First value is latitude to be used, 1713 * second value is longitude to be used, and the third value is the 1714 * method used to determine storm position values. 1715 */ 1716 private static double[] PickFinalLocation(int InputPositioningID, double ForecastLatitude, 1717 double ForecastLongitude, double RingSpiralLatitude, double RingSpiralLongitude, 1718 double RingSpiralScore, int RingSpiralSelectionIDValue) { 1719 1720 int PositioningMethodID = 0; 1721 int HistoryRecEyeScene; 1722 int HistoryRecCloudScene; 1723 double FinalLatitude = -99.99; 1724 double FinalLongitude = -99.99; 1725 double FinalScoreValue = 0.0; 1726 double HistoryRecFinalTno = 0.0; 1727 double MaximumHistoryRecCI = 0.0; 1728 double HistoryRecRawTno = 0.0; 1729 double HistoryRecCI = 0.0; 1730 boolean ForceAutoFixUseTF_forTesting = false; 1731 boolean FoundEyeSceneTF = false; 1732 1733 int CurDate = History.IRCurrentRecord.date; 1734 int CurTime = History.IRCurrentRecord.time; 1735 int CurCloudScene = History.IRCurrentRecord.cloudscene; 1736 int CurEyeScene = History.IRCurrentRecord.eyescene; 1737 double CurrentTime = Functions.calctime(CurDate, CurTime); 1738 1739 int NumRecsHistory = History.HistoryNumberOfRecords(); 1740 double InitStrengthValue = Env.InitRawTValue; 1741 boolean LandFlagTF = Env.LandFlagTF; 1742 1743 int EyeSceneCount = 0; 1744 1745 if ((Main.HistoryFileName != null) || (ForceAutoFixUseTF_forTesting)) { 1746 HistoryRecFinalTno = 9.0; 1747 MaximumHistoryRecCI = 9.0; 1748 HistoryRecRawTno = 9.0; 1749 } else if (NumRecsHistory == 0) { 1750 HistoryRecFinalTno = InitStrengthValue; 1751 MaximumHistoryRecCI = InitStrengthValue; 1752 HistoryRecRawTno = InitStrengthValue; 1753 } else { 1754 EyeSceneCount = 0; 1755 HistoryRecFinalTno = History.IRCurrentRecord.Traw; 1756 HistoryRecRawTno = History.IRCurrentRecord.Traw; 1757 MaximumHistoryRecCI = History.IRCurrentRecord.Traw; 1758 int XInc = 0; 1759 while (XInc < NumRecsHistory) { 1760 int RecDate = History.HistoryFile[XInc].date; 1761 int RecTime = History.HistoryFile[XInc].time; 1762 int RecLand = History.HistoryFile[XInc].land; 1763 double RecTnoRaw = History.HistoryFile[XInc].Traw; 1764 double HistoryRecTime = Functions.calctime(RecDate, RecTime); 1765 boolean LandCheckTF = true; 1766 if (((LandFlagTF) && (RecLand == 1)) || (RecTnoRaw < 1.0)) { 1767 LandCheckTF = false; 1768 } 1769 if ((HistoryRecTime < CurrentTime) && (LandCheckTF)) { 1770 HistoryRecCI = History.HistoryFile[XInc].CI; 1771 HistoryRecFinalTno = History.HistoryFile[XInc].Tfinal; 1772 HistoryRecRawTno = History.HistoryFile[XInc].Traw; 1773 HistoryRecEyeScene = History.HistoryFile[XInc].eyescene; 1774 HistoryRecCloudScene = History.HistoryFile[XInc].cloudscene; 1775 if (HistoryRecCI > MaximumHistoryRecCI) { 1776 MaximumHistoryRecCI = HistoryRecCI; 1777 } 1778 if ((HistoryRecEyeScene < 3) 1779 || ((HistoryRecCloudScene == 1) && (HistoryRecEyeScene == 3))) { 1780 EyeSceneCount = EyeSceneCount + 1; 1781 /* checking for eye or embedded center scene types */ 1782 if (EyeSceneCount == 3) { 1783 FoundEyeSceneTF = true; 1784 } 1785 } else { 1786 EyeSceneCount = 0; 1787 } 1788 } 1789 XInc++; 1790 } 1791 if (((CurEyeScene < 3) || ((CurCloudScene == 1) && (CurEyeScene == 3))) 1792 && (EyeSceneCount == 2)) { 1793 FoundEyeSceneTF = true; 1794 } 1795 } 1796 1797 /* System.out.printf("%d %d : ",CurDate,CurTime); */ 1798 /* 1799 * System.out.printf("RingSpiralScore=%f HistoryRecFinalTno=%f\n", 1800 * RingSpiralScore,HistoryRecFinalTno); 1801 */ 1802 /* 1803 * System.out.printf("HistoryRecRawTno=%f MaximumHistoryRecCI=%f\n", 1804 * HistoryRecRawTno,MaximumHistoryRecCI); 1805 */ 1806 1807 /* 1808 * check score for developing systems (MaximumHistoryRecCI<5.0) starting 1809 * at 3.0 or check score for weakeining systems 1810 * (MaximumHistoryRecCI>=5.0) only above 3.5 1811 */ 1812 if (((HistoryRecRawTno >= 3.0) && (MaximumHistoryRecCI < 5.0)) 1813 || ((HistoryRecRawTno >= 3.5) && (MaximumHistoryRecCI >= 5.0))) { 1814 1815 if ((HistoryRecFinalTno <= 4.5) && (RingSpiralScore < 1.0) && (!FoundEyeSceneTF)) 1816 RingSpiralScore = -99.99; 1817 1818 /* System.out.printf(" FINALRingSpiralScore=%f ",RingSpiralScore); */ 1819 1820 if (RingSpiralScore > 0.0) { 1821 /* use Spiral/Ring methodology for center point */ 1822 FinalLatitude = RingSpiralLatitude; 1823 FinalLongitude = RingSpiralLongitude; 1824 FinalScoreValue = RingSpiralScore; 1825 PositioningMethodID = RingSpiralSelectionIDValue; 1826 } else { 1827 /* CDO... can't find anything to focus on */ 1828 FinalLatitude = ForecastLatitude; 1829 FinalLongitude = ForecastLongitude; 1830 FinalScoreValue = 0.0; 1831 PositioningMethodID = InputPositioningID; 1832 } 1833 } else { 1834 /* 1835 * current Tfinal is less than 3.5 or current scene is not an eye or 1836 * embedded center WILL USE FORECAST POSITION FOR AODT ANALYSIS 1837 */ 1838 FinalLatitude = ForecastLatitude; 1839 FinalLongitude = ForecastLongitude; 1840 FinalScoreValue = 0.0; 1841 PositioningMethodID = InputPositioningID; 1842 } 1843 1844 return new double[] { FinalLatitude, FinalLongitude, FinalScoreValue, 1845 (double) PositioningMethodID }; 1846 1847 } 1848 1849 /** 1850 * Calls routines to setup transformation, transform, data move. 1851 * 1852 * Input data provided with global variable arrays containing original and 1853 * transformed arrays. 1854 */ 1855 private static void RemapData() { 1856 /* 1857 * The following routines were originally developed by Dave Santek of 1858 * UW/SSEC and were added to the ADT under permission. If executed, an 1859 * array of latitude and longitude position arrays will be remapped to a 1860 * rectilinear projection for Tony Wimmers routines 1861 */ 1862 1863 int NumberOfCorners = 0; 1864 int LineSplineValue = 3; 1865 int ElementSplineValue = LineSplineValue; 1866 1867 tiff_vars.in_elems = Data.IRData_NumberColumns; 1868 tiff_vars.in_lines = Data.IRData_NumberRows; 1869 1870 /* 1871 * System.out.printf("elems=%d lines=%d\n", 1872 * tiff_vars.in_elems,tiff_vars.in_lines); 1873 */ 1874 1875 /* Uinit( ); */ 1876 /* 1877 * dis_vars.xrectl = (double)tiff_vars.in_lines; dis_vars.xrecte = 1878 * (double)tiff_vars.in_elems; 1879 * System.out.printf("xrectl=%f xrecte=%f\n", 1880 * dis_vars.xrectl,dis_vars.xrecte); 1881 */ 1882 1883 DetermineDest(); 1884 /* System.out.printf("after determineDest\n"); */ 1885 1886 NumberOfCorners = Init(LineSplineValue, ElementSplineValue); 1887 1888 System.out.printf("number of corners=%d\n", NumberOfCorners); 1889 Corner(NumberOfCorners, LineSplineValue, ElementSplineValue); 1890 /* System.out.printf("after corners\n"); */ 1891 1892 if ((ElementSplineValue > 1) || (LineSplineValue > 1)) { 1893 DoMap(NumberOfCorners, LineSplineValue, ElementSplineValue); 1894 } 1895 1896 } 1897 1898 /** 1899 * Interpolate between two arrays of different size. 1900 */ 1901 private static void DetermineDest() { 1902 int XInc, YInc; 1903 int XSizeMax = Data.IRData_NumberColumns - 1; 1904 int YSizeMax = Data.IRData_NumberRows - 1; 1905 1906 double NWCornerLatitude = Data.IRData_Latitude[0][0]; 1907 double NWCornerLongitude = Data.IRData_Longitude[0][0]; 1908 double NECornerLatitude = Data.IRData_Latitude[0][XSizeMax]; 1909 double NECornerLongitude = Data.IRData_Longitude[0][XSizeMax]; 1910 double SWCornerLatitude = Data.IRData_Latitude[YSizeMax][0]; 1911 double SWCornerLongitude = Data.IRData_Longitude[YSizeMax][0]; 1912 double SECornerLatitude = Data.IRData_Latitude[YSizeMax][XSizeMax]; 1913 double SECornerLongitude = Data.IRData_Longitude[YSizeMax][XSizeMax]; 1914 1915 /* crosses dateline check */ 1916 if ((NWCornerLongitude < NECornerLongitude) || (SWCornerLongitude < SECornerLongitude)) { 1917 NWCornerLongitude = NWCornerLongitude + 360.0; 1918 if (SWCornerLongitude < SECornerLongitude) { 1919 SWCornerLongitude = SWCornerLongitude + 360.0; 1920 } 1921 /* System.out.printf("DATELINE CROSS\n"); */ 1922 for (XInc = 0; XInc < XSizeMax; XInc++) { 1923 for (YInc = 0; YInc < YSizeMax; YInc++) { 1924 double DataValue = Data.IRData_Longitude[YInc][XInc]; 1925 if (DataValue < 0.0) { 1926 Data.IRData_Longitude[YInc][XInc] = (float) (DataValue + 360.0); 1927 } 1928 } 1929 } 1930 NWCornerLatitude = Data.IRData_Latitude[0][0]; 1931 NWCornerLongitude = Data.IRData_Longitude[0][0]; 1932 NECornerLatitude = Data.IRData_Latitude[0][XSizeMax]; 1933 NECornerLongitude = Data.IRData_Longitude[0][XSizeMax]; 1934 SWCornerLatitude = Data.IRData_Latitude[YSizeMax][0]; 1935 SWCornerLongitude = Data.IRData_Longitude[YSizeMax][0]; 1936 SECornerLatitude = Data.IRData_Latitude[YSizeMax][XSizeMax]; 1937 SECornerLongitude = Data.IRData_Longitude[YSizeMax][XSizeMax]; 1938 } 1939 1940 double MaximumLatitudeValue = Math.min(NWCornerLatitude, NECornerLatitude); 1941 double MinimumLatitudeValue = Math.max(SWCornerLatitude, SECornerLatitude); 1942 double MaximumLongitudeValue = Math.max(NWCornerLongitude, SWCornerLongitude); 1943 double MinimumLongitudeValue = Math.min(NECornerLongitude, SECornerLongitude); 1944 1945 double LatitudeIncrement = (MaximumLatitudeValue - MinimumLatitudeValue) 1946 / (double) Data.IRData_NumberColumns; 1947 double LongitudeIncrement = (MaximumLongitudeValue - MinimumLongitudeValue) 1948 / (double) Data.IRData_NumberRows; 1949 double MaximumIncrementValue = Math.max(LatitudeIncrement, LongitudeIncrement); 1950 System.out.printf("REMAPPING INFO\n"); 1951 System.out.printf("Source Array Bounds\n"); 1952 System.out.printf(" NW Corner : %7.2f/%7.2f\n", NWCornerLatitude, NWCornerLongitude); 1953 System.out.printf(" NE Corner : %7.2f/%7.2f\n", NECornerLatitude, NECornerLongitude); 1954 System.out.printf(" SW Corner : %7.2f/%7.2f\n", SWCornerLatitude, SWCornerLongitude); 1955 System.out.printf(" SE Corner : %7.2f/%7.2f\n", SECornerLatitude, SECornerLongitude); 1956 System.out.printf("Destination Array Bounds\n"); 1957 System.out.printf(" Max Lat/Lon: %7.2f/%7.2f\n", MaximumLatitudeValue, 1958 MaximumLongitudeValue); 1959 System.out.printf(" Min Lat/Lon: %7.2f/%7.2f\n", MinimumLatitudeValue, 1960 MinimumLongitudeValue); 1961 System.out.printf(" Inc Lat/Lon: %5.3f/ %5.3f\n", MaximumIncrementValue, 1962 MaximumIncrementValue); 1963 1964 tiff_vars.out_lines = (int) ((MaximumLatitudeValue - MinimumLatitudeValue) / MaximumIncrementValue); 1965 tiff_vars.out_elems = (int) ((MaximumLongitudeValue - MinimumLongitudeValue) / MaximumIncrementValue); 1966 for (YInc = 0; YInc < tiff_vars.out_lines; YInc++) { 1967 double LatitudeValue = MaximumLatitudeValue - (YInc * MaximumIncrementValue); 1968 for (XInc = 0; XInc < tiff_vars.out_elems; XInc++) { 1969 double LongitudeValue = MaximumLongitudeValue - (XInc * MaximumIncrementValue); 1970 IRData_Remap_Latitude[YInc][XInc] = LatitudeValue; 1971 IRData_Remap_Longitude[YInc][XInc] = LongitudeValue; 1972 } 1973 } 1974 1975 IRData_Remap_NumberColumns = tiff_vars.out_elems; 1976 IRData_Remap_NumberRows = tiff_vars.out_lines; 1977 } 1978 1979 /** 1980 * Compute number of corners for transformation and block sizes. 1981 * 1982 * @param LineSplineInput 1983 * Spline function for line values. 1984 * @param ElementSplineInput 1985 * Spline function for element values. 1986 * 1987 * @return Total number of corners to interpolate. 1988 */ 1989 private static int Init(int LineSplineInput, int ElementSplineInput) { 1990 int NumberOfCorners = 0; 1991 1992 remap_vars.nspl = (tiff_vars.out_elems + ElementSplineInput - 1) / ElementSplineInput; 1993 remap_vars.nspe = (tiff_vars.out_lines + LineSplineInput - 1) / LineSplineInput; 1994 1995 remap_vars.ncl = remap_vars.nspl + 1; 1996 remap_vars.nce = remap_vars.nspe + 1; 1997 1998 if (((tiff_vars.out_elems + ElementSplineInput - 1) % ElementSplineInput) == 0) { 1999 remap_vars.ncl = remap_vars.nspl; 2000 } 2001 if (((tiff_vars.out_lines + LineSplineInput - 1) % LineSplineInput) == 0) { 2002 remap_vars.nce = remap_vars.nspe; 2003 } 2004 2005 NumberOfCorners = remap_vars.ncl * remap_vars.nce; 2006 2007 remap_vars.in_bfw = Math.max(MINBFW, MINBLKSIZ * tiff_vars.in_elems); 2008 remap_vars.out_bfw = Math.max(MINBFW, Math.max(LineSplineInput, MINBLKSIZ) 2009 * tiff_vars.out_elems); 2010 2011 remap_vars.slb = remap_vars.in_bfw / tiff_vars.in_elems; 2012 remap_vars.dlb = ((remap_vars.out_bfw / tiff_vars.out_elems) / LineSplineInput) 2013 * LineSplineInput; 2014 2015 return NumberOfCorners; 2016 } 2017 2018 /** 2019 * Compute transformations at corners. 2020 * 2021 * Operates on the {@link #LineCoordinateArray} and 2022 * {@link #ElementCoordinateArray}. 2023 * 2024 * @param NumberOfCornersInput 2025 * Total number of corners to interpolate. 2026 * @param LineSplineInput 2027 * Spline function for line values. 2028 * @param ElementSplineInput 2029 * Spline function for element values. 2030 */ 2031 private static void Corner(int NumberOfCornersInput, int LineSplineInput, int ElementSplineInput) { 2032 int XInc; 2033 int YInc; 2034 int ArrayInc; 2035 2036 /* initialize array of corners */ 2037 for (ArrayInc = 0; ArrayInc < NumberOfCornersInput; ArrayInc++) { 2038 LineCoordinateArray[ArrayInc] = -99.9; 2039 ElementCoordinateArray[ArrayInc] = -99.9; 2040 } 2041 2042 /* loop through destination file and record source coords */ 2043 ArrayInc = -1; 2044 int NumberLines = tiff_vars.out_lines + LineSplineInput - 1; 2045 int NumberElements = tiff_vars.out_elems + ElementSplineInput - 1; 2046 2047 for (YInc = 0; YInc < NumberLines; YInc = YInc + LineSplineInput) { 2048 int LineValue = YInc; 2049 for (XInc = 0; XInc < NumberElements; XInc = XInc + ElementSplineInput) { 2050 /* System.out.printf("yinc=%d xinc=%d... \n",YInc,XInc); */ 2051 int ElementValue = XInc; 2052 int[] UMapReturn = UMap(LineValue, ElementValue); 2053 int RetErr = UMapReturn[0]; 2054 int SplineLineValue = UMapReturn[1]; 2055 int SplineElementValue = UMapReturn[2]; 2056 /* 2057 * System.out.printf("spline line=%d element=%d \n ",SplineLineValue 2058 * ,SplineElementValue); 2059 */ 2060 if ((ElementSplineInput == 1) && (LineSplineInput == 1)) { 2061 IRData_Remap_Temperature[LineValue][ElementValue] = Data.IRData_Temperature[SplineLineValue][SplineElementValue]; 2062 } else { 2063 ++ArrayInc; 2064 if (RetErr == 0) { 2065 LineCoordinateArray[ArrayInc] = SplineLineValue; 2066 ElementCoordinateArray[ArrayInc] = SplineElementValue; 2067 } 2068 } 2069 } 2070 } 2071 2072 } 2073 2074 /** 2075 * Provide coordinates between original point and transformed point. 2076 * 2077 * @param LineValueInput 2078 * Original line coordinate. 2079 * @param ElementValueInput 2080 * Ooriginal element coordinate. 2081 * 2082 * @return Array containing three values: possible error code, interpolated 2083 * line coordinate, and interpolated element coordinate. 2084 */ 2085 private static int[] UMap(int LineValueInput, int ElementValueInput) { 2086 2087 /* dest array line position value */ 2088 double DestinationLatitude; 2089 2090 /* dest array element position value */ 2091 double DestinationLongitude; 2092 2093 /* 2094 * Convert destination LineValueInput, ElementValueInput to source 2095 * coords LineValue_Return, ElementValue_Return 2096 */ 2097 LineValueInput = Math.min(LineValueInput, tiff_vars.out_lines - 1); 2098 ElementValueInput = Math.min(ElementValueInput, tiff_vars.out_elems - 1); 2099 DestinationLatitude = IRData_Remap_Latitude[LineValueInput][ElementValueInput]; 2100 DestinationLongitude = IRData_Remap_Longitude[LineValueInput][ElementValueInput]; 2101 2102 int[] FindPointReturn = FindPoint(DestinationLatitude, DestinationLongitude); 2103 int RetErr = FindPointReturn[0]; 2104 int LineValueReturn = FindPointReturn[1]; 2105 int ElementValueReturn = FindPointReturn[2]; 2106 2107 return new int[] { RetErr, LineValueReturn, ElementValueReturn }; 2108 } 2109 2110 /** 2111 * Find specific lat/lon location in array and return index values. 2112 * 2113 * @param latitude 2114 * Latitude value. 2115 * @param longitude 2116 * Longitude value. 2117 * 2118 * @return Array of three values. The first value is the status (-1 for 2119 * error, 0 for ok), the second value is array line value of lat/lon 2120 * input, and the third value is array element value of lat/lon 2121 * input. 2122 */ 2123 private static int[] FindPoint(double latitude, double longitude) { 2124 int RetErr; 2125 int XInc; 2126 int IndexValue; 2127 int YLineValueReturn; 2128 int XElementValueReturn; 2129 2130 double[] CornerLatitudeArray = new double[4]; 2131 double[] CornerLongitudeArray = new double[4]; 2132 double[] NSEWDistanceValuesArray = new double[4]; 2133 2134 /* found point logical */ 2135 boolean FoundPointTF = false; 2136 2137 /* out of bounds flag logical */ 2138 boolean OutOfBoundsTF = false; 2139 2140 /* found latitude value logical */ 2141 boolean FoundLatitudeTF = false; 2142 2143 /* found longitude value logical */ 2144 boolean FoundLongitudeTF = false; 2145 2146 int NumberElements = tiff_vars.in_elems; 2147 int NumberLines = tiff_vars.in_lines; 2148 int ElementValue = 0; 2149 int LineValue = 0; 2150 double PreviousDistance = 9999.9; 2151 2152 for (XInc = 0; XInc < 4; XInc++) { 2153 NSEWDistanceValuesArray[XInc] = 0.0; 2154 CornerLatitudeArray[XInc] = 0.0; 2155 CornerLongitudeArray[XInc] = 0.0; 2156 } 2157 while ((!FoundPointTF) && (!OutOfBoundsTF)) { 2158 CornerLatitudeArray[0] = Data.IRData_Latitude[LineValue][ElementValue]; 2159 CornerLongitudeArray[0] = Data.IRData_Longitude[LineValue][ElementValue]; 2160 CornerLatitudeArray[3] = Data.IRData_Latitude[LineValue + 1][ElementValue + 1]; 2161 CornerLongitudeArray[3] = Data.IRData_Longitude[LineValue + 1][ElementValue + 1]; 2162 /* 2163 * System.out.printf( 2164 * "x=%d y=%d : CornerLatitudeArray0=%f CornerLongitudeArray0=%f " 2165 * , " CornerLatitudeArray3=%f CornerLongitudeArray3=%f\n", 2166 * ElementValue,LineValue, 2167 * CornerLatitudeArray[0],CornerLongitudeArray[0], 2168 * CornerLatitudeArray[3],CornerLongitudeArray[3]); 2169 */ 2170 if ((longitude > CornerLongitudeArray[0]) || (longitude < CornerLongitudeArray[3])) { 2171 FoundLongitudeTF = false; 2172 if (longitude < CornerLongitudeArray[3]) { 2173 ElementValue++; 2174 } else { 2175 if (longitude > CornerLongitudeArray[0]) { 2176 ElementValue--; 2177 } 2178 } 2179 } else { 2180 FoundLongitudeTF = true; 2181 } 2182 if ((latitude > CornerLatitudeArray[0]) || (latitude < CornerLatitudeArray[3])) { 2183 FoundLatitudeTF = false; 2184 if (latitude < CornerLatitudeArray[3]) { 2185 LineValue++; 2186 } else { 2187 if (latitude > CornerLatitudeArray[0]) { 2188 LineValue--; 2189 } 2190 } 2191 } else { 2192 FoundLatitudeTF = true; 2193 } 2194 double LocalValue1[] = Functions.distance_angle(latitude, longitude, 2195 CornerLatitudeArray[0], CornerLongitudeArray[0], 1); 2196 double DistanceValue = LocalValue1[0]; 2197 2198 /* 2199 * * System.out.printf("distance : latitude=%f longitude=%f",* 2200 * " CornerLatitudeArray0=%f CornerLongitudeArray0=%f",* 2201 * " dist=%f angle=%f\n",latitude,longitude,* 2202 * CornerLatitudeArray[0],CornerLongitudeArray[0],* 2203 * DistanceValue,AngleValue); 2204 */ 2205 if (FoundLatitudeTF && FoundLongitudeTF) { 2206 FoundPointTF = true; 2207 } 2208 if (PreviousDistance <= DistanceValue) { 2209 FoundPointTF = true; 2210 } 2211 if ((ElementValue < 0) || (ElementValue > (NumberElements - 2))) { 2212 OutOfBoundsTF = true; 2213 } 2214 if ((LineValue < 0) || (LineValue > (NumberLines - 2))) { 2215 OutOfBoundsTF = true; 2216 } 2217 PreviousDistance = DistanceValue; 2218 } 2219 if (FoundPointTF) { 2220 CornerLatitudeArray[1] = Data.IRData_Latitude[LineValue][ElementValue + 1]; 2221 CornerLongitudeArray[1] = Data.IRData_Longitude[LineValue][ElementValue + 1]; 2222 CornerLatitudeArray[2] = Data.IRData_Latitude[LineValue + 1][ElementValue]; 2223 CornerLongitudeArray[2] = Data.IRData_Longitude[LineValue + 1][ElementValue]; 2224 2225 double LocalValue2[] = Functions.distance_angle(latitude, longitude, 2226 CornerLatitudeArray[0], CornerLongitudeArray[0], 1); 2227 NSEWDistanceValuesArray[0] = LocalValue2[0]; 2228 double LocalValue3[] = Functions.distance_angle(latitude, longitude, 2229 CornerLatitudeArray[1], CornerLongitudeArray[1], 1); 2230 NSEWDistanceValuesArray[1] = LocalValue3[0]; 2231 IndexValue = (NSEWDistanceValuesArray[0] < NSEWDistanceValuesArray[1]) ? 0 : 1; 2232 double LocalValue4[] = Functions.distance_angle(latitude, longitude, 2233 CornerLatitudeArray[2], CornerLongitudeArray[2], 1); 2234 NSEWDistanceValuesArray[2] = LocalValue4[0]; 2235 IndexValue = (NSEWDistanceValuesArray[IndexValue] < NSEWDistanceValuesArray[2]) ? IndexValue 2236 : 2; 2237 double LocalValue5[] = Functions.distance_angle(latitude, longitude, 2238 CornerLatitudeArray[3], CornerLongitudeArray[3], 1); 2239 NSEWDistanceValuesArray[3] = LocalValue5[0]; 2240 IndexValue = (NSEWDistanceValuesArray[IndexValue] < NSEWDistanceValuesArray[3]) ? IndexValue 2241 : 3; 2242 2243 XElementValueReturn = ((IndexValue == 0) || (IndexValue == 2)) ? ElementValue 2244 : ElementValue + 1; 2245 YLineValueReturn = ((IndexValue == 0) || (IndexValue == 1)) ? LineValue : LineValue + 1; 2246 RetErr = 0; 2247 } else { 2248 XElementValueReturn = -1; 2249 YLineValueReturn = -1; 2250 RetErr = -1; 2251 } 2252 2253 return new int[] { RetErr, YLineValueReturn, XElementValueReturn }; 2254 } 2255 2256 private static int DoMap(int NumberOfCornersInput, int LineSplineInput, int ElementSplineInput) { 2257 /* LineCoordsArrayInput = LineCoordinateArray */ 2258 /* ElementCoordsArrayInput = ElementCoordinateArray */ 2259 2260 /* line increment */ 2261 int[] LineIncArray1 = { -2, -2, 0, 2, 2, 2, 0, -2 }; 2262 2263 /* line increment */ 2264 int[] LineIncArray2 = { -1, -1, 0, 1, 1, 1, 0, -1 }; 2265 2266 /* ele increment */ 2267 int[] ElementIncArray1 = { 0, 2, 2, 2, 0, -2, -2, -2 }; 2268 2269 /* ele increment */ 2270 int[] ElementIncArray2 = { 0, 1, 1, 1, 0, -1, -1, -1 }; 2271 2272 int BufferSize = tiff_vars.in_lines * tiff_vars.in_elems; 2273 2274 double[] SourceArray = new double[BufferSize]; 2275 double[] InterpArray = new double[BufferSize]; 2276 double[] TempCharArray = new double[NumberOfCornersInput]; 2277 2278 int Ret_ID; 2279 int XInc, YInc, IncVal; 2280 int ArrayIndexValue; 2281 int ArrayPointValue; 2282 int SplinePixelValue; 2283 int SourceBlockValue; 2284 int SourceBlockInc; 2285 int SplineInc; 2286 int OffsetValue; 2287 int OffsetValue0; 2288 int OffsetValue1; 2289 int OffsetValue2; 2290 int OffsetValueL = 0; 2291 int OffsetValueE = 0; 2292 int CornerPointer; 2293 int AccumulatedLines; 2294 int FirstCornerPointer; 2295 int FirstCornerPointerHold; 2296 int LastCornerPointer; 2297 int ValX; 2298 int LineMinimumFinal; 2299 int LineMaximumFinal; 2300 int ElementMinimumFinal; 2301 int ElementMaximumFinal; 2302 int EdgeFixID; 2303 int SourceLine = 0; 2304 int SourceElement = 0; 2305 int SourceOffsetValue; 2306 int MaximumSourceLine; 2307 2308 double LineValueA; 2309 double LineValueB; 2310 double LineValueC; 2311 double ElementValueA; 2312 double ElementValueB; 2313 double ElementValueC; 2314 double LineValueACounter; 2315 double LineValueBCounter; 2316 double LineValueCCounter; 2317 double ElementValueACounter; 2318 double ElementValueBCounter; 2319 double ElementValueCCounter; 2320 double LineValueACounter0; 2321 double ElementValueACounter0; 2322 double LineULValue; 2323 double LineURValue; 2324 double LineLLValue; 2325 double LineLRValue; 2326 double LineMinimum; 2327 double LineMaximum; 2328 double ElementULValue; 2329 double ElementURValue; 2330 double ElementLLValue; 2331 double ElementLRValue; 2332 double ElementMinimum; 2333 double ElementMaximum; 2334 boolean ReadAllTF = true; 2335 boolean SKIP; 2336 boolean LABEL30; 2337 boolean LABEL38; 2338 2339 if ((SourceArray == null) || (InterpArray == null) || (TempCharArray == null)) { 2340 Ret_ID = -1; 2341 } else { 2342 SplinePixelValue = LineSplineInput * ElementSplineInput; 2343 2344 for (YInc = 0; YInc < tiff_vars.in_lines; YInc++) { 2345 for (XInc = 0; XInc < tiff_vars.in_elems; XInc++) { 2346 ArrayIndexValue = (YInc * tiff_vars.in_elems) + XInc; 2347 SourceArray[ArrayIndexValue] = Data.IRData_Temperature[YInc][XInc]; 2348 } 2349 } 2350 2351 /* 2352 * System.out.printf("remap nce=%d ncl=%d\n",remap_vars.nce,remap_vars 2353 * .ncl); 2354 */ 2355 for (YInc = 1; YInc < remap_vars.nce + 1; YInc++) { 2356 for (XInc = 1; XInc < remap_vars.ncl + 1; XInc++) { 2357 ArrayPointValue = IND(YInc, XInc); 2358 2359 if (LineCoordinateArray[ArrayPointValue - 1] == (double) -99.9) { 2360 double LineSum = 0.0; 2361 double ElementSum = 0.0; 2362 int PointCounter = 0; 2363 for (IncVal = 0; IncVal < 8; IncVal++) { 2364 SKIP = false; 2365 int YInc1 = YInc + LineIncArray1[IncVal]; 2366 int YInc2 = YInc + LineIncArray2[IncVal]; 2367 if ((YInc1 < 1) || (YInc2 < 1)) { 2368 SKIP = true; 2369 } 2370 if ((YInc1 > remap_vars.nce) || (YInc2 > remap_vars.nce)) { 2371 SKIP = true; 2372 } 2373 2374 int XInc1 = XInc + ElementIncArray1[IncVal]; 2375 int XInc2 = XInc + ElementIncArray2[IncVal]; 2376 2377 if ((XInc1 < 1) || (XInc2 < 1)) { 2378 SKIP = true; 2379 } 2380 if ((XInc1 > remap_vars.ncl) || (XInc2 > remap_vars.ncl)) { 2381 SKIP = true; 2382 } 2383 2384 int ValueIND1 = IND(YInc1, XInc1); 2385 int ValueIND2 = IND(YInc2, XInc2); 2386 2387 if (!SKIP) { 2388 if ((LineCoordinateArray[ValueIND1 - 1] == -99.9) 2389 || (LineCoordinateArray[ValueIND2 - 1] == -99.9)) { 2390 SKIP = true; 2391 } 2392 if (TempCharArray[ValueIND1 - 1] != 0) { 2393 SKIP = true; 2394 } 2395 if (TempCharArray[ValueIND2 - 1] != 0) { 2396 SKIP = true; 2397 } 2398 if (!SKIP) { 2399 PointCounter = PointCounter + 1; 2400 LineSum = LineSum 2401 + ((double) (2.0 * LineCoordinateArray[ValueIND2 - 1])) 2402 - LineCoordinateArray[ValueIND1 - 1]; 2403 ElementSum = ElementSum 2404 + ((double) 2.0 * ElementCoordinateArray[ValueIND2 - 1]) 2405 - ElementCoordinateArray[ValueIND1 - 1]; 2406 } /* SKIP:; */ 2407 } /* SKIP:; */ 2408 } 2409 2410 if (PointCounter > 0) { 2411 LineCoordinateArray[ArrayPointValue - 1] = LineSum / PointCounter; 2412 ElementCoordinateArray[ArrayPointValue - 1] = ElementSum / PointCounter; 2413 TempCharArray[ArrayPointValue - 1] = 1; 2414 } 2415 } 2416 } 2417 } 2418 2419 /* Loop through by destination blocks */ 2420 int BlockCounter = 0; 2421 2422 for (IncVal = 1; IncVal < tiff_vars.out_lines + 1; IncVal = IncVal + remap_vars.dlb) { 2423 /* Accumulated lines/block */ 2424 AccumulatedLines = BlockCounter * remap_vars.dlb; 2425 2426 /* Pointer to first corner of splines for this dest block */ 2427 FirstCornerPointer = (BlockCounter * remap_vars.ncl * remap_vars.dlb) 2428 / LineSplineInput; 2429 FirstCornerPointerHold = FirstCornerPointer; 2430 2431 /* Pointer to last corner for this dest block */ 2432 LastCornerPointer = (((BlockCounter + 1) * remap_vars.ncl * remap_vars.dlb) / LineSplineInput) - 1; 2433 LastCornerPointer = Math.min(LastCornerPointer, NumberOfCornersInput 2434 - remap_vars.ncl); 2435 2436 /* For each destination block loop through entire source */ 2437 for (SourceBlockInc = 1; SourceBlockInc < tiff_vars.in_lines + 1; SourceBlockInc = SourceBlockInc 2438 + remap_vars.slb) { 2439 MaximumSourceLine = Math.min(tiff_vars.in_lines, (SourceBlockInc 2440 + remap_vars.slb - 1)); 2441 ReadAllTF = false; 2442 2443 /* Loop through splines and move any data */ 2444 2445 FirstCornerPointer = FirstCornerPointerHold; 2446 SourceBlockValue = 0; 2447 while (FirstCornerPointer < LastCornerPointer) { 2448 for (SplineInc = 1; SplineInc < remap_vars.nspl + 1; SplineInc++) { 2449 LABEL30 = false; 2450 OffsetValue0 = (SourceBlockValue / remap_vars.nspl) * LineSplineInput 2451 * tiff_vars.out_elems; 2452 OffsetValue1 = (SourceBlockValue % remap_vars.nspl) 2453 * ElementSplineInput; 2454 OffsetValue2 = OffsetValue1 + OffsetValue0 + 1; 2455 CornerPointer = FirstCornerPointer + SplineInc - 1; 2456 2457 /* 2458 * Get 4 corners in line space and check for out of 2459 * bounds 2460 */ 2461 2462 ValX = remap_vars.ncl; 2463 LineULValue = LineCoordinateArray[CornerPointer]; 2464 LineURValue = LineCoordinateArray[CornerPointer + 1]; 2465 LineLLValue = LineCoordinateArray[CornerPointer + ValX]; 2466 if ((CornerPointer + 1 + ValX) < NumberOfCornersInput) { 2467 LineLRValue = LineCoordinateArray[CornerPointer + 1 + ValX]; 2468 } else { 2469 LineLRValue = LineCoordinateArray[CornerPointer + ValX]; 2470 } 2471 LineMinimum = Math.min(LineULValue, LineURValue); 2472 LineMinimum = Math.min(LineMinimum, LineLLValue); 2473 LineMinimum = Math.min(LineMinimum, LineLRValue); 2474 2475 /* Test for the presence of a limb in the spline box */ 2476 2477 if (LineMinimum == -99.9) { 2478 LABEL30 = true; 2479 } 2480 2481 LineMinimumFinal = (int) (LineMinimum + 0.5); 2482 if (LineMinimumFinal > MaximumSourceLine) { 2483 LABEL30 = true; 2484 } 2485 LineMaximum = Math.max(LineULValue, LineURValue); 2486 LineMaximum = Math.max(LineMaximum, LineLLValue); 2487 LineMaximum = Math.max(LineMaximum, LineLRValue); 2488 LineMaximumFinal = (int) (LineMaximum + 0.5); 2489 if (LineMaximumFinal < SourceBlockInc) { 2490 LABEL30 = true; 2491 } 2492 2493 /* 2494 * Get 4 corners in elem space & check for out of 2495 * bound 2496 */ 2497 2498 ValX = remap_vars.ncl; 2499 ElementULValue = ElementCoordinateArray[CornerPointer]; 2500 ElementURValue = ElementCoordinateArray[CornerPointer + 1]; 2501 ElementLLValue = ElementCoordinateArray[CornerPointer + ValX]; 2502 if ((CornerPointer + 1 + ValX) < NumberOfCornersInput) { 2503 ElementLRValue = ElementCoordinateArray[CornerPointer + 1 + ValX]; 2504 } else { 2505 ElementLRValue = ElementCoordinateArray[CornerPointer + ValX]; 2506 } 2507 2508 ElementMaximum = Math.max(ElementULValue, ElementURValue); 2509 ElementMaximum = Math.max(ElementMaximum, ElementLLValue); 2510 ElementMaximum = Math.max(ElementMaximum, ElementLRValue); 2511 ElementMaximumFinal = (int) (ElementMaximum + 0.5); 2512 if (ElementMaximumFinal < 1) { 2513 LABEL30 = true; 2514 } 2515 2516 ElementMinimum = Math.min(ElementULValue, ElementURValue); 2517 ElementMinimum = Math.min(ElementMinimum, ElementLLValue); 2518 ElementMinimum = Math.min(ElementMinimum, ElementLRValue); 2519 ElementMinimumFinal = (int) (ElementMinimum + 0.5); 2520 2521 if (ElementMinimumFinal > tiff_vars.in_elems) { 2522 LABEL30 = true; 2523 } 2524 EdgeFixID = 0; 2525 2526 /* 2527 * If the max & min element fall off the 2528 * image...pitch it 2529 */ 2530 2531 if ((ElementMaximumFinal > tiff_vars.in_elems) 2532 && (ElementMinimumFinal < 1)) { 2533 LABEL30 = true; 2534 } 2535 2536 /* Fix if left & right edge should be continuous */ 2537 2538 if (!LABEL30) { 2539 if ((ElementMaximumFinal - ElementMinimumFinal) > (int) (.75 * tiff_vars.in_elems)) { 2540 if (ElementULValue < (tiff_vars.in_elems / 2)) { 2541 ElementULValue = ElementULValue + tiff_vars.in_elems; 2542 } 2543 if (ElementURValue < (tiff_vars.in_elems / 2)) { 2544 ElementURValue = ElementURValue + tiff_vars.in_elems; 2545 } 2546 if (ElementLLValue < (tiff_vars.in_elems / 2)) { 2547 ElementLLValue = ElementLLValue + tiff_vars.in_elems; 2548 } 2549 if (ElementLRValue < (tiff_vars.in_elems / 2)) { 2550 ElementLRValue = ElementLRValue + tiff_vars.in_elems; 2551 } 2552 EdgeFixID = 1; 2553 } 2554 2555 LineValueA = (LineURValue - LineULValue) / ElementSplineInput; 2556 LineValueB = (LineLLValue - LineULValue) / LineSplineInput; 2557 LineValueC = (LineLRValue + LineULValue - LineURValue - LineLLValue) 2558 / SplinePixelValue; 2559 ElementValueA = (ElementURValue - ElementULValue) 2560 / ElementSplineInput; 2561 ElementValueB = (ElementLLValue - ElementULValue) / LineSplineInput; 2562 ElementValueC = (ElementLRValue + ElementULValue - ElementLLValue - ElementURValue) 2563 / SplinePixelValue; 2564 2565 int ElementMiscValue = 0; 2566 LineValueBCounter = LineULValue + 0.5; 2567 LineValueCCounter = 0.0; 2568 ElementValueBCounter = ElementULValue + 0.5; 2569 ElementValueCCounter = 0.0; 2570 2571 if (ReadAllTF == false) { 2572 ReadAllTF = true; 2573 } 2574 2575 if ((SplineInc == remap_vars.nspl) 2576 || (ElementMinimumFinal < 1 || (ElementMaximumFinal > tiff_vars.in_elems)) 2577 || ((LineMinimumFinal < SourceBlockInc) || (LineMaximumFinal > MaximumSourceLine)) 2578 || ((FirstCornerPointer + (2 * remap_vars.ncl) - 1) > LastCornerPointer)) { 2579 for (YInc = 1; YInc < LineSplineInput + 1; YInc++) { 2580 LineValueACounter = LineValueCCounter + LineValueA; 2581 ElementValueACounter = ElementValueCCounter + ElementValueA; 2582 LineValueACounter0 = 0.0; 2583 ElementValueACounter0 = 0.0; 2584 for (XInc = 1; XInc < ElementSplineInput + 1; XInc++) { 2585 LABEL38 = false; 2586 SourceLine = (int) (LineValueBCounter + LineValueACounter0); 2587 if (SourceLine < SourceBlockInc) { 2588 LABEL38 = true; 2589 } 2590 if (SourceLine > MaximumSourceLine) { 2591 LABEL38 = true; 2592 } 2593 if (!LABEL38) 2594 SourceElement = (int) (ElementValueBCounter + ElementValueACounter0); 2595 if (SourceElement < 1) { 2596 LABEL38 = true; 2597 } 2598 if ((SourceElement > tiff_vars.in_elems) 2599 && (EdgeFixID == 0)) { 2600 LABEL38 = true; 2601 } 2602 if (!LABEL38) { 2603 if (SourceElement > tiff_vars.in_elems) { 2604 SourceElement = SourceElement 2605 - tiff_vars.in_elems; 2606 } 2607 SourceOffsetValue = ((SourceLine - SourceBlockInc) * tiff_vars.in_elems) 2608 + SourceElement; 2609 OffsetValueE = OffsetValue1 + XInc; 2610 if (OffsetValueE > tiff_vars.out_elems) { 2611 LABEL38 = true; 2612 } 2613 if (!LABEL38) 2614 OffsetValueL = OffsetValue0 + ElementMiscValue; 2615 if ((OffsetValueL / (tiff_vars.out_elems 2616 + AccumulatedLines - 1)) > tiff_vars.out_lines) { 2617 LABEL38 = true; 2618 } 2619 if (!LABEL38) { 2620 OffsetValue = OffsetValueL + OffsetValueE; 2621 InterpArray[OffsetValue - 1] = SourceArray[SourceOffsetValue - 1]; 2622 } /* LABEL38: */ 2623 } /* LABEL38: */ 2624 LineValueACounter0 = LineValueACounter0 2625 + LineValueACounter; 2626 ElementValueACounter0 = ElementValueACounter0 2627 + ElementValueACounter; 2628 } 2629 LineValueBCounter = LineValueBCounter + LineValueB; 2630 LineValueCCounter = LineValueCCounter + LineValueC; 2631 ElementValueBCounter = ElementValueBCounter + ElementValueB; 2632 ElementValueCCounter = ElementValueCCounter + ElementValueC; 2633 ElementMiscValue = ElementMiscValue + tiff_vars.out_elems; 2634 } 2635 } else { 2636 if (EdgeFixID == 0) { 2637 for (YInc = 1; YInc < LineSplineInput + 1; YInc++) { 2638 LineValueACounter = LineValueCCounter + LineValueA; 2639 ElementValueACounter = ElementValueCCounter 2640 + ElementValueA; 2641 LineValueACounter0 = 0.0; 2642 ElementValueACounter0 = 0.0; 2643 OffsetValue = OffsetValue2 + ElementMiscValue; 2644 for (XInc = 1; XInc < ElementSplineInput + 1; XInc++) { 2645 SourceLine = (int) (LineValueBCounter + LineValueACounter0); 2646 SourceElement = (int) (ElementValueBCounter + ElementValueACounter0); 2647 SourceOffsetValue = ((SourceLine - SourceBlockInc) * tiff_vars.in_elems) 2648 + SourceElement; 2649 InterpArray[OffsetValue - 1] = SourceArray[SourceOffsetValue - 1]; 2650 OffsetValue = OffsetValue + 1; 2651 LineValueACounter0 = LineValueACounter0 2652 + LineValueACounter; 2653 ElementValueACounter0 = ElementValueACounter0 2654 + ElementValueACounter; 2655 } 2656 LineValueBCounter = LineValueBCounter + LineValueB; 2657 LineValueCCounter = LineValueCCounter + LineValueC; 2658 ElementValueBCounter = ElementValueBCounter 2659 + ElementValueB; 2660 ElementValueCCounter = ElementValueCCounter 2661 + ElementValueC; 2662 ElementMiscValue = ElementMiscValue 2663 + tiff_vars.out_elems; 2664 } 2665 } else if (EdgeFixID == 1) { 2666 for (YInc = 1; YInc < LineSplineInput + 1; YInc++) { 2667 LineValueACounter = LineValueCCounter + LineValueA; 2668 ElementValueACounter = ElementValueCCounter 2669 + ElementValueA; 2670 LineValueACounter0 = 0.0; 2671 ElementValueACounter0 = 0.0; 2672 OffsetValue = OffsetValue2 + ElementMiscValue; 2673 for (XInc = 1; XInc < ElementSplineInput + 1; XInc++) { 2674 SourceLine = (int) (LineValueBCounter + LineValueACounter0); 2675 SourceElement = (int) (ElementValueBCounter + ElementValueACounter0); 2676 if (SourceElement > tiff_vars.in_elems) { 2677 SourceElement = SourceElement 2678 - tiff_vars.in_elems; 2679 } 2680 SourceOffsetValue = ((SourceLine - SourceBlockInc) * tiff_vars.in_elems) 2681 + SourceElement; 2682 InterpArray[OffsetValue - 1] = SourceArray[SourceOffsetValue - 1]; 2683 OffsetValue = OffsetValue + 1; 2684 LineValueACounter0 = LineValueACounter0 2685 + LineValueACounter; 2686 ElementValueACounter0 = ElementValueACounter0 2687 + ElementValueACounter; 2688 } 2689 LineValueBCounter = LineValueBCounter + LineValueB; 2690 LineValueCCounter = LineValueCCounter + LineValueC; 2691 ElementValueBCounter = ElementValueBCounter 2692 + ElementValueB; 2693 ElementValueCCounter = ElementValueCCounter 2694 + ElementValueC; 2695 ElementMiscValue = ElementMiscValue 2696 + tiff_vars.out_elems; 2697 } 2698 } 2699 } 2700 } /* LABEL30: */ 2701 SourceBlockValue = SourceBlockValue + 1; 2702 } 2703 FirstCornerPointer = FirstCornerPointer + remap_vars.ncl; 2704 } 2705 } 2706 BlockCounter = BlockCounter + 1; 2707 } 2708 2709 for (YInc = 0; YInc < tiff_vars.out_lines; YInc++) { 2710 for (XInc = 0; XInc < tiff_vars.out_elems; XInc++) { 2711 ArrayIndexValue = (YInc * tiff_vars.out_elems) + XInc; 2712 if ((ArrayIndexValue > 0) && (InterpArray[ArrayIndexValue] == 0.0)) { 2713 InterpArray[ArrayIndexValue] = InterpArray[ArrayIndexValue - 1]; 2714 } 2715 IRData_Remap_Temperature[YInc][XInc] = InterpArray[ArrayIndexValue]; 2716 } 2717 } 2718 Ret_ID = 0; 2719 } 2720 2721 return Ret_ID; 2722 2723 } 2724 2725 private static int IND(int y, int x) { 2726 /* #define IND(y,x) ((y-1)*remap_vars.ncl+x) */ 2727 return ((y - 1) * remap_vars.ncl + x); 2728 } 2729 2730}