00001 #include "FlowGeometry.h"
00002 #include "reverseBytes.h"
00003
00004 bool FlowGeometry::readFromFile(char* header, FILE* fp, bool bigEndian)
00005 {
00006 isFlipped = false;
00007
00008 sscanf(header,"SN4DB %d %d %d",&dim[0],&dim[1],&dim[2]);
00009 std::cout << "Dimensions: " << dim[0] << " x " << dim[1] << " x " << dim[2] << std::endl;
00010 if (dim[2]!=1)
00011 {
00012 std::cerr << "Invalid Z dimension value." << std::endl;
00013 return false;
00014 }
00015
00016 geometryData = new vec3[dim[0]*dim[1]];
00017
00018 int result = fread(geometryData,sizeof(vec3),dim[0]*dim[1],fp);
00019 if (result != dim[0]*dim[1])
00020 {
00021 std::cerr << "+ Error reading grid file." << std::endl << std::endl;
00022 return false;
00023 }
00024
00025 if (bigEndian)
00026 for (int j = 0; j < getDimX()*getDimY(); j++)
00027 for (int k = 0; k < 3; k++)
00028 geometryData[j][k] = reverseBytes<float>(geometryData[j][k]);
00029
00030
00031 boundaryMin = vec3(getPos(0));
00032
00033 boundaryMax = vec3(getPos((dim[0]*dim[1]) - 1));
00034 boundarySize = boundaryMax - boundaryMin;
00035
00036
00037 if (getPosY(dim[0]-1)>(boundaryMin[1] + boundaryMax[1])*0.5)
00038 {
00039
00040 isFlipped = true;
00041 std::cout << "Flipped Y and X dimensions." << std::endl;
00042 }
00043 else isFlipped = false;
00044
00045 std::cout << "X Boundaries: " << boundaryMin[0] << " ... " << boundaryMax[0] << std::endl;
00046 std::cout << "Y Boundaries: " << boundaryMin[1] << " ... " << boundaryMax[1] << std::endl;
00047
00048 return true;
00049 }
00050
00051 FlowGeometry::FlowGeometry()
00052 {
00053 geometryData = NULL;
00054 isFlipped = false;
00055 }
00056
00057 FlowGeometry::~FlowGeometry()
00058 {
00059 if (geometryData)
00060 delete[] geometryData;
00061 }
00062
00064 int FlowGeometry::getXYvtx(vec3 pos)
00065 {
00066 int i;
00067 int j;
00068
00069 for (i = 0; (i < dim[0])&&(getPosX(getVtx(i,0)) < pos[0]); i++);
00070
00071 for (j = 0; (j < dim[1])&&(getPosY(getVtx(0,j)) < pos[1]); j++);
00072
00073
00074 return getVtx((i<dim[0]) ? i : dim[0]-1, (j<dim[1]) ? j : dim[1]-1);
00075 }
00076
00077 inline vec3 FlowGeometry::getPos(int vtxID)
00078 {
00079 return geometryData[vtxID];
00080 }
00081
00082 inline float FlowGeometry::getPosX(int vtxID)
00083 {
00084 return geometryData[vtxID][0];
00085 }
00086
00087 inline float FlowGeometry::getPosY(int vtxID)
00088 {
00089 return geometryData[vtxID][1];
00090 }
00091
00092
00094 int FlowGeometry::getNearestVtx(vec3 pos)
00095 {
00096
00097 float dist = pos.dist2(getPos(0));
00098
00099 int closest = 0;
00100
00101 int next = getRightNeigh(0);
00102 float newd = pos.dist2(getPos(next));
00103
00104
00105
00106 while (newd <= dist)
00107 {
00108 closest = next;
00109 dist = newd;
00110 next = getRightNeigh(closest);
00111 newd = pos.dist2(getPos(next));
00112 }
00113
00114 next = getBottomNeigh(closest);
00115 newd = pos.dist2(getPos(next));
00116
00117
00118
00119 while (newd <= dist)
00120 {
00121 closest = next;
00122 dist = newd;
00123 next = getBottomNeigh(closest);
00124 newd = pos.dist2(getPos(next));
00125 }
00126 return closest;
00127 }
00128
00129 bool FlowGeometry::getInterpolationAt(vec3 pos, int* vtxID, float* coef)
00130 {
00131
00132
00133 if ((pos[0]<boundaryMin[0])||(pos[1]<boundaryMin[1])||(pos[0]>boundaryMax[0])||(pos[1]>boundaryMax[1]))
00134 return false;
00135
00136
00137 vtxID[0] = getNearestVtx(pos);
00138
00139
00140
00141
00142 float dX = pos[0] - getPos(vtxID[0])[0];
00143 float dY = pos[1] - getPos(vtxID[0])[1];
00144
00145 if (dX > 0) {
00146 vtxID[1] = getRightNeigh(vtxID[0]);
00147 } else {
00148 vtxID[1] = getLeftNeigh(vtxID[0]);
00149
00150 if (vtxID[1] < 0) vtxID[1] = getRightNeigh(vtxID[0]);
00151 }
00152
00153 if (dY > 0) {
00154 vtxID[2] = getBottomNeigh(vtxID[0]);
00155 vtxID[3] = getBottomNeigh(vtxID[1]);
00156 } else {
00157 vtxID[2] = getTopNeigh(vtxID[0]);
00158 vtxID[3] = getTopNeigh(vtxID[1]);
00159 if (vtxID[2] < 0) {
00160 vtxID[2] = getBottomNeigh(vtxID[0]);
00161 vtxID[3] = getBottomNeigh(vtxID[1]);
00162 }
00163 }
00164
00165
00166
00167 float gridDiffX = getPos(vtxID[1])[0]-getPos(vtxID[0])[0];
00168 float gridDiffY = getPos(vtxID[2])[1]-getPos(vtxID[0])[1];
00169
00170 float ratio[4];
00171 ratio[0] = (gridDiffX-dX)/gridDiffX;
00172 ratio[1] = dX/gridDiffX;
00173 ratio[2] = (gridDiffY-dY)/gridDiffY;
00174 ratio[3] = dY/gridDiffY;
00175 coef[0] = ratio[0]*ratio[2];
00176 coef[1] = ratio[1]*ratio[2];
00177 coef[2] = ratio[0]*ratio[3];
00178 coef[3] = ratio[1]*ratio[3];
00179
00180 return true;
00181 }
00182
00183 inline float FlowGeometry::getMinX()
00184 {
00185 return boundaryMin[0];
00186 }
00187
00188 inline float FlowGeometry::getMaxX()
00189 {
00190 return boundaryMax[0];
00191 }
00192
00193 inline float FlowGeometry::getMinY()
00194 {
00195 return boundaryMin[1];
00196 }
00197
00198 inline float FlowGeometry::getMaxY()
00199 {
00200 return boundaryMax[1];
00201 }
00202
00203 inline int FlowGeometry::getDimX()
00204 {
00205 return (isFlipped) ? dim[1] : dim[0];
00206 }
00207
00208 inline int FlowGeometry::getDimY()
00209 {
00210 return (isFlipped) ? dim[0] : dim[1];
00211 }
00212
00213 inline int FlowGeometry::getDimZ()
00214 {
00215 return dim[2];
00216 }
00217
00218 int FlowGeometry::getVtx(int x, int y)
00219 {
00220
00221 return (isFlipped)? (x*dim[1]) + y : (y*dim[0]) + x;
00222 }
00223
00224 int FlowGeometry::getVtxX(int vtxID)
00225 {
00226
00227 return (isFlipped)? vtxID / dim[1] : vtxID % dim[0];
00228 }
00229
00230 int FlowGeometry::getVtxY(int vtxID)
00231 {
00232
00233 return (isFlipped)? vtxID % dim[1] : vtxID / dim[0];
00234 }
00235
00236 int FlowGeometry::getRightNeigh(int vtxID)
00237 {
00238 int x = getVtxX(vtxID);
00239 return (x+1 < dim[0]) ? getVtx(x+1,getVtxY(vtxID)) : -1;
00240 }
00241
00242 int FlowGeometry::getTopNeigh(int vtxID)
00243 {
00244
00245
00246 int y = getVtxY(vtxID);
00247 return (y > 1) ? getVtx(getVtxX(vtxID), y-1) : -1;
00248 }
00249
00250 int FlowGeometry::getLeftNeigh(int vtxID)
00251 {
00252 int x = getVtxX(vtxID);
00253 return (x > 1) ? getVtx(x-1,getVtxY(vtxID)) : -1;
00254 }
00255
00256 int FlowGeometry::getBottomNeigh(int vtxID)
00257 {
00258
00259 int y = getVtxY(vtxID);
00260 return (y+1 < dim[1]) ? getVtx(getVtxX(vtxID), y+1) : -1;
00261 }
00262
00263 vec3 FlowGeometry::normalizeCoords(vec3 pos)
00264 {
00265 vec3 u(pos - boundaryMin);
00266
00267 u[0] /= boundarySize[0];
00268 u[1] /= boundarySize[1];
00269
00270 return u;
00271 }
00272
00273 vec3 FlowGeometry::unNormalizeCoords(vec3 pos)
00274 {
00275 vec3 u(pos);
00276
00277 u[0] *= boundarySize[0];
00278 u[1] *= boundarySize[1];
00279
00280 return u += boundaryMin;
00281 }
00282
00283 bool FlowGeometry::getIsFlipped()
00284 {
00285 return isFlipped;
00286 }