2.1.

Criterion for Voxel Size Selection

2.1.1.

Criterion based on geometrical approximation

In the VMC method, the curved interfaces are approximated by the boundary of the voxel. As shown in Fig. 1, the cross-section of a cylindrical blood vessel (vessel diameter $d=120\mu \mathrm{m}$) is approximated by regular voxels (red blocks in Fig. 1). If the voxel size is large, $d/6$ for example [Fig. 1(a)], then the approximate interface significantly deviates from the curve, which results in a deviation of the vessel area in the discretized domain. Consequently, a large error exists when the photons refract or reflect on the vessel wall. As the voxel size is reduced to $d/12$ [Fig. 1(b)] and $d/24$ [Fig. 1(c)], the serrated polygonal boundary is smoothed, and the cross-section area approaches the real vessel.

Fig. 1

Representation of blood vessel cross-section with a diameter of $120\mu \mathrm{m}$ by voxels when the voxel size is (a) $d/6$, (b) $d/12$, and (c) $d/24$.

For consistency, a dimensionless parameter is defined as $\lambda =d^{*}/\mathrm{\Delta }{w}^{*}$, where $d^{*}=d·{\mu }_{\mathrm{t}}$ is the dimensionless length and $\mathrm{\Delta }{w}^{*}=\mathrm{\Delta }w{\mu }_{\mathrm{t}}$ is the dimensionless voxel width. For a circular blood vessel cross-section, the ratio of the voxel-generated blood vessel area (red part in Fig. 1) to the real blood vessel area (area within dashed black lines in Fig. 1) is then related to $\lambda$ as indicated in Fig. 2. As can be seen, the area ratio quickly approaches 1 with increasing $\lambda$. If $\lambda$ is $>10$, then the area ratio changes very little ($<3\%$), and the relationship is true for all 3 vessel diameters (60, 120, and $180\mu \mathrm{m}$). Therefore, $\lambda =10$ can be used as a criterion from the geometrical point of view by which the dimensionless voxel width is given.

Eq. (1)

$$\mathrm{\Delta }{w}^{*}=0.1d^{*}.$$

Fig. 2

Relationship between $\lambda$ and the ratio of the generated blood vessel area to the real blood vessel area.

This criterion applies beyond the circular interface. For the irregular curved interface, the dimensionless diameter $d^{*}$ in Eq. (1) can be determined as $d^{*}=2{\mu }_{\mathrm{t}}/K_{\max }$, where $K_{\max }$ is the maximum curvature of the curved interface.

2.1.2.

Criterion based on geometrical approximation

A laser beam that propagates through a homogeneous medium with uniform absorption and scattering properties, as shown in Fig. 3, is considered. For the light with fluence of $F_{\mathrm{o}}$ at the incident point o, the analytical solution of light fluence $F(w)$ at the depth of $w$ within the medium can be assumed to be attenuated in the skin tissue as the consecutive black solid line in Fig. 3 and presented according to the Beer-Lambert law as follows:

Eq. (2)

$$F(w)=F_{\mathrm{o}}·\exp (-{\mu }_{\mathrm{t}}w),$$

where ${\mu }_{\mathrm{t}}$ is the attenuation coefficient of the medium. In a high absorbing medium, ${\mu }_{\mathrm{t}}={\mu }_{\mathrm{a}}$ (e.g., blood). In a high scattering medium defined by the diffusion theory, ${\mu }_{\mathrm{t}}={\{3{\mu }_{\mathrm{a}}\times [{\mu }_{\mathrm{a}}+(1-\mathrm{g}){\mu }_{\mathrm{s}}]\}}^{1/2}$ (like the situation in epidermis and dermis). ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}$, and $g$ are the absorption coefficient, scattering coefficient, and anisotropic factor, respectively.⁴

Fig. 3

Schematic of light attenuation and grid separation in homogeneous medium.

To calculate light attenuation with the MC method, a grid system has to be constructed to score the energy deposition in each grid cell (voxel) discretely. In this case, the medium can be discretized into cells of width $\mathrm{\Delta }w$ (assuming a uniform grid). The coordinates at the grid interface $w{s}_i$ ($i=1$ to $k$, $k$ is the node index of the last element) and the cell center $w_i$ ($i=1$ to $k$) can then be easily determined by

Eq. (3)

$$\mathrm{\Delta }w=w{s}_{i+1}-w{s}_i=w_{i+1}-w_i.$$

During the calculation, the photon weight absorbed within the cell (control volume) would be stored in each cell. After sufficient photon traveling, the calculated light fluence would be presented discretely as an array $\overline{F}(w_i)$ (zigzag red line shown in Fig. 3). $\bar{F}(w_i)$ can be calculated based on the energy (area under both lines) conservation with analytical solution over the domain of each cell $i$.

Eq. (4)

$$\overline{F}(w_i)\mathrm{\Delta }w={\int }_{w_i-\mathrm{\Delta }w/2}^{w_i+\mathrm{\Delta }w/2}F(w)\mathrm{d}w.$$

Integrating Eq. (2) into Eq. (4), we obtain

Eq. (5)

$$\overline{F}(w_i)=\frac{F_{\mathrm{o}}\times {\int }_{w_i-\mathrm{\Delta }w/2}^{w_i+\mathrm{\Delta }w/2}\exp (-{\mu }_{\mathrm{t}}w)\mathrm{d}w}{\mathrm{\Delta }w}=\frac{F_{\mathrm{o}}\times \exp (-{\mu }_{\mathrm{t}}{w}_i)[\exp ({\mu }_{\mathrm{t}}\mathrm{\Delta }w/2)-\exp (-{\mu }_{\mathrm{t}}\mathrm{\Delta }w/2)]}{\mathrm{\Delta }w·{\mu }_{\mathrm{t}}}.$$

However, if the light attenuation is large within the grid domain, then the scored light fluence may deviate from the analytical solution, particularly at the cell interface, although the energy conservation in each grid domain can be achieved. For example, the light attenuation is obvious in the first cell (Fig. 3). As a result, the scored light fluence in the first cell $\bar{F}(w_1)$ is much lower than the analytical solution $F(w)$ at the first cell interface $w{s}_1$, as shown in Fig. 3. To diminish such difference and increase the accuracy and sensitivity of the MC result, the grid cell width should ensure that the scored light fluence $\bar{F}(w_i)$ in the $i$’th cell is not significantly lower than the analytical solution of light fluence at the cell interface $w{s}_{i-1}F(w{s}_{i-1})$. An approximation factor $\delta$ ($0<\delta <1$) is then introduced and defined as

Eq. (6)

$$\delta =\frac{\overline{F}(w_i)}{F(w{s}_i)}=\frac{F\times \exp (-{\mu }_{\mathrm{t}}{w}_i)[\exp ({\mu }_{\mathrm{t}}\mathrm{\Delta }w/2)-\exp (-{\mu }_{\mathrm{t}}\mathrm{\Delta }w/2)]}{\mathrm{\Delta }w·{\mu }_{\mathrm{t}}}/\{F·\exp [-{\mu }_{\mathrm{t}}(w_i-\frac{w}{2})\left]\right\}=\frac{1-\exp (-\mathrm{\Delta }w{\mu }_{\mathrm{t}})}{\mathrm{\Delta }w{\mu }_{\mathrm{t}}}.$$

A dimensionless cell width is defined, $\mathrm{\Delta }{w}^{*}=\mathrm{\Delta }w{\mu }_{\mathrm{t}}$, which is the ratio of the cell width $\mathrm{\Delta }w$ to the light penetration depth $1/{\mu }_{\mathrm{t}}$ in the medium. Equation (6) then becomes

Eq. (7)

$$\delta =\frac{1-\exp (-\mathrm{\Delta }{w}^{*})}{\mathrm{\Delta }{w}^{*}}.$$

$\delta$ implies the approximation of the simulated result by the MC method to the analytical solution. A high value of $\delta$ must be achieved to ensure the computational accuracy of MC. For instance, $\delta \ge 0.9$ means that the computational error by MC simulation compared with that of the analytical solution is $<10\%$. As can be seen in Eq. (6), as the criterion, the corresponding $\mathrm{\Delta }{w}^{*}$ approaches 0.2, as

Eq. (8)

$$\mathrm{\Delta }{w}^{*}=0.2.$$

This criterion is introduced based on the light energy scoring independent of the photon movement. Therefore, this criterion is applicable to any MC method.

In summary, 2 factors affect the selection of voxel size in the VMC method: the discretization of irregular geometry [Eq. (1)] and the discretization of light propagation [Eq. (8)]. These 2 equations are plotted in Fig. 4. As can be seen, if the dimensionless diameter $d^{*}$ is $<2$, then the voxel size is dominated by the geometric approximation of the curved interface [Eq. (1)]. Otherwise, the voxel size is dominated by the optical property of the medium [Eq. (8)]. Finally, these 2 factors can be considered in the following unified expression:

Eq. (9)

$$\{\begin{array}{ll}\mathrm{\Delta }{w}^{*}=0.1d^{*}& d^{*}\le 2\\ \mathrm{\Delta }{w}^{*}=0.2& d^{*}\ge 2\end{array}.$$

Fig. 4

Relationship of $\mathrm{\Delta }{w}^{*}$ and $d^{*}$.

Notably, the attenuation coefficient of the medium is incorporated into the nondimensional voxel width, $\mathrm{\Delta }{w}^{*}=\mathrm{\Delta }w/(1/{\mu }_{\mathrm{t}})$. The attenuation coefficient significantly affects the real size of the voxel. In a high absorbing medium, ${\mu }_{\mathrm{t}}={\mu }_{\mathrm{a}}$ (e.g., blood). As a result, the voxel size is independent with scattering properties. On the other hand, in a high scattering medium, ${\mu }_{\mathrm{t}}={\{3{\mu }_{\mathrm{a}}\times [{\mu }_{\mathrm{a}}+(1-\mathrm{g}){\mu }_{\mathrm{s}}]\}}^{1/2}$ (like the situation in epidermis and dermis). The voxel size is highly related to scattering properties. The local voxel size should be small with high scattering coefficient and large with low scattering coefficient. As for the medium with nonuniform scattering properties, the nonuniform grid should be used.

2.2.

Criterion for the Determination of Photon Number

In MC simulation, the number of photons should be sufficiently large to satisfy the statistical requirement at the given voxel.⁸ The number of photons that arrived at each voxel depends on the total photon number in the simulation and the voxel size. For a coarse grid with a large cell size, a small number of photons may produce sufficient photons for each voxel. However, a small voxel should be used for a medium with large attenuation coefficients (either absorption or scattering), as shown in Eq. (9). The number of photons that arrive at each voxel will then become smaller if the same number of total photons is used, which may not meet the basic requirement of statistics and will lead to large simulation errors.

To demonstrate the effect of the photon number on the simulation result, a comparison of the light fluence distribution by the MC method with that by the analytical solution (Beer-Lambert law) was performed for a 1-D photon propagation problem within a homogeneous medium (blood) with uniform optical properties. The wavelength of the laser is 585 nm. The optical properties of the blood are listed in Table 2.⁸ As illustrated in Fig. 5, the MC method agrees well with the analytical solution when the total photon number is $>{10}^7$. A large error can be observed for the case with a photon number of ${10}^5$. The relative error of the MC method from the analytical solution can be plotted as a function of the photon number per voxel in Fig. 6. Interestingly, the error has a close relationship with the number of photon per voxel. Two regions show different dependences between the 2 variables. When the number of photons per voxel is $<5$, the MC method produces greater error. However, when the number is $>5$, the error produced by the MC method decreases. Further increasing the photon number per voxel will not significantly improve the accuracy of the MC method.

Table 2

Optical properties of skin components.8

Fig. 5

Comparison of analytical solution and Monte Carlo (MC) simulation result with various photon numbers.

Fig. 6

Relation of computation error by MC with photon number per voxel.

Apparently, 5 photons per voxel can control the selection of the total photons in the MC simulation. The 5 photons per voxel can be used as the basic requirement for the minimum photons required by the MC simulation, that is, the total photons required by the MC simulation will be 5 times that of the voxels in the simulation.

2.3.

Criterion Validation

The in-house VMC code is written by using MATLAB®, and the validation of the code can be seen in a previous study.¹³ The light propagation during laser treatment of PWS is taken as an example to validate the above criteria of voxel size and photon number.

PWS birthmarks are congenital vascular malformations that occur in $\sim 0.3\%$ of newborns. PWS is composed of ectatic venular capillary blood vessels, with diameters that range from 10 to $300\mu \mathrm{m}$, buried within healthy dermal tissue. PWS can be treated by pulsed dye laser with wavelength in the visible band (e.g., 585 nm) based on the theory of selective photothermolysis, with hemoglobin serving as the target chromophore.²³ According to the theory, the PWS blood vessels can be thermally damaged selectively because of their preferential absorption of laser energy compared with normal skin tissues, which are minimally affected.

The skin that contains PWS can be assumed as a 2-layer geometry composed of epidermis and dermis layers, as shown in Fig. 7. A single PWS blood vessel that comprises highly absorbing blood is considered, buried in the center of the dermis and parallel to the skin surface. The dimensions of the computational skin model are $x_k\times {\mathrm{y}}_k\times {\mathrm{z}}_k=1440\mu \mathrm{m}\times 1440\mu \mathrm{m}\times 1020\mu \mathrm{m}$. The thicknesses of the epidermis and dermis are 60 and $960\mu \mathrm{m}$, respectively. The diameter of the blood vessel is $120\mu \mathrm{m}$, and the depth of the vessel center to the skin surface is $250\mu \mathrm{m}$. A typical PWS treatment wavelength of 585 nm is used, as shown in Table 2.⁸ The circular laser beam is fired in the center of the skin surface, with a diameter of $1000\mu \mathrm{m}$.

Fig. 7

Schematic of the 2-layer skin model containing port wine stain.

For such skin tissue with PWS, the attention should be on the blood vessel. Based on the criterion, for the blood vessel with a diameter of $120\mu \mathrm{m}$, $d^{*}=120\mu \mathrm{m}\times 191/({10}^4\mu \mathrm{m})=2.292$. As a result, the voxel size for the blood vessel $\mathrm{\Delta }w=0.2/(191\times {10}^{-4})=10.47\mu \mathrm{m}$. To test for reasonability, the effect of the voxel size on the light attenuation within a blood vessel along the centerline of the laser spot is examined. The voxel size is selected as 5, 10, and $20\mu \mathrm{m}$, the cell numbers are $288\times 288\times 204$, $144\times 144\times 102$, and $72\times 72\times 51$, respectively, and the corresponding photon numbers are $5\times {10}^6$, $4\times {10}^7$, and $6.25\times {10}^5$, all of which meet the requirement for the total number of photons. As can be seen from Fig. 8, the highest energy density predicted with the voxel size of $20\mu \mathrm{m}$ ($136\mathrm{J}/{\mathrm{cm}}^3$) is significantly lower than those with the voxel size of either $10\mu \mathrm{m}$ ($167\mathrm{J}/{\mathrm{cm}}^3$) or $5\mu \mathrm{m}$ ($173\mathrm{J}/{\mathrm{cm}}^3$). At the same time, the calculation result with a cell size of $10\mu \mathrm{m}$ is $167\mathrm{J}/{\mathrm{cm}}^3$, which is close to the result with a cell size of $5\mu \mathrm{m}$ ($173\mathrm{J}/{\mathrm{cm}}^3$).

Fig. 8

Light attenuation within a blood vessel along the centerline of the laser spot with voxel size of 5, 10, and $20\mu \mathrm{m}$.

The energy deposition with cell size of $1\mu \mathrm{m}$ is the basic solution, and the simulated result is compared with that of a larger cell size from 2 to $80\mu \mathrm{m}$. The computation error of the maximum energy deposition within the blood vessel is illustrated in Fig. 9 as a function of the cell size. The computation error becomes smaller as the cell size decreases. The computation error is $>10\%$ when the voxel size is $>10\mu \mathrm{m}$. Therefore, the largest cell size for the $120\text{-}\mu \mathrm{m}$ blood vessel with a wavelength of 585 nm is $10\mu \mathrm{m}$, which is consistent with the above criterion.

Fig. 9

Computation error of the maximum light energy deposition within blood vessel as a function of grid size from 2 to $80\mu \mathrm{m}$.

The effect of the total photon number on the light density distribution in the blood vessel is shown in Fig. 10, with the voxel size of $20\mu \mathrm{m}$ [Fig. 10(a)] and $5\mu \mathrm{m}$ [Fig. 10(b)]. The minimum number of photons based on the requirement of 5 per voxel is $6.25\times {10}^5$ and $4\times {10}^7$ for 20 and $5\mu \mathrm{m}$, respectively. The results from different photon numbers (${10}^5$ to ${10}^7$ for $20\mu \mathrm{m}$ voxels and ${10}^6$ to ${10}^8$ for $5\mu \mathrm{m}$ voxels) are also plotted in the same figures. By comparison, if the total number of photons is less than the corresponding minimum number, then the calculated results show significant deviation from those results with more photons. For example, in the case of $20\mu \mathrm{m}$ voxels [Fig. 10(a)], 2 results with a photon number larger than the minimum number ($6.25\times {10}^5$) agree fairly well with that which uses the minimum photon number ($6.25\times {10}^5$), whereas the result for photon number of ${10}^5$ significantly deviates from the other results. Similar results can be observed in the case of $5\mu \mathrm{m}$ voxels [Fig. 10(b)]. All the results with photons less than the minimum required number ($4\times {10}^7$) show large fluctuations and significantly deviate from the results with large photons. A further increase in the photons beyond the minimum requirement makes no significant improvement. This comparison confirms the conclusion from the 1-D analysis, and the requirement for the number of photons for MC simulation, which is 5 photons per voxel, seems reasonable.

Fig. 10

Energy deposition within blood vessel along tissue depth at the laser spot center with voxel size of (a) $20\mu \mathrm{m}$ and (b) $5\mu \mathrm{m}$ with various photon numbers.

3.1.

Local Grid Refinement and Photon Splitting Technique

Based on the criterion in Eq. (9), the voxel size used in a 3-D VMC is greatly affected by the optical properties of the medium. For a system that consists of various tissue components with vastly different optical properties, the voxel size should be determined by the tissue component with the largest attenuation coefficient (such as blood vessel in PWS) if a uniform grid is used. For such a fine grid, a significant number of photons should also be used to satisfy the criterions for photon selection (5 times the number of voxels in the simulation) to guarantee computational accuracy. However, a coarse grid and fewer photons can be used for tissue components with less absorption, such as dermis.

Although the nonuniform grid is widely used in numerical simulations,^24,25 the uniform grid is usually used in MC simulation, at least for the simulation of laser–tissue interactions.^8,26–28 In this section, a nonuniform grid system is proposed to save computational cost. Aside from this nonuniform grid system, the photon splitting scheme is also developed to increase the photon number needed to satisfy the statistical accuracy requirement of 5 photons per voxel in the dense grid area. A 2-layer skin model that contains a single PWS blood vessel (Fig. 7) is taken as an example to illustrate the nonuniform grid system and photon splitting scheme.

First, the tissue interfaces should be defined to distinguish different grid densities. Two interfaces are in the 3 tissue components, that is, the epidermis–dermis interface and dermis–blood vessel interface. For simplicity, a rectangular cylinder is constructed to enclose the blood vessel as the dermis–blood vessel interface is curved. In the rectangular cylinder region that contains the blood vessel, a dense grid is determined based on the blood optical property by using the criterion in Eq. (9). The other domain of dermis can then be discretized by the coarse grid with a large voxel, with the width determined based on the dermis optical property. By contrast, the epidermis can have its own voxel size through the criterion. The schematic of the 2-D final grid for entire computational domain is shown in Fig. 11(a). The area within the dashed region is the blood vessel, which is represented by the discrete square voxels. For simplicity, the voxel size of the epidermis and dermis is the same.

Aside from the nonuniform grid system, a photon splitting scheme is suggested. For the nonuniform grid system shown in Fig. 11(a), the total number of photons used in the MC simulation can be determined based on the large voxels in the dermis, that is, 5 times the total number of large voxels. When a photon travels from the coarse grid into the fine one, the photon will be split into $n$ subphotons (either inside or outside the dense grid). As shown in Fig. 11(b), $n$ is determined based on the requirement of the photon numbers for statistical accuracy of the fine grid, that is, at least 5 photons per voxel. Then, all $n$ subphotons will be examined before a new photon is projected on the surface again.

Fig. 11

Schematic of (a) nonuniform grid system and (b) photon splitting.

3.2.

Effect of Local Grid Refinement and Photon Splitting Technique

The main objective of the nonuniform grid system and the photon splitting scheme is to reduce computational cost while maintaining accuracy. The single blood vessel skin model (Fig. 7) will be used to demonstrate the improvement of computation efficiency by the nonuniform grid system and the photon splitting scheme.

Based on Eq. (9), the maximum voxel size can be determined for each skin component. For the blood vessel with a diameter of $120\mu \mathrm{m}$, $d^{*}=d{\mu }_{\mathrm{t}}=120\mu \mathrm{m}\times 191/({10}^4\mu \mathrm{m})=2.23$; then the maximum voxel size should be $\mathrm{\Delta }w=0.2/{\mu }_{\mathrm{t}}=10.47\mu \mathrm{m}$. For the epidermis with planar interface, the maximum voxel size $\mathrm{\Delta }w=0.2/{\mu }_{\mathrm{t}}=0.2/79.38=25.19\mu \mathrm{m}$. For the dermis, similar to the epidermis, the maximum voxel size $\mathrm{\Delta }w=0.2/{\mu }_{\mathrm{t}}=0.2/13.9=143.9\mu \mathrm{m}$.

The light transportation in skin tissue and blood vessel was calculated with a uniform grid throughout the entire computational domain, with the voxel size determined by the requirement for the blood vessel because the finest voxels are needed in the blood vessel. The voxel size and the photon number are $10\mu \mathrm{m}$ and ${10}^7$, which has been used in the Ref. 8. The computation was completed in 17.5 h (CPU time) by HP xw6600 Workstation. The average energy density of the entire blood vessel area is $77.31\mathrm{J}/{\mathrm{cm}}^3$.

The result of the nonuniform grid was then checked. Voxel size is fixed for the epidermis ($20\mu \mathrm{m}$) and the rectangular dense grid area that contains the blood vessel ($10\mu \mathrm{m}$), but the dermis has different voxel sizes from 10 to $80\mu \mathrm{m}$ where the absorption is small and a large voxel size can be used, and the total number of photons is maintained at ${10}^7$. The computational time is reduced from 17.5 to 5.5 h when the voxel size in dermis was increased from 10 to $80\mu \mathrm{m}$ by the same workstation, as shown in Fig. 12. The figure also indicated that a further increase in voxel size in dermis would not further reduce the computational time. By contrast, increasing the voxel size in dermis also increases the computational error of the total energy density in the blood vessel, as shown in Fig. 13, although the overall error is within 3%. A coarse grid in the dermis shows no effect on the energy density distribution in the blood vessel, as shown in Fig. 14, which provides the light distributions within the blood vessel for 2 cases with dermis voxel sizes of $20\mu \mathrm{m}$ [Fig. 14(a)] and $80\mu \mathrm{m}$ [Fig. 14(b)], respectively.

Fig. 12

Relation between computation time and voxel size in dermis.

Fig. 13

Error of total energy density in the blood vessel as a function of voxel size in dermis.

Fig. 14

Energy density in refined area with dermis grid size of (a) $20\mu \mathrm{m}$ and (b) $80\mu \mathrm{m}$.

For the fixed voxel size of $20\mu \mathrm{m}$ in the epidermis, $80\mu \mathrm{m}$ voxel in the dermis, and $10\mu \mathrm{m}$ voxel for the rectangular dense grid area that contains the blood vessel, the computational time can be further reduced from 5.5 hours to 3.5 and 2.3 h by reducing the total photon number from ${10}^7$ to ${10}^6$ and ${10}^5$, which correspond to the splitting number $n$ from 1 to 10 and 100, as shown in Fig. 15. Both photon numbers satisfy the photon requirement for fine grid (blood vessel area).

Fig. 15

Relation between splitting number $n$ and computation time.

Figure 16 shows the effect of the splitting number $n$ on the accuracy of the simulation results, which is represented by the percentage of the difference from the average energy density calculated by using the finest uniform grid (voxel size of $10\mu \mathrm{m}$) for the entire computation domain and the total photon number ${10}^7$. The error introduced by photon splitting is relatively small ($<5\%$), which demonstrates that the present photon splitting technique is valid and maintains energy conservation. The light deposition within the grid refinement area with $n=10$ and $n=100$ are plotted in Figs. 17(a) and 17(b). Photon splitting in the fine grid area shows no effect on the energy density distribution in the blood vessel, as shown in Fig. 17.

Fig. 16

Error in the energy density in the blood vessel as a function of splitting number.

Fig. 17

Energy distribution within refined area with splitting number of (a) 10 and (b) 100 (voxel size: epidermis $20\mu \mathrm{m}$, unrefined dermis $80\mu \mathrm{m}$, and refined area $10\mu \mathrm{m}$).