def reconstruct(self, denoised=True, niter=10, verbose=False):
        """
        Reconstruct the original image using iterative method with
        L1 regularization, because the denoising violates the exact inverse
        procedure.

        arguments:
        * denoised: whether use the denoised coefficients
        * niter: number of iterations
        """
        if denoised:
            decomposition = self.denoised
        else:
            decomposition = self.decomposition
        # L1 regularization
        lbd = 1.0
        delta = lbd / (niter - 1)
        # initial solution
        solution = self.reconstruct_ivst(denoised=denoised,
                positive_project=True)
        #
        iuwt = IUWT(level=self.level)
        iuwt.calc_filters()
        # iterative reconstruction
        if verbose:
            print("Iteratively reconstructing (%d times): " % niter,
                    end="", flush=True, file=sys.stderr)
        for i in range(niter):
            if verbose:
                print("%d..." % i, end="", flush=True, file=sys.stderr)
            tempd = self.data.copy()
            solution_decomp = []
            for scale in range(1, self.level+1):
                approx, detail = iuwt.transform(tempd, scale)
                approx_sol, detail_sol = iuwt.transform(solution, scale)
                # Update coefficients according to the significant supports,
                # which are acquired during the denosing precodure with FDR.
                sig = self.sig_supports[scale]
                detail_sol[sig] = detail[sig]
                detail_sol = self.soft_threshold(detail_sol, threshold=lbd)
                #
                solution_decomp.append(detail_sol)
                tempd = approx.copy()
                solution = approx_sol.copy()
            # last approximation (the two are the same)
            solution_decomp.append(approx)
            # reconstruct
            solution = iuwt.reconstruct(decomposition=solution_decomp)
            # discard all negative values
            solution[solution < 0] = 0.0
            #
            lbd -= delta
        if verbose:
            print("DONE!", flush=True, file=sys.stderr)
        #
        self.reconstruction = solution
        return self.reconstruction