<<DiscreteMath`ComputationalGeometry`
A= ReadList["3d.dat", {Real,Real,Real}];
(* B=Partition[Flatten[Table[{i+Random[],j+Random[],Sin[i*j/30]+Random[]},{i,25},{j,25}]],3] *)

R=TriangularSurfacePlot[A,Background->GrayLevel[0.999],
  DisplayFunction->Identity,Boxed->False,Axes->False,AspectRatio->1]
P=Table[R[[1,k,1]],{k,Length[R[[1]]]}]

WriteString["surface.inc","union{\n"];
Do[WriteString["surface.inc","T(<",
                                 P[[i,1,1]],",",P[[i,1,2]],",",P[[i,1,3]],">,<",
                                 P[[i,2,1]],",",P[[i,2,2]],",",P[[i,2,3]],">,<",
                                 P[[i,3,1]],",",P[[i,3,2]],",",P[[i,3,3]],">)\n"],{i,Length[P]}];
WriteString["surface.inc","}\n"];





