diff --git a/cmath3d/TriangleMesh.cpp b/cmath3d/TriangleMesh.cpp index bcf0e6a4f59a9de8144d407574b0fcbfbdb94d6a..834e67678643aaa1dcd9db4c047256f4691598a8 100644 --- a/cmath3d/TriangleMesh.cpp +++ b/cmath3d/TriangleMesh.cpp @@ -2,11 +2,15 @@ #include <fstream> #include <map> #include <lapacke.h> + #include <i3d/draw.h> #include <i3d/morphology.h> +#include <i3d/DistanceTransform.h> +#include <i3d/filters.h> #include "TriangleMesh.h" #include "../cmath3d_v/TriangleMesh_v.h" +#include "../src/rnd_generators.h" #undef min #undef max @@ -834,6 +838,63 @@ void ActiveMesh::RenderMaskB(i3d::Image3d<i3d::GRAY16>& mask) } +void ActiveMesh::RenderOneTimeTexture(const i3d::Image3d<i3d::GRAY16>& mask, + i3d::Image3d<i3d::GRAY16>& texture) +{ + i3d::Image3d<float> dt; + i3d::GrayToFloat(mask,dt); + i3d::EDM(dt,0,100.f,false); + + dt.SaveImage("dt.ics"); + + float* p=dt.GetFirstVoxelAddr(); + float* const pL=p+dt.GetImageSize(); + int randomCounter=20; + while (p != pL) + { + //are we on a surface shell? + if (*p < 0.10f) + { + //invert the intensity levels (iso-lines) + if (*p > 0.f) *p=(0.12f-*p)*100.f; + + //introduce random bumbs? + if (randomCounter==0) + { + *p *= 2.f; + randomCounter=30; //something random here... + } + --randomCounter; + } + else *p=0.f; //else clear the pixel value + ++p; + } + + dt.SaveImage("dt2.ics"); + i3d::GaussIIR(dt,3.f,3.f,9.f); + p=dt.GetFirstVoxelAddr(); + for (; p != pL; ++p) + { + //uncertainty in the number of incoming photons + float noiseMean = sqrt(*p), // from statistics: shot noise = sqrt(signal) + noiseVar = noiseMean; // for Poisson distribution E(X) = D(X) + + *p=ceilf(*p + GetRandomPoisson(noiseMean) - noiseVar); + + //constants are parameters of Andor iXon camera provided from vendor: + //photon shot noise: dark current + *p+=GetRandomPoisson(0.06f); + + //read-out noise: + // variance up to 25.f (old camera on ILBIT) + // variance about 1.f (for the new camera on ILBIT) + //*p+=GetRandomGauss(700.f,90.f); + } + + i3d::FloatToGrayNoWeight(dt,texture); +} + + void ActiveMesh::CenterMesh(const Vector3F& newCentre) { //calc geom. centre diff --git a/cmath3d/TriangleMesh.h b/cmath3d/TriangleMesh.h index 25b2d1b83dfbb2f252116041dfe0e06a94823d27..f332791295842ec0f97211386abe2628f9ca88f7 100644 --- a/cmath3d/TriangleMesh.h +++ b/cmath3d/TriangleMesh.h @@ -109,6 +109,8 @@ class ActiveMesh void RenderMask(i3d::Image3d<i3d::GRAY16>& mask,const bool showTriangles=false); void RenderMaskB(i3d::Image3d<i3d::GRAY16>& mask); + void RenderOneTimeTexture(const i3d::Image3d<i3d::GRAY16>& mask, + i3d::Image3d<i3d::GRAY16>& texture); void CenterMesh(const Vector3F& newCentre); void ScaleMesh(const Vector3F& scale); diff --git a/src/main.cpp b/src/main.cpp index e05ae4ef1ce7a2b135a994bea601a2f2bed3b03b..32ff18dfe3494e595e20295e40bbff290bfbe6b1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -45,7 +45,8 @@ int main(void) //ISBI: - for (int i=0; i < 50; ++i) + //for (int i=0; i < 50; ++i) + for (int i=10; i < 11; ++i) { LoadNewMesh(i); RenderMesh(i); @@ -159,5 +160,11 @@ int RenderMesh(int fileNo) mask.SaveImage(fileName); std::cout << "rendered mesh #" << fileNo << "\n"; + i3d::Image3d<i3d::GRAY16> texture; + mesh.RenderOneTimeTexture(mask,texture); + sprintf(fileName,"text_t%03d.tif",fileNo); + texture.SaveImage(fileName); + std::cout << "textured mesh #" << fileNo << "\n"; + return(0); }