text:
G-0K8J8ZR168
검색
검색 팝업 닫기

Ex) Article Title, Author, Keywords

## Article

Curr. Opt. Photon. 2021; 5(5): 514-523

Published online October 25, 2021 https://doi.org/10.3807/COPP.2021.5.5.514

## Absolute Depth Estimation Based on a Sharpness-assessment Algorithm for a Camera with an Asymmetric Aperture

Beomjun Kim, Daerak Heo, Woonchan Moon, Joonku Hahn

School of Electronic and Electrical Engineering, Kyungpook National University, Daegu 41566, Korea

Corresponding author: jhahn@knu.ac.kr, ORCID 0000-0002-5038-7253

Received: June 14, 2021; Revised: July 19, 2021; Accepted: July 29, 2021

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Methods for absolute depth estimation have received lots of interest, and most algorithms are concerned about how to minimize the difference between an input defocused image and an estimated defocused image. These approaches may increase the complexity of the algorithms to calculate the defocused image from the estimation of the focused image. In this paper, we present a new method to recover depth of scene based on a sharpness-assessment algorithm. The proposed algorithm estimates the depth of scene by calculating the sharpness of deconvolved images with a specific point-spread function (PSF). While most depth estimation studies evaluate depth of the scene only behind a focal plane, the proposed method evaluates a broad depth range both nearer and farther than the focal plane. This is accomplished using an asymmetric aperture, so the PSF at a position nearer than the focal plane is different from that at a position farther than the focal plane. From the image taken with a focal plane of 160 cm, the depth of object over the broad range from 60 to 350 cm is estimated at 10 cm resolution. With an asymmetric aperture, we demonstrate the feasibility of the sharpness-assessment algorithm to recover absolute depth of scene from a single defocused image.

Keywords: Coded aperture, Depth estimation, Image reconstruction

OCIS codes: (100.2000) Digital image processing; (100.3008) Image recognition, algorithms and filters; (100.3020) Image reconstruction-restoration; (120.3940) Metrology; (170.1630) Coded aperture imaging

Depth estimation is one of the most important areas of study in three-dimensional (3D) metrology. Depth estimation occupies a crucial position in a variety of industries, as it has been actively used for self-driving cars, inspection of defects, 3D recognition, for example. Time-of-flight and structured-light illumination are regarded as representative technologies of depth estimation, where a specific light is emitted onto a target and the reflected light signal is collected by a detector. These technologies provide a depth map with high resolution and accuracy, but they usually require expensive and complicated optical instruments.

On the other hand, depth estimation from defocused images has the great advantage that it requires only one or two images from a single conventional camera and its system has a small form factor compared to previous ones. This method originates from the relation between level of defocus (LOD) and depth profile. When the object is placed outside of the depth of field, the captured image is defocused and LOD increases proportionally to the distance from the focal plane. Based on this fact, many studies of depth recovery have concentrated on obtaining LOD from defocused images [13]. Therefore, evaluation of LOD from the scene plays a significant role in various applications such as image deblurring, image segmentation and depth estimation.

In general, methods for depth recovery are classified into relative depth estimation (RDE) and absolute depth estimation (ADE) from the defocused images. In RDE studies, the LOD is evaluated by computing the standard deviation, which is inferred from edge-detection. The standard deviation is related to the size of the point-spread function (PSF) and can be obtained from the high-frequency content using the derivative operator [4, 5]. Pentland [6] proposed the framework of RDE to obtain the focal disparity map by utilizing the relation between high frequency content coming from the difference of two defocused images. But this approach has an inconvenience: It needs to capture several images from the same scene, with different camera parameters. Bae and Durand [7] recovered the depth map from a single defocused image by employing the method of Elder and Zucker [1] and Levin et al. [8]. Zhou and Terence [9] use a ratio of gradient between the defocused input image and the estimated image to calculate the depth map. The above methods evaluate LOD of a scene from the edge-detection of a defocused image modeled by convolution of a gaussian PSF and a sharp image. Thus they take little computational time and uncomplicated processing, but it is difficult to recover absolute depth without the relation between standard deviation and depth layer [10, 11].

Most methods of ADE are based on deconvolution to recover the depth from a defocused image. The features of the PSF are important to estimate the depth precisely, and it is necessary to specially design an optical stop with a coded pattern. Levin et al. [8] proposed the framework of ADE to recover absolute depth of scene and all-focus image using a camera with a centrosymmetric coded aperture. However, the centrosymmetric coded aperture makes it difficult to distinguish if the object is placed in front of or behind the focal plane. Zhou et al. [12] proposed a method using two asymmetric apertures that are optimized for reconstruction of defocused image. They recover the depth map over a broad range and obtain an all-focus image in high quality, but they need to take two defocused images with different coded apertures for the same scene.

In this paper, we propose a new ADE method to recover the depth of scene by using an asymmetric aperture based on the sharpness-assessment algorithm. We estimate the depth of scene in a broad range from 60 to 350 cm by using an asymmetric aperture. The proposed algorithm recovers the depth of scene by calculating the sharpness of the deconvolved images. In our experiments, we present the estimated depth and the difference in estimation compared to the ground truth. Also, we show the absolute depth map textured with an all-focus image for several objects placed at different distances. Therefore, we demonstrate that our algorithm provides a feasible solution to recover absolute depth of scene from a single defocused image.

This paper is organized as follows. In Section 2, the camera with asymmetric aperture is modeled. In Section 3, the sharpness-assessment algorithm is described along the entire procedure. In Section 4, we show the experimental results to demonstrate our proposed algorithm, and in Section 5 the conclusion is given.

### II. Modeling for a camera with an asymmetric aperture

When we take a picture, the scene apart from the focal plane is defocused, and this blurring becomes dominant as the scene gets farther away from the focal plane. Figure 1 shows a simple camera model representing a circle of confusion (CoC) and its radius when the half of the aperture size is 7 mm. The radius is related to the position of a point source from the focal plane. As the distance from the focal plane to the point source increases, the size of the CoC also increases. In this thin lens model, the radius of the blur circle σ is given by

Figure 1.Defocus model. (a) Geometry of the camera, and (b) the radius of the CoC when half of the aperture size is 7 mm.

σ=±A(1f/zobj1f/z11),

where A is the aperture radius, f is the focal length, zobj is the object’s distance, and z1 is the focal distance. When an object is located closer than the focal plane, the sign of Eq. (1) becomes negative.

Figure 1(b) shows the radius of the CoC when the distance ranges from 60 to 350 cm. Here the focal plane of the camera is at 160 cm. The radius of the CoC changes rapidly when the distance of the object is closer than the focal plane. Depth is better distinguished when the object is closer than the focal plane. On the other hand, when the object is farther than the focal plane, the variation of the CoC is very small and it is difficult to determine the depth of the object. So, we use these facts to determine the appropriate depth range.

In this study, the concept of ADE comes from the fact that the PSF differs according to capture distance. The captured CoC using a circular aperture has a simple symmetric shape, so it is impossible to distinguish whether the point source of light is located closer or farther than the focal plane. To solve this problem, we use a camera with an asymmetric aperture shaped like the numeral 7, as shown in Fig. 2(a). When the camera’s aperture is asymmetric, we get the PSF by taking a point light source with a camera. To easily move the point light source, it is displayed on a panel. The PSFs are obtained by shifting the display panel from 60 to 350 cm in 10 cm increments. Figure 2(b) shows several examples of captured PSFs when the focus plane is set to 160 cm.

Figure 2.Camera with an asymmetric aperture. (a) The shape of an aperture like a “7”, and (b) several examples of captured PSFs, from 60 to 350 cm.

The image of a simple planar scene is computed as the convolution of the scene x and the PSF of the camera P:

I=xP+n,

where is a convolution operation and n is the noise including the camera’s shot noise and aberration. Since the image in the scene is formed by convolution of the PSF and the depth layer, the defocused image is recovered by deconvolution using the corresponding PSF for the plane of the scene. On the other hand, the defocused image is not recovered clearly with a PSF at the wrong depth. Therefore, the depth of the defocused image is estimated as the depth of the PSF that produces the sharpest image.

### III. Absolute depth estimation based on the sharpness-assessment algorithm

In this section, we show a new method of ADE based on the sharpness-assessment algorithm. The process is shown in Fig. 3. First, the defocused image is deconvolved using the PSF set. Second, its high-frequency contents are obtained using edge-detection operators. Then, the defocused image is segmented with respect to the objects in the scene. To estimate the depths of the objects, the regions are segmented using a set of masks {Ok}, where the order of segmentation is indexed by k. Third, the depth from the pre-processed images is estimated using the sharpness-assessment algorithm. It is composed of the total-sum scaling normalization, denoising based on cumulative distribution function (CDF), and scoring the sharpness of the deconvolved image. The absolute depth map and all-focus image are obtained by applying this algorithm for every segmentation respectively.

Figure 3.Flow chart for absolute depth estimation (ADE) based on the sharpness-assessment algorithm.

### 3.1. Image Deconvolution

For deconvolution of a defocused image, we use the Richardson-Lucy method that is known as a non-linear iterative deconvolution algorithm [13, 14]. It is useful for retrieving a focused image when we know the PSF of the depth layer. In our experiments, the number of iterations is set to 20. The deconvolved image using the PSF is obtained by

Jt=deconvRL(I,Pt),

where t represents the order of the depth layer, ranging from 60 to 350 cm in 10 cm increments, and is an integer from 1 to 30.

Figure 4(a) shows the defocused image of the scene where a resolution chart is located 290 cm from the camera. As shown in Fig. 4(b), we compute the reconstructed image by deconvolution with the PSF at the corresponding depth.

Figure 4.Reconstruction of a defocused resolution chart. (a) Defocused image and (b) reconstructed image when the objects are located 290 cm from the camera.

### 3.2. Edge-detection

There are many gradient operators for edge-detection, such as Sobel, Prewitt, and Roberts, but these methods are too sensitive to choose proper parameters for various features of an image. Therefore, we use four derivative edge-detection operators in the x and y directions [15]. In the x-direction, the first and second derivative operators are defined by

xJt(i,j)=[Jt(i,j)Jt(i1,j)],

xxJt(i,j)=[Jt(i+2,j)Jt(i,j)][Jt(i,j)Jt(i2,j)]

For expanding the width of the edge, the window summations with respect to ∂xJt and ∂xxJt are defined by

Wx,t,k(i,j)=Ok(i,j) iwmi+wxJt(m,j),

Wxx,t,k(i,j)=Ok(i,j) iwmi+wxxJt(m,j).

Wx,t,k represents the window summation of ∂xJt from i-w to i + w when the x-axis coordinate i is given for each segment by masking Ok. Wy,t,k and Wyy,t,k are also computed in the same manner. After applying these window summations, the features of the edges stand out. In our experiments, we set the window size, w to 2.

Figure 5 shows the window summations of 1st and 2nd derivatives along the x-axis and y-axis respectively. Images are pre-processed for the region of interest of the resolution chart from the scene. The first derivative operator is used to compute the changes in image intensity, and the second derivative is used to localize the edges.

Figure 5.Window summations using (a) 1st derivative along the x-axis, (b) 1st derivative along the y-axis, (c) 2nd derivative along the x-axis, and (d) 2nd derivative along the y-axis.

### 3.3. Total-sum Scaling Normalization

The purpose of our depth estimation algorithm is to find the sharpest image among the defocused images. In general, a sharp edge has a greater energy density than a defocused edge under the same conditions, but unfortunately some defocused edges do have more energy than sharp edges. This comes from ringing artifacts, which are caused by deconvolution with an improper convolutional kernel having a different size. These ringing artifacts increase the energy of the defocused image, so that they obstruct accurate depth estimation. Therefore, we use total-sum scaling normalization to reduce the effect of these artifacts. The normalized window summations are defined by

NWx,t,k=Wx,t,kij Wx,t,k,

NWxx,t,k=Wxx,t,kij Wxx,t,k,

### 3.4. Denoising Based on the Cumulative Distribution Function

In the high-frequency parts of the edge image, there is a lot of noise that interferes with accurately assessing the sharpness. Therefore, we remove the noise contained in the high-frequency of each image. We discover the features of the noise from the histogram. Usually, the noise in edge images is located in the region of small magnitude, and there are many such pixels, as shown in Fig. 6(a). Figure 6(b) shows the CDF, which is useful to determine a threshold point. We divide the CDF graph into 10 sections along the y-axis, and then count the number of magnitudes of each section. The threshold for a section is determined in the section where the increment of the magnitude of the next section is at least twice the increment of the current section. The representative value for the section is selected as the maximum value of the range. In most cases, the threshold point is the magnitude where the CDF value is 0.9. Thus, the denoised normalized window summations are defined by

Figure 6.Denoising based on the cumulative distribution function (CDF). (a) Example of the histogram for an edge image, and (b) threshold point determined by the CDF.

DNWx,t,k(i,j)={NWx,t,k(i,j)forNWx,t,k>Threshold0otherwise,

DNWxx,t,k(i,j)={NWxx,t,k(i,j)forNWxx,t,k>Threshold0otherwise.

Figure 7 shows the denoised normalized window summations for 1st and 2nd derivatives with respect to x and y-axes. Compared to Figs. 5(a)5(d), the noise is eliminated successfully, except for the sharp components.

Figure 7.Denoised normalized window summations from (a) 1st derivative along the x-axis, (b) 1st derivative along the y-axis, (c) 2nd derivative along the x-axis, and (d) 2nd derivative along the y-axis.

### 3.5. Scoring the Sharpness of a Deconvolved Image

Commonly, the energy of a sharp edge is greater than that of a defocused edge. Based on this fact, the scores for sharpness are defined by

Sx,t,k=ij DNWx,t,k(i,j)Nx,t,k,

Sxx,t,k=ij DNWxx,t,k(i,j)Nxx,t,k,

where Nx,t,k and Nxx,t,k present the number of non-zero pixels in the kth segment.

The scores for sharpness are undesirably affected by the size of the PSF, so we need to compensate for this effect. The compensation factor λt is defined by

λt=1σt,

where σt is the radius of the CoC at the tth depth.

The sharpness of the deconvolved image is obtained by summing the overall scores. The total score for sharpness is defined by

St,ktotal=λt(Sx,t,k+Sy,t,k+Sxx,t,k+Syy,t,k).

We estimate the absolute depth of the object as the depth of the PSF for which the total score for sharpness has a maximum value within the depth range. The absolute depth for the kth segment is given by

Dk=argmaxt(St,ktotal).

Figure 8 plots the normalized total score for sharpness with respect to the depth, where each score is normalized by the highest score. The red circle at the peak represents the score of the target position. From the graph, the sharpness score increases as the depth of the PSF becomes close to the actual depth. The highest value occurs at 290 cm. Therefore, the depth of object is estimated as 290 cm, which corresponds to the ground-truth.

Figure 8.Normalized total score for sharpness, when the resolution chart is located at 290 cm.

Experimentally, we demonstrate the proposed algorithm for two examples. First, the depth of the target is estimated as the target moves from 80 to 350 cm in 30-cm increments. Since the PSFs are captured in 10-cm increments, the depth resolution is 10 cm. The focus of the camera is fixed at 160 cm in this study. We use a Canon EOS 650D DSLR (Canon, Tokyo, Japan) with a Nikon AF-S 50 mm f/1.8G lens (Nikon, Tokyo, Japan).

As shown in Fig. 9(a), the images of the resolution chart are taken by moving the target from 80 to 350 cm in 30-cm increments. The depth of the target is estimated using the sharpness-assessment algorithm. Figure 9(b) shows the normalized scores for different target depths. Figure 9(c) shows the difference between the estimated depth and the ground- truth from 80 to 350 cm. When the target is positioned at 170 cm, the depth difference is −20 cm; this means that the estimated depth is 150 cm. The depth difference is explained by the features of the PSF: As mentioned, the radius of the PSF is too small around the focal plane, and the change in PSF is slightly smaller behind the focal plane. Therefore, the accuracy of the depth estimation near and behind the focal plane is relatively low.

Figure 9.Depth estimation of the target. (a) Illustration of the experimental setup. (b) Normalized total scores for sharpness for several depths of target. (c) Depth difference, from 80 to 350 cm.

In Fig. 9(b), the values of two peaks are slightly different. The evaluated total score at 200 cm is only about 0.2% higher than that at 150 cm. Since our algorithm estimates the absolute depth of an object as the position with the highest total score of sharpness, the position of 200 cm is chosen. However, the double peak problem is a serious problem that can mislead with erroneous estimates. We think that this ambiguity of double peaks results from the “dead-zone” which is the region near the focal plane [8]. Near the focal plane it is relatively difficult to estimate the absolute depth precisely, due to the small variation in PSF. Therefore, the accuracy of depth estimation near the focal plane is relatively low. In addition, this small PSF brings about the ambiguity in the depth estimation, and it increases the opportunity for an additional peak to appear around 160 cm.

The second experiment is conducted on a defocused image containing a painting, a cup, a photo frame, and a post box positioned at different depths as shown in Fig. 10(a). They are positioned at 100, 160, 250, and 320 cm respectively. The focal plane of the camera is again set to 160 cm. Figure 10(b) shows the segmented image that comes from using the region merging algorithm proposed by Nock and Nielsen [16]. The regions of interest are numbered sequentially according to their depths. For each segmented region, the proposed depth estimation algorithm is applied. Figure 10(c) shows the combination of all four normalized window summations obtained from the deconvolved images for the corresponding real object depths. In this image each object is indicated by bright borders, which represent the segmented areas. When the boundary of an object overlaps with that of another object placed at a different depth, the segmentation process may distort the edge-sharpness of the original image; this is regarded as one of the factors that obstructs accurate depth estimation. For that reason, only an inner feature of a segmented area is used for depth estimation.

Figure 10.Image segmentation. (a) Captured image and (b) segmented image, by applying the region merging algorithm. (c) Combination of the normalized window summations of 1st and 2nd derivatives.

Figure 11(a) shows the normalized total scores for sharpness for the four regions. As shown in Fig. 11(b), the estimated depths are identical to the actual depths except for Region 3. In Region 3, the depth difference is −10 cm. Figure 11(c) shows the depth map textured with the all-focus image. The all-focus image is reconstructed by deconvolution using the PSF at the estimated depth for each region. On the other hand, the other regions having the table and wall are set to 350 cm. Therefore, the proposed algorithm provides a feasible solution to recover the absolute depth of the scene from a single defocused image.

Figure 11.Depth estimation for a defocused image. (a) Normalized total scores for sharpness, and (b) depth difference. (c) Depth map textured with the all-focus image.

In this paper, we have presented a new method to estimate the depth of scene by using an asymmetric aperture, based on the sharpness-assessment algorithm. The asymmetric aperture is used to distinguish whether the target is located closer or farther than the focal plane. The sharpness-assessment algorithm is composed of total-sum scaling normalization, denoising based on the CDF, and scoring of the sharpness of the deconvolved image. In our experiments we used an asymmetric aperture shaped like a “7”. With the proposed method, the depth of scene was estimated over the wide range from 60 to 350 cm and the depth difference was within −20 cm, even near and behind the focal plane. Therefore, we have demonstrated that our algorithm provides a feasible solution to recover the absolute depth of scene from a single defocused image. Even though the optimization of the features of an asymmetric aperture is very important to enhance the performance of depth estimation, we have focused on demonstration of the application of an asymmetric aperture and verification of the feasibility of our sharpness-assessment algorithm. In the future, we plan to replace the conventional camera lens with a multi-aperture lens to reduce the depth difference near and behind the focal plane. We also plan to optimize the coded aperture pattern for each aperture to enhance depth discrimination resolution and accuracy.

This research was supported by ‘The Cross-Ministry Giga Korea Project’ grant funded by Korea government (MSIT) (No. 1711116979, Development of Telecommunications Terminal with Digital Holographic Table-top Display).

1. J. H. Elder and S. W. Zucker, “Local scale control for edge detection and blur estimation,” IEEE Trans. Pattern Anal. Mach. Intell. 20, 699-716 (1998).
2. M. Subbarao, T.-C. Wei and G. Surya, “Focused image recovery from two defocused images recorded with different camera settings,” IEEE Trans. Image Process 4, 1613-1628 (1995).
3. C. Swain and T. Chen, “Defocus-based image segmentation,,” in Proc. International Conference on Acoustics, Speech, and Signal Processing-ICASSP, (Detroit, MI, USA, 1995). pp. 2403-2406.
4. M. Subbarao and G. Surya, “Depth from defocus: a spatial domain approach,” Int. J. Comput. Vis. 13, 271-294 (1994).
5. Z. Djemel and D. Francois, “Depth from defocus estimation in spatial domain,” Comput. Vis. Image Underst. 81, 143-165 (2001).
6. A. P. Pentland, “A new sense for depth of field,” IEEE Trans. Pattern Anal. Mach. Intell. 9, 523-531 (1987).
7. S. Bae and F. Durand, “Defocus magnification,” Comput. Graph. Forum 26, 571-579 (2007).
8. A. Levin, R. Fergus, F. Durand and W. T. Freeman, “Image and depth from a conventional camera with a coded aperture,” ACM Trans. Graph. 26, 70-es (2007).
9. S. Zhou and S. Terence, “Defocus map estimation from a single image,” Pattern Recognit. 44, 1852-1858 (2011).
10. S. H. Lai, C. W. Fu and S. Chang, “A generalized depth estimation algorithm with a single image,” IEEE Trans. Pattern Anal. Mach. Intell. 14, 405-411 (1992).
11. C. Chen and Y. Chen, “Recovering depth from a single image using spectral energy of the defocused step edge gradient,,” in Proc 18th IEEE International Conference on Image essing-ICIP, (Brussels, Belgium, 2011). pp. 1981-1984.
12. C. Zhou, S. Lin and S. K. Nayar, “Coded aperture pairs for depth from defocus,,” in Proc. IEEE 12th International Conference on Computer Vision-ICCV, (Kyoto, Japan, 2009). pp. 325-332.
13. W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. America 62, 55-59 (1972).
14. L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. 79, 745-754 (1974).
15. J. Kumar, F. Chen and D. Doermann, “Sharpness estimation for document and scene images,,” in Proc. 21st International Conference on Pattern Recognition-ICPR2012, (Tsukuba, Japan, 2012). pp. 3292-3295.
16. R. Nock and F. Nielsen, “Statistical region merging,” IEEE Trans. Pattern Anal. Mach. Intell. 26, 1452-1458 (2004).

### Article

#### Article

Curr. Opt. Photon. 2021; 5(5): 514-523

Published online October 25, 2021 https://doi.org/10.3807/COPP.2021.5.5.514

## Absolute Depth Estimation Based on a Sharpness-assessment Algorithm for a Camera with an Asymmetric Aperture

Beomjun Kim, Daerak Heo, Woonchan Moon, Joonku Hahn

School of Electronic and Electrical Engineering, Kyungpook National University, Daegu 41566, Korea

Correspondence to:jhahn@knu.ac.kr, ORCID 0000-0002-5038-7253

Received: June 14, 2021; Revised: July 19, 2021; Accepted: July 29, 2021

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

### Abstract

Methods for absolute depth estimation have received lots of interest, and most algorithms are concerned about how to minimize the difference between an input defocused image and an estimated defocused image. These approaches may increase the complexity of the algorithms to calculate the defocused image from the estimation of the focused image. In this paper, we present a new method to recover depth of scene based on a sharpness-assessment algorithm. The proposed algorithm estimates the depth of scene by calculating the sharpness of deconvolved images with a specific point-spread function (PSF). While most depth estimation studies evaluate depth of the scene only behind a focal plane, the proposed method evaluates a broad depth range both nearer and farther than the focal plane. This is accomplished using an asymmetric aperture, so the PSF at a position nearer than the focal plane is different from that at a position farther than the focal plane. From the image taken with a focal plane of 160 cm, the depth of object over the broad range from 60 to 350 cm is estimated at 10 cm resolution. With an asymmetric aperture, we demonstrate the feasibility of the sharpness-assessment algorithm to recover absolute depth of scene from a single defocused image.

Keywords: Coded aperture, Depth estimation, Image reconstruction

### I. INTRODUCTION

Depth estimation is one of the most important areas of study in three-dimensional (3D) metrology. Depth estimation occupies a crucial position in a variety of industries, as it has been actively used for self-driving cars, inspection of defects, 3D recognition, for example. Time-of-flight and structured-light illumination are regarded as representative technologies of depth estimation, where a specific light is emitted onto a target and the reflected light signal is collected by a detector. These technologies provide a depth map with high resolution and accuracy, but they usually require expensive and complicated optical instruments.

On the other hand, depth estimation from defocused images has the great advantage that it requires only one or two images from a single conventional camera and its system has a small form factor compared to previous ones. This method originates from the relation between level of defocus (LOD) and depth profile. When the object is placed outside of the depth of field, the captured image is defocused and LOD increases proportionally to the distance from the focal plane. Based on this fact, many studies of depth recovery have concentrated on obtaining LOD from defocused images [13]. Therefore, evaluation of LOD from the scene plays a significant role in various applications such as image deblurring, image segmentation and depth estimation.

In general, methods for depth recovery are classified into relative depth estimation (RDE) and absolute depth estimation (ADE) from the defocused images. In RDE studies, the LOD is evaluated by computing the standard deviation, which is inferred from edge-detection. The standard deviation is related to the size of the point-spread function (PSF) and can be obtained from the high-frequency content using the derivative operator [4, 5]. Pentland [6] proposed the framework of RDE to obtain the focal disparity map by utilizing the relation between high frequency content coming from the difference of two defocused images. But this approach has an inconvenience: It needs to capture several images from the same scene, with different camera parameters. Bae and Durand [7] recovered the depth map from a single defocused image by employing the method of Elder and Zucker [1] and Levin et al. [8]. Zhou and Terence [9] use a ratio of gradient between the defocused input image and the estimated image to calculate the depth map. The above methods evaluate LOD of a scene from the edge-detection of a defocused image modeled by convolution of a gaussian PSF and a sharp image. Thus they take little computational time and uncomplicated processing, but it is difficult to recover absolute depth without the relation between standard deviation and depth layer [10, 11].

Most methods of ADE are based on deconvolution to recover the depth from a defocused image. The features of the PSF are important to estimate the depth precisely, and it is necessary to specially design an optical stop with a coded pattern. Levin et al. [8] proposed the framework of ADE to recover absolute depth of scene and all-focus image using a camera with a centrosymmetric coded aperture. However, the centrosymmetric coded aperture makes it difficult to distinguish if the object is placed in front of or behind the focal plane. Zhou et al. [12] proposed a method using two asymmetric apertures that are optimized for reconstruction of defocused image. They recover the depth map over a broad range and obtain an all-focus image in high quality, but they need to take two defocused images with different coded apertures for the same scene.

In this paper, we propose a new ADE method to recover the depth of scene by using an asymmetric aperture based on the sharpness-assessment algorithm. We estimate the depth of scene in a broad range from 60 to 350 cm by using an asymmetric aperture. The proposed algorithm recovers the depth of scene by calculating the sharpness of the deconvolved images. In our experiments, we present the estimated depth and the difference in estimation compared to the ground truth. Also, we show the absolute depth map textured with an all-focus image for several objects placed at different distances. Therefore, we demonstrate that our algorithm provides a feasible solution to recover absolute depth of scene from a single defocused image.

This paper is organized as follows. In Section 2, the camera with asymmetric aperture is modeled. In Section 3, the sharpness-assessment algorithm is described along the entire procedure. In Section 4, we show the experimental results to demonstrate our proposed algorithm, and in Section 5 the conclusion is given.

### II. Modeling for a camera with an asymmetric aperture

When we take a picture, the scene apart from the focal plane is defocused, and this blurring becomes dominant as the scene gets farther away from the focal plane. Figure 1 shows a simple camera model representing a circle of confusion (CoC) and its radius when the half of the aperture size is 7 mm. The radius is related to the position of a point source from the focal plane. As the distance from the focal plane to the point source increases, the size of the CoC also increases. In this thin lens model, the radius of the blur circle σ is given by

Figure 1. Defocus model. (a) Geometry of the camera, and (b) the radius of the CoC when half of the aperture size is 7 mm.

$σ=±A(1−f/zobj1−f/z1−1),$

where A is the aperture radius, f is the focal length, zobj is the object’s distance, and z1 is the focal distance. When an object is located closer than the focal plane, the sign of Eq. (1) becomes negative.

Figure 1(b) shows the radius of the CoC when the distance ranges from 60 to 350 cm. Here the focal plane of the camera is at 160 cm. The radius of the CoC changes rapidly when the distance of the object is closer than the focal plane. Depth is better distinguished when the object is closer than the focal plane. On the other hand, when the object is farther than the focal plane, the variation of the CoC is very small and it is difficult to determine the depth of the object. So, we use these facts to determine the appropriate depth range.

In this study, the concept of ADE comes from the fact that the PSF differs according to capture distance. The captured CoC using a circular aperture has a simple symmetric shape, so it is impossible to distinguish whether the point source of light is located closer or farther than the focal plane. To solve this problem, we use a camera with an asymmetric aperture shaped like the numeral 7, as shown in Fig. 2(a). When the camera’s aperture is asymmetric, we get the PSF by taking a point light source with a camera. To easily move the point light source, it is displayed on a panel. The PSFs are obtained by shifting the display panel from 60 to 350 cm in 10 cm increments. Figure 2(b) shows several examples of captured PSFs when the focus plane is set to 160 cm.

Figure 2. Camera with an asymmetric aperture. (a) The shape of an aperture like a “7”, and (b) several examples of captured PSFs, from 60 to 350 cm.

The image of a simple planar scene is computed as the convolution of the scene x and the PSF of the camera P:

$I=x⊗P+n,$

where is a convolution operation and n is the noise including the camera’s shot noise and aberration. Since the image in the scene is formed by convolution of the PSF and the depth layer, the defocused image is recovered by deconvolution using the corresponding PSF for the plane of the scene. On the other hand, the defocused image is not recovered clearly with a PSF at the wrong depth. Therefore, the depth of the defocused image is estimated as the depth of the PSF that produces the sharpest image.

### III. Absolute depth estimation based on the sharpness-assessment algorithm

In this section, we show a new method of ADE based on the sharpness-assessment algorithm. The process is shown in Fig. 3. First, the defocused image is deconvolved using the PSF set. Second, its high-frequency contents are obtained using edge-detection operators. Then, the defocused image is segmented with respect to the objects in the scene. To estimate the depths of the objects, the regions are segmented using a set of masks {Ok}, where the order of segmentation is indexed by k. Third, the depth from the pre-processed images is estimated using the sharpness-assessment algorithm. It is composed of the total-sum scaling normalization, denoising based on cumulative distribution function (CDF), and scoring the sharpness of the deconvolved image. The absolute depth map and all-focus image are obtained by applying this algorithm for every segmentation respectively.

Figure 3. Flow chart for absolute depth estimation (ADE) based on the sharpness-assessment algorithm.

### 3.1. Image Deconvolution

For deconvolution of a defocused image, we use the Richardson-Lucy method that is known as a non-linear iterative deconvolution algorithm [13, 14]. It is useful for retrieving a focused image when we know the PSF of the depth layer. In our experiments, the number of iterations is set to 20. The deconvolved image using the PSF is obtained by

$Jt=deconvRL(I,Pt),$

where t represents the order of the depth layer, ranging from 60 to 350 cm in 10 cm increments, and is an integer from 1 to 30.

Figure 4(a) shows the defocused image of the scene where a resolution chart is located 290 cm from the camera. As shown in Fig. 4(b), we compute the reconstructed image by deconvolution with the PSF at the corresponding depth.

Figure 4. Reconstruction of a defocused resolution chart. (a) Defocused image and (b) reconstructed image when the objects are located 290 cm from the camera.

### 3.2. Edge-detection

There are many gradient operators for edge-detection, such as Sobel, Prewitt, and Roberts, but these methods are too sensitive to choose proper parameters for various features of an image. Therefore, we use four derivative edge-detection operators in the x and y directions [15]. In the x-direction, the first and second derivative operators are defined by

$∂xJt(i,j)=[Jt(i,j)−Jt(i−1,j)],$

$∂xxJt(i,j)=[Jt(i+2,j)−Jt(i,j)]−[Jt(i,j)−Jt(i−2,j)]$

For expanding the width of the edge, the window summations with respect to ∂xJt and ∂xxJt are defined by

$Wx,t,k(i,j)=Ok(i,j)∑ i−w≤m≤i+w∂xJt(m,j),$

$Wxx,t,k(i,j)=Ok(i,j)∑ i−w≤m≤i+w∂xxJt(m,j).$

Wx,t,k represents the window summation of ∂xJt from i-w to i + w when the x-axis coordinate i is given for each segment by masking Ok. Wy,t,k and Wyy,t,k are also computed in the same manner. After applying these window summations, the features of the edges stand out. In our experiments, we set the window size, w to 2.

Figure 5 shows the window summations of 1st and 2nd derivatives along the x-axis and y-axis respectively. Images are pre-processed for the region of interest of the resolution chart from the scene. The first derivative operator is used to compute the changes in image intensity, and the second derivative is used to localize the edges.

Figure 5. Window summations using (a) 1st derivative along the x-axis, (b) 1st derivative along the y-axis, (c) 2nd derivative along the x-axis, and (d) 2nd derivative along the y-axis.

### 3.3. Total-sum Scaling Normalization

The purpose of our depth estimation algorithm is to find the sharpest image among the defocused images. In general, a sharp edge has a greater energy density than a defocused edge under the same conditions, but unfortunately some defocused edges do have more energy than sharp edges. This comes from ringing artifacts, which are caused by deconvolution with an improper convolutional kernel having a different size. These ringing artifacts increase the energy of the defocused image, so that they obstruct accurate depth estimation. Therefore, we use total-sum scaling normalization to reduce the effect of these artifacts. The normalized window summations are defined by

$NWx,t,k=Wx,t,k∑i∑j Wx,t,k,$

$NWxx,t,k=Wxx,t,k∑i∑j Wxx,t,k,$

### 3.4. Denoising Based on the Cumulative Distribution Function

In the high-frequency parts of the edge image, there is a lot of noise that interferes with accurately assessing the sharpness. Therefore, we remove the noise contained in the high-frequency of each image. We discover the features of the noise from the histogram. Usually, the noise in edge images is located in the region of small magnitude, and there are many such pixels, as shown in Fig. 6(a). Figure 6(b) shows the CDF, which is useful to determine a threshold point. We divide the CDF graph into 10 sections along the y-axis, and then count the number of magnitudes of each section. The threshold for a section is determined in the section where the increment of the magnitude of the next section is at least twice the increment of the current section. The representative value for the section is selected as the maximum value of the range. In most cases, the threshold point is the magnitude where the CDF value is 0.9. Thus, the denoised normalized window summations are defined by

Figure 6. Denoising based on the cumulative distribution function (CDF). (a) Example of the histogram for an edge image, and (b) threshold point determined by the CDF.

$DNWx,t,k(i,j)={NWx,t,k(i,j)for NWx,t,k>Threshold0otherwise,$

$DNWxx,t,k(i,j)={NWxx,t,k(i,j)for NWxx,t,k>Threshold0otherwise.$

Figure 7 shows the denoised normalized window summations for 1st and 2nd derivatives with respect to x and y-axes. Compared to Figs. 5(a)5(d), the noise is eliminated successfully, except for the sharp components.

Figure 7. Denoised normalized window summations from (a) 1st derivative along the x-axis, (b) 1st derivative along the y-axis, (c) 2nd derivative along the x-axis, and (d) 2nd derivative along the y-axis.

### 3.5. Scoring the Sharpness of a Deconvolved Image

Commonly, the energy of a sharp edge is greater than that of a defocused edge. Based on this fact, the scores for sharpness are defined by

$Sx,t,k=∑i∑j DNWx,t,k(i,j)Nx,t,k,$

$Sxx,t,k=∑i∑j DNWxx,t,k(i,j)Nxx,t,k,$

where Nx,t,k and Nxx,t,k present the number of non-zero pixels in the kth segment.

The scores for sharpness are undesirably affected by the size of the PSF, so we need to compensate for this effect. The compensation factor λt is defined by

$λt=1−σt,$

where σt is the radius of the CoC at the tth depth.

The sharpness of the deconvolved image is obtained by summing the overall scores. The total score for sharpness is defined by

$St,ktotal=λt(Sx,t,k+Sy,t,k+Sxx,t,k+Syy,t,k).$

We estimate the absolute depth of the object as the depth of the PSF for which the total score for sharpness has a maximum value within the depth range. The absolute depth for the kth segment is given by

$Dk=arg maxt(St,ktotal).$

Figure 8 plots the normalized total score for sharpness with respect to the depth, where each score is normalized by the highest score. The red circle at the peak represents the score of the target position. From the graph, the sharpness score increases as the depth of the PSF becomes close to the actual depth. The highest value occurs at 290 cm. Therefore, the depth of object is estimated as 290 cm, which corresponds to the ground-truth.

Figure 8. Normalized total score for sharpness, when the resolution chart is located at 290 cm.

### IV. EXPERIMENTAL RESULTS

Experimentally, we demonstrate the proposed algorithm for two examples. First, the depth of the target is estimated as the target moves from 80 to 350 cm in 30-cm increments. Since the PSFs are captured in 10-cm increments, the depth resolution is 10 cm. The focus of the camera is fixed at 160 cm in this study. We use a Canon EOS 650D DSLR (Canon, Tokyo, Japan) with a Nikon AF-S 50 mm f/1.8G lens (Nikon, Tokyo, Japan).

As shown in Fig. 9(a), the images of the resolution chart are taken by moving the target from 80 to 350 cm in 30-cm increments. The depth of the target is estimated using the sharpness-assessment algorithm. Figure 9(b) shows the normalized scores for different target depths. Figure 9(c) shows the difference between the estimated depth and the ground- truth from 80 to 350 cm. When the target is positioned at 170 cm, the depth difference is −20 cm; this means that the estimated depth is 150 cm. The depth difference is explained by the features of the PSF: As mentioned, the radius of the PSF is too small around the focal plane, and the change in PSF is slightly smaller behind the focal plane. Therefore, the accuracy of the depth estimation near and behind the focal plane is relatively low.

Figure 9. Depth estimation of the target. (a) Illustration of the experimental setup. (b) Normalized total scores for sharpness for several depths of target. (c) Depth difference, from 80 to 350 cm.

In Fig. 9(b), the values of two peaks are slightly different. The evaluated total score at 200 cm is only about 0.2% higher than that at 150 cm. Since our algorithm estimates the absolute depth of an object as the position with the highest total score of sharpness, the position of 200 cm is chosen. However, the double peak problem is a serious problem that can mislead with erroneous estimates. We think that this ambiguity of double peaks results from the “dead-zone” which is the region near the focal plane [8]. Near the focal plane it is relatively difficult to estimate the absolute depth precisely, due to the small variation in PSF. Therefore, the accuracy of depth estimation near the focal plane is relatively low. In addition, this small PSF brings about the ambiguity in the depth estimation, and it increases the opportunity for an additional peak to appear around 160 cm.

The second experiment is conducted on a defocused image containing a painting, a cup, a photo frame, and a post box positioned at different depths as shown in Fig. 10(a). They are positioned at 100, 160, 250, and 320 cm respectively. The focal plane of the camera is again set to 160 cm. Figure 10(b) shows the segmented image that comes from using the region merging algorithm proposed by Nock and Nielsen [16]. The regions of interest are numbered sequentially according to their depths. For each segmented region, the proposed depth estimation algorithm is applied. Figure 10(c) shows the combination of all four normalized window summations obtained from the deconvolved images for the corresponding real object depths. In this image each object is indicated by bright borders, which represent the segmented areas. When the boundary of an object overlaps with that of another object placed at a different depth, the segmentation process may distort the edge-sharpness of the original image; this is regarded as one of the factors that obstructs accurate depth estimation. For that reason, only an inner feature of a segmented area is used for depth estimation.

Figure 10. Image segmentation. (a) Captured image and (b) segmented image, by applying the region merging algorithm. (c) Combination of the normalized window summations of 1st and 2nd derivatives.

Figure 11(a) shows the normalized total scores for sharpness for the four regions. As shown in Fig. 11(b), the estimated depths are identical to the actual depths except for Region 3. In Region 3, the depth difference is −10 cm. Figure 11(c) shows the depth map textured with the all-focus image. The all-focus image is reconstructed by deconvolution using the PSF at the estimated depth for each region. On the other hand, the other regions having the table and wall are set to 350 cm. Therefore, the proposed algorithm provides a feasible solution to recover the absolute depth of the scene from a single defocused image.

Figure 11. Depth estimation for a defocused image. (a) Normalized total scores for sharpness, and (b) depth difference. (c) Depth map textured with the all-focus image.

### V. CONCLUSION

In this paper, we have presented a new method to estimate the depth of scene by using an asymmetric aperture, based on the sharpness-assessment algorithm. The asymmetric aperture is used to distinguish whether the target is located closer or farther than the focal plane. The sharpness-assessment algorithm is composed of total-sum scaling normalization, denoising based on the CDF, and scoring of the sharpness of the deconvolved image. In our experiments we used an asymmetric aperture shaped like a “7”. With the proposed method, the depth of scene was estimated over the wide range from 60 to 350 cm and the depth difference was within −20 cm, even near and behind the focal plane. Therefore, we have demonstrated that our algorithm provides a feasible solution to recover the absolute depth of scene from a single defocused image. Even though the optimization of the features of an asymmetric aperture is very important to enhance the performance of depth estimation, we have focused on demonstration of the application of an asymmetric aperture and verification of the feasibility of our sharpness-assessment algorithm. In the future, we plan to replace the conventional camera lens with a multi-aperture lens to reduce the depth difference near and behind the focal plane. We also plan to optimize the coded aperture pattern for each aperture to enhance depth discrimination resolution and accuracy.

### ACKNOWLEDGEMENT

This research was supported by ‘The Cross-Ministry Giga Korea Project’ grant funded by Korea government (MSIT) (No. 1711116979, Development of Telecommunications Terminal with Digital Holographic Table-top Display).

### Fig 1.

Figure 1.Defocus model. (a) Geometry of the camera, and (b) the radius of the CoC when half of the aperture size is 7 mm.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 2.

Figure 2.Camera with an asymmetric aperture. (a) The shape of an aperture like a “7”, and (b) several examples of captured PSFs, from 60 to 350 cm.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 3.

Figure 3.Flow chart for absolute depth estimation (ADE) based on the sharpness-assessment algorithm.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 4.

Figure 4.Reconstruction of a defocused resolution chart. (a) Defocused image and (b) reconstructed image when the objects are located 290 cm from the camera.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 5.

Figure 5.Window summations using (a) 1st derivative along the x-axis, (b) 1st derivative along the y-axis, (c) 2nd derivative along the x-axis, and (d) 2nd derivative along the y-axis.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 6.

Figure 6.Denoising based on the cumulative distribution function (CDF). (a) Example of the histogram for an edge image, and (b) threshold point determined by the CDF.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 7.

Figure 7.Denoised normalized window summations from (a) 1st derivative along the x-axis, (b) 1st derivative along the y-axis, (c) 2nd derivative along the x-axis, and (d) 2nd derivative along the y-axis.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 8.

Figure 8.Normalized total score for sharpness, when the resolution chart is located at 290 cm.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 9.

Figure 9.Depth estimation of the target. (a) Illustration of the experimental setup. (b) Normalized total scores for sharpness for several depths of target. (c) Depth difference, from 80 to 350 cm.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 10.

Figure 10.Image segmentation. (a) Captured image and (b) segmented image, by applying the region merging algorithm. (c) Combination of the normalized window summations of 1st and 2nd derivatives.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### Fig 11.

Figure 11.Depth estimation for a defocused image. (a) Normalized total scores for sharpness, and (b) depth difference. (c) Depth map textured with the all-focus image.
Current Optics and Photonics 2021; 5: 514-523https://doi.org/10.3807/COPP.2021.5.5.514

### References

1. J. H. Elder and S. W. Zucker, “Local scale control for edge detection and blur estimation,” IEEE Trans. Pattern Anal. Mach. Intell. 20, 699-716 (1998).
2. M. Subbarao, T.-C. Wei and G. Surya, “Focused image recovery from two defocused images recorded with different camera settings,” IEEE Trans. Image Process 4, 1613-1628 (1995).
3. C. Swain and T. Chen, “Defocus-based image segmentation,,” in Proc. International Conference on Acoustics, Speech, and Signal Processing-ICASSP, (Detroit, MI, USA, 1995). pp. 2403-2406.
4. M. Subbarao and G. Surya, “Depth from defocus: a spatial domain approach,” Int. J. Comput. Vis. 13, 271-294 (1994).
5. Z. Djemel and D. Francois, “Depth from defocus estimation in spatial domain,” Comput. Vis. Image Underst. 81, 143-165 (2001).
6. A. P. Pentland, “A new sense for depth of field,” IEEE Trans. Pattern Anal. Mach. Intell. 9, 523-531 (1987).
7. S. Bae and F. Durand, “Defocus magnification,” Comput. Graph. Forum 26, 571-579 (2007).
8. A. Levin, R. Fergus, F. Durand and W. T. Freeman, “Image and depth from a conventional camera with a coded aperture,” ACM Trans. Graph. 26, 70-es (2007).
9. S. Zhou and S. Terence, “Defocus map estimation from a single image,” Pattern Recognit. 44, 1852-1858 (2011).
10. S. H. Lai, C. W. Fu and S. Chang, “A generalized depth estimation algorithm with a single image,” IEEE Trans. Pattern Anal. Mach. Intell. 14, 405-411 (1992).
11. C. Chen and Y. Chen, “Recovering depth from a single image using spectral energy of the defocused step edge gradient,,” in Proc 18th IEEE International Conference on Image essing-ICIP, (Brussels, Belgium, 2011). pp. 1981-1984.
12. C. Zhou, S. Lin and S. K. Nayar, “Coded aperture pairs for depth from defocus,,” in Proc. IEEE 12th International Conference on Computer Vision-ICCV, (Kyoto, Japan, 2009). pp. 325-332.
13. W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. America 62, 55-59 (1972).
14. L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. 79, 745-754 (1974).
15. J. Kumar, F. Chen and D. Doermann, “Sharpness estimation for document and scene images,,” in Proc. 21st International Conference on Pattern Recognition-ICPR2012, (Tsukuba, Japan, 2012). pp. 3292-3295.
16. R. Nock and F. Nielsen, “Statistical region merging,” IEEE Trans. Pattern Anal. Mach. Intell. 26, 1452-1458 (2004).

Wonshik Choi,
Editor-in-chief


text:
Volume 10 (2014) Article 5 pp. 107-131
Competing-Provers Protocols for Circuit Evaluation
by
Revised: May 18, 2014
Published: August 9, 2014
[PDF (336K)]    [PS (1451K)]    [PS.GZ (381K)]
[Source ZIP]
Keywords: communication complexity, competing provers, delegation, complexity theory
ACM Classification: H.4, F.0
AMS Classification: 68Q10, 68Q15

Abstract: [Plain Text Version]

Let $C$ be a (fan-in $2$) Boolean circuit of size $s$ and depth $d$, and let $x$ be an input for $C$. Assume that a verifier, that knows $C$ but does not know $x$, can access the low degree extension of $x$ at one random point. Two competing provers try to convince the verifier that $C(x)=0$ and $C(x)=1$, respectively, and it is assumed that one of the provers is honest.

For any $r\geq1$, we construct1 an $r$-round protocol with communication complexity $d^{1/r}\poly\log\left(s\right)$ that convinces the verifier of the correct value of $C(x)$ (with small probability of error). In particular, when we allow only one round, the protocol exchanges $d\cdot\poly\log\left(s\right)$ bits, and when we allow $r=O\left(\frac{\log\left(d\right)}{\log\log\left(s\right)}\right)$ rounds, the protocol exchanges only $\poly\log\left(s\right)$ bits. Moreover, the complexity of the verifier and honest prover in this protocol is $\poly(s)$, and if in addition the circuit is $\log(s)$-space uniform, the complexity of the verifier is $d^{1/r}\poly\log\left(s\right)$. The protocol is obtained by combining the delegation protocol of Goldwasser, Kalai, and Rothblum (STOC 2008), the competing-provers protocols of Feige and Kilian (STOC 1997), and some new techniques. We suggest two applications of these results:

Delegating computation to competing clouds: The main motivation behind the protocol of GKR'08 was delegating computation to a cloud. Using our new protocol, a verifier can delegate computation to two (or more) competing clouds. If at least one of the clouds is reliable the verifier can trust that the computation is correct (with high probability). The advantage over the protocol of GTR'08 is that the communication complexity and the number of rounds in our protocol are significantly lower.

Communication complexity with competing provers, and circuit lower bounds: Aaronson and Wigderson (2009) suggested the model of communication complexity with competing provers, where two competing provers try to convince two players that $f(x,y)=0$ and $f(x,y)=1$, respectively, where $x$ is an input held by the first player and $y$ is an input held by the second player. By scaling down the competing-provers protocols of FK'97, they showed that strong enough lower bounds for the communication complexity of $f$, in this model, imply lower bounds for the computational complexity of $f$.

Our results strengthen this connection. More precisely, we show that if $f$ can be computed by a Boolean circuit of size $s$ and depth $d$ then for any $r\geq1$ there is an $r$-round protocol for $f$, in this model, with communication complexity $d^{1/r}\poly\log\left(s\right)$. This can be viewed as a possible direction towards proving circuit lower bounds. For instance, in order to prove $f\notin\NC$, it suffices to show that any one round protocol for $f$, in this model, requires the exchange of $\omega\left(\poly\log\left(n\right)\right)$ bits. This gives a relatively simple combinatorial property that implies strong circuit lower bounds.

A conference version of this paper appeared in the Proceedings of the 4th Innovations in Theoretical Computer Science Conference (ITCS), 2013.

1 We have learned that a result similar to ours for the case $r=1$, as well as the motivation of delegating computation to competing clouds, was independently obtained by Canetti, Riva, and Rothblum (ICITS 2012).

format: plain


text:
# OpenStax_CNX

You are here: Home » Content » SHM equation

### Recently Viewed

This feature requires Javascript to be enabled.

# SHM equation

Module by: Sunil Kumar Singh. E-mail the author

Note: You are viewing an old version of this document. The latest version is available here.

Force is the cause of simple harmonic motion. However, it plays a very different role than that in translational or rotational motion. First thing that we should be aware that oscillating systems like block-spring arrangement or a pendulum are examples of systems in stable equilibrium. When we apply external force, it tends to destabilize or disturb the state of equilibrium by imparting acceleration to the object in accordance with Newton’s second law.

Secondly, role of destabilizing external force is one time act. The oscillation does not require destabilizing external force subsequently. It, however, does not mean that SHM is a non-accelerated motion. As a matter of fact, oscillating system generates a restoring mechanism or restoring force that takes over the role of external force once disturbed. The name “restoring” signifies that force on the oscillating object attempts to restore equilibrium (central position in the figure) – undoubtedly without success in SHM.

SHM is an accelerated motion in which object keeps changing its velocity all the time. The analysis of SHM involves consideration of “restoring force” – not the external force that initially starts the motion. Further, we need to understand that initial external force and hence restoring force are relatively small than the force required to cause translation or rotation. For example, if we displace pendulum bob by a large angle and release the same for oscillation, then the force on the system may not fulfill the requirement of SHM and as such the resulting motion may not be a SHM.

We conclude the discussion by enumerating requirements of SHM as :

• The object is in stable equilibrium before start of the motion.
• External destabilizing force is applied only once.
• The object accelerates and executes SHM under the action of restoring force.
• The magnitude of restoring force or displacement is relatively small.
• There is no dissipation of energy during motion (ideal reference assumption).

## Force equation

Here, we set out to figure out nature of restoring force that maintains SHM. For understanding purpose we consider the block-spring system and analyze “to and fro” motion of the block. Let the origin of reference coincides with the position of the block for the neutral length of spring. The block is moved right by a small displacement “x” and released to oscillate about neutral position or center of oscillation. The restoring spring force is given by (“k” is spring constant) :

F = - k x F = - k x

In the case of pendulum, we describe motion in terms of torque as it involves angular motion. Here, torque is (we shall study this relation later):

τ = - m g l θ τ = - m g l θ

In either case, we see that “cause” (whether force or torque) is proportional to negative of displacement – linear or angular as the case be. Alternatively, we can also state the nature of restoring force in terms of acceleration,

Linear acceleration, a = - k m x a = - k m x ------ for linear SHM of block-spring system

Angular acceleration, τ = - m g l I θ τ = - m g l I θ ------ for rotational SHM of pendulum

where “m” and “I” are mass and moment of inertia of the oscillating objects in two systems.

In order to understand the nature of cause, we focus on the block-spring system. When the block is to the right of origin, “x” is positive and restoring spring force is negative. This means that restoring force (resulting from elongation of the spring) is directed left towards the neutral position (center of oscillation). This force accelerates the block towards the center. As a result, the block picks up velocity till it reaches the center.

The plot, here, depicts nature of force about the center of oscillation bounded between maximum displacements on either side.

As the block moves past the center, “x” is negative and force is positive. This means that restoring force (resulting from compression of the spring) is directed right towards the center. The acceleration is positive, but opposite to direction of velocity. As such restoring force decelerates the block.

In the nutshell, after the block is released at one extreme, it moves, first, with acceleration up to the center and then moves beyond center towards left with deceleration till velocity becomes zero at the opposite extreme. It is clear that block has maximum velocity at the center and least at the extreme positions (zero).

From the discussion, the characterizing aspects of the restoring force responsible for SHM are :

• The restoring force is always directed towards the center of oscillation.
• The restoring force changes direction across the center.
• The restoring force first accelerates the object till it reaches the center and then decelerates the object till it reaches the other extreme.
• The process of “acceleration” and “deceleration” keeps alternating in each half of the motion.

## Differential form of SHM equation

We observed that acceleration of the object in SHM is proportional to negative of displacement. Here, we shall formulate the general equation for SHM in linear motion with the understanding that same can be extended to SHM along curved path. In that case, we only need to replace linear quantities with corresponding angular quantities.

a = - ω 2 x a = - ω 2 x

where “ ω 2 ω 2 ” is a constant. The constant “ω” turns out to be angular frequency of SHM. This equation is the basic equation for SHM. For block-spring system, it can be seen that :

ω = k m ω = k m

where “k” is the spring constant and “m” is the mass of the oscillating block. We can, now, write acceleration as differential,

2 x t 2 = - ω 2 x 2 x t 2 = - ω 2 x

2 x t 2 + ω 2 x = 0 2 x t 2 + ω 2 x = 0

This is the SHM equation in differential form for linear oscillation. A corresponding equation of motion in the context of angular SHM is :

2 θ t 2 + ω 2 θ = 0 2 θ t 2 + ω 2 θ = 0

where "θ" is the angular displacement.

## Solution of SHM differential equation

In order to solve the differential equation, we consider position of the oscillating object at an initial displacement x 0 x 0 at t =0. We need to emphasize that “ x 0 x 0 ” is initial position – not the extreme position, which is equal to amplitude “A”. Let

t = 0, x = x 0, v = v 0 t = 0, x = x 0, v = v 0

We shall solve this equation in two parts. We shall first solve equation of motion for the velocity as acceleration can be written as differentiation of velocity w.r.t to time. Once, we have the expression for velocity, we can solve velocity equation to obtain a relation for displacement as its derivative w.r.t time is equal to velocity.

### Velocity

We write SHM equation as differential of velocity :

a = v t = = - ω 2 x a = v t = = - ω 2 x

v x X x t = - ω 2 x v x X x t = - ω 2 x

v v x = - ω 2 x v v x = - ω 2 x

Arranging terms with same variable on either side, we have :

v v = - ω 2 x x v v = - ω 2 x x

Integrating on either side between interval, while keeping constant out of the integral sign :

v 0 v v v = - ω 2 x 0 x x x v 0 v v v = - ω 2 x 0 x x x

[ v 2 2 ] v 0 v = - ω 2 [ x 2 2 ] x 0 x [ v 2 2 ] v 0 v = - ω 2 [ x 2 2 ] x 0 x

v 2 v 0 2 = - ω 2 x 2 x 0 2 v 2 v 0 2 = - ω 2 x 2 x 0 2

v = { v 0 2 + ω 2 x 0 2 ω 2 x 2 } v = { v 0 2 + ω 2 x 0 2 ω 2 x 2 }

v = ω { v 0 2 ω 2 + x 0 2 x 2 } v = ω { v 0 2 ω 2 + x 0 2 x 2 }

We put v 0 2 ω 2 + x 0 2 = A 2 v 0 2 ω 2 + x 0 2 = A 2 . We shall see that “A” turns out to be the amplitude of SHM.

v = ω A 2 x 2 v = ω A 2 x 2

This is the equation of velocity. When x = A or -A,

v = ω A 2 A 2 = 0 v = ω A 2 A 2 = 0

when x = 0,

v max = ω A 2 0 2 = ω A v max = ω A 2 0 2 = ω A

### Displacement

We write velocity as differential of displacement :

v = x t = ω A 2 x 2 v = x t = ω A 2 x 2

Arranging terms with same variable on either side, we have :

x A 2 x 2 = ω t x A 2 x 2 = ω t

Integrating on either side between interval, while keeping constant out of the integral sign :

x 0 x x A 2 x 2 = ω 0 t t x 0 x x A 2 x 2 = ω 0 t t

[ sin 1 x A ] x 0 x = ω t [ sin 1 x A ] x 0 x = ω t

sin 1 x A sin 1 x 0 A = ω t sin 1 x A sin 1 x 0 A = ω t

Let sin 1 x 0 A = φ sin 1 x 0 A = φ . We shall see that “φ” turns out to be the phase constant of SHM.

sin 1 x A = ω t + φ sin 1 x A = ω t + φ

x = A sin ω t + φ x = A sin ω t + φ

This is one of solutions of the differential equation. We can check this by differentiating this equation twice with respect to time to yield equation of motion :

x t = A ω cos ω t + φ x t = A ω cos ω t + φ

2 x t 2 = A ω 2 sin ω t + φ = ω 2 x 2 x t 2 = A ω 2 sin ω t + φ = ω 2 x

2 x t 2 + ω 2 x = 0 2 x t 2 + ω 2 x = 0

Similarly, it is found that equation of displacement in cosine form, x = A cos ω t + φ x = A cos ω t + φ , also satisfies the equation of motion. As such, we can use either of two forms to represent displacement in SHM. Further, we can write general solution of the equation as :

x = A sin ω t + B cos ω t x = A sin ω t + B cos ω t

This equation can be reduced to single sine or cosine function with appropriate substitution.

## Example

Problem 1: Find the time taken by a particle executing SHM in going from mean position to half the amplitude. The time period of oscillation is 2 s.

Solution : Employing expression for displacement, we have :

x = A sin ω t x = A sin ω t

We have deliberately used sine function to represent displacement as we are required to determine time for displacement from mean position to a certain point. We could ofcourse stick with cosine function, but then we would need to add a phase constant “π/2” or “-π/2”. The two approach yields the same expression of displacement as above.

Now, according to question,

A 2 = A sin ω t A 2 = A sin ω t

sin ω t = 1 2 = sin π 6 sin ω t = 1 2 = sin π 6

ω t = π 6 ω t = π 6

t = π 6 ω = π T 6 X 2 π = T 12 = 2 12 = 1 6 s t = π 6 ω = π T 6 X 2 π = T 12 = 2 12 = 1 6 s

## Content actions

### Give feedback:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks


text:
1
0
Entering edit mode
12 months ago
renaegeier ▴ 10

I've been having issues getting snakemake to work. Here is the first rule that will not run. I want to take protein annotation files and make blast databases out of them.

STRAINS = ["A", "B", "C"]

rule make_db:
input:
"prokka_faa/{strain}.faa"
output:
"dbs/{strain}-db"
conda:
"envs/finding.yaml"
shell:
"makeblastdb -input_type prot -dbtype prot -in {input} -out {output}"

But when I run this:

$snakemake -n --snakefile snakefile.txt --cores 2 -p I get the output "Building DAG of jobs... WorkflowError: Target rules may not contain wildcards. Please specify concrete files or a rule without wildcards." I can get this to run with just one strain, but I need to run this on 185 files. This may be a simple fix but I have tried everything I can think of to fix this and I'm not great with python syntax. Any help is greatly appreciated, thank you! snakemake snakefile • 515 views ADD COMMENT 4 Entering edit mode 12 months ago Snakemake doesn't know what strains you want it to make dbs/*-db for. You could tell it on the command line, i.e., specify the targets concretely:$ snakemake -n --snakefile snakefile.txt --cores 2 -p dbs/A-db dbs/B-db dbs/C-db

or you can add an umbrella rule that fans out to all the strains. Make it the first rule in the snakefile as that is the default target rule used when you don't specify any on the command line:

rule make_db_all:
input:
expand("dbs/{strain}-db", strain = STRAINS)
0
Entering edit mode

Thank you!! This got the command to run! Do I need to make a rule like that before every rule that uses wildcards?


text:
# Source of a river

In this example we impose a variable flow rate for a “source” located inside the computation domain.

``````#include "saint-venant.h"
#include "discharge.h"``````

The domain is 10 metres squared, centered on the origin. Time is in seconds.

``````#define LEVEL 8

int main()
{
size (10.);
origin (- L0/2., - L0/2.);
G = 9.81;
N = 1 << LEVEL;
run();
}``````

## Boundary conditions

We create a new boundary for the source, with a Neumann condition for the normal velocity (i.e. an inflow).

``````bid source;
u.n[source] = neumann(0);``````

The flow rate varies in time and is set by computing the the elevation ${\eta }_{s}$ of the water surface necessary to match this flow rate.

``````double etas;

event inflow (i++) {
etas = eta_b (0.1*(1.1 - cos(4.*π*t)), source);
h[source] = max (etas - zb[], 0);
η[source] = max (etas - zb[], 0) + zb[];
}``````

## Initial conditions

The river bed is a single valley. The source is created by masking and is a narrow slot located at the head of the valley.

``````event init (i = 0)
{
mask (fabs(x) < 0.5 && fabs(y - 3.5) < Δ/2. ? source : none);``````

We start with a dry riverbed, so that the problem does not have a natural timescale the Saint-Venant solver can use. We set a maximum timestep to set this timescale.

``````  DT = 1e-2;

foreach()
zb[] = (- cos(x) + y)/2.;
boundary ({zb});
}``````

## Outputs

We compute the time-derivative of the total water volume (i.e. the net flow rate), and make a GIF movie.

``````event logfile (i++; t <= 2.) {
static double volo = 0., to = 0.;
double vol = statsf(h).sum;
if (i > 0)
fprintf (ferr, "%g %.6f %g %g\n", t, (vol - volo)/(t - to), vol, etas);
volo = vol, to = t;
}

event output (i += 5) {
static FILE * fp = popen ("ppm2gif > source.gif", "w");
output_ppm (h, fp, min = 0, max = 0.05);
}``````


text:
Gold processing | Britannica

Gold processing, preparation of the ore for use in various products.. For thousands of years the word gold has connoted something of beauty or value. These images are derived from two properties of gold, its colour and its chemical stability.The colour of gold is due to the electronic structure of the gold atom, which absorbs electromagnetic radiation with wavelengths less than 5600 angstroms ...

Ore processing – Waihi Gold

Typically, 80% of the ore particles are less than 70 microns in size (about the same consistency as fine beach sand). Such a fine particle size is required for gold liberation – the size required for the cyanide to be able to ‘see’ the gold in the leaching process.

Gold processing 101 - Mining Magazine

The hashing stage (corresponding to metal extraction and recovery stages) is a little more complex for gold ores, as the optimal process flowsheet selection choice is heavily dependent on a good understanding of two fundamental geometallurgical parameters, the gold mineralogical associations, and the gold particle size and liberation ...

mm size product lime thickener gold ore

lime gold ore ball mill pictures . Ore processingWaihi Gold Ore is stockpiled 1 at the processing plant and the process begins by feeding the ore into a hopper with a loader The ore is conveyed and lime is added 2 to raise the pH of the ore Following crushing through a jaw crusher 3 the ore is fed into the semi autogenous grinding SAG mill 4 along with water and steel balls

Processing, smelting and refining gold | World Gold Council

In cyanidation, metallic gold is oxidised and dissolved in an alkaline cyanide solution. When gold dissolution is complete, the gold-bearing solution is separated from the solids. With ores of higher gold content (greater than 20 grams of gold per tonne of ore), cyanidation is accomplished by vat leaching, which involves holding a slurry of ore ...

Gold processing - Mining and concentrating | Britannica

Gold processing - Gold processing - Mining and concentrating: The nature of the ore deposit determines the mining and mineral processing techniques applied. Oxide ore deposits are frequently of such low grade (e.g., 3 to 10 parts per million) that extensive mineral processing cannot economically be justified. In this case they are merely shattered by explosives and then piled into heaps for ...

Gold Processing,Extraction,Smelting Plant Design ...

2019/03/06 · Gold ore . Prominer maintains a team of senior gold processing engineers with expertise and global experience. These gold professionals are specifically in gold processing through various beneficiation technologies, for gold ore of different characteristics, such as flotation, cyanide leaching, gravity separation, etc., to achieve the processing plant of optimal and cost-efficient process designs.

Gold Processing Methods & Gold Ore Extraction

Of all the methods of extracting gold & processing it from its ore, I used a few to evaluate two principal flowsheets in this case study. The flowsheets utilized operations that involved flotation, cyanidation and gravity concentration. Tests that mirror each of these unit operations were utilized to evaluate the principal flowsheets. This page offers a comparative review of gold recovery ...

Production of Gold

Process Description Unit 100 – Size Reduction of Ore The BFD of the overall process is shown in Figure 1. The PFD for Unit 100, shown in Figure 2, is designed to reduce 41.5 tons/hr of gold ore from a feed range of 2-5” to 160 microns. The mined ore is fed using a Grizzly Feeder, F-101, into a Jaw

How Gold Is Mined | Gold Mining Process | World Gold Council

The lifecycle of a gold mine. People in hard hats working underground is what often comes to mind when thinking about how gold is mined. Yet mining the ore is just one stage in a long and complex gold mining process. Long before any gold can be extracted, significant exploration and development needs to take place, both to determine, as accurately as possible, the size of the deposit as well ...

Mining and mineral processing | The lithosphere | Siyavula

Processing the gold ore (ESBRM) For every ton of ore that is mined, only a very small amount (about $$\text{5}$$ $$\text{g}$$) of gold is extracted. A number of different methods can be used to separate gold from its ore, but one of the more common methods is called gold cyanidation .

The process design of gold leaching and carbon-in-pulp ...

The CIP process A block-flow diagram of a typical CIP plant for a non-refractory gold ore is shown in Figure 2. Table I and Table II illustrate the capital and operating cost breakdowns for a typical South African gold plant. These figures are not a standard but reflect the nature of the ore and the design basis for a particular situation.

Gold extraction - QueensMineDesignWiki

Through various pre-treatment options, the goal is to achieve sufficient gold recovery in the most economical way possible. This article serves its purpose in providing an overall guide to gold processing, and the associated processing options for each gold ore classification. References. Brooy, S. R. (1994). Review of gold extraction from ores.

Gold particle size distribution in Au deposits - Assaying ...

Hi there, I'm completely new to 911 Metallurgist, but maybe someone here can help me out! I'm looking for studies that characterize the gold particle size distribution in different geological deposit types (I'm a geologist). I would like to know how the size distribution of Au particles varies in

Mineral Processing - SlideShare

2014/05/15 · Used in gold ore grinding. CASCADE AND AEROFALL MILL- Have larger dia. Dia/Length > 3. The lumps fall through a large distance. The feed and discharge are through trunnions. The mill rotates fast. High reduction ratio. Used for iron and gold ores. 120 Feed+ air to cyclone Drive 115. 121 AEROFALL MILL- 116. 122 117.

Day 51, Procedures to obtain gold: Vat Leaching Process ...

2013/07/29 · The Vat Leaching Process Vat leaching can be described as leaching in a bathtub inside a building. Treatment of ore is site-specific as determined by metallurgical tests. Bench scale metallurgical testing of ores for their amenability to vat leaching is similar to that for heap leaching. First, bottle-roll tests of finely ground material are run…

Gold Processing at Boddington Mine - Blogger

Ore processing stage in gold processing plant at Boddington mine through the stages of crushing . The Boddington primary crushing circuit comprises two 60” x 113” XHD primary gyratory crushers located at the edge of the mine with a nominal capacity of 3,670 t/h each at an OSS (open side setting) of 175mm.

Gold Mining in South Africa | Africa Mining IQ

Over 50% of all gold reserves are found in South Africa, with the Witwatersrand Basin remaining the largest gold resource in the world.. In 1975, South Africa was responsible for producing 40% of the gold ever mined. By 2010, however, China affirmed its status as the world’s largest gold producer with production of 324 tonnes of gold, followed by Australia (222.8 tonnes), and South Africa ...


text:
# Convenient error analysis

Is there any convenient way to find out the error in evaluating an expression numerically, providing all errors of the arguments are given? Like

Sin[DataError[123, .1] + 3*DataError[456, .2]/DataError[789, .5]]^2

DataError[xxx, xx]

• You might consider using Interval[] for this application. – J. M. will be back soon Dec 12 '16 at 14:51
• What would be the form of the result you have in mind? – Daniel Lichtblau Dec 12 '16 at 15:51
• @DanielLichtblau , like what's written above, plug in all tuples of data and it's error, giving out a tuple of result and calculated error (according to error analysis and pde stuff) – ZisIsNotZis Dec 12 '16 at 15:54

It was not stated how DataError is to be defined, and in particular whether the second argument is to be taken as an absolute or relative error. I am assuming absolute.

Easiest is, as suggested by @J.M., to use Interval.

dataError1[a_, b_] := Interval[{a - b, a + b}]
val1 = Sin[
dataError1[123, .1] + 3*dataError1[456, .2]/dataError1[789, .5]]^2

(* Out[580]= Interval[{0.542555735765, 0.736481835669}] *)

Could use explicit optimization to avoid the "dependency problem" (which does not arise in this example).

dataError2[a_, b_, eps_] := a + b*eps
val2 = Sin[
dataError2[123, .1, eps1] +
3*dataError2[456, .2, eps2]/dataError2[789, .5, eps3]]^2;

evars = {eps1, eps2, eps3};
NMinimize[{val2, Thread[-1 <= evars <= 1]}, evars]
NMaximize[{val2, Thread[-1 <= evars <= 1]}, evars]

(* Out[584]= {0.542555735765, {eps1 -> 1., eps2 -> 1., eps3 -> -1.}}

Out[585]= {0.736481835669, {eps1 -> -1., eps2 -> -1., eps3 -> 1.}} *)

If the error offsets are fairly tight, could do this as a first-order estimate. Much better for speed, in complicated cases.

firstorder =
Normal[Series[val2 /. Thread[evars -> t*evars], {t, 0, 1}]] /. t -> 1
NMinimize[{firstorder, Thread[-1 <= evars <= 1]}, evars]
NMaximize[{firstorder, Thread[-1 <= evars <= 1]}, evars]

(* Out[586]= -0.095854777084 eps1 - 0.000728933666038 eps2 +
0.00105321214104 eps3 + Sin[32805/263]^2

Out[587]= {0.544829012526, {eps1 -> 1., eps2 -> 1., eps3 -> -1.}}

Out[588]= {0.740102858308, {eps1 -> -1., eps2 -> -1., eps3 -> 1.}} *)

One could use the built-in arbitrary precision numerics, I suppose. I don't know the accuracy of Accuracy[], but it seems to be close if an approximate error is what is desired and not hard bounds. The smaller the relative error, the better it will be. Accuracy is measure on a log-10 scale. See also SetAccuracy[].

Sin[SetAccuracy[123, 10^-.1] +
3*SetAccuracy[456, 10^-.2]/SetAccuracy[789, 10^-.5]]^2
% + {1, -1} Log10[Accuracy[%]]
(*
0.*10^-1
{0.549593, 0.735339}
*)

We can see from the 0.*10^-1, whose FullForm is 0.6424659354173337650.6153213200244898`, is that one drawback is that the error can overwhelm the computed value. If it does, Mathematica will return a value of zero and this method will fail.


text:
# Scaling up NTC measurement through Multi-Channel Ohm meter

For a test bed, we want to measure temperatures with NTC sensors (~10kOhm at RT) for the range 10°C to 70°C. Our test machine accepts auxiliary signal through CAN-bus, so we used a Raspberry with a "CAN-hat" to prototype the temperature measurement (three-point probing) and configure the test machine.

We now want to scale it up to 80 (phase 1) then 320 (phase 2) channels. We already have all the necessary NTC sensors (recovered from previous projects). So far we have 2 solutions which are not satisfying:

• An all inclusive solution by a provider for a hefty price

• Acquiring eight 10-channel Ohm meters + eight RS-232 to CAN. That's about $100 to$150 per channel. That would slash the price in 2 compared to solution 1.

In both cases, that's more than the cost of the RPi + the CAN-hat per channel. I would expect that scaling up a prototype would mean reducing the price per unit.

Is there any less resource-intensive solution(s) we are missing? Ideally measurement with 3 or 4-point method would be wished.

Thank you. edit: accuracy around 1K, resolution around 0.1K.

• What is your required resolution and accuracy? How much time do you have to take these readings? – Elliot Alderson Aug 13 '18 at 16:27
• Have you considered Analog Switch Mux? – Sunnyskyguy EE75 Aug 13 '18 at 18:04
• Thanks for your questions. Resolution should lie around 0.1K and accuracy around 1K. Acquisition frequency is typically 30s to few minutes. I have not heard of Analog Switch Mux. – Wli Aug 14 '18 at 7:29

At the quantities you desire, it would most likely be best to design a custom PCB with a micro with CAN bus, and ADC and an OPamp could be had for the 50$-70$ range depending on your requirements. (0.1C 0.5C or 1C accuracy?)

The circuit looks like this (doesn't matter if it's a PTC or NTC): Source: Here

An add on board with analog for an RPI would probably be in the same price range. If you did multiple channels on the same board the cost would go down even more. Depending on the accuracy and temperature resolution, maybe less. There would be an additional cost of time and software writing for the custom solution.

Another option that might work is the LTC2986 which has a demo board with an 'arudino' (actually linear's own flavor), if I remember right they have at least 4 channels (or 8) with 0.1C accuracy.

Source: Digikey.com

With a one of these solutions the NTC's may also reqire calibration (depending on your accuracy and resolution requirements) with a bath or some other method which might be more complexity

• Thanks for the great answer! Calibration is what scares me the most, for 320 channels it'll be very repetitive and long. But our other "all-in" solutions have the same requirements. Accuracy should lie around 1K, resolution 0.1K. – Wli Aug 14 '18 at 7:51


text:
# How can Lagrange Multipliers and Penalty Method be applied to optimize algorithms?

I have a programming assignment which I was told that is solvable with some DP algorithm. The question is a variant of LIS problem where at most $$k$$ exceptions (restarts) are allowed.

But I know that there is a better solution. My professor mentioned Lagrange Multipliers and giving a penalty for each restart. But after googling these terms I wasn't able to find out anything related to algorithms. Every article is about Calculus and function optimization. I also read about them on Wikipedia but I can't figure out how to use them.

Is there a keyword that can describe better what I want to read?

This technique is called Lagrangian relaxation.

The regular $$DP$$ approach, where $$DP[a][b]$$ represents the length of the longest increasing subsequence that ends in the $$a$$'th number and restarts at most $$b$$ times, is $$\mathcal{O}(nk \log n)$$. For convenience we'll assume the last number is the largest, and therefore $$DP[n][k]$$ is the value we are looking for (if this isn't the case, append $$\infty$$ and after calculating the answer decrement it by $$1$$).

To optimize this, we'll select some $$\lambda \in \mathbb{N}$$ which represents the cost of every exception, and compute $$DP'[a] = \max_{b} DP[a][b] - \lambda b$$. This can be done in $$\mathcal{O}(n \log n)$$: first sort the values, and keep a range maximum data structure over them, with all positions $$j$$ initialised to $$v_{j} = 0$$. Assume the value at position $$i$$ is the $$p_{i}$$th in the sorted list of values. Then $$DP'[i] = \max(1 + \max_{j < p_{i}} v_{j}, 1 - \lambda + \max_{j > p_{i}} v_{j})$$, and we set $$v_{p_{i}} = DP'[i]$$.

What use is computing $$\max_{b} DP[n][b] - \lambda b$$ to us? Notice that as we increase $$\lambda$$, the optimal $$b$$ in the maximum cannot increase. Similarly as we decrease $$\lambda$$, the optimal $$b$$ cannot decrease. If $$\lambda = 0$$, it is optimal to take all elements to our subsequence, and if $$\lambda = n$$, it is optimal to have $$0$$ exceptions. If we can find $$\lambda$$ for which the optimal $$b$$ can be $$k$$, then $$DP[n][b] = DP'[n] + \lambda b$$. Further, if such $$\lambda$$ exists, we can do binary search for it, yielding a $$\mathcal{O}(n \log^{2} n)$$ algorithm.

Note that we can modify the $$\mathcal{O}(n \log n)$$ algorithm to calculate the minimum and maximum values of $$b$$ that achieve the maximum value with the specific $$\lambda$$ with no increase in complexity. We can always find a $$\lambda$$ for which $$\min_{b} \leq k \leq \max_{b}$$, but this doesn't guarantee that there exists some subsequence with $$k$$ exceptions achieving the maximum. However, if we can show that $$DP[n][b]$$ is concave, i.e. $$DP[n][b+2] - DP[n][b+1] \leq DP[n][b+1] - DP[n][b]$$, we get this result, as we know that $$DP[n][\min_{b} + 1] \leq DP[n][\min_{b}] + \lambda$$ (due to maximality), therefore $$DP[n][\max_{b}] \leq DP[n][k] + (\max_{b} - k) \lambda$$, hence $$DP'[n] = DP[n][\max_{b}] - \lambda \max_{b} \leq DP[n][k] - \lambda k$$ and there exists a subsequence with $$k$$ exceptions achieving the maximum.

EDIT: Here is a C++ program for finding a maximal subsequence in $$\mathcal{O}(n \log^{2} n)$$. I use a segment tree for the range maximum data structure.

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
using ll = long long;
const int INF = 2 * (int)1e9;

pair<ll, pair<int, int>> combine(pair<ll, pair<int, int>> le, pair<ll, pair<int, int>> ri) {
if (le.first < ri.first) swap(le, ri);
if (ri.first == le.first) {
le.second.first = min(le.second.first, ri.second.first);
le.second.second = max(le.second.second, ri.second.second);
}
return le;
}

// Specialised range maximum segment tree
class SegTree {
private:
vector<pair<ll, pair<int, int>>> seg;
int h = 1;

pair<ll, pair<int, int>> recGet(int a, int b, int i, int le, int ri) const {
if (ri <= a || b <= le) return {-INF, {INF, -INF}};
else if (a <= le && ri <= b) return seg[i];
else return combine(recGet(a, b, 2*i, le, (le+ri)/2), recGet(a, b, 2*i+1, (le+ri)/2, ri));
}
public:
SegTree(int n) {
while(h < n) h *= 2;
seg.resize(2*h, {-INF, {INF, -INF}});
}
void set(int i, pair<ll, pair<int, int>> off) {
seg[i+h] = combine(seg[i+h], off);
for (i += h; i > 1; i /= 2) seg[i/2] = combine(seg[i], seg[i^1]);
}
pair<ll, pair<int, int>> get(int a, int b) const {
return recGet(a, b+1, 1, 0, h);
}
};

// Binary searches index of v from sorted vector
int bins(const vector<int>& vec, int v) {
int low = 0;
int high = (int)vec.size() - 1;
while(low != high) {
int mid = (low + high) / 2;
if (vec[mid] < v) low = mid + 1;
else high = mid;
}
return low;
}

// Finds longest strictly increasing subsequence with at most k exceptions in O(n log^2 n)
vector<int> lisExc(int k, vector<int> vec) {
// Compress values
vector<int> ord = vec;
sort(ord.begin(), ord.end());
ord.erase(unique(ord.begin(), ord.end()), ord.end());
for (auto& v : vec) v = bins(ord, v) + 1;

// Binary search lambda
int n = vec.size();
int m = ord.size() + 1;
int lambda_0 = 0;
int lambda_1 = n;
while(true) {
int lambda = (lambda_0 + lambda_1) / 2;
SegTree seg(m);
if (lambda > 0) seg.set(0, {0, {0, 0}});
else seg.set(0, {0, {0, INF}});

// Calculate DP
vector<pair<ll, pair<int, int>>> dp(n);
for (int i = 0; i < n; ++i) {
auto off0 = seg.get(0, vec[i]-1); // previous < this
off0.first += 1;

auto off1 = seg.get(vec[i], m-1); // previous >= this
off1.first += 1 - lambda;
off1.second.first += 1;
off1.second.second += 1;

dp[i] = combine(off0, off1);
seg.set(vec[i], dp[i]);
}

// Is min_b <= k <= max_b?
auto off = seg.get(0, m-1);
if (off.second.second < k) {
lambda_1 = lambda - 1;
} else if (off.second.first > k) {
lambda_0 = lambda + 1;
} else {
// Construct solution
ll r = off.first + 1;
int v = m;
int b = k;
vector<int> res;
for (int i = n-1; i >= 0; --i) {
if (vec[i] < v) {
if (r == dp[i].first + 1 && dp[i].second.first <= b && b <= dp[i].second.second) {
res.push_back(i);
r -= 1;
v = vec[i];
}
} else {
if (r == dp[i].first + 1 - lambda && dp[i].second.first <= b-1 && b-1 <= dp[i].second.second) {
res.push_back(i);
r -= 1 - lambda;
v = vec[i];
--b;
}
}
}
reverse(res.begin(), res.end());
return res;
}
}
}

int main() {
int n, k;
cin >> n >> k;

vector<int> vec(n);
for (int i = 0; i < n; ++i) cin >> vec[i];

vector<int> ans = lisExc(k, vec);
for (auto i : ans) cout << i+1 << ' ';
cout << '\n';
}

EDIT2: Thanks to Jaehyun Koo over at Codeforces I now know how to show concavity. The following is a modified version of his proof.

Consider the array partitioning problem. In it we are given values $$cost[A][B]$$ representing the cost of interval $$[a, b)$$, and wish to partition the array into intervals $$[0, x_{1}), [x_{1}, x_{2}), \dots, x_{k}, n)$$. Let $$DP[n][k]$$ denote the maximum sum $$\sum_{i = 0}^{k} cost[x_{i}][x_{i+1}]$$, where $$x_{0} = 0$$, $$x_{k+1} = n$$. We claim that $$DP[n][k]$$ is concave if the costs are Monge, that is, for all $$a \leq b \leq c \leq d$$ we have $$cost[a][d] + cost[b][c] \leq cost[a][c] + cost[b][d]$$.

First we'll show that our problem is an instance of the array partitioning problem with Monge costs. Set $$cost[a][b]$$ to be the length of the longest increasing subsequence in the interval $$[a, b)$$. Then $$DP[n][k]$$ for this instance of the array partitioning problem equals $$DP[n][k]$$ for our longest increasing subsequence problem. It remains to show that the costs are Monge.

Choose $$a \leq b \leq c \leq d$$, and take any LIS $$L_{a, d} = x_{1}, \dots, x_{cost[a][d]}$$ in the interval $$[a, d)$$, and any LIS $$L_{b, d} = y_{1}, \dots, y_{cost[b][c]}$$ in the interval $$[b, c)$$. We will combine them into two increasing subsequences in the intervals $$[a, c)$$ and $$[b, d)$$ of equal total length. To do this, let $$x_{i}$$ be the first $$x$$ and $$x_{j}$$ be the last $$x$$ in $$[b, c)$$. If $$x_{i} \leq y_{1}$$, set $$L_{a, c} = (x_{1}, \dots, x_{i-1}, y_{1}, \dots, y_{cost[b][c]})$$, $$L_{b, d} = (x_{i}, \dots, x_{cost[a][d]})$$. If $$x_{j} \geq y_{cost[b][c]}$$ do the same reversed. Otherwise, exists some $$t, h$$ s.t. $$y_{h} \leq x_{t} \leq y_{h+1}$$. Then set $$L_{a, c} = (x_{1}, \dots, x_{t}, y_{h+1}, \dots, y_{cost[b, c]})$$ and $$L_{b, d} = (y_{1}, \dots, y_{h}, x_{t+1}, \dots, x_{cost[a, d]})$$. Hence our cost array is Monge.

Now we'll show that the array partitioning problem with Monge cost is concave. Note that $$DP[n][k+2] - DP[n][k+1] \leq DP[n][k+1] - DP[n][k]$$ is the same inequality as $$DP[n][k+2] + DP[n][k] \leq 2 DP[n][k+1]$$. Take any partitions $$x_{0}, \dots, x_{k+3}$$ and $$y_{0}, \dots, y_{k+1}$$ with values $$DP[n][k+2]$$ and $$DP[n][k]$$ respectively. Take any $$0 \leq i \leq k$$ such that $$y_{i} \leq x_{i+1} \leq x_{i+2} \leq y_{i+1}$$. Such $$i$$ always exists, as some interval $$[y_{i}, y_{i+1}]$$ must be the first such that the last $$x$$ before the end of the interval, $$x_{j+2} \leq y_{i+1}$$ has $$j \geq i$$, thus $$x_{i+2} \leq x_{j+2} \leq y_{i+1}$$ and $$y_{i} \leq x_{i+1}$$ as otherwise the interval $$[y_{i-1}, y_{i}]$$ would contain $$x_{i+1}$$ contradicting the minimality of $$i$$.

We make the partitions $$y_{0}, \dots, y_{i}, x_{i+2}, \dots, x_{k+3}$$ and $$x_{0}, \dots, x_{i+1}, y_{i+1}, \dots, y_{k+1}$$, both of length $$k+1$$. What is the difference in total value? Most terms cancel, but in the sum of values of the original we have $$cost[x_{i+1}][x_{i+2}]$$ and $$cost[y_{i}][y_{i+1}]$$, while in the new one we have $$cost[y_{i}][x_{i+2}]$$ and $$cost[x_{i+1}][y_{i+1}]$$. But since $$y_{i} \leq x_{i+1} \leq x_{i+2} \leq y_{i+1}$$, by the Monge property $$cost[y_{i}][y_{i+1}] + cost[x_{i+1}][x_{i+2}] \leq cost[y_{i}][x_{i+2}] + cost[x_{i+1}][y_{i+1}]$$, hence the total value can only increase, and $$DP[n][k+2] + DP[n][k] \leq 2 DP[n][k+1]$$.

• Thanks for taking the time to give the algorithm too. Where can I find more about this? I'm just curious because searching "Lagrangian relaxation" doesn't give anything algorithm-related.
– user113632
Dec 25 '19 at 10:53
• In competitive programming this trick is called the "Aliens trick" after the problem "Aliens" from IOI 2016 where it appeared. Searching for it should find you many blogs on the topic. I've also seen this technique in an approximation algorithm for the k-median problem, see the paper "Approximation algorithms for metric facility location and k-median problems using the primal-dual schema and Lagrangian relaxation", Journal of the ACM 2001. Dec 25 '19 at 13:31
• How could I find the optimal $b$ in the end so I can calculate the answer?
– user113632
Dec 26 '19 at 19:09
• If for every position you know $DP'[i]$, $\min_{b}[i]$ and $\max_{b}[i]$, you can maintain the current endpoint and number of restarts, $i = n, b = k$, then build a maximum subsequence iterating backwards. For $j < i$, if $v_{j} \leq v_{i}$ and $DP'[j] = DP'[i] - 1$, $\min_{b}[j] \leq b \leq \max_{b}[j]$, include $j$ to your sequence and set $i = j$. If $v_{j} > v_{i}$ and $DP'[j] = DP'[i] - 1 + \lambda$, $\min_{b}[j] \leq b-1 \leq \max_{b}[j]$, include $j$, set $i = j$ and $b = b-1$. Dec 26 '19 at 20:16
• I added code to solve this problem to my answer, you can refer to it for implementation details. Dec 26 '19 at 22:30


text:
RandomBallCover  1.2.1 Hosted by GitHub
File List
Here is a list of all files with brief descriptions:
 2d_plot.cpp An application examining the resulting arrangement of the Random Ball Cover data structure on 2-D data 2d_plot_nn.cpp An application examining the resulting arrangement and operation of the Random Ball Cover data structure on 2-D data algorithms.cpp Defines classes that organize the execution of OpenCL kernels algorithms.hpp Declares classes that organize the execution of OpenCL kernels common.hpp Declares classes used by the OpenCL interface classes in cl_algo data_types.hpp Declarations of data types used by the Random Ball Cover data structure helper_funcs.cpp Helper functions for testing helper_funcs.hpp Declarations of helper functions for testing img_segmentation.cpp An application examining the potential of the Random Ball Cover data structure for image segmentation kinect_frame_grabber.cpp Grabs and stores in binary files an RGB and a Depth Kinect frame rbc_kernels.cl Kernels for building and accessing the Random Ball Cover data structure reduce_kernels.cl Kernels for performing the Reduce operation rgb_img_segmentation.cpp An application examining the potential of the Random Ball Cover data structure for image segmentation scan_kernels.cl Kernels for performing the Scan operation testsRBC.cpp Google Test Unit Tests for the kernels associated with the Random Ball Cover data structure testsRBC_RGBD.cpp Google Test Unit Tests for the kernels associated with the Random Ball Cover data structure testsReduce.cpp Google Test Unit Tests for the Reduce kernels testsScan.cpp Google Test Unit Tests for the Scan kernels