(*
  Makes a mesh from a black and white image with white parts being the interior.
The distances between nodes are approximately dx pixels.
*)
makemesh[inputFilename_,dx_] := Module[{img,region,mesh},
img=ColorNegate[Import[inputFilename]];
region = ImageMesh[img];
mesh=TriangulateMesh[region,MaxCellMeasure->dx]
]
getnormal[{i1_,i2_},coords_]:={{0,-1},{1,0}}.(coords[[i1]]-coords[[i2]])
getnormal[{i1_,i2_,i3_},coords_]:=Cross[coords[[i1]]-coords[[i2]],coords[[i2]]-coords[[i3]]]
(*
Save mesh in HDF5 format. 
Normals in nodes are computed as averages of normals of dim-1 "faces" bordering on the node.
Nodes closer than cutoffDx are cleared.
The mesh is saved in HDF5 format (you filename must end in .h5) ready for import into C++
*)
savemesh[mesh_, cuttofDx_,outputFilename_,outputGroupname_:"domain"]:=Module[{coords,dim,int,bnd,bndmap,types,normals,bndlines,rot,problematic,problematicBnd,remove,labeled,smallerThan,cnt,pidx,me},
coords=MeshCoordinates[mesh];
dim = Length[coords[[1]]];
Print["Generated "<>ToString[Length[coords]]<>" points."];
int=Map[First,MeshCells[mesh,{0,"Interior"}]];
bnd=Map[First,MeshCells[mesh,{0,"Boundary"}]];
bndmap = Table[-1,Length[coords]];
types=bndmap;
types[[int]]=1;
types[[bnd]]=-1;
Table[bndmap[[bnd[[i]]]]=i-1;,{i,1,Length[bnd]}];
normals=Table[Table[0,{i,dim}], {x,bnd}];
bndlines = MeshCells[mesh,{dim-1,"Boundary"}];
Table[
With[{n=getnormal[line,coords]},
Table[
normals[[bndmap[[idx]]+1]]+=n;
,{idx,line}]];
,{line,First/@bndlines}];
normals=Normalize/@normals;

(* clean too close points *)
labeled=Table[coords[[i]]->i, {i, Length[coords]}];
problematic = Nearest[labeled,coords,{1000, cuttofDx}];
remove =Union@Flatten@Rest@Reap[Table[me=p[[1]];Table[If[me < o, Sow[o],],{o, Rest[p]}],{p,problematic}];];
problematicBnd =Sort[Select[bndmap[[Flatten@remove]]+1,#≠0&]];
normals=Delete[normals,List/@problematicBnd];
smallerThan=Table[0,{i,Length[bndmap]}];
pidx=1; cnt=0;
Table[
If[pidx<=Length@problematicBnd&&problematicBnd[[pidx]]==i-1,cnt++;pidx++;,];
smallerThan[[i]]=cnt;
,{i,Length@bndmap}];
bndmap-=smallerThan;
bndmap[[int]]=-1;
bndmap=Delete[bndmap,List/@remove];
coords=Delete[coords,List/@remove];
types=Delete[types,List/@remove];
int = Flatten@Position[types,1];
bnd = Flatten@Position[types,-1];
Print["Exporing "<>ToString[Length[coords]]<>" points."];

Export[outputFilename,{bndmap,Transpose@coords,types,Transpose@normals},{"Datasets",{outputGroupname<>"/bmap",outputGroupname<>"/pos",outputGroupname<>"/types",outputGroupname<>"/normals"}}];
Print["Exported."];
{coords,int,bnd,bndmap,normals} 
]

mesh=makemesh["medusa.png",30]

{coords,int,bnd,bndmap,normals}=savemesh[mesh, 3, "medusa.h5"];

{Length[int]+Length[bnd],Length[coords],Length[bndmap],Length[normals],Length[bnd]}

arrows = Table[Arrow[{coords[[i]], coords[[i]]+30 normals[[bndmap[[i]]+1]]}],{i,bnd}];
Show[
ListPlot[{coords[[int]],coords[[bnd]]},AspectRatio->Automatic],
Graphics[{Arrowheads[Medium],arrows}]
]

model=Import["heatsink.stl","BoundaryMeshRegion"]
mesh = TriangulateMesh[model,MaxCellMeasure->1];
MeshCellCount[mesh,0]

{coords,int,bnd,bndmap,normals}=savemesh[mesh, 0.1, "heatsink.h5"];

arrows = Table[Arrow[{coords[[i]], coords[[i]]+2 normals[[bndmap[[i]]+1]]}],{i,bnd}];
Show[
ListPointPlot3D[{coords[[int]],coords[[bnd]]},AspectRatio->Automatic],
Graphics3D[{Arrowheads[Small],arrows}]
]