1 module dsurf.cartesian; 2 3 import std.range; 4 import std.stdio; 5 import std..string; 6 import std.algorithm; 7 import std.file; 8 import std.conv; 9 import std.math; 10 11 import mir.ndslice; 12 13 /*********************************** 14 * Class CartesianSurface represents regular rectangular-cell surface 15 * Authors: Dmitriy Linev 16 * License: MIT 17 * Date: 2019-2020 18 */ 19 class CartesianSurface { 20 21 /** 22 * Default constructor, doesn't allocate memory for height map 23 */ 24 this() pure { 25 m_xOrigin = 0; 26 m_yOrigin = 0; 27 m_dx = 0; 28 m_dy = 0; 29 m_nx = 0; 30 m_ny = 0; 31 } 32 33 /** 34 * Sets surface header and allocates memory for height map 35 */ 36 37 /** 38 * Sets surface header and allocates surface memory 39 * Params: 40 * nx = number of nodes along X coordinate 41 * ny = number of nodes along Y coordinate 42 * xOrigin = surface origin X coordinate 43 * yOrigin = surface origin Y coordinate 44 * dx = surface increment (cell size) along X coordinate 45 * dy = surface increment (cell size) along Y coordinate 46 */ 47 void setHeader(int nx, int ny, double xOrigin, double yOrigin, double dx, double dy) pure { 48 m_nx = nx; 49 m_ny = ny; 50 m_xOrigin = xOrigin; 51 m_yOrigin = yOrigin; 52 m_dx = dx; 53 m_dy = dy; 54 m_z = slice!double(nx, ny); 55 } 56 57 unittest { 58 CartesianSurface s1 = new CartesianSurface; 59 s1.setHeader(2, 2, 0, 0, 500, 500); 60 s1.m_z[] = 50; 61 for (int i = 0; i < s1.nx; i++) { 62 for (int j = 0; j < s1.ny; j++) { 63 assert(s1.z[i][j] == 50); 64 } 65 } 66 } 67 68 /// Copy constructor, returns the exact copy of the given surface 69 this(CartesianSurface surface) pure { 70 this.setHeader(surface.nx, surface.ny, surface.xOrigin, surface.yOrigin, surface.dx, surface.dy); 71 this.m_z[] = surface.m_z[]; 72 } 73 74 unittest { 75 CartesianSurface s1 = new CartesianSurface; 76 s1.setHeader(2, 2, 0, 0, 500, 500); 77 for (int i = 0; i < s1.nx; i++) { 78 for (int j = 0; j < s1.ny; j++) { 79 s1.z[i][j] = 50; 80 } 81 } 82 83 CartesianSurface s2 = new CartesianSurface(s1); 84 for (int i = 0; i < s1.nx; i++) { 85 for (int j = 0; j < s1.ny; j++) { 86 assert(s2.z[i][j] == 50); 87 } 88 } 89 90 for (int i = 0; i < s1.nx; i++) { 91 for (int j = 0; j < s1.ny; j++) { 92 s2.z[i][j] = 10; 93 } 94 } 95 96 for (int i = 0; i < s1.nx; i++) { 97 for (int j = 0; j < s1.ny; j++) { 98 assert(s2.z[i][j] == 10); 99 } 100 } 101 102 for (int i = 0; i < s1.nx; i++) { 103 for (int j = 0; j < s1.ny; j++) { 104 assert(s1.z[i][j] == 50); 105 } 106 } 107 } 108 109 /// Returns X coordinate of surface origin 110 @property pure double xOrigin() const @safe @nogc { return m_xOrigin; } 111 112 /// Returns Y coordinate of surface origin 113 @property pure double yOrigin() const @safe @nogc { return m_yOrigin; } 114 115 /// Returns surface increment (cell size) along X axis 116 @property pure double dx() const @safe @nogc { return m_dx; } 117 118 /// Returns surface increment (cell size) along Y axis 119 @property pure double dy() const @safe @nogc { return m_dy; } 120 121 /// Returns number of nodes along X axis 122 @property pure int nx() const @safe @nogc { return m_nx; } 123 124 /// Returns number of nodes along Y axis 125 @property pure int ny() const @safe @nogc { return m_ny; } 126 127 128 /** 129 Method to access height map. 130 Returns: `Slice!(double*, 2)` containing surface's height map with dimensions nx * ny 131 Example: 132 --- 133 for (int i = 0; i < surface.nx; i++) { 134 for (int j = 0; j < surface.ny; j++) { 135 surface.z[i][j] = 0; 136 } 137 } 138 --- 139 */ 140 @property Slice!(double*, 2) z() { return m_z; } 141 142 /// Returns `true` if point with given coordinates is inside surface boundaries, otherwise returns `false` 143 bool isInsideSurface(double x, double y) const @nogc { 144 if (x < m_xOrigin || y < m_yOrigin || x > m_xOrigin + m_dx * (m_nx - 1) || y > m_yOrigin + m_dy * (m_ny - 1)) 145 return false; 146 return true; 147 } 148 149 /// Returns cell index from a given `X` coordinate. Maximum number is `nx - 1` 150 /// Returns -1 if a given coordinate is outside the surface 151 int cellXIndex(double x) const { 152 if (x < m_xOrigin || x > m_xOrigin + m_dx * (m_nx - 1)) 153 return -1; //TODO throw? 154 return ((x - m_xOrigin) / m_dx).to!int; 155 } 156 157 /// Returns cell index from a given `Y` coordinate. Maximum number is `ny - 1` 158 /// Returns `-1` if a given coordinate is outside the surface 159 int cellYIndex(double y) const { 160 if (y < m_yOrigin || y > m_yOrigin + m_dy * (m_ny - 1)) 161 return -1; //TODO throw? 162 return ((y - m_yOrigin) / m_dy).to!int; 163 } 164 165 unittest { 166 CartesianSurface surface = new CartesianSurface; 167 surface.setHeader(10, 10, 0, 1000, 500, 500); 168 assert(surface.cellXIndex(-1000) == -1); 169 assert(surface.cellYIndex(-1000) == -1); 170 assert(surface.cellXIndex(0) == 0); 171 assert(surface.cellYIndex(1000) == 0); 172 assert(surface.cellXIndex(550) == 1); 173 assert(surface.cellYIndex(1550) == 1); 174 assert(surface.cellXIndex(5000) == -1); 175 assert(surface.cellYIndex(6000) == -1); 176 } 177 178 /// Returns z value of a point with given coordinate using bilinear interpolation 179 double getZ(double x, double y) const { 180 if (!isInsideSurface(x, y)) 181 return double.nan; 182 immutable int i = cellXIndex(x); 183 immutable int j = cellYIndex(y); 184 185 if (i == m_nx - 1 && j == m_ny - 1) //top right corner 186 return m_z[$ - 1][$ - 1]; 187 else if (i == m_nx - 1) { //right edge 188 return m_z[i][j] + (m_z[i][j + 1] - m_z[i][j]) / m_dy * (y - (m_yOrigin + j * m_dy)); 189 } 190 else if (j == m_ny - 1) { //top edge 191 return m_z[i][j] + (m_z[i + 1][j] - m_z[i][j]) / m_dx * (x - (m_xOrigin + i * m_dx)); 192 } 193 194 //z: top left, top right, bottom left, bottom right 195 immutable double ztl = m_z[i][j + 1]; 196 immutable double zbl = m_z[i][j]; 197 immutable double ztr = m_z[i + 1][j + 1]; 198 immutable double zbr = m_z[i + 1][j]; 199 200 //FIXME handle blanks 201 immutable double z1 = (m_xOrigin + (i + 1) * m_dx - x) / dx * zbl + 202 (x - (m_xOrigin + i * m_dx)) / dx * zbr; 203 immutable double z2 = (m_xOrigin + (i + 1) * m_dx - x) / dx * ztl + 204 (x - (m_xOrigin + i * m_dx)) / dx * ztr; 205 return (m_yOrigin + (j + 1) * m_dy - y) / m_dy * z1 + 206 (y - (m_yOrigin + j * m_dy)) / m_dy * z2; 207 } 208 209 /// Operators +=, -=, *=, /= overloading for two surfaces 210 CartesianSurface opOpAssign(string op)(CartesianSurface rhs) { 211 for (int i = 0; i < m_nx; i++) { 212 for (int j = 0; j < m_ny; j++) { 213 immutable double x = m_xOrigin + i * dx; 214 immutable double y = m_yOrigin + j * dy; 215 immutable double rhsz = rhs.getZ(x, y); 216 if (isNaN(rhsz)) 217 continue; 218 static if (op == "+") { 219 m_z[i][j] += rhsz; 220 } 221 else static if (op == "-") { 222 m_z[i][j] -= rhsz; 223 } 224 else static if (op == "*") { 225 m_z[i][j] *= rhsz; 226 } 227 else static if (op == "/") { 228 if (abs(rhsz) < 1e-9) 229 m_z[i][j] = math.nan; 230 else 231 m_z[i][j] -= rhsz; 232 } 233 else static assert(0, "Operator "~op~" not implemented"); 234 } 235 } 236 return this; 237 } 238 239 /// Operators +, -, *, / overloading for two surfaces 240 CartesianSurface opBinary(string op)(CartesianSurface rhs) { 241 CartesianSurface result = new CartesianSurface(this); 242 static if (op == "+") { 243 result += rhs; 244 } 245 else static if (op == "-") { 246 result -= rhs; 247 } 248 else static if (op == "*") { 249 result *= result; 250 } 251 else static if (op == "/") { 252 result /= rhs; 253 } 254 else static assert(0, "Operator "~op~" not implemented"); 255 return result; 256 } 257 258 /// Operators +=, -=, *=, /= overloading for a surface and a fixed value 259 CartesianSurface opOpAssign(string op)(double rhs) { 260 for (int i = 0; i < m_nx; i++) { 261 for (int j = 0; j < m_ny; j++) { 262 static if (op == "+") { 263 m_z[i][j] += rhs; 264 } 265 else static if (op == "-") { 266 m_z[i][j] -= rhs; 267 } 268 else static if (op == "*") { 269 m_z[i][j] *= rhs; 270 } 271 else static if (op == "/") { 272 m_z[i][j] -= rhs; 273 } 274 else static assert(0, "Operator "~op~" not implemented"); 275 } 276 } 277 return this; 278 279 } 280 281 /// Operators +, -, *, / overloading for a surface and a fixed value 282 CartesianSurface opBinary(string op)(double rhs) { 283 CartesianSurface result = new CartesianSurface(this); 284 static if (op == "+") { 285 result += rhs; 286 } 287 else static if (op == "-") { 288 result -= rhs; 289 } 290 else static if (op == "*") { 291 result *= result; 292 } 293 else static if (op == "/") { 294 result /= rhs; 295 } 296 else static assert(0, "Operator "~op~" not implemented"); 297 return result; 298 } 299 300 private: 301 Slice!(double*, 2) m_z; 302 double m_xOrigin; 303 double m_yOrigin; 304 double m_dx; 305 double m_dy; 306 int m_nx; 307 int m_ny; 308 } 309 310 /** 311 * Loads `surface` from file trying to detect format automatically 312 * Params: 313 * surface = `CartesianSurface` to load data to 314 * fileName = Path to file for loading 315 * Currently supported formats are CPS3 ASCII and ZMAP+ 316 */ 317 void loadFromFile(CartesianSurface surface, string fileName) { 318 immutable auto format = surfaceFormat(fileName); 319 if (format == "cps") 320 loadFromCps3Ascii(surface, fileName); 321 else if (format == "zmap") 322 loadFromZmap(surface, fileName); 323 else 324 throw new FileException("Unknown surface format"); 325 } 326 327 /** 328 * Loads `surface` from file of CPS3 ASCII format 329 * Params: 330 * surface = `CartesianSurface` to load data to 331 * fileName = Path to file for loading 332 */ 333 void loadFromCps3Ascii(CartesianSurface surface, string fileName) { 334 File file = File(fileName, "r"); 335 bool readingHeader = true; 336 int i = -1; 337 int j = -1; 338 double xOrigin = -1, yOrigin = -1; 339 double dx = -1, dy = -1; 340 int nx = -1, ny = -1; 341 double blank = 1e31; 342 while(!file.eof()) { 343 string line = file.readln().chomp().replace("\"", ""); 344 if (!line) 345 continue; 346 if (readingHeader) { 347 auto words = line.split(); 348 if (words[0] == "FSASCI") { 349 blank = words[$ - 1].to!double; 350 } 351 if (words[0] == "FSLIMI") { 352 xOrigin = words[1].to!double; 353 yOrigin = words[3].to!double; 354 } 355 else if (words[0] == "FSNROW") { 356 nx = words[2].to!int; 357 ny = words[1].to!int; 358 } 359 else if (words[0] == "FSXINC") { 360 dx = words[1].to!double; 361 dy = words[2].to!double; 362 readingHeader = false; 363 surface.setHeader(nx, ny, xOrigin, yOrigin, dx, dy); 364 i = 0; 365 j = ny - 1; 366 } 367 } 368 else { 369 if (i < 0 || j < 0) 370 throw new Exception("Invalid index"); //TODO add some information 371 auto words = line.split(); 372 if (words[0].startsWith("->")) //petrel specific 373 continue; 374 foreach (word; words) { 375 double value = 0; 376 try 377 value = to!double(word); 378 catch (ConvException e) 379 value = double.nan; 380 if (value == blank) 381 surface.z[i][j] = double.nan; 382 else 383 surface.z[i][j] = value; 384 385 j--; 386 if (j < 0) { 387 i++; 388 j = surface.ny() - 1; 389 } 390 } 391 } 392 } 393 file.close(); 394 } 395 396 /** 397 * Loads `surface` from file of ZMap+ ASCII format 398 * Params: 399 * surface = `CartesianSurface` to load data to 400 * fileName = Path to file for loading 401 */ 402 void loadFromZmap(CartesianSurface surface, string fileName) { 403 File file = File(fileName, "r"); 404 bool readingHeader = true; 405 int i = -1; 406 int j = -1; 407 double xOrigin = -1, yOrigin = -1; 408 double xMax = -1, yMax = -1; 409 double dx = -1, dy = -1; 410 int nx = -1, ny = -1; 411 double blank = double.nan; 412 while(!file.eof()) { 413 string line = file.readln().chomp().replace("\"", ""); 414 if (readingHeader) { 415 if (line.startsWith("!")) 416 continue; 417 string [] words = line.replace(",", "").split(); 418 if (isNaN(blank)) { 419 if (line.startsWith("@")) 420 continue; 421 blank = words[1].to!double; 422 } 423 else if (nx < 0) { 424 nx = words[1].to!int; 425 ny = words[0].to!int; 426 xOrigin = words[2].to!double; 427 yOrigin = words[4].to!double; 428 xMax = words[3].to!double; 429 yMax = words[5].to!double; 430 } 431 else if (dx < 0) { 432 dx = words[0].to!double; 433 dy = words[1].to!double; 434 if (approxEqual(0.0, dx, 0.0) || approxEqual(0.0, dy, 0.0)) { 435 dx = (xMax - xOrigin) / (nx - 1); 436 dy = (yMax - yOrigin) / (ny - 1); 437 } 438 } 439 else if (line.startsWith("@")) { 440 surface.setHeader(nx, ny, xOrigin, yOrigin, dx, dy); 441 readingHeader = false; 442 i = 0; 443 j = ny - 1; 444 } 445 } 446 else { 447 if (i < 0 || j < 0) 448 throw new Exception("Invalid index"); //TODO add some information 449 auto words = line.split(); 450 if (words.empty) 451 continue; 452 if (words[0].startsWith("+")) //RMS specific 453 continue; 454 foreach (word; words) { 455 double value = 0; 456 try 457 value = to!double(word); 458 catch (ConvException e) 459 value = double.nan; 460 if (value == blank) 461 surface.z[i][j] = double.nan; 462 else 463 surface.z[i][j] = value; 464 465 j--; 466 if (j < 0) { 467 i++; 468 j = surface.ny() - 1; 469 } 470 } 471 } 472 } 473 } 474 475 476 /** 477 * Saves `surface` to file using specified `format` 478 * Params: 479 * surface = surface to save 480 * fileName = path to output file 481 * format = format to save surface to. Currently supported formats for export are: CPS3 ASCII 482 */ 483 void saveToFile(CartesianSurface surface, string fileName, string format) { 484 saveToCps3Ascii(surface, fileName); 485 if (format.startsWith("cps".toLower)) 486 saveToCps3Ascii(surface, fileName); 487 else 488 throw new FileException(format ~ " cartesian surface format is not supported for export"); 489 } 490 491 /** 492 * Saves `surface` to file using CPS3 ASCII format 493 * Params: 494 * surface = surface to save 495 * fileName = path to output file 496 */ 497 void saveToCps3Ascii(CartesianSurface surface, string fileName) { 498 File file = File(fileName, "w"); 499 immutable double blank = 1e30; 500 file.writeln("FSASCI 0 1 COMPUTED 0 ", blank); 501 file.writeln("FSATTR 0 0"); 502 file.writeln("FSLIMI ", 503 surface.m_xOrigin, " ", 504 surface.m_xOrigin + (surface.m_nx) * surface.m_dx, " ", 505 surface.m_yOrigin, " ", 506 surface.m_yOrigin + (surface.m_ny) * surface.m_dy, " ", 507 surface.m_z[surface.m_z.minIndex], " ", 508 surface.m_z[surface.m_z.maxIndex]); 509 file.writeln("FSNROW ", surface.m_ny, " ", surface.m_nx); 510 file.writeln("FSXINC ", surface.m_dx, " ", surface.m_dy); 511 int n = 0; 512 for (int i = 0; i < surface.nx; i++) { 513 for (int j = surface.ny - 1; j >= 0; j--) { 514 if (isNaN(surface.z[i][j])) 515 file.write(blank); 516 else 517 file.write(surface.z[i][j]); 518 n++; 519 if (n > 5) { 520 n = 0; 521 file.write("\n"); 522 } 523 else { 524 file.write(" "); 525 } 526 } 527 } 528 } 529 530 531 532 unittest { 533 CartesianSurface surface = new CartesianSurface; 534 surface.loadFromZmap("./test/test_pet_rect.zmap"); 535 assert(surface.nx == 5); 536 assert(surface.ny == 3); 537 assert(surface.dx == 250); 538 assert(surface.dy == 500); 539 assert(surface.xOrigin == 5000); 540 assert(surface.yOrigin == 0); 541 assert(surface.z[0][0] == 1); 542 assert(surface.z[$ - 1][$ - 1] == 15); 543 544 //TODO test_pet_rect_blank.zmap 545 surface.loadFromZmap("./test/test_pet_sq.zmap"); 546 assert(surface.nx == 3); 547 assert(surface.ny == 3); 548 assert(surface.dx == 500); 549 assert(surface.dy == 500); 550 assert(surface.xOrigin == 5000); 551 assert(surface.yOrigin == 0); 552 assert(surface.z[0][0] == 0); 553 assert(surface.z[$ - 1][$ - 1] == 8); 554 555 surface.loadFromZmap("./test/test_rms_sq.zmap"); 556 assert(surface.nx == 3); 557 assert(surface.ny == 3); 558 assert(surface.dx == 500); 559 assert(surface.dy == 500); 560 assert(surface.xOrigin == 5000); 561 assert(surface.yOrigin == 0); 562 assert(surface.z[0][0] == 0); 563 assert(surface.z[$ - 1][$ - 1] == 8); 564 565 surface.loadFromZmap("./test/test_rms_rect.zmap"); 566 assert(surface.nx == 5); 567 assert(surface.ny == 3); 568 assert(surface.dx == 250); 569 assert(surface.dy == 500); 570 assert(surface.xOrigin == 5000); 571 assert(surface.yOrigin == 0); 572 assert(surface.z[0][0] == 1); 573 assert(surface.z[$ - 1][$ - 1] == 15); 574 } 575 576 577 /** 578 * Tries to detect surface format 579 * Params: 580 * fileName = 581 * Returns: string containing surface format. `cps` for CPS3 ASCII; 'zmap' for ZMAP+ ASCII; 'unknown' if format hasn't been detected. 582 */ 583 string surfaceFormat(string fileName) { 584 File file = File(fileName, "r"); 585 string format = "unknown"; 586 string str = file.readln(); 587 if (str.startsWith("FSASCI")) { 588 format = "cps"; 589 } 590 else { 591 if (str.startsWith("!")) { //probably zmap 592 while(str.startsWith("!")) { 593 str = file.readln(); 594 } 595 if (str.startsWith("@")) { 596 string [] words = str.replace(",", "").split(); 597 if (words.length >= 4 && icmp(words[2], "grid") == 0) 598 format = "zmap"; 599 } 600 } 601 } 602 file.close(); 603 return format; 604 } 605 606 unittest { 607 assert(surfaceFormat("./test/test_pet_rect.cps") == "cps"); 608 assert(surfaceFormat("./test/test_pet_rect_blank.cps") == "cps"); 609 assert(surfaceFormat("./test/test_pet_sq.cps") == "cps"); 610 assert(surfaceFormat("./test/test_rms_rect.cps") == "cps"); 611 assert(surfaceFormat("./test/test_rms_rect_blank.cps") == "cps"); 612 assert(surfaceFormat("./test/test_rms_sq.cps") == "cps"); 613 614 assert(surfaceFormat("./test/test_pet_rect.zmap") == "zmap"); 615 assert(surfaceFormat("./test/test_pet_rect_blank.zmap") == "zmap"); 616 assert(surfaceFormat("./test/test_pet_sq.zmap") == "zmap"); 617 assert(surfaceFormat("./test/test_rms_rect.zmap") == "zmap"); 618 assert(surfaceFormat("./test/test_rms_rect_blank.zmap") == "zmap"); 619 assert(surfaceFormat("./test/test_rms_sq.zmap") == "zmap"); 620 } 621 622 /** 623 * Samples height map to `surface` using height map from the given `source` 624 */ 625 void sampleFromSurface(CartesianSurface surface, CartesianSurface source) { 626 for (int i = 0; i < surface.nx; i++) { 627 for (int j = 0; j < surface.ny; j++) { 628 surface.z[i][j] = source.getZ(surface.xOrigin + i * surface.dx, surface.yOrigin + j * surface.dy); 629 } 630 } 631 } 632 633 /** 634 * Translates the given surface 635 * Params: 636 * surface = surface to translate 637 * dx = translation value along X direction 638 * dy = translation value along Y direction 639 * Returns: translated surface for more convenient call chaining 640 */ 641 CartesianSurface translate(CartesianSurface surface, double dx, double dy) { 642 surface.m_xOrigin += dx; 643 surface.m_yOrigin += dy; 644 return surface; 645 } 646 647 /** 648 * Scales the given surface around it's origin point 649 * Params: 650 * surface = surface to scale 651 * xf = scale factor along X direction 652 * xf = scale factor value along Y direction 653 * Returns: translated surface for more convenient call chaining 654 */ 655 CartesianSurface scale(CartesianSurface surface, double xf, double yf) { //scales around origin point 656 //TODO filter negative factors 657 surface.m_dx *= xf; 658 surface.m_dy *= yf; 659 return surface; 660 } 661 662 CartesianSurface normalize(CartesianSurface surface) { 663 immutable double zmax = surface.m_z[surface.m_z.maxIndex]; 664 immutable double zmin = surface.m_z[surface.m_z.minIndex]; 665 surface.m_z[] -= zmin; 666 surface.m_z[] /= zmax; 667 return surface; 668 } 669 670 unittest { 671 //TODO implement 672 } 673 674 //TODO flipAnlogI/flipAlongJ 675