1 module dsurf.cartesian; 2 3 import std.range; 4 import std..string; 5 import std.algorithm; 6 import std.conv; 7 import std.math; 8 9 /*********************************** 10 * Class CartesianSurface represents regular rectangular-cell surface 11 * Authors: Dmitriy Linev 12 * License: MIT 13 * Date: 2019-2021 14 */ 15 class CartesianSurface { 16 17 /** 18 * Default constructor, doesn`t allocate memory for height map 19 */ 20 this() pure { 21 m_xOrigin = 0; 22 m_yOrigin = 0; 23 m_dx = 0; 24 m_dy = 0; 25 m_nx = 0; 26 m_ny = 0; 27 } 28 29 /** 30 * Sets surface header and allocates surface memory 31 * Params: 32 * nx = number of nodes along X coordinate 33 * ny = number of nodes along Y coordinate 34 * xOrigin = surface origin X coordinate 35 * yOrigin = surface origin Y coordinate 36 * dx = surface increment (cell size) along X coordinate 37 * dy = surface increment (cell size) along Y coordinate 38 */ 39 void setHeader(int nx, int ny, double xOrigin, double yOrigin, double dx, double dy) pure { 40 if (nx <= 0 || ny <= 0) 41 throw new Exception("How do we suppose to set negative size of the surface?"); 42 m_nx = nx; 43 m_ny = ny; 44 m_xOrigin = xOrigin; 45 m_yOrigin = yOrigin; 46 m_dx = dx; 47 m_dy = dy; 48 m_zd = new double[](m_nx * m_ny); 49 m_z = m_zd.chunks(m_ny).array; 50 } 51 52 /** 53 * Assigns each surface note to a constant height value 54 * Params: 55 * z = height 56 */ 57 void assignConstant(double z) { 58 m_zd[] = z; 59 } 60 61 /// Copy constructor, returns the exact copy of the given surface 62 this(CartesianSurface surface) pure { 63 this.setHeader(surface.nx, surface.ny, surface.xOrigin, surface.yOrigin, surface.dx, surface.dy); 64 this.m_zd[] = surface.m_zd[]; 65 } 66 67 68 69 /// Returns X coordinate of surface origin 70 @property pure double xOrigin() const @safe @nogc { return m_xOrigin; } 71 72 /// Returns Y coordinate of surface origin 73 @property pure double yOrigin() const @safe @nogc { return m_yOrigin; } 74 75 /// Returns surface increment (cell size) along X axis 76 @property pure double dx() const @safe @nogc { return m_dx; } 77 78 /// Returns surface increment (cell size) along Y axis 79 @property pure double dy() const @safe @nogc { return m_dy; } 80 81 /// Returns number of nodes along X axis 82 @property pure int nx() const @safe @nogc { return m_nx; } 83 84 /// Returns number of nodes along Y axis 85 @property pure int ny() const @safe @nogc { return m_ny; } 86 87 /// Returns minimum z value 88 @property pure double zMin() const @safe @nogc { return m_zd[m_zd.minIndex]; } 89 90 /// Returns maximum z value 91 @property pure double zMax() const @safe @nogc { return m_zd[m_zd.maxIndex]; } 92 93 /** 94 Method to access height map. 95 Returns: `Slice!(double*, 2)` containing surface`s height map with dimensions nx * ny 96 Example: 97 --- 98 foreach (i; 0 .. surface.nx) { 99 foreach (j; 0 .. surface.ny) { 100 surface.z[i][j] = 0; 101 } 102 } 103 --- 104 */ 105 @property double[][] z() { return m_z; } 106 107 /// Returns `true` if point with given coordinates is inside surface boundaries, otherwise returns `false` 108 bool isInsideSurface(double x, double y) const @nogc { 109 if (x < m_xOrigin || y < m_yOrigin || x > m_xOrigin + m_dx * (m_nx - 1) || y > m_yOrigin + m_dy * (m_ny - 1)) 110 return false; 111 return true; 112 } 113 114 /// Returns cell index from a given `X` coordinate. Maximum number is `nx - 1` 115 /// Returns -1 if a given coordinate is outside the surface 116 int cellXIndex(double x) const { 117 if (x < m_xOrigin || x > m_xOrigin + m_dx * (m_nx - 1)) 118 return -1; //TODO throw? 119 return ((x - m_xOrigin) / m_dx).to!int; 120 } 121 122 /// Returns cell index from a given `Y` coordinate. Maximum number is `ny - 1` 123 /// Returns `-1` if a given coordinate is outside the surface 124 int cellYIndex(double y) const { 125 if (y < m_yOrigin || y > m_yOrigin + m_dy * (m_ny - 1)) 126 return -1; //TODO throw? 127 return ((y - m_yOrigin) / m_dy).to!int; 128 } 129 130 /// Returns z value of a point with given coordinate using bilinear interpolation 131 double getZ(double x, double y) { 132 if (!isInsideSurface(x, y)) 133 return double.nan; 134 immutable int i = cellXIndex(x); 135 immutable int j = cellYIndex(y); 136 137 if (i == m_nx - 1 && j == m_ny - 1) //top right corner 138 return m_z[nx - 1][ny - 1]; 139 else if (i == m_nx - 1) { //right edge 140 return m_z[i][j] + (m_z[i][j + 1] - m_z[i][j]) / m_dy * (y - (m_yOrigin + j * m_dy)); 141 } 142 else if (j == m_ny - 1) { //top edge 143 return m_z[i][j] + (m_z[i + 1][j] - m_z[i][j]) / m_dx * (x - (m_xOrigin + i * m_dx)); 144 } 145 146 //z: top left, top right, bottom left, bottom right 147 double [] zcells = [m_z[i][j + 1], m_z[i + 1][j + 1], m_z[i][j], m_z[i + 1][j]]; 148 149 immutable ulong blanks = zcells.count!isNaN; 150 if (blanks > 1) { // more than one node is blank - value in not defined 151 return double.nan; 152 } 153 else if (blanks == 1) { // one node is blank - tryng to reconstruct actual value 154 if (isNaN(zcells[0])) { 155 if (sqrt(pow(x - (m_xOrigin + i * m_dx), 2) + pow(y - (m_yOrigin + (j + 1) * m_dy), 2)) < 156 sqrt(pow(x - (m_xOrigin + (i + 1) * m_dx), 2) + pow(y - (m_yOrigin + j * m_dy), 2))) 157 // checking if given point is inside defined triangle 158 // by comparing a distance to the blank corner with a distance to the opposite corner 159 return double.nan; 160 // if it is defined, reconstruct the missing corner assuming the fixed gradient along two opposite edges 161 zcells[0] = zcells[2] + zcells[1] - zcells[3]; 162 } 163 else if (isNaN(zcells[1])) { 164 if (sqrt(pow(x - (m_xOrigin + (i + 1) * m_dx), 2) + pow(y - (m_yOrigin + (j + 1) * m_dy), 2)) < 165 sqrt(pow(x - (m_xOrigin + i * m_dx), 2) + pow(y - (m_yOrigin + j * m_dy), 2))) 166 return double.nan; 167 zcells[1] = zcells[3] + zcells[0] - zcells[2]; 168 } 169 else if (isNaN(zcells[2])) { 170 if (sqrt(pow(x - (m_xOrigin + i * m_dx), 2) + pow(y - (m_yOrigin + j * m_dy), 2)) < 171 sqrt(pow(x - (m_xOrigin + (i + 1) * m_dx), 2) + pow(y - (m_yOrigin + (j + 1) * m_dy), 2))) 172 return double.nan; 173 zcells[2] = zcells[0] + zcells[3] - zcells[1]; 174 } 175 else if (isNaN(zcells[3])) { 176 if (sqrt(pow(x - (m_xOrigin + (i + 1) * m_dx), 2) + pow(y - (m_yOrigin + j * m_dy), 2)) < 177 sqrt(pow(x - (m_xOrigin + i * m_dx), 2) + pow(y - (m_yOrigin + (j + 1) * m_dy), 2))) 178 return double.nan; 179 zcells[3] = zcells[1] + zcells[2] - zcells[0]; 180 } 181 } 182 // bilinear interpolation 183 immutable double z1 = (m_xOrigin + (i + 1) * m_dx - x) / dx * zcells[2] + 184 (x - (m_xOrigin + i * m_dx)) / dx * zcells[3]; 185 immutable double z2 = (m_xOrigin + (i + 1) * m_dx - x) / dx * zcells[0] + 186 (x - (m_xOrigin + i * m_dx)) / dx * zcells[1]; 187 return (m_yOrigin + (j + 1) * m_dy - y) / m_dy * z1 + 188 (y - (m_yOrigin + j * m_dy)) / m_dy * z2; 189 } 190 191 /// Returns a minimum height value in a surface 192 double min() const { 193 return m_zd[this.m_zd.minIndex]; 194 } 195 196 /// Returns a maximum height value in a surface 197 double max() const { 198 return m_zd[this.m_zd.maxIndex]; 199 } 200 201 /// Operators +=, -=, *=, /= overloading for two surfaces 202 CartesianSurface opOpAssign(string op)(CartesianSurface rhs) { 203 for (int i = 0; i < m_nx; i++) { 204 for (int j = 0; j < m_ny; j++) { 205 immutable double x = m_xOrigin + i * dx; 206 immutable double y = m_yOrigin + j * dy; 207 immutable double rhsz = rhs.getZ(x, y); 208 if (isNaN(rhsz)) 209 continue; 210 static if (op == "+") 211 m_z[i][j] += rhsz; 212 else static if (op == "-") 213 m_z[i][j] -= rhsz; 214 else static if (op == "*") 215 m_z[i][j] *= rhsz; 216 else static if (op == "/") { 217 if (abs(rhsz) < 1e-9) 218 m_z[i][j] = double.nan; 219 else 220 m_z[i][j] /= rhsz; 221 } 222 else static assert(0, "Operator " ~ op ~ " not implemented"); 223 } 224 } 225 return this; 226 } 227 228 229 230 /// Operators +, -, *, / overloading for two surfaces 231 CartesianSurface opBinary(string op)(CartesianSurface rhs) { 232 CartesianSurface result = new CartesianSurface(this); 233 static if (op == "+") 234 result += rhs; 235 else static if (op == "-") 236 result -= rhs; 237 else static if (op == "*") 238 result *= rhs; 239 else static if (op == "/") 240 result /= rhs; 241 else static assert(0, "Operator "~op~" not implemented"); 242 return result; 243 } 244 245 246 247 /// Operators +=, -=, *=, /= overloading for a surface and a fixed value 248 CartesianSurface opOpAssign(string op)(double rhs) { 249 static if (op == "+") 250 m_zd [] += rhs; 251 else static if (op == "-") 252 m_zd [] -= rhs; 253 else static if (op == "*") 254 m_zd [] *= rhs; 255 else static if (op == "/") 256 m_zd [] /= rhs; 257 else static assert(0, "Operator "~op~" not implemented"); 258 return this; 259 } 260 261 /// Operators +, -, *, / overloading for a surface and a fixed value 262 CartesianSurface opBinary(string op)(double rhs) { 263 CartesianSurface result = new CartesianSurface(this); 264 static if (op == "+") 265 result += rhs; 266 else static if (op == "-") 267 result -= rhs; 268 else static if (op == "*") 269 result *= rhs; 270 else static if (op == "/") 271 result /= rhs; 272 else static assert(0, "Operator "~op~" not implemented"); 273 return result; 274 } 275 276 277 278 package: 279 double[] m_zd; /// dense representation of Z values of the surface 280 double[][] m_z; /// Z chunks 281 double m_xOrigin; 282 double m_yOrigin; 283 double m_dx; 284 double m_dy; 285 int m_nx; 286 int m_ny; 287 } 288 289 /** 290 * Samples height map to `surface` using height map from the given `source` 291 */ 292 void sampleFromSurface(CartesianSurface surface, CartesianSurface source) { 293 foreach (i; 0 .. surface.nx) { 294 foreach (j; 0 .. surface.ny) { 295 surface.z[i][j] = source.getZ(surface.xOrigin + i * surface.dx, surface.yOrigin + j * surface.dy); 296 } 297 } 298 } 299 300 /** 301 * Translates the given surface (moves its origin point to the given vector) leaving increments untouched 302 * Params: 303 * surface = surface to translate 304 * dx = translation value along X direction 305 * dy = translation value along Y direction 306 * Returns: translated surface for more convenient call chaining 307 */ 308 CartesianSurface translate(CartesianSurface surface, double dx, double dy) { 309 surface.m_xOrigin += dx; 310 surface.m_yOrigin += dy; 311 return surface; 312 } 313 314 //TODO flipAnlogI/flipAlongJ 315