\(おまじない\)

<<Utilities`MemoryConserve`
$MemoryIncrement=1000000;
<< LinearAlgebra`MatrixManipulation`;
Off[General::spell1,MemoryConserve::start,MemoryConserve::end];

■  Normalization

\(\(imageNormalize = 
        Compile[{{data2D, _Real, 
              2}}, 
		Module[{minimum = 
                Min[data2D]}, 
			\((data2D - minimum)\)/
              Max[data2D - minimum]]];\)\)

■ Base Functions

daub2={{1.,1.},{1.,-1.}}/Sqrt[2.];
daub4={{1.+Sqrt[3.], 3.+Sqrt[3.],3.-Sqrt[3.], 1.-Sqrt[3.]},
	{1.-Sqrt[3.],-3.+Sqrt[3.],3.+Sqrt[3.],-1.-Sqrt[3.]}}/
	(4 Sqrt[2.]);
daub6={{ .332670552950, .806891509311, .459877502118,
		-.135011020010, -.085441273882, .035226291882},
	   { .035226291882, .085441273882,-.135011020010,
	    -.459877502118,.806891509311,-.332670552950}};
daub8={{ 0.325803428051, 1.010945715092, 0.892200138246,
	    -0.039575026236,-0.264507167369, 0.043616300475,
	     0.046503601071,-0.014986989330},
	   {-0.014986989330,-0.046503601071, 0.043616300475,
		 0.264507167369,-0.039575026236,-0.892200138246,
		 1.010945715092,-0.325803428051}}/Sqrt[2.];
daub10={{ .160102397974, .603829269797, .724308528438,
		  .138428145901,-.242294887066,-.032244869585,
		  .077571493840,-.006241490213,-.012580751999,
		  .003335725285},
	    { .003335725285, .012580751999,-.006241490213,
	     -.077571493840,-.032244869585, .242294887066,
	      .138428145901,-.724308528438, .603829269797,
	     -.160102397974}};
daub12={{ 0.111540743350, 0.494623890398, 0.751133908021,
		  0.315250351709,-0.226264693965,-0.129766867567,
		  0.097501605587, 0.027522865530,-0.031582039318,
		  0.000553842201, 0.004777257511,-0.001077301085},
		{-0.001077301085,-0.004777257511, 0.000553842201,
		  0.031582039318, 0.027522865530,-0.097501605587,
		 -0.129766867567, 0.226264693965, 0.315250351709,
		 -0.751133908021, 0.494623890398,-0.111540743350}};
daub14={{ .077852054085, .396539319482, .729132090846, 
 		  .469782287405,-.143906003929,-.224036184994,
 		  .071309219267, .080612609151,-.038029936935, 
		 -.016574541631, .012550998556, .000429577973, 
		 -.001801640704, .000353713800},
		{ .000353713800, .001801640704, .000429577973, 
		 -.012550998556,-.016574541631, .038029936935,
		  .080612609151,-.071309219267,-.224036184994,
		  .143906003929, .469782287405,-.729132090846,
		  .396539319482,-.077852054085}};
daub16={{ .054415842243, .312871590914, .675630736297, 
 		  .585354683654,-.015829105256,-.284015542962,
 		  .000472484574, .128747426620,-.017369301002, 
		 -.044088253931, .013981027917, .008746094047, 
		 -.004870352993,-.000391740373, .000675449406,
		 -.000117476784},
		{-.000117476784,-.000675449406,-.000391740373,
		  .004870352993, .008746094047,-.013981027917,
		 -.044088253931, .017369301002, .128747426620,
		 -.000472484574,-.284015542962, .015829105256,
		  .585354683654,-.675630736297, .312871590914,
		 -.054415842243}};
daub18={{ .038077947364, .243834674613, .604823123690, 
 		  .657288078051, .133197385825,-.293273783279, 
		 -.096840783223, .148540749338, .030725681479,
		 -.067632829061, .000250947115, .022361662124,
		 -.004723204758,-.004281503682, .001847646883,
          .000230385764,-.000251963189, .000039347320},
        { .000039347320, .000251963189, .000230385764,
         -.001847646883,-.004281503682, .004723204758,
          .022361662124,-.000250947115,-.067632829061,
         -.030725681479, .148540749338, .096840783223,
         -.293273783279,-.133197385825, .657288078051, 
         -.604823123690, .243834674613,-.038077947364}};
daub20={{ 0.026670057901, 0.188176800078, 0.527201188932,
		  0.688459039454, 0.281172343661,-0.249846424327,
		 -0.195946274377, 0.127369340336, 0.093057364604,
		 -0.071394147166,-0.029457536822, 0.033212674059,
		  0.003606553567,-0.010733175483, 0.001395351747,
		  0.001992405295,-0.000685856695,-0.000116466855,
		  0.000093588670,-0.000013264203},
		{-0.000013264203,-0.000093588670,-0.000116466855,
		  0.000685856695, 0.001992405295,-0.001395351747,
		 -0.010733175483,-0.003606553567, 0.033212674059,
		  0.029457536822,-0.071394147166,-0.093057364604,
		  0.127369340336, 0.195946274377,-0.249846424327,
		 -0.281172343661, 0.688459039454,-0.527201188932,
		  0.188176800078,-0.026670057901}};

■  Construction of Wavelet Functions

c[n_,baseType_]:=Module[{mat,i,result},
	result=mat=AppendRows[baseType,
		ZeroMatrix[2,n-Length[baseType[[1]]]]];
       	rotateRightofList[a_List]:=RotateRight/@ a;
        rotateRight[b_List,k_]:=
        Nest[rotateRightofList,b,k];
        Do[result=
        	AppendColumns[result,rotateRight[mat,2 i]],
            {i,n/2-1}];
	result]

cPrime[i_,n_,baseType_]:=BlockMatrix[{{c[n/2^i,baseType],
	ZeroMatrix[n/2^i,n-n/2^i]},{ZeroMatrix[n-n/2^i,n/2^i],
	IdentityMatrix[n-n/2^i]}}]

p[n_]:=Module[{mat,i},mat=IdentityMatrix[n];
	Do[mat[[i]]=RotateRight[mat[[i]],i-1],{i,2,n}];
    Do[mat[[i]]=RotateRight[mat[[i]]],{i,n/2+1,n}];
    mat];
pPrime[i_,n_,baseType_]:=BlockMatrix[{{p[n/2^i],
	ZeroMatrix[n/2^i,n-n/2^i]},{ZeroMatrix[n-n/2^i,n/2^i],
  	IdentityMatrix[n-n/2^i]}}];

waveletMatrix[n_,baseType_]:=
	Apply[Dot,Flatten[Table[{pPrime[j,n,baseType],
		cPrime[j,n,baseType]},{j,Log[2,n]-
		Ceiling[Log[2.,Length[baseType[[2]]]]],1,-1}],
		1]].p[n].c[n,baseType];

■   Outlines and Sketch

● 1st sample

\(\(sampleDataA = {{1, 1, 1, 1, 1, 1, 1, 1}, 
		{1, 1, 1, 0, 0, 0, 1, 
            1}, 
		{1, 1, 1, 0, 0, 0, 1, 1}, 
		{1, 1, 1, 0, 0, 0, 1, 
            1}, 
		{1, 1, 1, 0, 0, 0, 1, 1}, 
		{1, 1, 1, 0, 0, 0, 0, 
            0}, 
		{1, 1, 1, 1, 1, 0, 0, 0}, 
		{1, 1, 1, 1, 1, 1, 1, 
            1}};\)\)

\(\(sampleFigureA = 
        ListDensityPlot[sampleDataA, 
	
          PlotLabel -> "Sample figure A", 
	
          Mesh -> False];\)\)

weightenWaveletSegnetation[test_,m_]:=
	Module[{n,n2,m2,m3,m4,ms,me,mx,k,p,ps,pe,wMatrix,wTMatrix,
		partSpect,segment,out,gMax,maxAbs,minAbs,list={}},
			n=Length[test];n2=n/2;
			mx=m;Which[m<=0,mx=1,m>n,mx=n];m2=n/mx;
			wMatrix=waveletMatrix[n,daub2];
			wTMatrix=Transpose[wMatrix];
			out=ZeroMatrix[n];
			segment=wMatrix.test.wTMatrix;
			Do[segment[[i,j]]=0,{i,n2},{j,n2}];
			segment=wTMatrix.segment.wMatrix;
			gMax=Max[Abs[segment]];
			minAbs=N[gMax/(mx mx)];
		Do[
			m3=k*m2;m4=p*m2;
				ms=m3+1;me=ms+m2-1;
				ps=m4+1;pe=ps+m2-1;
				partSpect=ZeroMatrix[n];
			Do[partSpect[[i,j]]=segment[[i,j]],{i,ms,me},{j,ps,pe}];	
			maxAbs=Max[Abs[partSpect]];
			If[maxAbs>=minAbs,
				partSpect=Round[partSpect/maxAbs],
				partSpect=0. partSpect];
			Do[out[[i,j]]=(maxAbs/gMax) Abs[partSpect[[i,j]]];
				If[partSpect[[i,j]]>0,
      				list=Append[list,{RGBColor[1,1,0],
      					Point[{j,i}]}]],
      		{i,ms,me},{j,ps,pe}],
      		
      	{p,0,mx-1},{k,0,mx-1}];
		{out,Map[Graphics,list]}
	];

Timing[segmentA16=weightenWaveletSegnetation[sampleDataA,16];]

sketchA16=ListDensityPlot[1-Abs[segmentA16[[1]]],
	PlotLabel -> "Sketch of sample A",
	Mesh->False,DisplayFunction->Identity];
outlineA16=Show[segmentA16[[2]],
	PlotLabel -> "Outline of sample A",
	Background->GrayLevel[0],AspectRatio->1,Frame->True,
	DisplayFunction->Identity];
Show[GraphicsArray[{sampleFigureA,outlineA16,sketchA16}]];

● 2nd sample

\(\(sampleDataA = {{1, 1, 1, 1, 1, 1, 1, 1}, 
		{1, 1, 1, 0, 0, 1, 1, 
            0}, 
		{1, 1, 1, 1, 0, 1, 1, 0}, 
		{1, 1, 1, 1, 0, 0, 1, 
            1}, 
		{1, 1, 1, 1, 1, 1, 1, 1}, 
		{1, 1, 1, 0, 0, 1, 1, 
            0}, 
		{1, 1, 1, 1, 1, 1, 1, 0}, 
		{1, 1, 1, 1, 1, 1, 1, 
            1}};\)\)

sampleFigureA=ListDensityPlot[sampleDataA,
	PlotLabel->"Sample figure A",
	Mesh->False];
sampleFigureB=ListDensityPlot[sampleDataB,
	PlotLabel->"Sample figure B",
	DisplayFunction->Identity,Mesh->False];

Timing[segmentA16=weightenWaveletSegnetation[sampleDataA,16];]

sketchA16=ListDensityPlot[1-Abs[segmentA16[[1]]],
	PlotLabel -> "Sketch of sample A",
	Mesh->False,DisplayFunction->Identity];
outlineA16=Show[segmentA16[[2]],
	PlotLabel -> "Outline of sample A",
	Background->GrayLevel[0],AspectRatio->1,Frame->True,
	DisplayFunction->Identity];
Show[GraphicsArray[{sampleFigureA,outlineA16,sketchA16}]];

●  Photo sample

sampleFigure = 
	ListDensityPlot[sample,
		Mesh -> False,PlotLabel->"Sample data",
		Background->GrayLevel[1]];

Timing[segment1=weightenWaveletSegnetation[sample,1];]

Show[segment1[[2]],
	PlotLabel -> "1 Segmented outline",
	Frame->True,AspectRatio->1,Background->GrayLevel[0]];

ListDensityPlot[1-segment1[[1]],
	PlotLabel->"1 segment sketch",Mesh->False];

Timing[segment16=weightenWaveletSegnetation[sample,16];]

Show[segment16[[2]], 
	PlotLabel -> "16 by 16 Segments",
	Background->GrayLevel[0],AspectRatio->1,Frame->True];

ListDensityPlot[1-segment16[[1]],
	PlotLabel->"16 by 16 segmented sketch",Mesh->False];