00001 /*
00002 File: MatRad.cc
00003
00004 Function: See header file
00005
00006 Author(s): Andrew Willmott
00007
00008 Copyright: (c) 1997-2000, Andrew Willmott
00009
00010 Notes:
00011
00012 */
00013
00014 #include "MatRad.h"
00015 #include "cl/Timer.h"
00016 #include "gcl/Draw.h"
00017 #include <fstream.h>
00018 #include <stdio.h>
00019
00020 Matd MatRad::EmissionVectors(PatchList &patches)
00021 {
00022 Matd result(3, patches.NumItems());
00023 Int i, j;
00024
00025 for (i = 0; i < patches.NumItems(); i++)
00026 for (j = 0; j < 3; j++)
00027 result[j][i] = patches[i]->Emittance()[j];
00028
00029 return(result);
00030 }
00031
00032 SparseVecd MatRad::FormFactorToVector(Int j, PatchList &patches)
00033 // Find the vector of patch factors to patch j
00034 {
00035 const GCLReal epsilon = 1e-6;
00036 SparseVecd result(patches.NumItems());
00037 Int i;
00038 GCLReal vis;
00039
00040 #ifdef RAD_VIS
00041 if (gRadControl->showRays && gRadControl->visibility == vis_1)
00042 display->Clear().Draw(GetScene());
00043 #endif
00044
00045 if (len(patches[j]->Reflectance()) <= 1e-6)
00046 result.MakeUnit(j);
00047 else
00048 {
00049 result.Begin();
00050 for (i = 0; i < patches.NumItems(); i++)
00051 if (i != j)
00052 {
00053 vis = patches[i]->Visibility(patches[j]);
00054
00055 if (vis > 0.0)
00056 result.AddElt(i, patches[i]->EstFormFactor(patches[j])
00057 * vis);
00058 }
00059 result.End();
00060 }
00061
00062 #ifdef RAD_VIS
00063 if (gRadControl->showRays && gRadControl->visibility == vis_1)
00064 {
00065 display->Show();
00066 Pause();
00067 }
00068 #endif
00069
00070 return(result);
00071 }
00072
00073 Bool MatRad::Render()
00074 // Returns true if no interruptions.
00075 {
00076 Int i, j;
00077 Vec3d error;
00078 Int numPolys = patches.NumItems();
00079 Colour c;
00080
00081 RadMethod::Render();
00082 SparseVecd::SetFuzz(1e-10);
00083
00084 SparseVecd FFRow(numPolys); // Row of form factors
00085
00086 for (i = 0; i < 3; i++)
00087 A[i].SetSize(numPolys, numPolys);
00088
00089 E.SetSize(3, numPolys);
00090 B.SetSize(3, numPolys);
00091
00092 numPatches = 0;
00093
00094 // Flat colour scene according to current reflectances
00095 ColourMeshInitial();
00096
00097 if (Stage(1)) return(0);
00098
00099 E = EmissionVectors(patches);
00100
00101 if (Stage(2)) return(0);
00102
00103 // Initialise A's to unit matrices
00104 for (j = 0; j < 3; j++)
00105 A[j] = vl_I;
00106
00107 if (Stage(3)) return(0);
00108
00109 // Calculate solution matrices
00110 for (i = 0; i < patches.NumItems(); i++)
00111 {
00112 numPatches++;
00113 // find form factors to this patch
00114 FFRow = FormFactorToVector(i, patches);
00115
00116 for (j = 0; j < 3; j++)
00117 {
00118 A[j][i] -= FFRow * patches[i]->Reflectance()[j];
00119
00120 c[j] = patches[i]->Emittance()[j] +
00121 patches[i]->Reflectance()[j] * dot(FFRow, E[j]);
00122 }
00123 patches[i]->SetColour(c);
00124
00125 if (Stage(4)) return(0);
00126 }
00127
00128 if (Stage(5)) return(0);
00129
00130 error = vl_1;
00131
00132 while (len(error) > 0.001)
00133 {
00134 if (!gRadControl->useConjGrad)
00135 {
00136 for (i = 0; i < 3; i++) // Solve one more iteration
00137 error[i] = SolveOverRelax(A[i], B[i], E[i], 1e-6,
00138 gRadControl->alpha);
00139 }
00140 else
00141 for (i = 0; i < 3; i++) // Solve one more iteration
00142 error[i] = SolveConjGrad(A[i], B[i], E[i], 1e-6);
00143
00144 if (Stage(6)) return(0);
00145 }
00146
00147 if (Stage(7)) return(0);
00148
00149 return(1);
00150 }
00151
00152 #ifdef RAD_VIS
00153
00154 Void MatRad::DumpStats() {}
00155
00156 Int MatRad::Stage(Int stage)
00157 {
00158 Bool animate = gRadControl->animate || gRadControl->pause;
00159 Int i, j;
00160 Colour c;
00161
00162 switch (stage)
00163 {
00164 case 1: // pre setup
00165 Field(out3) << patches.NumItems() << " patches." << show;
00166 Field(out1) << "Forming emission vector E..." << show;
00167 break;
00168
00169 case 2: // post setup
00170 Field(out1) << "Flat-colouring scene..." << show;
00171
00172 for (i = 0; i < patches.NumItems(); i++) // distribute colour
00173 {
00174 for (j = 0; j < 3; j++)
00175 c[j] = patches[i]->Reflectance()[j]
00176 * Max((double) 0.2, E[j][i]);
00177
00178 patches[i]->SetColour(c);
00179 }
00180
00181 ColourVertices();
00182 display->Redraw(); // draw
00183
00184 Field(out1) << "1. Initialise F matrices..." << show;
00185 break;
00186
00187 case 3:
00188 Field(out1) << "2. Calculate F matrices..." << show;
00189 StartUpdate();
00190 break;
00191
00192 case 4: // After patch setup
00193 doUpdate = animate || Update();
00194
00195 if (doUpdate)
00196 {
00197 if (gRadControl->animate && numPatches < patches.NumItems())
00198 patches[numPatches]->SetHighlight(1);
00199 ColourVertices();
00200 display->Redraw();
00201 if (gRadControl->animate && numPatches < patches.NumItems())
00202 patches[numPatches]->SetHighlight(0);
00203
00204 RenderMatrix();
00205 Field(out3) << numPatches << " / "
00206 << patches.NumItems() << " patches." << show;
00207
00208 UpdateCont();
00209 }
00210 break;
00211
00212 case 5: // After FF calculation
00213 ColourVertices();
00214 display->Redraw();
00215 RenderMatrix();
00216 Field(out1) << "3. Solve equations..." << show;
00217 StartUpdate();
00218 solverIterations = 0;
00219 break;
00220
00221 case 6: // Middle of solve loop.
00222 solverIterations++;
00223 if (animate || Update())
00224 {
00225 // Colour scene according to current B vector
00226 for (i = 0; i < patches.NumItems(); i++)
00227 {
00228 for (j = 0; j < 3; j++)
00229 c[j] = B[j][i];
00230
00231 patches[i]->SetColour(c);
00232 }
00233
00234 ColourVertices();
00235 display->Redraw(); // Draw!
00236 Field(out3) << solverIterations << " iterations" << show;
00237 UpdateCont();
00238 }
00239 break;
00240
00241 case 7:
00242 ColourVertices();
00243 display->Redraw();
00244
00245 Field(out3) << "Solving the matrix took "
00246 << solverIterations << " iterations" << show;
00247 Field(out1) << "Finished!" << show;
00248 break;
00249 }
00250
00251 if (Idle()) return(1);
00252
00253 return(Pause());
00254 }
00255
00256 #else
00257
00258 Void MatRad::DumpMatrix()
00259 {
00260 Int i, j, k;
00261 Char *rgb = "rgb";
00262 String temp;
00263 FILE *file;
00264
00265 for (k = 0; k < 3; k++)
00266 {
00267 temp.Printf("%s.%c.mlab", gRadControl->outFile, rgb[k]);
00268 file = fopen(temp, "w");
00269 fprintf(file, "[\n");
00270 for (i = 0; i < A[k].Rows(); i++)
00271 {
00272 for (j = 0; j < A[k].Cols() - 1; j++)
00273 fprintf(file, "%.18g, ", A[k][i][j]);
00274 fprintf(file, "%.18g\n", A[k][i][j]);
00275 }
00276
00277 fprintf(file, "]\n");
00278 fclose(file);
00279 }
00280
00281 temp.Printf("%s.stats", gRadControl->outFile);
00282 file = fopen(temp, "w");
00283 fprintf(file, "minArea %.18g\n", stats.minArea);
00284 fprintf(file, "maxArea %.18g\n", stats.maxArea);
00285 fprintf(file, "minRefl [%.18g %.18g %.18g]\n",
00286 stats.minRefl[0], stats.minRefl[1], stats.minRefl[2]);
00287 fprintf(file, "maxRefl [%.18g %.18g %.18g]\n",
00288 stats.maxRefl[0], stats.maxRefl[1], stats.maxRefl[2]);
00289 fclose(file);
00290 }
00291
00292 Void MatRad::DumpStats()
00293 {
00294 GCLReal density, mem;
00295 Int i, j;
00296
00297 density = 0;
00298 for (i = 0; i < 3; i++)
00299 for (j = 0; j < A[i].Rows(); j++)
00300 density += A[i][j].NumItems();
00301
00302 mem = (density * (sizeof(GCLReal) + sizeof(Int))) +
00303 sizeof(GCLReal) * E.Rows() * E.Cols() +
00304 sizeof(GCLReal) * B.Rows() * B.Cols();
00305 mem /= 1024.0;
00306
00307 density /= 3 * A[0].Rows() * A[0].Cols();
00308
00309 cout << dumpID
00310 << ' ' << totTime
00311 << ' ' << gRadControl->stage
00312 << ' ' << numPatches
00313 << ' ' << patches.NumItems()
00314 << ' ' << gRadControl->rays
00315 << ' ' << density
00316 << ' ' << mem
00317 << ' ' << grid->MemoryUsage()
00318 << ' ' << TotalMemoryUse()
00319 << endl;
00320
00321 DumpScene();
00322 }
00323
00324 Int MatRad::Stage(Int stage)
00325 {
00326 Colour c;
00327
00328 if (CheckTime()) return(1);
00329 // CheckTime calls DumpStats if enough time has passed, and
00330 // also stops the clock.
00331
00332 gRadControl->stage = stage;
00333
00334 switch (stage)
00335 {
00336 case 1: // pre setup
00337 cout << "method mat" << endl;
00338 if (!gRadControl->useConjGrad)
00339 cout << "solver overrelaxation" << endl;
00340 else
00341 cout << "solver conjugate-gradient" << endl;
00342
00343 cout << "sub " << gRadControl->patchSubdivs << endl;
00344 cout << "alpha " << gRadControl->alpha << endl;
00345 cout << "scene " << sceneName << endl;
00346 cout << "polys " << gRadControl->numPolys << endl;
00347 cout << "format "
00348 << "ID "
00349 << "time "
00350 << "stage "
00351 << "procPatches "
00352 << "patches "
00353 << "rays "
00354 << "density "
00355 << "mem "
00356 << "rtMem "
00357 << "rss "
00358 << endl;
00359 cout << "----------------------------------------------------------"
00360 << endl;
00361
00362 DumpStats();
00363 break;
00364
00365 case 5: // Finished calculating A[k]
00366 if (gRadControl->drawMatrix)
00367 DumpMatrix();
00368 break;
00369
00370 case 6: // Middle of solve loop.
00371 #ifndef RAD_BEST_TIMING
00372 Int i, j;
00373
00374 // Colour scene according to current B vector
00375 for (i = 0; i < patches.NumItems(); i++)
00376 {
00377 for (j = 0; j < 3; j++)
00378 c[j] = B[j][i];
00379
00380 patches[i]->SetColour(c);
00381 }
00382 #endif
00383 break;
00384
00385 case 7:
00386 DumpStats();
00387 break;
00388 }
00389
00390 if (Idle()) return(0);
00391 timer.ContTimer();
00392 return(0);
00393 }
00394
00395 #endif
00396
00397 Void MatRad::RemoveDirect()
00398 {
00399 Int i;
00400
00401 }
00402
00403 Void MatRad::DrawMatrix(Renderer &r)
00404 {
00405 #ifdef RAD_VIS
00406 SVIterd iter0, iter1, iter2;
00407
00408 if (!gRadControl->drawMatrix || gRadControl->stop)
00409 return;
00410
00411 Int i, j, n;
00412 GCLReal x1, x2, y1, y2, w;
00413
00414 w = 2.0 / A[0].Cols();
00415
00416 y1 = 1 - w;
00417 y2 = 1;
00418 n = Min(A[0].Rows(), (Int) numPatches);
00419
00420 for (i = 0; i < n; i++)
00421 {
00422 x1 = -1;
00423 x2 = w - 1;
00424
00425 iter0.Begin(A[0][i]);
00426 iter1.Begin(A[1][i]);
00427 iter2.Begin(A[2][i]);
00428
00429 for (j = 0; j < A[0].Cols(); j++)
00430 {
00431 iter0.Inc(j);
00432 iter1.Inc(j);
00433 iter2.Inc(j);
00434
00435 if (patches[i]->highlight == 2 && patches[j]->highlight == 3)
00436 r.C(cPurple);
00437 else if (patches[i]->highlight == 2)
00438 r.C(cYellow);
00439 else if (patches[j]->highlight == 3)
00440 r.C(cGreen);
00441 else if (i == j)
00442 r.C(cBlack);
00443 else
00444 {
00445 Colour mClr = cBlack;
00446
00447 if (iter0.Exists(j)) mClr[0] = -iter0.Data();
00448 if (iter1.Exists(j)) mClr[1] = -iter1.Data();
00449 if (iter2.Exists(j)) mClr[2] = -iter2.Data();
00450
00451 r.C(A[0].Cols() * mClr);
00452 }
00453
00454 PaintRect(r, Coord(x1, y1), Coord(x2, y2));
00455
00456 x1 += w;
00457 x2 += w;
00458 }
00459
00460 y1 -= w;
00461 y2 -= w;
00462 }
00463 #endif
00464 }
00465
00466 RadElem *MatRad::NewMesh()
00467 {
00468 return(new RadGrid());
00469 }
